text:
using Test

using SpectralUnmixing

datafile = "../data/basic_endmember_library.csv"
classname = "Class"

@info "Basic Loading"
lib = load_data!(SpectralLibrary(datafile, classname))
@test isnothing(lib) == false

@time @test isnothing(lib.spectra) == false
@time @test length(size(lib.spectra)) == 2

@info "Filter By Class"
tmplib = deepcopy(lib)
filter_by_class!(tmplib)
@time @test length(lib.classes) == length(tmplib.classes)

tmplib =  load_data!(SpectralLibrary(datafile, classname, 2, 0, [unique(lib.classes)[1]]))
filter_by_class!(tmplib)
@time @test length(lib.classes) != length(tmplib.classes)
@time @test length(unique(tmplib.classes)) == 1
tmplib=nothing

@info "Remove wavelength region inplace"

tmplib =  load_data!(SpectralLibrary(datafile, classname))

remove_wavelength_region_inplace!(tmplib)
@time @test size(lib.spectra) != size(tmplib.spectra)

@info "Wavelength interpolation"
tmplib = deepcopy(lib)
new_wl = 400 .+ Vector{Float64}(1:200)*10.
interpolate_library_to_new_wavelengths!(tmplib, new_wl)
@time @test size(tmplib.spectra)[2] == length(new_wl)
@time @test size(lib.spectra) != size(tmplib.spectra)

@info "Scale Library"
tmplib = deepcopy(lib)
@time scale_library!(tmplib)
@test lib.spectra == tmplib.spectra

@time scale_library!(tmplib, 2)
@test lib.spectra != tmplib.spectra

tmplib = deepcopy(lib)
@time brightness_normalize!(tmplib)
@test lib.spectra != tmplib.spectra

@info "Simulate Mixture - class-even type"
remove_wavelength_region_inplace!(lib)
class_idx = []
for uc in lib.class_valid_keys
    push!(class_idx, (1:size(lib.classes)[1])[lib.classes .== uc])
