text:
function setindex_multiprocesses_worker(block::Array{T,N}, ba::BigArray{D,T},
                                        chunkGlobalRange::CartesianIndices{N}) where {D,T,N}
    C = get_encoding(ba)
    ba.kvStore[ cartesian_range2string(chunkGlobalRange) ] = encode( block, C)
end

"""
    put array in RAM to a BigArray
this version uses channel to control the number of asynchronized request
"""
function setindex_multiprocesses!( ba::BigArray{D,T}, buf::Array{T,N},
                       idxes::Union{UnitRange, Int, Colon} ... ) where {D,T,N}
    idxes = colon2unit_range(buf, idxes)
    # check alignment
    @assert all(map((x,y,z)->mod(first(x) - 1 - y, z), idxes, ba.offset.I, ba.chunkSize).==0) "the start of index should align with BigArray chunk size"
    t1 = time()
    baIter = ChunkIterator(idxes, ba.chunkSize; offset=ba.offset)
    @sync begin
        for (blockID, chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer) in baIter
            chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer =
                adjust_volume_boundary(ba, chunkGlobalRange, globalRange,
                                       rangeInChunk, rangeInBuffer)
            block = buf[rangeInBuffer]
            @async remotecall_fetch(setindex_multiprocesses_worker, WORKER_POOL,
                                       block, ba, chunkGlobalRange)
        end
    end
    elapsed = time() - t1 # sec
    println("saving speed: $(sizeof(buf)/1024/1024/elapsed) MB/s")
end

function getindex_multiprocesses_worker(ba::BigArray{D,T}, jobs::RemoteChannel,
                                results::RemoteChannel) where {D,T,N}
    baRange = CartesianIndices(ba)
    blockId, chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer = take!(jobs)
    if any(map((x,y)->x>y, first(globalRange).I, last(baRange).I)) || any(map((x,y)->x<y, last(globalRange).I, first(baRange).I))
        @warn("out of volume range, keep it as zeros")
        return
    end
    chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer = adjust_volume_boundary(ba, chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer)
    # finalize to avoid memory leak, see
    # https://discourse.julialang.org/t/understanding-distributed-memory-garbage-collection/8726/2
    #finalize(jobs)
    chunkSize = (last(chunkGlobalRange) - first(chunkGlobalRange) + one(CartesianIndex{N})).I
    #println("processing block in global range: $(cartesian_range2string(globalRange))")
    try
        data = ba.kvStore[ cartesian_range2string(chunkGlobalRange) ]
        chk = Codings.decode(data, get_encoding(ba))
        chk = reshape(reinterpret(T, chk), chunkSize)
        chk = chk[rangeInChunk]
        arr = OffsetArray(chk, cartesian_range2unit_range(globalRange)...)
        put!(results, arr)
    catch err
        if isa(err, KeyError) && ba.fillMissing
            println("no such key: $(err), will fill with zeros.")
        else
            println("catch an error while get index in remote worker: $err")
            @show typeof(err)
            @show stacktrace()
            rethrow()
        end
    end
end

function getindex_multiprocesses( ba::BigArray{D, T, N}, idxes::Union{UnitRange, Int}...) where {D,T,N}
    t1 = time()
    sz = map(length, idxes)
    ret = OffsetArray(zeros(T, sz), idxes...)

    channelSize = cld( nworkers(), 2 )
    jobs    = RemoteChannel(()->Channel{Tuple}( channelSize ));
    results = RemoteChannel(()->Channel{OffsetArray}( channelSize ));

    baIter = ChunkIterator(idxes, ba.chunkSize; offset=get_offset(ba))

    @sync begin
        @async begin
            for iter in baIter
                put!(jobs, iter)
            end
            close(jobs)
        end
        # control the number of concurrent requests here
        for iter in baIter
            @async remote_do(getindex_multiprocesses_worker, WORKER_POOL, ba, jobs, results)
        end

        @async begin
            for iter in baIter
                block = take!(results)
                ret[axes(block)...] = parent(block)
            end
            close(results)
        end
    end
    elapsed = time() - t1 # seconds
    println("cutout speed: $(sizeof(ret)/1024/1024/elapsed) MB/s")
    # handle single element indexing, return the single value
    ret
end


text:
# event based vs euler depending on event (i.e. spike) rate

using Pkg; Pkg.activate("./../../core/"); Pkg.instantiate()
push!(LOAD_PATH, string(pwd(),"./../../core/src/lifintegrator/core/"))
using EventBasedIntegrator, LinearAlgebra, ProgressMeter
using TimerOutputs, CPUTime, BenchmarkTools
using HDF5, PyPlot

dts = [0.01, 0.1] # ms
t_end = 1e3 # ms
in_current = 1.51 .+ collect(0:2:14) # 8

V_REST = -65. # * b2.mV
V_RESET = -65. # * b2.mV
FIRING_THRESHOLD = -50. # * b2.mV
MEMBRANE_RESISTANCE = 10. # * b2.Mohm
MEMBRANE_TIME_SCALE = 10. # * b2.ms
ABSOLUTE_REFRACTORY_PERIOD = 0. # * b2.ms

include("./eventbasedlifintegrator_defs/eventbasedlifintegrator_defs.jl")

function single_neuron(in_current; dt = 0.1, doEuler = false)
    network_julia, spike_monitor_julia, profiling_time_julia =
    simulate_LIF_neuron(10 * in_current; # should actually be voltage (R = 1 Ω)
                        simulation_time = t_end,
                        doEuler = doEuler,
                        dt = dt,
                        v_rest=V_REST,
                        v_reset=V_RESET,
                        firing_threshold=FIRING_THRESHOLD,
                        membrane_resistance=MEMBRANE_RESISTANCE,
                        membrane_time_scale=MEMBRANE_TIME_SCALE,
                        abs_refractory_period=ABSOLUTE_REFRACTORY_PERIOD)

    spike_times_julia = [spike[2] for spike in spike_monitor_julia.spikes]
    return profiling_time_julia, length(spike_times_julia)
end

## eventbased

ts_eventbased = zeros(length(in_current))
ns_eventbased = zeros(length(in_current))
for j in 1:length(in_current)
        ts_eventbased[j], ns_eventbased[j] = single_neuron(in_current[j]; doEuler = false)
end

## Euler
ts_euler = zeros(length(dts), length(in_current))
ns_euler = zeros(length(dts), length(in_current))
for i in 1:length(dts)
    for j in 1:length(in_current)
        ts_euler[i,j], ns_euler[i,j] = single_neuron(in_current[j]; dt = dts[i], doEuler = true)
    end
end

fid = h5open("./data/profiling_eventrate.h5", "w")
write(fid, "ts_eventbased", ts_eventbased)
write(fid, "ns_eventbased", ns_eventbased)
write(fid, "ts_euler", ts_euler)
write(fid, "ns_euler", ns_euler)
write(fid, "dts", dts)
close(fid)
#
# figure()
# plot(ns_eventbased, ts_eventbased, label = "eventbased")
# plot(ns_euler[1,:], ts_euler[1,:], label = string("euler, dt=",dts[1])
# plot(ns_euler[2,:], ts_euler[2,:], label = string("euler, dt=",dts[2])
# xlabel("# spikes (during 1 s)")
# ylabel("profiling time [s]")
# legend()
#
# savefig("./figures/profiling_eventrate.pdf")

ext: jl


text:
ϕ_vortex_32 =  generate_complex_potential_vortex(; pos = complex(0.0, 1.0)m, vorticity = 1.0m²/s / 2π)
ϕ_source_32 = generate_complex_potential_source(; pos = complex(3.0, 0.0)m, massflowout = 1.0m²/s)
ϕ_sink_32 = generate_complex_potential_source(; pos = -complex(3.0, 0.0)m, massflowout = -1.0m²/s)
ϕ_32(p) = ϕ_vortex_32(p) + ϕ_source_32(p) + ϕ_sink_32(p)

"""
    convolute_image_32(f_xy; physwidth = 10.0m, physheight = 4.0m)

Input:
    f_xy: (Q, Q)  → Q   Function taking two coordinates, outputs a quantity
    physwidth
    physheight
    cutoff

Output:
    Array{Complex{Float64},2}

where
    Q = Quantity

Outputs phase and amplitude information for every pixel in a line-integral-convolution rendering.

The input typically represents velocities.
"""
function convolute_image_32(f_xy; physwidth = 10.0m, physheight = 4.0m)
    xs, ys = x_y_iterators_at_pixels(;physwidth, physheight)
    # Make a continuous function, interpolated between pixels. This is supposed to
    # be faster than the original function (which it may not be)
    fxy = if f_xy isa Extrapolation
        fxy
    else
        function_to_interpolated_function(f_xy; physwidth = physwidth, physheight = physheight)
    end
    convolute_image_32(fxy, xs, ys)
end

"""
    convolute_image_32(f_xy::Extrapolation, n_xy::Extrapolation, xs, ys)

Input:
    f_xy: (Q, Q)  → (Q, Q)   Function taking two coordinates, outputs a tuple of quantities
    xs    Iterator{Q}        Arguments for f_xy, to be combined with ys
    ys    Iterator{Q}        Arguments for f_xy, to be combined with xs

Output:
    Array{Complex{Float64},2}

where
    Q = Quantity

"""
function convolute_image_32(f_xy::Extrapolation, n_xy::Extrapolation, xs, ys)
    # Image processing convention: a column has a horizontal line of pixels
    M = Array{Complex{Float64}}(undef, length(ys), length(xs))
    @assert eltype(xs) <: Quantity{Float64}
    @assert eltype(ys) <: Quantity{Float64}
    @assert n_xy(xs[1], ys[1]) isa Float64
    # Length of a streamline in time, including forward and back
    Δt = 2.0s
    # Number of waves over one streamline
    wave_per_streamline = 2
    # Frequency of interest
    freq_0 = wave_per_streamline / Δt
    # Number of sample points per set (1 to 9 from tracking backward, 10 (start) to including 20 forward)
    n = 20
    # Sampling frequency.
    freq_s = n / Δt
    rowsy = length(ys)
    for (i::Int64, x::Quantity{Float64}) in enumerate(xs), (j::Int64, y::Quantity{Float64}) in enumerate(ys)
        # Find the phase and magnitude for one pixel by projecting noise on 20 points on the streamline passing through it
        pv  = line_integral_convolution_complex(f_xy, n_xy, x, y, freq_s, freq_0)
        # Find the original indexes and update the image matrix
        M[rowsy + 1 - j, i] = pv
    end
    M
end

"""
    convolute_image_32(f_xy, xs, ys)

Input:
    f_xy: (Q, Q)  → (Q, Q)   Function taking two coordinates, outputs a tuple of quantities
    xs    Iterator{Q}        Arguments for f_xy, to be combined with ys
    ys    Iterator{Q}        Arguments for f_xy, to be combined with xs

Output:
    Array{Complex{Float64},2}

Where: Q is a quantity

Outputs phase and amplitude information for every pixel in a line-integral-convolution rendering.

The input f_xy typically represents velocities, the output all frames in a visualization.
"""
function convolute_image_32(f_xy, xs, ys)
    # Prepare a noise function:
    # Simplex based continuous noise function
    nxy = noise_for_lic(f_xy, xs, ys)
    # Convolute forward and back for every pixel
    convolute_image_32(f_xy, nxy, xs, ys)
end

"""
    convolute_image_32(matrix::Array{Quantity{Complex{Float64}, D, U}, 2}) where {D, U}

Assuming the image processing convention: a column has a horizontal line of pixels.
"""
function convolute_image_32(matrix::Array{Quantity{Complex{Float64}, D, U}, 2}) where {D, U}
    # Make a continuous function, interpolated between pixels.
    fxy = matrix_to_function(matrix)
    ny, nx = size(matrix)
    # We assume every element correspond to a position, unit of length
    physwidth = nx * scale_pt_to_unit(m)
    physheight = ny * scale_pt_to_unit(m)
    # Iterators for position - linear integration between points.
    xs = range(-physwidth / 2, stop = physwidth / 2, length = nx)
    ys = range(-physheight / 2, stop = physheight / 2, length = ny)
    convolute_image_32(fxy, xs, ys)
end

"""
    phase_histog(z, n_bins)
    -> binends, binwidth, relative_bin_counts

z is a collection of complex numbers
n_bins is the number of bins for the histogram
"""
function phase_histog(z, n_bins)
    binwidth = 2π / n_bins
    binends = range(binwidth - π, π, length = n_bins)
    binno(θ) = 1 + round(Int, (θ + π)/ (binwidth ), RoundDown)
    bincounts = [0 for bin in binends]
    for x in z
        θ = angle(x)
        bincounts[binno(θ)] += 1
    end
    binends, binwidth, normalize_datarange(bincounts / sum(bincounts))
end


text:
function check_error(err, error_msg)
    if err != 0
        error(error_msg * " with error: $(err)")
    end
end

const FREE_FONT_LIBRARY = FT_Library[C_NULL]

function ft_init()
    FREE_FONT_LIBRARY[1] != C_NULL && error("Freetype already initalized. init() called two times?")
    err = FT_Init_FreeType(FREE_FONT_LIBRARY)
    return err == 0
end

function ft_done()
    FREE_FONT_LIBRARY[1] == C_NULL && error("Library == CNULL. FreeTypeAbstraction.done() called before init(), or done called two times?")
    err = FT_Done_FreeType(FREE_FONT_LIBRARY[1])
    FREE_FONT_LIBRARY[1] = C_NULL
    return err == 0
end

function newface(facename, faceindex::Real=0, ftlib=FREE_FONT_LIBRARY)
    face = Ref{FT_Face}()
    err = FT_New_Face(ftlib[1], facename, Int32(faceindex), face)
    check_error(err, "Couldn't load font $facename")
    return face[]
end

