text:
@mix @columns struct AlleleFrequency{AF}
    # Field              | Default  | Flat  | Bounds      | Description
    alleleFrequency::AF  | 0.5      | true  | (0.0, 1.0)  | "Allele frequency"
end

@mix @columns struct DeviationPhenotype{DP}
    # Field              | Default  | Flat  | Bounds      | Description
    deviationPhenotype::DP  | 0.5      | true  | (-1.0, 1.0)  | "Deviation of homozygous to the average phenotype"
end

@mix @columns struct DominanceDegree{DD}
    # Field              | Default  | Flat  | Bounds      | Description
    dominanceDegree::DD  | 0.5      | true  | (-1.0, 1.0)  | "AdditDegree of dominance between two alleles"
end

isprobvec(p::AbstractVector{<:Real}) =
    all(x -> x ≥ zero(x), p) && isapprox(sum(p), one(eltype(p)))

# struct AlleleFrequency{T<:AbstractVector{<:Real}}
#     p::T
#     AlleleFrequency{T}(p) where {T<:AbstractVector{<:Real}} = !isprobvec(p) ? throw(ArgumentError("p = $p is not a probability vector.")) : new{T}(p)
# end
#
# AlleleFrequency(1.0)
# AlleleFrequency([0.1, 0.9])
# AlleleFrequency([0.000001, 0.999999])
# AlleleFrequency([1])
# AlleleFrequency([0.1,0.1,0.1,0.1,0.1,0.2,0.2,0.1])

abstract type DeltaAlleleFrequencyRule{R,W} <: CellRule{R,W} end

DynamicGrids.precalcrules(rule::DeltaAlleleFrequencyRule, data) = rule

abstract type DeltaAlleleFrequencyMapRule{R,W} <: DeltaAlleleFrequencyRule{R,W} end

DynamicGrids.precalcrules(rule::DeltaAlleleFrequencyMapRule, data) = begin
    precalclayer(layer(rule, data), rule, data)
end

@Exposure @LC50 @HillCoefficient @DeviationPhenotype @DominanceDegree struct DeltaAlleleFrequencySurv{R,W} <: DeltaAlleleFrequencyRule{R,W} end

@inline applyrule(data, rule::DeltaAlleleFrequencySurv, alleleFrequency, args...) = begin
    alleleFrequency > zero(alleleFrequency) || return zero(alleleFrequency)
    @fastmath alleleFrequency +
        alleleFrequency * (1-alleleFrequency)*(rule.deviationPhenotype + rule.dominanceDegree*(1-2*alleleFrequency))*
        rule.hillcoefficient*(rule.exposure/rule.LC50)^rule.hillcoefficient / ( rule.LC50 *((rule.exposure/rule.LC50)^rule.hillcoefficient +1))
end

@Layers @LC50 @HillCoefficient @DominanceDegree @DeviationPhenotype struct DeltaAlleleFrequencySurvMap{R,W} <: DeltaAlleleFrequencyMapRule{R,W} end

@inline applyrule(data, rule::DeltaAlleleFrequencySurvMap, alleleFrequency, index, args...) = begin
    alleleFrequency > zero(alleleFrequency) || return zero(alleleFrequency)
    exposure = layer(rule, data, index)
    @fastmath alleleFrequency +
        alleleFrequency * (1-alleleFrequency)*(rule.deviationPhenotype + rule.dominanceDegree*(1-2*alleleFrequency))*
        rule.hillcoefficient*(exposure/rule.LC50)^rule.hillcoefficient / ( rule.LC50 *((exposure/rule.LC50)^rule.hillcoefficient +1))
end

@Layers @LC50 @HillCoefficient @DeviationPhenotype @DominanceDegree struct DeltaAlleleFrequencySurvMap_noFastmath{R,W} <: DeltaAlleleFrequencyMapRule{R,W} end

@inline applyrule(data, rule::DeltaAlleleFrequencySurvMap_noFastmath, alleleFrequency, index, args...) = begin
    alleleFrequency > zero(alleleFrequency) || return zero(alleleFrequency)
    exposure = layer(rule, data, index)
    alleleFrequency +
        alleleFrequency * (1-alleleFrequency)*(rule.deviationPhenotype + rule.dominanceDegree*(1-2*alleleFrequency))*
        rule.hillcoefficient*(exposure/rule.LC50)^rule.hillcoefficient / ( rule.LC50 *((exposure/rule.LC50)^rule.hillcoefficient +1))
end


text:
dupeelements(qq) = [qq[ceil(i)] for i in (1:(2*length(qq)))/2]

function update_statepos!(specs::Vector{EulerSpec})
	i = 1
	for spec in specs
		for j = 1:length(spec.statepos)
			spec.statepos[j] = i
			i += 1
		end
		#@show spec
	end
end

# given count_totals, which is unscaled sizes, compute a size which
# depends on the shape type, along with bounding box info.
# lb/ub for x/y params are half the maximum extent
# lb/ub for other parameters are usually [0,1]
function compute_shape_sizes(specs, count_totals, sizesum)
	@assert length(specs) == length(count_totals)
	shape_sizes = similar(count_totals, Float64)
	bounding_boxes = Float64[] #similar(count_totals, Float64)

	# rescale count_totals (areas) to sum to sizesum
	areas = sizesum * count_totals / sum(count_totals)

	for i in 1:length(specs)
		if specs[i].shape == :circle
			# area = pi r^2; r = sqrt(area/pi)
			shape_sizes[i] = sqrt(areas[i] / pi)
			append!(bounding_boxes, [shape_sizes[i], shape_sizes[i]])
		elseif specs[i].shape == :square
			# area = side^2; halfside = sqrt(area) / 2
			shape_sizes[i] = sqrt(areas[i]) / 2
			append!(bounding_boxes, [shape_sizes[i], shape_sizes[i]])
		elseif specs[i].shape == :triangle
			# area = side^2 * sqrt(3)/4; dist to vertex (bounding circle radius) = side/2 / cos(pi/6)
			side = sqrt(4 * areas[i] / sqrt(3))
			shape_sizes[i] = side / (2cos(pi/6))
			append!(bounding_boxes, [shape_sizes[i], shape_sizes[i]]) # I think side/2 is right, but get errors?
		elseif specs[i].shape == :rectangle
			# for rectangles, one of the parameters is the eccentricity, which
			# varies from 0 to 1. This maps to a height/width ratio that varies
			# from 1:4 to 4:1 (capped so that boundaries can be drawn). The
			# conversion formula is hw = 16 ^ (ecc - 0.5).
			# From there, height = sqrt(area) * hw
			# and width = sqrt(area) / hw. This means that shape_size for a
			# rectangle is sqrt(area)/2.
			shape_sizes[i] = sqrt(areas[i])/2
			append!(bounding_boxes, [2shape_sizes[i], 2shape_sizes[i], 0])
				# to give room depending on ecc; no need to bound eccentricity
		else
			error("unknown shape: ", specs[i].shape)
		end
	end
	shape_sizes, bounding_boxes
end

function make_euler_object(labels, counts, specs::Vector{EulerSpec}; sizesum = 1)
	target = DisjointSet(counts, labels)
	count_totals = vec(sum(counts,1))

	@assert length(labels) == length(count_totals)

	# use the spec to figure out how to translate to and from a state
	nparams = sum([length(spec.clamp) for spec in specs])

	update_statepos!(specs) # this lets us look the up parameters in the state vector

	# we have non-normalized areas, but need to compute a per-shape size (e.g., radius)
	# that indicates the maximum distance from the object center to the farthest edge
	shape_sizes, bounding_boxes = compute_shape_sizes(specs, count_totals, sizesum)
	# then, use those sizes to generate bounds

	# this forces the centers to not overlap with the boundary
	# and inserts any non-NaN clamps into the vector
	clamp_vec = [[sp.clamp for sp in specs]...]
	#@show clamp_vec
	#@show shape_sizes
	lb = ifelse(isnan(clamp_vec), bounding_boxes, clamp_vec)
	ub = ifelse(isnan(clamp_vec), 1 .- bounding_boxes, clamp_vec)
	@assert all(lb .<= ub)

	# return: state vector, state object (with bounds, closure, etc)
	eo = EulerObject(nparams, labels, lb, ub, shape_sizes, target, specs, identity)
	eo.evalfn = (x,g) -> begin
		@assert length(g) == 0
		cost = VennEuler.eval_euler_state(eo, x)
		#println(x, " (", cost, ")")
		cost
	end
	eo
end
make_euler_object(labels, counts, spec::EulerSpec; q...) =
	make_euler_object(labels, counts, EulerSpec[deepcopy(spec)::EulerSpec for x in 1:length(labels)]; q...)

function random_state(eo::EulerObject)
	rand(eo.nparams) .* (eo.ub .- eo.lb) .+ eo.lb
end

ext: jl


text:
@testset "1106.parsing-a-boolean-expression.jl" begin
    @test parse_bool_expr("|(&(t,f,t),!(t))") == false
    @test parse_bool_expr("&(t,f)") == false
    @test parse_bool_expr("|(f,t)") == true
    @test parse_bool_expr("!(f)") == true
end


text:
"""
Retrieve static data. Multiple tickers and multiple fields can be specified.
A data frame will be returned with securities as rows and fields as columns.

Bulk data fields return data frames and cannot be combined with other fields.
If multiple securities are given, a dictionary ticker => DataFrame is returned.
These special fields are not detected automatically, a few are predefined but
you may need to update the variable Bloomie.bulk_fields in file src/response.jl

#Examples
```julia
bdp(["AAPL US Equity","MSFT US Equity"],["PX_LAST","VOLUME"])
bdp("AAPL US Equity","BEST_PE_RATIO",overrides=Dict("BEST_FPERIOD_OVERRIDE"=>"1BF"))
bdp(["FB US Equity","TWTR US Equity"],"ERN_ANN_DT_AND_PER")
bdp("TSLA US Equity","BEST_ANALYST_RECS_BULK",overrides=Dict("END_DATE_OVERRIDE"=>"20151231"))
bdp(["INDU Index","TRAN Index"],"INDX_MEMBERS")
bdp("INDU Index","INDX_MWEIGHT_HIST",overrides=Dict("END_DT"=>"20101231"))
```
"""
function bdp(tickers,fields;overrides=Dict())
    global session
    global services
    single_security = false
    if !isa(tickers,Vector); tickers=[tickers]; single_security=true; end
    if !isa(fields,Vector); fields=[fields]; end
    service = services["refdata"]
    request = @with_pointer request blpapi_Service_createRequest(service,request,"ReferenceDataRequest")
    elements = blpapi_Request_elements(request)
    try
        securities_ = @with_pointer securities_ blpapi_Element_getElement(elements,securities_,"securities",C_NULL)
        foreach(ticker -> blpapi_Element_setValueString(securities_,ticker,-1), tickers)
        fields_ = @with_pointer fields_ blpapi_Element_getElement(elements,fields_,"fields",C_NULL)
        foreach(field -> blpapi_Element_setValueString(fields_,field,-1), fields)
        if length(overrides) > 0
            overrides_ = @with_pointer overrides_ blpapi_Element_getElement(elements,overrides_,"overrides",C_NULL)
            for (k,v) in overrides
                new_ = @with_pointer new_ blpapi_Element_appendElement(overrides_,new_)
                blpapi_Element_setElementString(new_, "fieldId", C_NULL, string(k))
                blpapi_Element_setElementString(new_, "value", C_NULL, string(v))
            end
        end
        correlation_id = blpapi_CorrelationId(0,0)
        blpapi_Session_sendRequest(session,request,Ref(correlation_id),C_NULL,C_NULL,C_NULL,0)
    finally
        blpapi_Request_destroy(request)
    end
    out=process_tree(retrieve_response())
    if typeof(out[1]) <: Dict
        data = [row for row in filter(x->typeof(x)<:Dict,out)]
        nodata = filter(x->!(typeof(x)<:Dict),out)
        map(x->warn(Printf.@sprintf "No data retrieved for %s" x[1]),nodata)
        if length(data) > 0
            DataFrames.DataFrame(data)
        else
            nothing
        end
    elseif typeof(out) <: Array
        out = vcat(out...)
        data = filter(x->x[2]!=nothing,out)
        nodata = filter(x->x[2]==nothing,out)
        map(x-> @warn Printf.@sprintf "No data retrieved for %s" x[1] , nodata)
        if length(data) > 0
            Dict(data)
            if single_security
                data[1][2]
            else
                Dict(data)
            end
        else
            nothing
        end
    else
        error("this should not happen")
    end
end


text:
using TetGen
using TetGen: TetgenIO
using GeometryBasics
using GeometryBasics: Mesh, QuadFace
# Construct a cube out of Quads
points = Point{3, Float64}[
    (0.0, 0.0, 0.0), (2.0, 0.0, 0.0),
    (2.0, 2.0, 0.0), (0.0, 2.0, 0.0),
    (0.0, 0.0, 12.0), (2.0, 0.0, 12.0),
    (2.0, 2.0, 12.0), (0.0, 2.0, 12.0)
]

# using Plots

# plot(points, markershape = :o)

facets = QuadFace{Cint}[
    1:4,
    5:8,
    [1,5,6,2],
    [2,6,7,3],
    [3, 7, 8, 4],
    [4, 8, 5, 1]
]
# markers = Cint[-1, -2, 0, 0, 0, 0]
# attach some additional information to our faces!
mesh = Mesh(points, facets)
# result = tetrahedralize(mesh, "vpq1.414a0.1")
# using FileIO, MeshIO
# # and save it in your favorite format
# FileIO.save("column.ply", mesh)
# # using GLMakie, AbstractPlotting
# # GLMakie.mesh(normal_mesh(result), color=(:blue, 0.1), transparency=true)
# # GLMakie.wireframe!(result)

#try a triangular mesh

points = Point{3, Float64}[
    (0.0, 0.0, 0.0), (0.0, 1.0, 0.0),
    (0.0, 0.0, 1.0)]

# using Plots

# plot(points, markershape = :o)

facets = TriangleFace{Cint}[
    [1, 2, 3]]
# markers = Cint[-1, -2, 0, 0, 0, 0]
# attach some additional information to our faces!
mesh = GeometryBasics.Mesh(points, facets)


