text:
function Base.modf(x::Double64)
    frc_hi, int_hi = modf(Hi(x))
    frc_lo, int_lo = modf(Lo(x))
    frc_part = Double64((frc_hi, frc_lo))
    int_part = Double64((int_hi, int_lo))
    return (frc_part, int_part)
end

fmod(parts::Tuple{Double64, Double64}) = fmod(parts...)
function fmod(frc_part::Double64, int_part::Double64)
    return int_part + frc_part
end

const log10of2 = log10(2)
const log2of10 = log2(10)

@inline digits2bits(x::Int) = floor(Int, log2of10*x)
@inline bits2digits(x::Int) = ceil(Int, log10of2*x)

# fractionalbits(x) is the smallest number of bits such that
# frac(x) = modf(x)[1]
# isinteger(frac(x) * 2.0^fractionalbits(frac(x)))
# fix for 0.0 and powers of 2
function fractionalbits(x::Float64)
    fr = modf(x)[1]
    nbits = 53 - trailing_zeros((reinterpret(UInt64, fr) & 0x000f_ffff_ffff_ffff))
    return nbits
end

# ~2^(-105) q.v. Error bounds from extra-precise iterative refinement
# ACM Transactions on Mathematical Software Vol. 32, No. 2 pg 337
# .. IEEE-754 double precision as working precision (εw = 2−53) and double-double as residual precision (εr ≈ 2−105)."

Base.eps(::Type{Double64})   = 2.465190328815662e-32
Base.eps(x::Double64) = iszero(Lo(x)) ? Hi(x) * 2.465190328815662e-32 : eps(Lo(x))