end
num_components = 3
seed = 13
@time simulated_rfl, true_complete_fractions, true_mixture_fractions = simulate_pixel(lib, num_components, "class-even", seed)
@test sum(isnothing.(simulated_rfl)) + sum(isnan.(simulated_rfl)) == 0
println(lib.spectra' * true_complete_fractions)
println(simulated_rfl[:])
@test sum(lib.spectra' * true_complete_fractions .≈ simulated_rfl[:]) == length(simulated_rfl[:])

for mode in ["sma", "sma-best", "mesma", "mesma-best"]
    n_mc = 2
    mode = "sma-best"
    num_endmembers=[3]
    normalization="brightness"
    optimization="bvls"

    max_combinations=10000
    combination_type="class_even"

    unmixing_library = lib

    options = []

    class_idx = []
    if combination_type == "class-even"
        for uc in library.class_valid_keys
            push!(class_idx, (1:size(unmixing_library.classes)[1])[unmixing_library.classes .== uc])
        end
    end

    @info "Unmix Pixel - Mode: " * mode
    @time mr, mv, cfr, cfv = unmix_pixel(unmixing_library, simulated_rfl, nothing, class_idx, options, mode, n_mc,
            num_endmembers, normalization, optimization, max_combinations, combination_type)
    @test sum(mr[1:end-1]) ≈ 1
    @test size(mr) == size(mv)
    @test sum(cfr[1:end-1]) ≈ 1
    @test size(cfr) == size(cfv)
end

max_stars_repo_name: fotxoa-geo/SpectralUnmixing

max_issues_repo_name: fotxoa-geo/SpectralUnmixing

max_forks_repo_name: fotxoa-geo/SpectralUnmixing


text:
# Copyright 2014 Mitchell Kember. Subject to the MIT License.
# Google Code Jam Round 1A 2014
# Problem B: ???

function mindel(r, t)
	d = 0
	stack = Int[]
	push!(stack, nothing)
	top = r
	while top !== nothing
		if top
	end
end

function main()
	n = int(readline())
	for i = 1:n
		nn = int(readline())
		tree = Array((Int, Int), nn-1)
		for j = 1:nn-1
			con = int(split(readline()))
			tree[j] = minmax(con[1], con[2])
		end
		println("Case #$i: ", mindel(tree))
	end
end

main()

ext: jl


text:
# Run tests on functions in src/by_reference_doc/playlists.jl

@testset verbose = true "GET-request endpoints for playlists" begin

    # Input arguments from composite types defined in src/types.jl
    playlist_id = SpPlaylistId()
    user_id = SpUserId()
    category_id = SpCategoryId()

    # Cycle through different input keywords for testing
    countries = ["US", "NL", "DE"]
    offsets = [10, 35]

    # Start test blocks
    @test ~isempty(Spotify.playlist_get(playlist_id)[1])

    # Does not run first time
    @test_throws BoundsError ~isempty(Spotify.playlist_get_current_user(; offset = 0)[1])

    @testset "For offset = $(num_offset)" for num_offset in offsets

        @test ~isempty(Spotify.playlist_get_tracks(playlist_id; offset = num_offset)[1])
        @test ~isempty(Spotify.playlist_get_user(user_id; offset = num_offset)[1])

        # Works okay second time
        @test ~isempty(Spotify.playlist_get_current_user(; offset = num_offset)[1])
    end

    @testset "For country = $(country_id)" for country_id in countries

        @test ~isempty(Spotify.playlist_get_featured(country = country_id)[1])
        @test ~isempty(Spotify.playlist_get_category(category_id; country = country_id)[1])

    end

    @test ~isempty(Spotify.playlist_get_cover_image(playlist_id)[1])

end

# Test if 404 exception handling works by using an incorrect playlist ID

@testset "Test 404 Not Found exception" begin

    @test_logs (:info, "404 Not Found - Check if input arguments are okay") match_mode=:any Spotify.playlist_get_tracks("37i9dQZF1E4vUblDJzzkV3")[1]

end


text:
using Test, Distributions
@testset "Normal" begin
    @test isa(convert(Normal{Float64}, Float16(0), Float16(1)),
              Normal{Float64})
    @test Inf === logpdf(Normal(0, 0), 0)
    @test -Inf === logpdf(Normal(), Inf)
    @test iszero(logcdf(Normal(0, 0), 0))
    @test iszero(logcdf(Normal(), Inf))
    @test -Inf === logccdf(Normal(0, 0), 0)
    @test iszero(logccdf(Normal(eps(), 0), 0))
    @test -Inf === quantile(Normal(), 0)
    @test iszero(quantile(Normal(), 0.5))
    @test Inf === quantile(Normal(), 1)
    @test -Inf === quantile(Normal(0, 0), 0)
    @test iszero(quantile(Normal(0, 0), 0.75))
    @test Inf === quantile(Normal(0, 0), 1)
    @test -Inf === quantile(Normal(0.25, 0), 0)
    @test 0.25 == quantile(Normal(0.25, 0), 0.95)
    @test Inf === quantile(Normal(0.25, 0), 1)
    @test Inf === cquantile(Normal(), 0)
    @test iszero(cquantile(Normal(), 0.5))
    @test -Inf === cquantile(Normal(), 1)
    @test Inf === cquantile(Normal(0, 0), 0)
    @test iszero(cquantile(Normal(0, 0), 0.75))
    @test -Inf === cquantile(Normal(0, 0), 1)
    @test Inf === cquantile(Normal(0.25, 0), 0)
    @test 0.25 == cquantile(Normal(0.25, 0), 0.95)
    @test -Inf === cquantile(Normal(0.25, 0), 1)
    @test -Inf === invlogcdf(Normal(), -Inf)
    @test isnan(invlogcdf(Normal(), NaN))
    @test Inf === invlogccdf(Normal(), -Inf)
    @test isnan(invlogccdf(Normal(), NaN))
    # test for #996 being fixed
    let d = Normal(0, 1), x = 1.0, ∂x = 2.0
        @inferred cdf(d, ForwardDiff.Dual(x, ∂x)) ≈ ForwardDiff.Dual(cdf(d, x), ∂x * pdf(d, x))
    end
end


text:
# Created by Alex Dowling (alexdowling.net) while at the University of Wisconsin-Madison

##### This example solves an optimization in parallel using pmap()
# More info: http://julia.readthedocs.org/en/latest/manual/parallel-computing/
# Use "julia -p n" to start Julia with n

### Import optimization problem definition/solve in a function

push!(LOAD_PATH,pwd())
using DummyModule

### Test function

my = solveOptProb2((1,1))

### Iterate over possible values of a and b
aOpt = 0:5
bOpt = 1:3

aN = length(aOpt)
bN = length(bOpt)

# Assemble array of input parameters for the optimization problem
P = Array(Tuple,aN*bN)

for i in 1:length(aOpt)
	for j in 1:length(bOpt)
		P[(i-1)*bN + j] = (aOpt[i],bOpt[j])
	end
end

# Display P
print(P)

# Solve the optimization problem in parallel
# Note: results is an array of tuples, with the same dimensions as P
results = pmap(solveOptProb2, P)

for i = 1:length(results)
	DummyModule.print(results[i])
end


text:
#=

Code to implement piecewise linear interpolation for an arbitrary number
of dimensions

=#

function piecewise_linear_nodes(n::S,domain = [1.0,-1.0]) where {S <: Integer}

  if n <= 0
    error("The number of nodes must be positive.")
  end

  nodes = zeros(n)

  if isodd(n)
    nodes[Int((n-1)/2)+1] = (domain[1]+domain[2])/2.0
  end

  if n == 1

    return nodes

  elseif n == 2

    return [domain[2], domain[1]]

  else

    nodes[1] = domain[2]
    nodes[n] = domain[1]

    inc = (domain[1]-domain[2])/(n-1)
    for i = 2:div(n,2)
      nodes[i]     = domain[2] + (i-1)*inc
      nodes[n-i+1] = domain[1] - (i-1)*inc
    end

    return nodes

  end

end

const linear_nodes = piecewise_linear_nodes

function bracket_nodes(x::Array{T,1},point::R) where {T <: AbstractFloat,R<:Number}

  if real(point) <= x[1] # Real is used because complex numbers are occasionally used in NLboxsolve.jl
    return (1,2)
  elseif real(point) >= x[end]
    return (length(x) - 1, length(x))
  else
    y = 0
    for i in x
      if i < real(point)
        y += 1
      else
        break
      end
    end
    return (y, y+1)
  end

end

function bracket_nodes(x::Union{NTuple{N,Array{T,1}},Array{Array{T,1},1}},point::Array{R,1}) where {T <: AbstractFloat, N,R<:Number}

  bracketing_nodes = Array{Int64,2}(undef,2,length(point))
  for i = 1:length(point)
    bracketing_nodes[:,i] .= bracket_nodes(x[i],point[i])
  end

  return bracketing_nodes

end

function piecewise_linear_weight(x::Array{T,1},point::R) where {T<:AbstractFloat,R<:Number}

    bracketing_nodes = bracket_nodes(x,point)
    weight = (point-x[bracketing_nodes[1]])/(x[bracketing_nodes[2]]-x[bracketing_nodes[1]])

    return weight

end

function piecewise_linear_weights(x::Union{NTuple{N,Array{T,1}},Array{Array{T,1},1}},point::Array{R,1}) where {T<:AbstractFloat,N,R<:Number}

    weights = Array{R,1}(undef,length(point))
    for i = 1:length(point)
        weights[i] = piecewise_linear_weight(x[i],point[i])
    end

    return weights

end

function select_bracketing_nodes(bounds::Array{S,2} where {S <: Integer})

  d = size(bounds,2)

  number_relevant_grid_points = 2^d
  bracketing_grid_points      = Array{Int64}(undef,number_relevant_grid_points,d)

  for i = 1:d
    bracketing_grid_points[:,i] .= repeat(repeat(bounds[:,i],inner=[2^(d-i),1]),2^(i-1))[:]
  end

  return bracketing_grid_points

end

function select_relevant_data(y::AbstractArray{T,N},bracketing_grid_points::Array{S,2}) where {T <: AbstractFloat, S <: Integer, N}

  data = zeros(size(bracketing_grid_points,1))
  for i = 1:length(data)
    data[i] = y[CartesianIndex(tuple(bracketing_grid_points[i,:]...))]
  end

  return data

end

function piecewise_linear_evaluate(y::AbstractArray{T,N},x::Union{NTuple{N,Array{T,1}},Array{Array{T,1},1}},point::Union{R,Array{R,1}}) where {T<:AbstractFloat,N,R<:Number}

  b = bracket_nodes(x,point)
  w = piecewise_linear_weights(x,point)

  d = size(b,2)

  relevant_points = select_bracketing_nodes(b)
  data = select_relevant_data(y,relevant_points)

  for j = d:-1:1

    new_data = zeros(R,div(length(data),2))
    for i = 1:length(new_data)
      new_data[i] = data[2*(i-1)+1] + w[j]*(data[2*i]-data[2*(i-1)+1])
    end

    data = copy(new_data)

  end

  return data[1]

end

function piecewise_linear_evaluate(y::AbstractArray{T,N},x::Union{NTuple{N,Array{T,1}},Array{Array{T,1},1}}) where {T <: AbstractFloat, N}

  function approximating_function(point::Union{R,Array{R,1}}) where {R<:Number}

    return piecewise_linear_evaluate(y,x,point)

  end

  return approximating_function

end

function piecewise_linear_evaluate(y::AbstractArray{T,N},x::Union{NTuple{N,Array{T,1}},Array{Array{T,1},1}},point::Union{R,Array{R,1}},integrals::Union{T,Array{T,1}}) where {T<:AbstractFloat,N,R<:Number}

  # This function is only needed to facilitate compatibility with SolveDSGE

  b = bracket_nodes(x,point)
  w = piecewise_linear_weights(x,point)
  w = w.*integrals

  d = size(b,2)

  relevant_points = select_bracketing_nodes(b)
  data = select_relevant_data(y,relevant_points)

  for j = d:-1:1

    new_data = zeros(R,Int(length(data)/2))
    for i = 1:length(new_data)
      new_data[i] = data[2*(i-1)+1] + w[j]*(data[2*i]-data[2*(i-1)+1])
    end

    data = copy(new_data)

  end

  return data[1]

end

# Functions to handle the 1-D case

function piecewise_linear_evaluate(y::Array{T,1},x::Array{T,1},point::R) where {T<:AbstractFloat,R<:Number}

    b = bracket_nodes(x,point)
    w = piecewise_linear_weight(x,point)

    y_estimate = y[b[1]] + w*(y[b[2]]-y[b[1]])

    return y_estimate

end

function piecewise_linear_evaluate(y::Array{T,1},nodes::Array{T,1}) where {T <: AbstractFloat}

  function approximating_function(point::R) where {R<:Number}

      return piecewise_linear_evaluate(y,nodes,point)

  end

  return approximating_function

end

# Function to transform from one uniform grid to another

function grid_reshape(f::AbstractArray{T,N},grid::NTuple{N,Array{T,1}}) where {T<:AbstractFloat,N}

  #1. Construct the old grid

  old_grid = Array{Array{T,1},1}(undef,N)
  for i = 1:N
    piecewise_linear_nodes
    old_grid[i] = piecewise_linear_nodes(size(f,i),[grid[i][end],grid[i][1]])
  end

  #2. Initialize the new y

  y = Array{T,N}(undef,tuple(length.(grid)...))

  #3. Interpolate to fill y

  point = Array{T,1}(undef,N)
  for i in CartesianIndices(y)
    for j = 1:N
      point[j] = grid[j][i[j]]
    end
    y[i] = piecewise_linear_evaluate(f,old_grid,point)
  end

  return y

end


text:
using GeneralizedDynamicsFromData
using OrderedCollections
using JLD2
using FileIO

experiment_name = "selkov_statistics_MSE"
repetitions = 100

weight_decay = 1e-4
η = 1e-2
η_decay_rate = 0.5
η_decay_step = 800
η_limit = 1e-4

construct_optimiser() = Flux.Optimiser(WeightDecay(weight_decay),
                                       ExpDecay(η,η_decay_rate,η_decay_step,η_limit),
                                       ADAM())
construct_loss(θ, y, predict) = mse_loss(θ, y, predict)

net_config = OrderedDict([
    :inputs => 2,
    :outputs => 1,
    :neurons => 16,
    :layers => 1,
    :non_lin => tanh,
    :initialization => Flux.glorot_normal
])

problem_oscillation = Dict([:equation => selkov,
                           :parameters => Float32[0.1, 0.6],
                           :u0 => Float32[1.0, 1.0],
                           :tspan => (0.0f0, 15.0f0),
                           :ts => 0.1,
                           :solver => Tsit5,
                           :optimizer => construct_optimiser,
                           :max_iter => 3000,
                           :loss => construct_loss]
                      )

for noise in [1e-2, 5e-3, 1e-3, 1e-4]
  summary, callbacks = repeat_experiment(problem_oscillation,
                                        net_config,
                                        repetitions;
                                        ε = noise,
                                        progress=false)
  filename_cb = joinpath("./data", experiment_name*"_oscillation_"*string(noise)*"_all.h5")
  filename_sum = joinpath("./data", experiment_name*"_oscillation_"*string(noise)*".jld2")
  callbacks_to_hdf5(callbacks, filename_cb)
  save(filename_sum, summary)
end

problem_steadystate = Dict([:equation => selkov,
                           :parameters => Float32[0.1, 0.15],
                           :u0 => Float32[1.0, 1.0],
                           :tspan => (0.0f0, 15.0f0),
                           :ts => 0.1,
                           :solver => Tsit5,
                           :optimizer => construct_optimiser,
                           :max_iter => 3000,
                           :loss => construct_loss]
                         )

for noise in [1e-2, 5e-3, 1e-3, 1e-4]
  summary, callbacks = repeat_experiment(problem_steadystate,
                                        net_config,
                                        repetitions;
                                        ε = noise,
                                        progress=false)
  filename_cb = joinpath("./data", experiment_name*"_steadystate_"*string(noise)*"_all.h5")
  filename_sum = joinpath("./data", experiment_name*"_steadystate_"*string(noise)*".jld2")
  callbacks_to_hdf5(callbacks, filename_cb)
  save(filename_sum, summary)
end

max_stars_repo_name: pnieters/GeneralizedDynamicsFromData

max_issues_repo_name: pnieters/GeneralizedDynamicsFromData

max_forks_repo_name: pnieters/GeneralizedDynamicsFromData


text:
using Random
import Base: length, size, iterate, eltype, IteratorSize, IteratorEltype, haslength, @propagate_inbounds, repeat, rand, tail
import .Iterators: cycle, Cycle

mutable struct Data{T}; x; y; batchsize; length; partial; imax; indices; shuffle; xsize; ysize; xtype; ytype; end

"""
    minibatch(x, [y], batchsize; shuffle, partial, xtype, ytype, xsize, ysize)

Return an iterator of minibatches [(xi,yi)...] given data tensors x, y and batchsize.

The last dimension of x and y give the number of instances and should be equal. `y` is
optional, if omitted a sequence of `xi` will be generated rather than `(xi,yi)` tuples.  Use
`repeat(d,n)` for multiple epochs, `Iterators.take(d,n)` for a partial epoch, and
`Iterators.cycle(d)` to cycle through the data forever (this can be used with `converge`).
If you need the iterator to continue from its last position when stopped early (e.g. by a
break in a for loop), use `Iterators.Stateful(d)` (by default the iterator would restart
from the beginning).

Keyword arguments:

- `shuffle=false`: Shuffle the instances every epoch.
- `partial=false`: If true include the last partial minibatch < batchsize.
- `xtype=typeof(x)`: Convert xi in minibatches to this type.
- `ytype=typeof(y)`: Convert yi in minibatches to this type.
- `xsize=size(x)`: Convert xi in minibatches to this shape.
- `ysize=size(y)`: Convert yi in minibatches to this shape.
"""
minibatch, Data

function minibatch(x,y,batchsize; shuffle=false,partial=false,xtype=typeof(x),ytype=typeof(y),xsize=size(x), ysize=size(y))
    nx = size(x)[end]
    if nx != size(y)[end]; throw(DimensionMismatch()); end
    x2 = reshape(x, :, nx)
    y2 = reshape(y, :, nx)
    imax = partial ? nx : nx - batchsize + 1
    # xtype,ytype may be underspecified, here we infer the exact types from the first batch:
    ids = 1:min(nx,batchsize)
    xt = typeof(convert(xtype, reshape(x2[:,ids],xsize[1:end-1]...,length(ids))))
    yt = typeof(convert(ytype, reshape(y2[:,ids],ysize[1:end-1]...,length(ids))))
    Data{Tuple{xt,yt}}(x2,y2,batchsize,nx,partial,imax,1:nx,shuffle,xsize,ysize,xtype,ytype)
end

function minibatch(x,batchsize; shuffle=false,partial=false,xtype=typeof(x),xsize=size(x))
    nx = size(x)[end]
    x2 = reshape(x, :, nx)
    imax = partial ? nx : nx - batchsize + 1
    # xtype may be underspecified, here we infer the exact types from the first batch:
    ids = 1:min(nx,batchsize)
    xt = typeof(convert(xtype, reshape(x2[:,ids],xsize[1:end-1]...,length(ids))))
    Data{xt}(x2,nothing,batchsize,nx,partial,imax,1:nx,shuffle,xsize,nothing,xtype,nothing)
end

@propagate_inbounds function iterate(d::Data, i=0)     # returns data in d.indices[i+1:i+batchsize]
    if i >= d.imax
        return nothing
    end
    if d.shuffle && i == 0
        d.indices = randperm(d.length)
    end
    nexti = min(i + d.batchsize, d.length)
    ids = d.indices[i+1:nexti]
    xbatch = convert(d.xtype, reshape(d.x[:,ids],d.xsize[1:end-1]...,length(ids)))
    if d.y == nothing
        return (xbatch,nexti)
    else
        ybatch = convert(d.ytype, reshape(d.y[:,ids],d.ysize[1:end-1]...,length(ids)))
        return ((xbatch,ybatch),nexti)
    end
end

eltype(::Type{Data{T}}) where T = T

function length(d::Data)
    n = d.length / d.batchsize
    d.partial ? ceil(Int,n) : floor(Int,n)
end

function rand(d::Data)
    i = rand(0:(d.length-d.batchsize))
    return iterate(d, i)[1]
end

# IterTools.ncycle(data,n) for multiple epochs
# Base.Iterators.cycle(data) to go forever
# Base.Iterators.take(data,n) for partial epochs
# IterTools.takenth(itr,n) to report every n iterations

struct Repeat; data::Data; n::Int; end

function repeat(d::Data, n::Int)
    @warn "repeat(d::Data,n) is deprecated, use IterTools.ncycle instead." maxlog=1
    @assert n >= 0
    Repeat(d,n)
end

length(r::Repeat) = r.n * length(r.data)
eltype(r::Repeat) = eltype(r.data)
eltype(c::Cycle{Data}) = eltype(c.xs)
eltype(c::Cycle{Repeat}) = eltype(c.xs)

@propagate_inbounds function iterate(r::Repeat, s=(1,))
    epoch, state = s[1], tail(s)
    epoch > r.n && return nothing
    next = iterate(r.data, state...)
    next === nothing && return iterate(r, (epoch+1,))
    (next[1], (epoch, next[2]))
end

# Give length info in summary:
Base.summary(d::Data) = "$(length(d))-element $(typeof(d))"
Base.show(io::IO, d::Data) = print(IOContext(io,:compact=>true), summary(d))
Base.show(io::IO, ::MIME"text/plain", d::Data) = show(io, d)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

# flip arguments so that points always come first
evaluate(d::PreMetric, g::Geometry, p::Point) = evaluate(d, p, g)

"""
    evaluate(Euclidean(), point, line)

Evaluate the Euclidean distance between `point` and `line`.
"""
function evaluate(::Euclidean, p::Point, l::Line)
  a, b = points(l)
  u = p - a
  v = b - a
  α = (u ⋅ v) / (v ⋅ v)
  norm(u - α*v)
end

"""
    evaluate(::PreMetric, point1, point2)

Evaluate pre-metric between coordinates of `point2` and `point2`.
"""
evaluate(d::PreMetric, p1::Point, p2::Point) =
  evaluate(d, coordinates(p1), coordinates(p2))


text:
"""
Uniformed Matrix for CTMC
"""

macro unif(Q, ufact)
    expr = quote
        qv = maximum(abs.(spdiag($Q))) * $ufact
        if iszero(qv)
            qv = 1.0e-12
        end
        P = $Q / qv
        d = spdiag(P)
        d .+= 1
        (P, qv)
    end
    esc(expr)
end

"""
unif(Q::AbstractSparseM{Tv,Ti}, ufact::Tv = 1.01)
unif(Q::Matrix{Tv}, ufact::Tv = 1.01)

Get an uniformed transition probability matrix from a CTMC kernel.

   P = I + Q / qv
   qv = max(abs(diag(Q))) * ufact

Parameters:
- Q: CTMC Kernel
- ufact: uniformization factor
Return value:
A tuple of
- P: The uniformed transition probability matrix
- qv: The maximum event rate

"""

function unif(Q::SparseMatrixCSC{Tv,Ti}, ufact::Tv = 1.01) where {Tv, Ti}
    @unif(Q, ufact)
end

function unif(Q::SparseCSR{Tv,Ti}, ufact::Tv = 1.01) where {Tv, Ti}
    @unif(Q, ufact)
end

function unif(Q::SparseCSC{Tv,Ti}, ufact::Tv = 1.01) where {Tv, Ti}
    @unif(Q, ufact)
end

function unif(Q::SparseCOO{Tv,Ti}, ufact::Tv = 1.01) where {Tv, Ti}
    @unif(Q, ufact)
end

function unif(Q::Matrix{Tv}, ufact::Tv = 1.01) where {Tv, Ti}
    @unif(Q, ufact)
end


text:
################################################################################
#
#   (8,3)b
#
################################################################################

# Implementation
# - implementation 1
# - labels <: Any
function getUnitcell_8_3_b(
			unitcell_type  :: Type{U},
			implementation :: Val{1}
		) :: U where {LS,LB,S<:AbstractSite{LS,3},B<:AbstractBond{LB,3},U<:AbstractUnitcell{S,B}}

	# return a new Unitcell
	return newUnitcell(
		# Type of the unitcell
		U,
		# lattice vectors
		Vector{Float64}[
			Float64[1/2., 1/(2*sqrt(3)), sqrt(2)/(5*sqrt(3))],
			Float64[0, 1/sqrt(3), (2*sqrt(2))/(5*sqrt(3))],
			Float64[0, 0, sqrt(6)/5],
		],
		# sites
		S[
			newSite(S, Float64[1/10., 1/(2*sqrt(3)), sqrt(2)/(5*sqrt(3))], getDefaultLabelN(LS,1)),
			newSite(S, Float64[1/5., sqrt(3)/5, sqrt(6)/5], getDefaultLabelN(LS,2)),
			newSite(S, Float64[3/10., 11/(10*sqrt(3)), (4*sqrt(2))/(5*sqrt(3))], getDefaultLabelN(LS,3)),
			newSite(S, Float64[1/5., 2/(5*sqrt(3)), (2*sqrt(2))/(5*sqrt(3))], getDefaultLabelN(LS,4)),
			newSite(S, Float64[3/10., (3*sqrt(3))/10., sqrt(6)/5], getDefaultLabelN(LS,5)),
			newSite(S, Float64[2/5., 1/sqrt(3), sqrt(2)/sqrt(3)], getDefaultLabelN(LS,6)),
		],
		# bonds
		B[
			newBond(B, 1, 4, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 4, 2, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 2, 5, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 5, 3, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 3, 6, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 4, 1, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 2, 4, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 5, 2, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 3, 5, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 6, 3, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 6, 1, getDefaultLabel(LB), (1, 0, 1)),
			newBond(B, 1, 6, getDefaultLabel(LB), (-1, 0, -1)),
			newBond(B, 4, 3, getDefaultLabel(LB), (0, -1, 0)),
			newBond(B, 3, 4, getDefaultLabel(LB), (0, 1, 0)),
			newBond(B, 1, 2, getDefaultLabel(LB), (0, 0, -1)),
			newBond(B, 2, 1, getDefaultLabel(LB), (0, 0, 1)),
			newBond(B, 5, 6, getDefaultLabel(LB), (0, 0, -1)),
			newBond(B, 6, 5, getDefaultLabel(LB), (0, 0, 1)),
		]
	)
end

# Implementation
# - implementation 4 (Kitaev)
# - labels <: Any
function getUnitcell_8_3_b(
			unitcell_type  :: Type{U},
			implementation :: Val{4}
		) :: U where {LS,LB,S<:AbstractSite{LS,3},B<:AbstractBond{LB,3},U<:AbstractUnitcell{S,B}}

	# return a new Unitcell
	return newUnitcell(
		# Type of the unitcell
		U,
		# lattice vectors
		Vector{Float64}[
			Float64[1/2., 1/(2*sqrt(3)), sqrt(2)/(5*sqrt(3))],
			Float64[0, 1/sqrt(3), (2*sqrt(2))/(5*sqrt(3))],
			Float64[0, 0, sqrt(6)/5],
		],
		# sites
		S[
			newSite(S, Float64[1/10., 1/(2*sqrt(3)), sqrt(2)/(5*sqrt(3))], getDefaultLabelN(LS,1)),
			newSite(S, Float64[1/5., sqrt(3)/5, sqrt(6)/5], getDefaultLabelN(LS,2)),
			newSite(S, Float64[3/10., 11/(10*sqrt(3)), (4*sqrt(2))/(5*sqrt(3))], getDefaultLabelN(LS,3)),
			newSite(S, Float64[1/5., 2/(5*sqrt(3)), (2*sqrt(2))/(5*sqrt(3))], getDefaultLabelN(LS,4)),
			newSite(S, Float64[3/10., (3*sqrt(3))/10., sqrt(6)/5], getDefaultLabelN(LS,5)),
			newSite(S, Float64[2/5., 1/sqrt(3), sqrt(2)/sqrt(3)], getDefaultLabelN(LS,6)),
		],
		# bonds
		B[
			newBond(B, 1, 4, getDefaultLabelZ(LB), (0, 0, 0)),
			newBond(B, 4, 2, getDefaultLabelY(LB), (0, 0, 0)),
			newBond(B, 2, 5, getDefaultLabelZ(LB), (0, 0, 0)),
			newBond(B, 5, 3, getDefaultLabelY(LB), (0, 0, 0)),
			newBond(B, 3, 6, getDefaultLabelZ(LB), (0, 0, 0)),
			newBond(B, 4, 1, getDefaultLabelZ(LB), (0, 0, 0)),
			newBond(B, 2, 4, getDefaultLabelY(LB), (0, 0, 0)),
			newBond(B, 5, 2, getDefaultLabelZ(LB), (0, 0, 0)),
			newBond(B, 3, 5, getDefaultLabelY(LB), (0, 0, 0)),
			newBond(B, 6, 3, getDefaultLabelZ(LB), (0, 0, 0)),
			newBond(B, 6, 1, getDefaultLabelY(LB), (1, 0, 1)),
			newBond(B, 1, 6, getDefaultLabelY(LB), (-1, 0, -1)),
			newBond(B, 4, 3, getDefaultLabelX(LB), (0, -1, 0)),
			newBond(B, 3, 4, getDefaultLabelX(LB), (0, 1, 0)),
			newBond(B, 1, 2, getDefaultLabelX(LB), (0, 0, -1)),
			newBond(B, 2, 1, getDefaultLabelX(LB), (0, 0, 1)),
			newBond(B, 5, 6, getDefaultLabelX(LB), (0, 0, -1)),
			newBond(B, 6, 5, getDefaultLabelX(LB), (0, 0, 1)),
		]
	)
end


text:
module TestKDTree
	using Base.Test
	using Distance
	using NearestNeighbors

	# Empty tree
	t = KDTree()
	@test typeof(t.root) == NearestNeighbors.EmptyKDTree

	t = KDTree(Euclidean())
	@test typeof(t.root) <: NearestNeighbors.EmptyKDTree
	@test typeof(t.metric) == Euclidean

	# With one element
	x = [1, 2, 3]
	t = KDTree(x, 1, Euclidean())
	@test typeof(t.root) <: NearestNeighbors.KDTreeNode
	@test t.root.k == x
	@test t.root.i == 1

	# With multiple points
	X = [1 60 29 7 86 44 23 54 12]
	t = KDTree(X)
	@test typeof(t.root) <: NearestNeighbors.KDTreeNode

	# Search
	#Distance.evaluate(m::Metric, a::Int, b::Float64) = evaluate(m, float(a) ,b)
	ind, dist = nearest(t, [3], 2)
	@test length(ind) == 2
	@test length(dist) ==  2
	@test 4 in ind && 1 in ind
	@test_approx_eq(4.0, dist[1])
	@test_approx_eq(2.0, dist[2])

	X = readdlm(Pkg.dir("NearestNeighbors", "test", "iris.csv"), ',')
	t = KDTree(X)
	v = X[:, 84]

	@test sort(nearest(t, v, 1)[1])[1] == 84
	@test sort(nearest(t, v, 2)[1])[1] == 84
	@test sort(nearest(t, v, 3)[1])[1] == 84
end


text:
## Occupation vector indexing.

Base.getindex(basis::AbstractSzbasis, i::Int) = @view basis.vectors[:, i]


text:
@testset "jitter.jl" begin
    println("             > jitter.jl")
    m_out = 5
    model = BinomialGridModel(
        m_out,
        1:5,
        LinRange(0.05,0.95,5),
        LinRange(0.1,2,5),
        LinRange(0.05,2,5),
        LinRange(0.05,2,5)
    )

    for i in 1:100
        model_old = deepcopy(model)
        jitter!(model, 12)

        # check whether new indices are in the allowed range of indices
        @test all(1 .<= model.Nind .<= length(model.Nrng))
        @test all(1 .<= model.pind .<= length(model.prng))
        @test all(1 .<= model.qind .<= length(model.qrng))
        @test all(1 .<= model.σind .<= length(model.σrng))
        @test all(1 .<= model.τind .<= length(model.τrng))

        # check whether new indices are different by at most one from old indices
        @test all(abs.(model.Nind .- model_old.Nind) .<= 1)
        @test all(abs.(model.pind .- model_old.pind) .<= 1)
        @test all(abs.(model.qind .- model_old.qind) .<= 1)
        @test all(abs.(model.σind .- model_old.σind) .<= 1)
        @test all(abs.(model.τind .- model_old.τind) .<= 1)
    end

    m_out = 1024
    model = BinomialGridModel(
        m_out,
        1:5,
        LinRange(0.05,0.95,5),
        LinRange(0.1,2,5),
        LinRange(0.05,2,5),
        LinRange(0.05,2,5)
    )

    if RUN_BENCHMARKS
        println("")
        println("Benchmarking function jitter!: should take about 80μs")
        display(@benchmark CUDA.@sync jitter!($model, 12))
        println("")
    end
end


text:
using Distributions, Test

@testset "truncated uniform" begin
    # just test equivalence of truncation results
    u = Uniform(1, 2)
    @test truncated(u) === u
    @test truncated(u; lower=0) == u
    @test truncated(u; lower=0, upper=Inf) == u
    @test truncated(u; upper=2.1) == u
    @test truncated(u; lower=-Inf, upper=2.1) == u
    @test truncated(u; lower=1.1) == Uniform(1.1, 2)
    @test truncated(u; lower=1.1, upper=Inf) == Uniform(1.1, 2)
    @test truncated(u, 1.1, 2.1) == Uniform(1.1, 2)
    @test truncated(u, 1.1, 1.9) == Uniform(1.1, 1.9)
end


text:
####################################################

#  This file is used for development and testing of PRM.jl.
#  It runs the PRM once and plots the results.
#  nouyang 2017

####################################################
# MAIN
####################################################
include("PRM.jl")

function main()
    numSamples = 25
    connectRadius =10
    param = algT.AlgParameters(numSamples, connectRadius)

    print("\n ---- Running PRM ------ \n")
    ## Define obstacles
    obs1 = HyperRectangle(Vec(8,3.), Vec(2,2.)) #Todo
    obs2 = HyperRectangle(Vec(4,4.), Vec(2,10.)) #Todo

    obstacles = Vector{HyperRectangle}()
    push!(obstacles, obs1, obs2)

    # Define walls
    walls = Vector{LineSegment}()
    w,h  = 20,20
    perimeter = HyperRectangle(Vec(0.,0), Vec(w,h))
    roomPerimeter = algfxn.decompRect(perimeter)
    #internalwalls =  (linesegs
    for l in roomPerimeter
        push!(walls, l)
    end

    r = algT.Room(w,h,walls,obstacles)
    plotfxn.plotRoom(r)

    ## Run preprocessing
    nodeslist, edgeslist = preprocessPRM(r, param)

    ## Query created RM
    startstate = Point(1.,1)
    goalstate = Point(18.,18)

    pathcost, isPathFound, solPath = queryPRM(startstate, goalstate, nodeslist, edgeslist, obstacles)

    ## Plot path found
    title = "PRM with # samples =$numSamples, \nPathfound = $isPathFound, \npathcost = $pathcost"
    roadmap = algT.roadmap(startstate, goalstate, nodeslist, edgeslist)

    plot = plotfxn.plotPRM(roadmap, solPath, title::String)

####
#startGoal = algT.GraphNode(0, Point(0,0))
#endNode = algT.GraphNode(0, Point(0,0))

end

function collTest()

####
    function decompRect(r::HyperRectangle) #GeometryTypes.HyperRectangle{2,Float64}
        corners = decompose(Point{2, Float64}, r)
        corners = [Point(pt) for pt in corners]
        lineBottom  = LineSegment(corners[1], corners[2])
        lineTop     = LineSegment(corners[3], corners[4])
        lineLeft    = LineSegment(corners[1], corners[3])
        lineRight   = LineSegment(corners[2], corners[4])
        lines = Vector{LineSegment}()
        push!(lines, lineTop, lineRight, lineBottom, lineLeft)
        return lines
    end

    function isCollidingEdge(line::LineSegment, obsList::Vector{HyperRectangle})

        print("\n -------- \n")
		@show line
        print("\n -------- \n")
        for obs in obsList
            rectLines = decompRect(obs)
            for rectline in rectLines
                @show intersects(line, rectline)[1]
				@show rectline
        print("\n -------- \n")
                if intersects(line, rectline)[1]
                    return true
                end
            end
        end
        return false
    end

function ccw(A,B,C)
    # determines direction of lines formed by three points
	return (C[2]-A[2]) * (B[1]-A[1]) > (B[2]-A[2]) * (C[1]-A[1])
	#return (C.y-A.y) * (B.x-A.x) > (B.y-A.y) * (C.x-A.x)
end

function intersectLineSeg(line1, line2) #no ":" at the end!
	A, B = line1[1], line1[2]
	C, D = line2[1], line2[2]
 	return ( (ccw(A, C, D) != ccw(B, C, D)) && ccw(A, B, C) != ccw(A, B, D))
end
####
	#obs1 = GeometryTypes.HyperRectangle{2,Float64}([7.67199, 10.8904], [1.67166, 2.51185])
	obs1 = GeometryTypes.HyperRectangle{2,Float64}([7,10], [2,3])

    obstacles = Vector{HyperRectangle}()
    push!(obstacles, obs1)

	#nodeslist = algT.GraphNode[algT.GraphNode(1, [18.7716, 16.8013]), algT.GraphNode(2, [12.8554, 11.3548]), algT.GraphNode(3, [13.4626, 13.6564]), algT.GraphNode(4, [14.5381, 7.52624]), algT.GraphNode(5, [15.5076, 17.8328]), algT.GraphNode(6, [4.65531, 14.1306]), algT.GraphNode(7, [15.3949, 17.6171]), algT.GraphNode(8, [1.34579, 4.52833]), algT.GraphNode(9, [15.2216, 11.153]), algT.GraphNode(10, [6.29213, 11.7306])]

    #anEdge = algT.Edge(2, 10)

	#startnode = nodeslist[2]
	startnode = algT.GraphNode(1, [6,11])
	#endnode = nodeslist[10]
	endnode = algT.GraphNode(2, [12,11])

	l = LineSegment(Point(6.0,11), Point(12.0,11))

    coll = isCollidingEdge(l, obstacles)
    @show coll
    #isCollidingNode

#(intersects(line, rectline))[1] = true
#r1 = GeometryTypes.Point{2,Float64}[[9.0, 10.0], [9.0, 13.0]]
r1 =LineSegment( Point(9.0,10.0), Point(9.0, 13.0))

 #--------
#(intersects(line, rectline))[1] = false
#r2 = GeometryTypes.Point{2,Float64}[[9.34365, 10.8904], [9.34365, 13.4022]]
r2 =LineSegment( Point(9.34365, 10.890), Point(9.34365, 13.402))

@show f= intersects(l, r1)
@show intersects(l, r2)
@show intersectLineSeg(l, r1)
@show intersectLineSeg(l, r2)

#https://github.com/JuliaGeometry/GeometryTypes.jl/blob/master/src/lines.jl

end

function rectUnionAreaTest()
    #obs1 = HyperRectangle(Vec(8,3.), Vec(4,4.)) #Todo
    #obs2 = HyperRectangle(Vec(10,5.), Vec(4,4.)) #Todo
    c = HyperRectangle(Vec(0,0),Vec(2,2))
    d = HyperRectangle(Vec(1,1),Vec(2,2))

    function areaRect(r::HyperRectangle)
        w, h = widths(r)
        area = w*h
    end

    function unionAreaRects(r1::HyperRectangle, r2::HyperRectangle)
        area1 = areaRect(r1)
        area2 = areaRect(r2)
        overlap = areaRect(intersect(r1, r2))
        area = area1 + area2 - overlap
        return area
    end

    bool = intersect(c,d) == HyperRectangle(Vec(1,1), Vec(1,1))
	@show bool

    areaRect(c)
    areaRect(d)
    unionAreaRects(c,d)

end

########################################
# Call main() function
########################################

main()
#collTest()
#rectUnionAreaTest()

#
 # EDGE REMOVED: algT.GraphNode(3, [13.4626, 13.6564]) WITH algT.GraphNode(10, [6.29213, 11.7306])
 # EDGE REMOVED: algT.GraphNode(9, [15.2216, 11.153]) WITH algT.GraphNode(10, [6.29213, 11.7306])
 # EDGE REMOVED: algT.GraphNode(10, [6.29213, 11.7306]) WITH algT.GraphNode(3, [13.4626, 13.6564])
 # EDGE REMOVED: algT.GraphNode(10, [6.29213, 11.7306]) WITH algT.GraphNode(9, [15.2216, 11.153])
 # ---- pre-processed ---
 # ---- queried ----
 # ---Solution Path-----
# solPath = algT.GraphNode[algT.GraphNode(0, [1.0, 1.0]), algT.GraphNode(8, [1.34579, 4.52833]), algT.GraphNode(10, [6.29213, 11.7306]), algT.GraphNode(2, [12.8554, 11.3548]), algT.GraphNode(1, [18.7716, 16.8013]), algT.GraphNode(0, [18.0, 18.0])]
 # --------
 # --- Time ---
# timestamp = 2017-10-28T03:19:26.972
 # --- Obstacles ---
# obstacles = GeometryTypes.HyperRectangle[GeometryTypes.HyperRectangle{2,Float64}([7.67199, 10.8904], [1.67166, 2.51185]), GeometryTypes.HyperRectangle{2,Float64}([1.27599, 14.9952], [1.25427, 3.5936]), GeometryTypes.HyperRectangle{2,Float64}([4.64947, 3.96388], [8.07466, 4.21098])]
 # --- Nodes ---
# nodeslist = algT.GraphNode[algT.GraphNode(1, [18.7716, 16.8013]), algT.GraphNode(2, [12.8554, 11.3548]), algT.GraphNode(3, [13.4626, 13.6564]), algT.GraphNode(4, [14.5381, 7.52624]), algT.GraphNode(5, [15.5076, 17.8328]), algT.GraphNode(6, [4.65531, 14.1306]), algT.GraphNode(7, [15.3949, 17.6171]), algT.GraphNode(8, [1.34579, 4.52833]), algT.GraphNode(9, [15.2216, 11.153]), algT.GraphNode(10, [6.29213, 11.7306])]
#
 # --- Edges ---
# edgeslist = algT.Edge[algT.Edge(1, 1), algT.Edge(1, 2), algT.Edge(1, 3), algT.Edge(1, 5), algT.Edge(1, 7), algT.Edge(1, 9), algT.Edge(2, 1), algT.Edge(2, 2), algT.Edge(2, 3), algT.Edge(2, 4), algT.Edge(2, 5), algT.Edge(2, 6), algT.Edge(2, 7), algT.Edge(2, 9), algT.Edge(2, 10), algT.Edge(3, 1), algT.Edge(3, 2), algT.Edge(3, 3), algT.Edge(3, 4), algT.Edge(3, 5), algT.Edge(3, 6), algT.Edge(3, 7), algT.Edge(3, 9), algT.Edge(4, 2), algT.Edge(4, 3), algT.Edge(4, 4), algT.Edge(4, 9), algT.Edge(4, 10), algT.Edge(5, 1), algT.Edge(5, 2), algT.Edge(5, 3), algT.Edge(5, 5), algT.Edge(5, 7), algT.Edge(5, 9), algT.Edge(6, 2), algT.Edge(6, 3), algT.Edge(6, 6), algT.Edge(6, 10), algT.Edge(7, 1), algT.Edge(7, 2), algT.Edge(7, 3), algT.Edge(7, 5), algT.Edge(7, 7), algT.Edge(7, 9), algT.Edge(8, 8), algT.Edge(8, 10), algT.Edge(9, 1), algT.Edge(9, 2), algT.Edge(9, 3), algT.Edge(9, 4), algT.Edge(9, 5), algT.Edge(9, 7), algT.Edge(9, 9), algT.Edge(10, 2), algT.Edge(10, 4), algT.Edge(10, 6), algT.Edge(10, 8), algT.Edge(10, 10)]


text:
module DiffEqBayes
using DiffEqBase, Distributions, Turing, MacroTools, Mamba
using ParameterizedFunctions, RecursiveArrayTools
using Parameters, Distributions, Optim, Requires
using Distances, ApproxBayes, DocStringExtensions, Random

STANDARD_PROB_GENERATOR(prob,p) = remake(prob;u0=eltype(p).(prob.u0),p=p)
STANDARD_PROB_GENERATOR(prob::EnsembleProblem,p) = EnsembleProblem(remake(prob.prob;u0=eltype(p).(prob.prob.u0),p=p))

include("turing_inference.jl")
include("abc_inference.jl")

function __init__()
    @require CmdStan="593b3428-ca2f-500c-ae53-031589ec8ddd" begin
        using CmdStan
        include("stan_inference.jl")
        include("stan_string.jl")
        export stan_inference, stan_string
    end

    @require DynamicHMC="bbc10e6e-7c05-544b-b16e-64fede858acb" begin
        using DynamicHMC, TransformVariables, LogDensityProblems
        include("dynamichmc_inference.jl")
        export dynamichmc_inference
    end

    @require StatsPlots="f3b207a7-027a-5e70-b257-86293d7955fd" begin
        include("utils.jl")
        export plot_chain
    end
end

export turing_inference, abc_inference

end # module


text:
using ADMIXTURE

println("running chr18.uniqueSNPs.bed with 26 populations and 10 threads")

P, Q = admixture("chr18.uniqueSNPs.bed", 26, j=10)

max_stars_repo_name: Hua-Zhou/MendelImpute


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

"""
    TransformedDensity <: AbstractDensity

Abstract type for transformed densities.

In addition to the [`AbstractDensity`](@ref) interface, subtypes must support
the functions `Base.parent` and [`trafoof`](@ref):

```julia
parent(d::SomeTransformedDensity)::AbstractDensity
trafoof(d::SomeTransformedDensity)::AbstractVariateTransform
```
"""
struct AbstractTransformedDensity end
export AbstractTransformedDensity

"""
    trafoof(d::AbstractTransformedDensity)::AbstractDensity

Get the transform from `parent(d)` to `d`, so that

```julia
trafoof(d)(parent(d)) == d
```
"""
function trafoof end
export trafoof

abstract type TDVolCorr end
struct TDNoCorr <: TDVolCorr end
struct TDLADJCorr <: TDVolCorr end

"""
    TransformedDensity

*BAT-internal, not part of stable public API.*
"""
struct TransformedDensity{D<:AbstractDensity,FT<:VariateTransform,VC<:TDVolCorr} <: AbstractDensity
    orig::D
    trafo::FT
    volcorr::VC

    # ToDo: Check varshape(orig) == varshape(trafo) in ctor
end

@inline function (trafo::VariateTransform)(density::AbstractDensity; volcorr::Val{vc} = Val(true)) where vc
    if vc
        TransformedDensity(density, trafo, TDLADJCorr())
    else
        TransformedDensity(density, trafo, TDNoCorr())
    end
end

Base.parent(density::TransformedDensity) = density.orig
trafoof(density::TransformedDensity) = density.trafo

# var_bounds(density::TransformedDensity) = ...

ValueShapes.varshape(density::TransformedDensity{<:AbstractDensity,<:DistributionTransform}) = varshape(density.trafo.target_dist)

# ToDo: Should not be neccessary, improve default implementation of
# ValueShapes.totalndof(density::AbstractDensity):
ValueShapes.totalndof(density::TransformedDensity) = totalndof(varshape(density))

var_bounds(density::TransformedDensity) = _trafo_var_bounds(density.trafo)

_trafo_var_bounds(trafo::VariateTransform) = missing

function _trafo_var_bounds(trafo::DistributionTransform{<:Any,<:Any,<:Union{StandardUvUniform,StandardMvUniform}})
    n = totalndof(trafo.target_varshape)
    HyperRectBounds(fill(_default_PT(0), n), fill(_default_PT(1), n))
end

function _trafo_var_bounds(trafo::DistributionTransform{<:Any,<:Any,<:Union{StandardUvNormal,StandardMvNormal}})
    n = totalndof(trafo.target_varshape)
    HyperRectBounds(fill(_default_PT(-Inf), n), fill(_default_PT(+Inf), n))
end

function eval_logval(density::TransformedDensity{D,FT,TDNoCorr}, v::Any, T::Type{<:Real}) where {D,FT}
    v_shaped = fixup_variate(varshape(density), v)
    v_orig = inv(density.trafo)(v_shaped)
    eval_logval(parent(density), v_orig, T)
end

function eval_logval(density::TransformedDensity{D,FT,TDLADJCorr}, v::Any, T::Type{<:Real}) where {D,FT,}
    v_shaped = fixup_variate(varshape(density), v)
    R = density_logval_type(v_shaped, T)
    r = inv(density.trafo)(v_shaped, 0)
    v_orig = r.v
    ldaj = r.ladj
    logd_orig = eval_logval(parent(density), v_orig, R)

    logd_result = logd_orig + ldaj
    R = typeof(logd_result)

    if isnan(logd_result) && logd_orig == -Inf && ldaj == +Inf
        # Zero density wins against infinite volume:
        R(-Inf)
    elseif isfinite(logd_orig) && (ldaj == -Inf)
        # Maybe  also for (logd_orig == -Inf) && isfinite(ldaj) ?
        # Return constant -Inf to prevent problems with ForwardDiff:
        #R(-Inf)
        near_neg_inf(R) # Avoids AdvancedHMC warnings
    else
        logd_result
    end
end

function eval_logval_unchecked(density::TransformedDensity, v::Any)
    eval_logval(density, v, default_dlt())
end


text:
function ChiSquareTest()
    return ChiSquareTest(())
end

function chi_square(obj::ChiSquareTest, arg0::Vector{Vector{jlong}})
    return jcall(obj, "chiSquare", jdouble, (Vector{Vector{jlong}},), arg0)
end

function chi_square(obj::ChiSquareTest, arg0::Vector{jdouble}, arg1::Vector{jlong})
    return jcall(obj, "chiSquare", jdouble, (Vector{jdouble}, Vector{jlong}), arg0, arg1)
end

function chi_square_data_sets_comparison(obj::ChiSquareTest, arg0::Vector{jlong}, arg1::Vector{jlong})
    return jcall(obj, "chiSquareDataSetsComparison", jdouble, (Vector{jlong}, Vector{jlong}), arg0, arg1)
end

function chi_square_test(obj::ChiSquareTest, arg0::Vector{Vector{jlong}})
    return jcall(obj, "chiSquareTest", jdouble, (Vector{Vector{jlong}},), arg0)
end

function chi_square_test(obj::ChiSquareTest, arg0::Vector{Vector{jlong}}, arg1::jdouble)
    return jcall(obj, "chiSquareTest", jboolean, (Vector{Vector{jlong}}, jdouble), arg0, arg1)
end

function chi_square_test(obj::ChiSquareTest, arg0::Vector{jdouble}, arg1::Vector{jlong})
    return jcall(obj, "chiSquareTest", jdouble, (Vector{jdouble}, Vector{jlong}), arg0, arg1)
end

function chi_square_test(obj::ChiSquareTest, arg0::Vector{jdouble}, arg1::Vector{jlong}, arg2::jdouble)
    return jcall(obj, "chiSquareTest", jboolean, (Vector{jdouble}, Vector{jlong}, jdouble), arg0, arg1, arg2)
end

function chi_square_test_data_sets_comparison(obj::ChiSquareTest, arg0::Vector{jlong}, arg1::Vector{jlong})
    return jcall(obj, "chiSquareTestDataSetsComparison", jdouble, (Vector{jlong}, Vector{jlong}), arg0, arg1)
end

function chi_square_test_data_sets_comparison(obj::ChiSquareTest, arg0::Vector{jlong}, arg1::Vector{jlong}, arg2::jdouble)
    return jcall(obj, "chiSquareTestDataSetsComparison", jboolean, (Vector{jlong}, Vector{jlong}, jdouble), arg0, arg1, arg2)
end


text:
# run examples
const EXAMPLES_SRC = joinpath(@__DIR__, "..", "examples")
const EXAMPLES_OUT = joinpath(@__DIR__, "src", "examples")
const LITERATEJL = joinpath(@__DIR__, "literate.jl")

# Always rerun examples
ispath(EXAMPLES_OUT) && rm(EXAMPLES_OUT; recursive=true)
mkpath(EXAMPLES_OUT)

# Run examples asynchronously
processes = map(filter!(isdir, readdir(EXAMPLES_SRC; join=true))) do example
    scriptjl = joinpath(example, "script.jl")
    return run(
        pipeline(
            `$(Base.julia_cmd()) $LITERATEJL $scriptjl $EXAMPLES_OUT`;
            stdin=devnull,
            stdout=devnull,
            stderr=stderr,
        );
        wait=false,
    )::Base.Process
end

# Check that all examples were run successfully
isempty(processes) || success(processes) || error("some examples were not run successfully")

# Build documentation
using KernelFunctions
using Documenter

# Print `@debug` statements (https://github.com/JuliaDocs/Documenter.jl/issues/955)
if haskey(ENV, "GITHUB_ACTIONS")
    ENV["JULIA_DEBUG"] = "Documenter"
end

# Doctest setup
DocMeta.setdocmeta!(
    KernelFunctions,
    :DocTestSetup,
    quote
        using KernelFunctions
        using LinearAlgebra
        using Random
        using PDMats: PDMats
    end;
    recursive=true,
)

makedocs(;
    sitename="KernelFunctions.jl",
    format=Documenter.HTML(),
    modules=[KernelFunctions],
    pages=[
        "Home" => "index.md",
        "userguide.md",
        "kernels.md",
        "transform.md",
        "metrics.md",
        "create_kernel.md",
        "API" => "api.md",
        "Design" => "design.md",
        "Examples" =>
            map(filter!(filename -> endswith(filename, ".md"), readdir(EXAMPLES_OUT))) do x
                return joinpath("examples", x)
            end,
    ],
    strict=true,
    checkdocs=:exports,
    doctestfilters=[
        r"{([a-zA-Z0-9]+,\s?)+[a-zA-Z0-9]+}",
        r"(Array{[a-zA-Z0-9]+,\s?1}|Vector{[a-zA-Z0-9]+})",
        r"(Array{[a-zA-Z0-9]+,\s?2}|Matrix{[a-zA-Z0-9]+})",
    ],
)

deploydocs(;
    repo="github.com/JuliaGaussianProcesses/KernelFunctions.jl.git", push_preview=true
)


text:
using Documenter, CPIMC

makedocs(
    modules = [CPIMC],
    sitename = "CPIMC.jl",
    authors = "Kai Hunger, Arif Yilmaz, Paul Hamann",
    pages = [
        "Home" => "index.md",
        "Manual" => "manual.md",
        "API reference" => "api.md"
            ])

deploydocs(
    repo = "github.com/CPIMC/CPIMC.jl.git",
)


text:
using LowRankTensorApproximation
using Test

@testset "LowRankTensorApproximation.jl" begin
    # Write your own tests here.
end


text:
module TestExamples2DEulerMulticomponent

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d")

@testset "Compressible Euler Multicomponent" begin
  @testset "elixir_eulermulti_shock_bubble.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_shock_bubble.jl"),
      l2   = [6.03088526e-02, 2.88034228e-03, 2.41548900e+01, 7.83393342e-02, 1.91342787e+04],
      linf = [7.37128632e-01, 1.05064677e-01, 1.99311122e+02, 2.01502716e+00, 1.57477455e+05],
      tspan = (0.0, 0.0001))
  end

  @testset "elixir_eulermulti_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_ec.jl"),
      l2   = [1.23457293e-02, 4.93829171e-02, 5.01948074e-02, 5.02023248e-02, 2.25886833e-01],
      linf = [5.96271450e-02, 2.38508580e-01, 3.06937711e-01, 3.06807092e-01, 1.06295287e+00])
  end

  @testset "elixir_eulermulti_es.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_es.jl"),
      l2   = [1.22783570e-02, 4.91134280e-02, 4.97009915e-02, 4.97029892e-02, 2.24364821e-01],
      linf = [6.33058769e-02, 2.53223508e-01, 2.52097945e-01, 2.48924663e-01, 9.32036428e-01])
  end

  @testset "elixir_eulermulti_eoc_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_eoc_ec.jl"),
      l2   = [3.41443064e-06, 1.36577226e-05, 1.44713621e-05, 1.44713621e-05, 2.28669565e-05],
      linf = [1.81596016e-05, 7.26384063e-05, 7.36760417e-05, 7.36760417e-05, 1.21009804e-04])
  end

  @testset "elixir_eulermulti_eoc_es.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_eoc_es.jl"),
      l2   = [5.54801268e-07, 2.21920507e-06, 1.99274552e-06, 1.99274552e-06, 6.05125237e-06],
      linf = [4.60991759e-06, 1.84396703e-05, 1.83105421e-05, 1.83105421e-05, 5.41312537e-05])
  end