text:
#Copyright (c) 2015 Rex Computing and Isaac Yonemoto

#see LICENSE.txt

#this work was supported in part by DARPA Contract D15PC00135

#profiling zero tests.

include ("../../unum.jl")

#run in unum{4,6} environment.  30% of the unums will have a ubit set.
#10% of the unums will be zero, 10% of them will be negative zero.
#10% of the unums will be a subnormal.

#calculates a random number between 0..i
function randint(i)
end

function randU(T)
  x = convert(T, randn())
  dieroll = rand()
  if (dieroll < 0.1)
    zero(T)
  elseif (dieroll < 0.2)
    -zero(T)
  elseif (dieroll < 0.3)
    T(z16, z16,
  else
    convert(T, rand())
  end
end

count = 100000

function profile(T)
  #generate the array.
  array = [randU() for i = 1:count]

  println("naive ordering")
  tic()
    map(iszero_naive, array)
  toc()

  println("optimized ordering")
  tic()
    map(iszero_naive, array)
  toc()
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
#=
Examples of conditional statements.
=#

x = 10
y = 17.2

# An if-else statement
if x > y
    println("the x > y")
else
    println("the x < y or x and y are equal")
end

# An if-elseif-else statement
if x > y
    println("the x > y")
elseif x < y
    println("the x < y")
else
    println("x and y are equal")
end

# Using the ternary operator for an if-else statement
(x > y) ? println("the x > y") : println("the x < y")

# Can also write the above as follows
println((x > y) ? "the x > y" : "the x < y")

max_stars_repo_name: wigging/julia-computing

max_issues_repo_name: wigging/julia-computing

max_forks_repo_name: wigging/julia-computing


text:
const julia_sind = @cfunction(sind, Cdouble, (Cdouble,))

# Mutually recursive fibonacci implementation which bounces back and
# forth between Julia and C.
function mutual_fibonacci(n::Cint, callback::Ptr{Cvoid})::Cint
    if n <= 1
        return n
    end

    return ccall(callback, Cint, (Cint,), n - 1) +
           ccall(callback, Cint, (Cint,), n - 2)
end

const mutual_fibonacci_ptr = @cfunction(mutual_fibonacci, Cint, (Cint, Ptr{Cvoid}))

# Example of receiving an array of data. This function adds 1 to the
# first element, 2 to the second, and so on.
function add_element_number(x::Ptr, n::Integer)
    for i = 1:n
        element = unsafe_load(x, i)
        unsafe_store!(x, element + i, i)
    end
    return
end

# @cfunction can instantiate pointers to multiple specializations of
# the function for different argument types.
const add_element_number_int = @cfunction(add_element_number, Cvoid,
                                          (Ptr{Cint}, Cint))
const add_element_number_float = @cfunction(add_element_number, Cvoid,
                                            (Ptr{Cfloat}, Cint))


text:
# push!(LOAD_PATH,"E:/GitHub/ScTenifold.jl/src/");
using ScTenifold                          # obviously the tests use the ScTenifoldNet module...
using Test                                # and the Base.Test module...
tests = ["code_test1", "code_test2"]      # the test file names are stored as strings...
for t in tests
  include("$(t).jl")                      # ... so that they can be evaluated in a loop
end

max_stars_repo_name: jamesjcai/ScTenifold

max_issues_repo_name: jamesjcai/ScTenifold

max_forks_repo_name: jamesjcai/ScTenifold


text:
using Revise

# using BSplines
# using BayesNegativeBinomial
using BNPRegressionGGA2021
using Distributions
using Gadfly
using Random
# using Statistics
# using Test
# using LinearAlgebra

#=================================#
# Example 1 - 1 trivial predictor #
#=================================#

function simulate_sample_01(rng, N0, N1)
    dy = MixtureModel(LogNormal, [(0.0, 1.0), (1.0, 0.25)], [0.6, 0.4]);
    X0 = ones(N0, 1)
    X1 = ones(N1, 1)
    z0 = rand(rng, dy, N0)
    ỹ0 = deepcopy(z0)
    ỹ1 = rand(rng, Exponential(10), N0)
    ỹ0 = min.(ỹ0, ỹ1)
    event = ỹ0 .< ỹ1
    y1 = LinRange(0, 7, N1) |> collect
    return dy, event, z0, ỹ0, X0, y1, X1
end

N0, N1 = 500, 50;
rng = MersenneTwister(2);
dy, event, z0, ỹ0, X0, y1, X1 = simulate_sample_01(rng, N0, N1);
mean(event)
smpl = BNPRegressionGGA2021.ErlangSampler(; rng, ỹ0, X0, y1, X1);
chainf, chainβ = BNPRegressionGGA2021.sample(rng, smpl; mcmcsize = 20000, burnin = 10000);
plot(
    layer(x = y1, y = mean(chainf), Geom.line, color=["bnp"]),
    layer(x = y1, y = pdf.(dy, y1), Geom.line, color=["true"]),
    layer(x = z0, Geom.histogram(density = true), color=["histogram"])
)
# smpl.φ
# smpl.λ
# unique(smpl.d)
# maximum(smpl.n)
# dy = Erlang(ceil(Int, smpl.φ[3]), 1.0 / smpl.λ[])
# mean(dy)
# std(dy)
# mean(dy)
# std(dy)

#=====================================#
# Example 2 - 1 non-trivial predictor #
#=====================================#

function simulate_sample_01(rng, N0, N1)
    dy1 = LogNormal(0, 0.5)
    dy2 = MixtureModel(LogNormal, [(0.0, 1.0), (1.0, 0.25)], [0.6, 0.4]);
    x0 = [ones(N0) rand(rng, [1, 2], N0)
    X1 = kron(ones(N1), [1 1; 1 2])
    z0 = [X0[i, 2] == 2 ? rand(rng, dy2) : rand(rng, dy1) for i in 1:N0]
    y1 = LinRange(0, 6, N1) |> collect  |> x -> repeat(x, 2)
    ỹ0 = deepcopy(z0)
    ỹ1 = rand(rng, Exponential(30), N0)
    event = ỹ0 .< ỹ1
    return dy1, dy2, event, z0, ỹ0, X0, y1, X1
end

N0, N1 = 1000, 50;
rng = MersenneTwister(1);
dy1, dy2, event, z0, ỹ0, X0, y1, X1 = simulate_sample_01(rng, N0, N1);
mean(event)
grid = LinRange(0, 6, N1) |> collect;
plot(
    layer(x = grid, y = pdf.(dy2, grid), Geom.line, color=["kdensity2"]),
    layer(x = grid, y = pdf.(dy1, grid), Geom.line, color=["kdensity1"]),
)
smpl = BNPRegressionGGA2021.ErlangSampler(; ỹ0, X0, y1, X1);
chainf, chainβ = BNPRegressionGGA2021.sample(rng, smpl; mcmcsize = 20000, burnin = 10000);
plot(x = y1, y = mean(chainf), color = string.(X1[:, 2]), Geom.line)
plot(x = z0, color = string.(X0[:, 2]), Geom.histogram(density=true, bincount = 50))

#===================================================#
# Example 3 - 1 trivial and 1 non-trivial predictor #
#===================================================#

function simulate_sample_01(rng, N0, N1)
    dy1 = LogNormal(0, 0.5)
    dy2 = MixtureModel(LogNormal, [(0.0, 1.0), (1.0, 0.25)], [0.6, 0.4]);
    X0 = [ones(N0) rand(rng, [1, 2], N0, 2)]
    X1 = kron(ones(N1), [1 1 1; 1 2 1; 1 1 2; 1 2 2])
    z0 = [X0[i, 2] == 2 ? rand(rng, dy2) : rand(rng, dy1) for i in 1:N0]
    ỹ0 = deepcopy(z0)
    ỹ1 = rand(rng, Exponential(30), N0)
    event = ỹ0 .< ỹ1
    y1 = LinRange(0, 6, N1) |> collect  |> x -> repeat(x, 4)
    return dy1, dy2, event, z0, ỹ0, X0, y1, X1
end

N0, N1 = 1000, 50;
rng = MersenneTwister(2);
dy1, dy2, event, z0, ỹ0, X0, y1, X1 = simulate_sample_01(rng, N0, N1);
ygrid = LinRange(0, 6, N1) |> collect;
plot(
    layer(x = ygrid, y = pdf.(dy2, ygrid), Geom.line, color=["pdf2"]),
    layer(x = ygrid, y = pdf.(dy1, ygrid), Geom.line, color=["pdf1"]),
)
plot(
    layer(x = ygrid, y = pdf.(dy2, ygrid) ./ (1.0 .- cdf.(dy2, ygrid)), Geom.line, color=["pdf2"]),
    layer(x = ygrid, y = pdf.(dy1, ygrid) ./ (1.0 .- cdf.(dy1, ygrid)), Geom.line, color=["pdf1"]),
)
smpl = BNPRegressionGGA2021.ErlangSampler(; rng, ỹ0, X0, y1, X1);
chainf, chainβ = BNPRegressionGGA2021.sample(rng, smpl; mcmcsize = 20000, burnin = 10000);
plot(x = y1, y = mean(chainf), color = string.(X1[:, 2]) .* " + " .* string.(X1[:, 3]), Geom.line)

#===================================================#
# Example 4 - 1 trivial and 4 non-trivial predictor #
#===================================================#

function simulate_sample_01(rng, N0, N1)
    dy1 = LogNormal(0, 0.5)
    dy2 = MixtureModel(LogNormal, [(0.0, 1.0), (1.0, 0.25)], [0.6, 0.4]);
    X0 = [ones(N0) rand(rng, [1, 2], N0, 5)]
    X1 = kron(ones(N1), [1 1 1 ones(1, 3); 1 2 1 ones(1, 3); 1 1 2 ones(1, 3); 1 2 2 ones(1, 3)])
    z0 = [X0[i, 2] == 2 ? rand(rng, dy2) : rand(rng, dy1) for i in 1:N0]
    ỹ0 = deepcopy(z0)
    ỹ1 = rand(rng, Exponential(30), N0)
    event = ỹ0 .< ỹ1
    y1 = LinRange(0, 6, N1) |> collect  |> x -> repeat(x, size(X1, 1) ÷ N1)
    return dy1, dy2, event, z0, ỹ0, X0, y1, X1
end

N0, N1 = 1000, 50;
rng = MersenneTwister(2);
dy1, dy2, event, z0, ỹ0, X0, y1, X1 = simulate_sample_01(rng, N0, N1);
ygrid = LinRange(0, 6, N1) |> collect;
plot(
    layer(x = ygrid, y = pdf.(dy2, ygrid), Geom.line, color=["pdf2"]),
    layer(x = ygrid, y = pdf.(dy1, ygrid), Geom.line, color=["pdf1"]),
)
plot(
    layer(x = ygrid, y = pdf.(dy2, ygrid) ./ (1.0 .- cdf.(dy2, ygrid)), Geom.line, color=["pdf2"]),
    layer(x = ygrid, y = pdf.(dy1, ygrid) ./ (1.0 .- cdf.(dy1, ygrid)), Geom.line, color=["pdf1"]),
)
smpl = BNPRegressionGGA2021.ErlangSampler(; rng, ỹ0, X0, y1, X1);
chainf, chainβ = BNPRegressionGGA2021.sample(rng, smpl; mcmcsize = 20000, burnin = 10000);
plot(x = y1, y = mean(chainf), color = string.(X1[:, 2]) .* " + " .* string.(X1[:, 3]), Geom.line)


text:
# Analysis tools

export
    visualize,
    displacements,
    distances,
    rdf,
    energy,
    kinetic_energy,
    potential_energy,
    velocity_autocorr

"""
    visualize(coord_logger, box_size, out_filepath; <keyword arguments>)

Visualize a simulation as an animation.
This function is only available when GLMakie is imported.
It can take a while to run, depending on the length and size of the simulation.

# Arguments
- `connections=Tuple{Int, Int}[]`: pairs of atoms indices to link with bonds.
- `connection_frames`: the frames in which bonds are shown. Is a list of the
    same length as the number of frames, where each item is a list of `Bool`s of
    the same length as `connections`. Defaults to always `true`.
- `trails::Integer=0`: the number of preceding frames to show as transparent
    trails.
- `framerate::Integer=30`: the frame rate of the animation.
- `color=:purple`: the color of the atoms. Can be a single color or a list of
    colors of the same length as the number of atoms.
- `connection_color=:orange`: the color of the bonds. Can be a single color or a
    list of colors of the same length as `connections`.
- `markersize=20.0`: the size of the atom markers.
- `linewidth=2.0`: the width of the bond lines.
- `transparency=true`: whether transparency is active on the plot.
- `kwargs...`: other keyword arguments are passed to the plotting function.
"""
function visualize end

"""
    displacements(coords, box_size)

Get the pairwise vector displacements of a set of coordinates, accounting for
the periodic boundary conditions.
"""
function displacements(coords, box_size)
    n_atoms = length(coords)
    coords_rep = repeat(reshape(coords, n_atoms, 1), 1, n_atoms)
    diffs = vector.(coords_rep, permutedims(coords_rep, (2, 1)), (box_size,))
    return diffs
end

"""
    distances(coords, box_size)

Get the pairwise distances of a set of coordinates, accounting for the periodic
boundary conditions.
"""
distances(coords, box_size) = norm.(displacements(coords, box_size))

"""
    rdf(coords, box_size; npoints=200)

Get the radial distribution function of a set of coordinates.
This describes how density varies as a function of distance from each atom.
Returns a list of distance bin centres and a list of the corresponding
densities.
"""
function rdf(coords, box_size; npoints::Integer=200)
    n_atoms = length(coords)
    dims = length(first(coords))
    dists = distances(coords, box_size)
    dists_vec = [dists[i, j] for i in 1:n_atoms, j in 1:n_atoms if j > i]
    dist_unit = unit(first(dists_vec))
    kd = kde(ustrip.(dists_vec), npoints=npoints)
    ρ = n_atoms / reduce(*, box_size)
    if dims == 3
        normalizing_factor = 4π .* ρ .* step(kd.x) .* kd.x .^ 2 .* dist_unit .^ 3
    elseif dims == 2
        normalizing_factor = 2π .* ρ .* step(kd.x) .* kd.x .* dist_unit .^ 2
    end
    bin_centres = collect(kd.x) .* dist_unit
    density_weighted = kd.density ./ normalizing_factor
    return bin_centres, density_weighted
end

"""
    energy(s, neighbors=nothing)

Compute the total energy of the system.
"""
energy(s, neighbors=nothing) = kinetic_energy(s) + potential_energy(s, neighbors)

"""
    kinetic_energy(s)

Compute the kinetic energy of the system.
"""
function kinetic_energy(s::System)
    ke = sum(i -> s.atoms[i].mass * dot(s.velocities[i], s.velocities[i]) / 2, axes(s.atoms, 1))
    # Convert energy to per mol if required
    if dimension(s.energy_unit) == u"𝐋^2 * 𝐌 * 𝐍^-1 * 𝐓^-2"
        T = typeof(ustrip(ke))
        return uconvert(s.energy_unit, ke * T(Unitful.Na))
    else
        return uconvert(s.energy_unit, ke)
    end
end

"""
    potential_energy(s, neighbors=nothing)

Compute the potential energy of the system.
"""
function potential_energy(s::System, neighbors=nothing)
    n_atoms = length(s)
    potential = zero(ustrip(s.box_size[1])) * s.energy_unit

    for inter in values(s.general_inters)
        if inter.nl_only
            @inbounds for ni in 1:neighbors.n
                i, j, weight_14 = neighbors.list[ni]
                if weight_14
                    potential += potential_energy(inter, s, i, j, true)
                else
                    potential += potential_energy(inter, s, i, j)
                end
            end
        else
            for i in 1:n_atoms
                for j in (i + 1):n_atoms
                    potential += potential_energy(inter, s, i, j)
                end
            end
        end
    end

    for inter_list in values(s.specific_inter_lists)
        potential += potential_energy(inter_list, s.coords, s.box_size)
    end

    return uconvert(s.energy_unit, potential)
end

@views function potential_energy(inter_list::InteractionList2Atoms, coords, box_size)
    return sum(potential_energy.(inter_list.inters, coords[inter_list.is], coords[inter_list.js],
                                    (box_size,)))
end

@views function potential_energy(inter_list::InteractionList3Atoms, coords, box_size)
    return sum(potential_energy.(inter_list.inters, coords[inter_list.is], coords[inter_list.js],
                                    coords[inter_list.ks], (box_size,)))
end

@views function potential_energy(inter_list::InteractionList4Atoms, coords, box_size)
    return sum(potential_energy.(inter_list.inters, coords[inter_list.is], coords[inter_list.js],
                                    coords[inter_list.ks], coords[inter_list.ls], (box_size,)))
end

"""
    velocity_autocorr(vl, first_ind, last_ind)

Calculates the autocorrelation function of velocity from the velocity logger.
This helps characterize the similarity between velocities observed at different
time instances.
"""
function velocity_autocorr(vl::VelocityLogger, first_ind::Integer=1, last_ind::Integer=length(vl.velocities))
    n_atoms = length(first(vl.velocities))
    return dot(vl.velocities[first_ind], vl.velocities[last_ind]) / n_atoms
end


text:
using RobotDynamics
using Test
using StaticArrays

struct MyModel <: AbstractModel end
model = MyModel()
@test_throws ErrorException state_dim(model)
@test_throws ErrorException control_dim(model)

model = Cartpole()
x,u = zeros(model)
@test sum(x) == 0
@test sum(u) == 0
x,u = rand(model)
@test sum(x) != 0
@test sum(u) != 0
xdot = dynamics(model, x, u)
@test sum(x) != 0

x,u = zeros(Float32,model)
@test x isa SVector{4,Float32}
@test u isa SVector{1,Float32}

x,u = rand(Float32,model)
@test x isa SVector{4,Float32}
@test u isa SVector{1,Float32}

x,u = fill(model, 10.0)
@test x ≈ @SVector fill(10,4)
@test u ≈ @SVector fill(10,1)

n,m = size(model)
dt = 0.1
F = zeros(n,n+m)
z = KnotPoint(x,u,dt)
jacobian!(F, model, z)
@test sum(F) != 0

D = RobotDynamics.DynamicsJacobian(n,m)
jacobian!(D, model, z)
@test D.A == F[:,1:n]
@test D.B ≈ F[:,n .+ (1:m)]

@test discrete_dynamics(RK3, model, x, u, 0.0, dt) ≈
    discrete_dynamics(RK3, model, z)
@test discrete_dynamics(RK3, model, z) ≈ discrete_dynamics(model, z)

z_ = KnotPoint(zero(x),zero(m),dt)
RobotDynamics.propagate_dynamics(RK3, model, z_, z)
@test state(z_) ≈ discrete_dynamics(model, z)

F = zeros(n,n+m)
discrete_jacobian!(RK3, F, model, z)

tmp = [RobotDynamics.DynamicsJacobian(n,m) for k = 1:3]
jacobian!(RK3, D, model, z, tmp)
@test D.A ≈ F[1:n,1:n]
@test D.B ≈ F[1:n,n .+ (1:m)]
@test sum(F) != 0
@test F[1] == 1

# Error state
x0 = rand(model)[1]
@test RobotDynamics.state_diff(model, x, x0) ≈ x - x0
@test state_diff_size(model) == state_dim(model)
G = zeros(n,n)
RobotDynamics.state_diff_jacobian!(G, model, z)
@test G ≈ I(n)
G .= 0
RobotDynamics.state_diff_jacobian!(G, model, x)
@test G ≈ I(n)
@test RobotDynamics.state_diff_jacobian(model, x) ≈ G

# @btime discrete_jacobian!($RK3, $F, $model, $z)
# @btime jacobian!($RK3, $D, $model, $z, $tmp)


text:
import SortingAlgorithms: uint_mapping
import Base.Ordering
# sorts vs using radix sort but skip the last few bits as defined by skipbits
function sort32!(vs::AbstractVector{T}, lo::Int=1, hi::Int=length(vs), o::Ordering=Base.Forward, ts=similar(vs); skipbits = 32, RADIX_SIZE = 11, RADIX_MASK = 0x000007ff) where T
    # Input checking
    if lo >= hi;  return vs;  end

    # Make sure we're sorting a bits type
    # T = Base.Order.ordtype(o, vs)

    # Init
    iters = ceil(Integer, sizeof(T)*8/RADIX_SIZE)
    bin = zeros(UInt32, 2^RADIX_SIZE, iters)
    if lo > 1;  bin[1,:] = lo-1;  end

    # Histogram for each element, radix
    for i = lo:hi
        v = uint_mapping(o, vs[i])
        for j = skipbits÷RADIX_SIZE+1:iters
            idx = Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK) + 1
            @inbounds bin[idx,j] += 1
        end
    end

    # Sort!
    swaps = 0
    len = hi-lo+1
    for j = skipbits÷RADIX_SIZE+1:iters
        # Unroll first data iteration, check for degenerate case
        v = uint_mapping(o, vs[hi])
        idx = Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK) + 1

        # are all values the same at this radix?
        if bin[idx,j] == len;  continue;  end

        cbin = cumsum(bin[:,j])
        ci = cbin[idx]
        ts[ci] = vs[hi]
        cbin[idx] -= 1

        # Finish the loop...
        @inbounds for i in hi-1:-1:lo
            v = uint_mapping(o, vs[i])
            idx = Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK) + 1
            ci = cbin[idx]
            ts[ci] = vs[i]
            cbin[idx] -= 1
        end
        vs,ts = ts,vs
        swaps += 1
    end

    if isodd(swaps)
        vs,ts = ts,vs
        for i = lo:hi
            vs[i] = ts[i]
        end
    end
    vs