struct FontExtent{T}
    vertical_bearing::Vec{2, T}
    horizontal_bearing::Vec{2, T}

    advance::Vec{2, T}
    scale::Vec{2, T}
end

hadvance(ext::FontExtent) = ext.advance[1]
inkwidth(ext::FontExtent) = ext.scale[1]
inkheight(ext::FontExtent) = ext.scale[2]
hbearing_ori_to_left(ext::FontExtent) = ext.horizontal_bearing[1]
hbearing_ori_to_top(ext::FontExtent) = ext.horizontal_bearing[2]
leftinkbound(ext::FontExtent) = hbearing_ori_to_left(ext)
rightinkbound(ext::FontExtent) = leftinkbound(ext) + inkwidth(ext)
bottominkbound(ext::FontExtent) = hbearing_ori_to_top(ext) - inkheight(ext)
topinkbound(ext::FontExtent) = hbearing_ori_to_top(ext)

BroadcastStyle(::Type{<: FontExtent}) = Style{FontExtent}()
BroadcastStyle(::Style{FontExtent}, x) = Style{FontExtent}()
BroadcastStyle(x, ::Style{FontExtent}) = Style{FontExtent}()

function broadcasted(op::Function, f::FontExtent, scaling::StaticVector)
    return FontExtent(
        op.(f.vertical_bearing, scaling[1]),
        op.(f.horizontal_bearing, scaling[2]),
        op.(f.advance, scaling),
        op.(f.scale, scaling),
    )
end

function broadcasted(op::Function, f::FontExtent)
    return FontExtent(
        op.(f.vertical_bearing),
        op.(f.horizontal_bearing),
        op.(f.advance),
        op.(f.scale),
    )
end

function broadcasted(op::Function, ::Type{T}, f::FontExtent) where T
    return FontExtent(
        map(x-> op(T, x), f.vertical_bearing),
        map(x-> op(T, x), f.horizontal_bearing),
        map(x-> op(T, x), f.advance),
        map(x-> op(T, x), f.scale),
    )
end

function FontExtent(fontmetric::FreeType.FT_Glyph_Metrics, scale::T = 64.0) where T <: AbstractFloat
    return FontExtent(
        Vec{2, T}(fontmetric.vertBearingX, fontmetric.vertBearingY) ./ scale,
        Vec{2, T}(fontmetric.horiBearingX, fontmetric.horiBearingY) ./ scale,
        Vec{2, T}(fontmetric.horiAdvance, fontmetric.vertAdvance) ./ scale,
        Vec{2, T}(fontmetric.width, fontmetric.height) ./ scale
    )
end

function ==(x::FontExtent, y::FontExtent)
    return (x.vertical_bearing == y.vertical_bearing &&
            x.horizontal_bearing == y.horizontal_bearing &&
            x.advance == y.advance &&
            x.scale == y.scale)
end

function FontExtent(fontmetric::FreeType.FT_Glyph_Metrics, scale::Integer)
    return FontExtent(
        div.(Vec{2, Int}(fontmetric.vertBearingX, fontmetric.vertBearingY), scale),
        div.(Vec{2, Int}(fontmetric.horiBearingX, fontmetric.horiBearingY), scale),
        div.(Vec{2, Int}(fontmetric.horiAdvance, fontmetric.vertAdvance), scale),
        div.(Vec{2, Int}(fontmetric.width, fontmetric.height), scale)
    )
end

function bearing(extent::FontExtent)
    return Vec2f0(extent.horizontal_bearing[1],
                  -(extent.scale[2] - extent.horizontal_bearing[2]))
end

function safe_free(face)
    ptr = getfield(face, :ft_ptr)
    if ptr != C_NULL && FREE_FONT_LIBRARY[1] != C_NULL
        FT_Done_Face(face)
    end
end

function boundingbox(extent::FontExtent)
    mini = bearing(extent)
    return Rect2D(mini, Vec2f0(extent.scale))
end

mutable struct FTFont
    ft_ptr::FreeType.FT_Face
    use_cache::Bool
    extent_cache::Dict{Char, FontExtent{Float32}}
    function FTFont(ft_ptr::FreeType.FT_Face, use_cache::Bool=true)
        extent_cache = Dict{Tuple{Int, Char}, FontExtent{Float32}}()
        face = new(ft_ptr, use_cache, extent_cache)
        finalizer(safe_free, face)
        return face
    end
end

use_cache(face::FTFont) = getfield(face, :use_cache)
get_cache(face::FTFont) = getfield(face, :extent_cache)

function FTFont(path::String)
    return FTFont(newface(path))
end

# C interop
function Base.cconvert(::Type{FreeType.FT_Face}, font::FTFont)
    return font
end

function Base.unsafe_convert(::Type{FreeType.FT_Face}, font::FTFont)
    return getfield(font, :ft_ptr)
end

function Base.propertynames(font::FTFont)
    return fieldnames(FreeType.FT_FaceRec)
end

function Base.getproperty(font::FTFont, fieldname::Symbol)
    fontrect = unsafe_load(getfield(font, :ft_ptr))
    field = getfield(fontrect, fieldname)
    if field isa Ptr{FT_String}
        field == C_NULL && return ""
        return unsafe_string(field)
    # Some fields segfault with unsafe_load...Lets find out which another day :D
    elseif field isa Ptr{FreeType.LibFreeType.FT_GlyphSlotRec}
        return unsafe_load(field)
    else
        return field
    end
end

# Allow broadcasting over fonts
Base.Broadcast.broadcastable(ft::FTFont) = Ref(ft)

function set_pixelsize(face::FTFont, size::Integer)
    err = FT_Set_Pixel_Sizes(face, size, size)
    check_error(err, "Couldn't set pixelsize")
    return size
end

function kerning(c1::Char, c2::Char, face::FTFont)
    i1 = FT_Get_Char_Index(face, c1)
    i2 = FT_Get_Char_Index(face, c2)
    kerning2d = Ref{FreeType.FT_Vector}()
    err = FT_Get_Kerning(face, i1, i2, FreeType.FT_KERNING_DEFAULT, kerning2d)
    # Can error if font has no kerning! Since that's somewhat expected, we just return 0
    err != 0 && return Vec2f0(0)
    # 64 since metrics are in 1/64 units (units to 26.6 fractional pixels)
    divisor = 64
    return Vec2f0(kerning2d[].x / divisor, kerning2d[].y / divisor)
end

function get_extent(face::FTFont, char::Char)
    if use_cache(face)
        get!(get_cache(face), char) do
            return internal_get_extent(face, char)
        end
    else
        return internal_get_extent(face, char)
    end
end

function internal_get_extent(face::FTFont, char::Char)
    #=
    Load chars without scaling. This leaves all glyph metrics that can be
    retrieved in font units, which can be normalized by dividing with the
    font's units_per_EM. This is more robust than relying on extents
    that are only valid with a specific pixelsize, because a font's
    pixelsize can be silently changed by third parties, such as Cairo.
    If that happens, all glyph metrics are incorrect. We avoid this by using the normalized space.
    =#
    err = FT_Load_Char(face, char, FT_LOAD_NO_SCALE)
    check_error(err, "Could not load char to get extent.")
    # This gives us the font metrics in normalized units (0, 1), with negative
    # numbers interpreted as an offset
    return FontExtent(face.glyph.metrics, Float32(face.units_per_EM))
end

descender(font) = font.descender / font.units_per_EM
ascender(font) = font.ascender / font.units_per_EM


text:
export hermite, HermiteBasis

"""
    y = hermite(n::Int, x::Real)

Evaluate the ``n``-degree Hermite polynomial at ``x``.

```math
n!\\sum_{m=0}^{[n/2]} \\frac{(-1)^m}{m!(n-2m)!} \\frac{x^{n-2m}}{2^m}
```
"""
function hermite(n::Int, x::Real)
  if n == 0
    1
  else
    N = Int(floor(n/2))
    fac = factorial
    s = reduce((s, m) -> s + (-1)^m * x^(n-2m) / (fac(m) * fac(n-2m) * 2^m), 0:N; init=0)
    fac(n) * s
  end
end

"""
    basis = HermiteBasis(vectors)

Encode an ordered basis (`basis`) of polynomials ``(ψ_j)_{j\\in \\mathcal J_i}`` of the form

```math
ψ_{(j_1,\\ldots,j_n)}(θ_1, \\ldots, θ_n) = \\prod_{k=1}^n φ_{j_k}(θ_k),
```

where ``φ_{j_k}`` is the Hermite polynomial of degree ``j_k``.
The indices ``\\mathcal J_i`` of the polynomials is provided by `vectors` in the natural
way. For example,

```julia
HermiteBasis([[1, 2, 3], [0, 1, 1, 2]])
```

will encode the ordered basis ``( ψ_{(1,2,3)}, ψ_{(0,1,1,2)} )``.
"""
struct HermiteBasis <: OrderedBasis
  vectors::AbstractArray{<:AbstractArray{<:Int}}
  function HermiteBasis(vectors::AbstractArray{<:AbstractArray{<:Int}})
    if all(vec -> all(x -> x >= 0, vec), vectors)
      new(vectors)
    else
      throw(
        ArgumentError("A `HermiteBasis` must only include vectors with positive components.")
       )
    end
  end
end

Base.length(basis::HermiteBasis) = length(basis.vectors)

function Base.union(basis::HermiteBasis, bases...)
  HermiteBasis(union(basis.vectors, map(b -> b.vectors, bases)...))
end

function (basis::HermiteBasis)(theta::AbstractArray{<:Real})
  if any(vec -> length(vec) > length(theta), basis.vectors)
    throw(ArgumentError("$theta is not in the domain of $basis"))
  else
    map(vec -> prod(map(p -> hermite(p[2], theta[p[1]]), enumerate(vec))), basis.vectors)
  end
end

max_stars_repo_name: mvarble/gen-transport-mcmc

max_issues_repo_name: mvarble/gen-transport-mcmc

max_forks_repo_name: mvarble/gen-transport-mcmc


text:
module ITensorParallel

using Distributed
using MPI
using ITensors

import ITensors: product, position!, noiseterm, lproj, rproj

include("default_in_partition.jl")
include("opsum_sum.jl")
include("threaded_projmposum.jl")
include("distributed_projmposum.jl")
include("mpi_projmposum.jl")

export ThreadedProjMPOSum, opsum_sum, DistributedSum, MPISum

end

max_stars_repo_name: mtfishman/ITensorParallel

max_issues_repo_name: mtfishman/ITensorParallel

max_forks_repo_name: mtfishman/ITensorParallel


text:
import SCS
import Convex
import InsarTimeseries

import PyPlot
plt = PyPlot

huber_fit = InsarTimeseries.huber_fit
huber = InsarTimeseries.huber

obj_l1(z) = norm(z, 1)
obj_l2(z) = norm(z, 2)
obj_huber(z) = sum(huber(z)) / 2

function l1_fit(A, b)
    x = Convex.Variable(size(A, 2))
    problem = Convex.minimize(norm(A * x - b, 1))
    Convex.solve!(problem, SCS.SCSSolver(verbose = false))
    return x.value
end

# Fit a line with L1, huber, and L2
function compare_l1_huber(m = 500, m_true = 1, b_true = -5, sigma = 3; with_outlier = true)
    # slope of line is m_true, y intercept is b_true
    x_true = collect(range(0, 10, length = m))
    data_true = (m_true .* x_true) .+ b_true
    @show m_true, b_true
    # b = x_true .+ (sigma .* randn(m))  # Gaussian noise
    b = data_true .+ (2sigma .* rand(m) .- sigma)  # Uniform noise [-sigma, sigma]
    if with_outlier
        extra = maximum(b) - minimum(b)
        b[1] += 2extra
        b[end] -= 2extra
    end

    n = 1
    A = hcat(x_true, ones(m))

    huber_sol = huber_fit(A, b)
    l1_sol = l1_fit(A, b)
    l2_sol = A \ b

    plt.figure()
    plt.plot(x_true, data_true, label = "truth")
    plt.scatter(x_true, b, label = "data")
    plt.plot(x_true, A * huber_sol, label = "Huber")
    plt.plot(x_true, A * l1_sol, label = "L1")
    plt.plot(x_true, A * l2_sol, label = "L2")
    plt.legend()
    plt.show(block = false)

    res_huber = A * huber_sol .- b
    res_l1 = A * l1_sol .- b
    res_l2 = A * l2_sol .- b
    return huber_sol, l1_sol, l2_sol, huber_sol - l1_sol, res_huber, res_l1, res_l2
end

function plot_resid(rh, r1, r2, bins = nothing)
    vm = max(maximum(abs.(rh)), maximum(abs.(r1)))
    fig, axes = plt.subplots(1, 3)
    n1, _, _ = axes[1].hist(rh, range = (-vm, vm), bins = bins)
    n2, _, _ = axes[2].hist(r1, range = (-vm, vm), bins = bins)
    n3, _, _ = axes[3].hist(r2, range = (-vm, vm), bins = bins)
    axes[1].set_title("Huber residuals")
    axes[2].set_title("L1 (CVX) residuals")
    axes[3].set_title("L2 least sq.")
    nmax = max(maximum(n1), maximum(n2), maximum(n3))
    axes[1].set_ylim(0, nmax + 2)
    axes[2].set_ylim(0, nmax + 2)
    axes[3].set_ylim(0, nmax + 2)
    return fig, axes
end

points = 40
sh, s1, s2, dif, rh, r1, r2 = compare_l1_huber(points, with_outlier = true)

nbins = 35
fig, axes = plot_resid(rh, r1, r2, nbins)