end

end # module


text:
using StaticArrays
using BenchmarkTools
using Statistics

mutable struct Multidual{N, T}
    val::T
    derivs::SVector{N, T}
end

Base.:+(f::Multidual, g::Multidual) = Multidual(f.val + g.val, f.derivs + g.derivs)
Base.:+(f::Multidual, g::Number) = Multidual(f.val + g, f.derivs)
Base.:+(f::Number, g::Multidual) = g + f

Base.:-(f::Multidual, g::Multidual) = Multidual(f.val - g.val, f.derivs - g.derivs)
Base.:-(f::Multidual, g::Number) = Multidual(f.val - g, f.derivs)
Base.:-(f::Number, g::Multidual) = g - f
Base.:-(f::Multidual) = Multidual(-f.val, -1 .* f.derivs)

Base.:*(f::Multidual, g::Multidual) = Multidual(f.val*g.val, f.val .* g.derivs + g.val .* f.derivs)
Base.:*(α::Number, f::Multidual) = Multidual(f.val * α, α .* f.derivs)
Base.:*(f::Multidual, α::Number) = α * f

function f(u, p)
    α, β, γ, δ = p

    du1 = α*u[1] - β*u[1]*u[2]
    du2 = -γ*u[2] + δ*u[1]*u[2]

    @SVector [du1, du2]