end

function load_uint2(s::String, skipbytes)
    ss = sizeof(s)
    if skipbytes > ss
        return UInt(0)
    elseif skipbytes + 4 > ss
        s1 = ((s |> pointer |> Ptr{UInt32}) + skipbytes) |> unsafe_load
        s2 = Base.zext_int(UInt64, s1) |> ntoh
        extra_bits_to_shift = 8(skipbytes + 4 - ss)
        return (s2 >> (32+extra_bits_to_shift)) << (32+extra_bits_to_shift)
    else
        s1 = ((s |> pointer |> Ptr{UInt32}) + skipbytes) |> unsafe_load
        return Base.zext_int(UInt64, s1) |> ntoh
    end
end

# this version is smaller
function fsortperm2(svec::AbstractVector{String})
    strlen = maximum(sizeof, svec)
    strlen = max(strlen-4, 0)
    underly_bits = load_uint2.(svec, strlen) .| collect(UInt32(1):UInt32(length(svec)))
    underly_bits = sort32!(underly_bits, 1, length(underly_bits),  Base.Forward)
    ss = (underly_bits .<< 32) .>> 32

    while strlen > 0
        strlen = max(strlen-4, 0)
        underly_bits = load_uint2.(@view(svec[ss]), strlen) .| ss
        underly_bits = sort32!(underly_bits, 1, length(underly_bits),  Base.Forward)
        ss = (underly_bits .<< 32) .>> 32
    end
    return ss
end

# too slow
function fsortperm4(svec::AbstractVector{String}, ::Type{T} = UInt) where T<:Unsigned
    strlen = maximum(sizeof, svec)
    strlen = max(strlen-sizeof(T), 0)

    l = length(svec)

    pairs = collect(enumerate(svec))

    sort!(pairs, by=x->load_uint(T, x[2], strlen), alg = RadixSort)

    while strlen > 0
        strlen = max(strlen-sizeof(T), 0)
        sort!(pairs, by=x->load_uint(T, x[2], strlen), alg = RadixSort)
    end
    return [x[1] for x in pairs]
end

# sort pointer but need to keep sizeof, so a bit slow
function fsortperm6(svec::AbstractVector{String}, ::Type{T} = UInt) where T<:Unsigned
    strlen = maximum(sizeof, svec)
    strlen = max(strlen-sizeof(T), 0)

    l = length(svec)

    pairs = Vector{Tuple{T, Int, Ptr{UInt}}}(l)
    for (i, svec1) = enumerate(svec)
        pairs[i] = tuple(i, sizeof(svec1), pointer(svec1) |> Ptr{UInt})
    end

    sort!(pairs, by=x->load_uint(x[3],x[2],strlen), alg = RadixSort)

    while strlen > 0
        strlen = max(strlen-sizeof(T), 0)
        sort!(pairs, by=x->load_uint(x[3],x[2],strlen), alg = RadixSort)
    end
    return [x[1] for x in pairs]
end


text:
using NearestNeighbors
using StaticArrays

if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

import Distances: Metric, evaluate
immutable CustomMetric1 <: Metric end
evaluate(::CustomMetric1, a::AbstractVector, b::AbstractVector) = maximum(abs(a - b))
function NearestNeighbors.interpolate{V <: AbstractVector}(::CustomMetric1,
                                                           a::V,
                                                           b::V,
                                                           x,
                                                           d,
                                                           ab)
    idx = (abs(b - a) .>= d - x)
    c = copy(a)
    c[idx] = (1 - x / d) .* a[idx] + (x / d) .* b[idx]
    return c, true
end
immutable CustomMetric2 <: Metric end
evaluate(::CustomMetric2, a::AbstractVector, b::AbstractVector) = norm(a - b) / (norm(a) + norm(b))

# TODO: Cityblock()
const metrics = [Chebyshev(), Euclidean(), Minkowski(3.5)]
const fullmetrics = [metrics; Hamming(); CustomMetric1(); CustomMetric2()]
const trees = [KDTree, BallTree]
const trees_with_brute = [BruteTree; trees]

include("test_knn.jl")
include("test_inrange.jl")
include("test_monkey.jl")
include("test_datafreetree.jl")


text:
module test_01_OneScalarSignal

using ModiaResult
using OrderedCollections
ModiaResult.@usingModiaPlot

t = range(0.0, stop=10.0, length=100)

result = OrderedDict{String,Any}()

result["time"] = t
result["phi"]  = sin.(t)

println("\n... test_01_OneScalarSignal.jl:\n")
ModiaResult.printResultInfo(result)

plot(result, "phi", heading="sine(time)")

end

max_stars_repo_name: ModiaSim/ModiaResult

max_issues_repo_name: ModiaSim/ModiaResult

max_forks_repo_name: ModiaSim/ModiaResult


text:
module StringSplitting

# split string on comma returning first element
fruit_string = "apple, banana, carrot"
first_fruit = strip(split(fruit_string, ",")[1])

# split string on comma returning all elements
fruit_string = "apple, banana, carrot"
fruits = map(strip, split(fruit_string, ","))

println(fruit_string)
println(first_fruit)
println(fruits)

end # module

max_stars_repo_name: curtisalexander/saspyrilia-sandboxes

max_issues_repo_name: curtisalexander/saspyrilia-sandboxes

max_forks_repo_name: curtisalexander/saspyrilia-sandboxes


text:
using Test
using WavePropBase
using WavePropBase.Geometry
using WavePropBase.Interpolation
using StaticArrays

@testset "Generic nodes" begin
    cheb_nodes = (n) -> [cos((2*k+1)/2n * π) for k in 0:n-1]
    # 1d
    nodes = cheb_nodes(10)
    f = x->cos(x)
    vals = map(f,nodes)
    p = TensorLagInterp(vals,nodes)
    xtest = 0.1
    @test f(xtest) ≈ p(xtest)
    # 2d
    nx = 10
    ny = 12
    x   = [0.5+0.5cos((2k-1)*π/2nx) for k in 1:nx] # Chebyshev nodes
    y   = [0.5+0.5cos((2k-1)*π/2ny) for k in 1:ny] # Chebyshev nodes
    f   = (x) -> cos(x[1]*x[2])
    vals = [f((x,y)) for x in x, y in y]
    p   = TensorLagInterp(vals,(x,y))
    xtest = (0.1,0.2)
    @test p(xtest) ≈ f(xtest)
    # 3d
    nx,ny,nz = 10, 12, 8
    x   = [0.5+0.5cos((2k-1)*π/2nx) for k in 1:nx] # Chebyshev nodes
    y   = [0.5+0.5cos((2k-1)*π/2ny) for k in 1:ny] # Chebyshev nodes
    z   = [0.5+0.5cos((2k-1)*π/2nz) for k in 1:nz] # Chebyshev nodes
    f   = (x) -> cos(x[1]*x[2])*x[3]
    vals = [f((x,y,z)) for x in x, y in y, z in z]
    p   = TensorLagInterp(vals,(x,y,z))
    xtest = (0.1,0.2,-0.1)
    @test p(xtest) ≈ f(xtest)