text:
@testset "README examples" begin

    # Don't forget to update README.md! TODO: Use some real tool to handle this instead...

    @testset "First example" begin
        using NaiveNASlib
        using Test
        in1 = inputvertex("in1", 1)
        in2 = inputvertex("in2", 1)

        # Create a new vertex which computes the sum of in1 and in2
        computation = in1 + in2
        @test typeof(computation) == MutationVertex

        # CompGraph helps evaluating the whole graph as a function
        graph = CompGraph([in1, in2], computation);

        # Evaluate the function represented by graph
        @test graph(2,3) == 5
    end

    @testset "Second and third examples" begin

        # First we need something to mutate. Batteries excluded, remember?
        mutable struct SimpleLayer
            W
            SimpleLayer(W) = new(W)
            SimpleLayer(nin, nout) = new(ones(Int, nin,nout))
        end
        (l::SimpleLayer)(x) = x * l.W

        # Helper function which creates a mutable layer.
        layer(in, outsize) = absorbvertex(SimpleLayer(nout(in), outsize), outsize, in, mutation=IoSize)

        invertex = inputvertex("input", 3)
        layer1 = layer(invertex, 4);
        layer2 = layer(layer1, 5);

        @test [nout(layer1)] == nin(layer2) == [4]

        # Lets change the output size of layer1:
        Δnout(layer1, -2);

        @test [nout(layer1)] == nin(layer2) == [2]

        ### Third example ###
        # First a few "normal" layers
        invertex = inputvertex("input", 6);
        start = layer(invertex, 6);
        split = layer(start, div(nout(invertex) , 3));

        # When multiplying with a scalar, the output size is the same as the input size.
        # This vertex type is said to be SizeInvariant (in lack of better words).
        scalarmult(v, f::Integer) = invariantvertex(x -> x .* f, v)

        # Concatenation is a third type of vertex, called SizeStack where the output size is the sum of input sizes
        joined = conc(scalarmult(split, 2), scalarmult(split,3), scalarmult(split,5), dims=2);
        @test trait(joined) == SizeStack()

        # Elementwise addition is SizeInvariant
        out = start + joined;
        @test trait(out) == SizeInvariant()

        @test [nout(invertex)] == nin(start) == nin(split) == [3 * nout(split)] == [sum(nin(joined))] == [nout(out)] == [6]
        @test [nout(start), nout(joined)] == nin(out) == [6, 6]

        graph = CompGraph(invertex, out)
        @test graph((ones(Int, 1,6))) == [78  78  114  114  186  186]

        # Ok, lets try to change the size of the vertex "out".
        # First we need to realize that we can only change it by integer multiples of 3
        # This is because it is connected to "split" through three paths which require nin==nout
        # Therefore, any size change to nout of "split" will result in 3 times the change of nin of "out".
        # Equivalently, nout of "split" is nin of "out" divided by 3 and nin/nout must be integers.

        # We need this information from the layer. Some layers have other requirements
        NaiveNASlib.minΔnoutfactor(::SimpleLayer) = 1
        NaiveNASlib.minΔninfactor(::SimpleLayer) = 1

        @test minΔnoutfactor(out) == minΔninfactor(out) == 3

        # Next, we need to define how to mutate our SimpleLayer
        NaiveNASlib.mutate_inputs(l::SimpleLayer, newInSize) = l.W = ones(Int, newInSize, size(l.W,2))
        NaiveNASlib.mutate_outputs(l::SimpleLayer, newOutSize) = l.W = ones(Int, size(l.W,1), newOutSize)

        # In some cases it is useful to hold on to the old graph before mutating
        # To do so, we need to define the clone operation for our SimpleLayer
        NaiveNASlib.clone(l::SimpleLayer) = SimpleLayer(l.W)
        parentgraph = copy(graph)

        Δnout(out, 3)

        # We didn't touch the input when mutating...
        @test [nout(invertex)] == nin(start) == [6]
        # Start and joined must have the same size due to elementwise op.
        # All three scalarmult vertices are transparent and propagate the size change to split
        @test [nout(start)] == nin(split) == [3 * nout(split)] == [sum(nin(joined))] == [nout(out)] == [9]
        @test [nout(start), nout(joined)] == nin(out) == [9, 9]

        # However, this only updated the mutation metadata, not the actual layer.
        # Some reasons for this are shown in the pruning example below
        @test graph((ones(Int, 1,6))) == [78  78  114  114  186  186]

        # To mutate the graph, we need to apply the mutation:
        apply_mutation(graph);

        @test graph((ones(Int, 1,6))) == [114  114  114  168  168  168  276  276  276]

        # Copy is still intact
        @test parentgraph((ones(Int, 1,6))) == [78  78  114  114  186  186]

        @testset "Add layers example" begin

            invertex = inputvertex("input", 3)
            layer1 = layer(invertex, 5)
            graph = CompGraph(invertex, layer1)

            @test nv(graph) == 2
            @test graph(ones(Int, 1, 3)) == [3 3 3 3 3]

            # Insert a layer between invertex and layer1
            insert!(invertex, vertex -> layer(vertex, nout(vertex)))

            @test nv(graph) == 3
            @test graph(ones(Int, 1, 3)) == [9 9 9 9 9]
        end

        @testset "Remove layers example" begin
            invertex = inputvertex("input", 3)
            layer1 = layer(invertex, 5)
            layer2 = layer(layer1, 4)
            graph = CompGraph(invertex, layer2)

            @test nv(graph) == 3
            @test graph(ones(Int, 1, 3)) == [15 15 15 15]

            # Remove layer1 and change nin of layer2 from 5 to 3
            # Would perhaps have been better to increase nout of invertex, but it is immutable
            remove!(layer1)
            apply_mutation(graph)

            @test nv(graph) == 2
            @test graph(ones(Int, 1, 3)) == [3 3 3 3]
        end

        @testset "Add edge example" begin
            invertices = inputvertex.(["input1", "input2"], [3, 2])
            layer1 = layer(invertices[1], 4)
            layer2 = layer(invertices[2], 4)
            add = layer1 + layer2
            out = layer(add, 5)
            graph = CompGraph(invertices, out)

            @test nin(add) == [4, 4]
            # Two inputs this time, remember?
            @test graph(ones(Int, 1, 3), ones(Int, 1, 2)) == [20 20 20 20 20]

            # This graph is not interesting enough for there to be a good showcase for adding a new edge.
            # Lets create a new layer which has a different output size just to see how things change
            # The only vertex which support more than one input is add
            layer3 = layer(invertices[2], 6)
            create_edge!(layer3, add)
            apply_mutation(graph)

            # By default, NaiveNASlib will try to increase the size in case of a mismatch
            @test nin(add) == [6, 6, 6]
            @test graph(ones(Int, 1, 3), ones(Int, 1, 2)) == [42 42 42 42 42]
        end

        @testset "Remove edge example" begin
            invertex = inputvertex("input", 4)
            layer1 = layer(invertex, 3)
            layer2 = layer(invertex, 5)
            merged = conc(layer1, layer2, layer1, dims=2)
            out = layer(merged, 3)
            graph = CompGraph(invertex, out)

            @test nin(merged) == [3, 5, 3]
            @test graph(ones(Int, 1, 4)) == [44 44 44]

            remove_edge!(layer1, merged)
            apply_mutation(graph)

            @test nin(merged) == [5, 3]
            @test graph(ones(Int, 1, 4)) == [32 32 32]
        end

        @testset "Pruning example" begin
            # Some mockup 'batteries' for this example

            # First, how to select or add rows or columns to a matrix
            # Negative values in selected indicate rows/cols insertion at that index
            function select_params(W, selected, dim)
                Wsize = collect(size(W))
                indskeep = repeat(Any[Colon()], 2)
                newmap = repeat(Any[Colon()], 2)

                # The selected indices
                indskeep[dim] = filter(ind -> ind > 0, selected)
                # Where they are 'placed', others will be zero
                newmap[dim] = selected .> 0
                Wsize[dim] = length(newmap[dim])

                newmat = zeros(Int64, Wsize...)
                newmat[newmap...] = W[indskeep...]
                return newmat
            end

            NaiveNASlib.mutate_inputs(l::SimpleLayer, selected::Vector{<:Integer}) = l.W = select_params(l.W, selected, 1)
            NaiveNASlib.mutate_outputs(l::SimpleLayer, selected::Vector{<:Integer}) = l.W = select_params(l.W, selected, 2)

            # Return layer just so we can easiliy look at it
            function prunablelayer(in, outsize)
                l = SimpleLayer(reshape(1: nout(in) * outsize, nout(in), :))
                return absorbvertex(l, outsize, in), l
            end

            # Ok, now lets get down to business!
            invertices = inputvertex.(["in1", "in2"], [3,4])
            v1, l1 = prunablelayer(invertices[1], 4)
            v2, l2 = prunablelayer(invertices[2], 3)
            merged = conc(v1, v2, dims=2)
            v3, l3 = prunablelayer(merged, 2)
            graph = CompGraph(invertices, v3)

            # These weights are of course complete nonsense from a neural network perspective.
            # They are just to make it easier to spot what has changed after pruning below.
            @test l1.W ==
            [ 1  4  7  10 ;
            2  5  8  11 ;
            3  6  9  12 ]

            @test l2.W ==
            [ 1  5   9 ;
            2  6  10 ;
            3  7  11 ;
            4  8  12 ]

            @test l3.W ==
            [ 1   8 ;
            2   9 ;
            3  10 ;
            4  11 ;
            5  12 ;
            6  13 ;
            7  14 ]

            # A limitation in current implementation is that one must change the size before pruning
            # See https://github.com/DrChainsaw/NaiveNASlib.jl/issues/40
            Δnin(v3, -3)

            # What did that do?
            @test nout(v1) == 2
            @test nout(v2) == 2

            # Doing this however makes it possible to do several mutations without throwing away
            # more information than needed.
            # For example, if we had first applied the previous mutation we would have thrown away
            # weights for v2 which would then just be replaced by 0s when doing this:
            Δnout(v2, 2)

            # What did that do?
            @test nout(v1) == 2
            @test nout(v2) == 4
            @test nin(v3) == [6]
            # Net result is that v1 shall decrease output size by 1 and v2 shall increase its output size by 1

            # Now, we need a utility/value metric per neuron in order to determine which neurons to keep
            # Give high utility to neurons 1 and 3 of v1, same for all others...
            utility(v) = v == v1 ? [10, 1, 10, 1] : ones(nout_org(v))
            # Then select the neurons.
            Δoutputs(graph, utility)
            # And apply it to the actual weights
            apply_mutation(graph)

            @test l1.W ==
            [ 1  7 ;
            2  8 ;
            3  9 ]

            # Note how column 3 was not replaced by zeros: We increased the target size before pruning
            @test l2.W ==
            [ 1  5   9  0;
            2  6  10  0;
            3  7  11  0;
            4  8  12  0]

            @test l3.W ==
            [ 1   8 ;
            3  10 ;
            5  12 ;
            6  13 ;
            7  14 ;
            0   0]
        end

        @testset "Advanced usage" begin

            @testset "Strategies" begin
                # A simple graph where one vertex has a constraint for changing the size.
                invertex = inputvertex("in", 3)
                layer1 = layer(invertex, 4)
                joined = conc(scalarmult(layer1, 2), scalarmult(layer1, 3), dims=2)

                # joined can only change in steps of 2
                @test minΔnoutfactor(joined) == 2

                # all_in_graph finds all vertices in the same graph as the given vertex
                verts = all_in_graph(joined)

                # Strategy to try to change it by one and throw an error when not successful
                exactOrFail = ΔNout{Exact}(joined, 1, ΔSizeFailError("Size change failed!!"))

                # Note that we now call Δsize instead of Δnout as the direction is given by the strategy
                @test_throws ErrorException Δsize(exactOrFail, verts)

                # No change was made
                @test nout(joined) == 2*nout(layer1) == 8

                # Try to change by one and fail silently when not successful
                exactOrNoop = ΔNout{Exact}(joined, 1, ΔSizeFailNoOp())

                Δsize(exactOrNoop, verts)

                # No change was made
                @test nout(joined) == 2*nout(layer1) == 8

                # In many cases it is ok to not get the exact change which was requested
                relaxedOrFail = ΔNout{Relaxed}(joined, 1, ΔSizeFailError("This should not happen!!"))

                Δsize(relaxedOrFail, verts)

                # Changed by two as this was the smallest possible change
                @test nout(joined) == 2*nout(layer1) == 10

                # Logging when fallback is applied is also possible
                using Logging
                # Yeah, this is not easy on the eyes, but it gets the job done...
                exactOrLogThenRelax = ΔNout{Exact}(joined, 1, LogΔSizeExec(Logging.Info, "Exact failed, relaxing", ΔNout{Relaxed}(joined, 1, ΔSizeFailError("This should not happen!!"))))

                @test_logs (:info, "Exact failed, relaxing") Δsize(exactOrLogThenRelax, verts)

                @test nout(joined) == 2*nout(layer1) == 12
            end

            @testset "Traits" begin
                noname = layer(inputvertex("in", 2), 2)
                @test name(noname) == "MutationVertex::SizeAbsorb"

                # Naming vertices is so useful for logging and debugging I almost made it mandatory
                named = absorbvertex(SimpleLayer(2, 3), 3, inputvertex("in", 2), traitdecoration = t -> NamedTrait(t, "named layer"))
                @test name(named) == "named layer"

                # Speaking of logging...
                layer1 = absorbvertex(SimpleLayer(2, 3), 3, inputvertex("in", 2), traitdecoration = t -> SizeChangeLogger(NamedTrait(t, "layer1")))

                # What info is shown can be controlled by supplying an extra argument to SizeChangeLogger
                nameonly = NameInfoStr()
                layer2 = absorbvertex(SimpleLayer(nout(layer1), 4), 4, layer1, traitdecoration = t -> SizeChangeLogger(nameonly, NamedTrait(t, "layer2")))

                @test_logs(
                (:info, "Change nout of layer1, inputs=[in], outputs=[layer2], nin=[2], nout=[3], SizeAbsorb() by 1"),
                (:info, "Change nin of layer2 by 1"), # Note: less verbose compared to layer1 due to NameInfoStr
                 Δnout(layer1, 1))

                 # traitdecoration works exactly the same for conc and invariantvertex as well, no need for an example

                 # Use the >> operator when creating SizeInvariant vertices using arithmetic operators:
                 add = "addvertex" >> inputvertex("in1", 1) + inputvertex("in2", 1)
                 @test name(add) == "addvertex"

                 # For more elaborate traits one can use traitconf
                 add2 = traitconf(t -> SizeChangeLogger(NamedTrait(t, "layer1 + layer2"))) >> layer1 + layer2
                 @test name(add2) == "layer1 + layer2"

                 @test_logs(
                 (:info, "Change nout of layer1, inputs=[in], outputs=[layer2, layer1 + layer2], nin=[2], nout=[4], SizeAbsorb() by 1"),
                 (:info, "Change nin of layer2 by 1"),
                 (:info, "Change nout of layer2 by 1"),
                 (:info, "Change nin of layer1 + layer2, inputs=[layer1, layer2], outputs=[], nin=[4, 4], nout=[4], SizeInvariant() by 1, 1"),
                 (:info, "Change nout of layer1 + layer2, inputs=[layer1, layer2], outputs=[], nin=[5, 5], nout=[4], SizeInvariant() by 1"),
                  Δnout(layer1, 1))

                  # When creating own trait wrappers, remember to subtype DecoratingTrait or else there will be pain!

                  # Wrong!! Not a subtype of DecoratingTrait
                  struct PainfulTrait{T<:MutationTrait} <: MutationTrait
                      base::T
                  end
                  painlayer = absorbvertex(SimpleLayer(2, 3), 3, inputvertex("in", 2), traitdecoration = PainfulTrait)

                  # Now one must implement a lot of methods for PainfulTrait...
                  @test_throws MethodError Δnout(painlayer, 1)

                  # Right! Is a subtype of DecoratingTrait
                  struct SmoothSailingTrait{T<:MutationTrait} <: DecoratingTrait
                      base::T
                  end
                  # Just implement base and all will be fine
                  NaiveNASlib.base(t::SmoothSailingTrait) = t.base

                  smoothlayer = absorbvertex(SimpleLayer(2, 3), 3, inputvertex("in", 2), traitdecoration = SmoothSailingTrait)

                  Δnout(smoothlayer, 1)
                  @test nout(smoothlayer) == 4
            end

            @testset "Graph instrumentation and modification" begin
                invertex = inputvertex("in", 2)
                layer1 = layer(invertex, 3)
                layer2 = layer(layer1, 4)

                graph = CompGraph(invertex, layer2)

                @test name.(vertices(graph)) == ["in", "MutationVertex::SizeAbsorb", "MutationVertex::SizeAbsorb"]

                # Ok, lets add names to layer1 and layer2 and change the name of invertex

                # Lets first define the default: Fallback to "clone"
                # clone is the existing function to copy things in this manner as I did not want to override Base.copy
                copyfun(args...;cf) = clone(args...;cf=cf) # Keyword argument cf is the function to use for copying all fields of the input

                # Add a name to layer1 and layer2
                function copyfun(v::MutationVertex,args...;cf)
                    # This is probably not practical to do in a real graph, so make sure you have names when first creating it...
                    name = v == layer1 ? "layer1" : "layer2"
                    addname(args...;cf) = clone(args...;cf=cf)
                    addname(t::SizeAbsorb;cf) = NamedTrait(t, name) # SizeAbsorb has no fields, otherwise we would have had to run cf for each one of them...
                    clone(v, args...;cf=addname)
                end

                # Change name of invertex
                # Here we can assume that invertex name is unique in the whole graph or else we would have had to use the above way
                copyfun(s::String; cf) = s == name(invertex) ? "in changed" : s

                # Now supply copyfun when copying the graph.
                # I must admit that thinking about what this does makes me a bit dizzy...
                namedgraph = copy(graph, copyfun)

                @test name.(vertices(namedgraph)) == ["in changed", "layer1", "layer2"]
            end

        end
    end