end

function dopri5!(u_vals, p, t_span, f, derivs)
    # Values from the Dormand Prince Tableau
    a21  =  (1.0/5.0)
    a31  = (3.0/40.0)
    a32  = (9.0/40.0)
    a41  = (44.0/45.0)
    a42  = (-56.0/15.0)
    a43  = (32.0/9.0)
    a51  = (19372.0/6561.0)
    a52  = (-25360.0/2187.0)
    a53  = (64448.0/6561.0)
    a54  = (-212.0/729.0)
    a61  = (9017.0/3168.0)
    a62  = (-355.0/33.0)
    a63  = (46732.0/5247.0)
    a64  = (49.0/176.0)
    a65  = (-5103.0/18656.0)
    a71  = (35.0/384.0)
    a72  = (0.0)
    a73  = (500.0/1113.0)
    a74  = (125.0/192.0)
    a75  = (-2187.0/6784.0)
    a76  = (11.0/84.0)

    b1   = (35.0/384.0)
    b2   = (0.0)
    b3   = (500.0/1113.0)
    b4   = (125.0/192.0)
    b5   = (-2187.0/6784.0)
    b6   = (11.0/84.0)
    b7   = (0.0)

    dt = 0.25
    uₚ = SVector{2, Float64}(u_vals[1, :])
    @inbounds for i in 1:Int64(t_span/dt)
        k1 = f(uₚ, p)
        k2 = f(uₚ + dt.*(a21.*k1), p)
        k3 = f(uₚ + dt.*(a31.*k1 + a32.*k2), p)
        k4 = f(uₚ + dt.*(a41.*k1 + a42.*k2 + a43.*k3), p)
        k5 = f(uₚ + dt.*(a51.*k1 + a52.*k2 + a53.*k3 + a54.*k4), p)
        k6 = f(uₚ + dt.*(a61.*k1 + a62.*k2 + a63.*k3 + a64.*k4 + a65.*k5), p)
        u = uₚ + dt.*(b1.*k1 + b3.*k3 + b4.*k4 + b5.*k5 + b6.*k6)
        u_vals[i, 1] = u[1].val
        u_vals[i, 2] = u[2].val
        derivs[i, 1:2] = u[1].derivs
        derivs[i, 3:4] = u[2].derivs
        uₚ = SVector{2, Float64}(u_vals[i, :])
    end