#=
   faster 2.0^p
   for p in trunc(Int,log2(nextfloat(0.0))) .. trunc(Int, log2(prevfloat(floatmax(Float64), 354))
   for smaller p --> 0.0 for larger p --> Inf
=#
ldexp2pow(p::Float64) = ldexp(1.0, p)
#=
    ufp(x) is "unit in the first place"
    ufp(0) = 0
    ufp(x!=0) = 2.0^floor(log2(abs(x)))

    unsafe_ufp(x) errors on x in (NaN, Inf)
=#
@inline unsafe_ufp(x::Float64) = !iszero(x) ? ldexp(1.0, exponent(x)) : x
unsafe_ufp(x::Double64) = unsafe_ufp(Hi(x))

ufp(x::Float64) = isfinite(x) ? unsafe_ufp(x) : x
ufp(x::Double64) = ufp(Hi(x))

#=
Formal Verification of a Floating-Point Expansion Renormalization Algorithm
by Sylvie Boldo, Mioara Joldes, Jean-Michel Muller, Valentina Popescu
8th International Conference on Interactive Theorem Proving (ITP’2017),
Sep 2017, Brasilia, Brazil

Definition 21 A normal binary precision-p floating-point (FP) number has the
form x = M_x · 2^(e_x - p + 1) with 2^(p−1) ≤ |Mx| ≤ 2^(p)−1.
The integer e_x is the exponent of x, and M_x · 2^(-p + 1) is the significand of x.

We denote
ulp(x) = 2^(e_x−p+1) # unit in the last place
uls(x) = ulp(x) * 2^(z_x) # unit in last significant place
  where z_x is the number of trailing zeros at the end of M_x

Algorithms for triple-word arithmetic
by Nicolas Fabiano and Jean-Michel Muller and Joris Picot

[We denote]
ufp(x) = 2^(floor(log2(abs(x)))) # unit in the first place
ulp(x) = ufp(x) * 2^(-p+1)       # for Float64, p=53
ulp(x) = ufp(x) * 2^(-52)
uls(x) = the largest 2^k such that isinteger(x/2^k)

ufp(x) is the weight of of its most significant bit
ulp(x) is the weight of of its least significant bit
uls(x) is the weight of its rightmost nonzero bit

u = 2^(-p) == ulp(1.0)/2 # roundoff error unit
=#

ulp(x::Float64) = ldexp(ufp(x), -52) # or eps(x)

function uls(x::Float64)
    k =  54 - trailing_zeros(reinterpret(UInt64, x))
    return ldexp(ufp(x), -k)
end

# -------------------------------------------------------

# significant_significand_bits(x::Float64)
ssbits(x::Float64) =
    exponent(x) - (x >= 1.0) +
    (54 - trailing_zeros(reinterpret(UInt64, x)))

# --------------------------------------------------------

function nearestint(x::T) where {T<:Union{Float64, Double64}}
    s, absx = signbit(x), abs(x)
    absx = absx + 0.5
    absx = trunc(absx)
    return s ? -absx : absx
end

function accurate_c1c2(value::Real, ::Type{T}, p=significant_bits(T)) where {T<:Real}
    r1 = one(typeof(value)) / value
    r = T(r1)
    ir = one(typeof(r)) / r
    irulp = T(4 * eps(ir))
    c1a = nearestint(inv(r * irulp))
    c1 = T(c1a * irulp)
    p2 = T(2.0)^(-p+4)
    scaledulp = T(p2 * eps(c1))
    c2a = nearestint( (value - c1) / scaledulp )
    c2 = T(c2a * scaledulp)
    return c1, c2
end


text: println(sum(1:100)^2 - sum(i^2 for i in 1:100))

ext: jl

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
using Flux
using Flux: @epochs
using Flux: params
using Flux: throttle
using Plots
using Distributions
using DataFrames
using CSV
using StaticArrays

@inline phi(x) = max(x^3,0)
@inline f(x) = 1

#build a block
block_layer = SkipConnection(Chain(Dense(10,10, phi),Dense(10,10, phi)), +)

#construct network skeleton
M = Chain(Dense(2,10),
block_layer,
block_layer,
block_layer,
block_layer,
Dense(10,1))

#check number of params
params(M)

#inner loss function
function inner_loss(N, space)
    a = rand(Uniform(-1,1),2,N)
    loss = 0
    wt = 4*ones(N)/N
    l11 = [1, 0]
    l21 = [2, 0]
    l12 = [0, 1]
    l22 = [0, 2]

    @inbounds loss += ((0.5*((-1/12*M(a .+ space.*l21) .+ 2/3*M(a .+ space.*l11) .- 2/3*M(a .- space.*l11) .+ 1/12*M(a .- space.*l21))./space).^2) *wt)[1]
    @inbounds loss += ((0.5*((-1/12*M(a .+ space.*l22) .+ 2/3*M(a .+ space.*l12) .- 2/3*M(a .- space.*l12) .+ 1/12*M(a .- space.*l22))./space).^2) *wt)[1]
    @inbounds loss += (-1 .*M(a)*wt)[1]

    return loss
end

#boundary condition loss function
#add [0,1) as a boundary
function bc_loss(N, β)
    dist = Uniform(-1,1)
    a = hcat(vcat(rand(dist, 1,floor(Int32, N/5)), ones(1,floor(Int32, N/5))),vcat(rand(dist, 1,floor(Int32, N/5)), ones(1,floor(Int32, N/5)).*-1),
    vcat(ones(1,floor(Int32, N/5)), rand(dist, 1,floor(Int32, N/5))), vcat(ones(1,floor(Int32, N/5)).*-1, rand(dist, 1,floor(Int32, N/5))),
    vcat(zeros(1,floor(Int32, N/5)), rand(1,floor(Int32, N/5))))
    out_original = M(a)
    loss = (β/N)*sum(abs2, out_original)
    return loss
end

#add losses
function total_loss(N, space = 0.01, β = 750)
    return inner_loss(N, space) + bc_loss(floor(Int32, N/2), β)
end

#train model
function train_model(N = 100, train_time = 20)
    df = DataFrame(time = Int[], in_loss = Float64[],bound_loss = Float64[])
    c = 1
    n = floor(Int32, N/2)
    evalcb = function()
        @show(total_loss(N))
        in_loss = inner_loss(N, 0.01)
        bound_loss = bc_loss(n, 750)
        push!(df, [c, in_loss, bound_loss])
        c += 1
        if c == train_time
            CSV.write("$(c)_2d_cutsq.csv",df)
        end
    end
    opt = ADAM()
    dataset = [(N) for i=1:train_time]
    Flux.train!(total_loss, params(M), zip(dataset), opt, cb=throttle(evalcb, 10))
end

@epochs 1000 train_model(4000)

#plot the solution over a mesh in the domain
function plot_surface()
    data = zeros(201,201)
    grid_space = [i for i in -100:100]
    for i in grid_space
        for j in grid_space
            point = [i/100 j/100]
            @inbounds data[i+101,j+101] = M(point')[1]
        end
    end

    gr()
    heatmap(1:size(data,1),
        1:size(data,2), data,
        c = cgrad([:blue, :white,:red, :yellow]),
        xlabel = "X", ylabel = "Y",
        title = "2D Cut Square Fitted Solution")

    savefig("2d_cutsq_fitted.png")
end

plot_surface()


text:
##################################################################
# homogeneous mixture: arbitrary number of the same distribution #
##################################################################

"""
    HomogeneousMixture(distribution::Distribution, dims::Vector{Int})

Define a new distribution that is a mixture of some number of instances of single base distributions.

The first argument defines the base distribution of each component in the mixture.

The second argument must have length equal
to the number of arguments taken by the base distribution.  A value of 0
at a position in the vector an indicates that the corresponding argument to the
base distribution is a scalar, and integer values of i for i >= 1 indicate that
the corresponding argument is an i-dimensional array.

Example:

```julia
mixture_of_normals = HomogeneousMixture(normal, [0, 0])
```

The resulting distribution (e.g. `mixture_of_normals` above) can then be used like the built-in distribution values like `normal`.
The distribution takes `n+1` arguments where `n` is the number of arguments taken by the base distribution.
The first argument to the distribution is a vector of non-negative mixture weights, which must sum to 1.0.
The remaining arguments to the distribution correspond to the arguments of the base distribution, but have a different type:
If an argument to the base distribution is a scalar of type `T`, then the corresponding argument to the mixture distribution is a `Vector{T}`, where each element of this vector is the argument to the corresponding mixture component.
If an argument to the base distribution is an `Array{T,N}` for some `N`, then the corresponding argument to the mixture distribution is of the form `arr::Array{T,N+1}`, where each slice of the array of the form `arr[:,:,...,i]` is the argument for the `i`th mixture component.

Example:

```julia
mixture_of_normals = HomogeneousMixture(normal, [0, 0])
mixture_of_mvnormals = HomogeneousMixture(mvnormal, [1, 2])

@gen function foo()
    # mixture of two normal distributions
    # with means -1.0 and 1.0
    # and standard deviations 0.1 and 10.0
    # the first normal distribution has weight 0.4; the second has weight 0.6
    x ~ mixture_of_normals([0.4, 0.6], [-1.0, 1.0], [0.1, 10.0])

    # mixture of two multivariate normal distributions
    # with means: [0.0, 0.0] and [1.0, 1.0]
    # and covariance matrices: [1.0 0.0; 0.0 1.0] and [10.0 0.0; 0.0 10.0]
    # the first multivariate normal distribution has weight 0.4;
    # the second has weight 0.6
    means = [0.0 1.0; 0.0 1.0] # or, cat([0.0, 0.0], [1.0, 1.0], dims=2)
    covs = cat([1.0 0.0; 0.0 1.0], [10.0 0.0; 0.0 10.0], dims=3)
    y ~ mixture_of_mvnormals([0.4, 0.6], means, covs)
end
```
"""
struct HomogeneousMixture{T} <: Distribution{T}
    base_dist::Distribution{T}
    dims::Vector{Int}
end

(dist::HomogeneousMixture)(args...) = random(dist, args...)

Gen.has_output_grad(dist::HomogeneousMixture) = has_output_grad(dist.base_dist)
Gen.has_argument_grads(dist::HomogeneousMixture) = (true, has_argument_grads(dist.base_dist)...)
Gen.is_discrete(dist::HomogeneousMixture) = is_discrete(dist.base_dist)

function args_for_component(dist::HomogeneousMixture, k::Int, args)
    # returns a generator
    return (arg[fill(Colon(), dim)..., k]
            for (arg, dim) in zip(args, dist.dims))
end

function Gen.random(dist::HomogeneousMixture, weights, args...)
    k = categorical(weights)
    return random(dist.base_dist, args_for_component(dist, k, args)...)
end

function Gen.logpdf(dist::HomogeneousMixture, x, weights, args...)
    K = length(weights)
    log_densities = [Gen.logpdf(dist.base_dist, x, args_for_component(dist, k, args)...) for k in 1:K]
    log_densities = log_densities .+ log.(weights)
    return logsumexp(log_densities)
end

function Gen.logpdf_grad(dist::HomogeneousMixture, x, weights, args...)
    K = length(weights)
    log_densities = [Gen.logpdf(dist.base_dist, x, args_for_component(dist, k, args)...) for k in 1:K]
    log_weighted_densities = log_densities .+ log.(weights)
    relative_weighted_densities = exp.(log_weighted_densities .- logsumexp(log_weighted_densities))

    # log_grads[k] contains the gradients for the k'th component
    log_grads = [Gen.logpdf_grad(dist.base_dist, x, args_for_component(dist, k, args)...) for k in 1:K]

    # compute gradient with respect to x
    log_grads_x = [log_grad[1] for log_grad in log_grads]
    x_grad = if has_output_grad(dist.base_dist)
        sum(log_grads_x .* relative_weighted_densities)
    else
        nothing
    end

    # compute gradients with respect to the weights
    weights_grad = exp.(log_densities .- logsumexp(log_weighted_densities))

    # compute gradients with respect to each argument
    arg_grads = Any[]
    for (i, (has_grad, arg, dim)) in enumerate(zip(has_argument_grads(dist)[2:end], args, dist.dims))
        if has_grad
            if dim == 0
                grads = [log_grad[i+1] for log_grad in log_grads]
                grad_weights = relative_weighted_densities
            else
                grads = cat(
                    [log_grad[i+1] for log_grad in log_grads]...,
                    dims=dist.dims[i]+1)
                grad_weights = reshape(
                    relative_weighted_densities,
                    (1 for d in 1:dist.dims[i])..., K)
            end
            push!(arg_grads, grads .* grad_weights)
        else
            push!(arg_grads, nothing)
        end
    end

    return (x_grad, weights_grad, arg_grads...)
end

export HomogeneousMixture

##############################################################################
# heterogeneous mixture: fixed number of potentially different distributions #
##############################################################################

"""
    HeterogeneousMixture(distributions::Vector{Distribution{T}}) where {T}

Define a new distribution that is a mixture of a given list of base distributions.

The argument is the vector of base distributions, one for each mixture component.

Note that the base distributions must have the same output type.

Example:
```julia
uniform_beta_mixture = HeterogeneousMixture([uniform, beta])
```

The resulting mixture distribution takes `n+1` arguments, where `n` is the sum of the number of arguments taken by each distribution in the list.
The first argument to the mixture distribution is a vector of non-negative mixture weights, which must sum to 1.0.
The remaining arguments are the arguments to each mixture component distribution, in order in which the distributions are passed into the constructor.

Example:
```julia
@gen function foo()
    # mixure of a uniform distribution on the interval [`lower`, `upper`]
    # and a beta distribution with alpha parameter `a` and beta parameter `b`
    # the uniform as weight 0.4 and the beta has weight 0.6
    x ~ uniform_beta_mixture([0.4, 0.6], lower, upper, a, b)
end
```
"""
struct HeterogeneousMixture{T} <: Distribution{T}
    K::Int
    distributions::Vector{Distribution{T}}
    has_output_grad::Bool
    has_argument_grads::Tuple
    is_discrete::Bool
    num_args::Vector{Int}
    starting_args::Vector{Int}
end

(dist::HeterogeneousMixture)(args...) = random(dist, args...)

Gen.has_output_grad(dist::HeterogeneousMixture) = dist.has_output_grad
Gen.has_argument_grads(dist::HeterogeneousMixture) = dist.has_argument_grads
Gen.is_discrete(dist::HeterogeneousMixture) = dist.is_discrete

const MIXTURE_WRONG_NUM_COMPONENTS_ERR = "the length of the weights vector does not match the number of mixture components"

function HeterogeneousMixture(distributions::Vector{Distribution{T}}) where {T}
    _has_output_grad = true
    _has_argument_grads = Bool[true] # weights
    _is_discrete = true
    for dist in distributions
        _has_output_grad = _has_output_grad && has_output_grad(dist)
        for has_arg_grad in has_argument_grads(dist)
            push!(_has_argument_grads, has_arg_grad)
        end
        _is_discrete = _is_discrete && is_discrete(dist)
    end
    num_args = Int[]
    starting_args = Int[]
    for dist in distributions
        push!(starting_args, sum(num_args) + 1)
        push!(num_args, length(has_argument_grads(dist)))
    end
    K = length(distributions)
    return HeterogeneousMixture{T}(
        K, distributions,
        _has_output_grad,
        tuple(_has_argument_grads...),
        _is_discrete,
        num_args,
        starting_args)
end

function extract_args_for_component(dist::HeterogeneousMixture, component_args_flat, k::Int)
    start_arg = dist.starting_args[k]
    n = dist.num_args[k]
    return component_args_flat[start_arg:start_arg+n-1]
end

function Gen.random(dist::HeterogeneousMixture{T}, weights, component_args_flat...) where {T}
    (length(weights) != dist.K) && error(MIXTURE_WRONG_NUM_COMPONENTS_ERR)
    k = categorical(weights)
    value::T = random(
        dist.distributions[k],
        extract_args_for_component(dist, component_args_flat, k)...)
    return value
end

function Gen.logpdf(dist::HeterogeneousMixture, x, weights, component_args_flat...)
    (length(weights) != dist.K) && error(MIXTURE_WRONG_NUM_COMPONENTS_ERR)
    log_densities = [Gen.logpdf(
            dist.distributions[k], x,
            extract_args_for_component(dist, component_args_flat, k)...)
        for k in 1:dist.K]
    log_densities = log_densities .+ log.(weights)
    return logsumexp(log_densities)
end

function Gen.logpdf_grad(dist::HeterogeneousMixture, x, weights, component_args_flat...)
    (length(weights) != dist.K) && error(MIXTURE_WRONG_NUM_COMPONENTS_ERR)
    log_densities = [Gen.logpdf(
            dist.distributions[k], x,
            extract_args_for_component(dist, component_args_flat, k)...)
        for k in 1:dist.K]
    log_weighted_densities = log_densities .+ log.(weights)
    relative_weighted_densities = exp.(log_weighted_densities .- logsumexp(log_weighted_densities))

    # log_grads[k] contains the gradients for that k in the mixture
    log_grads = [Gen.logpdf_grad(
            dist.distributions[k], x,
            extract_args_for_component(dist, component_args_flat, k)...)
        for k in 1:dist.K]

    # gradient with respect to x
    log_grads_x = [log_grad[1] for log_grad in log_grads]
    x_grad = if has_output_grad(dist)
        sum(log_grads_x .* relative_weighted_densities)
    else
        nothing
    end

    # gradients with respect to the weights
    weights_grad = exp.(log_densities .- logsumexp(log_weighted_densities))

    # gradients with respect to each argument of each component
    component_arg_grads = Any[]
    cur = 1
    for k in 1:dist.K
        for i in 1:dist.num_args[k]
            if dist.has_argument_grads[cur]
                @assert log_grads[k][i+1] != nothing
                push!(component_arg_grads, relative_weighted_densities[k] * log_grads[k][i+1])
            else
                @assert log_grads[k][i+1] == nothing
                push!(component_arg_grads, nothing)
            end
            cur += 1
        end
    end

    return (x_grad, weights_grad, component_arg_grads...)
end

export HeterogeneousMixture


text:
export DataIterator
export fit!, evaluate
using Random
using ProgressMeter

mutable struct DataIterator
    dataset
    batchsize
    shuffle
	device::Int
end

function DataIterator(dataset; batchsize::Int, shuffle::Bool, device::Int)
    DataIterator(dataset, batchsize, shuffle, device)
end

Base.length(iter::DataIterator) = length(iter.dataset)

function Base.iterate(iter::DataIterator)
	perm = shuffle ? randperm(length(dataset)) : collect(1:length(dataset))
    res = []
    for i = 1:batchsize:length(dataset)
        j = min(i+batchsize-1, length(dataset))
        data = dataset[perm[i:j]]
        push!(res, f(data))
        update!(prog, j)
    end
    res
end
function Base.iterate(iter::DataIterator, state)
	return i >= iter.length ? nothing : (el, (el + it.n, i + 1))
end

function eachbatch(f, dataset, batchsize::Int, shuffle::Bool)
    perm = shuffle ? randperm(length(dataset)) : collect(1:length(dataset))
    prog = Progress(length(dataset))
    res = []
    for i = 1:batchsize:length(dataset)
        j = min(i+batchsize-1, length(dataset))
        data = dataset[perm[i:j]]
        push!(res, f(data))
        update!(prog, j)
    end
    res
end

function fit!(lossfun, model, dataset, opt; batchsize, shuffle, device=-1)
    settrain(true)
    params = collect(Iterators.flatten(parameters.(graphs(model))))
    res = eachbatch(dataset, batchsize, shuffle) do data
        y = lossfun((model,data))
		loss = sum(Array(y.data))
        gradient!(y)
        foreach(opt, params)
		device >= 0 && CUDA.synchronize()
        loss
    end
    loss = sum(res) / length(dataset)
    loss
end

function evaluate(f, model, dataset; batchsize, device=-1)
    settrain(false)
    eachbatch(dataset, batchsize, false) do data
        y = f((model,data))
		device >= 0 && CUDA.synchronize()
		y
    end
end


text:
function drawTagDetection(vis::Visualizer, tagname, Q, T, bTc, bP2t; posename=:test)
  # draw tag triad
  setobject!(vis[currtag], Triad(0.2))
  settransform!(vis[currtag], bTt)

  # draw ray to tag
  v = vis[posename][:lines][tagname]
  geometry = PointCloud(
  GeometryTypes.Point.([bTt.translation[1];0],
                       [bTt.translation[2];0],
                       [0.0;0])
  )
  setobject!(v, LineSegments(geometry, LineBasicMaterial()))
  settransform!(v, Translation(0.0,0,0))

  # draw orientation tail for tag
  v = vis[posename][:orient][tagname]
  geometry = PointCloud(
  GeometryTypes.Point.([0.1;0],
                       [0.0;0],
                       [0.0;0])
  )
  setobject!(v, LineSegments(geometry, LineBasicMaterial()))
  settransform!(v, bP2t) #
  nothing
end

function drawThickLine!(image, startpoint, endpoint, colour, thickness)
    (row, col) = size(image)
    x1 = startpoint.x
    y1 = startpoint.y
    x2 = endpoint.x
    y2 = endpoint.y

    draw!(image, LineSegment(x1,y1,x2,y2), colour)

    if x1 != x2 && (y2-y1)/(x2-x1) < 1
        for tn = 1:thickness
            i = tn ÷ 2

            x1mi = x1-i
            x1pi = x1+i
            y1mi = y1-i
            y1pi = y1+i
            x2mi = x2-i
            x2pi = x2+i
            y2mi = y2-i
            y2pi = y2+i
            #clip to protect bounds
            (x1mi < 1) && (x1mi = 1)
            (y1mi < 1) && (y1mi = 1)
            (x1pi > col) && (x1pi = col)
            (y1pi > row) && (y1pi = row)
            (x2mi < 1) && (x2mi = 1)
            (y2mi < 1) && (y2mi = 1)
            (x2pi > col) && (x2pi = col)
            (y2pi > row) && (y2pi = row)

            iseven(tn) && draw!(image, LineSegment(x1,y1mi,x2,y2mi), colour)
            isodd(tn)  && draw!(image, LineSegment(x1,y1pi,x2,y2pi), colour)
        end
    else
        for tn = 1:thickness
            i = tn ÷ 2

            x1mi = x1-i
            x1pi = x1+i
            y1mi = y1-i
            y1pi = y1+i
            x2mi = x2-i
            x2pi = x2+i
            y2mi = y2-i
            y2pi = y2+i
            #clip to protect bounds
            (x1mi < 1) && (x1mi = 1)
            (y1mi < 1) && (y1mi = 1)
            (x1pi > col) && (x1pi = col)
            (y1pi > row) && (y1pi = row)
            (x2mi < 1) && (x2mi = 1)
            (y2mi < 1) && (y2mi = 1)
            (x2pi > col) && (x2pi = col)
            (y2pi > row) && (y2pi = row)

            iseven(tn) && draw!(image, LineSegment(x1mi,y1,x2mi,y2), colour)
            isodd(tn)  && draw!(image, LineSegment(x1pi,y1,x2pi,y2), colour)
        end
    end
end

function drawTagLine!(imgl, tag_detection)
  # naive
  # ch, cw = round(Int, height/2), round(Int, width/2)
  @warn "Using hard-coded camera calibration parameters in drawTagLine!"
  cw, ch = 330.4173, 196.32587  # from ZED driver config
  focal = 340.97913
  imheight = 376

  tagsize = 0.172

  bearing0 = tag_detection[:bearing][1]

  # convert bearing to position on image
  tag_u_coord = cw + focal*tan(-bearing0)
  tag_u_coord = round(Int, tag_u_coord)

  c1 = Point(tag_u_coord, 1)
  c2 = Point(tag_u_coord, imheight)
  drawThickLine!(imgl, c1, c2, RGB{N0f8}(0.8,0.0,0.4), 10)

  nothing
end

# tan(bear) = u/f

# L = X + DX  =>  X\L = DX

function imageFactor(fg, fsym, im, cfg)
  img = Gray.(im)

  zmeas = getData(getVert(fg, fsym,nt=:fnc)).fnc.usrfnc!.Zpose.z.μ

  # @show IIF.lsf(fg, fsym)
  XX, LL = (KDE.getKDEMax.(IIF.getBelief.(fg, IIF.lsf(fg, fsym)))...)
  @show xyt = se2vee(SE2(XX[1:3]) \ SE2([LL[1:2];0.0]))
  bear= TU.wrapRad(atan2(-xyt[2],xyt[1]) - XX[3])
  # bear = 0.0
  focal= cfg[:intrinsics][:fx]
  cu = cfg[:intrinsics][:cx]
  #
  @show bu = tan(bear)*focal  +  cu
  intens = pdf.(Distributions.Normal(bu,30.0), collect(1:size(img,2)))
  intens ./= maximum(intens)

  zbear= TU.wrapRad(atan2(-zmeas[2],zmeas[1]) - 0.0)

  @show zbu = tan(zbear)*focal  +  cu
  zintens = pdf.(Distributions.Normal(zbu,30.0), collect(1:size(img,2)))
  zintens ./= maximum(zintens)

  for i in 1:round(Int,0.3*size(img, 1))
    img[i, :] .*= intens
  end
  for i in round(Int,0.3*size(img, 1)):size(img,1)
    img[i, :] .*= zintens
  end

  img
end

function drawAllBearingImgs(fg, IMGS, cfg)
  BRIMGS = Dict{Symbol,Dict{Symbol, Any}}()
  XX, LL = ls(fg)

  idx = 0
  for x in XX
    idx += 1
    BRIMGS[x] = Dict{Symbol,Any}()
    for l in LL
      try
        fsym = Symbol("$(x)$(l)f1")
        im = imageFactor(fg, fsym, IMGS[idx], cfg)
        BRIMGS[x][l] = im
      catch Er
        #
      end
    end
  end
  return BRIMGS
end

function genGifLandm(fg, sym, IMGS; show=true, delay=50, dir="")
  BIM = drawAllBearingImgs(fg, IMGS, cfg)
  TEMP = []
  for ss in split.(string.(ls(fg,sym)),string(sym))
    if haskey(BIM, Symbol(ss[1]))
      if haskey(BIM[Symbol(ss[1])], sym)
       push!(TEMP, BIM[Symbol(ss[1])][sym])
      end
    end
  end

  # mkdir("/tmp/animate_tags")
  # Base.rm("/tmp/animate_tags")
  run(`rm -fr /tmp/animate_tags`)
  mkdir("/tmp/animate_tags")
  for i in 1:length(TEMP)
    save("/tmp/animate_tags/im$(i).png", TEMP[i])
  end

  run(`convert -delay $(delay) /tmp/animate_tags/im*.png $(dir)$(sym).gif`)

  !show ? nothing : @async run(`eog $(sym).gif`)
  return "$(sym).gif"
end

function getAllLandmGifs(fg, IMGS; show=false, dir="")
  info("Dropping gifs here $(Base.pwd())")
  for l in ls(fg)[2]
    genGifLandm(fg, l, IMGS, show=show, dir=dir)
  end
  nothing
end

function plotPoseVelAsMax(fg, poserange)

  xx = Symbol.(string.("x",collect(poserange)))

  len = length(xx)
  VV = zeros(len, 2)

  i = 0
  for state in KDE.getKDEMax.(getBelief.(fg, xx))
    i += 1
    VV[i,:] = state[4:5]
  end

  Gadfly.plot(
    Gadfly.layer(x=poserange, y=VV[:,1], Geom.line, Theme(default_color=colorant"red")),
    Gadfly.layer(x=poserange, y=VV[:,2], Geom.line, Theme(default_color=colorant"green"))
  )
end

# # save factor graph for later testing and evaluation
# fg, = loadjld(file=resultsdir*"/racecar_fg_x200_presolve.jld2")
# fg, = loadjld(file=resultsdir*"/racecar_fg_??.jld2")
# @time tree = batchSolve!(fg, N=N, drawpdf=true, show=true, recursive=true)
# 0
#
# IIF.savejld(fg, file=resultsdir*"/racecar_fg_final_resolve.jld2")
# # fgr, = loadjld(file=resultsdir*"/racecar_fg_final_resolve.jld2")
# results2csv(fg; dir=resultsdir, filename="results_resolve.csv")
#
#
# # pl = plotKDE(fg, :x1, levels=1, dims=[1;2]);
#
#
#
# #,xmin=-3,xmax=6,ymin=-5,ymax=2);
# Gadfly.push_theme(:default)
# pl = drawPosesLandms(fg, spscale=0.1, drawhist=false, meanmax=:max);
# # Gadfly.set(:default_theme)
# Gadfly.draw(PNG(joinpath(resultsdir,"images","final.png"),15cm, 10cm),pl);
# Gadfly.draw(SVG(joinpath(resultsdir,"images","final.svg"),15cm, 10cm),pl);
#
# # pl = drawPosesLandms(fg, spscale=0.1, meanmax=:mean) # ,xmin=-3,xmax=3,ymin=-2,ymax=2);
# # Gadfly.draw(PNG(joinpath(resultsdir,"hist_final.png"),15cm, 10cm),pl)
#
# 0

# debugging

# using Profile
#
# using Logging
# global_logger(NullLogger())
# @profile println("test")
# Profile.clear()
# @profile tree1 = IIF.wipeBuildNewTree!(fg)
#
# Juno.profiler()

# vars = lsRear(fg, 5)
#
#
# num_edges(fg.g)
#
# subgraphFromVerts(fg, vars)
#
# ensureAllInitialized!(fg)
# isInitialized(fg, :x220)
#
# ls(fg, :l14)

#


text:
module TestMeasures
    using Test
    using BioEnergeticFoodWebs
    using LinearAlgebra

    # Internal functions
    i = ones(10)
    @test BioEnergeticFoodWebs.shannon(i) ≈ 1.0

    i = zeros(100).+0.001
    i[1] = 1.0
    @test BioEnergeticFoodWebs.shannon(i) ≈ 0.0 atol=0.2

    @test isnan(BioEnergeticFoodWebs.shannon(vec([1.0])))
    @test isnan(BioEnergeticFoodWebs.shannon(vec([-1.0])))

    i = ones(5)
    @test BioEnergeticFoodWebs.coefficient_of_variation(i) == 0.0

    i = collect(range(0.0, stop = 1.0, length = 3))
    @test BioEnergeticFoodWebs.coefficient_of_variation(i) ≈ 1+1/(4*length(i))

    # Test the total biomass thing
    B = Matrix{Float64}(I, 10, 10)
    A = Matrix{Float64}(I, 10, 10)
    fake_params = Dict{Symbol, Any}(:A => A)
    s = Dict{Symbol, Any}(:B => B, :p => fake_params)
    @test total_biomass(s, last=10) == 1.0
    @test_throws AssertionError total_biomass(s, last=1000)
    @test population_biomass(s, last=10)[1] == 0.1

    # Population stability
    @test isnan(population_stability(s, last=1))
    @test population_stability(s, last=2) ≈ -1.59099  atol=0.01
    @test species_richness(s, last=1) == 1.0
    @test species_persistence(s, last=1) == 0.1

    # Foodweb evenness
    B = repeat([0 .5 0 .5], 5)
    s = Dict{Symbol, Any}(:B => B)
    @test foodweb_evenness(s, last = 2) == 1.00

    # Test when the total biomass is 0
    empty_p = Dict{Symbol, Any}(:B => zeros(10, 10))
    @test isnan(population_biomass(empty_p, last=2))
    @test isnan(population_stability(empty_p, last=2))
    @test isnan(total_biomass(empty_p, last=2))
    @test isnan(foodweb_evenness(empty_p, last=2))
    @test isnan(species_richness(empty_p, last = 2))

end

module TestSave
    using Test
    using BioEnergeticFoodWebs
    using LinearAlgebra
    using JLD2
    using JSON

    A = [0 0 0 ; 0 0 0 ; 0 0 0]
    b = rand(3)
    p = model_parameters(A)
    #default variable name and extension
    def_vname = "befwm_simul"
    def_ext = :jld
    #test default arguments
    cd(tempdir())
    s = simulate(p,b)
    # default file name
    def_fname = "befwm_" * string(hash(s)) * ".jld2"
    save(s, as = def_ext)

    # Test if the file is saved (under the default name)
    @test isfile(def_fname)
    # Test if the content is the same
    #@load def_fname sim
    #@test sim == s
    rm(def_fname)

    # Test with .json
    ext = :json
    fname = "befwm_" * string(hash(s)) * ".json"
    save(s, as = ext)
    @test isfile(fname)
    rm(fname)

    # Test with .JSON
    ext = :JSON
    save(s, as = ext)
    @test isfile(fname)
    rm(fname)

end

module TestPGrowth
    using Test
    using BioEnergeticFoodWebs
    import Statistics.mean, Statistics.std

    A = zeros(Int64, (3, 3))
    n = ones(3)

    # System productivity
    p = model_parameters(A, productivity=:system)
    s = simulate(p, n, start=0, stop=500)
    pgall = producer_growth(s, last = 100, out_type = :all)
    pgstd = producer_growth(s, last = 100, out_type = :std)
    meanpg = mean(pgall, dims = 1)
    stdpg = std(pgall, dims = 1)
    pg = producer_growth(s, last = 100, out_type = :mean)
    @test pg[1] == meanpg[1]
    @test pg[2] == meanpg[2]
    @test pg[3] == meanpg[3]
    @test pgstd[1] == stdpg[1]
    @test pgstd[2] == stdpg[2]
    @test pgstd[3] == stdpg[3]
    K = p[:K][1]
    B = population_biomass(s, last = 100)
    @test B[1] == B[2] == B[3]
    @test B[1] ≈ K/3 atol = 1e-4
    pg_calc = 0
    @test pg[1] == pg[2] == pg[3]
    @test pg[1] ≈ pg_calc atol = 1e-4

    # Nutrient intake
    A = [0 0 ; 0 0]
    b0 = [0.5, 0.5]
    c0 = [2.0, 2.0]
    k1 = [0.2, 0.1]
    k2 = [0.1, 0.2]
    p = model_parameters(A, productivity = :nutrients, K1 = k1, K2 = k2)
    s = simulate(p, b0, start=0, stop=500)
    pg = producer_growth(s, last = 100, out_type = :mean)
    @test pg[1] ≈ 0 atol = 1e-4
    @test pg[2] ≈ 0 atol = 1e-4

end


text:
# File: normalplot.jl
# Author:
# Date: 29-June-2020

# import XL file
[~,dates,returns] = getXLData('hfma_matlab_data.xlsx','CTA Index');

figure; # create figure
normplot(returns); # plot normal probability

#Plotting univariate distributions
using Distributions
using StatsPlots

norm = fit(Normal, rand(100))
plot(norm)

# plotting a multivariate normal Distribution

mvnorm = fit(MvNormal, [rand(0.0:100.0, 100) rand(0.0:100.0, 100)]')

Z = [pdf(mvnorm,[i,j]) for i in 0:100, j in 0:100]
plot(0:100,0:100,Z,st=:surface)

# I was able to hack together something that looks decent.
# Below is the code as well as the output for some data I’m working with:
X = 3
Y = 5
x = 0:maximum(X)
y = 0:(maximum(Y)/length(x))+1:maximum(Y) # because x and y are of different lengths
z = [pdf(mvnorm, [i, j]) for i in x, j in y]
using PyPlot
pyplot()

plot(x, y, z, linetype=:surface, legend=false, color=:blues)
plot!(x, y, z, linetype=:wireframe, legend=false, color=:black, width=0.1)


text:
mutable struct Mesh
    type  :: Symbol
    n     ::Ref{Int64}
    x     ::Vector{Float64}
    Mesh() = new( )
end

# Create structure with default missing fields
mesh = Mesh()

mesh.n  = 101
mesh.x  = [1.1, 1.5, 1.9]

println( mesh.n )
println( mesh.x )


text:
# imports
using FMI
using FMIFlux
using FMIZoo
using Flux
using DifferentialEquations: Tsit5
import Plots

# set seed
import Random
Random.seed!(1234);

tStart = 0.0
tStep = 0.01
tStop = 5.0
tSave = tStart:tStep:tStop

referenceFMU = fmiLoad("SpringPendulumExtForce1D", "Dymola", "2022x")
fmiInfo(referenceFMU)

param = Dict("mass_s0" => 1.3, "mass.v" => 0.0)   # increase amplitude, invert phase
vrs = ["mass.s", "mass.v", "mass.a"]
referenceSimData = fmiSimulate(referenceFMU, tStart, tStop; parameters=param, recordValues=vrs, saveat=tSave)
fmiPlot(referenceSimData)

posReference = fmi2GetSolutionValue(referenceSimData, vrs[1])
velReference = fmi2GetSolutionValue(referenceSimData, vrs[2])
accReference = fmi2GetSolutionValue(referenceSimData, vrs[3])

defaultFMU = referenceFMU
param = Dict("mass_s0" => 0.5, "mass.v" => 0.0)

defaultSimData = fmiSimulate(defaultFMU, tStart, tStop; parameters=param, recordValues=vrs, saveat=tSave)
fmiPlot(defaultSimData)

posDefault = fmi2GetSolutionValue(defaultSimData, vrs[1])
velDefault = fmi2GetSolutionValue(defaultSimData, vrs[2])
accDefault = fmi2GetSolutionValue(defaultSimData, vrs[3])

function extForce(t)
    return [0.0]
end

# loss function for training
function lossSum()
    solution = csNeuralFMU(extForce, tStep)

    accNet = fmi2GetSolutionValue(solution, 1; isIndex=true)

    Flux.Losses.mse(accReference, accNet)
end

# callback function for training
global counter = 0
function callb()
    global counter += 1

    if counter % 20 == 1
        avgLoss = lossSum()
        @info "Loss [$counter]: $(round(avgLoss, digits=5))"
    end
end

# NeuralFMU setup
numInputs = length(defaultFMU.modelDescription.inputValueReferences)
numOutputs = length(defaultFMU.modelDescription.outputValueReferences)

net = Chain(inputs -> fmiInputDoStepCSOutput(defaultFMU, tStep, inputs),
            Dense(numOutputs, 16, tanh),
            Dense(16, 16, tanh),
            Dense(16, numOutputs))

csNeuralFMU = CS_NeuralFMU(defaultFMU, net, (tStart, tStop); saveat=tSave);

solutionBefore = csNeuralFMU(extForce, tStep)
accNeuralFMU = fmi2GetSolutionValue(solutionBefore, 1; isIndex=true)
Plots.plot(tSave, accNeuralFMU, label="acc CS-NeuralFMU", linewidth=2)

# train
paramsNet = Flux.params(csNeuralFMU)

optim = ADAM()
Flux.train!(lossSum, paramsNet, Iterators.repeated((), 300), optim; cb=callb)

# plot results mass.a
solutionAfter = csNeuralFMU(extForce, tStep)

fig = Plots.plot(xlabel="t [s]", ylabel="mass acceleration [m/s^2]", linewidth=2,
                 xtickfontsize=12, ytickfontsize=12,
                 xguidefontsize=12, yguidefontsize=12,
                 legendfontsize=8, legend=:topright)

accNeuralFMU = fmi2GetSolutionValue(solutionAfter, 1; isIndex=true)

Plots.plot!(fig, tSave, accDefault, label="defaultFMU", linewidth=2)
Plots.plot!(fig, tSave, accReference, label="referenceFMU", linewidth=2)
Plots.plot!(fig, tSave, accNeuralFMU, label="CS-NeuralFMU (300 eps.)", linewidth=2)
fig

fmiUnload(defaultFMU)


text:
# ps4000closeunit.jl

export ps4000closeunit

function ps4000closeunit(handle::UInt16)
  ps_status = ccall((:ps4000CloseUnit, ps4000driver),
                  PicoStatus,
                  (UInt16,),
                  handle)
  checkps_status(ps_status)
  return nothing
end


text:
@inbounds usage within string code #35842 #Open @KristofferC opened this issue on May 11, 2020 · 0 comments
#https://github.com/JuliaLang/julia/issues/35842
#=
Is there a True Significance for @inbounds? Really?
=#
#=
Comments
@KristofferC

Sponsor Member KristofferC commented on May 11, 2020 •

The string code is filled with @inbounds, @propagate_inbounds, @checkbounds etc but it isn't (at least to me) completely clear what these really mean and it doesn't really feel documented.

Let's take string indexing, here @inbounds is only used to turn off checking that you are not outside the bounds of the string but you still need to index on a valid index:
=#

f(a, i) = @inbounds a[i]
f(0, 1)
f(1, 10)
#f (generic function with 1 method)

f("aα", 2) # if string, it checks it , at a valid index
f("aα", 3)
# ERROR: StringIndexError("aα", 3) #erronr annotations been update meaningfully # `valid nearby index` [2]=> `a`
#useful in
s = "aα"
@time begin
    try
        @inbounds for i in size(f(s, size(f(s, 1))))
        end
    catch
        if (!OutOfMemoryError)
            println("caught ;)")
        end
    end
end

#=
However, you can also add an @inbounds to the SubString constructor (how would you even know that) which then turns off checking for a valid index:
=#
SubString("fα", 1, 3)  #same problem index 3 higher than max size of string , if you add

SubString("fα", 1, 2) #successfully returns fa

#=
ERROR: StringIndexError("fα", 3) =#

g(a, i, j) = @inbounds SubString(a, i, j) #works
# g (generic function with 1 method)

g("fα", 1, 2) #valid around max length (string)
# no error #it doesn't works
# "fα"
#=
So to me,
`it feels a bit inconsistent``.
Many string functions that take some kind of string index use @checkbounds but
what do you actually promise when you use @inbounds and
How do you find out which functions
`that are useful to pass it to`.

I can only comment on changes right
I guess there are more things yet to be uncovered
#

=#

#---testing

#-- speed
using BenchmarkTools

#=ensures value stays withing bounds  =#
f(a, i) = @inbounds a[i]
@benchmark f(1, 8^2) # rand(1,1) is unpredictable
f(a, i) = a[i] #redo without
@benchmark f(1, 8^2) #ERROR bounds error!
#=
 Memory estimate: 0 bytes, allocs estimate: 0.

BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
 Range (min … max):  1.300 ns … 21.300 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     1.400 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.505 ns ±  1.169 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

                              █                            ▁
  ▅▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▂
  1.3 ns         Histogram: frequency by time         1.5 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

 =#
f(1, 10) #bounds Error

#without @inbounds

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
using BlockBandedMatrices

@testset "Adj/Trans" begin
    A = BandedBlockBandedMatrix(randn(ComplexF64,10,14), 1:4,2:5, (1,2), (2,1))

    @test A'[Block(1,1)] == A[Block(1,1)]'
    @test A'[Block(2,3)] == A[Block(3,2)]'
    @test transpose(A)[Block(1,1)] == transpose(A[Block(1,1)])
    @test transpose(A)[Block(2,3)] == transpose(A[Block(3,2)])

    @test blockbandwidths(A') == blockbandwidths(transpose(A)) == (2,1)
    @test subblockbandwidths(A') == subblockbandwidths(transpose(A)) == (1,2)

    @test BandedBlockBandedMatrix(A') == A'
    @test BandedBlockBandedMatrix(transpose(A)) == transpose(A)
end


text:
function is_path_valid(path_to_file::String)::PSResult{Bool}

    # the config file should be a json file, and should be reachable -
    # TODO: fill me in ...
    return PSResult{Bool}(true)
end

function is_string_empty(raw_string::String)::PSResult{Bool}

    # if we have an empty string - return true
    if isempty(raw_string) == true
        return PSResult{Bool}(true)
    end

    # default Return is flase -
    return PSResult{Bool}(false)
end

function check_missing_api_key(user_model::PSUserModel)::(Union{T, Nothing} where T<:Any)

    # do we have the alpha_vantage_api_key -
    if (hasfield(PSUserModel, :alphavantage_api_key) == false)
        # throw -
        return PSResult{PSError}(PSError("user model is missing API key information"))
    end

    # get the key -
    api_key = user_model.alphavantage_api_key

    # check -
    if (isempty(api_key) == true)

        # formulate an error message -
        error_message = "the API key is empty in the user model"

        # throw -
        return PSResult{PSError}(PSError(error_message))
    end

    # ok -
    return nothing
end

function check_missing_symbol(stock_symbol::String)::(Union{T, Nothing} where T<:Any)

    if (isempty(stock_symbol) == true)

        # formulate an error message -
        error_message = "missing stock symbol"

        # throw -
        return PSResult{PSError}(PSError(error_message))
    end

    # return nothing -
    return nothing
end

function check_json_api_return_data(api_call_raw_data::String)::(Union{T, Nothing} where T<:Any)

    # well formed JSON?
    if is_valid_json(api_call_raw_data).value == false
        return PSError("invalid JSON $(api_call_raw_data)")
    end

    # need to check to see if legit data is coming back from the service -
    api_data_dictionary = JSON.parse(api_call_raw_data)
    if (haskey(api_data_dictionary,"Error Message") == true)

        # grab the error mesage -
        error_message = api_data_dictionary["Error Message"]

        # throw -
        return PSResult{PSError}(PSError(error_message))
    end

    # need to check - are we hitting the API call limit?
    if (haskey(api_data_dictionary,"Note") == true)

        # grab the error mesage -
        error_message = api_data_dictionary["Note"]

        # throw -
        return PSResult{PSError}(PSError(error_message))
    end

    # default -
    return nothing
end

function check_user_model(user_model::PSUserModel)::(Union{T, Nothing} where T<:Any)
    return PSResult{Bool}(true)
end

"""
    _is_valid_json(raw_string)->Bool

Checks to see if the input `string` is a valid JSON structure.
Returns `true` indicating valid JSON, `false` otherwise.
"""
function is_valid_json(raw_string::String)::PSResult{Bool}

    # check: do we have an empty string?
    if (is_string_empty(raw_string).value == true)
        return PSResult{Bool}(false)
    end

    # otherwise, to check to see if the string is valid JSON, try to
    # parse it.
    try
        JSON.parse(raw_string)
        return PSResult{Bool}(true)
    catch
        return PSResult{Bool}(false)
    end
end


text:
module REPLPromptTimes

function __init__()
	println("Paste the following into your Julia REPL:")
	println("""include(joinpath(dirname(pathof(REPLPromptTimes)), "include.jl"))""")
end

end #module


text:
module DEBTrait

using Parameters, FieldMetadata
import FieldMetadata: @default, @description, @units, @bounds, @logscaled, @flattenable, @plottable, @selectable,
                      default, description, units, bounds, logscaled, flattenable, plottable, selectable
using Unitful
Unitful.register(@__MODULE__);
@unit bp "bp" Basepair 1u"m" true;
using Unitful: °C, K, Pa, kPa, MPa, J, kJ, W, L, g, kg, cm, m, s, hr, d, mol, mmol, μmol, σ, R

using Roots
using LinearAlgebra

include("metabolism.jl")
include("cell_composition.jl")
include("assimilation.jl")
include("supeca.jl")
include("soil_properties.jl")
include("ThermoStoichWizard.jl")
include("models.jl")

export AbstractMetabolism, PMetabolismC, PMetabolismCN, growth!, growth_production
export AbstractGenomicData, BaseGenome, PCellComposition
export AbstractMonomer, AbstractAssimilation, BaseMonomer, PAssimilation
export AbstractSoil, BaseSoil
export AbstractParameterSpace, SetupBatchC, ParamsBatchC, rhs_BatchC!

# Write your package code here.

end


text:
include("dependencies_for_runtests.jl")

using Oceananigans.Solvers: solve!
using Statistics

function identity_operator!(b, x)
    parent(b) .= parent(x)
    return nothing
end

function run_identity_operator_test(grid)
    b = CenterField(grid)
    solver = PreconditionedConjugateGradientSolver(identity_operator!, template_field = b, reltol=0, abstol=10*sqrt(eps(eltype(grid))))
    initial_guess = solution = similar(b)
    set!(initial_guess, (x, y, z) -> rand())

    solve!(initial_guess, solver, b)

    @test norm(solution) .< solver.abstol
end

function run_poisson_equation_test(grid)
    arch = architecture(grid)
    # Solve ∇²ϕ = r
    ϕ_truth = CenterField(grid)

    # Initialize zero-mean "truth" solution with random numbers
    set!(ϕ_truth, (x, y, z) -> rand())
    parent(ϕ_truth) .-= mean(ϕ_truth)
    fill_halo_regions!(ϕ_truth)

    # Calculate Laplacian of "truth"
    ∇²ϕ = r = CenterField(grid)
    compute_∇²!(∇²ϕ, ϕ_truth, arch, grid)

    solver = PreconditionedConjugateGradientSolver(compute_∇²!, template_field=ϕ_truth, reltol=eps(eltype(grid)))

    # Solve Poisson equation
    ϕ_solution = CenterField(grid)
    solve!(ϕ_solution, solver, r, arch, grid)
    fill_halo_regions!(ϕ_solution)

    # Diagnose Laplacian of solution
    ∇²ϕ_solution = CenterField(grid)
    compute_∇²!(∇²ϕ_solution, ϕ_solution, arch, grid)

    # Test
    extrema_tolerance = 1e-12
    std_tolerance = 1e-13

    CUDA.@allowscalar begin
        @test minimum(abs, interior(∇²ϕ_solution) .- interior(∇²ϕ)) < extrema_tolerance
        @test maximum(abs, interior(∇²ϕ_solution) .- interior(∇²ϕ)) < extrema_tolerance
        @test          std(interior(∇²ϕ_solution) .- interior(∇²ϕ)) < std_tolerance

        @test   minimum(abs, interior(ϕ_solution) .- interior(ϕ_truth)) < extrema_tolerance
        @test   maximum(abs, interior(ϕ_solution) .- interior(ϕ_truth)) < extrema_tolerance
        @test            std(interior(ϕ_solution) .- interior(ϕ_truth)) < std_tolerance
    end

    return nothing
end

@testset "PreconditionedConjugateGradientSolver" begin
    for arch in archs
        @info "Testing PreconditionedConjugateGradientSolver [$(typeof(arch))]..."
        grid = RectilinearGrid(arch, size=(4, 8, 4), extent=(1, 3, 1))
        run_identity_operator_test(grid)
        run_poisson_equation_test(grid)
    end
end


text:
#!/usr/bin/env julia

using ArgParse: ArgParseSettings, @add_arg_table, parse_args, parse_item
import ArgParse

using AkazeFeatures
# include("akaze-config.jl")
# include("akaze-descriptor.jl")
# include("akaze.jl")
# include("fed.jl")
# include("nonlinear-diffusion.jl")
# include("utils.jl")

function ArgParse.parse_item(::Type{DIFFUSIVITY_TYPE}, str::AbstractString)
    return DIFFUSIVITY_TYPE(parse(Int, str))
end

s = ArgParseSettings()
@add_arg_table s begin
    "--nsublevels"
    help = "Number of sublevels per scale level"
    arg_type = Int
    default = 4

    "--omax"
    help = "Maximum octave evolution of the image 2^sigma (coarsest scale sigma units)"
    arg_type = Int
    default = 4

    "--soffset"
    help="Base scale offset (sigma units)"
    arg_type = Float32
    default = 1.6f0

    "--derivative_factor"
    help="Factor for the multiscale derivatives"
    arg_type = Float32
    default = 1.5f0

    "--sderivatives"
    help="Smoothing factor for the derivatives"
    arg_type=Float32
    default = 1f0

    "--diffusivity"
    # arg_type=Int
    # default = Int(PM_G2)
    arg_type=DIFFUSIVITY_TYPE
    default = PM_G2
    help="Diffusivity type"

    "--dthreshold"
    help="Detector response threshold to accept point"
    arg_type=Float32
    default = 1f-3
    # default = 3f-5

    "--min_dthreshold"
    help="Minimum detector threshold to accept a point"
    arg_type=Float32
    default = 1f-5

    "--output"
    help="output filename"
    arg_type=String
    default = "keypoints.txt"

    "imagename"
    help = "Image file to be analyzed"
    required = true
end

parsed_args = parse_args(ARGS, s)

img = load_image_as_grayscale(parsed_args["imagename"])

img_height, img_width = size(img)

opt = AKAZEOptions(
    omin = 3,
    omax = parsed_args["omax"],
    nsublevels = parsed_args["nsublevels"],
    img_width = img_width,
    img_height = img_height,
    diffusivity = parsed_args["diffusivity"],
    dthreshold = parsed_args["dthreshold"],
)
akaze = AKAZE(opt)

Create_Nonlinear_Scale_Space(akaze, img)
kpts = Feature_Detection(akaze)
Compute_Main_Orientation.([akaze], kpts)
println("Extracted $(length(kpts)) points")
desc = Compute_Descriptors(akaze, kpts)

open(parsed_args["output"],"w") do io
    for row in dump_keypoints_text(kpts, desc)
        println(io, row)
    end
end


text:
using PiAnalysis
using Test

# Check that chud() is calculating pi accurately by comparing the
# result of chud() to that of Julia's pi then check piAccurate() by
# comparing the difference to the result of piAccurate()
#
# Either they're both right or they're both wrong

@testset "Test of chud() and piAccurate()" begin

  # Check the accuracy of chud() and pi from 1024 bits
  # to 102400 bits

  for i in 1:10

    # Set the precision

    setprecision(BigFloat,1024*i)

    # Compute the difference between pi and chud() and then
    # compare to piAccurate()

    p1 = BigFloat(pi)
    p2 = chud(1024*i)
    diff1 = abs(p1-p2)
    diff2 = piAccurate(1024*i)

    a_temp = 10^(BigInt(300*i))
    a = BigFloat(1/a_temp)

    @test diff1 < a
    @test diff2 == diff1

  end

end

# Check that the first 10 digits of numDig() are correct and then
# check that the results of numDig() and digCompare() don't vary
# outside statistical probability

@testset "Test of numDig()" begin

  # Check that the first ten digits are correct

  oneTo10_temp = numDig(1,10,1024)
  oneTo10 = oneTo10_temp[:,2]
  correctOneTo10 = [0,2,1,1,1,3,1,0,0,1]

  @test oneTo10-correctOneTo10 == [0,0,0,0,0,0,0,0,0,0]

  # Check that the results of numDig() and digCompare() don't vary
  # outside statistical probability
  #
  # The digits of pi are randomly distributed, so the count for each
  # digit shouldn't vary outside (+/-)sqrt(N) of 1/10th of N
  # where N is the length of the interval. By the same reasoning,
  # the result for each digit in digCompare() should be within
  # (+/-)sqrt(2)*sqrt(N) since errors add in quadrature.

  for i in 1:100

    ary1 = numDig(i,5*i,2048)
    ary2 = digCompare(i,2*i,2*i+1,3*i+1,2048)
    intvl = 5*i-i+1
    intvl2 = 2*i-i+1

    for j in 1:10
      @test ary1[j,2] < (intvl/10) + sqrt(intvl) # Test of numDig()
      @test ary1[j,2] > (intvl/10) - sqrt(intvl)
      @test ary2[j,2] < sqrt(2) * sqrt(intvl2)         # Test of digCompare()
      @test ary2[j,2] > -sqrt(2) * sqrt(intvl2)
    end

  end
end


text:
@testset "Power Series Differentiation" begin

    T, t = Nemo.PowerSeriesRing(Nemo.QQ, 300, "t"; model=:capped_absolute)

    @test ps_diff(zero(T)) == zero(T)

    @test ps_diff(one(T)) == zero(T)

    @test ps_diff(gen(T)) == one(T)

    @test ps_diff(t^3 - t^4 * (1 // 8) + 5 * t^10) == 3 * t^2 - t^3 * (1 // 2) + 50 * t^9

    @test ps_diff( (1 + t)^1000 ) == truncate(1000 * (1 + t)^999, 299)

end


text:
module priors

using ..gpflow
import ..gpflow: instantiate!

export
    Gaussian,
    LogNormal,
    Gamma,
    Laplace,
    Uniform,
    instantiate!,
    logp,
    sample

function logp(obj::T, x) where T <: Prior
    return  obj.o.logp(x)
end

function sample(obj::T; shape=(1,)) where T <: Prior
    return  obj.o.sample(;shape=shape)
end

mutable struct Gaussian{T1,T2} <: Prior
    mu::T1
    var::T2
    o::Union{PyObject, Nothing}
end

function Gaussian(mu, var)
    out = Gaussian(mu, var, nothing)
    instantiate!(out)
    out
end

function instantiate!(o::Union{Gaussian, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.priors.Gaussian(o.mu, o.var)
    return o.o
end

mutable struct LogNormal{T1,T2} <: Prior
    mu::T1
    var::T2
    o::Union{PyObject, Nothing}
end

function LogNormal(mu, var)
    out = LogNormal(mu, var, nothing)
    instantiate!(out)
    out
end

function instantiate!(o::Union{LogNormal, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.priors.LogNormal(o.mu, o.var)
    return o.o
end

mutable struct Gamma{T1,T2} <: Prior
    shape::T1
    scale::T2
    o::Union{PyObject, Nothing}
end

function Gamma(shape, scale)
    out = Gamma(shape, scale, nothing)
    instantiate!(out)
    out
end

function instantiate!(o::Union{Gamma, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.priors.Gamma(o.shape, o.scale)
    return o.o
end

mutable struct Laplace{T1,T2} <: Prior
    mu::T1
    sigma::T2
    o::Union{PyObject, Nothing}
end

function Laplace(mu, sigma)
    out = Laplace(mu, sigma, nothing)
    instantiate!(out)
    out
end

function instantiate!(o::Union{Laplace, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.priors.Laplace(o.mu, o.sigma)
    return o.o
end

mutable struct Uniform{T1,T2} <: Prior
    lower::T1
    upper::T2
    o::Union{PyObject, Nothing}
end

function Uniform(;lower=0.0, upper=1.0)
    out = Uniform(lower, upper, nothing)
    instantiate!(out)
    out
end

function instantiate!(o::Union{Uniform, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.priors.Uniform(;lower=o.lower, upper=o.upper)
    return o.o
end

end # module


text:
"""
    PyList{T=Py}([x])

Wraps the Python list `x` (or anything satisfying the sequence interface) as an `AbstractVector{T}`.

If `x` is not a Python object, it is converted to one using `pylist`.
"""
struct PyList{T} <: AbstractVector{T}
    py :: Py
    PyList{T}(::Val{:new}, py::Py) where {T} = new{T}(py)
end
export PyList

PyList{T}(x=pylist()) where {T} = PyList{T}(Val(:new), ispy(x) ? Py(x) : pylist(x))
PyList(x=pylist()) = PyList{Py}(x)

ispy(::PyList) = true
getpy(x::PyList) = x.py
pydel!(x::PyList) = pydel!(x.py)

pyconvert_rule_sequence(::Type{T}, x::Py, ::Type{PyList{V}}=Utils._type_ub(T)) where {T<:PyList,V} =
    if PyList{Py} <: T
        pyconvert_return(PyList{Py}(x))
    else
        pyconvert_return(PyList{V}(x))
    end

Base.length(x::PyList) = Int(pylen(x))

Base.size(x::PyList) = (length(x),)

Base.@propagate_inbounds function Base.getindex(x::PyList{T}, i::Int) where {T}
    @boundscheck checkbounds(x, i)
    return pyconvert_and_del(T, @py x[@jl(i-1)])
end

Base.@propagate_inbounds function Base.setindex!(x::PyList{T}, v, i::Int) where {T}
    @boundscheck checkbounds(x, i)
    pysetitem(x, i-1, convert(T, v))
    return x
end

Base.@propagate_inbounds function Base.insert!(x::PyList{T}, i::Integer, v) where {T}
    @boundscheck (i==length(x)+1 || checkbounds(x, i))
    pydel!(@py x.insert(@jl(i-1), @jl(convert(T, v))))
    return x
end

function Base.push!(x::PyList{T}, v) where {T}
    pydel!(@py x.append(@jl(convert(T, v))))
    return x
end

function Base.pushfirst!(x::PyList, v)
    return @inbounds Base.insert!(x, 1, v)
end

function Base.append!(x::PyList, vs)
    for v in vs
        push!(x, v)
    end
    return x
end

function Base.push!(x::PyList, v1, v2, vs...)
    push!(x, v1)
    push!(x, v2, vs...)
end

Base.@propagate_inbounds function Base.pop!(x::PyList{T}) where {T}
    @boundscheck (isempty(x) && throw(BoundsError(x)))
    return pyconvert_and_del(T, @py x.pop())
end

Base.@propagate_inbounds function Base.popat!(x::PyList{T}, i::Integer) where {T}
    @boundscheck checkbounds(x, i)
    return pyconvert_and_del(T, @py x.pop(@jl(i-1)))
end

Base.@propagate_inbounds Base.popfirst!(x::PyList) = popat!(x, 1)

function Base.reverse!(x::PyList)
    pydel!(@py x.reverse())
    return x
end

function Base.empty!(x::PyList)
    pydel!(@py x.clear())
    return x
end

function Base.copy(x::PyList{T}) where {T}
    PyList{T}(Val(:new), @py x.copy())
end


text:
module HierGaussLikelihood

using AutoDiff.Backward

# This is a standard hierarchical likelihood.  We imagine that we have a set of
# objects, each of which has a parameter `x`.  Collectively, the `x` are
# distributed according to a Gaussian with mean `mu` and s.d. `sigma`.  We have
# some number of observations of each x, `xobs` with a known uncertainty
# (assumed Gaussian).  We want to be able to take gradients with respect to
# `mu`, `sigma`, and the `x`.

function make_hglikelihood(xobs, sigmaobs)
    nx, no = size(xobs)

    function loglike(mu, sigma, xs...)
        ll = sum(-0.5.*(xs.-mu).*(xs.-mu)./(sigma*sigma)) - float(nx)*log(sigma)

        ll2 = sum(-0.5.*(xobs .- xs).*(xobs .- mu)./(sigmaobs.*sigmaobs) .- log.(sigmaobs))

        ll + ll2
    end

    loglike
end

nx = 1000
no = 5

mu = 0.0
sigma = 1.0

xs = mu .+ sigma.*randn(nx)
sigmaobs = 1.0 .+ 1.0.*rand(nx, no)
xobs = xs .+ sigmaobs.*randn(nx, no)

loglike = make_hglikelihood(xobs, sigmaobs)
g = gradient(loglike)

g(mu, sigma, xs...)

end


text:
# make_polaron.jl

# Example: MAPI: make_polaron(4.5, 24.1, 2.25e12, 0.12)

function make_polaron(ϵ_optic, ϵ_static, phonon_freq, m_eff; temp = 300.0, efield_freq = 0.0, coupling = 0.0, verbose = true)

    # Collect data.
    ω = 2 * π * phonon_freq
    Ω = efield_freq
    T = temp
    if coupling == 0.0
        α = frohlich_α(ϵ_optic, ϵ_static, phonon_freq, m_eff)
    elseif coupling > 0.0
        α = coupling
    end

    # Prepare empty arrays for different temperatures.
    β = [] # Reduced thermodynamic beta (unitless)
    v = [] # Variational parameter v (1 / s, Hz)
    w = [] # Variational parameter w (1 / s, Hz)
    κ = [] # Spring constant (kg / s^2)
    M = [] # Fictitious mass (kg)
    F = [] # Free energy (meV)
    μ = [] # Mobility (cm^2 / Vs)
    Γ = [] # Optical absorption (1 / cm)

    # Initialise variation parameters.
    v_t, w_t = 0.0, 0.0
    if verbose
        print("\n")
    end

    for t in T # Iterate over temperatures.
        if verbose
            print("\e[2K", "Working on temperature: $(t) K / $(temp[end]) K.\n")
        end

        if t == 0.0 # If T = 0
            append!(β, Inf)  # set β = Inf

            # Evaluate variational parameters.
            v_t, w_t = variation(α; v = v_t, w = w_t)
            append!(v, v_t)
            append!(w, w_t)

            # Evaluate fictitious particle mass and spring constant.
            κ_t = (v_t^2 - w_t^2) # spring constant
            M_t = ((v_t^2 - w_t^2) / w_t^2) # mass
            append!(κ, κ_t)
            append!(M, M_t)

            # Evaluate free energy at zero temperature. NB: Enthalpy.
            F_t = free_energy(v_t, w_t, α) * 1000 * ħ * ω / eV
            append!(F, F_t)

            # Broadcast data.
            if verbose
                println("\e[2K", "α: ", round(α, digits = 3), " | β: ", Inf, " | v: ", round(v_t, digits = 3), " s^-1 | w: ", round(w_t, digits = 3), " s^-1")
                println("\e[2K", "κ: ", round(κ_t, digits = 3), " m_e kg/s^2 | M: ", round(M_t, digits = 3), " m_e kg")
                println("\e[2K", "Free Energy (Enthalpy): ", round(F_t, digits = 3), " meV")
            end

            # Prepare empty arrays for different frequencies.
            μ_t = [] # Mobilities
            Γ_t = [] # Optical absorptions

            for f in Ω # Iterate over frequencies

                if f == 0.0 # If Ω = 0 at T = 0

                    # Evaluate AC mobilities. NB: Ω = 0 is DC mobility.
                    μ_f = polaron_mobility_dc(1e3, α, v_t, w_t)
                    append!(μ_t, μ_f)

                    # Evaluate frequency-dependent optical absorptions.
                    Γ_f = optical_absorption_dc(1e3, α, v_t, w_t)
                    append!(Γ_t, Γ_f)

                    # Broadcast data.
                    if verbose
                        println("\e[2K", "Working on Frequency: $(f) Hz / $(efield_freq[end]) Hz")
                        println("\e[2K", "DC Mobility: ", round(μ_f, digits = 3), " cm^2/Vs")
                        println("\e[2K", "Optical absorption: ", round(Γ_f, digits = 3), " cm^-1")
                    end

                elseif f > 0.0 # If Ω > 0 at T = 0

                    # Evaluate AC mobilities.
                    μ_f = 100^2 * eV * polaron_mobility_ac(f, α, v_t, w_t) / (ω * m_eff * m_e)
                    append!(μ_t, μ_f)

                    # Evaluate optical absorptions.
                    Γ_f = optical_absorption_ac(f, α, v_t, w_t) / (100 * c * ϵ_0 * sqrt(ϵ_optic))
                    append!(Γ_t, Γ_f)

                    # Broadcast data.
                    if verbose
                        println("\e[2K", "Working on Frequency: $(f) Hz / $(efield_freq[end]) Hz")
                        println("\e[2K", "AC Mobility: ", round(μ_f, digits = 3), " cm^2/Vs")
                        println("\e[2K", "Optical absorption: ", round(Γ_f, digits = 3), " cm^-1")
                    end
                end

                print("\033[F"^3) # Move cursor back
            end

            # Add freq-dependent mobilities and absorptions for T = 0.
            append!(μ, [μ_t])
            append!(Γ, [Γ_t])
            if verbose
                print("\033[F"^4) # Move cursor back
            end

        elseif t > 0.0 # If T > 0

            # Evaluate reduced thermodynamic beta.
            β_t = ħ * ω / (k_B * t)
            append!(β, β_t)

            # Evaluate variational parameters.
            v_t, w_t = variation(α, β_t; v = v_t, w = w_t)
            append!(v, v_t)
            append!(w, w_t)

            # Evaluate fictitious particle mass and spring constant.
            κ_t = (v_t^2 - w_t^2) # spring constant
            M_t = ((v_t^2 - w_t^2) / w_t^2) # mass
            append!(κ, κ_t)
            append!(M, M_t)

            # Evaluate free energy at finite temperature.
            F_t = free_energy(v_t, w_t, α, β_t) * 1000 * ħ * ω / eV
            append!(F, F_t)

            # Broadcast data.
            if verbose
                println("\e[2K", "α: ", round(α, digits = 3), " | β: ", round(β_t, digits = 3), " | v: ", round(v_t, digits = 3), " s^-1 | w: ", round(w_t, digits = 3), " s^-1")
                println("\e[2K", "κ: ", round(κ_t, digits = 3), " m_e kg/s^2 | M: ", round(M_t, digits = 3), " m_e kg")
                println("\e[2K", "Free Energy: ", round(F_t, digits = 3), " meV")
            end

            # Prepare empty arrays for different frequencies.
            μ_t = [] # Mobilities
            Γ_t = [] # Optical absorptions

            for f in Ω # Iterate over frequencies

                if f == 0.0 # If Ω = 0 at T > 0

                    # Evaluate DC mobility.
                    μ_f = 100^2 * eV * polaron_mobility_dc(β_t, α, v_t, w_t) / (ω * m_eff * m_e)
                    append!(μ_t, μ_f)

                    # Evaluate DC optical absorption.
                    Γ_f = optical_absorption_dc(β_t, α, v_t, w_t) / (100 * c * ϵ_0 * sqrt(ϵ_optic))
                    append!(Γ_t, Γ_f)

                    # Broadcast data.
                    if verbose
                        println("\e[2K", "Working on Frequency: $(f) Hz / $(efield_freq[end]) Hz")
                        println("\e[2K", "DC Mobility: ", round(μ_f, digits = 3), " cm^2/Vs")
                        println("\e[2K", "Optical absorption: ", round(Γ_f, digits = 3), " cm^-1")
                    end

                elseif f > 0.0 # If Ω > 0 at T > 0

                    # Evaluate AC mobilities.
                    μ_f = 100^2 * eV * polaron_mobility_ac(f, β_t, α, v_t, w_t) / (ω * m_eff * m_e)
                    append!(μ_t, μ_f)

                    # Evaluate optical absorptions.
                    Γ_f = optical_absorption_ac(f, β_t, α, v_t, w_t) / (100 * c * ϵ_0 * sqrt(ϵ_optic))
                    append!(Γ_t, Γ_f)

                    # Broadcast data.
                    if verbose
                        println("\e[2K", "Working on Frequency: $(f) Hz / $(efield_freq[end]) Hz")
                        println("\e[2K", "DC Mobility: ", round(μ_f, digits = 3), " cm^2/Vs")
                        println("\e[2K", "Optical absorption: ", round(Γ_f, digits = 3), " cm^-1")
                    end
                end

                if verbose
                    print("\033[F"^3) # Move cursor back
                end
            end

            # Add freq-dependent mobilities and absorptions for T > 0.
            append!(μ, [μ_t])
            append!(Γ, [Γ_t])
            if verbose
                print("\033[F"^4) # Move cursor back
            end

        else # Negative temperatures are unphysical!
            println("Temperature must be either zero or positive.")
        end
    end

    if verbose
        print("\n"^7) # Clear prints
    end

    # Return Polaron mutable struct with evaluated data.
    return Polaron(α, T, β, v, w, κ, M, F, Ω, μ, Γ)
end


text:
import FireflyModule.Firefly

function inversion_mutation(firefly::Firefly, distance_matrix, r::Int)::Firefly
    """Reverses the random part of given firefly path with
    given difference.
    """
    length_of_mutation = rand(2:r);
    max_len = length(firefly.path)-length_of_mutation
    index1 = rand(1:max_len)
    index2 = index1 + length_of_mutation
    mutated_path = reverse(firefly.path, index1, index2)
    new_f = copy(firefly)
    new_f.path = mutated_path
    new_f.cost = path_cost(new_f, distance_matrix)
    return new_f
end

max_stars_repo_name: zahidkizmaz/firefly-algorithm-tsp

max_issues_repo_name: zahidkizmaz/firefly-algorithm-tsp

max_forks_repo_name: zahidkizmaz/firefly-algorithm-tsp


text:
using Revise

using RCall

const N = Int(1e8); const K = 100;

R"""
memory.limit(2e9)
library(data.table)
N=$N; K=$K
set.seed(1)
DT <- data.table(
  id1 = sample(sprintf("id%010d",1:(N/K)), N, TRUE),      # large groups (char)
  v3 =  sample(round(runif(100,max=100),4), N, TRUE) # numeric e.g. 23.5749
)
cat("GB =", round(sum(gc()[,2])/1024, 3), "\n")
pt = proc.time()
DT[, sum(v3), keyby=id1]
x = timetaken(pt)
rm(DT); gc()
"""

@rget x;

using FastGroupBy, BenchmarkTools, DataFrames

srand(1);
df = DataFrame(id_str_large = rand(["id"*dec(k,10) for k = 1:Int(N/K)],N) , val = rand(round.(rand(K)*100,4), N));
x1 = @elapsed sumby!(df, :id_str_large, :val);

# this is way too slow
# id_str_large = rand(["id"*dec(k,10) for k = 1:Int(N/K)],N);
# val = rand(round.(rand(K)*100,4), N);
# x2 = @elapsed fastby!(sum, id_str_large, val)

x
x1
# x2
x1/parse(Float64, x[1:end-3])

srand(1);
df = DataFrame(id_str_large = rand(["id"*dec(k,10) for k = 1:Int(N/K)],N) , val = rand(round.(rand(K)*100,4), N));
using StatsBase
x3 = @elapsed countmap(df[:id_str_large], weights(df[:val]));

x3/parse(Float64, x[1:end-3])


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

const StdUvDist = Union{StandardUvUniform, StandardUvNormal}
const StdMvDist = Union{StandardMvUniform, StandardMvNormal}

struct DistributionTransform{
    VT <: AbstractValueShape,
    VF <: AbstractValueShape,
    DT <: ContinuousDistribution,
    DF <: ContinuousDistribution
} <: VariateTransform{VT,VF}
    target_dist::DT
    source_dist::DF
    target_varshape::VT
    source_varshape::VF
end

# ToDo: Add specialized dist trafo types able to cache relevant quantities, etc.

function _distrafo_ctor_impl(target_dist::Distribution, source_dist::Distribution)
    @argcheck eff_totalndof(target_dist) == eff_totalndof(source_dist)
    target_varshape = varshape(target_dist)
    source_varshape = varshape(source_dist)
    DistributionTransform(target_dist, source_dist, target_varshape, source_varshape)
end

DistributionTransform(target_dist::Distribution{VF,Continuous}, source_dist::Distribution{VF,Continuous}) where VF =
    _distrafo_ctor_impl(target_dist, source_dist)

DistributionTransform(target_dist::Distribution{Multivariate,Continuous}, source_dist::Distribution{VF,Continuous}) where VF =
    _distrafo_ctor_impl(target_dist, source_dist)

DistributionTransform(target_dist::Distribution{VF,Continuous}, source_dist::Distribution{Multivariate,Continuous}) where VF =
    _distrafo_ctor_impl(target_dist, source_dist)

DistributionTransform(target_dist::Distribution{Multivariate,Continuous}, source_dist::Distribution{Multivariate,Continuous}) =
    _distrafo_ctor_impl(target_dist, source_dist)

# apply_dist_trafo(trg_d, src_d, src_v, prev_ladj)
function apply_dist_trafo end

import Base.∘
function ∘(a::DistributionTransform, b::DistributionTransform)
    @argcheck a.source_dist == b.target_dist
    DistributionTransform(a.target_dist, b.source_dist)
end

Base.inv(trafo::DistributionTransform) = DistributionTransform(trafo.source_dist, trafo.target_dist)

ValueShapes.varshape(trafo::DistributionTransform) = trafo.source_varshape

function apply_vartrafo_impl(trafo::DistributionTransform, v::Any, prev_ladj::OptionalLADJ)
    apply_dist_trafo(trafo.target_dist, trafo.source_dist, v, prev_ladj)
end

const _StdDistType = Union{Uniform, Normal}

_trg_disttype(::Type{Uniform}, ::Type{Univariate}) = StandardUvUniform
_trg_disttype(::Type{Uniform}, ::Type{Multivariate}) = StandardMvUniform
_trg_disttype(::Type{Normal}, ::Type{Univariate}) = StandardUvNormal
_trg_disttype(::Type{Normal}, ::Type{Multivariate}) = StandardMvNormal

function _trg_dist(disttype::Type{<:_StdDistType}, source_dist::Distribution{Univariate,Continuous})
    trg_dt = _trg_disttype(disttype, Univariate)
    trg_dt()
end

function _trg_dist(disttype::Type{<:_StdDistType}, source_dist::Distribution{Multivariate,Continuous})
    trg_dt = _trg_disttype(disttype, Multivariate)
    trg_dt(eff_totalndof(source_dist))
end

function _trg_dist(disttype::Type{<:_StdDistType}, source_dist::ContinuousDistribution)
    trg_dt = _trg_disttype(disttype, Multivariate)
    trg_dt(eff_totalndof(source_dist))
end

function DistributionTransform(disttype::Type{<:_StdDistType}, source_dist::ContinuousDistribution)
    trg_d = _trg_dist(disttype, source_dist)
    DistributionTransform(trg_d, source_dist)
end

function std_dist_from(src_d::Distribution)
    throw(ArgumentError("No standard intermediate distribution defined to transform from $(typeof(src_d).name)"))
end

function std_dist_to(trg_d::Distribution)
    throw(ArgumentError("No standard intermediate distribution defined to transform into $(typeof(trg_d).name)"))
end

@inline function _intermediate_std_dist(trg_d::Distribution, src_d::Distribution)
    _select_intermediate_dist(std_dist_to(trg_d), std_dist_from(src_d))
end

@inline _intermediate_std_dist(::Union{StdUvDist,StdMvDist}, src_d::Distribution) = std_dist_from(src_d)

@inline _intermediate_std_dist(trg_d::Distribution, ::Union{StdUvDist,StdMvDist}) = std_dist_to(trg_d)

function _intermediate_std_dist(::Union{StdUvDist,StdMvDist}, ::Union{StdUvDist,StdMvDist})
    throw(ArgumentError("Direct conversions must be used between standard intermediate distributions"))
end

@inline _select_intermediate_dist(a::D, ::D) where D<:Union{StdUvDist,StdMvDist} = a
@inline _select_intermediate_dist(a::D, ::D) where D<:Union{StandardUvUniform,StandardMvUniform} = a
@inline _select_intermediate_dist(a::Union{StandardUvUniform,StandardMvUniform}, ::Union{StdUvDist,StdMvDist}) = a
@inline _select_intermediate_dist(::Union{StdUvDist,StdMvDist}, b::Union{StandardUvUniform,StandardMvUniform}) = b

_check_conv_eff_totalndof(trg_d::Uniform, src_d::Uniform) = nothing

function _check_conv_eff_totalndof(trg_d::Distribution, src_d::Distribution)
    trg_d_n = eff_totalndof(trg_d)
    src_d_n = eff_totalndof(src_d)
    if trg_d_n != src_d_n
        throw(ArgumentError("Can't convert to $(typeof(trg_d).name) with $(trg_d_n) eff. DOF from $(typeof(src_d).name) with $(src_d_n) eff. DOF"))
    end
    nothing
end

function apply_dist_trafo(trg_d::Distribution, src_d::Distribution, src_v::Any, prev_ladj::OptionalLADJ)
    _check_conv_eff_totalndof(trg_d, src_d)
    intermediate_d = _intermediate_std_dist(trg_d, src_d)
    intermediate_v, intermediate_ladj = apply_dist_trafo(intermediate_d, src_d, src_v, prev_ladj)
    apply_dist_trafo(trg_d, intermediate_d, intermediate_v, intermediate_ladj)
end

function apply_dist_trafo(trg_d::DT, src_d::DT, src_v::Real, prev_ladj::OptionalLADJ) where {DT <: StdUvDist}
    (v = src_v, ladj = prev_ladj)
end

function apply_dist_trafo(trg_d::DT, src_d::DT, src_v::AbstractVector{<:Real}, prev_ladj::OptionalLADJ) where {DT <: StdMvDist}
    @argcheck length(trg_d) == length(src_d) == length(eachindex(src_v))
    (v = src_v, ladj = prev_ladj)
end

function apply_dist_trafo(trg_d::Distribution{Univariate}, src_d::StdMvDist, src_v::AbstractVector{<:Real}, prev_ladj::OptionalLADJ)
    @argcheck length(src_d) == length(eachindex(src_v)) == 1
    apply_dist_trafo(trg_d, view(src_d, 1), first(src_v), prev_ladj)
end

function apply_dist_trafo(trg_d::StdMvDist, src_d::Distribution{Univariate}, src_v::Real, prev_ladj::OptionalLADJ)
    @argcheck length(trg_d) == 1
    r = apply_dist_trafo(view(trg_d, 1), src_d, first(src_v), prev_ladj)
    (v = unshaped(r.v), ladj = r.ladj)
end

@inline _trafo_cdf(d::Distribution{Univariate,Continuous}, x::Real) = _trafo_cdf_impl(params(d), d, x)

@inline _trafo_cdf_impl(::NTuple, d::Distribution{Univariate,Continuous}, x::Real) = cdf(d, x)

@inline function _trafo_cdf_impl(::NTuple{N,Union{Integer,AbstractFloat}}, d::Distribution{Univariate,Continuous}, x::ForwardDiff.Dual{TAG}) where {N,TAG}
    x_v = ForwardDiff.value(x)
    u = cdf(d, x_v)
    dudx = pdf(d, x_v)
    ForwardDiff.Dual{TAG}(u, dudx * ForwardDiff.partials(x))
end

@inline _trafo_quantile(d::Distribution{Univariate,Continuous}, u::Real) = _trafo_quantile_impl(params(d), d, u)

@inline _trafo_quantile_impl(::NTuple, d::Distribution{Univariate,Continuous}, u::Real) = _trafo_quantile_impl_generic(d, u)

@inline function _trafo_quantile_impl(::NTuple{N,Union{Integer,AbstractFloat}}, d::Distribution{Univariate,Continuous}, u::ForwardDiff.Dual{TAG}) where {N,TAG}
    x = _trafo_quantile_impl_generic(d, ForwardDiff.value(u))
    dxdu = inv(pdf(d, x))
    ForwardDiff.Dual{TAG}(x, dxdu * ForwardDiff.partials(u))
end

@inline _trafo_quantile_impl_generic(d::Distribution{Univariate,Continuous}, u::Real) = quantile(d, u)

# Workaround for rounding errors that can result in quantile values outside of support of Truncated:
@inline function _trafo_quantile_impl_generic(d::Truncated{<:Distribution{Univariate,Continuous}}, u::Real)
    x = quantile(d, u)
    T = typeof(x)
    min_x = T(minimum(d))
    max_x = T(maximum(d))
    if x < min_x && isapprox(x, min_x, atol = 4 * eps(T))
        min_x
    elseif x > max_x && isapprox(x, max_x, atol = 4 * eps(T))
        max_x
    else
        x
    end
end

@inline function _eval_dist_trafo_func(f::typeof(_trafo_cdf), d::Distribution{Univariate,Continuous}, src_v::Real, prev_ladj::OptionalLADJ)
    R_V = float(typeof(src_v))
    R_LADJ = !ismissing(prev_ladj) ? float(promote_type(typeof(src_v), typeof(prev_ladj))) : Missing
    if insupport(d, src_v)
        trg_v = f(d, src_v)
        trafo_ladj = !ismissing(prev_ladj) ? + logpdf(d, src_v) : missing
        var_trafo_result(convert(R_V, trg_v), src_v, convert(R_LADJ, trafo_ladj), prev_ladj)
    else
        var_trafo_result(convert(R_V, NaN), src_v, convert(R_LADJ, NaN), prev_ladj)
    end
end

@inline function _eval_dist_trafo_func(f::typeof(_trafo_quantile), d::Distribution{Univariate,Continuous}, src_v::Real, prev_ladj::OptionalLADJ)
    R_V = float(typeof(src_v))
    R_LADJ = !ismissing(prev_ladj) ? float(promote_type(typeof(src_v), typeof(prev_ladj))) : Missing
    if 0 <= src_v <= 1
        trg_v = f(d, src_v)
        trafo_ladj = !ismissing(prev_ladj) ? - logpdf(d, trg_v) : missing
        var_trafo_result(convert(R_V, trg_v), src_v, convert(R_LADJ, trafo_ladj), prev_ladj)
    else
        var_trafo_result(convert(R_V, NaN), src_v, convert(R_LADJ, NaN), prev_ladj)
    end
end

std_dist_from(src_d::Distribution{Univariate,Continuous}) = StandardUvUniform()

function apply_dist_trafo(::StandardUvUniform, src_d::Distribution{Univariate,Continuous}, src_v::Real, prev_ladj::OptionalLADJ)
    _eval_dist_trafo_func(_trafo_cdf, src_d, src_v, prev_ladj)
end

std_dist_to(trg_d::Distribution{Univariate,Continuous}) = StandardUvUniform()

function apply_dist_trafo(trg_d::Distribution{Univariate,Continuous}, ::StandardUvUniform, src_v::Real, prev_ladj::OptionalLADJ)
    TV = float(typeof(src_v))
    # Avoid src_v ≈ 0 and src_v ≈ 1 to avoid infinite variate values for target distributions with infinite support:
    mod_src_v = ifelse(src_v == 0, zero(TV) + eps(TV), ifelse(src_v == 1, one(TV) - eps(TV), convert(TV, src_v)))
    trg_v, ladj = _eval_dist_trafo_func(_trafo_quantile, trg_d, mod_src_v, prev_ladj)
    (v = trg_v, ladj = ladj)
end

function _dist_trafo_rescale_impl(trg_d, src_d, src_v::Real, prev_ladj::OptionalLADJ)
    R = float(typeof(src_v))
    trg_offs, trg_scale = location(trg_d), scale(trg_d)
    src_offs, src_scale = location(src_d), scale(src_d)
    rescale_factor = trg_scale / src_scale
    trg_v = (src_v - src_offs) * rescale_factor + trg_offs
    trafo_ladj = !ismissing(prev_ladj) ? log(rescale_factor) : missing
    var_trafo_result(trg_v, src_v, trafo_ladj, prev_ladj)
end

@inline apply_dist_trafo(trg_d::Uniform, src_d::Uniform, src_v::Real, prev_ladj::OptionalLADJ) = _dist_trafo_rescale_impl(trg_d, src_d, src_v, prev_ladj)
@inline apply_dist_trafo(trg_d::StandardUvUniform, src_d::Uniform, src_v::Real, prev_ladj::OptionalLADJ) = _dist_trafo_rescale_impl(trg_d, src_d, src_v, prev_ladj)
@inline apply_dist_trafo(trg_d::Uniform, src_d::StandardUvUniform, src_v::Real, prev_ladj::OptionalLADJ) = _dist_trafo_rescale_impl(trg_d, src_d, src_v, prev_ladj)

# ToDo: Use StandardUvNormal as standard intermediate dist for Normal? Would
# be useful if StandardUvNormal would be a better standard intermediate than
# StandardUvUniform for some other uniform distributions as well.
#
#     std_dist_from(src_d::Normal) = StandardUvNormal()
#     std_dist_to(trg_d::Normal) = StandardUvNormal()

@inline apply_dist_trafo(trg_d::Normal, src_d::Normal, src_v::Real, prev_ladj::OptionalLADJ) = _dist_trafo_rescale_impl(trg_d, src_d, src_v, prev_ladj)
@inline apply_dist_trafo(trg_d::StandardUvNormal, src_d::Normal, src_v::Real, prev_ladj::OptionalLADJ) = _dist_trafo_rescale_impl(trg_d, src_d, src_v, prev_ladj)
@inline apply_dist_trafo(trg_d::Normal, src_d::StandardUvNormal, src_v::Real, prev_ladj::OptionalLADJ) = _dist_trafo_rescale_impl(trg_d, src_d, src_v, prev_ladj)

# ToDo: Optimized implementation for Distributions.Truncated <-> StandardUvUniform

@inline function apply_dist_trafo(trg_d::StandardMvUniform, src_d::StandardMvNormal, src_v::AbstractVector{<:Real}, prev_ladj::OptionalLADJ)
    apply_dist_trafo(convert(Product, trg_d), convert(Product, src_d), src_v, prev_ladj)
end

@inline function apply_dist_trafo(trg_d::StandardMvNormal, src_d::StandardMvUniform, src_v::AbstractVector{<:Real}, prev_ladj::OptionalLADJ)
    apply_dist_trafo(convert(Product, trg_d), convert(Product, src_d), src_v, prev_ladj)
end

std_dist_from(src_d::MvNormal) = StandardMvNormal(length(src_d))

function apply_dist_trafo(trg_d::StandardMvNormal, src_d::MvNormal, src_v::AbstractVector{<:Real}, prev_ladj::OptionalLADJ)
    @argcheck length(trg_d) == length(src_d)
    A = cholesky(src_d.Σ).U
    trg_v = transpose(A) \ (src_v - src_d.μ)
    trafo_ladj = !ismissing(prev_ladj) ? - logabsdet(A)[1] : missing
    var_trafo_result(trg_v, src_v, trafo_ladj, prev_ladj)
end

std_dist_to(trg_d::MvNormal) = StandardMvNormal(length(trg_d))

function apply_dist_trafo(trg_d::MvNormal, src_d::StandardMvNormal, src_v::AbstractVector{<:Real}, prev_ladj::OptionalLADJ)
    @argcheck length(trg_d) == length(src_d)
    A = cholesky(trg_d.Σ).U
    trg_v = transpose(A) * src_v + trg_d.μ
    trafo_ladj = !ismissing(prev_ladj) ? + logabsdet(A)[1] : missing
    var_trafo_result(trg_v, src_v, trafo_ladj, prev_ladj)
end

eff_totalndof(d::Dirichlet) = length(d) - 1

std_dist_to(trg_d::Dirichlet) = StandardMvUniform(eff_totalndof(trg_d))

function apply_dist_trafo(trg_d::Dirichlet, src_d::StandardMvUniform, src_v::AbstractVector{<:Real}, prev_ladj::OptionalLADJ)
    # See https://arxiv.org/abs/1010.3436
    len_t = length(trg_d)
    @argcheck len_t == length(src_d) + 1
    in_d = product_distribution([Beta(sum(trg_d.alpha[i+1:end]),trg_d.alpha[i]) for i in 1:len_t-1])
    in_v, in_ladj = apply_dist_trafo(in_d, src_d, src_v, missing)
    trg_v = [prod(in_v[1:i-1]) * (i < len_t ? 1-in_v[i] : 1) for i in 1:len_t]
    trafo_ladj = !ismissing(prev_ladj) ? - logpdf(trg_d, trg_v) : missing
    var_trafo_result(trg_v, src_v, trafo_ladj, prev_ladj)
end

function apply_dist_trafo(trg_d::Distributions.Product, src_d::Distributions.Product, src_v::AbstractVector{<:Real}, prev_ladj::OptionalLADJ)
    rs = apply_dist_trafo.(trg_d.v, src_d.v, src_v, zero(Float32))
    trg_v = broadcast(r -> r.v, rs)
    trafo_ladj = !ismissing(prev_ladj) ? sum(map(r -> r.ladj, rs)) : missing
    var_trafo_result(trg_v, src_v, trafo_ladj, prev_ladj)
end

function apply_dist_trafo(trg_d::StdMvDist, src_d::Distributions.Product, src_v::AbstractVector{<:Real}, prev_ladj::OptionalLADJ)
    apply_dist_trafo(convert(Distributions.Product, trg_d), src_d, src_v, prev_ladj)
end

function apply_dist_trafo(trg_d::Distributions.Product, src_d::StdMvDist, src_v::AbstractVector{<:Real}, prev_ladj::OptionalLADJ)
    apply_dist_trafo(trg_d, convert(Distributions.Product, src_d), src_v, prev_ladj)
end

function _ntdistelem_to_stdmv(trg_d::StdMvDist, sd::Distribution, src_v_unshaped::AbstractVector{<:Real}, trg_acc::ValueAccessor, init_ladj::OptionalLADJ)
    td = view(trg_d, ValueShapes.view_range(Base.OneTo(length(trg_d)), trg_acc))
    sv = stripscalar(view(src_v_unshaped, trg_acc))
    apply_dist_trafo(td, sd, sv, init_ladj)
end

function _ntdistelem_to_stdmv(trg_d::StdMvDist, sd::ConstValueDist, src_v_unshaped::AbstractVector{<:Real}, trg_acc::ValueAccessor, init_ladj::OptionalLADJ)
    (v = Bool[], ladj = init_ladj)
end

_transformed_ntd_elshape(d::Distribution{Univariate}) = varshape(d)
_transformed_ntd_elshape(d::Distribution{Multivariate}) = ArrayShape{Real}(eff_totalndof(d))
function _transformed_ntd_elshape(d::Distribution)
    vs = varshape(d)
    @argcheck totalndof(vs) == eff_totalndof(d)
    vs
end

function _transformed_ntd_accessors(d::NamedTupleDist{names}) where names
    shapes = map(_transformed_ntd_elshape, values(d))
    vs = NamedTupleShape(NamedTuple{names}(shapes))
    values(vs)
end

function apply_dist_trafo(trg_d::StdMvDist, src_d::ValueShapes.UnshapedNTD, src_v::AbstractVector{<:Real}, prev_ladj::OptionalLADJ)
    src_vs = varshape(src_d.shaped)
    @argcheck length(src_d) == length(eachindex(src_v))
    trg_accessors = _transformed_ntd_accessors(src_d.shaped)
    init_ladj = ismissing(prev_ladj) ? missing : zero(Float32)
    rs = map((acc, sd) -> _ntdistelem_to_stdmv(trg_d, sd, src_v, acc, init_ladj), trg_accessors, values(src_d.shaped))
    trg_v = vcat(map(r -> r.v, rs)...)
    trafo_ladj = !ismissing(prev_ladj) ? sum(map(r -> r.ladj, rs)) : missing
    var_trafo_result(trg_v, src_v, trafo_ladj, prev_ladj)
end

function apply_dist_trafo(trg_d::StdMvDist, src_d::NamedTupleDist, src_v::Union{NamedTuple,ShapedAsNT}, prev_ladj::OptionalLADJ)
    src_v_unshaped = unshaped(src_v, varshape(src_d))
    apply_dist_trafo(trg_d, unshaped(src_d), src_v_unshaped, prev_ladj)
end

function _stdmv_to_ntdistelem(td::Distribution, src_d::StdMvDist, src_v::AbstractVector{<:Real}, src_acc::ValueAccessor, init_ladj::OptionalLADJ)
    sd = view(src_d, ValueShapes.view_range(Base.OneTo(length(src_d)), src_acc))
    sv = view(src_v, ValueShapes.view_range(axes(src_v, 1), src_acc))
    apply_dist_trafo(td, sd, sv, init_ladj)
end

function _stdmv_to_ntdistelem(td::ConstValueDist, src_d::StdMvDist, src_v::AbstractVector{<:Real}, src_acc::ValueAccessor, init_ladj::OptionalLADJ)
    (v = Bool[], ladj = init_ladj)
end

function apply_dist_trafo(trg_d::ValueShapes.UnshapedNTD, src_d::StdMvDist, src_v::AbstractVector{<:Real}, prev_ladj::OptionalLADJ)
    trg_vs = varshape(trg_d.shaped)
    @argcheck length(src_d) == length(eachindex(src_v))
    src_accessors = _transformed_ntd_accessors(trg_d.shaped)
    init_ladj = ismissing(prev_ladj) ? missing : zero(Float32)
    rs = map((acc, td) -> _stdmv_to_ntdistelem(td, src_d, src_v, acc, init_ladj), src_accessors, values(trg_d.shaped))
    trg_v_unshaped = vcat(map(r -> unshaped(r.v), rs)...)
    trafo_ladj = !ismissing(prev_ladj) ? sum(map(r -> r.ladj, rs)) : missing
    var_trafo_result(trg_v_unshaped, src_v, trafo_ladj, prev_ladj)
end

function apply_dist_trafo(trg_d::NamedTupleDist, src_d::StdMvDist, src_v::AbstractVector{<:Real}, prev_ladj::OptionalLADJ)
    unshaped_result = apply_dist_trafo(unshaped(trg_d), src_d, src_v, prev_ladj)
    (v = varshape(trg_d)(unshaped_result.v), ladj = unshaped_result.ladj)
end

const AnyReshapedDist = Union{ReshapedDist,MatrixReshaped}

function apply_dist_trafo(trg_d::Distribution{Multivariate}, src_d::AnyReshapedDist, src_v::Any, prev_ladj::OptionalLADJ)
    src_vs = varshape(src_d)
    @argcheck length(trg_d) == totalndof(src_vs)
    apply_dist_trafo(trg_d, unshaped(src_d), unshaped(src_v, src_vs), prev_ladj)
end

function apply_dist_trafo(trg_d::AnyReshapedDist, src_d::Distribution{Multivariate}, src_v::AbstractVector{<:Real}, prev_ladj::OptionalLADJ)
    trg_vs = varshape(trg_d)
    @argcheck totalndof(trg_vs) == length(src_d)
    r = apply_dist_trafo(unshaped(trg_d), src_d, src_v, prev_ladj)
    (v = trg_vs(r.v), ladj = r.ladj)
end

function apply_dist_trafo(trg_d::AnyReshapedDist, src_d::AnyReshapedDist, src_v::AbstractVector{<:Real}, prev_ladj::OptionalLADJ)
    trg_vs = varshape(trg_d)
    src_vs = varshape(src_d)
    @argcheck totalndof(trg_vs) == totalndof(src_vs)
    r = apply_dist_trafo(unshaped(trg_d), unshaped(src_d), unshaped(src_v, src_vs), prev_ladj)
    (v = trg_vs(r.v), ladj = r.ladj)
end

function apply_dist_trafo(trg_d::StdMvDist, src_d::UnshapedHDist, src_v::AbstractVector{<:Real}, prev_ladj::OptionalLADJ)
    src_v_primary, src_v_secondary = _hd_split(src_d, src_v)
    trg_d_primary = typeof(trg_d)(length(eachindex(src_v_primary)))
    trg_d_secondary = typeof(trg_d)(length(eachindex(src_v_secondary)))
    trg_v_primary, ladj_primary = apply_dist_trafo(trg_d_primary, _hd_pridist(src_d), src_v_primary, prev_ladj)
    trg_v_secondary, ladj = apply_dist_trafo(trg_d_secondary, _hd_secdist(src_d, src_v_primary), src_v_secondary, ladj_primary)
    trg_v = vcat(trg_v_primary, trg_v_secondary)
    (v = trg_v, ladj = ladj)
end

function apply_dist_trafo(trg_d::StdMvDist, src_d::HierarchicalDistribution, src_v::Any, prev_ladj::OptionalLADJ)
    src_v_unshaped = unshaped(src_v, varshape(src_d))
    apply_dist_trafo(trg_d, unshaped(src_d), src_v_unshaped, prev_ladj)
end

function apply_dist_trafo(trg_d::UnshapedHDist, src_d::StdMvDist, src_v::AbstractVector{<:Real}, prev_ladj::OptionalLADJ)
    src_v_primary, src_v_secondary = _hd_split(trg_d, src_v)
    src_d_primary = typeof(src_d)(length(eachindex(src_v_primary)))
    src_d_secondary = typeof(src_d)(length(eachindex(src_v_secondary)))
    trg_v_primary, ladj_primary = apply_dist_trafo(_hd_pridist(trg_d), src_d_primary, src_v_primary, prev_ladj)
    trg_v_secondary, ladj = apply_dist_trafo(_hd_secdist(trg_d, trg_v_primary), src_d_secondary, src_v_secondary, ladj_primary)
    trg_v = vcat(trg_v_primary, trg_v_secondary)
    (v = trg_v, ladj = ladj)
end

function apply_dist_trafo(trg_d::HierarchicalDistribution, src_d::StdMvDist, src_v::AbstractVector{<:Real}, prev_ladj::OptionalLADJ)
    unshaped_result = apply_dist_trafo(unshaped(trg_d), src_d, src_v, prev_ladj)
    (v = varshape(trg_d)(unshaped_result.v), ladj = unshaped_result.ladj)
end


text:
struct Neurons{T₁ <: Number,
               T₂ <: ActivationFunction}
    n_neurons::Int
    weights::Matrix{T₁}
    activation_function::T₂

    function Neurons(::Type{T₁},
                     n_neurons::Int,
                     n_features::Int,
                     activation_function::T₂) where {T₁ <: Number,
                                                     T₂ <: ActivationFunction}
        weights = gaussian_projection_matrix(T₁, n_neurons, n_features)
        new{T₁,T₂}(n_neurons, weights, activation_function)
    end

    function Neurons(n_neurons::Int,
                     weights::AbstractMatrix{T₁},
                     activation_function::T₂) where {T₁ <: Number,
                                                     T₂ <: ActivationFunction}
        new{T₁,T₂}(n_neurons, weights, activation_function)
    end
end

function Base.copy(ns::Neurons)
    fields = [deepcopy(getfield(ns, k)) for k ∈ fieldnames(Neurons)]
    Neurons(fields...)
end

function gaussian_projection_matrix(::Type{T},
                                    n_neurons::Int,
                                    n_features::Int) where {T <: Number}
    randn(T, n_neurons, n_features) ./ sqrt(n_features)
end


text:
# ps3000_get_values.jl

export ps3000_get_values!, ValueBuffer

immutable ValueBuffer
  a :: Array{Int16,1}
  b :: Array{Int16,1}
  c :: Array{Int16,1}
  d :: Array{Int16,1}
  ValueBuffer(length) = new(Array(Int16,length),Array(Int16,length),
                            Array(Int16,length),Array(Int16,length))
end

function ps3000_get_values!(handle       :: Int16,
                            buffer       :: ValueBuffer,
                            no_of_values)
  overflow = Ref{Int16}()  # bitfield, see doc
  values_per_channel_returned = ccall((:ps3000_get_values, ps3000driver),Int32,
    (Int16,
      Ptr{Int16}, Ptr{Int16}, Ptr{Int16}, Ptr{Int16},
      Ref{Int16}, Int32),
     handle,
     buffer.a, buffer.b, buffer.c, buffer.d,
     overflow, no_of_values)
  if values_per_channel_returned == 0
    error("ps3000 Error: Parameter Out Of Range")
  end
  return (values_per_channel_returned, overflow[])
end


text:
for (mod, f) in ((Base, :startswith), (Base, :endswith), (Compat, :contains))
    # add Symbol overloads
    @eval $mod.$f(haystack::Symbol, needle::Symbol) = $f(String(haystack), String(needle))
    @eval $mod.$f(haystack::Symbol, needle) = $f(String(haystack), needle)
    if f === :contains
        # `AbstractString` constraint needed to avoid ambiguities with Base in v1.5+.
        @eval $mod.$f(haystack::AbstractString, needle::Symbol) = $f(haystack, String(needle))
    else
        @eval $mod.$f(haystack, needle::Symbol) = $f(haystack, String(needle))
    end

    # add _any varients
    @eval begin
        """
            $($f)_any(haystack, needles)

        Given a list of `needles`, this returns `true`
        if `$($f)` would return true for any needle in that list.
        """
        function $(Symbol(f, :_any))(haystack, needles)
            return any($f(haystack, needle) for needle in needles)
        end
    end
end

# add curried varients
for f in (:startswith_any, :endswith_any, :contains_any)
    @eval begin
        """
            $($f)(needles) -> Function

        Curried form of `$($f)(haystack, needles)`
        """
        $f(needles) = Base.Fix2($f, needles)
    end
end