end


text:
#!/usr/bin/env julia

include("logger.jl")

"""
Do k-means clustering for inputs, each of which is a single data sample.
    @k the number of clusters
    @inputs the input data vector, each element is a vector in feature space.

    return the centroids, each of which is a cluster center vector in feature space.
        the cluster indicator array, each element is the cluster id which
        contains the data example indicated by the array index.
"""
function kmeans{T <: Number}(k::Int, inputs::Vector{Vector{T}}; lowerbound=0.00001, maxiter=1000)
    assert(k > 1)

    # initialize k centroids from inputs
    centroids = init_centroids(k, inputs)
    #dbglog("init: centroids=$(formatnum(centroids))")

    iter, sumdiff = 0, 1
    while sumdiff > lowerbound && iter < maxiter
        nearest_centroids = find_nearest_centroids(inputs, centroids)
        new_centroids = find_new_centroids(k, inputs, nearest_centroids)

        sumdiff = sum([dot(v, v) for v in new_centroids - centroids])
        iter += 1

        centroids = new_centroids

        #dbglog("======\n$iter: diff=$sumdiff centroids=$(formatnum(centroids))")
    end

    # the final cluster assignment
    nearest_centroids = find_nearest_centroids(inputs, centroids)

    return centroids, nearest_centroids
end

"""
Do k-means clustering for inputs
    @k the number of clusters
    @inputs the input data matrix, each column is a vector in feature space.

    return the centroids, each of which is a cluster center vector in feature space.
        the cluster indicator array, each element is the cluster id which
        contains the data example indicated by the array index.
"""
function kmeans{T <: Number}(k::Int, inputs::Matrix{T}; lowerbound=0.00001, maxiter=1000)
    inputs = Vector{T}[inputs[:, j] for j = 1:(size(inputs)[2])]
    kmeans(k, inputs, lowerbound=lowerbound, maxiter=maxiter)
end

"""
Find the new centroids
"""
function find_new_centroids{T <: Number}(k::Int, inputs::Vector{Vector{T}}, nearest_centroids)
    new_centroids = Vector{Float64}[ begin
        vsum, count = 0, 0
        for (i, min_c) in filter(x -> x[2] == c, nearest_centroids)
            vsum += inputs[i]
            count += 1
        end
        float(vsum) / count
    end for c = 1:k ]
end

"""
Find the nearest centroids for all the data samples.
    @k the number of clusters
    @inputs
    @centroids vector of length k

    return an array, each element of which is a pair of array index and id of the
    nearest cluster.
"""
function find_nearest_centroids{T <: Number}(inputs::Vector{Vector{T}}, centroids::Vector{Vector{T}})
    # find the nearest centroid for each input
    k = length(centroids)
    belonging = map(enumerate(inputs)) do x
        i, v = x
        min_c = indmin([(d = v - centroids[c]; dot(d, d)) for c in 1:k])
        i, min_c
    end
end

"""
Initialize the centroids.
    @k the number of cluster needed
    @inputs

return an array of centroids, each is a vector in feature space.
"""
function init_centroids{T <: Number}(k::Int, inputs::Vector{Vector{T}})
    assert(length(inputs) > k)

    centroids = Vector{Vector{T}}()

    for c = 1:k
        r = rand(inputs)
        while r in centroids
            r = rand(inputs)
        end
        push!(centroids, r)
    end

    return centroids
end

max_stars_repo_name: zxteloiv/julia-ex

max_issues_repo_name: zxteloiv/julia-ex

max_forks_repo_name: zxteloiv/julia-ex


text:
using BinDeps

@BinDeps.setup

somoclu_version="1.7.5"

link = "https://github.com/peterwittek/somoclu/releases/download/$(somoclu_version)/somoclu-$(somoclu_version).tar.gz"

if !is_windows()
	libsomoclu = library_dependency("libsomoclu",
								 aliases=["libsomoclu", "libsomoclu.so"], os=:Unix)
	provides(Sources, Dict([URI(link) => libsomoclu]))
	provides(BuildProcess,
			 Autotools(libtarget=joinpath("src", "libsomoclu.so"),
			           configure_options=[AbstractString("--without-mpi")]), libsomoclu, os=:Unix)
	@BinDeps.install Dict([:libsomoclu => :libsomoclu])
else
    using WinRPM
	using Base.Filesystem
	using Base.Libdl

	WinRPM.install("gcc", yes=true)
	WinRPM.install("binutils", yes=true)
	WinRPM.install("gcc-c++", yes=true)
	WinRPM.install("libstdc++6", yes=true)
	WinRPM.install("libwinpthread1", yes=true)
	WinRPM.install("win_iconv", yes=true)
	WinRPM.install("zlib1", yes=true)

	GCCROOT = joinpath(Pkg.dir("WinRPM"),"deps")

	libsomoclu = library_dependency("libsomoclu",
								 aliases=["libsomoclu", "libsomoclu.dll"],
								 os=:Windows)
	provides(Sources, Dict([URI(link) => libsomoclu]))

	makedeplnk = "https://sourceforge.net/projects/gnuwin32/files/make/3.81/make-3.81-dep.zip/download"
	makebinlnk = "https://sourceforge.net/projects/gnuwin32/files/make/3.81/make-3.81-bin.zip/download"

	makedepzip 	= joinpath(BinDeps.downloadsdir(libsomoclu), "make-3.81-dep.zip")
	makebinzip   = joinpath(BinDeps.downloadsdir(libsomoclu), "make-3.81-bin.zip")
	makebuilddir = joinpath(BinDeps.builddir(libsomoclu),  "make")

	make_s, make_d =
		joinpath(makebuilddir, "bin", "make.exe"),
		joinpath(BinDeps.bindir(libsomoclu), "make.exe")
	icnv_s, icnv_d =
		joinpath(makebuilddir, "bin", "libiconv2.dll"),
		joinpath(BinDeps.bindir(libsomoclu), "libiconv2.dll")
	intl_s, intl_d =
		joinpath(makebuilddir, "bin", "libintl3.dll"),
		joinpath(BinDeps.bindir(libsomoclu), "libintl3.dll")

	somoclusrcdir   = joinpath(BinDeps.srcdir(libsomoclu),   "somoclu-$somoclu_version")
	somoclubuilddir = joinpath(BinDeps.builddir(libsomoclu), "somoclu-$somoclu_version")
	somoclumakefile = joinpath(somoclusrcdir, "src", "Makefile.libsomoclu.mingw")

	provides(SimpleBuild,
		(@build_steps begin
			CreateDirectory(BinDeps.bindir(libsomoclu))
			FileRule(joinpath(BinDeps.bindir(libsomoclu), "make.exe"), @build_steps begin
				FileDownloader(makedeplnk, makedepzip)
				FileDownloader(makebinlnk, makebinzip)
				FileUnpacker(makedepzip, makebuilddir, joinpath("bin", "libiconv2.dll"))
				FileUnpacker(makebinzip, makebuilddir, joinpath("bin", "make.exe"))
				()-> begin
					cp(make_s, make_d, remove_destination=true)
					cp(icnv_s, icnv_d, remove_destination=true)
					cp(intl_s, intl_d, remove_destination=true)
				end
			end)
			GetSources(libsomoclu)
			CreateDirectory(somoclubuilddir)
			@build_steps begin
				ChangeDirectory(somoclubuilddir)
				FileRule("Makefile", ()-> begin
					println("Copying... Makefile")
					cp(somoclumakefile, joinpath(somoclubuilddir, "Makefile"),
					   remove_destination=true)
					end)
				`../../usr/bin/make.exe ARCH=$(Sys.ARCH) version=$(somoclu_version)`
				CreateDirectory(BinDeps.libdir(libsomoclu))
				FileRule(joinpath(BinDeps.usrdir(libsomoclu),"lib","libsomoclu.dll"), ()->
					cp("libsomoclu.dll", joinpath(BinDeps.libdir(libsomoclu), "libsomoclu.dll"),
					   remove_destination=true))
			end
		end), libsomoclu, os=:Windows)

    push!(BinDeps.defaults, SimpleBuild)
	@BinDeps.install Dict([:libsomoclu => :libsomoclu])
    pop!(BinDeps.defaults)
end


text:
module Interop

using ..LLVM
import ..LLVM: API, ref

const jlctx = Ref{LLVM.Context}()

include("interop/base.jl")
include("interop/asmcall.jl")
include("interop/passes.jl")

end


text:
function smoothing_within_interval(X::Vector{Float64}, Y::Function, interval::Tuple{Float64, Float64})
  @assert !isempty(X)
  @assert !isempty(interval)
  @assert length(interval) == 2
  @assert interval[1] ≠ interval[2]

  σˣ = sign(X[end] - X[1])
  @assert issorted(X; rev = σˣ < 0)

  σⁱ = sign(interval[2] - interval[1])

  _X = if σˣ > 0 X else reverse(X) end
  _Y = if σˣ > 0 Y else x->Y(x - X[end]) end
  _interval = if σⁱ > 0 interval else reverse(interval) end
end


text:
using Test
using SparseIR