end

# Generating the dataset for parameter estimation
tp1 = Multidual(1.5, SVector(0.0, 0.0))
tp2 = Multidual(1.0, SVector(0.0, 0.0))
tp3 = Multidual(3.0, SVector(0.0, 0.0))
tp4 = Multidual(1.0, SVector(0.0, 0.0))

tp = (tp1, tp2, tp3, tp4)

dataset = Array{Float64}(undef, 40, 2)
dataset[1, :] = [1.0, 1.0]
t_derivs = Array{Float64}(undef, 40, 4)

t_span = 10.0
dopri5!(dataset, tp, t_span, f, t_derivs)

function cost_fun(derivs, ds, up, m)
    diff = abs.(ds .- up)
    dcdu = SVector{2, Float64}((1/m) .* sum(diff, dims=1))
    derivs[:, 1:2] = dcdu[1] .* derivs[:, 1:2] # dx/dp
    derivs[:, 3:4] = dcdu[2] .* derivs[:, 3:4] # dy/dp
    dcdp = SVector{4, Float64}(mean(derivs, dims=1))
    cost = SVector{2, Float64}((1/(2*m)) .* sum((diff .^ 2), dims=1))
    return @SVector[cost, dcdp]
end

function estim_params!(p, ds, t_span, f)
    m = 40
    up = Array{Float64}(undef, 40, 2)
    derivs = Array{Float64}(undef, 40, 4)
    cost = SVector{2, Float64}(1.0, 1.0)
    diff = 1.0
    tol = 10e-6
    lr = 0.001
    count = 0

    while !isapprox(diff, 0.0, atol=tol)
        up[1, :] = [1.0, 1.0]
        dopri5!(up, p, t_span, f, derivs)
        prev = cost
        cost, dcdp = cost_fun(derivs, ds, up, m)
        p[1].val = p[1].val - lr*dcdp[1]
        p[2].val = p[2].val - lr*dcdp[2]
        p[3].val = p[3].val - lr*dcdp[3]
        p[4].val = p[4].val - lr*dcdp[4]
        diff = sum(abs.(cost .- prev))
        count += 1
        if count%250 == 0
            println("Training pass $(count), cost = $(cost), $(p[1].val), $(p[2].val), $(p[3].val), $(p[4].val)")
        end
    end
    println("Finished Training.\nThe estimated parameters are: \nα = $(p[1].val)\ttrue value = 1.5\nβ = $(p[2].val)\ttrue value = 1.0\nγ = $(p[3].val)\ttrue value = 3.0\nδ = $(p[4].val)\ttrue value = 1.0\n")