end

@testset "Cheb nodes" begin
    # 1,2,3 dimensions
    for d in 1:3
        p = ntuple(i->10+i,d)
        vals   = zeros(p)
        nodes1d   = map(p->cheb1nodes(p,-1,1),p)
        weights1d = map(p->cheb1weights(p),p)
        poly    = TensorLagInterp(vals,nodes1d,weights1d)
        f     = x->cos(prod(x))
        for I in CartesianIndices(poly)
            vals[I] = f(interpolation_nodes(poly,I))
        end
        xtest = ntuple(i->0.1,d)
        @test f(xtest) ≈ poly(xtest)
    end
    # 1,2,3 dimensions
    for d in 1:3
        p = ntuple(i->10+i,d)
        vals   = zeros(p)
        nodes1d   = map(p->cheb2nodes(p,-1.3,1.5),p)
        weights1d = map(p->cheb2weights(p),p)
        poly    = TensorLagInterp(vals,nodes1d,weights1d)
        f     = x->cos(prod(x))
        for I in CartesianIndices(poly)
            vals[I] = f(interpolation_nodes(poly,I))
        end
        xtest = ntuple(i->0.1,d)
        @test f(xtest) ≈ poly(xtest)
        # todo add test when xtest is one of the interpolation points
    end
end

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
#### Pressure model kernels

function perturbation_pressure(bl::AtmosModel{FT}, args, env, buoy) where {FT}
    dpdz = vuntuple(n_updrafts(turbconv_model(bl))) do i
        perturbation_pressure(bl, turbconv_model(bl).pressure, args, env, buoy, i)
    end
    return dpdz
end

"""
    perturbation_pressure(
        m::AtmosModel{FT},
        press::PressureModel,
        args,
        env,
        buoy,
        i,
    ) where {FT}

Returns the value of perturbation pressure gradient
for updraft i following He et al. (JAMES, 2020), given:

 - `m`, an `AtmosModel`
 - `press`, a `PressureModel`
 - `args`, top-level arguments
 - `env`, NamedTuple of environment variables
 - `buoy`, NamedTuple of environment and updraft buoyancies
 - `i`, index of the updraft
"""
function perturbation_pressure(
    m::AtmosModel{FT},
    press::PressureModel,
    args,
    env,
    buoy,
    i,
) where {FT}
    @unpack state, diffusive, aux = args
    # Alias convention:
    up = state.turbconv.updraft
    up_aux = aux.turbconv.updraft
    up_dif = diffusive.turbconv.updraft

    w_up_i = fix_void_up(up[i].ρa, up[i].ρaw / up[i].ρa)

    nh_press_buoy = press.α_b * buoy.up[i]
    nh_pressure_adv = -press.α_a * w_up_i * up_dif[i].∇w[3]
    # TO DO: Add updraft height dependency (non-local)
    nh_pressure_drag =
        press.α_d * (w_up_i - env.w) * abs(w_up_i - env.w) / press.H_up_min

    dpdz = nh_press_buoy + nh_pressure_adv + nh_pressure_drag

    return dpdz
end;


text:
#!/usr/bin/julia
#
# Definitions of parameters for different simulation methods and types.
export SimParams, SimMethod, SimVariation, SplitStep, Lanczos;

"""
Simulation parameters for simulation completing the time evolution of a system.

These parameters are common to all the simulation methods.
"""
struct SimParams{T<:Number}
  "Time step."
  Δt::T;

  "Total number of time steps."
  steps::UInt;
  "Error tolerance parameter."
  ϵ::Number;

  function SimParams{T}(Δt::T, steps::Real, ϵ::Number=1e-8) where T
    cΔt = ustrip(Δt);
    # Error check the passed arguments.
    if imag(Δt) != 0
      if (real(cΔt) < 0 || imag(cΔt) < 0)
        throw(DomainError(Δt, "timestep must be positive"));
      elseif (real(cΔt) > 0 && imag(cΔt) > 0)
        throw(DomainError(Δt, "timestep must either be real or imaginary"));
      elseif abs(cΔt) == 0
        throw(DomainError(Δt, "timestep must not be zero"));
      end
    elseif cΔt <= 0
      throw(DomainError(Δt, "timestep must be positive"));
    end

    steps <= 0 && throw(DomainError(steps, "number of timesteps must be positive"));
    steps % 1 !== 0 && throw(DomainError(steps, "number of timesteps must be an integer"));

    ustrip(ϵ) < 0 && throw(DomainError(ϵ, "error must be positive"));

    new(Δt, steps, ϵ);
  end
end
SimParams(Δt::T, steps::Real, ϵ::Number=1e-8) where T<:Number = SimParams{T}(Δt, steps, ϵ);

"""
Simulation methods.

Concrete instances of SimMethod should implement all of the parameters
specific to the simulation method.
"""
abstract type SimMethod end

"Split step simulation method."
struct SplitStep <: SimMethod end

"Short iterative Lanczos method."
struct Lanczos <: SimMethod
  "Number of basis functions."
  basisno::UInt;

  "Number of steps before re-generating basis."
  basissteps::UInt;
end

"""
Simulation variations

Modifications to the SimMethod types.

Any variation specific to a SimMethod should inherit from an abstract type
with name '<SimMethod Name>Variation'.

See also: [`SimMethod`](@ref), [`AdaptiveΔt`](@ref), [`LanczosVariation`](@ref)
"""
abstract type SimVariation end;


text:
module RchannelExperiments
using Reexport, Parameters
import Images
import Base.Dates: Time
import DSP, PyPlot
const P=PyPlot

export RcExp, RcRes, rho, g
const rho = 1000
const g = 9.81
const mu = 0.001792 # dynamic viscosity of water at 0C https://www.thermexcel.com/english/tables/eau_atm.htm
const nu = mu/rho  # kinematic viscosity

include("helpers.jl")
# The experiment setup, such as geometry, measurements, etc.
include("setup/RcSetup.jl")
@reexport using .RcSetup
# Image processing to get diameter, scalloping, etc
include("images/RcImages.jl")
@reexport using .RcImages
# All measurements which were recorded by the computer
include("labview/RcLabView.jl")
@reexport using .RcLabView
# volume
include("volume/RcVolume.jl")
@reexport using .RcVolume
# hand-recorded temperature
include("temperature/RcTemperature.jl")
@reexport using .RcTemperature

"""
This is one Experiment.
"""
@with_kw struct RcExp
    name::String # probably just the ISO date: 2017-10-06
    # start-time of experiment.  Plotting will be relative to this time
    experiment_start::DateTime
    # The time periods when water is flowing
    experiment_running::Vector{Tuple{Time,Time}}=Tuple{Time,Time}[]
    # Time of events
    events::Dict{Any,Any}=Dict{Any,Any}()
    setup::Setup
    # sub exps
    exi::Union{ExpImgs,Void}
    exl::Union{ExpLabView,Void}
    exv::ExpVol=ExpVol()
    ext::ExpTemp=ExpTemp(notebook_pdf="")
    extras::Dict{Symbol}=Dict{Symbol,Any}()
end

"Holds all results"
@with_kw struct RcRes
    ex::RcExp
    resi::Union{ExpImgsResults,Void}
    resl::Union{ExpLabViewResults,Void}
    rest::ExpTempResults=ExpTempResults()
    resv::ExpVolResults=ExpVolResults()
    extras::Dict{Symbol}=Dict{Symbol,Any}()
end

include("plotting.jl")

include("proc.jl")
include("uncertainty.jl")
# get relevant data
function get_Q(res::RcRes)
    return t, Q
end
function get_D(res::RcRes)
    return t, D
end
function get_p(res::RcRes)
    return t, p
end
function calc_f(res::RcRes)
end
function calc_heat_flux(res::RcRes)
end

end


text:
using TestParticle, Meshes, OrdinaryDiffEq, StaticArrays, Random
using Test

"Initial state perturbation for EnsembleProblem."
function prob_func(prob, i, repeat)
   remake(prob, u0=rand(MersenneTwister(i))*prob.u0)
end

@testset "TestParticle.jl" begin
   @testset "numerical field" begin
      x = range(-10, 10, length=15)
      y = range(-10, 10, length=20)
      z = range(-10, 10, length=25)
      B = fill(0.0, 3, length(x), length(y), length(z)) # [T]
      E = fill(0.0, 3, length(x), length(y), length(z)) # [V/m]

      B[3,:,:,:] .= 10e-9
      E[3,:,:,:] .= 5e-10

      Δx = x[2] - x[1]
      Δy = y[2] - y[1]
      Δz = z[2] - z[1]

      mesh = CartesianGrid((length(x)-1, length(y)-1, length(z)-1),
         (x[1], y[1], z[1]),
         (Δx, Δy, Δz))

      x0 = [0.0, 0.0, 0.0] # initial position, [m]
      u0 = [1.0, 0.0, 0.0] # initial velocity, [m/s]
      stateinit = [x0..., u0...]

      param = prepare(mesh, E, B)
      tspan = (0.0, 1.0)
      trace! = trace_numeric!

      prob = ODEProblem(trace!, stateinit, tspan, param)

      sol = solve(prob, Tsit5(); save_idxs=[1,2,3])

      x = getindex.(sol.u, 1)
      y = getindex.(sol.u, 2)
      z = getindex.(sol.u, 3)

      @test length(x) == 8 && x[end] ≈ 0.8540967226885379

      trajectories = 10
      prob = ODEProblem(trace!, stateinit, tspan, param)
      ensemble_prob = EnsembleProblem(prob, prob_func=prob_func)
      sol = solve(ensemble_prob, Tsit5(), EnsembleThreads();
         trajectories=trajectories, save_idxs=[1,2,3])

      x = getindex.(sol.u[10].u, 1)
      y = getindex.(sol.u[10].u, 2)
      z = getindex.(sol.u[10].u, 3)

      @test x[7] ≈ 0.09615629718624641 rtol=1e-6

      stateinit = SA[x0..., u0...]

      param = prepare(mesh, E, B)
      tspan = (0.0, 1.0)
      trace = trace_numeric

      prob = ODEProblem(trace, stateinit, tspan, param)

      sol = solve(prob, Tsit5(); save_idxs=[1,2,3])

      x = getindex.(sol.u, 1)
      y = getindex.(sol.u, 2)
      z = getindex.(sol.u, 3)

      @test length(x) == 8 && x[end] ≈ 0.8540967226885379
   end

   @testset "analytical field" begin
      # in-place version
      Ek = 5e7 # [eV]

      m = TestParticle.mᵢ
      q = TestParticle.qᵢ
      c = TestParticle.c
      Rₑ = TestParticle.Rₑ

      # initial velocity, [m/s]
      v₀ = TestParticle.sph2cart(c*sqrt(1-1/(1+Ek*q/(m*c^2))^2), 0.0, π/4)
      # initial position, [m]
      r₀ = TestParticle.sph2cart(2.5*Rₑ, 0.0, π/2)
      stateinit = [r₀..., v₀...]

      param = prepare(TestParticle.getE_dipole, TestParticle.getB_dipole)
      tspan = (0.0, 1.0)

      trace! = trace_analytic!
      prob = ODEProblem(trace!, stateinit, tspan, param)

      sol = solve(prob, Tsit5(); save_idxs=[1,2,3])

      x = getindex.(sol.u, 1)

      @test x[300] ≈ 1.2592311352654776e7 rtol=1e-6

      # static array version (results not identical with above: maybe some bugs?)
      stateinit = SA[r₀..., v₀...]

      trace = trace_analytic
      prob = ODEProblem(trace, stateinit, tspan, param)

      sol = solve(prob, Tsit5(); save_idxs=[1,2,3])

      x = getindex.(sol.u, 1)

      @test x[306] ≈ 1.2588844644203672e7 rtol=1e-6

   end
end


text:
function white_noise(N::Int, pix_size::Float64, wnl::Float64)

    rng = MersenneTwister()
    white_noise_map = randn(rng, N, N) .* (wnl / pix_size)

    return white_noise_map

end

function atmospheric_noise(N::Int, pix_size::Float64, anl::Float64)

    rng = MersenneTwister()
    random_realization = randn(rng, N, N)

    inds = ((0:(N-1)) .+ 0.5 .- (N/2))
    X = reshape(repeat(inds, N), N ,N)
    Y = X'
    R = ((X .^2 .+ Y .^2).^0.5) .* (pix_size / 60.0)

    mag_k = ((2.0 * π) ./ (R .+ 0.01 )).^(5.0 / 3.0)
    atm_rand = fft(random_realization)
    atm_map = ifft(atm_rand .* fftshift(mag_k))
    atm_map_norm = atm_map .* (anl / pix_size)

    return real.(atm_map_norm)

end

function one_over_f(N::Int, pix_size::Float64, oofnl::Float64)

    rng = MersenneTwister()
    random_realization = randn(rng, N, N)

    inds = ((0:(N-1)) .+ 0.5 .- (N/2))
    X = reshape(repeat(inds, N), N ,N)' .* (pix_size / 60.0)
    kx = (2.0 * π) ./ (X .+ 0.01)
    rand_ft = fft(random_realization)
    oofn_map = ifft(rand_ft .* fftshift(kx))
    oofn_map_norm = oofn_map .* (oofnl / pix_size)

    return real.(oofn_map_norm)

end

function one_over_f2(N::Int, pix_size::Float64, oofnl::Float64)

    rng = MersenneTwister()
    random_realization = randn(rng, N, N)

    inds = ((0:(N-1)) .+ 0.5 .- (N/2))
    X = reshape(repeat(inds, N), N ,N)' .* (pix_size / 60.0)
    kx = ((2.0 * π) ./ (X .+ 0.01)).^(3)
    rand_ft = fft(random_realization)
    oofn_map = ifft(rand_ft .* fftshift(kx))
    oofn_map_norm = oofn_map .* (oofnl / pix_size)

    return real.(oofn_map_norm)