@testset "poly.jl" begin
    @testset "typestability" begin
        @test typestable(
            SparseIR._evaluate, [SparseIR.PiecewiseLegendrePoly{Float64}, Float64];
            checkonlyany=true,
        )
        @test typestable(
            SparseIR._evaluate,
            [SparseIR.PiecewiseLegendrePoly{Float64}, Vector{Float64}];
            checkonlyany=true,
        )
    end

    @testset "shape" begin
        u, s, v = sve_logistic[42]
        l = length(s)
        @test size(u) == (l,)

        @test size(u[4]) == ()
        @test size(u[3:5]) == (3,)
    end

    @testset "knots" begin
        u, s, v = sve_logistic[42]
        @test first(u[1].knots) == -1.0
        @test last(u[1].knots) == 1.0
    end

    @testset "slice" begin
        sve_result = sve_logistic[42]

        basis = DimensionlessBasis(fermion, 42.0; sve_result)
        @test length(basis[1:5]) == 5

        basis = FiniteTempBasis(fermion, 4.2, 10.0; sve_result)
        @test length(basis[1:4]) == 4
    end

    @testset "eval" begin
        u, s, v = sve_logistic[42]
        l = length(s)

        # Evaluate
        @test u(0.4) == [u[i](0.4) for i in 1:l]
        @test u.([0.4, -0.2]) == [[u[i](x) for i in 1:l] for x in (0.4, -0.2)]
    end

    @testset "matrix_hat" begin
        u, s, v = sve_logistic[42]
        uhat = SparseIR.hat.(u, :odd)

        n = [1, 3, 5, -1, -3, 5]
        result1 = uhat[1](n)
        result = uhat(reshape(n, (3, 2)))
        result_iter = reshape(uhat(n), (:, 3, 2))
        @test size(result1) == (length(n),)
        @test size(result) == size(result_iter)
        @test result == result_iter
    end

    @testset "overlap with Λ = $Λ" for Λ in (10, 42, 10_000)
        atol = 1e-13
        u, s, v = sve_logistic[Λ]

        # Keep only even number of polynomials
        u, s, v = u[1:(end - end % 2)],
        s[1:(end - end % 2)],
        v[1:(end - end % 2)]

        @test overlap(u[1], u[1]) ≈ 1 rtol = 0 atol = atol
        @test overlap(u[1], u[2]) ≈ 0 rtol = 0 atol = atol

        ref = float.(eachindex(s) .== 1)
        @test all(isapprox.(overlap(u[1], u), ref; rtol=0, atol))
    end

    @testset "eval unique" begin
        u, s, v = sve_logistic[42]
        û = SparseIR.hat.(u, :odd)

        # evaluate
        res1 = û([1, 3, 3, 1])
        idx = [1, 2, 2, 1]
        res2 = û([1, 3])[:, idx]
        @test res1 == res2
    end
end


text:
# Programs to test writing to destinations.

# Write to hardware pointers.
@asap4asm function pointer_write_1()
    MOVI(set_pointer[0], 1)
    MOVI(set_pointer[1], 2)
    MOVI(set_pointer[2], 3)
    MOVI(set_pointer[3], 4)
end

# Write to address generators.
#
# Start is the lower 8 bits and stop is the upper 16 bits.
@asap4asm function address_generator_write_1()
    # Set start to 10, stop to 20 and stride to 3
    MOVI(ag_start[0], ((20 << 8) | 10))
    MOVI(ag_stride[0], 3)

    # Set start to 0, stop to 255, and stride to "-1"
    MOVI(ag_start[1], 255 << 8)
    MOVI(ag_stride[1], 0xFF)

    # Set start to 100, stride to 1
    MOVI(ag_start[2], 100)
    MOVI(ag_stride[2], 1)
end

# Write to CX Mask
@asap4asm function cxmask_write_1()
    MOVI(cxmask[0], 0xFFFF)
end

# Write to OBUF mask - just write all 1's.
@asap4asm function obufmask_write_1()
    MOVI(obuf_mask, 0xFF)
end


text:
function fifth_power_digit_sum(n::Int)

	a=Int[];

	while(n>=1)
		push!(a,n%10);
		n=n-n%10;
		n/=10;
		n=Int(n)
	end
	ans=0;
	for i in a
		ans+=i^5
	end
	return ans
end

function compute()
	ans=0
	for i=2:1000000
		if i==fifth_power_digit_sum(i)
			ans+=i
		end
	end
	return ans
end

println(compute())


text:
"""
# Search for the symbol of an economic data series that matches the search text.

    search_symbol(search_text, must_contain)

## Arguments

    - `search_title` : A string containing the words with which you expect to find the id of the variable.
    - `must_contain` : A string like "Germany" to explicitly search for a variable which contains "Germany" in its title

## Description

The resulting DataFrame is sorted by popularity.

## Example

```jldoctest
julia> search_symbol("Exports", "Mexico")
julia> search_symbol("GDP", "Germany")
julia> search_symbol("G", "France")
julia> x = search_symbol("Exports", "Mexico")
julia> println(x)
```
"""
function search_symbol(search_title::String, must_contain::String = "")

    search_text = replace(search_title, " " => "+")
    must_contain = replace(must_contain, " " => "+")
    url = "https://api.stlouisfed.org/fred/series/search"

    parameters = Dict("api_key" => ENV["API_FRED"],
                      "file_type" => "json",
                      "search_text" => search_text)

    response = HTTP.request("GET", url; query = parameters)
    results = JSON.parse(String(response.body))["seriess"]

    x = DataFrame(
        popularity = zeros(Int, eachindex(results)),
        title = ["" for i ∈ eachindex(results)],
        id = ["" for i ∈ eachindex(results)]
    )

    for i ∈ eachindex(results)
        x[i, 1] = results[i]["popularity"]
        x[i, 2] = results[i]["id"]
        x[i, 3] = results[i]["title"]
    end

    i = findall([occursin(must_contain, x[i, 3]) for i ∈ 1:nrow(x)] .== true)
    x = sort!(x[i, :], rev = true)

    return(x)
end


text:
function DAVIS_LSTM(w,h,c,x)
    x = relu(conv4(w[1],x;padding=1).+w[2])
    x = relu(conv4(w[3],x;padding=1).+w[4])
    x = pool(x)
    x = relu(conv4(w[5],x;padding=1).+w[6])
    x = relu(conv4(w[7],x;padding=1).+w[8])
    x = pool(x)
    x = relu(conv4(w[9],x;padding=1).+w[10])
    x = relu(conv4(w[11],x;padding=1).+w[12])
    x = relu(conv4(w[13],x;padding=1).+w[14])
    x = relu(conv4(w[15],x;padding=1).+w[16])
    x = pool(x)
    x = relu(conv4(w[17],x;padding=1).+w[18])
    x = relu(conv4(w[19],x;padding=1).+w[20])
    x = relu(conv4(w[21],x;padding=1).+w[22])
    x = relu(conv4(w[23],x;padding=1).+w[24])
    #28,28,512,1
    i = sigm(conv4(w[25],x;padding=1) .+ conv4(w[26],h;padding=1) .+ w[27])
    f = sigm(conv4(w[28],x;padding=1) .+ conv4(w[29],h;padding=1) .+ w[30])
    o = sigm(conv4(w[31],x;padding=1) .+ conv4(w[32],h;padding=1) .+ w[33])
    j = tanh(conv4(w[34],x;padding=1) .+ conv4(w[35],h;padding=1) .+ w[36])
    c = c.*f + i.*j
    h = tanh(c).*o

    x = relu(conv4(w[37],h).+w[38])
    x = deconv4(w[39],x;padding=1,stride=6).+w[40]
end

function DAVIS_LSTM_w()
    w = Array[
        #Conv64
        xavier(Float32,3,3,3,64),        zeros(Float32,1,1,64,1),
        xavier(Float32,3,3,64,64),       zeros(Float32,1,1,64,1),
        #Conv128
        xavier(Float32,3,3,64,128),      zeros(Float32,1,1,128,1),
        xavier(Float32,3,3,128,128),     zeros(Float32,1,1,128,1),
        #Conv256
        xavier(Float32,3,3,128,256),     zeros(Float32,1,1,256,1),
        xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
        xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
        xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
        #Conv512
        xavier(Float32,3,3,256,512),     zeros(Float32,1,1,512,1),
        xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
        xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
        xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
        #LSTM
        xavier(Float32,3,3,512,512),     xavier(Float32,3,3,512,512),   zeros(Float32,1,1,512,1),
        xavier(Float32,3,3,512,512),     xavier(Float32,3,3,512,512),   zeros(Float32,1,1,512,1),
        xavier(Float32,3,3,512,512),     xavier(Float32,3,3,512,512),   zeros(Float32,1,1,512,1),
        xavier(Float32,3,3,512,512),     xavier(Float32,3,3,512,512),   zeros(Float32,1,1,512,1),
        #Scores
        xavier(Float32,1,1,512,2),       zeros(Float32,1,1,2,1),
        #Deconv
        bilinear(Float32,32,32,2,2),     zeros(Float32,1,1,2,1),
    ]
    w[1:24] = load(pretrainedPath("DAVIS_Conv.jld"),"w")[1:24]
    return w
end

function DAVIS_LSTM_loss(w,h,c,x,y,t,D)
    total = 0.0
    for i in max(t-D,1):t
        ypred = DAVIS_LSTM(w,h,c,KnetArray(x[:,:,:,i:i]))
        total += pixelwiseSoftloss(Array(ypred),y[:,:,:,i:i])
    end
    return total
end

max_stars_repo_name: eeyrdkl/MScThesis

max_issues_repo_name: eeyrdkl/MScThesis

max_forks_repo_name: eeyrdkl/MScThesis


text:
using PyPlot
#lines
x = LinRange(0, 10, 100)
y1 = sin.(x)
y2 = cos.(x)

subplot(221)
plot(x, y1, "r")
subplot(222)
plot(x, y2, "b")

subplot(2,2, (3,4))
plot(x, y1, "r", label = "sin")
plot(x, y2, "b", label = "cos")

if !isinteractive()
    PyPlot.show()
end

#fig = Figure()
#lines(fig[1, 1], x, y1, color = :red)
#lines(fig[1, 2], x, y2, color = :blue)

#lines(fig[2, 1:2], x, y1, color = :red, label = "sin")
#lines!(fig[2, 1:2], x, y2, color = :blue, label = "cos")
#axislegend()
#current_figure()
#fig

max_stars_repo_name: mkitti/GaussianBeamLens

max_issues_repo_name: mkitti/GaussianBeamLens

max_forks_repo_name: mkitti/GaussianBeamLens