end

α = Multidual(1.2, SVector(1.0, 0.0))
β = Multidual(0.8, SVector(0.0, 1.0))
γ = Multidual(2.8, SVector(1.0, 0.0))
δ = Multidual(0.8, SVector(0.0, 1.0))

p = (α, β, γ, δ)

estim_params!(p, dataset, t_span, f)


text:
using Turing, Random, MacroTools, Distributions, Test
import Turing.translate_tilde!

dir = splitdir(splitdir(pathof(Turing))[1])[1]
include(dir*"/test/test_utils/AllUtils.jl")

Random.seed!(129)

priors = 0 # See "new grammar" test.

@turing_testset "compiler.jl" begin
    @testset "assume" begin
        @model test_assume() = begin
            x ~ Bernoulli(1)
            y ~ Bernoulli(x / 2)
            x, y
        end

        smc = SMC()
        pg = PG(10)

        res1 = sample(test_assume(), smc, 1000)
        res2 = sample(test_assume(), pg, 1000)

        check_numerical(res1, [:y], [0.5], atol=0.1)
        check_numerical(res2, [:y], [0.5], atol=0.1)

        # Check that all xs are 1.
        @test all(isone, res1[:x].value)
        @test all(isone, res2[:x].value)
    end
    @testset "beta binomial" begin
        prior = Beta(2,2)
        obs = [0,1,0,1,1,1,1,1,1,1]
        exact = Beta(prior.α + sum(obs), prior.β + length(obs) - sum(obs))
        meanp = exact.α / (exact.α + exact.β)

        @model testbb(obs) = begin
            p ~ Beta(2,2)
            x ~ Bernoulli(p)
            for i = 1:length(obs)
                obs[i] ~ Bernoulli(p)
            end
            p, x
        end

        smc = SMC()
        pg = PG(10)
        gibbs = Gibbs(HMC(0.2, 3, :p), PG(10, :x))

        chn_s = sample(testbb(obs), smc, 1000)
        chn_p = sample(testbb(obs), pg, 2000)
        chn_g = sample(testbb(obs), gibbs, 1500)

        check_numerical(chn_s, [:p], [meanp], atol=0.05)
        check_numerical(chn_p, [:x], [meanp], atol=0.1)
        check_numerical(chn_g, [:x], [meanp], atol=0.1)
    end
    @testset "forbid global" begin
        xs = [1.5 2.0]
        # xx = 1

        @model fggibbstest(xs) = begin
            s ~ InverseGamma(2,3)
            m ~ Normal(0,sqrt(s))
            # xx ~ Normal(m, sqrt(s)) # this is illegal

            for i = 1:length(xs)
                xs[i] ~ Normal(m, sqrt(s))
                # for xx in xs
                # xx ~ Normal(m, sqrt(s))
            end
            s, m
        end

        gibbs = Gibbs(PG(10, :s), HMC(0.4, 8, :m))
        chain = sample(fggibbstest(xs), gibbs, 2);
    end
    @testset "model macro" begin
        model_info = Dict(:main_body_names => Dict(:vi => :vi, :sampler => :sampler))
        # unit test model macro
        expr = Turing.generate_observe(:x, :y, model_info)
        @test expr.head == :block
        @test :(vi.logp += Turing.observe(sampler, y, x, vi)) in expr.args

        @model testmodel_comp(x, y) = begin
            s ~ InverseGamma(2,3)
            m ~ Normal(0,sqrt(s))

            x ~ Normal(m, sqrt(s))
            y ~ Normal(m, sqrt(s))

            return x, y
        end
        testmodel_comp(1.0, 1.2)

        # check if drawing from the prior works
        @model testmodel0(x) = begin
            x ~ Normal()
            return x
        end
        f0_mm = testmodel0()
        @test mean(f0_mm() for _ in 1:1000) ≈ 0. atol=0.1

        # Test #544
        @model testmodel0(x = Vector{Float64}(undef, 2)) = begin
            x[1] ~ Normal()
            x[2] ~ Normal()
            return x
        end
        f0_mm = testmodel0()
        @test all(x -> isapprox(x, 0; atol = 0.1), mean(f0_mm() for _ in 1:1000))

        @model testmodel01(x) = begin
            x ~ Bernoulli(0.5)
            return x
        end
        f01_mm = testmodel01()
        @test mean(f01_mm() for _ in 1:1000) ≈ 0.5 atol=0.1

        # test if we get the correct return values
        @model testmodel1(x1, x2) = begin
            s ~ InverseGamma(2,3)
            m ~ Normal(0,sqrt(s))

            x1 ~ Normal(m, sqrt(s))
            x2 ~ Normal(m, sqrt(s))

            return x1, x2
        end
        f1_mm = testmodel1(1., 10.)
        @test f1_mm() == (1, 10)
        f1_mm = testmodel1(x1=1., x2=10.)
        @test f1_mm() == (1, 10)

        @info "Testing the compiler's ability to catch bad models..."

        # Test for assertions in observe statements.
        @model brokentestmodel_observe1(x1, x2) = begin
            s ~ InverseGamma(2,3)
            m ~ Normal(0,sqrt(s))

            x1 ~ Normal(m, sqrt(s))
            x2 ~ x1 + 2

            return x1, x2
        end

        btest = brokentestmodel_observe1(1., 2.)
        @test_throws ArgumentError btest()

        @model brokentestmodel_observe2(x) = begin
            s ~ InverseGamma(2,3)
            m ~ Normal(0,sqrt(s))

            x = Vector{Float64}(undef, 2)
            x ~ [Normal(m, sqrt(s)), 2.0]

            return x
        end

        btest = brokentestmodel_observe2([1., 2.])
        @test_throws ArgumentError btest()

        # Test for assertions in assume statements.
        @model brokentestmodel_assume1() = begin
            s ~ InverseGamma(2,3)
            m ~ Normal(0,sqrt(s))

            x1 ~ Normal(m, sqrt(s))
            x2 ~ x1 + 2

            return x1, x2
        end

        btest = brokentestmodel_assume1()
        @test_throws ArgumentError btest()

        @model brokentestmodel_assume2() = begin
            s ~ InverseGamma(2,3)
            m ~ Normal(0,sqrt(s))

            x = Vector{Float64}(undef, 2)
            x ~ [Normal(m, sqrt(s)), 2.0]

            return x
        end

        btest = brokentestmodel_assume2()
        @test_throws ArgumentError btest()
    end
    @testset "new grammar" begin
        x = Float64[1 2]

        @model gauss(x) = begin
            priors = TArray{Float64}(2)
            priors[1] ~ InverseGamma(2,3)         # s
            priors[2] ~ Normal(0, sqrt(priors[1])) # m
            for i in 1:length(x)
                x[i] ~ Normal(priors[2], sqrt(priors[1]))
            end
            priors
        end

        chain = sample(gauss(x), PG(10), 10)
        chain = sample(gauss(x), SMC(), 10)

        @model gauss2(x, ::Type{TV}=Vector{Float64}) where {TV} = begin
            priors = TV(undef, 2)
            priors[1] ~ InverseGamma(2,3)         # s
            priors[2] ~ Normal(0, sqrt(priors[1])) # m
            for i in 1:length(x)
                x[i] ~ Normal(priors[2], sqrt(priors[1]))
            end
            priors
        end

        chain = sample(gauss2(x), PG(10), 10)
        chain = sample(gauss2(x=x, TV=Vector{Float64}), PG(10), 10)
        chain = sample(gauss2(x), SMC(), 10)
        chain = sample(gauss2(x=x, TV=Vector{Float64}), SMC(), 10)
    end
    @testset "new interface" begin
        obs = [0, 1, 0, 1, 1, 1, 1, 1, 1, 1]

        @model newinterface(obs) = begin
          p ~ Beta(2,2)
          for i = 1:length(obs)
            obs[i] ~ Bernoulli(p)
          end
          p
        end

        chain = sample(
            newinterface(obs),
            HMC{Turing.ForwardDiffAD{2}}(0.75, 3, :p, :x),
            100)
    end
    @testset "no return" begin
        @model noreturn(x) = begin
            s ~ InverseGamma(2,3)
            m ~ Normal(0, sqrt(s))
            for i in 1:length(x)
                x[i] ~ Normal(m, sqrt(s))
            end
        end

        chain = sample(noreturn([1.5 2.0]), HMC(0.15, 6), 1000)
        check_numerical(chain, [:s, :m], [49/24, 7/6])
    end
    @testset "observe" begin
        @model test() = begin
          z ~ Normal(0,1)
          x ~ Bernoulli(1)
          1 ~ Bernoulli(x / 2)
          0 ~ Bernoulli(x / 2)
          x
        end

        is  = IS()
        smc = SMC()
        pg  = PG(10)

        res_is = sample(test(), is, 10000)
        res_smc = sample(test(), smc, 1000)
        res_pg = sample(test(), pg, 100)

        @test all(isone, res_is[:x].value)
        @test res_is.logevidence ≈ 2 * log(0.5)

        @test all(isone, res_smc[:x].value)
        @test res_smc.logevidence ≈ 2 * log(0.5)

        @test all(isone, res_pg[:x].value)
    end
    @testset "sample" begin
        alg = Gibbs(HMC(0.2, 3, :m), PG(10, :s))
        chn = sample(gdemo_default, alg, 1000);
    end
    @testset "vectorization" begin
        @model vdemo1(x) = begin
            s ~ InverseGamma(2,3)
            m ~ Normal(0, sqrt(s))
            x ~ [Normal(m, sqrt(s))]
            return s, m
        end

        alg = HMC(0.01, 5)
        x = randn(100)
        res = sample(vdemo1(x), alg, 250)

        D = 2
        @model vdemo2(x) = begin
            μ ~ MvNormal(zeros(D), ones(D))
            x ~ [MvNormal(μ, ones(D))]
        end

        alg = HMC(0.01, 5)
        res = sample(vdemo2(randn(D,100)), alg, 250)

        # Vector assumptions
        N = 10
        setchunksize(N)
        alg = HMC(0.2, 4)

        @model vdemo3() = begin
            x = Vector{Real}(undef, N)
            for i = 1:N
                x[i] ~ Normal(0, sqrt(4))
            end
        end

        t_loop = @elapsed res = sample(vdemo3(), alg, 1000)

        # Test for vectorize UnivariateDistribution
        @model vdemo4() = begin
          x = Vector{Real}(undef, N)
          x ~ [Normal(0, 2)]
        end

        t_vec = @elapsed res = sample(vdemo4(), alg, 1000)

        @model vdemo5() = begin
            x ~ MvNormal(zeros(N), 2 * ones(N))
        end

        t_mv = @elapsed res = sample(vdemo5(), alg, 1000)

        println("Time for")
        println("  Loop : $t_loop")
        println("  Vec  : $t_vec")
        println("  Mv   : $t_mv")

        # Transformed test
        @model vdemo6() = begin
            x = Vector{Real}(undef, N)
            x ~ [InverseGamma(2, 3)]
        end

        sample(vdemo6(), alg, 1000)
    end
    @testset "Type parameters" begin
        N = 10
        setchunksize(N)
        alg = HMC(0.01, 5)
        x = randn(1000)
        @model vdemo1(::Type{T}=Float64) where {T} = begin
            x = Vector{T}(undef, N)
            for i = 1:N
                x[i] ~ Normal(0, sqrt(4))
            end
        end

        t_loop = @elapsed res = sample(vdemo1(), alg, 250)
        t_loop = @elapsed res = sample(vdemo1(Float64), alg, 250)
        t_loop = @elapsed res = sample(vdemo1(T=Float64), alg, 250)

        @model vdemo2(::Type{T}=Float64) where {T <: Real} = begin
            x = Vector{T}(undef, N)
            x ~ [Normal(0, 2)]
        end

        t_vec = @elapsed res = sample(vdemo2(), alg, 250)
        t_vec = @elapsed res = sample(vdemo2(Float64), alg, 250)
        t_vec = @elapsed res = sample(vdemo2(T=Float64), alg, 250)

        @model vdemo3(::Type{TV}=Vector{Float64}) where {TV <: AbstractVector} = begin
            x = TV(undef, N)
            x ~ [InverseGamma(2, 3)]
        end

        sample(vdemo3(), alg, 250)
        sample(vdemo3(Vector{Float64}), alg, 250)
        sample(vdemo3(TV=Vector{Float64}), alg, 250)
    end
    @testset "tilde" begin
        model_info = Dict(
            :name => "model",
            :main_body_names => Dict(:model => :model,
                                    :vi => :vi,
                                    :sampler => :sampler),
            :arg_syms => [],
            :tent_pvars_list => [])

        ex = :(y ~ Normal(1,1))
        model_info[:main_body] = ex
        translate_tilde!(model_info)
        res = model_info[:main_body]
        Base.@assert res.head == :block

        ex = quote
            x = 1
            y = rand()
            y ~ Normal(0,1)
        end

        model_info[:main_body] = ex
        translate_tilde!(model_info)
        res = model_info[:main_body]
        Base.@assert res.head == :block
    end