end


text:
# using Printf

println("... test threading 2 ...")

# println("... test channel 4 example parallel ...")
# A = []
# ch1 = Channel(32)
# ch2 = Channel(32)
#
# function simple1(c::Clock, input::Channel, output::Channel, name, id, op)
#     token = take!(input)         # take something from the input
#     push!(A, (tau(c), name, id, token))
# #    println(@sprintf("%5.2f: %s %d took %d", tau(c), name, id, token))
#     d = delay!(c, rand())           # after a delay
# #    println(@sprintf("%5.2f: %s %d now after a delay ...", tau(c), name, id))
#     put!(output, op(token, id))  # put it out with some op applied
# #    println(@sprintf("%5.2f: %s %d sent %d", tau(c), name, id, op(token, id)))
# end
#
# clk = Clock(0.001)
# Random.seed!(123)
#
# for i in 1:2:8    # create, register and start 8 SimProcesses SP
#     process!(clk, SP(i, simple1, clk, ch1, ch2, "foo", i, +), spawn=true)
#     process!(clk, SP(i+1, simple1, clk, ch2, ch1, "bar", i+1, *), spawn=true)
# end
#
# @test length(clk.processes) == 8
# for p in values(clk.processes)
#     @test p.state == Simulate.Idle()
#     @test istaskstarted(p.task)
# end
#
# put!(ch1, 1)
# sleep(0.1)
# run!(clk, 10)
#
# # @test length(A) > 20
# # p = [i[3] for i in A]
# # for i in 1:8
# #     @test i ∈ p  # all processes did something
# # end


text:
using ForwardDiff, Roots, Plots, LinearAlgebra, Optim,Statistics

function strong_backtracking(f, ∇, x, d; α=1, β=0.1, σ=0.4)
    y0, g0, y_prev, α_prev = f(x), ∇(x)⋅d, NaN, 0
    αlo, αhi = NaN, NaN
    # bracket phase
    while true
        y = f(x + α*d)
        if y > y0 + β*α*g0 || (!isnan(y_prev) && y ≥ y_prev)
            αlo, αhi = α_prev, α
            break
        end
        g = ∇(x + α*d)⋅d
        if abs(g) ≤ -σ*g0
            return α
        elseif g ≥ 0
            αlo, αhi = α, α_prev
            break
        end
        y_prev, α_prev, α = y, α, 2α
    end

    # zoom phase
    ylo = f(x + αlo*d)
    while true
        α = (αlo + αhi)/2
        y = f(x + α*d)
        if y > y0 + β*α*g0 || y ≥ ylo
            αhi = α
        else
            g = ∇(x + α*d)⋅d
            if abs(g) ≤ -σ*g0
                return α
            elseif g*(αhi - αlo) ≥ 0
                αhi = αlo
            end
            αlo = α
        end
    end
end