text:
mutable struct gpuCTPF <: GPUTopicModel
	K::Int
	M::Int
	V::Int
	U::Int
	N::Vector{Int}
	C::Vector{Int}
	R::Vector{Int}
	B::Int
	corp::Corpus
	batches::Vector{UnitRange{Int}}
	topics::VectorList{Int}
	scores::Matrix{Float32}
	libs::VectorList{Int}
	drecs::VectorList{Int}
	urecs::VectorList{Int}
	a::Float32
	b::Float32
	c::Float32
	d::Float32
	e::Float32
	f::Float32
	g::Float32
	h::Float32
	alef::Matrix{Float32}
	bet::Vector{Float32}
	gimel::VectorList{Float32}
	dalet::Vector{Float32}
	he::Matrix{Float32}
	vav::Vector{Float32}
	zayin::VectorList{Float32}
	het::Vector{Float32}
	phi::MatrixList{Float32}
	xi::MatrixList{Float32}
	newalef::Void
	newhe::Void
	Npsums::VectorList{Int}
	Jpsums::VectorList{Int}
	Rpsums::VectorList{Int}
	Ypsums::VectorList{Int}
	terms::VectorList{Int}
	counts::VectorList{Int}
	words::VectorList{Int}
	readers::VectorList{Int}
	ratings::VectorList{Int}
	views::VectorList{Int}
	device::cl.Device
	context::cl.Context
	queue::cl.CmdQueue
	Npsumsbuf::cl.Buffer{Int}
	Jpsumsbuf::cl.Buffer{Int}
	Rpsumsbuf::cl.Buffer{Int}
	Ypsumsbuf::cl.Buffer{Int}
	termsbuf::cl.Buffer{Int}
	countsbuf::cl.Buffer{Int}
	wordsbuf::cl.Buffer{Int}
	readersbuf::cl.Buffer{Int}
	ratingsbuf::cl.Buffer{Int}
	viewsbuf::cl.Buffer{Int}
	alefkern::cl.Kernel
	newalefkern::cl.Kernel
	betkern::cl.Kernel
	gimelkern::cl.Kernel
	daletkern::cl.Kernel
	hekern::cl.Kernel
	newhekern::cl.Kernel
	vavkern::cl.Kernel
	zayinkern::cl.Kernel
	hetkern::cl.Kernel
	phikern::cl.Kernel
	phinormkern::cl.Kernel
	xikern::cl.Kernel
	xinormkern::cl.Kernel
	alefbuf::cl.Buffer{Float32}
	newalefbuf::cl.Buffer{Float32}
	betbuf::cl.Buffer{Float32}
	gimelbuf::cl.Buffer{Float32}
	daletbuf::cl.Buffer{Float32}
	hebuf::cl.Buffer{Float32}
	newhebuf::cl.Buffer{Float32}
	vavbuf::cl.Buffer{Float32}
	zayinbuf::cl.Buffer{Float32}
	hetbuf::cl.Buffer{Float32}
	phibuf::cl.Buffer{Float32}
	xibuf::cl.Buffer{Float32}
	elbo::Float32
	newelbo::Float32

	function gpuCTPF(corp::Corpus, K::Integer, batchsize::Integer=length(corp), basemodel::Union{Void, BaseTopicModel}=nothing)
		@assert !isempty(corp)
		@assert all(ispositive.([K, batchsize]))
		checkcorp(corp)

		M, V, U = size(corp)
		N = [length(doc) for doc in corp]
		C = [size(doc) for doc in corp]
		R = [length(doc.readers) for doc in corp]

		batches = partition(1:M, batchsize)
		B = length(batches)

		@assert ispositive(U)
		@assert isequal(collect(1:U), sort(collect(keys(corp.users))))
		libs = [Int[] for _ in 1:U]
		for u in 1:U, d in 1:M
			u in corp[d].readers && push!(libs[u], d)
		end

		a, b, c, d, e, f, g, h = fill(0.1, 8)

		if isa(basemodel, Union{AbstractLDA, AbstractCTM})
			@assert isequal(size(basemodel.beta), (K, V))
			alef = exp.(basemodel.beta)
			topics = basemodel.topics
		elseif isa(basemodel, Union{AbstractfLDA, AbstractfCTM})
			@assert isequal(size(basemodel.fbeta), (K, V))
			alef = exp.(basemodel.fbeta)
			topics = basemodel.topics
		else
			alef = exp.(rand(Dirichlet(V, 1.0), K)' - 0.5)
			topics = [collect(1:V) for _ in 1:K]
		end

		bet = ones(K)
		gimel = [ones(K) for _ in 1:M]
		dalet = ones(K)
		he = ones(K, U)
		vav = ones(K)
		zayin = [ones(K) for _ in 1:M]
		het = ones(K)
		phi = [ones(K, N[d]) / K for d in batches[1]]
		xi = [ones(2K, R[d]) / 2K for d in batches[1]]

		model = new(K, M, V, U, N, C, R, B, copy(corp), batches, topics, zeros(M, U), libs, Vector[], Vector[], a, b, c, d, e, f, g, h, alef, bet, gimel, dalet, he, vav, zayin, het, phi, xi)
		fixmodel!(model, check=false)

		for (b, batch) in enumerate(batches)
			model.phi = [ones(K, N[d]) / K for d in batch]
			model.xi = [ones(2K, R[d]) / 2K for d in batch]
			updateNewELBO!(model, b)
		end
		model.phi = [ones(K, N[d]) / K for d in batches[1]]
		model.xi = [ones(2K, R[d]) / 2K for d in batches[1]]
		updateELBO!(model)
		return model
	end
end

gpuCTPF(corp::Corpus, K::Int, basemodel::Union{Void, BaseTopicModel}) = gpuCTPF(corp, K, length(corp), basemodel)

function Elogpya(model::gpuCTPF, d::Int, m::Int)
	x = 0
	readers, ratings = model.corp[d].readers, model.corp[d].ratings
	for (u, (re, ra)) in enumerate(zip(readers, ratings)), i in 1:model.K
		binom = Binomial(ra, model.xi[m][i,u])
		x += (ra * model.xi[m][i,u] * (digamma(model.gimel[d][i]) - log(model.dalet[i]) + digamma(model.he[i,re]) - log(model.vav[i]))
				- (model.gimel[d][i] / model.dalet[i]) * (model.he[i,re] / model.vav[i]) - sum([pdf(binom, y) * lgamma(y + 1) for y in 0:ra]))
	end
	return x
end

function Elogpyb(model::gpuCTPF, d::Int, m::Int)
	x = 0
	readers, ratings = model.corp[d].readers, model.corp[d].ratings
	for (u, (re, ra)) in enumerate(zip(readers, ratings)), i in 1:model.K
		binom = Binomial(ra, model.xi[m][model.K+i,u])
		x += (ra * model.xi[m][model.K+i,u] * (digamma(model.zayin[d][i]) - log(model.het[i]) + digamma(model.he[i,re]) - log(model.vav[i]))
				- (model.zayin[d][i] / model.het[i]) * (model.he[i,re] / model.vav[i]) - sum([pdf(binom, y) * lgamma(y + 1) for y in 0:ra]))
	end
	return x
end

function Elogpz(model::gpuCTPF, d::Int, m::Int)
	x = 0
	terms, counts = model.corp[d].terms, model.corp[d].counts
	for (n, (j, c)) in enumerate(zip(terms, counts)), i in 1:model.K
		binom = Binomial(c, model.phi[m][i,n])
		x += (c * model.phi[m][i,n] * (digamma(model.gimel[d][i]) - log(model.dalet[i]) + digamma(model.alef[i,j]) - log(model.bet[i]))
				- (model.gimel[d][i] / model.dalet[i]) * (model.alef[i,j] / model.bet[i]) - sum([pdf(binom, z) * lgamma(z + 1) for z in 0:c]))
	end
	return x
end

function Elogpbeta(model::gpuCTPF)
	x = model.V * model.K * (model.a * log(model.b) - lgamma(model.a))
	for j in 1:model.V, i in 1:model.K
		x += (model.a - 1) * (digamma(model.alef[i,j]) - log(model.bet[i])) - model.b * model.alef[i,j] / model.bet[i]
	end
	return x
end

function Elogptheta(model::gpuCTPF, d::Int)
	x = model.K * (model.c * log(model.d) - lgamma(model.c))
	for i in 1:model.K
		x += (model.c - 1) * (digamma(model.gimel[d][i]) - log(model.dalet[i])) - model.d * model.gimel[d][i] / model.dalet[i]
	end
	return x
end

function Elogpeta(model::gpuCTPF)
	x = model.U * model.K * (model.e * log(model.f) - lgamma(model.e))
	for u in 1:model.U, i in 1:model.K
		x += (model.e - 1) * (digamma(model.he[i,u]) - log(model.vav[i])) - model.f * model.he[i,u] / model.vav[i]
	end
	return x
end

function Elogpepsilon(model::gpuCTPF, d::Int)
	x = model.K * (model.g * log(model.h) - lgamma(model.g))
	for i in 1:model.K
		x += (model.g - 1) * (digamma(model.zayin[d][i]) - log(model.het[i])) - model.h * model.zayin[d][i] / model.het[i]
	end
	return x
end

function Elogqy(model::gpuCTPF, d::Int, m::Int)
	x = 0
	for (u, ra) in enumerate(model.corp[d].ratings)
		x -= entropy(Multinomial(ra, model.xi[m][:,u]))
	end
	return x
end

function Elogqz(model::gpuCTPF, d::Int, m::Int)
	x = 0
	for (n, c) in enumerate(model.corp[d].counts)
		x -= entropy(Multinomial(c, model.phi[m][:,n]))
	end
	return x
end

function Elogqbeta(model::gpuCTPF)
	x = 0
	for j in 1:model.V, i in 1:model.K
		x -= entropy(Gamma(model.alef[i,j], 1 / model.bet[i]))
	end
	return x
end

function Elogqtheta(model::gpuCTPF, d::Int)
	x = 0
	for i in 1:model.K
		x -= entropy(Gamma(model.gimel[d][i], 1 / model.dalet[i]))
	end
	return x
end

function Elogqeta(model::gpuCTPF)
	x = 0
	for u in 1:model.U, i in 1:model.K
		x -= entropy(Gamma(model.he[i,u], 1 / model.vav[i]))
	end
	return x
end

function Elogqepsilon(model::gpuCTPF, d::Int)
	x = 0
	for i in 1:model.K
		x -= entropy(Gamma(model.zayin[d][i], 1 / model.het[i]))
	end
	return x
end

function updateELBO!(model::gpuCTPF)
	model.elbo = model.newelbo + Elogpbeta(model) + Elogpeta(model) - Elogqbeta(model) - Elogqeta(model)
	model.newelbo = 0
	return model.elbo
end

function updateNewELBO!(model::gpuCTPF, b::Int)
	batch = model.batches[b]
	for (m, d) in enumerate(batch)
		model.newelbo += (Elogpya(model, d, m)
						+ Elogpyb(model, d, m)
						+ Elogpz(model, d, m)
						+ Elogptheta(model, d)
						+ Elogpepsilon(model, d)
						- Elogqy(model, d, m)
						- Elogqz(model, d, m)
						- Elogqtheta(model, d)
						- Elogqepsilon(model, d))
	end
end

const CTPF_ALEF_cpp =
"""
kernel void
updateAlef(long K,
			float a,
			global float *newalef,
			global float *alef)

			{
			long i = get_global_id(0);
			long j = get_global_id(1);

			alef[K * j + i] = newalef[K * j + i];
			newalef[K * j + i] = a;
			}
			"""

function updateAlef!(model::gpuCTPF)
	model.queue(model.alefkern, (model.K, model.V), nothing, model.K, model.a, model.newalefbuf, model.alefbuf)
end

const CTPF_NEWALEF_cpp =
"""
kernel void
updateNewalef(long K,
				const global long *Jpsums,
				const global long *counts,
				const global long *words,
				const global float *phi,
				global float *newalef)

				{
				long i = get_global_id(0);
				long j = get_global_id(1);

				float acc = 0.0f;

				for (long w=Jpsums[j]; w<Jpsums[j+1]; w++)
					acc += counts[words[w]] * phi[K * words[w] + i];

				newalef[K * j + i] += acc;
				}
				"""

function updateNewalef!(model::gpuCTPF)
	model.queue(model.newalefkern, (model.K, model.V), nothing, model.K, model.Jpsumsbuf, model.countsbuf, model.wordsbuf, model.phibuf, model.newalefbuf)
end

const CTPF_BET_cpp =
"""
kernel void
updateBet(long K,
			long M,
			float b,
			const global float *alef,
			const global float *gimel,
			const global float *dalet,
			global float *bet)

			{
			long i = get_global_id(0);

			float acc = 0.0f;

			for (long d=0; d<M; d++)
				acc += gimel[K * d + i];

			bet[i] = b + acc / dalet[i];
			}
			"""

function updateBet!(model::gpuCTPF)
	model.queue(model.betkern, model.K, nothing, model.K, model.M, model.b, model.alefbuf, model.gimelbuf, model.daletbuf, model.betbuf)
end

const CTPF_GIMEL_cpp =
"""
kernel void
updateGimel(long F,
			long K,
			float c,
			const global long *Npsums,
			const global long *Rpsums,
			const global long *counts,
			const global long *ratings,
			const global float *phi,
			const global float *xi,
			global float *gimel)

			{
			long i = get_global_id(0);
			long d = get_global_id(1);

			float accphi = 0.0f;
			float accxi = 0.0f;

			for (long n=Npsums[d]; n<Npsums[d+1]; n++)
				accphi += phi[K * n + i] * counts[n];

			for (long r=Rpsums[d]; r<Rpsums[d+1]; r++)
				accxi += xi[2 * K * r + i] * ratings[r];

			gimel[K * (F + d) + i] = c + accphi + accxi;
			}
			"""

function updateGimel!(model::gpuCTPF, b::Int)
	batch = model.batches[b]
	model.queue(model.gimelkern, (model.K, length(batch)), nothing, batch[1] - 1, model.K, model.c, model.Npsumsbuf, model.Rpsumsbuf, model.countsbuf, model.ratingsbuf, model.phibuf, model.xibuf, model.gimelbuf)
end

const CTPF_DALET_cpp =
"""
kernel void
updateDalet(long K,
			long V,
			long U,
			float d,
			const global float *alef,
			const global float *bet,
			const global float *he,
			const global float *vav,
			global float *dalet)

			{
			long i = get_global_id(0);

			float accalef = 0.0f;
			float acche = 0.0f;

			for (long j=0; j<V; j++)
				accalef += alef[K * j + i];

			for (long u=0; u<U; u++)
				acche += he[K * u + i];

			dalet[i] = d + accalef / bet[i] + acche / vav[i];
			}
			"""

function updateDalet!(model::gpuCTPF)
	model.queue(model.daletkern, model.K, nothing, model.K, model.V, model.U, model.d, model.alefbuf, model.betbuf, model.hebuf, model.vavbuf, model.daletbuf)
end

const CTPF_HE_cpp =
"""
kernel void
updateHe(long K,
			float e,
			global float *newhe,
			global float *he)

			{
			long i = get_global_id(0);
			long u = get_global_id(1);

			he[K * u + i] = newhe[K * u + i];
			newhe[K * u + i] = e;
			}
			"""

function updateHe!(model::gpuCTPF)
	model.queue(model.hekern, (model.K, model.U), nothing, model.K, model.e, model.newhebuf, model.hebuf)
end

const CTPF_NEWHE_cpp =
"""
kernel void
updateNewhe(long K,
			const global long *Ypsums,
			const global long *ratings,
			const global long *views,
			const global float *xi,
			global float *newhe)

			{
			long i = get_global_id(0);
			long u = get_global_id(1);

			float acc = 0.0f;

			for (long r=Ypsums[u]; r<Ypsums[u+1]; r++)
				acc += ratings[views[r]] * (xi[2 * K * views[r] + i] + xi[K * (2 * views[r] + 1) + i]);

			newhe[K * u + i] += acc;
			}
			"""

function updateNewhe!(model::gpuCTPF)
	model.queue(model.newhekern, (model.K, model.U), nothing, model.K, model.Ypsumsbuf, model.ratingsbuf, model.viewsbuf, model.xibuf, model.newhebuf)
end

const CTPF_VAV_cpp =
"""
kernel void
updateVav(long K,
			long M,
			float f,
			const global float *gimel,
			const global float *dalet,
			const global float *zayin,
			const global float *het,
			global float *vav)

			{
			long i = get_global_id(0);

			float accgimel = 0.0f;
			float acczayin = 0.0f;

			for (long d=0; d<M; d++)
			{
				accgimel += gimel[K * d + i];
				acczayin += zayin[K * d + i];
			}

			vav[i] = f + accgimel / dalet[i] + acczayin / het[i];
			}
			"""

function updateVav!(model::gpuCTPF)
	model.queue(model.vavkern, model.K, nothing, model.K, model.M, model.f, model.gimelbuf, model.daletbuf, model.zayinbuf, model.hetbuf, model.vavbuf)
end

const CTPF_ZAYIN_cpp =
"""
kernel void
updateZayin(long F,
			long K,
			float g,
			const global long *Rpsums,
			const global long *ratings,
			const global float *xi,
			global float *zayin)

			{
			long i = get_global_id(0);
			long d = get_global_id(1);

			float acc = 0.0f;

			for (long r=Rpsums[d]; r<Rpsums[d+1]; r++)
				acc += xi[K * (2 * r + 1) + i] * ratings[r];

			zayin[K * (F + d) + i] = g + acc;
			}
			"""

function updateZayin!(model::gpuCTPF, b::Int)
	batch = model.batches[b]
	model.queue(model.zayinkern, (model.K, length(batch)), nothing, batch[1] - 1, model.K, model.g, model.Rpsumsbuf, model.ratingsbuf, model.xibuf, model.zayinbuf)
end

const CTPF_HET_cpp =
"""
kernel void
updateHet(long K,
			long U,
			float h,
			const global float *he,
			const global float *vav,
			global float *het)

			{
			long i = get_global_id(0);

			float acc = 0.0f;

			for (long u=0; u<U; u++)
				acc += he[K * u + i];

			het[i] = h + acc / vav[i];
			}
			"""

function updateHet!(model::gpuCTPF)
	model.queue(model.hetkern, model.K, nothing, model.K, model.U, model.h, model.hebuf, model.vavbuf, model.hetbuf)
end

const CTPF_PHI_cpp =
"""
$(DIGAMMA_cpp)

kernel void
updatePhi(long F,
			long K,
			const global long *Npsums,
			const global long *terms,
			const global float *alef,
			const global float *bet,
			const global float *gimel,
			const global float *dalet,
			global float *phi)

			{

			long i = get_global_id(0);
			long d = get_global_id(1);

			float gdb = digamma(gimel[K * (F + d) + i]) - log(dalet[i]) - log(bet[i]);

			for (long n=Npsums[d]; n<Npsums[d+1]; n++)
				phi[K * n + i] = exp(gdb + digamma(alef[K * terms[n] + i]));
			}
			"""

const CTPF_PHI_NORM_cpp =
"""
kernel void
normalizePhi(long K,
				global float *phi)

				{
				long dn = get_global_id(0);

				float normalizer = 0.0f;

				for (long i=0; i<K; i++)
					normalizer += phi[K * dn + i];

				for (long i=0; i<K; i++)
					phi[K * dn + i] /= normalizer;
				}
				"""

function updatePhi!(model::gpuCTPF, b::Int)
	batch = model.batches[b]
	model.queue(model.phikern, (model.K, length(batch)), nothing, batch[1] - 1, model.K, model.Npsumsbuf, model.termsbuf, model.alefbuf, model.betbuf, model.gimelbuf, model.daletbuf, model.phibuf)
	model.queue(model.phinormkern, sum(model.N[batch]), nothing, model.K, model.phibuf)
end

const CTPF_XI_cpp =
"""
$(DIGAMMA_cpp)

kernel void
updateXi(long F,
			long K,
			const global long *Rpsums,
			const global long *readers,
			const global float *bet,
			const global float *gimel,
			const global float *dalet,
			const global float *he,
			const global float *vav,
			const global float *zayin,
			const global float *het,
			global float *xi)

			{
			long i = get_global_id(0);
			long d = get_global_id(1);

			float gdv = digamma(gimel[K * (F + d) + i]) - log(dalet[i]) - log(bet[i]);
			float zhv = digamma(zayin[K * (F + d) + i]) - log(het[i]) - log(vav[i]);

			for (long r=Rpsums[d]; r<Rpsums[d+1]; r++)
			{
				xi[2 * K * r + i] = exp(gdv + digamma(he[K * readers[r] + i]));
				xi[K * (2 * r + 1) + i] = exp(zhv + digamma(he[K * readers[r] + i]));
			}
			}
			"""

const CTPF_XI_NORM_cpp =
"""
kernel void
normalizeXi(long K,
			global float *xi)

				{
				long dn = get_global_id(0);

				float normalizer = 0.0f;

				for (long i=0; i<2*K; i++)
					normalizer += xi[2 * K * dn + i];

				for (long i=0; i<2*K; i++)
					xi[2 * K * dn + i] /= normalizer;
				}
				"""

function updateXi!(model::gpuCTPF, b::Int)
	batch = model.batches[b]
	model.queue(model.xikern, (model.K, length(batch)), nothing, batch[1] - 1, model.K, model.Rpsumsbuf, model.readersbuf, model.betbuf, model.gimelbuf, model.daletbuf, model.hebuf, model.vavbuf, model.zayinbuf, model.hetbuf, model.xibuf)
	model.queue(model.xinormkern, sum(model.R[batch]), nothing, model.K, model.xibuf)
end

function train!(model::gpuCTPF; iter::Int=150, tol::Real=1.0, viter::Int=10, vtol::Real=1/model.K^2, chkelbo::Int=1)
	@assert all(.!isnegative.([tol, vtol]))
	@assert all(ispositive.([iter, viter, chkelbo]))
	lowVRAM = model.B > 1

	for k in 1:iter
		chk = (k % chkelbo == 0)
		for b in 1:model.B
			for _ in 1:viter
				oldgimel = @host b model.gimelbuf
				updatePhi!(model, b)
				updateXi!(model, b)
				updateGimel!(model, b)
				updateZayin!(model, b)
				gimel = @host b model.gimelbuf
				if sum([norm(diff) for diff in oldgimel - gimel]) < length(model.batches[b]) * vtol
					break
				end
			end
			updateNewalef!(model)
			updateNewhe!(model)
			if chk
				updateHost!(model, b)
				updateNewELBO!(model, b)
			end
			lowVRAM && updateBuf!(model, b)
		end
		if checkELBO!(model, k, chk, tol)
			updateDalet!(model)
			updateHet!(model)
			updateAlef!(model)
			updateBet!(model)
			updateHe!(model)
			updateVav!(model)
			break
		end
		updateDalet!(model)
		updateHet!(model)
		updateAlef!(model)
		updateBet!(model)
		updateHe!(model)
		updateVav!(model)
	end
	updateHost!(model, 1)
	Ebeta = model.alef ./ model.bet
	model.topics = [reverse(sortperm(vec(Ebeta[i,:]))) for i in 1:model.K]

	Eeta = model.he ./ model.vav
	for d in 1:model.M
		Etheta = model.gimel[d] ./ model.dalet
		Eepsilon = model.zayin[d] ./ model.het
		model.scores[d,:] = sum(Eeta .* (Etheta + Eepsilon), 1)
	end

	model.drecs = Vector{Int}[]
	for d in 1:model.M
		nr = setdiff(keys(model.corp.users), model.corp[d].readers)
		push!(model.drecs, nr[reverse(sortperm(vec(model.scores[d,nr])))])
	end

	model.urecs = Vector{Int}[]
	for u in 1:model.U
		ur = filter(d -> !(u in model.corp[d].readers), collect(1:model.M))
		push!(model.urecs, ur[reverse(sortperm(model.scores[ur,u]))])
	end
	nothing
end


text:
using NORBiP
using Test

import SCIP
using JuMP
using LinearAlgebra
import Polyhedra
import CDDLib

test_bp() = NORBiP.BilevelProblem(
    [1.,0.], [0.],
    [-.1 0.], -ones(1,1), [-1.],
    [-1.],
    [-.1 0.], ones(1,1), ones(1,)
)

function test_bp2()
    G = ones(2,1); G[1,1] = -1
    H = 2 .* ones(2,1); H[1,1] = 4
    A = zeros(2,1) .+ [-2, 5]
    q = [11.0, 13.0]
    cx = [1.0]
    cy = [-10.0]
    B = zeros(2,1) .+ [-1, -4]
    b = [-5.0, 30.0]
    d = [1.0]
    return NORBiP.BilevelProblem(cx, cy, G, H, q, d, A, B, b)
end

@testset "Basic tests model 1" begin
    bp = test_bp()
    m = Model(() -> SCIP.Optimizer(display_verblevel = 0))
    (m, x, v) = NORBiP.high_point_relaxation(m, bp)
    optimize!(m)
    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test JuMP.value(x[1]) ≈ 0.0
    m2 = Model(() -> SCIP.Optimizer(display_verblevel = 0))
    (_, x, v) = NORBiP.bilevel_optimality(m2, bp)
    optimize!(m2)
    @test JuMP.termination_status(m2) == MOI.OPTIMAL
    @test JuMP.value(x[1]) ≈ 0.0
    @test JuMP.value(v[1]) ≈ 1.0
end

@testset "Basic tests model 2" begin
    bp = test_bp2()
    m = Model(() -> SCIP.Optimizer(display_verblevel = 0))
    (m, x, v) = NORBiP.high_point_relaxation(m, bp)
    NORBiP.add_upper_level(m, bp, x, v)
    optimize!(m)
    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test JuMP.value(x[1]) ≈ 5.0
    @test JuMP.value(v[1]) ≈ 4.0
    m = Model(() -> SCIP.Optimizer(display_verblevel = 0))
    (_, x, v) = NORBiP.bilevel_optimality(m, bp)
    @constraint(m, bp.G * x .+ bp.H .<= bp.q)
    @objective(m, Min, dot(bp.cx, x) + dot(bp.cy, v))
    optimize!(m)
    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test JuMP.value(x[1]) ≈ 1.0
    @test JuMP.value(v[1]) ≈ 3.0
    vertex_list = [
        [(zeros(2,), 4.0)], #k = 1
        [(zeros(2,), 2.0)], #k = 2
    ]
    m = Model(() -> SCIP.Optimizer(display_verblevel = 0))
    (m, x, v, λ, σ, s, upperfeas, lowerfeas, kkt, kkt2_var, kkt2_bounds, w) = NORBiP.solve_near_optimal(m, bp, NORBiP.ExtendedMethod(), 0.5, vertex_list)
    NORBiP.add_upper_level(m, bp, x, v)
    optimize!(m)
    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test JuMP.value(x[1]) ≈ 1.222 atol = 10^-3
    @test JuMP.value(v[1]) ≈ 2.555 atol = 10^-3
    empty!(m)
    (m, x, v, λ, σ, s, upperfeas, lowerfeas, kkt, kkt2_var, kkt2_bounds, w) = NORBiP.solve_near_optimal(m, bp, NORBiP.ExtendedMethod(), 1.0, vertex_list)
    NORBiP.add_upper_level(m, bp, x, v)
    optimize!(m)
    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test JuMP.value(x[1]) ≈ 1.444 atol = 10^-3
    @test JuMP.value(v[1]) ≈ 2.111 atol = 10^-3
    empty!(m)
    (m, x, v, λ, σ, s, upperfeas, lowerfeas, kkt, kkt2_var, kkt2_bounds, w) = NORBiP.solve_near_optimal(m, bp, NORBiP.ExtendedMethod(), 1.0, vertex_list, resvar=true)
    NORBiP.add_upper_level(m, bp, x, v)
    optimize!(m)
    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test JuMP.value(x[1]) ≈ 1.444 atol = 10^-3
    @test JuMP.value(v[1]) ≈ 2.111 atol = 10^-3
end

@testset "Tests model 2 with lazy constraints" begin
    bp = test_bp2()
    m = Model(() -> SCIP.Optimizer(display_verblevel = 0))
    (st, problem_ref) = NORBiP.solve_near_optimal(m, bp, NORBiP.LazyExtended(), 0.5, optimizer = () -> SCIP.Optimizer(display_verblevel = 0), poly_lib = CDDLib.Library(:exact))
    (x, v, w) = (problem_ref[:x], problem_ref[:v], problem_ref[:w])
    @test JuMP.termination_status(problem_ref[:model]) == MOI.OPTIMAL
    @test st === :OPTIMAL
    @test JuMP.value(x[1]) ≈ 1.222 atol = 10^-3
    @test JuMP.value(v[1]) ≈ 2.555 atol = 10^-3
    for k in eachindex(w)
        if !isempty(w[k])
            @test sum(JuMP.value.(w[k])) == 1
        end
    end
end

@testset "Tests model 2 with lazy constraints & batch size 1 and 5" begin
    bp = test_bp2()
    for Γ in (1, 2, 3, 10)
        m = Model(() -> SCIP.Optimizer(display_verblevel = 0))
        (st, problem_ref) = NORBiP.solve_near_optimal(m, bp, NORBiP.LazyBatched(Γ), 0.5, optimizer = () -> SCIP.Optimizer(display_verblevel = 0), poly_lib = CDDLib.Library(:exact))
        (x, v, w) = (problem_ref[:x], problem_ref[:v], problem_ref[:w])
        @test JuMP.termination_status(problem_ref[:model]) == MOI.OPTIMAL
        @test st === :OPTIMAL
        @test JuMP.value(x[1]) ≈ 1.222 atol = 10^-3
        @test JuMP.value(v[1]) ≈ 2.555 atol = 10^-3
        for k in eachindex(w)
            if !isempty(w[k])
                @test sum(JuMP.value.(w[k])) == 1
            end
        end
    end
end

@testset "Bound computations" begin
    vertices = [
        ([1,2],3),
        ([-1,-1],4),
    ]
    (α_bounds, β_bounds) = NORBiP.find_hypercube(vertices)
    @test α_bounds[1] == (-1,1)
    @test α_bounds[2] == (-1,2)
    @test β_bounds == (3,4)
end

@testset "Single-vertex heuristic" begin
    bp = test_bp2()
    for bsize in (1, 2, 5, 10)
        res = NORBiP.solve_near_optimal(bp, NORBiP.SingleVertexHeuristic(bsize), 0.5, () -> SCIP.Optimizer(display_verblevel = 0))
        @test length(res) == 5
        (st, model, x, v, C) = res
        @test st == MOI.OPTIMAL
        @test JuMP.value(x[1]) ≈ 1.222 atol = 10^-3
        @test JuMP.value(v[1]) ≈ 2.555 atol = 10^-3
    end
end

@testset "Randomized single-vertex heuristic" begin
    bp = test_bp2()
    for bsize in (1, 2, 5, 10)
        res = NORBiP.solve_near_optimal(bp, NORBiP.RandomizedSingleVertexHeuristic(bsize), 0.5, () -> SCIP.Optimizer(display_verblevel = 0))
        @test length(res) == 5
        (st, model, x, v, C) = res
        @test st == MOI.OPTIMAL
        @test JuMP.value(x[1]) ≈ 1.222 atol = 10^-3
        @test JuMP.value(v[1]) ≈ 2.555 atol = 10^-3
    end
end

@testset "Tests models with bilinear terms" begin
    bp = test_bp2()
    m = Model(() -> SCIP.Optimizer(display_verblevel = 0))
    problem_ref = NORBiP.solve_near_optimal(m, bp, NORBiP.BilinearTerms(), 0.5)
    x, v = (problem_ref[:x], problem_ref[:v])
    optimize!(m)
    @test JuMP.termination_status(problem_ref[:model]) == MOI.OPTIMAL
    @test JuMP.value(x[1]) ≈ 1.222 atol = 10^-3
    @test JuMP.value(v[1]) ≈ 2.555 atol = 10^-3

    m = Model(() -> SCIP.Optimizer(display_verblevel = 0))
    problem_ref = NORBiP.solve_near_optimal(m, bp, NORBiP.BilinearTerms(), 0.0)
    x, v = (problem_ref[:x], problem_ref[:v])
    optimize!(m)
    @test JuMP.termination_status(problem_ref[:model]) == MOI.OPTIMAL
    @test JuMP.value(x[1]) ≈ 1.0 atol = 10^-3
    @test JuMP.value(v[1]) ≈ 3.0 atol = 10^-3

    bp = test_bp()
    m = Model(() -> SCIP.Optimizer(display_verblevel = 0))
    problem_ref = NORBiP.solve_near_optimal(m, bp, NORBiP.BilinearTerms(), 0.1)
    x, v = (problem_ref[:x], problem_ref[:v])
    optimize!(m)
    @test JuMP.termination_status(problem_ref[:model]) == MOI.OPTIMAL
    @test JuMP.value(x[1]) ≈ 0.5 atol = 10^-3
    @test JuMP.value(v[1]) ≈ 1.05 atol = 10^-3
end

@testset "Adversarial test model 2" begin
    bp = test_bp2()
    m = Model(() -> SCIP.Optimizer(display_verblevel = 0))
    (_, x, v) = NORBiP.bilevel_optimality(m, bp)
    NORBiP.add_upper_level(m, bp, x, v)
    optimize!(m)
    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test JuMP.value(x[1]) ≈ 1.0
    @test JuMP.value(v[1]) ≈ 3.0
    ks = NORBiP.near_optimal_robust_violated(bp, JuMP.value.(x), JuMP.value.(v), 0.5; optimizer=SCIP.Optimizer)
    @test length(ks) == 1
end

@testset "Valid inequality" begin
    # Basic tests model 2
    bp = test_bp2()
    vertex_list = [
        [(zeros(2,), 4.0)], #k = 1
        [(zeros(2,), 2.0)], #k = 2
    ]
    m = Model(() -> SCIP.Optimizer(display_verblevel = 0))
    (m, x, v, λ, σ, s, upperfeas, lowerfeas, kkt, kkt2_var, kkt2_bounds, w) = NORBiP.solve_near_optimal(m, bp, NORBiP.ExtendedMethod(), 0.5, vertex_list)
    NORBiP.add_upper_level(m, bp, x, v)
    optimize!(m)
    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test JuMP.value(x[1]) ≈ 1.222 atol = 10^-3
    @test JuMP.value(v[1]) ≈ 2.555 atol = 10^-3
    obj1 = objective_value(m)
    NORBiP.add_upper_level(m, bp, x, v)
    optimize!(m)
    @test obj1 ≈ objective_value(m) atol = 10^-3
    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test JuMP.value(x[1]) ≈ 1.222 atol = 10^-3
    @test JuMP.value(v[1]) ≈ 2.555 atol = 10^-3
    (st, Abounds) = NORBiP.compute_bounds(bp; optimizer=() -> SCIP.Optimizer(display_verblevel = 0))
    @test st == MOI.OPTIMAL
    m = Model(() -> SCIP.Optimizer(display_verblevel = 0))
    (m, x, v, λ, σ, s, upperfeas, lowerfeas, kkt, kkt2_var, kkt2_bounds, w) = NORBiP.solve_near_optimal(m, bp, NORBiP.ExtendedMethod(), 0.5, vertex_list)
    NORBiP.add_upper_level(m, bp, x, v)
    NORBiP.add_bound_inequalities(m, bp, x, v, λ, Abounds)
    optimize!(m)
    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test obj1 ≈ objective_value(m) atol = 10^-3
    @test JuMP.value(x[1]) ≈ 1.222 atol = 10^-3
    @test JuMP.value(v[1]) ≈ 2.555 atol = 10^-3
    m = Model(() -> SCIP.Optimizer(display_verblevel = 0))
    (m, x, v, λ, σ, s, upperfeas, lowerfeas, kkt, kkt2_var, kkt2_bounds, w) = NORBiP.solve_near_optimal(m, bp, NORBiP.ExtendedMethod(), 1.0, vertex_list)
    NORBiP.add_upper_level(m, bp, x, v)
    NORBiP.add_bound_inequalities(m, bp, x, v, λ, Abounds)
    optimize!(m)
    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test JuMP.value(x[1]) ≈ 1.444 atol = 10^-3
    @test JuMP.value(v[1]) ≈ 2.111 atol = 10^-3
end

include("test_inf.jl")


text:
module Apply

export apply, offset_apply
using Base: tail

apply(f::F, a::Tuple{A,Vararg}, args...) where {F,A} =
    (f(a[1], args...), apply(f, tail(a), args...)...)
apply(f::F, a::Tuple{A,Vararg}, b::Tuple{B,Vararg}, args...) where {F,A,B} =
    (f(a[1], b[1], args...), apply(f, tail(a), tail(b), args...)...)
apply(f::F, a::Tuple{A,Vararg}, b::Tuple{B,Vararg}, c::Tuple{C,Vararg}, args...) where {F,A,B,C} =
    (f(a[1], b[1], c[1], args...), apply(f, tail(a), tail(b), tail(c), args...)...)

apply(f, ::Tuple{}, ::Tuple{}, ::Tuple{}, args...) = ()
apply(f, a, ::Tuple{}, ::Tuple{}, args...) = ()
apply(f, ::Tuple{}, b, ::Tuple{}, args...) = ()
apply(f, a, b, ::Tuple{}, args...) = ()
apply(f, ::Tuple{}, ::Tuple{}, args...) = ()
apply(f, a, ::Tuple{}, args...) = ()
apply(f, ::Tuple{}, args...) = ()

offset_apply(f, o::Tuple{O,Vararg}, a::AbstractArray, offset::Int, args...) where O = begin
    offset = f(o[1], a, offset, args...)
    offset_apply(f, tail(o), a, offset::Int, args...)
    nothing
end
offset_apply(f, a::AbstractArray, o::Tuple{O,Vararg}, offset::Int, args...) where O = begin
    offset = f(a, o[1], offset, args...)
    offset_apply(f, a, tail(o), offset::Int, args...)
    nothing
end
offset_apply(f, o::Tuple{O,Vararg}, offset::Int, args...) where O = begin
    offset = f(o[1], offset, args...)
    offset_apply(f, tail(o), offset::Int, args...)
    nothing
end
offset_apply(f, o::Tuple{}, offset::Int, args...) = nothing
offset_apply(f, a::AbstractArray, o::Tuple{}, offset::Int, args...) = nothing
offset_apply(f, o::Tuple{}, a::AbstractArray, offset::Int, args...) = nothing

end # module


text:
using Distributions
using Random

struct SpreadingParams
  dist::Truncated{Pareto{Float64},Continuous,Float64}
  spreading::Vector{Float32}
end

function SpreadingParams(rng::AbstractRNG, N::Integer; alpha::Real, x0::Real=1, truncation::Real=Inf)
  dist=truncated(Pareto(alpha, x0), x0, truncation)
  spreading = Vector{Float32}(undef, N)
  rand!(rng, dist, spreading)
  SpreadingParams(dist, spreading)
end

function saveparams(dict, p::SpreadingParams, prefix::AbstractString="")
  dict[prefix*"dist"] = params(p.dist)
  dict[prefix*"values"] = p.spreading
  nothing
end

function show(io::IO, p::SpreadingParams)
  print(io, "Spreading params for ", numindividuals(p), " individuals, distribution=", p.dist)
end

spreading(p::SpreadingParams, source_id::Integer) = Float64(p.spreading[source_id])


text:
# Chemfiles.jl, a modern library for chemistry file reading and writing
# Copyright (C) Guillaume Fraux and contributors -- BSD license

export lengths, set_lengths!, angles, set_angles!
export shape, set_shape!, CellShape, volume, matrix, wrap!

__ptr(cell::UnitCell) = __ptr(cell.__handle)
__const_ptr(cell::UnitCell) = __const_ptr(cell.__handle)

"""
The possible shape for an unit cell are:

- ``Chemfiles.Orthorhombic`` for unit cells with the three angles are 90°.
- ``Chemfiles.Triclinic`` for unit cells where the three angles may not be 90°.
- ``Chemfiles.Infinite`` for unit cells without boundaries.
"""
@enum CellShape begin
    Infinite = lib.CHFL_CELL_INFINITE
    Orthorhombic = lib.CHFL_CELL_ORTHORHOMBIC
    Triclinic = lib.CHFL_CELL_TRICLINIC
end

"""
Create an ``UnitCell`` from the three lenghts, with all the angles equal to 90°.
"""
function UnitCell(a::Number, b::Number, c::Number)
    ptr = @__check_ptr(lib.chfl_cell(Float64[a, b, c]))
    return UnitCell(CxxPointer(ptr, is_const=false))
end

"""
Create an ``UnitCell`` from the three lenghts and three angles.
"""
function UnitCell(a::Number, b::Number, c::Number, α::Number, β::Number, γ::Number)
    ptr = @__check_ptr(lib.chfl_cell_triclinic(Float64[a, b, c], Float64[α, β, γ]))
    return UnitCell(CxxPointer(ptr, is_const=false))
end

"""
Get a copy of the ``UnitCell`` of a ``frame``.
"""
function UnitCell(frame::Frame)
    ptr = @__check_ptr(lib.chfl_cell_from_frame(__const_ptr(frame)))
    cell = UnitCell(CxxPointer(ptr, is_const=false))
    copy = deepcopy(cell)
    finalize(cell)
    return copy
end

"""
Get the unit ``cell`` volume.
"""
function volume(cell::UnitCell)
    result = Ref{Float64}(0)
    __check(lib.chfl_cell_volume(__const_ptr(cell), result))
    return result[]
end

"""
Get the three ``cell`` lengths (a, b and c) in angstroms.
"""
function lengths(cell::UnitCell)
    result = Float64[0, 0, 0]
    __check(lib.chfl_cell_lengths(__const_ptr(cell), result))
    return result
end

"""
Set the ``cell`` lengths to ``a``, ``b`` and ``c``.

``a``, ``b`` and ``c`` should be in angstroms.
"""
function set_lengths!(cell::UnitCell, a::Real, b::Real, c::Real)
    __check(lib.chfl_cell_set_lengths(__ptr(cell), Float64[a, b, c]))
    return nothing
end

"""
Get the three ``cell`` angles (alpha, beta and gamma) in degrees.
"""
function angles(cell::UnitCell)
    result = Float64[0, 0, 0]
    __check(lib.chfl_cell_angles(__const_ptr(cell), result))
    return result
end

"""
Set the `cell` angles to ``α``, ``β`` and ``γ``.

``α``, ``β`` and ``γ`` should be in degrees.
"""
function set_angles!(cell::UnitCell, α::Real, β::Real, γ::Real)
    __check(lib.chfl_cell_set_angles(__ptr(cell), Float64[α, β, γ]))
    return nothing
end

"""
Get the ``cell`` matricial representation, *i.e.* the representation of the
three base vectors as::

        | a_x   b_x   c_x |
        |  0    b_y   c_y |
        |  0     0    c_z |
"""
function matrix(cell::UnitCell)
    result = Array{Float64}(undef, 3, 3)
    __check(lib.chfl_cell_matrix(__const_ptr(cell), pointer(result)))
    return result
end

"""
Get the ``cell`` shape, as a ``CellShape`` value.
"""
function shape(cell::UnitCell)
    result = Ref{lib.chfl_cellshape}(0)
    __check(lib.chfl_cell_shape(__const_ptr(cell), result))
    return CellShape(result[])
end

"""
Set the ``cell`` shape to the given ``shape``.
"""
function set_shape!(cell::UnitCell, shape::CellShape)
    __check(lib.chfl_cell_set_shape(__ptr(cell), lib.chfl_cellshape(shape)))
    return nothing
end

"""
Wrap a `vector` in the unit `cell`.
"""
function wrap!(cell::UnitCell, vector::Vector{Float64})
    if length(vector) != 3
        throw(ChemfilesError("Can only wrap vectors containing 3 elements in an unit cell"))
    end
    __check(lib.chfl_cell_wrap(__const_ptr(cell), vector))
    return vector
end

"""
Make a deep copy of a ``cell``.
"""
function Base.deepcopy(cell::UnitCell)
    ptr = lib.chfl_cell_copy(__const_ptr(cell))
    return UnitCell(CxxPointer(ptr, is_const=false))
end


text:
function to_prior(pos,prior)
    return quantile(prior, .5+.5*pos)
end

function to_unit_ball(pos,prior)
    return (cdf(prior,pos)-.5)/.5
end

function box_bound!(pos,box)
    if !(all(box[:,1].<pos.<box[:,2]))
        pos[pos.<box[:,1]].=box[:,1][pos.<box[:,1]]
        pos[pos.>box[:,2]].=box[:,2][pos.>box[:,2]]
    end
end