end


text:
# TODO: Implement the log-likelihood, which marginalizes over lam, gam.
# This will be used for parallel tempering.
"""
Basing this off of conditional DIC:
p(y_in, m_in | theta, beta) \\propto
    p(y_obs | y_mis, theta) * p(m | y_obs, y_mis, beta)
We don't have this exactly when we marginalize out the indicators,
but we can have

p(y_complete, m | theta)

where y_complete is the (y_obs, y_mis), for a given iteration.
"""
function compute_marg_loglike(s::State, c::Constants, d::Data, temper::Float64;
                              y=nothing, compute_prob_miss=true)
  # Boolean: Whether y was supplied
  y_supplied = (y != nothing)

  if !y_supplied
    y = s.y_imputed
  end

  # log-likelihood
  ll = 0.0

  # Dimensions
  J = d.J
  K = c.K
  L = c.L

  # Precompute / reshapes
  mus0 = reshape(mus(false, s, c, d), 1, 1, L[0])
  mus1 = reshape(mus(true, s, c, d), 1, 1, L[1])
  sig = sqrt.(s.sig2)

  for i in 1:d.I
    # Dimensions
    Ni = if y_supplied
      size(y[i], 1)
    else
      d.N[i]
    end
    Z = reshape(s.Z, 1, J, K)
    eta0i = s.eta[false][i:i, :, :]
    eta1i = s.eta[true][i:i, :, :]

    # Ni x J x 1
    yi = reshape(y[i], Ni, J, 1)

    # Ni x J x 1
    kernel0 = MCMC.lpdf_gmm(yi, mus0, sig[i], eta0i, temper, dims=3)
    kernel1 = MCMC.lpdf_gmm(yi, mus1, sig[i], eta1i, temper, dims=3)

    # (1 x J x K) .* (Ni x 1 x K) -> Ni x K
    Z_mix = let
      tmp = sum(Z .* kernel1 + (1 .- Z) .* kernel0, dims=2)
      dropdims(tmp, dims=2)
    end

    # (Ni x K) .+ (1 x K) --> (Ni x K)
    f = Z_mix .+ log.(s.W[i:i, :])

    # Ni-dim
    lli = MCMC.logsumexp(f, dims=2)

    # Increment loglikelihood
    ll += sum(lli)

    if compute_prob_miss
      # Add probability of missing for only imputed values only, because for
      # observed values, the probability evaluates to a constant across models.
      mi = Bool.(d.m[i])
      pm = [Model.prob_miss(y_inj, c.beta[:, i]) for y_inj in s.y_imputed[i][mi]]

      # Increment loglikelihood with missing prob stuff
      ll += sum(log.(pm))
    end
  end  # i loop

  return ll
end


text:
using FastRationals
using FastRationals: IsReduced, MayReduce

using Test
using LinearAlgebra

# adapted from Base/test/rational.jl

const FR = FastRational