function bfgs(f,grad_f,x0;tol=1e-6,maxitr=10000)

    #init data
    ite = 0;x=x0;conv=false;err=0;
    xh = zeros(length(x0),maxitr+1);xh[:,1] = x0;

    #init bfgs parameters
    m = length(x)
    D = Matrix(1.0I,m,m)
    g = grad_f(x)

    #iteractions
    while ite < maxitr
        #line search
        d = -D*g

        #Inexact line search (strongWolf)
        alpha = strong_backtracking(f,grad_f,x,d;α=1, β=1e-4, σ=0.1)

        #update x
        s = alpha*d;
        x += s

        gl = grad_f(x)
        y = gl - g
        D = D - (s*y'*D + D*y*s')/(s'*y) + (1+(y'*D*y)/(s'*y))[1]*(s*s')/(s'*y)

        #update bfgs parameters
        g = gl

        #update ite and x history
        ite += 1
        xh[:,ite+1] = x

        #local minimum point
        err = norm(s);
        if err < tol ; conv=true; break; end
    end

    println(ite+1)
    #return x history. convergence and error
    return xh[:,1:ite+1], conv, err
end

#valores dos parametros
ka = 9.0
kb = 2.0
La = 10.0
Lb = 10.0

#Caso 1
#F1 = 2.0; F2 = 4.0;

#Caso 2
#F1 = 0.0 ;F2 = 4.0; #nao ta convergindo

#Caso 3
F1 = 0.0; F2 = 400.0;

#function
F(u) = 0.5e0 * ka * (sqrt(u[1] ^ 2 + (La - u[2]) ^ 2) - La) ^ 2 + 0.5e0 * kb * (sqrt(u[1] ^ 2 + (Lb + u[2]) ^ 2) - Lb) ^ 2 - F1 * u[1] - F2 * u[2];
grad_F(u) = ForwardDiff.gradient(F,u)

#tolerancia
tol = 1e-8

#numero maximo de iteracoes
iterac_max = 10000

#ponto inicial
#u0 = [-0.01,0.0]
u0 = [ 0.00,0.0]

#u*
@time u_est = bfgs(F,grad_F,u0;tol=tol,maxitr=iterac_max)
u_est = u_est[1][:,end]
println("u*: ",u_est)
println("f(u*): ",F(u_est))

#gradiente no ponto de minimo
println("∇f(u*): ",grad_F(u_est))

#verificar a hessiana no ponto de minimo
hess_F(u) = ForwardDiff.hessian(F,u)
println("H(u*): ",hess_F(u_est))
println("λ: ",eigvals(hess_F(u_est)))


text:
using Test
using LinearAlgebra

using problems.matrix

@testset "matrix rotation 1" begin

    @testset "1x1" begin
        mat = Int8[1;]
        actual = rotate_matrix1(mat)
        expected = Int8[1;]

        @test actual == expected
    end

    @testset "2x2" begin
        mat = Int8[1 2;
                   3 4]

        actual = rotate_matrix1(mat)

        expected = Int8[3 1;
                        4 2]

        @test actual == expected
    end

    @testset "3x3" begin
        mat = Int8[1 2 3;
                   4 5 6;
                   7 8 9]

        actual = rotate_matrix1(mat)

        expected = Int8[7 4 1;
                        8 5 2;
                        9 6 3]

        @test actual == expected
    end

    @testset "4x4" begin

        mat = Int8[1 2 3 4;
                   5 6 7 8;
                   9 10 11 12;
                   13 14 15 16]

        actual = rotate_matrix1(mat)

        expected = Int8[13 9 5 1;
                        14 10 6 2;
                        15 11 7 3;
                        16 12 8 4]

        @test actual == expected
    end

    @testset "5x5" begin

        mat = Int8[1 2 3 4 5;
                   6 7 8 9 10;
                   11 12 13 14 15;
                   16 17 18 19 20;
                   21 22 23 24 25]

        actual = rotate_matrix1(mat)

        expected = Int8[21 16 11 6 1;
                        22 17 12 7 2;
                        23 18 13 8 3;
                        24 19 14 9 4;
                        25 20 15 10 5]

        @test actual == expected
    end
end

@testset "matrix rotation 2" begin

    @testset "1x1" begin
        mat = Int8[1;]
        actual = rotate_matrix2(mat)
        expected = Int8[1;]

        @test actual == expected
    end

    @testset "2x2" begin
        mat = Int8[1 2;
                   3 4]

        actual = rotate_matrix2(mat)

        expected = Int8[3 1;
                        4 2]

        @test actual == expected
    end

    @testset "3x3" begin
        mat = Int8[1 2 3;
                   4 5 6;
                   7 8 9]

        actual = rotate_matrix2(mat)

        expected = Int8[7 4 1;
                        8 5 2;
                        9 6 3]

        @test actual == expected
    end

    @testset "4x4" begin

        mat = Int8[1 2 3 4;
                   5 6 7 8;
                   9 10 11 12;
                   13 14 15 16]

        actual = rotate_matrix2(mat)

        expected = Int8[13 9 5 1;
                        14 10 6 2;
                        15 11 7 3;
                        16 12 8 4]

        @test actual == expected
    end

    @testset "5x5" begin

        mat = Int8[1 2 3 4 5;
                   6 7 8 9 10;
                   11 12 13 14 15;
                   16 17 18 19 20;
                   21 22 23 24 25]

        actual = rotate_matrix2(mat)

        expected = Int8[21 16 11 6 1;
                        22 17 12 7 2;
                        23 18 13 8 3;
                        24 19 14 9 4;
                        25 20 15 10 5]

        @test actual == expected
    end

    @testset "zero matrix 1" begin

        @testset "2x2" begin
            mat = Int8[1 2;
                       0 4]

            actual = zero_matrix1(mat)

            expected = Int8[0 2;
                            0 0]

            @test actual == expected
        end

    end

    @testset "3x3" begin
        mat = Int8[1 2 3;
                   4 0 6;
                   7 8 9]

        actual = zero_matrix1(mat)

        expected = Int8[1 0 3;
                        0 0 0;
                        7 0 9]

        @test actual == expected
    end

    @testset "4x4" begin

        mat = Int8[0 2 3 4;
                   5 6 7 8;
                   9 10 0 0;
                   13 14 15 16]

        actual = zero_matrix1(mat)

        expected = Int8[0 0 0 0;
                        0 6 0 0;
                        0 0 0 0;
                        0 14 0 0]

        @test actual == expected
    end

    @testset "5x5" begin

        mat = Int8[1 2 3 4 5;
                   6 7 8 9 10;
                   0 12 0 14 15;
                   16 17 18 19 20;
                   0 22 0 24 25]

        actual = zero_matrix1(mat)

        expected = Int8[0 2 0 4 5;
                        0 7 0 9 10;
                        0 0 0 0 0;
                        0 17 0 19 20;
                        0 0 0 0 0]

        @test actual == expected
    end

end

max_stars_repo_name: bmoretz/Daily-Coding-Problem

max_issues_repo_name: bmoretz/Daily-Coding-Problem

max_forks_repo_name: bmoretz/Daily-Coding-Problem


text:
export knn

import ..KNNModel

using mlpack._Internal.io

import mlpack_jll
const knnLibrary = mlpack_jll.libmlpack_julia_knn

# Call the C binding of the mlpack knn binding.
function knn_mlpackMain()
  success = ccall((:knn, knnLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module knn_internal
  import ..knnLibrary

import ...KNNModel

# Get the value of a model pointer parameter of type KNNModel.
function IOGetParamKNNModel(paramName::String)::KNNModel
  KNNModel(ccall((:IO_GetParamKNNModelPtr, knnLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type KNNModel.
function IOSetParamKNNModel(paramName::String, model::KNNModel)
  ccall((:IO_SetParamKNNModelPtr, knnLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializeKNNModel(stream::IO, model::KNNModel)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializeKNNModelPtr, knnLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializeKNNModel(stream::IO)::KNNModel
  buffer = read(stream)
  KNNModel(ccall((:DeserializeKNNModelPtr, knnLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    knn(; [algorithm, epsilon, input_model, k, leaf_size, query, random_basis, reference, rho, seed, tau, tree_type, true_distances, true_neighbors, verbose])

This program will calculate the k-nearest-neighbors of a set of points using
kd-trees or cover trees (cover tree support is experimental and may be slow).
You may specify a separate set of reference points and query points, or just a
reference set which will be used as both the reference and query set.

For example, the following command will calculate the 5 nearest neighbors of
each point in `input` and store the distances in `distances` and the neighbors
in `neighbors`:

```julia
julia> using CSV
julia> input = CSV.read("input.csv")
julia> distances, neighbors, _ = knn(k=5, reference=input)
```

The output is organized such that row i and column j in the neighbors output
matrix corresponds to the index of the point in the reference set which is the
j'th nearest neighbor from the point in the query set with index i.  Row j and
column i in the distances output matrix corresponds to the distance between
those two points.

# Arguments

 - `algorithm::String`: Type of neighbor search: 'naive', 'single_tree',
      'dual_tree', 'greedy'.  Default value `dual_tree`.

 - `epsilon::Float64`: If specified, will do approximate nearest neighbor
      search with given relative error.  Default value `0`.

 - `input_model::KNNModel`: Pre-trained kNN model.
 - `k::Int`: Number of nearest neighbors to find.  Default value `0`.

 - `leaf_size::Int`: Leaf size for tree building (used for kd-trees, vp
      trees, random projection trees, UB trees, R trees, R* trees, X trees,
      Hilbert R trees, R+ trees, R++ trees, spill trees, and octrees).  Default
      value `20`.

 - `query::Array{Float64, 2}`: Matrix containing query points (optional).
 - `random_basis::Bool`: Before tree-building, project the data onto a
      random orthogonal basis.  Default value `false`.

 - `reference::Array{Float64, 2}`: Matrix containing the reference
      dataset.
 - `rho::Float64`: Balance threshold (only valid for spill trees).
      Default value `0.7`.

 - `seed::Int`: Random seed (if 0, std::time(NULL) is used).  Default
      value `0`.

 - `tau::Float64`: Overlapping size (only valid for spill trees).  Default
      value `0`.

 - `tree_type::String`: Type of tree to use: 'kd', 'vp', 'rp', 'max-rp',
      'ub', 'cover', 'r', 'r-star', 'x', 'ball', 'hilbert-r', 'r-plus',
      'r-plus-plus', 'spill', 'oct'.  Default value `kd`.

 - `true_distances::Array{Float64, 2}`: Matrix of true distances to
      compute the effective error (average relative error) (it is printed when
      -v is specified).
 - `true_neighbors::Array{Int, 2}`: Matrix of true neighbors to compute
      the recall (it is printed when -v is specified).
 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `distances::Array{Float64, 2}`: Matrix to output distances into.
 - `neighbors::Array{Int, 2}`: Matrix to output neighbors into.
 - `output_model::KNNModel`: If specified, the kNN model will be output
      here.

"""
function knn(;
             algorithm::Union{String, Missing} = missing,
             epsilon::Union{Float64, Missing} = missing,
             input_model::Union{KNNModel, Missing} = missing,
             k::Union{Int, Missing} = missing,
             leaf_size::Union{Int, Missing} = missing,
             query = missing,
             random_basis::Union{Bool, Missing} = missing,
             reference = missing,
             rho::Union{Float64, Missing} = missing,
             seed::Union{Int, Missing} = missing,
             tau::Union{Float64, Missing} = missing,
             tree_type::Union{String, Missing} = missing,
             true_distances = missing,
             true_neighbors = missing,
             verbose::Union{Bool, Missing} = missing,
             points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, knnLibrary), Nothing, ());

  IORestoreSettings("k-Nearest-Neighbors Search")

  # Process each input argument before calling mlpackMain().
  if !ismissing(algorithm)
    IOSetParam("algorithm", convert(String, algorithm))
  end
  if !ismissing(epsilon)
    IOSetParam("epsilon", convert(Float64, epsilon))
  end
  if !ismissing(input_model)
    knn_internal.IOSetParamKNNModel("input_model", convert(KNNModel, input_model))
  end
  if !ismissing(k)
    IOSetParam("k", convert(Int, k))
  end
  if !ismissing(leaf_size)
    IOSetParam("leaf_size", convert(Int, leaf_size))
  end
  if !ismissing(query)
    IOSetParamMat("query", query, points_are_rows)
  end
  if !ismissing(random_basis)
    IOSetParam("random_basis", convert(Bool, random_basis))
  end
  if !ismissing(reference)
    IOSetParamMat("reference", reference, points_are_rows)
  end
  if !ismissing(rho)
    IOSetParam("rho", convert(Float64, rho))
  end
  if !ismissing(seed)
    IOSetParam("seed", convert(Int, seed))
  end
  if !ismissing(tau)
    IOSetParam("tau", convert(Float64, tau))
  end
  if !ismissing(tree_type)
    IOSetParam("tree_type", convert(String, tree_type))
  end
  if !ismissing(true_distances)
    IOSetParamMat("true_distances", true_distances, points_are_rows)
  end
  if !ismissing(true_neighbors)
    IOSetParamUMat("true_neighbors", true_neighbors, points_are_rows)
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("distances")
  IOSetPassed("neighbors")
  IOSetPassed("output_model")
  # Call the program.
  knn_mlpackMain()

  return IOGetParamMat("distances", points_are_rows),
         IOGetParamUMat("neighbors", points_are_rows),
         knn_internal.IOGetParamKNNModel("output_model")
end


text:
struct PlsrAgg
    fm::Plsr
    nlv
    wagg
    w
end

"""
    plsr_agg(X, Y, weights = ones(size(X, 1)); nlv, wagg = "unif")
Aggregation of PLSR models with different numbers of LVs.
* `X` : X-data.
* `Y` : Y-data. Must be univariate if `wagg` != "unif".
* weights : Weights of the observations.
* `nlv` : A character string such as "5:20" defining the range of the numbers of LVs
    to consider ("5:20": the predictions of models with nb LVS = 5, 6, ..., 20
    are averaged). Syntax such as "10" is also allowed ("10": correponds to
    the single model with 10 LVs).
* `wagg` : Type of averaging.

Ensemblist method where the predictions are calculated by "averaging"
the predictions of a set of models built with different numbers of
latent variables (LVs).

For instance, if argument `nlv` is set to `nlv = "5:10"`, the prediction for
a new observation is the average (eventually weighted) of the predictions
returned by the models with 5 LVS, 6 LVs, ... 10 LVs, respectively.

Dependending argument `wagg`, the average is the simple mean (`wagg` = "unif") or
a weighted mean using AIC weights computed from the models (see function `aicplsr`):
* "aic" : Usual AIC weights (exp(-delta/2)).
* "sqrt" : Sqrt(delta) is used in place of delta in AIC weights.
* "inv" : Weights are computed by 1 / AIC.
* "fair" : A decreasing "fair" function is applied to delta.
* "shenk" : Shenk et al. (1997) weights. See function `wshenk'.
"""
function plsr_agg(X, Y, weights = ones(size(X, 1)); nlv, wagg = "unif")
    plsr_agg!(copy(X), copy(Y), weights; nlv = nlv, wagg = wagg)
end

function plsr_agg!(X, Y, weights = ones(size(X, 1)); nlv, wagg = "unif")
    n = size(X, 1)
    p = size(X, 2)
    nlv = eval(Meta.parse(nlv))
    nlvmax = maximum(nlv)
    nlv = (max(minimum(nlv), 0):min(nlvmax, n, p))
    if isequal(wagg, "unif")
        w = ones(nlvmax + 1)
    elseif isequal(wagg, "aic")
        w = aicplsr(X, Y; nlv = nlvmax).w.aic
    elseif isequal(wagg, "sqrt")
        d = aicplsr(X, Y; nlv = nlvmax).delta.aic
        w = exp.(-sqrt.(d) / 2)
    elseif isequal(wagg, "inv")
        w = 1 ./ aicplsr(X, Y; nlv = nlvmax).crit.aic
    elseif isequal(wagg, "fair")
        d = aicplsr(X, Y; nlv = nlvmax).delta.aic
        d = d / maximum(d[isnan.(d) .== 0])
        w = 1 ./ (1 .+ d).^2
    elseif isequal(wagg, "shenk")
        fm = plskern(X, Y, weights; nlv = (nlvmax))
        w = vec(sum(Jchemo.wshenk(fm, X).W, dims = 1))
        w = [0 ; w]
    end
    w[isnan.(w)] .= 0
    #w = vec(mavg(w', 3))
    w = w[collect(nlv) .+ 1]
    w ./= sum(w)
    fm = plskern!(X, Y, weights; nlv = nlvmax)
    PlsrAgg(fm, nlv, wagg, w)
end

"""
    predict(object::PlsrAgg, X)
Compute Y-predictions from a fitted model.
* `object` : The fitted model.
* `X` : X-data for which predictions are computed.
"""
function predict(object::PlsrAgg, X)
    nlv = object.nlv
    le_nlv = length(nlv)
    zpred = predict(object.fm, X; nlv = nlv).pred
    if(le_nlv == 1)
        pred = zpred
    else
        acc = object.w[1] * copy(zpred[1])
        @inbounds for i = 2:le_nlv
            acc .+= object.w[i] * zpred[i]
        end
        pred = acc
    end
    (pred = pred, predlv = zpred)
end


text:
Base.@kwdef mutable struct EvaluateDQNAgent{M, AP<:AbstractValuePolicy, Φ, IR<:HistImageReplay} <: AbstractAgent
    model::M
    ap::AP
    replay::IR
    prev_s::Φ
    hist::Int
    action::Int = 0
    INFO::Dict{Symbol, Any} = Dict{Symbol, Any}()
end

EvaluateDQNAgent(model, ϵ, hist, action_set, img_manip, img_norm, img_size) =
    EvaluateDQNAgent(
        model = model,
        ap = ϵGreedy(0.1, action_set),
        replay = HistImageReplay(1, img_size, img_manip, img_norm, hist, 1),
        prev_s = zeros(Int, hist),
        hist = hist
    )

function RLCore.start!(agent::EvaluateDQNAgent,
                       env_s_tp1,
                       rng::AbstractRNG;
                       kwargs...)

    agent.prev_s .= add!(agent.replay, env_s_tp1)
    state = agent.replay.img_norm(
        reshape(getindex(agent.replay.image_buffer, agent.prev_s),
                agent.replay.image_buffer.img_size..., agent.replay.hist,
                1))
    agent.action = sample(agent.ap,
                          agent.model(state),
                          rng)#::eltype(agent.ap)
    return agent.action
end

function RLCore.step!(agent::EvaluateDQNAgent,
                      env_s_tp1,
                      r,
                      terminal,
                      rng::AbstractRNG;
                      kwargs...)

    add_ret = add!(agent.replay,
                   (agent.prev_s,
                    findfirst((a)->a==agent.action,
                              agent.ap.action_set)::Int,
                    env_s_tp1,
                    r,
                    terminal))
    agent.prev_s .= add_ret
    prev_s =
        agent.replay.img_norm(
            reshape(getindex(agent.replay.image_buffer, agent.prev_s),
                    agent.replay.image_buffer.img_size..., agent.replay.hist,
                    1))

    agent.action = sample(agent.ap,
                          agent.model(prev_s),
                          rng)
    return agent.action
end


text:
"""
An evaluator specialised to the travelling salesman problem.
"""
type TSPEvaluator <: Evaluator
  stage::AbstractString
  cities::Int
  threads::Int
  distance::Array{Float, 2}

  TSPEvaluator(st::AbstractString, n::Int, t::Int, d::Array{Float, 2}) =
    new(st, n, t, d)
end

"""
Contains a definition of a TSP evaluator.
"""
type TSPEvaluatorDefinition <: EvaluatorDefinition
  stage::AbstractString
  file::AbstractString
  threads::Int

  TSPEvaluatorDefinition() = new("", "", 1)
  TSPEvaluatorDefinition(s::AbstractString, f::AbstractString, t::Int) =
    new(s, f, t)
end

"""
Composes a TSP evaluator from its definition.
"""
function compose!(ev::TSPEvaluatorDefinition, p::Population)
  # Find the name of the stage that is used for evaluation if none is provided.
  if ev.stage == ""
    ev.stage = p.demes[1].species.genotype.label
  end

  # Load the contents of the file and convert into a list of co-ordinates.
  coords = open(ev.file, "r") do f
    c = Vector{Float}[]
    for ln in eachline(f)
      !isempty(ln) && push!(c, Float[float(x) for x in split(ln)])
    end
    c
  end
  num = length(coords)

  # Generate a distance matrix from the set of co-ordinates.
  matrix = Array(Float, num, num)
  for i = 1:num
    for j = 1:num
      matrix[i, j] = sqrt(sum((coords[i] - coords[j]) .^ 2))
    end
  end

  # Build and return the evaluator.
  return TSPEvaluator(ev.stage, num, ev.threads, matrix)
end

"""
Evaluates the fitness of a TSP tour for a pre-determined geometric problem.

**Positional Arguments**:

* `file::AbstractString`, the path to the file containing the co-ordinates of
  each of the cities for the geometric TSP instance being solved, relative to
  the current working directory.

**Keyword Arguments**:

* `threads::Int`, the number of threads that the evaluation process should be
  split across. Defaults to `1`.
"""
tsp(file::AbstractString; stage::AbstractString = "", threads = 1) =
  TSPEvaluatorDefinition(stage, file, max(1, threads))

"""
Computes the length of a provided tour and transforms it into a fitness value
according to a provided fitness scheme.
"""
function evaluate!(e::TSPEvaluator, sch::FitnessScheme, tour::Vector{Int})
  length = zero(Float)
  for i in 1:(e.cities - 1)
    length += e.distance[tour[i], tour[i+1]]
  end
  length += e.distance[tour[end], tour[1]]
  assign(sch, length)
end


text:
@testset "RCF" begin

  @testset "Various examples" begin

    l=collect(Hecke.quadratic_extensions(50))
    @test length(l)==30

    fields=Hecke.D5_extensions(3*fmpz(10)^3, l)
    @test length(fields)==1

    fields=Hecke.Dn_extensions(3, 2*fmpz(10)^4, l)
    @test length(fields)==1

    fields=Hecke.C3xD5_extensions(fmpz(2000000))
    @test length(fields)==0

    fields=Hecke.S3xC5_extensions(fmpz(353517))
    @test length(fields)==0

    fields1=Hecke.C9semiC4(fmpz(10)^30)
    @test fields1==1

  end

  @testset "abelian extensions && related examples" begin

    Qx,x=PolynomialRing(FlintQQ, "x")
    K,a=NumberField(x-1,"a")
    O=maximal_order(K)
    l=Hecke.abelian_normal_extensions(O,[2,2], fmpz(10)^4)
    @test length(l)==47
    l1=collect(Hecke.C22_extensions(10^4))
    @test length(l1)==47

    K,a=NumberField(x^2+1, "a")
    O=maximal_order(K)
    l=Hecke.abelian_normal_extensions(O,[2], fmpz(10)^5, with_autos=Val{true})
    @test length(l)==41
    for x in l[1:5]
      K, autos=Hecke._from_relative_to_abs(x)
      @test length(autos)==2
      y=small_generating_set(closure(autos, *))
      @test length(y)==1 || length(y)==2
      if length(y)==2
        @test y[1]*y[2]==y[2]*y[1]
      else
        @test y[1]^2!=Hecke.NfToNfMor(K,K,gen(K))
      end
    end
  end

end


text:
module ANOVA

using GLM, DataFrames, Distributions, CategoricalArrays, ArgCheck, StatsModels, LinearAlgebra
import StatsBase

struct AnovaObject
  Source::Array{String,1}
  DF::Array{AbstractFloat,1}
  SS::Array{AbstractFloat,1}
  MSS::Array{AbstractFloat,1}
  F::Array{AbstractFloat,1}
  p::Array{AbstractFloat,1}
end

struct AnovaDataFrameRegressionModel{M} <: StatsModels.RegressionModel
  model::M
  mf::ModelFrame
  mm::ModelMatrix
  anova::AnovaObject
end

struct AnovaOptions
  anovatype::Int
end

## helper functions ##

function merge_bool_array(a, b)
  return ifelse.(b,b,a)
end

# calculate effects for type I ANOVA
function effects(mod)
    return (mod.pp.X / cholesky!(mod.pp).U)' * mod.rr.y
end

# this will drop variables from the matrix using the mask argument then fit a new linear model
function droptermbymask(mod, mask)
  matrix = mod.pp.X[:,.!mask]
  y = mod.rr.y
  return lm(matrix,y)
end

# calculate ANOVA

"""
    anova(mod, mf, mm[, anovatype = 3])

Compute an ANOVA table for the given linear model.
"""
function anova(mod::LinearModel, mf::ModelFrame, mm::ModelMatrix; anovatype = 2)

  eff = effects(mod) # get effects

  T = eltype(mod.pp.X) # get types
  response = mf.terms.eterms[1]
  terms = mf.terms.terms

  ## calculate variables for residuals (this is the full modell) ##
  DFres = dof_residual(mod)           # degrees of freedom
  #SSres = sum(abs2.(residuals(mod))) # used deviance() instead
  SSres = deviance(mod)               # calculate sum of squares of weighted residuals
  MSSres = SSres / DFres              # calculate mean sum of squares for residuals

  # check if there is an intercept
  v = all(mod.pp.X[:,1] .== 1.0) ? 1 : 0

  # get assigned numbers for factors
  factors_assig = unique(mm.assign)

  ## calculate some variables for computation ##
  n = length(predict(mod))              # sample size
  Nfactors = length(factors_assig) - v  # number of factors

  ## create arrays for ANOVA factors ##
  DF = zeros(T,Nfactors)   # will contain degrees of freedoms for factor
  SS = zeros(T,Nfactors)   # will contain sum of squares for factor
  MSS = zeros(T,Nfactors)  # will contain mean sum of squares for factor
  RSS = zeros(T,Nfactors)  #
  F = zeros(T,Nfactors)    # will contain F values for factors
  p = zeros(T,Nfactors)    # will contain pvalue for factors

for i in 1:Nfactors
  mask = mm.assign .== factors_assig[i+v]
  if anovatype == 3 # if this is a type III ANOVA, we use a similar procedure to R's drop1() function
    fullmod = mod
    newmod = droptermbymask(fullmod, mask)
    RSS[i] = deviance(newmod)
    SS[i] = RSS[i] - deviance(fullmod)
    DF[i] = sum(mask)
    MSS[i] = SS[i]/DF[i]
    F[i] = ftest(fullmod,newmod).fstat[1]
  elseif anovatype == 2 # if this is a type II ANOVA
    full_model_mask = falses(length(mm.assign))
    for j in 1:Nfactors
      term = terms[j]
      if(typeof(term) == Expr && typeof(terms[i])==Symbol) # this is not an interaction term
        if(terms[i] in term.args)
          full_model_mask = merge_bool_array(full_model_mask, mm.assign .== factors_assig[j+v])
        end
       elseif(typeof(term) == Expr && typeof(terms[i])==Expr) # this is an interaction term
        if(all(in.(terms[i].args, (term.args,) )) && term != terms[i])
          full_model_mask = merge_bool_array(full_model_mask, mm.assign .== factors_assig[j+v])
        end
      end
    end
    fullmod = droptermbymask(mod, full_model_mask) # we have to drop all interaction terms containing the current term in the full modell
    newmod = droptermbymask(mod, merge_bool_array(full_model_mask, mask))
    RSS[i] = deviance(newmod)
    SS[i] = RSS[i] - deviance(fullmod)
    DF[i] = sum(mask)
    MSS[i] = SS[i]/DF[i]
    F[i] = MSS[i]/MSSres
  else # type I ANOVA
    SS[i] = sum(abs2.(eff[mask]))
    RSS[i] = 0
    DF[i] = sum(mask)
    MSS[i] = SS[i]/DF[i]
    F[i] = MSS[i]/MSSres
  end

  p[i] = ccdf(FDist(DF[i], DFres), F[i])
end

AnovaObject(
          vcat(string.(terms),"Residuals"),
          vcat(DF,DFres),
          vcat(SS, SSres),
          vcat(MSS, MSSres),
          vcat(F, 0),
          vcat(p, 0)
          )
end

anova(mod::StatsModels.DataFrameRegressionModel; anovatype = 3) = anova(mod.model, mod.mf, mod.mm, anovatype = anovatype)

function Base.isapprox(a::ANOVA.AnovaObject, b::ANOVA.AnovaObject)
  all( [a.Source == b.Source,
  all(a.DF .≈ b.DF),
  all(a.SS .≈ b.SS),
  all(a.MSS .≈ b.MSS),
  all(a.F .≈ b.F),
  all( a.p .≈ b.p)])
end

function Base.show(io::IO, m::AnovaDataFrameRegressionModel)
  println("ANOVA table for linear model")
  print(m.anova)
end

function Base.show(io::IO, x::AnovaObject)
  print(DataFrame(
          Source = x.Source,
          DF = x.DF,
          SS = x.SS,
          MSS = x.MSS,
          F = x.F,
          p = x.p
          ))
end

export AnovaDataFrameRegressionModel, AnovaOptions, anova, AnovaObject

end # module


text:
using Pkg3Dev
using Base.Test
using Pkg3

@testset "Pkg3Dev" begin
    # Make a depot
    depotpath = tempname()
    Pkg3Dev.create_depot(depotpath)
    @test isdir(depotpath)
    @test isdir(joinpath(depotpath, "registries"))
    @test Pkg3Dev.DEFAULT_DEPOT[] == depotpath
    @test depotpath in Pkg3.depots()

    # Make a registry
    Pkg3Dev.create_registry("TestRegistry", depot = depotpath)
    @test isdir(joinpath(depotpath, "registries", "TestRegistry"))
    @test isfile(joinpath(depotpath, "registries", "TestRegistry", "registry.toml"))
    @test Pkg3Dev.registry_depot("TestRegistry") == depotpath

    # TODO Make a new, registered package
    Pkg3Dev.create_package("TestPackage", "TestRegistry"; license = "MIT")
    @test isdir(joinpath(depotpath, "development", "TestPackage"))
    @test isfile(joinpath(depotpath, "development", "TestPackage", "README.md"))
    @test isfile(joinpath(depotpath, "development", "TestPackage", "LICENSE.md"))
    @test isfile(joinpath(depotpath, "development", "TestPackage", "src", "TestPackage.jl"))
    @test isfile(joinpath(depotpath, "development", "TestPackage", "test", "runtests.jl"))

    # TODO Add some dependencies to package and make a version of the package

    # TODO Try adding and `using` the package and check that Pkg3 understands the registry.

end


text:
module MLP
    using Random
    using Statistics
    using Plots
    using ProgressMeter

    export
        SimpleMLP,
        fit!,
        adaptive_fit!,
        predict,
        save_figures

    abstract type NeuralNetwork end
    abstract type FeedForwardNN <: NeuralNetwork end
    abstract type MultiLayerPerceptron <: FeedForwardNN end

    error(ŷ, y) = ŷ - y
    σ(x) = 1 / (1 + exp(-x))
    σ′(x) = x * (1 - x)

    mutable struct SimpleMLP{T} <: MultiLayerPerceptron where T <: AbstractFloat
        train_data
        test_data

        loss
        activation
        activation_gradient
        n_epochs::Integer
        η::T

        wL₁::Matrix{T}
        wL₂::Matrix{T}
        loss_history::Vector{T}
        train_accuracy_history::Vector{T}
        test_accuracy_history::Vector{T}

        function SimpleMLP{T}(
                data_X::Matrix{T},
                data_Y::Matrix{T},
                ;
                train_test_ratio::AbstractFloat=0.8,
                loss=error,
                activation=σ,
                activation_gradient=σ′,
                n_hidden_layers::Integer=size(data_X, 2),
                n_epochs::Integer=100,
                η::T=convert(T, 1e-1),
        ) where T <: AbstractFloat
            data = shuffle(collect(zip(eachrow(data_X), eachrow(data_Y))))
            n_train_data = round(Integer, length(data) * train_test_ratio)
            new(data[1:n_train_data], data[n_train_data+1:end],
                loss, activation, activation_gradient, n_epochs, η,
                rand(n_hidden_layers, size(data_X, 2)), rand(size(data_Y, 2), n_hidden_layers), T[], T[], T[])
        end
    end

    function fit!(mlp::MultiLayerPerceptron)
        last_up_wL₂ = zeros(size(mlp.wL₂))
        new_up_wL₂ = zeros(size(mlp.wL₂))
        last_up_wL₁ = zeros(size(mlp.wL₁))
        new_up_wL₁ = zeros(size(mlp.wL₁))
        @showprogress for c ∈ 1:mlp.n_epochs
            for (x, y) ∈ mlp.train_data
                x = hcat(x)
                y = hcat(y)
                aL₁ = mlp.activation.(mlp.wL₁ * x)
                aL₂ = mlp.activation.(mlp.wL₂ * aL₁)
                e_nn = mlp.loss(aL₂, y)
                dL₂ = mlp.activation_gradient.(aL₂) .* e_nn
                dL₁ = mlp.activation_gradient.(aL₁) .* (mlp.wL₂' * dL₂)

                ΔL₂ = dL₂ * aL₁'
                ΔL₁ = dL₁ * x'
                ΔL₂ *= mlp.η
                ΔL₁ *= mlp.η
                mlp.wL₂ -= ΔL₂
                mlp.wL₁ -= ΔL₁
                push!(mlp.loss_history, mean(e_nn))
            end
            check_accuracy!(mlp)
        end
    end

    function adaptive_fit!(mlp::MultiLayerPerceptron)
        ϵ = eps()
        GwL₂ = zeros(size(mlp.wL₂))
        GwL₁ = zeros(size(mlp.wL₁))
        @showprogress for c ∈ 1:mlp.n_epochs
            for (x, y) ∈ mlp.train_data
                x = hcat(x)
                y = hcat(y)
                aL₁ = mlp.activation.(mlp.wL₁ * x)
                aL₂ = mlp.activation.(mlp.wL₂ * aL₁)
                e_nn = mlp.loss(aL₂, y)
                dL₂ = mlp.activation_gradient.(aL₂) .* e_nn
                dL₁ = mlp.activation_gradient.(aL₁) .* (mlp.wL₂' * dL₂)

                ΔL₂ = dL₂ * aL₁'
                ΔL₁ = dL₁ * x'
                GwL₂ += ΔL₂ .^ 2
                GwL₁ += ΔL₁ .^ 2
                @. ΔL₂ *= mlp.η / (√GwL₂ + ϵ)
                @. ΔL₁ *= mlp.η / (√GwL₁ + ϵ)
                mlp.wL₂ -= ΔL₂
                mlp.wL₁ -= ΔL₁
                push!(mlp.loss_history, mean(e_nn))
            end
            check_accuracy!(mlp)
        end
    end

    function predict(mlp::MultiLayerPerceptron, x::Matrix)
        mlp.activation.(mlp.wL₂ * mlp.activation.(mlp.wL₁ * x))
    end

    function save_figures(mlp::MultiLayerPerceptron, fn::AbstractString="figs/")
        savefig(plot(
            mlp.loss_history,
            label="loss",
            title="loss",
            legend=:outertopright,
            size=(1024, 512),
        ), "$(fn)loss_fig.png")
        savefig(plot(
            [mlp.train_accuracy_history mlp.test_accuracy_history],
            label=["train accuracy" "test accuracy"],
            title="accuracy",
            legend=:outertopright,
            size=(1024, 512),
        ), "$(fn)accuracy_fig.png")
    end

    function check_accuracy!(mlp::MultiLayerPerceptron)
        train_true_prediction = 0
        for (x, y) ∈ mlp.train_data
            x = hcat(x)
            y = hcat(y)
            out = predict(mlp, x)
            train_true_prediction += argmax(out) == argmax(y)
        end
        train_acc = train_true_prediction / length(mlp.train_data)
        push!(mlp.train_accuracy_history, train_acc)

        test_true_prediction = 0
        for (x, y) ∈ mlp.test_data
            x = hcat(x)
            y = hcat(y)
            out = predict(mlp, x)
            test_true_prediction += argmax(out) == argmax(y)
        end
        test_acc = test_true_prediction / length(mlp.test_data)
        push!(mlp.test_accuracy_history, test_acc)
    end
end


text:
module Simulations

    using ProgressMeter

    using ..Equations
    using ..Integrators
    using ..Solutions

    export Simulation, run!

    include("simulations/simulation.jl")

end


text:
using Base.Test
using SQLAlchemy

engine = createengine("sqlite:///:memory:")
metadata = MetaData()
users = Table("Users", metadata,
            Column("name", SQLString()),
            Column("age", SQLFloat()))
createall(metadata, engine)
db = connect(engine)
db(insert(users), name="Alice", age=27)
db(insert(users), name="Bob", age=31)
res = db(select([users])) |> fetchall
@test res[1][:name]|>get == "Alice"
@test res[2][2]|>get == 31

include("orm.jl")


text:
module TestReadablestring

using Orihime
using Orihime.Internal: format_with_thousands_separator
using Test

macro ex_str(str)
    expr = esc(Meta.parse(str))
    :($expr => $str)
end

macro tsep_str(str)
    value = parse(Int, replace(str, "," => ""))
    :($value => $str)
end

function test_int_expr()
    @testset "$str" for (value, str) in [
        ex"0",
        ex"512",
        ex"2^10",
        ex"2^12",
        ex"-2^10",
        ex"-2^12",
        ex"10^3",
        ex"10^4",
        ex"-10^3",
        ex"-10^4",
        tsep"1,234",
        tsep"12,345",
        tsep"123,456",
        tsep"1,234,567",
        tsep"1,234,567,890",
        tsep"10,101,010,101",
        tsep"101,010,101,010",
        tsep"-1,234",
        tsep"-12,345",
        tsep"-123,456",
        tsep"-1,234,567",
        tsep"-1,234,567,890",
        tsep"-101,010,101,010",
    ]
        @test Orihime.readablestring(value) == str
    end
end

function test_no_thousands_separator()
    @testset "$value" for value in [1:10; 0; -1:-1:-10; 100; -100]
        @test format_with_thousands_separator(value) == string(value)
    end
end

end  # module


text:
struct SDDEProblem{uType,tType,lType,lType2,isinplace,P,NP,F,G,H,K,ND} <:
                          AbstractSDDEProblem{uType,tType,lType,isinplace,ND}
  f::F
  g::G
  u0::uType
  h::H
  tspan::tType
  p::P
  noise::NP
  constant_lags::lType
  dependent_lags::lType2
  kwargs::K
  noise_rate_prototype::ND
  seed::UInt64
  neutral::Bool
  order_discontinuity_t0::Rational{Int}

  @add_kwonly function SDDEProblem{iip}(f::AbstractSDDEFunction{iip}, g, u0, h, tspan, p = NullParameters();
                                       noise_rate_prototype = nothing, noise= nothing, seed = UInt64(0),
                                       constant_lags = (), dependent_lags = (),
                                       neutral = f.mass_matrix !== I && det(f.mass_matrix) != 1,
                                       order_discontinuity_t0 = 0//1,
                                       kwargs...) where {iip}
    _tspan = promote_tspan(tspan)
    new{typeof(u0),typeof(_tspan),typeof(constant_lags),typeof(dependent_lags),isinplace(f),
        typeof(p),typeof(noise),typeof(f),typeof(g),typeof(h),typeof(kwargs),typeof(noise_rate_prototype)}(
        f, g, u0, h, _tspan, p, noise, constant_lags, dependent_lags, kwargs, noise_rate_prototype, seed, neutral, order_discontinuity_t0)
  end

  function SDDEProblem{iip}(f::AbstractSDDEFunction{iip}, g , h, tspan::Tuple, p = NullParameters();
                           order_discontinuity_t0 = 1//1, kwargs...) where iip
    SDDEProblem{iip}(f, g, h(p, first(tspan)), h, tspan, p;
                    order_discontinuity_t0 = max(1//1, order_discontinuity_t0), kwargs...)
  end

  function SDDEProblem{iip}(f,g, args...; kwargs...) where iip
    SDDEProblem{iip}(convert(SDDEFunction{iip}, f,g),g, args...; kwargs...)
  end
end

SDDEProblem(f, g, args...; kwargs...) =
  SDDEProblem(convert(SDDEFunction, f, g), g, args...; kwargs...)

SDDEProblem(f::AbstractSDDEFunction, args...; kwargs...) =
  SDDEProblem{isinplace(f)}(f, args...; kwargs...)


text:
# Tests for Dirichlet distribution

using Distributions
using Base.Test

d = Dirichlet(3, 2.0)

@test dim(d) == 3
@test d.alpha == [2.0, 2.0, 2.0]
@test d.alpha0 == 6.0

@test_approx_eq mean(d) fill(1.0/3, 3)
@test_approx_eq cov(d) [8 -4 -4; -4 8 -4; -4 -4 8] / (36 * 7)
@test_approx_eq var(d) diag(cov(d))

@test_approx_eq pdf(d, [0.2, 0.3, 0.5]) 3.6
@test_approx_eq pdf(d, [0.4, 0.5, 0.1]) 2.4
@test_approx_eq logpdf(d, [0.2, 0.3, 0.5]) log(3.6)
@test_approx_eq logpdf(d, [0.4, 0.5, 0.1]) log(2.4)

x = rand(d, 100)
p = pdf(d, x)
lp = logpdf(d, x)
for i in 1 : size(x, 2)
	@test_approx_eq lp[i] logpdf(d, x[:,i])
	@test_approx_eq p[i] pdf(d, x[:,i])
end

v = [2.0, 1.0, 3.0]
d = Dirichlet(v)

@test dim(d) == length(v)
@test d.alpha == v
@test d.alpha0 == sum(v)

@test_approx_eq mean(d) v / sum(v)
@test_approx_eq cov(d) [8 -2 -6; -2 5 -3; -6 -3 9] / (36 * 7)
@test_approx_eq var(d) diag(cov(d))

@test_approx_eq pdf(d, [0.2, 0.3, 0.5]) 3.0
@test_approx_eq pdf(d, [0.4, 0.5, 0.1]) 0.24
@test_approx_eq logpdf(d, [0.2, 0.3, 0.5]) log(3.0)
@test_approx_eq logpdf(d, [0.4, 0.5, 0.1]) log(0.24)

x = rand(d, 100)
p = pdf(d, x)
lp = logpdf(d, x)
for i in 1 : size(x, 2)
	@test_approx_eq p[i] pdf(d, x[:,i])
	@test_approx_eq lp[i] logpdf(d, x[:,i])
end

# Sampling

x = rand(d)
@test isa(x, Vector{Float64})
@test length(x) == 3

x = rand(d, 10)
@test isa(x, Matrix{Float64})
@test size(x) == (3, 10)

# Test MLE

n = 10000
x = rand(d, n)
x = x ./ sum(x, 1)

r = fit_mle(Dirichlet, x)
@test_approx_eq_eps r.alpha d.alpha 0.25

r = fit_mle(Dirichlet, x, fill(2.0, n))
@test_approx_eq_eps r.alpha d.alpha 0.25


text:
using IterableTables
using NamedTuples
using DataFrames
using DataStreams
using CSV
using SQLite
using Feather
using Base.Test

@testset "DataStreams" begin

csv_source = CSV.Source(joinpath(dirname(@__FILE__),"testdata.csv"))
df = DataFrame(csv_source)

@test size(df) == (3,3)
@test isa(df[:a], Vector{Union{Int,Missing}})
@test isa(df[:b], Vector{Union{Float64,Missing}})
@test isa(df[:c], Vector{Union{String,Missing}})
@test df[:a] == [1,2,3]
@test df[:b] == [1.,2.,3.]
@test df[:c] == ["A","B","C"]

feather_source = Feather.Source(joinpath(Pkg.dir("Feather"),"test", "data", "CO2.feather"))
df_feather = DataFrame(feather_source)

@test size(df_feather) == (84,5)
@test eltype(df_feather[:Plant]) <: String
@test eltype(df_feather[:Type]) <:String
@test eltype(df_feather[:Treatment]) <:String
@test eltype(df_feather[:conc]) <:Float64
@test eltype(df_feather[:uptake]) <:Float64
@test df_feather[1,:Plant] == "Qn1"
@test df_feather[1,:Type] == "Quebec"
@test df_feather[1,:Treatment] == "nonchilled"
@test df_feather[1,:conc] == 95.0
@test df_feather[1,:uptake] == 16.

sqlite_source = SQLite.Source(SQLite.DB(joinpath(Pkg.dir("SQLite"), "test", "Chinook_Sqlite.sqlite")), "SELECT * FROM Employee;")
df_sqlite = DataFrame(sqlite_source)
@test size(df_sqlite) == (8,15)
@test isa(df_sqlite[:EmployeeId], Vector{Union{Int64,Missing}})
@test df_sqlite[1, :EmployeeId] == 1

end


text:
function decompose(obj::MatrixDecomposer, arg0::RealMatrix)
    return jcall(obj, "decompose", DecompositionSolver, (RealMatrix,), arg0)
end


text:
# Automatically generated using Clang.jl

const ssize_t = Cptrdiff_t
const iio_context = Cvoid
const iio_device = Cvoid
const iio_channel = Cvoid
const iio_buffer = Cvoid
const iio_context_info = Cvoid
const iio_scan_context = Cvoid
const iio_scan_block = Cvoid

@cenum iio_chan_type::UInt32 begin
    IIO_VOLTAGE = 0
    IIO_CURRENT = 1
    IIO_POWER = 2
    IIO_ACCEL = 3
    IIO_ANGL_VEL = 4
    IIO_MAGN = 5
    IIO_LIGHT = 6
    IIO_INTENSITY = 7
    IIO_PROXIMITY = 8
    IIO_TEMP = 9
    IIO_INCLI = 10
    IIO_ROT = 11
    IIO_ANGL = 12
    IIO_TIMESTAMP = 13
    IIO_CAPACITANCE = 14
    IIO_ALTVOLTAGE = 15
    IIO_CCT = 16
    IIO_PRESSURE = 17
    IIO_HUMIDITYRELATIVE = 18
    IIO_ACTIVITY = 19
    IIO_STEPS = 20
    IIO_ENERGY = 21
    IIO_DISTANCE = 22
    IIO_VELOCITY = 23
    IIO_CONCENTRATION = 24
    IIO_RESISTANCE = 25
    IIO_PH = 26
    IIO_UVINDEX = 27
    IIO_ELECTRICALCONDUCTIVITY = 28
    IIO_COUNT = 29
    IIO_INDEX = 30
    IIO_GRAVITY = 31
    IIO_POSITIONRELATIVE = 32
    IIO_PHASE = 33
    IIO_MASSCONCENTRATION = 34
    IIO_CHAN_TYPE_UNKNOWN = 2147483647
end

@cenum iio_modifier::UInt32 begin
    IIO_NO_MOD = 0
    IIO_MOD_X = 1
    IIO_MOD_Y = 2
    IIO_MOD_Z = 3
    IIO_MOD_X_AND_Y = 4
    IIO_MOD_X_AND_Z = 5
    IIO_MOD_Y_AND_Z = 6
    IIO_MOD_X_AND_Y_AND_Z = 7
    IIO_MOD_X_OR_Y = 8
    IIO_MOD_X_OR_Z = 9
    IIO_MOD_Y_OR_Z = 10
    IIO_MOD_X_OR_Y_OR_Z = 11
    IIO_MOD_LIGHT_BOTH = 12
    IIO_MOD_LIGHT_IR = 13
    IIO_MOD_ROOT_SUM_SQUARED_X_Y = 14
    IIO_MOD_SUM_SQUARED_X_Y_Z = 15
    IIO_MOD_LIGHT_CLEAR = 16
    IIO_MOD_LIGHT_RED = 17
    IIO_MOD_LIGHT_GREEN = 18
    IIO_MOD_LIGHT_BLUE = 19
    IIO_MOD_QUATERNION = 20
    IIO_MOD_TEMP_AMBIENT = 21
    IIO_MOD_TEMP_OBJECT = 22
    IIO_MOD_NORTH_MAGN = 23
    IIO_MOD_NORTH_TRUE = 24
    IIO_MOD_NORTH_MAGN_TILT_COMP = 25
    IIO_MOD_NORTH_TRUE_TILT_COMP = 26
    IIO_MOD_RUNNING = 27
    IIO_MOD_JOGGING = 28
    IIO_MOD_WALKING = 29
    IIO_MOD_STILL = 30
    IIO_MOD_ROOT_SUM_SQUARED_X_Y_Z = 31
    IIO_MOD_I = 32
    IIO_MOD_Q = 33
    IIO_MOD_CO2 = 34
    IIO_MOD_VOC = 35
    IIO_MOD_LIGHT_UV = 36
    IIO_MOD_LIGHT_DUV = 37
    IIO_MOD_PM1 = 38
    IIO_MOD_PM2P5 = 39
    IIO_MOD_PM4 = 40
    IIO_MOD_PM10 = 41
    IIO_MOD_ETHANOL = 42
    IIO_MOD_H2 = 43
    IIO_MOD_O2 = 44
end

struct iio_data_format
    length::UInt32
    bits::UInt32
    shift::UInt32
    is_signed::Bool
    is_fully_defined::Bool
    is_be::Bool
    with_scale::Bool
    scale::Cdouble
    repeat::UInt32
end


text:
module Tissue

import Base.Threads: @spawn, @threads, threadid
using MLStyle

export Graph, start, stop, wait_until_done
export @graph, @calculator, @bindstreams

"""
    process(calculator, ...; graph::Graph)

Process the input streams into an output.

The first argument must be the associated calculator. Any other argument defines an input stream to the calculator. An input stream will be named with the exact same name as the corresponding argument. At all times, data coming from mutliple input streams into a `process()` method is guaranteed to be derived from the same datum generated by the *source* calculator.

Optionally, if you need a reference to the `graph` object in the `process()` method, you can add it as a keyword argument named `graph` and Tissue.jl will supply a reference to the active graph object.

The value returned from the method is the value sent to the calculator's output stream. Hence, any calculator that binds one of their input streams with [`@bindstreams`](@ref) to this calculator's output stream will receive the return value of this method as an argument. `nothing` is treated specially: it means that the calculator will not process the given input. Any downstream calculator that depends on other calculators too will simply drop the value taken from their output stream.

It is important to remember that each calculator runs in its own task. Therefore, best practices about concurrency must be applied. Notably, the stream arguments of `process()` could possibly be accessed by multiple threads simultaneously. Therefore, it is highly recommended to make a [`deepcopy`](https://docs.julialang.org/en/v1/base/base/#Base.deepcopy) of an argument before mutating it. This will also make your calculator more reusable across graphs.

Multiple `process()` methods can be defined per calculator, allowing the calculator to be used in different graph topologies.

# Arguments
- `calculator`: the associated calculator
- `...`: the input streams

# Examples
```julia
using Base: deepcopy

function process(calc::MyCalculator, in_stream; graph)
    if some_condition
        stop(graph)
    end

    in_stream_copy = deepcopy(in_stream)

    mutate(in_stream_copy)

    return in_stream_copy
end
```
"""
function process(calculator)
    throw("Called unimplemented process(). This shouldn't happen.")
end

"""
    close(calculator)

Perform cleanup for the calculator. Optional.

You can define a method for your calculator type to perform any necessary cleanup. Called by [`wait_until_done(graph)`](@ref) on each calculator to perform cleanup.

# Examples
```julia
struct GoofyCalculator
    resource
    function GoofyCalculator()
        resource = acquire_resource()
        new(resource)
    end
end

function process(c::GoofyCalculator, some_stream)
    use_resource(c.resource)
end

function close(c::GoofyCalculator)
    release_resource(c.resource)
end
```
"""
function close(calculator)
    # nothing
end

include("graph.jl")
include("core.jl")
include("macros.jl")

end # module Tissue


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

julia_versions = [v"1.6.3", v"1.7"]

name = "XyceWrapper"
version = v"0.3.0"

# Collection of sources required to complete build
sources = [
    DirectorySource("./src"),
]

# See https://github.com/JuliaLang/Pkg.jl/issues/2942
# Once this Pkg issue is resolved, this must be removed
uuid = Base.UUID("a83860b7-747b-57cf-bf1f-3e79990d037f")
delete!(Pkg.Types.get_last_stdlibs(v"1.6.3"), uuid)

# Bash recipe for building across all platforms
script = raw"""
mkdir -p build
cd build
cmake \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_FIND_ROOT_PATH=${prefix} \
    -DJulia_PREFIX=${prefix} \
    -DCMAKE_BUILD_TYPE=Release \
    ..
cmake --build . --config Release --target install -- -j${nproc}
install_license /usr/share/licenses/MIT
"""

include("../../L/libjulia/common.jl")
platforms = vcat(libjulia_platforms.(julia_versions)...)
platforms = expand_cxxstring_abis(platforms)
filter!(p -> arch(p) != "armv6l", platforms)
filter!(p -> !(arch(p) == "aarch64" && Sys.isapple(p)), platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libxycelib", :xycelib),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Xyce_jll"; compat="^7.4.0"),
    Dependency("libcxxwrap_julia_jll"),
    BuildDependency("libjulia_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat = "1.6", preferred_gcc_version=v"7")