#=
@testset "Rationals" begin
    @test FR(1//1) == 1
    @test FR(2//2) == 1
    @test FR(1//1) == 1//1
    @test FR(2//2) == 1//1
    @test FR(2//4) == FR(8//16)
    @test FR(1//2) + FR(1//2) == 1
    @test FR((-1)//3) == FR(-(1//3)) == -(1//3)
    @test FR(1//2) + FR(3//4) == FR(5//4) == 5//4
    @test FR(1//3) * FR(3//4) == FR(1//4) == 1//4
    @test FR(1//2) / FR(3//4) == FR(2//3)

    @test_throws OverflowError -FR(0x01//0x0f)
    @test_throws OverflowError -FR(typemin(Int)//1)
    @test_throws OverflowError FR(typemax(Int)//3) + 1
    @test_throws OverflowError FR(typemax(Int)//3) * 2
    @test FR(typemax(Int)//1) * FR(1//typemax(Int)) == 1
    @test FR(typemax(Int)//1) / FR(typemax(Int)//1) == 1
    @test FR(1//typemax(Int)) / FR(1//typemax(Int)) == 1
    @test_throws OverflowError FR(1//2)^63

    for a = -5:5, b = -5:5
        if a == b == 0; continue; end
        if ispow2(b)
            @test FR(a//b) == FR(a/b)
            @test convert(FastRational,a/b) == FR(a//b)
        end
        # @test rationalize(a/b) == a//b
        @test FR(a//b) == FR(a//b)
        if b == 0
            @test_throws DivideError round(Integer,FR(a//b)) == round(Integer,a/b)
        else
            @test round(Integer,FR(a//b)) == round(Integer,a/b)
        end
        for c = -5:5
            @test (FR(a//b) == c) == (a/b == c)
            @test (FR(a//b) != c) == (a/b != c)
            @test (FR(a//b) <= c) == (a/b <= c)
            @test (FR(a//b) <  c) == (a/b <  c)
            @test (FR(a//b) >= c) == (a/b >= c)
            @test (FR(a//b) >  c) == (a/b >  c)
            for d = -5:5
                if b == d == 0; continue; end
                if c == d == 0; continue; end
                @test (FR(a//b) == FR(c//d)) == (a/b == c/d)
                @test (FR(a//b) != FR(c//d)) == (a/b != c/d)
                @test (FR(a//b) <= FR(c//d)) == (a/b <= c/d)
                @test (FR(a//b) <  FR(c//d)) == (a/b <  c/d)
                @test (FR(a//b) >= FR(c//d)) == (a/b >= c/d)
                @test (FR(a//b) >  FR(c//d)) == (a/b >  c/d)
            end
        end
    end

    @test 0.5 == FR(1//2)
    @test 0.1 != FR(1//10)
    @test 0.1 == FR(3602879701896397//36028797018963968)

    @test 1/3 < FR(1//3)
    @test !(FR(1//3) < 1/3)
    @test -1/3 < FR(1//3)
    @test -1/3 > FR(-1//3)
    @test 1/3 > FR(-1//3)
    @test 1/5 > FR(1//5)
    @test !(1//3 < NaN)
    @test !(1//3 == NaN)
    @test !(1//3 > NaN)

    # PR 29561
    @test abs(one(FastRational{UInt,IsReduced})) === one(FastRational{UInt,IsReduced})
    @test abs(one(FastRational{Int,IsReduced})) === one(FastRational{Int,IsReduced})
    @test abs(-one(FastRational{Int,MayReduce})) === one(FastRational{Int,MayReduce})
end
=#

@testset "tests" begin include("tests.jl") end


text:
include("utils.jl")

function task_sample()
    parent = current_task()
    task = @task begin
        ptr1 = sample_stackpointer()
        #=
        bt_data = ccall(
            :jl_backtrace_from_here,
            Ref{Base.SimpleVector},
            (Cint, Cint),
            true,
            0,
        )
        =#
        yieldto(parent, ptr1)
        # Suspend the task using `yieldto` to avoid the stack buffer to be GC'ed
    end
    return task, yieldto(task)
end

"""
    guess_task_layout(task, ptr1) -> offset_bytes

Return the guessed offset of the hidden field `stkbuf` of `jl_task_t` in bytes.

`ptr1` is a frame address of `task` "reasonably early" in the call stack.
"""
function guess_task_layout(task, ptr1)
    GC.@preserve task begin
        offset_base = (sizeof(Task) ÷ sizeof(UInt)) * sizeof(UInt)
        ptr0 = pointer_from_objref(task) + offset_base
        for offset in 0:64
            pstkbuf = Ptr{Ptr{Cvoid}}(ptr0 + sizeof(UInt) * offset)
            pbufsz = Ptr{Csize_t}(pstkbuf + sizeof(UInt))
            stkbuf = unsafe_load(pstkbuf)
            bufsz = unsafe_load(pbufsz)
            stkbase = stkbuf + bufsz
            if stkbase - ptr1 < 2000
                return offset_base + offset * sizeof(UInt)
            end
        end
    end
    return nothing
end

if !@isdefined(only)
    function only(xs)
        xs = collect(xs)
        @assert length(xs) == 1
        return xs[1]
    end
end

function guess_task_layout()
    samples = [task_sample() for _ in 1:100]
    guesses = [guess_task_layout(args...) for args in samples]
    return only(unique(guesses))
end

if realpath(PROGRAM_FILE) == realpath(@__FILE__)
    sample_stackpointer()
    println(guess_task_layout())
end


text:
function instantiate_network(network::Type{N}, sys::PowerSystems.PowerSystem; args...) where N <: AbstractDCPowerModel

    d_netinjection_p =  JumpAffineExpressionArray(undef, length(sys.buses), sys.time_periods)

    ts_active = active_timeseries_netinjection(sys)

    netinjection = (var_active = d_netinjection_p, var_reactive = nothing, timeseries_active = ts_active, timeseries_reactive = nothing)

    return netinjection

end

function instantiate_network(network::Type{N}, sys::PowerSystems.PowerSystem; args...) where N <: AbstractACPowerModel

    d_netinjection_p =  JumpAffineExpressionArray(undef, length(sys.buses), sys.time_periods)

    ts_active = active_timeseries_netinjection(sys)

    d_netinjection_q =  JumpAffineExpressionArray(undef, length(sys.buses), sys.time_periods)

    ts_reactive = reactive_timeseries_netinjection(sys)

    netinjection = (var_active = d_netinjection_p, var_reactive = d_netinjection_q, timeseries_active = ts_active, timeseries_reactive = ts_reactive)

    return netinjection

end

function instantiate_network(network::Type{N}, sys::PowerSystems.PowerSystem; args...) where N <: CopperPlatePowerModel

    return instantiate_network(AbstractDCPowerModel, sys)

end

function instantiate_network(network::Type{N}, sys::PowerSystems.PowerSystem; args...) where N <: AbstractFlowForm

    return instantiate_network(AbstractDCPowerModel, sys)

end


text:
import Pkg
import TOML

function is_manifest_filename(file::AbstractString)
    is_manifest = any(startswith.(Ref(file), first.(splitext.(Base.manifest_names))))
    is_toml = endswith(lowercase(file), ".toml")
    return is_manifest && is_toml
end

function _shorten_manifest_version_master(; manifest_filename::AbstractString)
    julia_version_key = "julia_version"
    manifest_dict = TOML.parsefile(manifest_filename)
    old_julia_version = get(manifest_dict, julia_version_key, nothing)
    if old_julia_version !== nothing
        m = match(
            r"^([0-9]*?)\.([0-9]*?)\.([0-9]*?)\-DEV\.([0-9]*?)$",
            old_julia_version,
        )
        if m !== nothing
            major = m[1]
            minor = m[2]
            patch = m[3]
            new_julia_version = "$(major).$(minor).$(patch)-DEV"
            manifest_dict[julia_version_key] = new_julia_version
            @info "Fixing manifest file at $(manifest_filename)" old_julia_version new_julia_version
            Pkg.Types.write_manifest(manifest_dict, manifest_filename)
        end
    end
    return nothing
end

function shorten_manifest_version_master(dir::AbstractString)
    for (root, dirs, files) in walkdir(abspath(dir))
        for file in files
            if is_manifest_filename(file)
                manifest_filename = joinpath(root, file)
                @info "Checking manifest file at $(manifest_filename)"
                _shorten_manifest_version_master(; manifest_filename)
            end
        end
    end
    return nothing
end


text:
module Gettext

using PyCall

@pyimport gettext as gt

bindtextdomain(domain, localedir=nothing) = gt.bindtextdomain(domain, localedir)
textdomain(domain=nothing) = gt.textdomain(domain)

gettext(message) = gt.gettext(message)
dgettext(domain, message) = gt.dgettext(domain, message)
ngettext(singular, plural, n) = gt.ngettext(singular, plural, n)
dngettext(domain, singular, plural, n) = gt.dngettext(domain, singular, plural, n)

macro __str(s)
    gettext(s)
end

macro N__str(s)
    s
end

export bindtextdomain, textdomain, gettext, lgettext, dgettext, dlgettext, ngettext, lngettext, dngettext, ldngettext, @__str, @N__str

end # module


text:
export Q_update!

"""
update the Q matrix in SparseMatrixCSC or CDS format. Q=sum_{i=1}^p AtA[i]rho[i]
if any of the rho[i] change
"""
function Q_update!(
                  Q            ::Union{Array{TF,2},SparseMatrixCSC{TF,Int64},joAbstractLinearOperator{TF, TF}},
                  AtA          ::Union{Vector{Array{TF, 2}}, Vector{SparseMatrixCSC{TF, TI}}, Vector{joAbstractLinearOperator{TF, TF}}, Vector{Union{Array{TF, 2}, SparseMatrixCSC{TF, TI}, joAbstractLinearOperator{TF, TF}}}},
                  set_Prop,
                  rho          ::Vector{TF},
                  ind_updated  ::Vector{TI},
                  log_PARSDMM,
                  i            ::Integer,
                  Q_offsets=[] ::Vector{TI}
                  ) where {TF<:Real,TI<:Integer}

  if isempty(ind_updated) == false
    if typeof(Q) <: joAbstractLinearOperator #some of the AtA are JOLI operators
      #JOLI operators cannot be defined recursively for many iterations for many operators
      #because this causes memory issues and slow performance -> re-assemble Q from scratch
      # for ii=1:length(ind_updated)
      #   if typeof(AtA[ind_updated[ii]]) <: joAbstractLinearOperator
      #     Q = Q + ( AtA[ind_updated[ii]] )*( rho[ind_updated[ii]]-log_PARSDMM.rho[i,ind_updated[ii]] );
      #   else
      #     Q = Q + ( joMatrix(AtA[ind_updated[ii]]) )*( rho[ind_updated[ii]]-log_PARSDMM.rho[i,ind_updated[ii]] );
      #   end
      # end
      if typeof(AtA[1]) <: joAbstractLinearOperator
        Q = rho[1]*AtA[1]
      else
        Q = rho[1]*joMatrix(AtA[1])
      end
      for i=2:length(AtA)
        if typeof(AtA[i]) <: joAbstractLinearOperator
          Q = Q + rho[i]*AtA[i]
        else
          Q = Q + rho[i]*joMatrix(AtA[i])
        end
      end
    elseif typeof(Q) == SparseMatrixCSC{TF,Int64} #all AtA are sparse arrays
      for ii=1:length(ind_updated)
          Q = Q + ( AtA[ind_updated[ii]] )*( rho[ind_updated[ii]]-log_PARSDMM.rho[i,ind_updated[ii]] );
      end
    else #all AtA are in CDS format
      for ii=1:length(ind_updated)
        CDS_scaled_add!(Q,AtA[ind_updated[ii]],Q_offsets,set_Prop.AtA_offsets[ind_updated[ii]],rho[ind_updated[ii]] .- log_PARSDMM.rho[i,ind_updated[ii]] )
      end
    end
  end

return Q
end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
using SRKGenerator
#Standard Driver Script

if isempty(ARGS)

  ## Called from test. Use a default setup
  imin = -12
  jmax = 1

  ## Parallelize correctly
  host = gethostname()
  if host == "ChrisRackBed"
    #numCores = 4
    #addprocs(numCores-length(workers()))
    ptx_str = joinpath(Pkg.dir("SRKGenerator"),"deps","integrationWin.ptx")
    cudaCores=[1664,640]
  end

  if host == "crackauc2"
    #numCores = 16
    #addprocs(numCores-length(workers()))
    ptx_str = joinpath(Pkg.dir("SRKGenerator"),"deps","integrationLinux.ptx")
    cudaCores = [2816]
  end

  if host == "ChrisRackTV"
    #numCores = 4
    #addprocs(numCores-length(workers()))
    ptx_str = joinpath(Pkg.dir("SRKGenerator"),"deps","integrationWin.ptx")
    cudaCores=[2560]
  end

else

  ## Called from commandline. Assume Linux.
  ptx_str = joinpath(Pkg.dir("SRKGenerator"),"deps","integrationLinux.ptx")
  ## Arg Form:
  cudaCores = [parse(Int,ARGS[1])]
  imin = parse(Int,ARGS[2])
  jmax = parse(Int,ARGS[3]); jmin = -jmax

end

resString = srk_optimize(:LN_AUGLAG_EQ,1/100,0,imin,1,-jmax,jmax,
                         rand_minmax = rand(1:3),
                         NLoptRandSeed=rand(1:Int(1e8)),gpuEnabled = true,
                         cudaCores=cudaCores,ptx_str=ptx_str)


text:
#
# backtrace_list() = [ccall(:jl_lookup_code_address, Any, (Ptr{Void}, Int32), b, 0) for b in backtrace()]
#
# const IGNORED_METHODS = Set{Symbol}([:LOG, :publish, :schedule_do])
#
# function backtracestring()
#
#   # get the backtrace and find the position of the "backtracestring()" call
#   btlist = backtrace_list()
#   btpos = find(x->x[1] == :backtracestring, btlist)[1]
#   # println(btlist)
#
#   # find the first non-LOG, non-julia method and return that string
#   for i in btpos+1:length(btlist)
#
#     # skip julia internals and LOG calls
#     methodsym = btlist[i][1]
#     if methodsym in IGNORED_METHODS
#       continue
#     end
#
#     methodname, pathname, filenum, tmp1, tmp2 = map(string, btlist[i])
#     if methodname[1:3] == "jl_"
#       continue
#     end
#
#     # get the filename from the full path, then return the string
#     filename = split(pathname, "/")[end]
#     if filename == "FastAnonymous.jl"
#       continue
#     end
#
#     return "[$filename:$filenum]"
#   end
#   "[]"
# end
#
#
#
# # --------------------------------------------------------
#
#
# immutable LogSeverity
#   val::Int
# end
#
# const DebugSeverity = LogSeverity(0)
# const InfoSeverity = LogSeverity(1)
# const ErrorSeverity = LogSeverity(2)
#
#
# Base.string(sev::LogSeverity) = (sev == DebugSeverity ? "Debug" : (sev == InfoSeverity ? "Info" : "Error"))
# Base.print(io::IO, sev::LogSeverity) = print(io, string(sev))
# Base.show(io::IO, sev::LogSeverity) = print(io, string(sev))
# Base.isless(sev1::LogSeverity, sev2::LogSeverity) = sev1.val < sev2.val
#
# # --------------------------------------------------------
#
# type SevObj
#   sev::LogSeverity
#   io::IO
# end
#
# const LOG_SEVERITY = SevObj(InfoSeverity, STDOUT)
#
# log_severity() = LOG_SEVERITY.sev
# log_severity!(sev::LogSeverity) = (LOG_SEVERITY.sev = sev; nothing)
#
# log_io() = LOG_SEVERITY.io
# log_io!(io::IO) = (LOG_SEVERITY.io = io; nothing)
#
#
# # --------------------------------------------------------
#
# LOG(args...) = LOG(InfoSeverity, args...)
# DEBUG(args...) = LOG(DebugSeverity, args...)
# ERROR(args...) = LOG(ErrorSeverity, args...)
#
# function LOG(sev::LogSeverity, args...)
#   if sev >= log_severity()
#     io = log_io()
#     print(io, "$(now()) [$sev]: ")
#     for arg in args
#       print(io, arg, " ")
#     end
#     println(io, backtracestring())
#       # join(vcat(map(string,args)..., backtracestring()), " "))
#   end
# end
#
# # note: the macro version give "x: xval" for "@LOG x"
#
#
# # TODO: maybe change the macros a little:
#   # can we get/hardcode the file/linenumber within the macro body?
#   # can we wrap the call with an "if log_severity() >= DebugSeverity; DEBUG(...); end" for speed?
#
#
# # default to InfoSeverity
# macro LOG(symbols...)
#   expr = :(LOG(InfoSeverity))
#   for s in symbols
#     push!(expr.args, "$s:")
#     push!(expr.args, esc(s))
#   end
#   expr
# end
#
# macro ERROR(symbols...)
#   expr = :(LOG(ErrorSeverity))
#   for s in symbols
#     push!(expr.args, "$s:")
#     push!(expr.args, esc(s))
#   end
#   expr
# end
#
# macro DEBUG(symbols...)
#   expr = :(LOG(DebugSeverity))
#   for s in symbols
#     push!(expr.args, "$s:")
#     push!(expr.args, esc(s))
#   end
#   expr
# end