text:
function init_pa_shim()
    libdir = joinpath(dirname(@__FILE__), "..", "deps", "usr", "lib")
    libsuffix = ""
    basename = "pa_shim"
    @static if is_linux() && Sys.ARCH == :x86_64
        libsuffix = "x86_64-linux-gnu"
    elseif is_linux() && Sys.ARCH == :i686
        libsuffix = "i686-linux-gnu"
    elseif is_apple() && Sys.ARCH == :x86_64
        libsuffix = "x86_64-apple-darwin14"
    elseif is_windows() && Sys.ARCH == :x86_64
        libsuffix = "x86_64-w64-mingw32"
    elseif is_windows() && Sys.ARCH == :i686
        libsuffix = "i686-w64-mingw32"
    elseif !any(
            (sfx) -> isfile(joinpath(libdir, "$basename.$sfx")),
            ("so", "dll", "dylib"))
        error("Unsupported platform $(Sys.MACHINE). You can build your own library by running `make` from $(joinpath(@__FILE__, "..", "deps", "src"))")
    end
    # if there's a suffix-less library, it was built natively on this machine,
    # so load that one first, otherwise load the pre-built one
    global const libpa_shim = Base.Libdl.find_library(
            [basename, "$(basename)_$libsuffix"],
            [libdir])
    libpa_shim == "" && error("Could not load $basename library, please file an issue at https://github.com/JuliaAudio/RingBuffers.jl/issues with your `versioninfo()` output")
    shim_dlib = Libdl.dlopen(libpa_shim)
    # pointer to the shim's process callback
    global const shim_processcb_c = Libdl.dlsym(shim_dlib, :pa_shim_processcb)
    if shim_processcb_c == C_NULL
        error("Got NULL pointer loading `pa_shim_processcb`")
    end
end

const pa_shim_errmsg_t = Cint
const PA_SHIM_ERRMSG_OVERFLOW = Cint(0) # input overflow
const PA_SHIM_ERRMSG_UNDERFLOW = Cint(1) # output underflow
const PA_SHIM_ERRMSG_ERR_OVERFLOW = Cint(2) # error buffer overflowed

# This struct is shared with pa_shim.c
mutable struct pa_shim_info_t
    inputbuf::Ptr{PaUtilRingBuffer} # ringbuffer for input
    outputbuf::Ptr{PaUtilRingBuffer} # ringbuffer for output
    errorbuf::Ptr{PaUtilRingBuffer} # ringbuffer to send error notifications
    sync::Cint # keep input/output ring buffers synchronized (0/1)
    notifycb::Ptr{Void} # Julia callback to notify on updates (called from audio thread)
    inputhandle::Ptr{Void} # condition to notify on new input data
    outputhandle::Ptr{Void} # condition to notify when ready for output
    errorhandle::Ptr{Void} # condition to notify on new errors
end

"""
    PortAudio.shimhash()

Return the sha256 hash(as a string) of the source file used to build the shim.
We may use this sometime to verify that the distributed binary stays in sync
with the rest of the package.
"""
shimhash() = unsafe_string(
        ccall((:pa_shim_getsourcehash, libpa_shim), Cstring, ()))
Base.unsafe_convert(::Type{Ptr{Void}}, info::pa_shim_info_t) = pointer_from_objref(info)


text:
using Kwant, PyPlot, PyCall
kwant = Kwant; pyplot = PyPlot

####
include("band_structure.jl")

lead = make_lead().finalized()
kwant.plotter.bands(lead, show=false)
pyplot.xlabel("momentum [(lattice constant)^-1]")
pyplot.ylabel("energy [t]")
pyplot.gcf()

####

include("closed_system.jl")

syst = make_system()

# Check that the system looks as intended.
kwant.plot(syst)

# Finalize the system.
system = syst.finalized()

# We should observe energy levels that flow towards Landau
# level energies with increasing magnetic field.
plot_spectrum(system, [iB * 0.002 for iB in range(0,length=(100-1))])

# Plot an eigenmode of a circular dot. Here we create a larger system for
# better spatial resolution.
system = make_system(r=20).finalized()

plot_wave_function(system)
plot_current(system)

ext: jl


text:
"""
    space(::SiteType"Qudit";
          dim = 2,
          conserve_qns = false,
          conserve_number = false,
          qnname_number = "Number")

Create the Hilbert space for a site of type "Qudit".

Optionally specify the conserved symmetries and their quantum number labels.
"""
function ITensors.space(
  ::SiteType"Qudit";
  dim=2,
  conserve_qns=false,
  conserve_number=conserve_qns,
  qnname_number="Number",
)
  if conserve_number
    return [QN(qnname_number, n - 1) => 1 for n in 1:dim]
  end
  return dim
end

function ITensors.val(::ValName{N}, ::SiteType"Qudit") where {N}
  return parse(Int, String(N)) + 1
end

function ITensors.state(::StateName{N}, ::SiteType"Qudit", s::Index) where {N}
  n = parse(Int, String(N))
  st = zeros(dim(s))
  st[n + 1] = 1.0
  return itensor(st, s)
end

# one-body operators
function _op(::OpName"Id", ::SiteType"Qudit"; dim::Tuple=(2,))
  d = dim[1]
  mat = zeros(d, d)
  for k in 1:d
    mat[k, k] = 1.0
  end
  return mat
end

function _op(::OpName"Adag", ::SiteType"Qudit"; dim::Tuple=(2,))
  d = dim[1]
  mat = zeros(d, d)
  for k in 1:(d - 1)
    mat[k + 1, k] = √k
  end
  return mat
end
_op(::OpName"adag", st::SiteType"Qudit"; kwargs...) = _op(OpName"Adag"(), st; kwargs...)
_op(::OpName"a†", st::SiteType"Qudit"; kwargs...) = _op(OpName"Adag"(), st; kwargs...)

function _op(::OpName"A", ::SiteType"Qudit"; dim::Tuple=(2,))
  d = dim[1]
  mat = zeros(d, d)
  for k in 1:(d - 1)
    mat[k, k + 1] = √k
  end
  return mat
end
_op(::OpName"a", st::SiteType"Qudit"; kwargs...) = _op(OpName"A"(), st; kwargs...)

function _op(::OpName"N", ::SiteType"Qudit"; dim::Tuple=(2,))
  d = dim[1]
  mat = zeros(d, d)
  for k in 1:d
    mat[k, k] = k - 1
  end
  return mat
end
_op(::OpName"n", st::SiteType"Qudit"; kwargs...) = _op(OpName"N"(), st; kwargs...)

# two-body operators
function _op(::OpName"ab", st::SiteType"Qudit"; dim::Tuple=(2, 2))
  return kron(_op(OpName("a"), st; dim=(dim[1],)), _op(OpName("a"), st; dim=(dim[2],)))
end

function _op(::OpName"a†b", st::SiteType"Qudit"; dim::Tuple=(2, 2))
  return kron(_op(OpName("a†"), st; dim=(dim[1],)), _op(OpName("a"), st; dim=(dim[2],)))
end

function _op(::OpName"ab†", st::SiteType"Qudit"; dim::Tuple=(2, 2))
  return kron(_op(OpName("a"), st; dim=(dim[1],)), _op(OpName("a†"), st; dim=(dim[2],)))
end

function _op(::OpName"a†b†", st::SiteType"Qudit"; dim::Tuple=(2, 2))
  return kron(_op(OpName("a†"), st; dim=(dim[1],)), _op(OpName("a†"), st; dim=(dim[2],)))
end

# interface
function op(on::OpName, st::SiteType"Qudit", s::Index...)
  rs = reverse([s...])
  d⃗ = dim.(Tuple(rs))
  opmat = _op(on, st; dim=d⃗)
  return ITensors.itensor(opmat, prime.(rs)..., dag.(rs)...)
end

# Zygote
@non_differentiable _op(::OpName"ab", ::SiteType"Qudit")
@non_differentiable _op(::OpName"a†b", ::SiteType"Qudit")
@non_differentiable _op(::OpName"ab†", ::SiteType"Qudit")
@non_differentiable _op(::OpName"a†b†", ::SiteType"Qudit")
@non_differentiable _op(::OpName"a", ::SiteType"Qudit")
@non_differentiable _op(::OpName"a†", ::SiteType"Qudit")
@non_differentiable _op(::OpName"N", ::SiteType"Qudit")


text:
# Trait to determine if an Index, Index collection, Tensor, or ITensor
# has symmetries
abstract type SymmetryStyle end

symmetrystyle(T) = error("No SymmetryStyle defined for the specified type")

symmetrystyle(T, S, U, V...)::SymmetryStyle = (
  Base.@_inline_meta; symmetrystyle(symmetrystyle(T), symmetrystyle(S, U, V...))
)

symmetrystyle(T, S)::SymmetryStyle = symmetrystyle(symmetrystyle(T), symmetrystyle(S))

# Rules for basic collections
symmetrystyle(inds::Tuple) = symmetrystyle(inds...)
symmetrystyle(inds::AbstractVector)::SymmetryStyle = symmetrystyle(inds...)


text:
function __init__()
    ScientificTypesBase.set_convention(DefaultConvention())
    ScientificTypesBase.TRAIT_FUNCTION_GIVEN_NAME[:table] = Tables.istable
end


text:
using Base

mutable struct PriceLevels
    levels::Array{PriceLevel, 1}
    ask::UInt128
    bid::UInt128
    max_price::UInt128
    min_price::UInt128

    function PriceLevels(max_price::UInt128, min_price::UInt128)
        n_prices = max_price - min_price + 1
        new(Array{PriceLevel, 1}(undef, n_prices), max_price,
            min_price, max_price, min_price)::PriceLevels
    end
end

function Base.size(pls::PriceLevels)::Tuple{UInt128}
    return size(pls.levels)
end

function price_to_index(price::UInt128, min_price::UInt128)::Int
    return Int(price - min_price + 1)
end

function index_to_price(index::UInt128, min_price::UInt128)::UInt128
    return index + min_price - 1
end

function insert_order!(pls::PriceLevels, order::LimitOrder)

    price_index = price_to_index(order.price, pls.min_price)
    if !isassigned(pls.levels, Int(price_index))
        pls.levels[price_index] = PriceLevel()
    end

    if side(order) == BUY

        # Assert price lower than ask, might be removed
        @assert order.price < pls.ask "Inserted buy order with larger price than current ask"

        add!(get_level(pls, order.price), order)

        if order.price > pls.bid
            pls.bid = order.price
        end
    else
        # Assert price higher than bid, might be removed
        @assert order.price > pls.bid "Inserted sell order with lower price than current bid"
        add!(get_level(pls, order.price), order)

        if order.price < pls.ask
            pls.ask = order.price
        end

    end

end

function get_level(pls::PriceLevels, price::UInt128)::PriceLevel
    price_index = price_to_index(price, pls.min_price)
    return pls.levels[price_index]
end

function update_ask!(pls::PriceLevels)::UInt128
    ask = pls.ask + 1
    while ask < pls.max_price
        ask_index = price_to_index(ask, pls.min_price)
        if isassigned(pls.levels, ask_index)
            if pls.levels[ask_index].size > 0
                pls.ask = ask
                return ask
            end
        end
        ask = ask + 1
    end
    pls.ask = ask
    return pls.ask
end

function update_bid!(pls::PriceLevels)::UInt128
    bid = pls.bid - 1
    while bid > pls.min_price
        bid_index = price_to_index(bid, pls.min_price)
        if isassigned(pls.levels, bid_index)
            if pls.levels[bid_index].size > 0
                pls.bid = bid
                return bid

            end
        end
        bid = bid - 1
    end
    pls.bid = bid
    return pls.bid
end


text:
module geometry

export is_in_pipe, is_in_half_space, is_in_sphere

# Returns distance from cylinder center, squared
# query: DenseVector for which to compute cylindrical norm
# d₀: Direction of the cylinder
# center: Center of cylinder goes through this point
function cylindrical_norm_squared{T}(query::DenseVector{T}, d₀::DenseVector{T}, center::DenseVector{T})
    @assert length(d₀) == length(center) == length(query)
    @assert norm(d₀, 1) > 1e-8
    q₀ = query - center
    q⟂ = q₀ - (q₀⋅d₀)d₀ / (d₀⋅d₀)
    q⟂⋅q⟂
end

# True if query point is inside cylinder
# query: DenseVector for which to compute cylindrical norm
# d₀: Direction of the cylinder
# center: Center of cylinder goes through this point
# radius: of the cylinder
function is_in_pipe{T}(query::DenseVector{T}, d₀::DenseVector{T}, center::DenseVector{T}, radius::T)
    cylindrical_norm_squared(query, d₀, center) <= (radius^2)
end

# True if in half-space defined by normal and offset
# query: DenseVector for which to compute condition
# d₀: Normal to the plane defining the half-space
# origin: Origin of the plane defining the half-space
function is_in_half_space{T}(query::DenseVector{T}, d₀::DenseVector{T}, origin::DenseVector{T})
    @assert length(d₀) == length(origin) == length(query)
    @assert norm(d₀, 1) > 1e-8
    (query - origin)⋅d₀ >= 0
end

# True if in sphere defined by center and radius
# query: DenseVector for which to compute condition
function is_in_sphere{T}(query::DenseVector{T}, center::DenseVector{T}, radius::T)
    @assert length(center) == length(query)
    q₀ = center - query
    q₀⋅q₀ <= radius^2
end

end # module geometry


text:
using MultiChannelFactorAnalysis
using Test
using Statistics
using LinearAlgebra
#using CSV

const datadir = joinpath(dirname(@__FILE__), "..", "data")

# @testset "MCFMSim" begin
#     fc, fus = simFactor(3, [4, 5, 6])
#     @test length(fc) == 3
#     @test length(fus)== 3
#     @test length(fus[1]) == 4
#     Ls = [5, 3, 4]
#     nc = 3
#     nus = [2, 2, 3]
#     Hs, Gs = randomLoadings(Ls, nc, nus)
#     @test length(Hs) == length(Ls)
#     @test size(Hs[1])[2] == nc
#     @test sum([size(H)[1] for H in Hs]) == sum(Ls)
# end

# @testset "MCFM" begin
#     Ls = [20, 15, 10]
#     N  = 100
#     p = 2
#     ps = [4, 3, 2]
#     lf = (l, p) ->  randomLoadingNorm(l, p; target_tr=0.3)
#     Hs, Gs = randomLoadings(Ls, p, ps, lf)
#     edists = vcat([repeat([Normal(0, sqrt(0.4 / L))], L) for L in Ls]...)
# end

@testset "All tests" begin

    @testset "layout" begin
        # Channel Layout
        clay = MultiChannelLayout([5, 3, 1])
        @test clay.nobs == 9
        @test clay.nchannels == 3
        @test_throws ArgumentError MultiChannelLayout([-2, 0, 1])
        # Factor Layout
        flay = MultiChannelFactorLayout(3, [2, 1, 1])
        @test flay.nchannels == 3
        @test flay.nchannelspecific == 4
        @test flay.nfactors == 7
    end
    @testset "utils" begin
        m = [1 2 3; 4 5 6; 7 8 9]
        b1 = [[1 2 3], [4 5 6; 7 8 9]]
        b2 = [[1 2; 4 5], reshape([9], 1, 1)]
        @test extract_horizontal_blocks(m, [1, 2]) == b1
        @test extract_diagonal_blocks(m, [2, 1], [2, 1]) == b2
        m[1, 1] = 300
        @test b1[1][1] == 1
        @test extract_horizontal_blocks(m, [1, 2]; get_view=true)[1][1] == 300
    end

    @testset "MCFMStructs" begin
        c = 3
        n = 10
        nobs_per_chan = [2, 2, 1]
        vals = [[zeros(p) for p in nobs_per_chan] for i in 1:n]
        channellayout = MultiChannelLayout(nobs_per_chan)
        @test MultiChannelData(vals, channellayout, n).values == vals
        @test_throws ArgumentError MultiChannelData(vals[1:(end-1), :], channellayout, n)
        @test_throws ArgumentError MultiChannelData(vals, MultiChannelLayout([2, 1, 1]), n)
        @test_throws ArgumentError MultiChannelData(vals, channellayout, 0 )
        @test_throws ArgumentError MultiChannelData(vals, channellayout, n+1)
    end

    @testset "data" begin
        c = 2
        n = 2
        nobs_per_chan = [2, 2]
        vals = [[j.* ones(p) .+ i for (j, p) in enumerate(nobs_per_chan)] for i in 1:n]
        channellayout = MultiChannelLayout(nobs_per_chan)
        data = MultiChannelData(vals, channellayout, n)
        @test stack(data) == [2.0 3.0; 2.0 3.0; 3.0 4.0; 3.0 4.0]
        @test unstack(stack(data), channellayout).values == data.values
        @test (data + data).values[1][1][1] == 4.0

        #from long, no missing
        df = DataFrame(sample_id = [1,  1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],
                channel = ["a", "a", "b", "b", "a", "a", "b", "b", "a", "a", "b", "b"],
                observation_id = [1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2 ],
                value = [0.5, 0.3, 0.2, 0.1, 0.3, 0.4, 0.6, 0.2, 0.3, 0.5, 0.2, 0.1])
        data = data_fromlong(df)
        @test data.values[1][1][1] == 0.5
        @test data.values[2][2][2] == 0.2

        #from long, some missing to trim
        df2 = DataFrame(sample_id = [1, 1, 1, 1, 2, 2, 3, 3],
                channel = ["a", "a", "b", "b", "a", "b", "a", "b"],
                observation_id = [1, 2, 1, 2, 1, 1, 1, 1],
                value = [0.5, 0.3, 0.2, 0.1, 0.3, 0.4, 0.6, 0.3])
        @test dat2.channellayout.nobs== 6

        # from long, some missing to interpolate
        df3 = DataFrame(sample_id = [1,  1, 1, 1, 2, 2, 2, 2, 3, 3],
                channel = ["a", "a", "b", "b", "a", "a", "b", "b", "a", "b"],
                observation_id = [1, 2, 1, 2, 1, 2, 1, 2, 1, 1 ],
                value = [0.5, 0.3, 0.2, 0.1, 0.3, 0.4, 0.6, 0.2, 0.3, 0.2])
        dat3 = data_fromlong(df3; max_missingpcent =0.5)
        @test df3.values[3][1][2] == 0.5
        @test dat3.nobs == 12

    end

    @testset "simulation" begin
        factorlayout = MultiChannelFactorLayout(1, [1, 1])
        channellayout = MultiChannelLayout([5, 5])
        nsamples = 100
        @test periodic_factors(nsamples, factorlayout).factors[1, :] == sin.((1:nsamples)/nsamples * 2 * pi)
        gauss_fac = classical_factors(nsamples, factorlayout)
        @test abs(mean(gauss_fac.factors[1, :])) <= 0.2
        @test size(gauss_fac.factors) == (3, 100)
        @test std(gauss_fac.factors[1, :]) <= 1.2
        @test std(gauss_fac.factors[1, :]) >= 0.8
        errors = indep_gauss_errors(nsamples, channellayout)
        @test size(errors.values)[1] == nsamples

        #Simulation
        H = reshape([1; 0.5; -0.5; 0.2; -0.2; -1; 0.2; 2.0; 1.0; 0.2], (10, 1))
        G = [0.5 0; 0.2 0; 0.1 0; -0.1 0; 0.1 0; 0 0.5; 0 0.3; 0 0.2; 0 -0.1; 0 -1]
        Sig = diagm(repeat([0.1], 10))
        params = MCFMParams(factorlayout, channellayout, H, G, Sig)
        simdata = simulate(nsamples, params, classical_factors, indep_gauss_errors)
        @test abs(mean(stack(simdata)[1, :])) <= 0.25
        @test abs(std(stack(simdata)[1, :]) - sqrt(1.1)) <= 0.35
    end

    @testset "factortest" begin
        M1 = [3; 2; 1; 2; 3; 4; 2; 1.0; 3.0; 2.0] * transpose([3; 2; 1; 2; 3; 4; 2; 1.0; 3.0; 2.0])
        M2 = M +  [5.0; -3.0; .2; 0.6; 0.3; 4.1; 3.2; -2.6; 4.0; 5.0] * transpose([5.0; -3.0; .2; 0.6; 0.3; 4.1; 3.2; -2.6; 4.0; 5.0])
        M3b = randn(100, 30)
        M3  = M3b * transpose(M3b)
        @test eigenvalue_ratio_test(reverse(eigvals(M1))) == 1
        @test eigenvalue_ratio_test(reverse(eigvals(M2))) == 2
        @test eigenvalue_ratio_test(reverse(eigvals(M3)); threshold=5) = 0

    @testset "fitting" begin
        factorlayout = MultiChannelFactorLayout(1, [1, 1])
        channellayout = MultiChannelLayout([5, 5])
        nsamples = 1000
        H = reshape([1; 0.5; -0.5; 0.2; -0.2; -1; 0.2; 2.0; 1.0; 0.2], (10, 1))
        G = [0.5 0; 0.2 0; 0.1 0; -0.1 0; 0.1 0; 0 0.5; 0 0.3; 0 0.2; 0 -0.1; 0 -1]
        Sig = diagm(repeat([0.1], 10))
        params = MCFMParams(factorlayout, channellayout, H, G, Sig)
        simdata = simulate(nsamples, params, classical_factors, indep_gauss_errors)
        fitted = fit(simdata,  factorlayout; maxiter=1e2)
        fitH = fitted.finalparams.H
        fitG = fitted.finalparams.G
        fitSig = fitted.finalparams.Σ
        @test norm(transpose(fitH) * fitH - transpose(H) * H)^2/norm(transpose(H) * H)^2 <= 1
        @test norm(transpose(fitG) * fitG - transpose(G) * G)^2/norm(transpose(G) * G)^2 <= 1
        #@test norm(fitSig - Sig)^2/norm(Sig)^2 <= 1

        # 0 channel-specific factors in some channels
        factorlayout = MultiChannelFactorLayout(1, [2, 0])
        channellayout = MultiChannelLayout([5, 5])
        nsamples = 1000
        H = reshape([1; 0.5; -0.5; 0.2; -0.2; -1; 0.2; 2.0; 1.0; 0.2], (10, 1))
        G = [0.5 0; 0.2 -0.1; 0.1 0.1; -0.1 2.0; 0.1 -2.0; 0 0; 0 0; 0 0; 0.0 0.0; 0.0 0.0]
        Sig = diagm(repeat([0.1], 10))
        params = MCFMParams(factorlayout, channellayout, H, G, Sig)
        simdata = simulate(nsamples, params, classical_factors, indep_gauss_errors)
        fitted = fit(simdata,  factorlayout; maxiter=1e2)
        fitH = fitted.finalparams.H
        fitG = fitted.finalparams.G
        fitSig = fitted.finalparams.Σ
        @test norm(transpose(fitH) * fitH - transpose(H) * H)^2/norm(transpose(H) * H)^2 <= 1
        @test norm(transpose(fitG) * fitG - transpose(G) * G)^2/norm(transpose(G) * G)^2 <= 1

        # Common factors only
        factorlayout = MultiChannelFactorLayout(2, [0, 0])
        channellayout = MultiChannelLayout([5, 5])
        nsamples = 1000
        H = [1 0.0; 0.5 0.2; -0.5 -0.1; 0.2 -0.3; -0.2 0.4; -1 -0.2; 0.2 0.6; 2.0 -0.2; 1.0 0.4; 0.2 -0.2]
        G = zeros(10, 0)
        Sig = diagm(repeat([0.1], 10))
        params = MCFMParams(factorlayout, channellayout, H, G, Sig)
        simdata = simulate(nsamples, params, classical_factors, indep_gauss_errors)
        fitted = fit(simdata,  factorlayout; maxiter=1e2)
        fitH = fitted.finalparams.H
        fitG = fitted.finalparams.G
        fitSig = fitted.finalparams.Σ
        facs = extract_factors(simdata, fitted.finalparams)
        #@test norm(fitSig - Sig)^2/norm(Sig)^2 <= 1
    end

    @testset "extraction" begin
        factorlayout = MultiChannelFactorLayout(1, [1, 1])
        channellayout = MultiChannelLayout([5, 5])
        nsamples = 100
        H = reshape([1; 0.5; -0.5; 0.2; -0.2; -1; 0.2; 2.0; 1.0; 0.2], (10, 1))
        G = [0.5 0; 0.2 0; 0.1 0; -0.1 0; 0.1 0; 0 0.5; 0 0.3; 0 0.2; 0 -0.1; 0 -1]
        Sig = diagm(repeat([0.1], 10))
        params = MCFMParams(factorlayout, channellayout, H, G, Sig)
        simdata = simulate(nsamples, params, periodic_factors, indep_gauss_errors)
        true_fac = periodic_factors(nsamples, factorlayout )
        fitted_fac = extract_factors(simdata, params)
        @test norm(true_fac.factors - fitted_fac.factors)^2/norm(true_fac.factors)^2 <= 1

        errorless_preds   = predict(true_fac, params)
        preds = predict(fitted_fac, params)
        @test norm(errorless_preds - stack(simdata))^2/norm(stack(simdata)) <= 1
        @test norm(preds - stack(simdata))^2/norm(stack(simdata)) <= 1
    end

    @testset "writeout" begin
        factorlayout = MultiChannelFactorLayout(1, [1, 1])
        channellayout = MultiChannelLayout([5, 5])
        nsamples = 100
        H = reshape([1; 0.5; -0.5; 0.2; -0.2; -1; 0.2; 2.0; 1.0; 0.2], (10, 1))
        G = [0.5 0; 0.2 0; 0.1 0; -0.1 0; 0.1 0; 0 0.5; 0 0.3; 0 0.2; 0 -0.1; 0 -1]
        Sig = diagm(repeat([0.1], 10))
        params = MCFMParams(factorlayout, channellayout, H, G, Sig)
        param_long = to_long(params)
        @test size(parm_long) == (20, 4)
        fac = periodic_factors(nsamples, factorlayout )
        @test size(to_long(fac)) == (100, 3)

        df = DataFrame(sample_id = [1,  1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],
                channel = ["a", "a", "b", "b", "a", "a", "b", "b", "a", "a", "b", "b"],
                observation_id = [1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2 ],
                value = [0.5, 0.3, 0.2, 0.1, 0.3, 0.4, 0.6, 0.2, 0.3, 0.5, 0.2, 0.1])
        data = data_fromlong(df)
end


text:
"""
    set_dual_model_sense!(dual_model::MOI.ModelLike, model::MOI.ModelLike)

Set the dual model objective sense.
"""
function set_dual_model_sense(dual_model::MOI.ModelLike, primal_model::MOI.ModelLike)::Nothing where T
    # Get model sense
    primal_sense = MOI.get(primal_model, MOI.ObjectiveSense())
    if primal_sense == MOI.FEASIBILITY_SENSE
        error(primal_sense, " is not supported") # Feasibility should be supported?
    end
    # Set dual model sense
    dual_sense = (primal_sense == MOI.MIN_SENSE) ? MOI.MAX_SENSE : MOI.MIN_SENSE
    MOI.set(dual_model, MOI.ObjectiveSense(), dual_sense)
    return
end

function _scalar_quadratic_function(func::MOI.ScalarQuadraticFunction{T}, ::Type{T}) where T
    return MOIU.canonical(func)
end
function _scalar_quadratic_function(func::MOI.ScalarAffineFunction{T}, ::Type{T}) where T
    return _scalar_quadratic_function(
        SQF{T}(func.terms, MOI.ScalarQuadraticTerm{T}[], func.constant), T)
end
function _scalar_quadratic_function(func::MOI.SingleVariable, T::Type)
    return _scalar_quadratic_function(SAF{T}(func), T)
end

# Primals
"""
    PrimalObjective{T}

Primal objective is defined as a `MOI.ScalarAffineFunction`
"""
mutable struct PrimalObjective{T}
    obj::SQF{T}
    quad_cross_parameters::Dict{VI, Vector{MOI.ScalarAffineTerm{T}}}
    obj_parametric::Union{SQF{T},Nothing}

    function PrimalObjective{T}(obj) where T
        canonical_obj = _scalar_quadratic_function(obj, T)
        # if isempty(canonical_obj.terms)
        #     error("Dualization does not support models with no variables in the objective function.")
        # end
        # This was commented for now, because the current understanding is that
        # problems like {min 0*x} are well defined and have well defined dual problems.
        # Therefore, they present no issue to dualization as opposed to problems
        # with FEASIBILITY_SENSE that do not have a well defined dual problem.
        # Moreover, JuMP and MOI default is FEASIBILITY_SENSE, if a MIN_SENSE
        # is in the problem, it is because the user set it explicitly.
        # For more on the original discussion, see:
        # https://github.com/JuliaOpt/Dualization.jl/pull/64#discussion_r347484642
        quad_cross_parameters = Dict{VI, Vector{MOI.ScalarAffineTerm{T}}}()
        return new(canonical_obj,quad_cross_parameters,nothing)
    end
end

# Duals
"""
    DualObjective{T}

Dual objective is defined as a `MOI.ScalarAffineFunction`.
"""
mutable struct DualObjective{T}
    obj::SQF{T}

end

const Objective{T} = Union{PrimalObjective{T}, DualObjective{T}}

function get_raw_obj(objective::Objective{T}) where T
    return objective.obj
end
function get_affine_terms(objective::Objective{T}) where T
    return objective.obj.affine_terms
end

function get_primal_objective(primal_model::MOI.ModelLike, T::Type=Float64)
    F = MOI.get(primal_model, MOI.ObjectiveFunctionType())
    return _get_primal_objective(MOI.get(primal_model, MOI.ObjectiveFunction{F}()), T)
end

function _get_primal_objective(obj_fun, T::Type)
    return PrimalObjective{T}(obj_fun)
end

# allow removing variables from objective function
function get_primal_objective(primal_model::MOI.ModelLike, variable_parameters::Vector{VI}, T::Type)
    p_obj = get_primal_objective(primal_model, T)
    if length(variable_parameters) > 0
        vars_func, quad_cross_params, params_func = split_variables(p_obj.obj, variable_parameters)
        p_obj.obj = vars_func
        p_obj.quad_cross_parameters = quad_cross_params
        p_obj.obj_parametric = params_func
    end
    return p_obj
end

function split_variables(func::MOI.ScalarQuadraticFunction{T},
    variable_parameters::Vector{VI}) where T

    # linear part
    lin_params = MOI.ScalarAffineTerm{T}[]
    lin_vars = MOI.ScalarAffineTerm{T}[]
    for term in func.affine_terms
        if term.variable_index in variable_parameters
            push!(lin_params, term)
        else
            push!(lin_vars, term)
        end
    end

    # Quadratic part
    quad_params = MOI.ScalarQuadraticTerm{T}[]
    quad_vars = MOI.ScalarQuadraticTerm{T}[]
    quad_cross_params = Dict{VI, Vector{MOI.ScalarAffineTerm{T}}}()
    for term in func.quadratic_terms
        is_param_1 = term.variable_index_1 in variable_parameters
        is_param_2 = term.variable_index_2 in variable_parameters
        if is_param_1 && is_param_2
            push!(quad_params, term)
        elseif is_param_1
            push_affine_term(quad_cross_params, term, false)
        elseif is_param_2
            push_affine_term(quad_cross_params, term, true)
        else
            push!(quad_vars, term)
        end
    end

    variables_func = MOI.ScalarQuadraticFunction{T}(lin_vars, quad_vars, func.constant)
    parameters_func = MOI.ScalarQuadraticFunction{T}(lin_params, quad_params, zero(T))

    return variables_func, quad_cross_params, parameters_func
end

function push_affine_term(dic, term::MOI.ScalarQuadraticTerm{T}, var_is_first::Bool) where T
    variable = var_is_first ? term.variable_index_1 : term.variable_index_2
    parameter = var_is_first ? term.variable_index_2 : term.variable_index_1
    if haskey(dic, variable)
        push!(dic[variable],
            MOI.ScalarAffineTerm{T}(term.coefficient, parameter))
    else
        dic[variable] =
            [MOI.ScalarAffineTerm{T}(term.coefficient, parameter)]
    end
end

"""
    set_dual_objective(dual_model::MOI.ModelLike, dual_objective::DualObjective{T})::Nothing where T

Add the objective function to the dual model.
"""
function set_dual_objective(dual_model::MOI.ModelLike,
    dual_objective::DualObjective{T})::Nothing where T
    # Set dual model objective function
    raw_obj = get_raw_obj(dual_objective)
    if MOIU.number_of_quadratic_terms(T, raw_obj) > 0
        MOI.set(dual_model, MOI.ObjectiveFunction{SQF{T}}(), raw_obj)
    else
        MOI.set(dual_model, MOI.ObjectiveFunction{SAF{T}}(),
            SAF{T}(raw_obj.affine_terms, raw_obj.constant))
    end
    return
end

"""
    get_dual_objective(dual_model::MOI.ModelLike, dual_obj_affine_terms::Dict,
                       primal_objective::PrimalObjective{T})::DualObjective{T} where T

build the dual model objective function from the primal model.
"""
function get_dual_objective(dual_problem, dual_obj_affine_terms::Dict,
                            primal_objective::PrimalObjective{T},
                            con_types, scalar_affine_terms,
                            variable_parameters)::DualObjective{T} where T

    dual_model = dual_problem.dual_model
    map = dual_problem.primal_dual_map
    sense_change = MOI.get(dual_model, MOI.ObjectiveSense()) == MOI.MAX_SENSE ? -one(T) : one(T)

    # standard linear part
    num_objective_terms = length(dual_obj_affine_terms)
    lin_terms = MOI.ScalarAffineTerm{T}[]
    sizehint!(lin_terms, num_objective_terms)
    for var in keys(dual_obj_affine_terms) # Number of constraints of the primal model
        coef = dual_obj_affine_terms[var]
        push!(lin_terms, MOI.ScalarAffineTerm{T}(
            # Add positive terms bi if dual model sense is max
            sense_change * coef,
            # Variable index associated with term bi
            var
            ))
    end

    # standard quadratic part
    quad_terms = MOI.ScalarQuadraticTerm{T}[]
    sizehint!(quad_terms, length(primal_objective.obj.quadratic_terms))
    for term in primal_objective.obj.quadratic_terms
        push!(quad_terms, MOI.ScalarQuadraticTerm{T}(
            - MOI.coefficient(term),
            map.primal_var_dual_quad_slack[term.variable_index_1],
            map.primal_var_dual_quad_slack[term.variable_index_2]
            ))
    end

    # parametric part

    if nothing !== primal_objective.obj_parametric

        # linear
        for term in primal_objective.obj_parametric.affine_terms
            push!(lin_terms, MOI.ScalarAffineTerm{T}(
                MOI.coefficient(term),
                map.primal_parameter[term.variable_index]
                ))
        end

        # quadratic
        for term in primal_objective.obj_parametric.quadratic_terms
            push!(quad_terms, MOI.ScalarQuadraticTerm{T}(
                MOI.coefficient(term),
                map.primal_parameter[term.variable_index_1],
                map.primal_parameter[term.variable_index_2],
                ))
        end

        # crossed
        # TODO? set_dot
        for vi in variable_parameters
            param = map.primal_parameter[vi]
            for term in scalar_affine_terms[vi]
                push!(quad_terms, MOI.ScalarQuadraticTerm{T}(
                    sense_change * MOI.coefficient(term),
                    param,
                    term.variable_index,
                    ))
            end
        end
    end

    saf_dual_objective = MOI.ScalarQuadraticFunction{T}(
        lin_terms, quad_terms, MOI.constant(get_raw_obj(primal_objective)))
    return DualObjective{T}(saf_dual_objective)
end


text:
function compute_dodge(x, i, n; width, space = 0.2)
    unique_x = unique(sort(x))
    width === automatic && (width = minimum(diff(unique_x))*(1-space))
    w = width/n
    x′ = x .+ i*w .- w*(n+1)/2
    return x′, w
end


text:
using SparseArrays
using LinearAlgebra
using Makie
using LidJul
using DataStructures #for DefaultDict
using PrettyTables
using StaticArrays
using Interpolations
using SIMD

#MIT18086_NAVIERSTOKES
#    Solves the incompressible Navier-Stokes equations in a
#    rectangular domain with prescribed velocities along the
#    boundary. The solution method is finite differencing on
#    a staggered grid with implicit diffusion and a Chorin
#    projection method for the pressure.
#    Visualization is done by a colormap-isoline plot for
#    pressure and normalized quiver and streamline plot for
#    the velocity field.
#    The standard setup solves a lid driven cavity problem.
# 07/2007 by Benjamin Seibold
#            http://www-math.mit.edu/~seibold/
# Feel free to modify for teaching and learning.
const   ARRAY=Array
const   REAL=Float64
function mit18086_navierstokes()
   # scene=Scene(resolution=(1000,1000),scale=false)
   # display(scene)
   #-----------------------------------------------------------------------
   # CuArrays.allowscalar(false)

   mte=DefaultDict{String,Float64}(0.0)
   mte["initialization"]+=@elapsed   begin

   Re = 1.e4     # Reynolds number
   dt = 0.5e-2;    # time step
   tf = 50e-0;    # final time
   lx = 1.0;       # width of box
   ly = 1.0;       # height of box
   nx = 128 ;      # number of x-gridpoints
   ny = 128 ;      # number of y-gridpoints
   #nx = 41;
   #ny = 41;
   nsteps = 120;  # number of steps with graphic output
   #-----------------------------------------------------------------------
   nt = ceil(tf/dt); #dt = tf/nt;
   @show nt,dt

   x = range(0,stop=lx,length=nx+1); hx = lx/nx;
   y = range(0,stop=ly,length=ny+1); hy = ly/ny;

   println("initialization")

   L2dp=Laplacian2D(nx,ny,1,1,neumann,neumann,neumann,neumann)
   gmgp=PoissonGMG(L2dp,GSSmoother)
   Tlp=TensorialOperator(ARRAY,REAL,nx,ny,hx,hy,3,3,1.0,1.0,1.e-5,false)
   Tlu=TensorialOperator(ARRAY,REAL,nx-1,ny,hx,hy,2,3,dt/Re,dt/Re,1)
   Tlv=TensorialOperator(ARRAY,REAL,nx,ny-1,hx,hy,3,2,dt/Re,dt/Re,1)
   Tlq=TensorialOperator(ARRAY,REAL,nx-1,ny-1,hx,hy,2,2,1,1,0)

   Ue=zeros(nx+2,ny+2)
   Ve=zeros(nx+2,ny+2)

   uN = collect(x*0.0 .+ 1.0);  vN=zeros(nx) #vN = avg(x)*0;
   uS = collect(x*0.0);  vS=zeros(nx) #vS = avg(x)*0;
   uW = zeros(ny); vW = collect(y*0.0);
   uE = zeros(ny); vE = collect(y*0.0);

   Ue[1,2:ny+1].=uW
   Ue[nx+1,2:ny+1].=uE
   Ue[1:nx+1,1].=2.0.*uS.-Ue[1:nx+1,2]
   Ue[1:nx+1,end].=2.0.*uN.-Ue[1:nx+1,end-1]

   Ve[2:end-1,1].=vS
   Ve[2:end-1,end].=vN
   Ve[1,1:ny+1].=2.0.*vW.-Ve[2,1:ny+1]
   Ve[end,1:ny+1].=2.0.*vE.-Ve[end-1,1:ny+1]

   U=view(Ue,2:nx,2:ny+1)
   V=view(Ve,2:nx+1,2:ny)

   UVx=zeros(nx,ny)
   UVy=zeros(nx,ny)
   UVx2=zeros(nx,ny)
   UVy2=zeros(nx,ny)

   vUVy=view(UVy,2:nx,1:ny)
   vUVx=view(UVx,1:nx,2:ny)
   vUVx2=view(UVx2,1:nx-1,1:ny)
   vUVy2=view(UVy2,1:nx,1:ny-1)

   Bx = zeros(nx+1,ny); By = zeros(nx,ny+1); Bxy=zeros(nx,ny)

   Bx[1,:] .= uW
   Bx[end,:] .= uE
   By[:,1] .= vS
   By[:,end] .= vN

   bU=zeros(nx-1,ny)
   bV=zeros(nx,ny-1)
   Qbxy=zeros(nx-1,ny-1)
   P=zeros(nx,ny) ;  Q = zeros(nx-1,ny-1)

   end

   mte["init Makie"]+=@elapsed   begin

   scene,qn,un,vn,uvn=plotflow_makie(nx,ny,hx,hy,x,y,U,V,Ue,Ve,P,Q,Tlq)

   end
   display(scene)
   record(scene,"test.gif",framerate=6) do io

   hxm1=1/hx
   hym1=1/hy

      @inbounds for k = 1:nt

      # treat nonlinear terms
      mte["gamma"]+=@elapsed begin
         # mu=maxabs(U)
         # mv=maxabs(V)
         mu=vmaxabs(Ue,SIMD.Vec{16,Float64})
         mv=vmaxabs(Ve,SIMD.Vec{16,Float64})

         gamma = min(1.2*dt*max(mu/hx,mv/hy),1.0);
      end

      mte["Ue boundary"]+=@elapsed begin
         @simd for i=1:nx+1
            Ue[i,1]=2.0uS[i]-Ue[i,2]
            Ue[i,end]=2.0uN[i]-Ue[i,end-1]
         end
      end

      mte["Ve boundary"]+=@elapsed begin
         @simd for j=1:ny+1
            Ve[1,j]=2.0vW[j]-Ve[2,j]
            Ve[end,j]=2.0vE[j]-Ve[end-1,j]
         end
      end

      compute_derivative!(UVx, UVx2, UVy, UVy2,Ue,Ve,nx,ny,gamma,hxm1,hym1,mte)

      mte["bu"]+=@elapsed begin
         for j=1:ny
            for i=1:nx-1
               bU[i,j]=U[i,j]-dt*(vUVy[i,j]+vUVx2[i,j]);
            end
         end
         border=(dt/Re)*2/(hx^2)
         for i=1:nx-1
            bU[i,ny] +=  border
         end
      end

      mte["bv"]+=@elapsed   begin
         for j=1:ny-1
            for i=1:nx
               bV[i,j]=V[i,j]-dt*(vUVx[i,j]+vUVy2[i,j]);
            end
         end

      end
      mte["solve TT"]+=@elapsed solve!(U,bU,Tlu)
      mte["solve TT"]+=@elapsed solve!(V,bV,Tlv)

         # pressure correction
      mte["Bxy"]+=@elapsed   begin
         j=1
         i=1
         Bxy[i,j]=-((U[i,j]-Bx[i,j])*hxm1+(V[i,j]-By[i,j])*hym1)
         for i=2:nx
            Bxy[i,j]=-((U[i,j]-U[i-1,j])*hxm1+(V[i,j]-By[i,j])*hym1)
         end
         i=1
         for j=2:ny
               Bxy[i,j]=-((U[i,j]-Bx[i,j])*hxm1+(V[i,j]-V[i,j-1])*hym1)
         end
         for j=2:ny
            for i=2:nx

               Bxy[i,j]=-((U[i,j]-U[i-1,j])*hxm1+(V[i,j]-V[i,j-1])*hym1)
            end
         end
      end

      mte["solveMG"]+=@elapsed LidJul.solve!(P,Bxy,gmgp)

      # @time begin
      mte["U V"]+=@elapsed begin
         @simd for j=1:ny
            for i=1:nx-1

               U[i,j] -= (P[i+1,j]-P[i,j])*hxm1
            end
         end
         @inbounds for j=1:ny-1
            for i=1:nx
               V[i,j] -= (P[i,j+1]-P[i,j])*hym1
            end
         end
      end

      # visualization
      mte["Makie Frames"]+=@elapsed begin
      if floor(25*k/nt)>floor(25*(k-1)/nt)
         print('.')
      end
      #
      if k==1 || floor(nsteps*k/nt)>floor(nsteps*(k-1)/nt)
         mte["solve Qbxy"]+=@elapsed begin
         @inbounds for j=1:ny-1
            for i=1:nx-1
               Qbxy[i,j]=(U[i,j+1]-U[i,j])*hym1-(V[i+1,j]-V[i,j])*hxm1
            end
         end
         end
         mte["solve TTQ"]+=@elapsed solve!(Q,Qbxy,Tlq)

         Qb=zeros(nx+1,ny+1)
         Qb[2:end-1,2:end-1].=Q
         qn[] = copy!(to_value(qn),Qb)
         uvn[] = copy!(to_value(uvn),[Ue,Ve])
         recordframe!(io) # record a new frame

      end
      end
   end
end
    sort_and_dump(mte)
   return
end

# function make_grad_qn(vq)

#    function grad_qn(x,y)
#       # @show x,y
#       nxq,nyq=size(vq)
#       # @show nxq,nyq
#       xi=1+x*(nxq-1)
#       yi=1+y*(nyq-1)
#       # @show xi,yi
#        f_qn=interpolate(vq,BSpline(Linear()))
#        Point2f0(Interpolations.gradient(f_qn,xi,yi))
#    end
#    return grad_qn
# end

function make_velocity(uv)
      f_vx=interpolate(uv[1],BSpline(Quadratic(Reflect(OnCell()))))
      f_vy=interpolate(uv[2],BSpline(Quadratic(Reflect(OnCell()))))
   function velocity(x,y)

      nx,ny=size(uv[1])
      xi=1+x*(nx-1)
      yi=1+y*(ny-1)

      Point2f0(f_vx(xi,yi),f_vy(xi,yi))
   end
   return velocity
end

function plotflow_makie(nx,ny,hx,hy,x,y,U,V,Ue,Ve,P,Q,Tlq)
   scene=Scene(resolution=(300,300),scale=false)
   limits = FRect(0.2,0.2, 0.72, 0.72)

   sp=2
   stp=1
   endp=2

   qn=Node(rand(nx+1,ny+1))

   uvn=Node([rand(size(Ue)...),rand(size(Ve)...)])
   un=Node(rand(size(Ue)...))
   vn=Node(rand(size(Ve)...))
   # lift(a->println(size(a),norm(a)),qn)
   as=2

   nlevels=10
   contour!(scene,x,y,lift(a->a,qn),linewidth=1,levels=200,fillrange=false,limits=limits,colormap = :plasma, axis = (showticks = false, showaxis = false,), alpha=0.95)
   streamplot!(scene,lift(a->make_velocity(to_value(a)),uvn),x,y, colormap = :magma, arrow_size=0.02, gridsize=(20,20))

   axis = scene[Axis]
   axis[:names][:axisnames] = (" ", " ")

   scene,qn,un,vn,uvn
end

function maxabs(U)
   minu=typemax(eltype(U))
   maxu=typemin(eltype(U))
   for uij in U
      minu > uij && (minu=uij)
      maxu < uij && (maxu=uij)
   end
   max(abs(minu),abs(maxu))
end

# function vmaxabs(xin, ::Type{SIMD.Vec{N,T}}=SIMD.Vec{4,Float64}) where {N, T}
#LP: needs more tests
function vmaxabs(xin, ::Type{SIMD.Vec{N,T}}) where {N, T}
   l=length(xin)
   xs=reshape(xin,(l,))
   lane = VecRange{N}(0)

   xmv= SIMD.Vec{N,T}(typemin(T))
   nb=div(l,N)
   is=1
   @inbounds for ib in 1:nb
      xmv=max(xmv,abs(xs[lane + is]))
      is+=N
   end
   xm=maximum(xmv)
   @inbounds for i=is:l
      xm=max(xm,abs(xs[is]))
   end
   xm
end

function sort_and_dump(dict)
   s=reverse(sort(collect(dict), by=x->x[2]))
   println("")
   t=collect(x[2] for x in s)
   # @show t
   total=sum(t)
   @show total
   for si in s
      si[2]>0.001total && println(si)
   end
end

const BS=4

#LP: The ugly part which should be enhanced via V. Churavy techniques
function compute_derivative!(UVx, UVx2, UVy, UVy2,Ue,Ve,nx,ny,gamma,hxm1,hym1,mte)

      njb=div(ny,BS)
      nib=div(nx,BS)

   mte["UVx static"]+=@elapsed begin
      for jb=1:njb
         jd=BS*(jb-1)

         @inbounds for ib=1:nib
            id=BS*(ib-1)
            # @show id,jd

            # @inbounds ue00=SMatrix{BS,BS,Float64}((Ue[id+1
            # [Ue[i,j] for i=ia+1:ia+BS, j=ja+1:ja+BS]
            ue00=@inbounds @SMatrix [Ue[ic+id+0,jc+jd+0] for ic=1:BS,jc=1:BS]
            ue10=@inbounds @SMatrix [Ue[ic+id+1,jc+jd+0] for ic=1:BS,jc=1:BS]
            ue01=@inbounds @SMatrix [Ue[ic+id+0,jc+jd+1] for ic=1:BS,jc=1:BS]
            ue11=@inbounds @SMatrix [Ue[ic+id+1,jc+jd+1] for ic=1:BS,jc=1:BS]
            ue21=@inbounds @SMatrix [Ue[ic+id+2,jc+jd+1] for ic=1:BS,jc=1:BS]
            ue02=@inbounds @SMatrix [Ue[ic+id+0,jc+jd+2] for ic=1:BS,jc=1:BS]

            ve00=@inbounds @SMatrix [Ve[ic+id+0,jc+jd+0] for ic=1:BS,jc=1:BS]
            ve10=@inbounds @SMatrix [Ve[ic+id+1,jc+jd+0] for ic=1:BS,jc=1:BS]
            ve20=@inbounds @SMatrix [Ve[ic+id+2,jc+jd+0] for ic=1:BS,jc=1:BS]
            ve01=@inbounds @SMatrix [Ve[ic+id+0,jc+jd+1] for ic=1:BS,jc=1:BS]
            ve11=@inbounds @SMatrix [Ve[ic+id+1,jc+jd+1] for ic=1:BS,jc=1:BS]
            ve12=@inbounds @SMatrix [Ve[ic+id+1,jc+jd+2] for ic=1:BS,jc=1:BS]

            uay00=0.5.*(ue01.+ue00)
            udy00=0.5.*(ue01.-ue00)
            uay01=0.5.*(ue02.+ue01)
            udy01=0.5.*(ue02.-ue01)
            uay10=0.5.*(ue11.+ue10)
            uax00=0.5.*(ue11.+ue01)
            udx00=0.5.*(ue11.-ue01)
            uax10=0.5.*(ue21.+ue11)
            udx10=0.5.*(ue21.-ue11)
         #
            vax00=0.5.*(ve10.+ve00)
            vdx00=0.5.*(ve10.-ve00)
            vax10=0.5.*(ve20.+ve10)
            vdx10=0.5.*(ve20.-ve10)
            vax01=0.5.*(ve11.+ve01)
            vay10=0.5.*(ve11.+ve10)
            vdy10=0.5.*(ve11.-ve10)
            vay11=0.5.*(ve12.+ve11)
            vdy11=0.5.*(ve12.-ve11)
            #
            uva100=uay00.*vax00.-gamma.*abs.(uay00).*vdx00
            uva110=uay10.*vax10.-gamma.*abs.(uay10).*vdx10
            uvx=hxm1.*(uva110.-uva100)

            uva2200=vay10.*vay10.-gamma.*abs.(vay10).*vdy10
            uva2201=vay11.*vay11.-gamma.*abs.(vay11).*vdy11
            uvy2=hym1.*(uva2201.-uva2200)

            uva200=uay00.*vax00.-gamma.*abs.(vax00).*udy00
            uva201=uay01.*vax01.-gamma.*abs.(vax01).*udy01
            uvy=hym1.*(uva201.-uva200)

            uva1200=uax00.*uax00.-gamma.*abs.(uax00).*udx00
            uva1210=uax10.*uax10.-gamma.*abs.(uax10).*udx10
            uvx2=hxm1.*(uva1210.-uva1200)

            # # #
            for jc=1:BS

               @simd for ic=1:BS
                  @inbounds UVx[id+ic,jd+jc]=uvx[ic,jc]
               end
               @simd for ic=1:BS
                  @inbounds UVy[id+ic,jd+jc]=uvy[ic,jc]
               end
               @simd for ic=1:BS
                  @inbounds UVx2[id+ic,jd+jc]=uvx2[ic,jc]
               end
               @simd for ic=1:BS
                  @inbounds UVy2[id+ic,jd+jc]=uvy2[ic,jc]
               end
            end

         end
      end
   end

end

@time mit18086_navierstokes()


text:
block_eltype(b) = eltype(b)
block_getindex(b, args...) = getindex(b, args...)

struct BlockMap{T, N, B, R<:AbstractArray{B,N}, BS<:NTuple{N,AbstractUnitRange{Int}}} <: AbstractBlockArray{T,N}
    blocks::R
    axes::BS
    BlockMap(blocks::R, block_sizes::BS) where {N, B, R<:AbstractArray{B,N}, BS} =
        new{block_eltype(B), N, B, R, BS}(blocks, block_sizes)
end

@inline Base.axes(block_array::BlockMap) = block_array.axes

isblockzero(A::BlockMap, args...) = false

isblockentryzero(b, args...) = false
isblockentryzero(b::Diagonal, i, j) = i ≠ j
isblockentryzero(b::UpperTriangular, i, j) = i > j
isblockentryzero(b::LowerTriangular, i, j) = i < j

zero_map(::Type{B}, sizes...) where B =
    zeros(eltype(B), sizes...)

@inline function getblock(block_arr::BlockMap{T,N,B}, block::Vararg{Integer, N}) where {T,N,B}
    @boundscheck blockcheckbounds(block_arr, block...)
    if isblockzero(block_arr, block...)
        zero_map(B,blocksize(block_arr.block_sizes, block)...)
    else
        block_arr.blocks[block...]
    end
end

@inline function Base.getindex(block_arr::BlockMap{T,N}, blockindex::BlockIndex{N}) where {T,N}
    @boundscheck blockcheckbounds(block_arr, Block(blockindex.I))
    @inbounds block = getblock(block_arr, blockindex.I...)
    @boundscheck checkbounds(block, blockindex.α...)
    @inbounds v = block[blockindex.α...]
    return v
end

function _check_setblock!(block_arr::BlockMap{T, N}, v, block::NTuple{N, Integer}) where {T,N}
    bls = blocklengths.(axes(block_arr))
    for i in 1:N
        if size(v, i) != bls[i][block[i]]
            throw(DimensionMismatch(string("tried to assign $(size(v)) array to ", getindex.(bls, block), " block")))
        end
    end
end

@inline function setblock!(block_arr::BlockMap{T, N}, v, block::Vararg{Integer, N}) where {T,N}
    @boundscheck blockcheckbounds(block_arr, block...)
    @boundscheck _check_setblock!(block_arr, v, block)
    @inbounds block_arr.blocks[block...] = v
    return block_arr
end

Base.@propagate_inbounds function Base.setindex!(block_array::BlockMap{T, N}, v, block_index::BlockIndex{N}) where {T,N}
    getblock(block_array, block_index.I...)[block_index.α...] = v
end

# * Display

## structured matrix methods ##
function Base.replace_in_print_matrix(A::BlockMap, i::Integer, j::Integer, s::AbstractString)
    I,J = findblockindex.(axes(A), (i,j))
    bi,bj = Int.(block.((I,J)))
    !isblockzero(A,I,J) && !isblockentryzero(A.blocks[bi,bj], blockindex.((I,J))...)  ? s : Base.replace_with_centered_mark(s)
end

block_spy_block(b) = "■"
block_spy_block(::Diagonal) = "⟍"
block_spy_block(::LowerTriangular) = "◣"
block_spy_block(::UpperTriangular) = "◥"

function block_spy(io::IO, A::BlockMap{<:Any,2})
    dm,dn = blocklengths.(axes(A))
    pxl = reduce(gcd, vcat(dm, dn))

    for (I,i) in enumerate(dm)
        for i′ in div(i, pxl)
            for (J,j) in enumerate(dn)
                s = if isblockzero(A, I, J)
                    Base.replace_with_centered_mark(" ")
                else
                    block_spy_block(A.blocks[I,J])
                end
                for j′ in div(j, pxl)
                    write(io, " $s")
                end
            end
            write(io, "\n")
        end
    end
end

block_spy(A) = block_spy(stdout, A)

export BlockMap


text:
using Test

@time @testset "Materials" begin
    include("test_materials.jl")
    include("test_materials2.jl")
    include("test_materials3.jl")
    include("test_materials4.jl")
    include("test_materials5.jl")
    include("test_materials6.jl")
end

@time @testset "Operations" begin
    include("test_operations.jl")
end

true


text:
module TestDfExtensions

using Base.Test
using DataArrays
using DataFrames
using Dates

using TimeData
## using TimeData

println("Running DataFrame extension tests")

######################
## composeDataFrame ##
######################

df = DataFrame()
df[:a] = @data([4, 5])
df[:b] = @data([3, 8])

vals = [4 3; 5 8]
df2 = TimeData.composeDataFrame(vals, [:a, :b])
@test isequal(df, df2)

#############################
## conversion to DataFrame ##
#############################

df = DataFrame()
df[:a] = @data([4, 5])
df[:b] = @data([3, 8])
dats = [Date(2014,1,1):Date(2014,1,2)]
tn = TimeData.Timenum(df, dats)

df2 = DataFrame()
df2[:idx] = dats
df2[:a] = @data([4, 5])
df2[:b] = @data([3, 8])

@test isequal(df2, TimeData.convert(DataFrame, tn))

###########
## round ##
###########

df = DataFrame()
df[:a] = @data([3.44444, 5., 2.88888, NA])
df[:b] = @data(["he", "ll", NA, "o"])

@test isequal(round(df), round(df, 2))

df2 = DataFrame()
df2[:a] = @data([3.44, 5., 2.89, NA])
df2[:b] = @data(["he", "ll", NA, "o"])

@test isequal(round(df, 2), df2)

###############
## setDfRow! ##
###############

## with Array
##-----------

df = DataFrame()
df[:a] = @data([4, 5])
df[:b] = @data([3, 8])

df2 = DataFrame()
df2[:a] = @data([4, 10])
df2[:b] = @data([3, 10])

@test isequal(TimeData.setDfRow!(df, Any[10 10], 2), df2)

df2 = DataFrame()
df2[:a] = @data([4, NA])
df2[:b] = @data([3, 10])

@test isequal(TimeData.setDfRow!(df, Any[NA 10], 2), df2)
@test_throws ErrorException TimeData.setDfRow!(df, [1 2 3], 2)

## with DataFrame
##---------------

df = DataFrame()
df[:a] = @data([4, 5])
df[:b] = @data([3, 8])

df2 = DataFrame()
df2[:a] = @data([4, 10])
df2[:b] = @data([3, 10])

@test isequal(TimeData.setDfRow!(df, df2[2, :], 2), df2)

df2 = DataFrame()
df2[:a] = @data([4, NA])
df2[:b] = @data([3, 10])

@test isequal(TimeData.setDfRow!(df, df2[2, :], 2), df2)

## with DataFrameRow
##------------------

indCounter = 1
for row in eachrow(df2)
    TimeData.setDfRow!(df, row, indCounter)
    indCounter += 1
end
@test isequal(df, df2)

df2 = DataFrame()
df2[:a] = @data([4, NA])
df2[:b] = @data([3, 10])

indCounter = 1
for row in eachrow(df2)
    TimeData.setDfRow!(df, row, indCounter)
    indCounter += 1
end
@test isequal(df, df2)

## preallocateDf
##--------------

tm = TimeData.testcase(TimeData.Timematr, 1)
df = TimeData.preallocateDf(tm)
@test size(tm) == size(df)
@test names(tm) == names(df)

## preallocateBoolDf
##------------------

tm = TimeData.testcase(TimeData.Timematr, 1)
df = TimeData.preallocateBoolDf(names(tm), size(tm, 1))
@test size(tm) == size(df)
@test names(tm) == names(df)
@test all(eltypes(df) .== Bool)

end


text:
# sample_file.jl

# Compute some data
x = [1:10]
y = 2x

# Plot out some results
println("This is some output!")
println("Number of x values is ", length(x))

# Let's plot our results
using PGFPlots
p = Plots.Linear(x,y)
a = Axis(p)
save("sample_file_plot.pdf", a)


text:
module NahaGraphs

using Markdown

include("graphs.jl")
include("transform.jl")
include("generic_dot.jl")
include("dot_style.jl")

end


text:
# See https://udel.edu/~arthij/MD.pdf for information on forces
# See https://arxiv.org/pdf/1401.1181.pdf for applying forces to atoms
# See Gromacs manual for other aspects of forces

export
    force,
    accelerations,
    LennardJones,
    SoftSphere,
    Mie,
    Coulomb,
    Gravity,
    HarmonicBond,
    HarmonicAngle,
    Torsion,
    FinnisSinclair

"""
    force(inter, coord_i, coord_j, atom_i, atom_j, box_size)

Calculate the force between a pair of atoms due to a given interation type.
Custom interaction types should implement this function.
"""
function force end

# Allow 2D broadcasting whilst eliminating the diagonal corresponding to self interaction
@inline @inbounds function force(inter, coord_i, coord_j, atom_i, atom_j, box_size, self_interaction)
    if isone(self_interaction)
        return zero(coord_i)
    else
        return force(inter, coord_i, coord_j, atom_i, atom_j, box_size)
    end
end

@inline @inbounds function force!(forces, inter, s::Simulation, i::Integer, j::Integer)
    fdr = force(inter, s.coords[i], s.coords[j], s.atoms[i], s.atoms[j], s.box_size)
    forces[i] -= fdr
    forces[j] += fdr
    return nothing
end

"""
    accelerations(simulation; parallel=true)

Calculate the accelerations of all atoms using the general and specific
interactions and Newton's second law.
"""
function accelerations(s::Simulation; parallel::Bool=true)
    n_atoms = length(s.coords)

    if parallel && nthreads() > 1 && n_atoms >= 100
        forces_threads = [zero(s.coords) for i in 1:nthreads()]

        # Loop over interactions and calculate the acceleration due to each
        for inter in values(s.general_inters)

            if typeof(inter) <: GlueInteraction
                update_glue_densities!(inter, s.coords, s)
            end

            if inter.nl_only
                neighbors = s.neighbors
                @threads for ni in 1:length(neighbors)
                    i, j = neighbors[ni]
                    force!(forces_threads[threadid()], inter, s, i, j)
                end
            else
                @threads for i in 1:n_atoms
                    for j in 1:(i - 1)
                        force!(forces_threads[threadid()], inter, s, i, j)
                    end
                end
            end
        end

        forces = sum(forces_threads)
    else
        forces = zero(s.coords)

        for inter in values(s.general_inters)

            if typeof(inter) <: GlueInteraction
                update_glue_densities!(inter, s.coords, s)
            end

            if inter.nl_only
                neighbors = s.neighbors
                for ni in 1:length(neighbors)
                    i, j = neighbors[ni]
                    force!(forces, inter, s, i, j)
                end
            else
                for i in 1:n_atoms
                    for j in 1:(i - 1)
                        force!(forces, inter, s, i, j)
                    end
                end
            end
        end
    end

    for inter_list in values(s.specific_inter_lists)
        sparse_forces = force.(inter_list, (s.coords,), (s,))
        ge1, ge2 = getindex.(sparse_forces, 1), getindex.(sparse_forces, 2)
        sparse_vec = SparseVector(n_atoms, reduce(vcat, ge1), reduce(vcat, ge2))
        forces += Array(sparse_vec)
    end

    for i in 1:n_atoms
        forces[i] /= s.atoms[i].mass
        s.forces[i] = forces[i]
    end
    return forces
end

function accelerations(s::Simulation, coords, coords_is, coords_js, atoms_is, atoms_js, self_interactions)
    n_atoms = length(coords)
    forces = zero(coords)

    for inter in values(s.general_inters)
        # Currently the neighbor list is not used for this implementation
        forces -= reshape(sum(force.((inter,), coords_is, coords_js, atoms_is, atoms_js,
                                        s.box_size, self_interactions), dims=2), n_atoms)
    end

    for inter_list in values(s.specific_inter_lists)
        # Take coords off the GPU if they are on there
        coords_cpu = Array(coords)
        sparse_forces = force.(inter_list, (coords_cpu,), (s,))
        ge1, ge2 = getindex.(sparse_forces, 1), getindex.(sparse_forces, 2)
        sparse_vec = SparseVector(n_atoms, reduce(vcat, ge1), reduce(vcat, ge2))
        # Move back to GPU if required
        forces += convert(typeof(coords), Array(sparse_vec))
    end

    mass_i = findfirst(x -> x == :mass, fieldnames(eltype(atoms_is)))
    return forces ./ getfield.(s.atoms, mass_i)
end

include("interactions/lennard_jones.jl")
include("interactions/soft_sphere.jl")
include("interactions/mie.jl")
include("interactions/coulomb.jl")
include("interactions/gravity.jl")
include("interactions/harmonic_bond.jl")
include("interactions/harmonic_angle.jl")
include("interactions/torsion.jl")
include("interactions/glue_fs.jl")


text:
function matgen(n)
  tmp = 1.0 / n / n
  [ float64(tmp * (i - j) * (i + j - 2)) for i=1:n, j=1:n ]
end

function mul_line(n, linea, lineb)
  s = float64(0)
  for i=1:n
    s += linea[i] * lineb[i]
  end
  s
end

function mul(a, b)
  m = size(a)[1]
  n = size(a)[2]
  p = size(b)[2]

  b2 = [ b[j, i] for i=1:p, j=1:n ]
  [ mul_line(n, a[i,:], b2[j,:]) for i=1:m, j=1:p ]
end

function main()
  n = 100
  if length(ARGS) >= 1
    n = int(ARGS[1])
  end
  t = time()
  n = int(n / 2 * 2)
  a = matgen(n)
  b = matgen(n)
  c = mul(a, b)
  v = int(n/2) + 1
  println(c[v, v])
  println(time() - t)
end

main()

max_stars_repo_name: miloyip/benchmarks

max_issues_repo_name: miloyip/benchmarks

max_forks_repo_name: miloyip/benchmarks


text:
# ============================================================================ #
# [Variable length integer]                                                    #
# ============================================================================ #

primitive type VarUInt <: Unsigned 64 end

VarUInt(x::Number) = reinterpret(VarUInt, UInt64(x))
Base.UInt64(x::VarUInt) = reinterpret(UInt64, x)
Base.show(io::IO, x::VarUInt) = print(io, UInt64(x))


text:
export StaticLineSearch, Backtracking

# ================================================

"""
`type StaticLineSearch`: does nothing, i.e., just returns the step
proposal.
"""
type StaticLineSearch end

linesearch!(ls::StaticLineSearch, F, f0, df0, x, p, α; f0_goal=f0) = α, 0, nothing

"""
`type Backtracking`

kw-parameters:
* `c1`: Armijo constant
* `order`: order of interpolation during backtracking; currently only 2 is
   implemented
* `minα`: smallest allowed step (throws exception otherwise)
* `mindecfact`: minimal decrease factor
"""
@with_kw type Backtracking
   c1::Float64 = 0.1
   order::Int = 2
   minα::Float64 = 1e-8
   mindecfact::Float64 = 0.25
end

function linesearch!(ls::Backtracking, F, f0, df0, x, p, α; f0_goal=f0, condition=iter->false)
   # read parameters, initialise variables
   @unpack c1, order, minα, mindecfact = ls
   @assert f0_goal >= f0
   fα = F(x + α * p)
   numE = 1
   # start backtracking loop
   while (fα > f0_goal + c1 * α * df0) && !condition(numE)
      # check which form of backtracking to use
      if order == 2
         α1 = - 0.5 * (df0 * α) / ((fα - f0)/α - df0)
         α = max(α1, α * mindecfact)     # avoid miniscule steps
      else
         error("currently only quadratic backtracking is implemented")
      end
      # check whether α is too small now; probably something has gone wrong
      if α < minα
         return NaN, numE, fα
      end
      # evaluate f again and loop
      fα = F(x + α * p)
      numE += 1
   end
   return α, numE, fα
end


text:
struct PimcParameters
    steps_equil::Int
    steps_production::Int
    step_size_1::StepSize
    step_size_j::StepSize
end

function adjustment_factor(acceptance_ratio)
    if acceptance_ratio < 0.2
        0.9
    elseif 0.2 <= acceptance_ratio < 0.3
        0.95
    elseif 0.3 <= acceptance_ratio < 0.4
        0.99
    elseif 0.4 == acceptance_ratio
        1.0
    elseif 0.4 < acceptance_ratio <= 0.5
        1.01
    elseif 0.5 < acceptance_ratio <= 0.6
        1.05
    elseif 0.6 < acceptance_ratio
        1.1
    end
end

function pimc(sys::PimcSystem, ip::PimcParameters, st::State{P}) where {P}
    sss = [ip.step_size_1; [ip.step_size_j for _ in 2:P]]

    nums_total = [0 for _ in 1:P]
    nums_accept = [0 for _ in 1:P]

    for i in 1:ip.steps_equil
        if i % 100 == 0
            for j in 1:P
                sss[j] *= adjustment_factor(nums_accept[j] / nums_total[j])
                nums_total[j] = 0
                nums_accept[j] = 0
            end
        end

        accept, st = step_1(sys, st, sss[1])
        nums_total[1] += 1
        accept && (nums_accept[1] += 1)

        for j in 2:P
            accept, st = step_j(sys, st, j, sss[j])
            nums_total[j] += 1
            accept && (nums_accept[j] += 1)
        end
    end

    num_total_1 = 0
    num_total_j = 0
    num_accept_1 = 0
    num_accept_j = 0

    derivs1 = zeros(ip.steps_production)
    derivs2 = zeros(ip.steps_production)
    dist = zeros(P, ip.steps_production)

    for i in 1:ip.steps_production
        accept, st = step_1(sys, st, sss[1])
        num_total_1 += 1
        accept && (num_accept_1 += 1)

        for j in 2:P
            accept, st = step_j(sys, st, j, sss[j])
            num_total_j += 1
            accept && (num_accept_j += 1)
        end

        dlogJ, chvar = deriv_estimator_pre(st)

        # Free energy derivative estimator 1.
        derivs1[i] += dlogJ

        for j in 1:P
            derivs1[i] += deriv_estimator_1(sys, chvar, st, j) * sys.beta / P
        end

        # Free energy derivative estimator 2.
        derivs2[i] += dlogJ
        derivs2[i] += deriv_estimator_1(sys, chvar, st, 1) * sys.beta / P
        derivs2[i] -= deriv_estimator_2(chvar, st) * sys.mass * P / (hbar^2 * sys.beta)

        dist[:, i] .= path_xi(st)
    end

    hists = NTuple{2,Vector{Float64}}[]

    for j in 1:P
        h = normalize(fit(Histogram, dist[j, :], nbins=50))
        edges = only(h.edges)
        mids = (edges[begin:(end-1)] .+ edges[(begin+1):end]) ./ 2
        push!(hists, (mids, h.weights))
    end

    (aggregate(derivs1), aggregate(derivs2), num_accept_1 / num_total_1,
        num_accept_j / num_total_j, hists)
end


text:
# Discrete univariate

convert(::Type{Binomial}, d::Bernoulli) = Binomial(1, d.p)

# Continuous univariate

convert(::Type{Gamma}, d::Exponential) = Gamma(1.0, d.θ)
convert(::Type{Gamma}, d::Erlang) = Gamma(d.α, d.θ)

convert(::Type{GeneralizedGamma}, d::Gamma) = GeneralizedGamma(d.θ, d.α, 1.0)
convert(::Type{GeneralizedGamma}, d::Exponential) = GeneralizedGamma(d.θ, 1.0, 1.0)
convert(::Type{GeneralizedGamma}, d::Weibull) = GeneralizedGamma(d.θ, d.α, d.α)


text:
"Dirichlet Process Mixture Model Hyperparameters"
immutable DPMHyperparam <: AbstractHyperparam

  γ_a::Float64
  γ_b::Float64

  "default values"
  DPMHyperparam() = new(1.0, 1.0)

end

"Dirichlet Process Mixture Model Data Object"
type DPMData <: AbstractModelData

  # Energy
  energy::Float64

  # Dirichlet concentration parameter
  α::Float64

  # Distributions
  distributions::Array{ConjugatePostDistribution}

  # Assignments
  assignments::Array{Int}

  # Weights
  weights::Array{Float64}

end

type DPMBuffer <: AbstractModelBuffer

  # ids used for random access
  ids::Array{Int}

  # dimensionality of data
  D::Int

  # number of samples
  N::Int

  # samples
  X::AbstractArray

  # assignments
  Z::Array{Int}

  # number of samples per cluster
  C::Array{Int}

  # number of active cluster
  K::Int

  # distributions
  G::Array{ConjugatePostDistribution}

  # base distribution
  G0::ConjugatePostDistribution

  # concentration parameter
  alpha::Float64
end

@doc doc"""
Initialization of dirichlet process mixture model using K-Means.
""" ->
function init_kmeans_dpmm{T <: Real}(X::AbstractMatrix{T}, G0::ConjugatePostDistribution; k = 2, maxiterations = 10)

  (N, D) = size(X)

	if issparse(X)
    if !(T == Float64)
		  R = Clustering.kmeans(float(full(X')), k; maxiter = maxiterations)
    else
      R = Clustering.kmeans(full(X'), k; maxiter = maxiterations)
    end
	else
		R = Clustering.kmeans(X', k; maxiter = maxiterations)
	end

  Z = assignments(R)

  G = Array{ConjugatePostDistribution}(N)

  for c in 1:N
    idx = find(Z .== c)

    if length(idx) > 0
      G[c] = add_data(G0, X[idx,:])
    else
      G[c] = deepcopy(G0)
    end
  end

  return (Z, G)
end

@doc doc"""
Initialization of dirichlet process mixture model using random assignments.
""" ->
function init_random_dpmm(X::AbstractArray, G0::ConjugatePostDistribution; k = 2)

  (N, D) = size(X)

  Z = 1 + (randperm(N) % k)

  G = Array{ConjugatePostDistribution}(N)

  for c in 1:N
    idx = find(Z .== c)

    if length(idx) > 0
      G[c] = add_data(G0, X[idx,:])
    else
      G[c] = deepcopy(G0)
    end
  end

  return (Z, G)
end

@doc doc"""
Initialization of dirichlet process mixture model using precomputed assignments.
""" ->
function init_precomputed_dpmm(X::AbstractArray, G0::ConjugatePostDistribution, Z::Array{Int})

  (N, D) = size(X)

  G = Array{ConjugatePostDistribution}(N)

  for c in 1:N
    idx = find(Z .== c)

    if length(idx) > 0
      G[c] = add_data(G0, X[idx,:])
    else
      G[c] = deepcopy(G0)
    end
  end

  return (Z, G)
end

"Single iteration of collabsed Gibbs sampling using CRP."
function cgibbs_crp_dpmm!(B::DPMBuffer)

  # randomize data
  shuffle!(B.ids)

  z = -1
  k = -1

  for index in B.ids

    x = @view B.X[index, :]

    # get assignment
    z = B.Z[index]

    # remove sample from cluster
    remove_data!(B.G[z], x)

    # remove cluster assignment
    B.C[z] -= 1

    # udpate number of active clusters if necessary
    if B.C[z] < 1
      B.K -= 1
    end

    # remove cluster assignment
    z = -1

    # compute posterior predictive
    # compute priors using chinese restaurant process
    # see: Samuel J. Gershman and David M. Blei, A tutorial on Bayesian nonparametric models.
    # In Journal of Mathematical Psychology (2012)
    p = ones(B.K + 1) * -Inf
    k2id = ones(Int, B.K + 1)

    j = 1
    for i in 1:length(B.C)
      if B.C[i] >= 1

        llh = logpred( B.G[i], x )[1]
        crp = log( B.C[i] / (B.N + B.alpha - 1) )

        p[j] = llh + crp
        k2id[j] = i

        j += 1

      end
    end

    k2id[B.K + 1] = k2id[B.K] + 1

    p[B.K + 1] = logpred(B.G0, x)[1] + log( B.alpha / (B.N + B.alpha - 1) )
    p = exp(p - maximum(p))

    k = k2id[rand_indices(p)]

    if k > length(B.G)
      # add new cluster
      Gk = add_data(B.G0, x)
      B.G = cat(1, B.G, Gk)
      B.C = cat(1, B.C, 0)
    else
      # add to cluster
      add_data!(B.G[k], x)
    end

    if k > B.K
      B.K += 1
    end

    B.C[k] += 1
    B.Z[index] = k
  end

  B
end

"Compute Energy of model for given data"
function compute_energy!(B::DPMData, X::AbstractArray)

  E = 0.00001

  for xi in 1:size(X, 1)

    pp = 0.0
    c = 0

    for i = 1:length(B.weights)
      p = exp( logpred( B.distributions[i], X[xi, :] )[1] ) * B.weights[i]

      # only sum over actual values (excluding nans)
      if p == p
        pp += p
        c += 1
      end
    end

    E += pp

  end

  B.energy = log( E / size(X, 1) )

end

"Training Dirichlet Process Mixture Model using collabsed Gibbs sampling"
function train_cgibbs_dpmm(X::AbstractArray, G0::ConjugatePostDistribution, Z::Array{Int}, G::Array{ConjugatePostDistribution}, hyper::DPMHyperparam; alpha = 1.0, burnin = 0, thinout = 1, maxiter = 100)

  (N, D) = size(X)

  results = DPMData[]

  # construct C
  zz = unique(Z)

  C = zeros(Int, N)
  for i = 1:N
    C[i] = sum(Z .== i)
  end

  K = length(zz)

  # construct Buffer
  B = DPMBuffer(
    collect(1:N), # ids
    D,
    N,
    X,
    Z,
    C,
    K,
    G,
    G0,
    alpha)

  # burn-in phase
  for iter in 1:burnin
    # run one gibbs iteration using chinese restaurant process
    B = cgibbs_crp_dpmm!(B)

    # sample concentration parameter from G(1, 1)
    B.alpha = random_concentration_parameter(B.alpha, hyper.γ_a, hyper.γ_b, B.N, B.K)
  end

  # Gibbs sweeps

  for iter in 1:maxiter

    #println("   # [DPMM]: Iteration $(iter) of $(maxiter)")

    for t in 1:thinout
      # run one gibbs iteration using chinese restaurant process
      B = cgibbs_crp_dpmm!(B)

      # sample concentration parameter from G(1, 1)
      B.alpha = random_concentration_parameter(B.alpha, hyper.γ_a, hyper.γ_b, B.N, B.K)
    end

    # record model
    model = DPMData(0.0, deepcopy(B.alpha), deepcopy(B.G), deepcopy(B.Z), map(i -> B.C[i] / (B.N + B.alpha - 1), 1:B.K))

    # compute energy
    compute_energy!(model, X)

    # record results
    push!(results, model)
  end

  return results

end


text:
# cubic is default

struct Linear end
struct Cubic end

struct CompletePolnomial{order} end
struct Smolyak end

# incomplete type inference:
DecisionRule(exo_grid, endo_grid, i::Int64) = DecisionRule(exo_grid, endo_grid, Val{i})
DecisionRule(exo_grid, endo_grid, i::Int64, interp_type) = DecisionRule(exo_grid, endo_grid, Val{i}, interp_type)

# cubic interpolation
function DecisionRule(exo_grid::EmptyGrid, endo_grid::CartesianGrid, ::Type{Val{nx}}, ::Type{Cubic}) where nx
    CubicDR(exo_grid, endo_grid, Val{nx})
end

function DecisionRule(exo_grid::UnstructuredGrid, endo_grid::CartesianGrid, ::Type{Val{nx}}, ::Type{Cubic}) where nx
    CubicDR(exo_grid, endo_grid, Val{nx})
end

function DecisionRule(exo_grid::CartesianGrid, endo_grid::CartesianGrid, ::Type{Val{nx}}, ::Type{Cubic}) where nx
    CubicDR(exo_grid, endo_grid, Val{nx})
end

# Smolyak interpolation
function DecisionRule(
        exo_grid::S, endo_grid::SmolyakGrid, ::Type{Val{nx}}, ::Type{Smolyak}
    ) where S <: Union{EmptyGrid,<:UnstructuredGrid} where nx
    SmolyakDR(exo_grid, endo_grid, Val{nx})
end

# Complete polynomials
function DecisionRule(
        exo_grid::S, endo_grid::Grid, ::Type{Val{nx}}, ::Type{CompletePolnomial{order}}
    ) where S <: Union{EmptyGrid,<:UnstructuredGrid} where nx where order
    CompletePolyDR(exo_grid, endo_grid, Val{nx}, order)
end

# default to order 3
function DecisionRule(
        exo_grid::S, endo_grid::Grid, ::Type{Val{nx}}, ::Type{<:CompletePolnomial}
    ) where S <: Union{EmptyGrid,<:UnstructuredGrid} where nx
    CompletePolyDR(exo_grid, endo_grid, Val{nx})
end

# generic constructors

function DecisionRule(exo_grid::Grid, endo_grid::Grid, vals::Vector{ListOfPoints{n_x}}) where n_x
    dr = DecisionRule(exo_grid, endo_grid, Val{n_x}, Cubic)
    set_values!(dr, vals)
    dr
end

function DecisionRule(exo_grid::Grid, endo_grid::Grid, vals::Vector{ListOfPoints{n_x}}, tt) where n_x
    dr = DecisionRule(exo_grid, endo_grid, Val{n_x}, tt)
    set_values!(dr, vals)
    dr
end

# #####
# ##### Cached Decision Rules (do we really need them ?)
# #####

mutable struct CachedDecisionRule{T,S}
    dr::T
    process::S
end

CachedDecisionRule(process::AbstractDiscretizedProcess, grid::Grid, n_x::Int) =
    CachedDecisionRule(DecisionRule(process.grid, grid, n_x), process)

CachedDecisionRule(process::AbstractDiscretizedProcess, grid::Grid, values) =
    CachedDecisionRule(DecisionRule(process.grid, grid, values), process)

set_values!(cdr::CachedDecisionRule, v) = set_values!(cdr.dr, v)

# defaults

(cdr::CachedDecisionRule)(v::Union{<:Point,<:ListOfPoints}) = evaluate(cdr.dr,v)
(cdr::CachedDecisionRule)(m::Union{<:Point,<:ListOfPoints}, s::Union{<:Point,<:ListOfPoints}) = evaluate(cdr.dr, m, s)

# Interpolating rules
(cdr::CachedDecisionRule{<:AbstractDecisionRule{T,<:Grid}, <:AbstractDiscretizedProcess})(i::Int, s::Union{<:Point,<:ListOfPoints}) where T<:Union{CartesianGrid, RandomGrid} = evaluate(cdr.dr,node(Point, cdr.process, i), s)

(cdr::CachedDecisionRule{<:AbstractDecisionRule{T,<:Grid}, <:AbstractDiscretizedProcess})(i::Int, j::Int, s::Union{<:Point,<:ListOfPoints}) where T<:Union{CartesianGrid, RandomGrid} = evaluate(cdr.dr,inode(Point, cdr.process, i, j), s)

(cdr::CachedDecisionRule{<:AbstractDecisionRule{<:EmptyGrid,<:Grid}, <:AbstractDiscretizedProcess})(i::Int, j::Int, s::Union{<:Point,<:ListOfPoints}) = evaluate(cdr.dr, s)

(cdr::CachedDecisionRule{<:AbstractDecisionRule{<:EmptyGrid,<:Grid}, <:AbstractDiscretizedProcess})(i::Int, s::Union{<:Point,<:ListOfPoints}) = evaluate(cdr.dr, s)

# maybe keep only the i,j,s calls.

(cdr::CachedDecisionRule{<:AbstractDecisionRule{T,<:Grid}, <:DiscreteMarkovProcess})(i::Int, s::Union{<:Point,<:ListOfPoints}) where T<:Union{CartesianGrid, RandomGrid} = evaluate(cdr.dr,node(Point, cdr.process, i), s)
(cdr::CachedDecisionRule{<:AbstractDecisionRule{T,<:Grid}, <:DiscreteMarkovProcess})(i::Int, j::Int, s::Union{<:Point,<:ListOfPoints}) where T<:Union{CartesianGrid, RandomGrid} = evaluate(cdr.dr,inode(Point, cdr.process, i, j), s)

(cdr::CachedDecisionRule{<:AbstractDecisionRule{<:UnstructuredGrid,<:Grid}, <:DiscreteMarkovProcess})(i::Int, s::Union{<:Point,<:ListOfPoints}) = evaluate(cdr.dr, i, s)
(cdr::CachedDecisionRule{<:AbstractDecisionRule{<:UnstructuredGrid,<:Grid}, <:DiscreteMarkovProcess})(i::Int, j::Int, s::Union{<:Point,<:ListOfPoints}) = evaluate(cdr.dr, j, s)

#
# # We don't do compat for CachedDR cuz it's a nigthmare.
#
# (cdr::CachedDecisionRule{<:AbstractDecisionRule{<:UnstructuredGrid,<:Grid}, <:DiscreteMarkovProcess})(i::Int, s::Union{Vector{Float64},Matrix{Float64}}) = evaluate(cdr.dr, i, s)
# (cdr::CachedDecisionRule{<:AbstractDecisionRule{<:UnstructuredGrid,<:Grid}, <:DiscreteMarkovProcess})(i::Int, j::Int, s::Union{Vector{Float64},Matrix{Float64}}) = evaluate(cdr.dr, j, s)
#
#
# (cdr::CachedDecisionRule)(v::Union{Vector{Float64},Matrix{Float64}}) = evaluate(cdr.dr,v)
# (cdr::CachedDecisionRule)(m::Union{Vector{Float64},Matrix{Float64}}, s::Union{Vector{Float64},Matrix{Float64}}) = evaluate(cdr.dr, m, s)
# (cdr::CachedDecisionRule)(i::Int, s::Union{Vector{Float64},Matrix{Float64}}) = evaluate(cdr.dr,node(cdr.process, i), s)
# (cdr::CachedDecisionRule)(i::Int, j::Int, s::Union{Vector{Float64},Matrix{Float64}}) = evaluate(cdr.dr,inode(cdr.process, i, j), s)


text:
"""
    mutable struct Linear <: Block

Applies a linear transformation to the incoming data: y = w * x .+ b
"""
mutable struct Linear <: Block
    w::VarOrNil # input to hidden weights
    b::VarOrNil # bias of hidden units
    function Linear(inputSize::Int, hiddenSize::Int; type::Type=Array{Float32})
        T = eltype(type)
        A = sqrt(T(1/hiddenSize))
        w = randn(T, hiddenSize, inputSize) .* A
        b = randn(T, hiddenSize,         1) .* A
        new(Variable{type}(w,true,true,true),
            Variable{type}(b,true,true,true))
    end
    function Linear()
        new(nothing, nothing)
    end
end

function clone(this::Linear; type::Type=Array{Float32})
    cloned = Linear()
    cloned.w = clone(this.w, type=type)
    cloned.b = clone(this.b, type=type)
    return cloned
end

# pretty show
function Base.show(io::IO, m::Linear)
    SIZE = size(m.w)
    TYPE = typeof(m.w.value)
    print(io, "Linear($(SIZE[2]), $(SIZE[1]); type=$TYPE)")
end

"""
    unbiasedof(m::Linear)

unbiased weights of `Linear` block
"""
function unbiasedof(m::Linear)
    weights = Vector(undef, 1)
    weights[1] = m.w.value
    return weights
end

function weightsof(m::Linear)
    weights = Vector(undef, 2)
    weights[1] = m.w.value
    weights[2] = m.b.value
    return weights
end

function gradsof(m::Linear)
    grads = Vector(undef, 2)
    grads[1] = m.w.delta
    grads[2] = m.b.delta
    return grads
end

function zerograds!(m::Linear)
    for v in gradsof(m)
        v .= 0.0
    end
end

function paramsof(m::Linear)
    params = Vector{Variable}(undef,2)
    params[1] = m.w
    params[2] = m.b
    return params
end

function xparamsof(m::Linear)
    xparams = Vector{XVariable}(undef,2)
    xparams[1] = ('w', m.w)
    xparams[2] = ('b', m.b)
    return xparams
end

function nparamsof(m::Linear)
    lw = length(m.w)
    lb = length(m.b)
    return (lw + lb)
end

function bytesof(model::Linear, unit::String="MB")
    n = nparamsof(model) * elsizeof(model.w)
    return blocksize(n, uppercase(unit))
end

function forward(m::Linear, x::Variable)
    w = m.w
    b = m.b
    return matAddVec(w * x, b)
end

function predict(m::Linear, x)
    w = m.w.value
    b = m.b.value
    return (w * x .+ b)
end

function to(type::Type, m::Linear)
    m.w = to(type, m.w)
    m.b = to(type, m.b)
    return m
end

function to!(type::Type, m::Linear)
    m.w = to(type, m.w)
    m.b = to(type, m.b)
    return nothing
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    logfactorial

"""
    logfactorial(n::Int) -> Float64

Return the approximate value of ln ``n!`` by Ramanujan method.
"""
function logfactorial(n::Int)
    n >= 0 || throw(DomainError(n, "n should be 0 or more"))
    n == 0 && return 0.0
    return n * log(n) - n + 1/6 * log(8*n^3 + 4*n^2 + n + 1/30) + 1/2 * log(pi)
end


text:
using Documenter, HTgnss

makedocs(;
    modules=[HTgnss],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/htyeim/HTgnss.jl/blob/{commit}{path}#L{line}",
    sitename="HTgnss.jl",
    authors="htyeim <htyeim@gmail.com>",
    assets=String[],
)

deploydocs(;
    repo="github.com/htyeim/HTgnss.jl",
)


text:
"""
    abstract type AbstractCorrespondence

Method for determining the correspondence matching between blobs and measurements.

Supports functions `assign(c::AbstractCorrespondence, blobs, coordinates)`, `too_far(c::AbstractCorrespondence, blob, coordinate)`

# Subtypes
- `HungarianCorrespondence`
- `NearestNeighborCorrespondence`
"""
abstract type AbstractCorrespondence end

dist_th(c::AbstractCorrespondence) = c.dist_th

"""
    HungarianCorrespondence <: AbstractCorrespondence

Use the Hungarian algorithm to assign measuements to blobs. Each measurement is assigned to one blob only.
Parameter `p > 0` influences how eager the assignement is on a spectrum between `p ≈ 0` corresponding to nearest neighbor matching, `p = 1` corresponding to minimizing the earth-movers distance, `p → ∞` corresponding to minimizing the maximum error. The default is `p = 1`.

# Parameters
- `p=1`: the exponent of the cost matrix
- `dist_th=2`: maximum allowed Mahalanobis distance between a blob and a measurement
"""
Base.@kwdef struct HungarianCorrespondence <: AbstractCorrespondence
    p::Float64 = 1
    dist_th::Float64 = 2
end

"""
    NearestNeighborCorrespondence <: AbstractCorrespondence

Assign each blob the measurement that appears closest to the blob. This method can assign the same measurement to multiple blobs.

# Parameters
- `dist_th=2`: maximum allowed Mahalanobis distance between a blob and a measurement
"""
Base.@kwdef struct NearestNeighborCorrespondence <: AbstractCorrespondence
    dist_th::Float64 = 2
end

"""
    MCCorrespondence <: AbstractCorrespondence

Assigns blobs to measurement by approximately integrating over the posterior distribution over blobs and performing the assignment using the inner assignment

# Parameters
- `inner`: inner assignment object
- `num_samples::Int = 20` number of Monte Carlo samples to draw. The inner assignment routine will be called this many times so it can get expensive to set this too high.
"""
Base.@kwdef struct MCCorrespondence{T<:AbstractCorrespondence} <: AbstractCorrespondence
    inner::T = HungarianCorrespondence()
    num_samples::Int = 20
end
dist_th(c::MCCorrespondence) = dist_th(c.inner)


text:
### Day 20
## Firewall rules
## Author: Thanasis Georgiou
workspace()

# Read ruleset file
f = open("rules.txt")
input = readlines(f)
inputLen = length(input)
close(f)

# Organize array
rules = Array(Array{UInt32}, inputLen)
for i = 1:inputLen
  rules[i] = map(str -> parse(UInt32, str), split(input[i], '-'))
end

# Sort rules
sort!(rules, lt = (a, b) -> a[1] < b[1])

# Lower valid IP
validIp = 0x1
for (lower, upper) in rules
  if validIp >= lower && validIp <= upper
    validIp = upper + 1
  end
end
println("Lower valid address: $validIp")

# Count allowed IPs
blocked = 0
count = 0
for (lower, upper) in rules
  if lower - 1 > blocked
    count += lower - blocked - 1
  end
  blocked = max(blocked, upper)
end
count += typemax(UInt32) - blocked

println("Allowed IPs: $count")


text:
"""
    mincost_flow(graph, node_demand, edge_capacity, edge_cost, solver,<keyword arguments>)

Find a flow satisfying the `node_demand` at each node and `edge_capacity` constraints for each edge
while minimizing the `sum(edge_cost.*flow)`.

-`node_demand` is a vector of nodal values, which should be positive for sources and
 negative at sink nodes, any node with zero demand is a transport node

- The problem can be seen as a linear programming problem and uses a LP
solver under the hood. We use Clp in the examples and tests.

Returns a flow matrix, flow[i,j] corresponds to the flow on the (i,j) arc.
# Arguments
- `edge_demand::AbstractMatrix`: demand a minimum flow for a given edge.
- `edge_demand_exact=true`: changes the capacity of a non zero edge to the demanded value
- `source_nodes::AbstractVector` Sources at which the nodal netflow is allowed to be greater than nodal demand **
- `sink_nodes::AbstractVector` Sinks at which the nodal netflow is allowed to be less than nodal demand **
   ** source_nodes & sink_nodes are only needed when nodal flow are not explictly set in node_demand

### Usage Example:

```julia
julia> using Clp: ClpSolver # use your favorite LP solver here
julia> g = lg.DiGraph(6) # Create a flow-graph
julia> lg.add_edge!(g, 5, 1)
julia> lg.add_edge!(g, 5, 2)
julia> lg.add_edge!(g, 3, 6)
julia> lg.add_edge!(g, 4, 6)
julia> lg.add_edge!(g, 1, 3)
julia> lg.add_edge!(g, 1, 4)
julia> lg.add_edge!(g, 2, 3)
julia> lg.add_edge!(g, 2, 4)
julia> w = zeros(6,6)
julia> w[1,3] = 10.
julia> w[1,4] = 5.
julia> w[2,3] = 2.
julia> w[2,4] = 2.
julia> demand = spzeros(6)
julia> demand[5] = 2
julia> demand[6] = -2
julia> capacity = ones(6,6)
julia> flow = mincost_flow(g, demand, capacity, w, ClpSolver())
```
"""

function mincost_flow(g::lg.DiGraph,
		node_demand::AbstractVector,
		edge_capacity::AbstractMatrix,
		edge_cost::AbstractMatrix,
		solver::AbstractMathProgSolver;
		edge_demand::Union{Nothing,AbstractMatrix} = nothing,
		edge_demand_exact::Bool = false, #If true changes capacity of that node to the value in edge demand
		source_nodes::AbstractVector{<:Integer} = Vector{Integer}(), #Source nodes at which to allow a netflow greater than nodal demand
		sink_nodes::AbstractVector{<:Integer} = Vector{Integer}()	 #Sink nodes at which to allow a netflow less than nodal demand
		)
	T = eltype(g)
	VT = promote_type(eltype(edge_capacity),eltype(node_demand),eltype(edge_cost))

	nv = lg.nv(g)
	ne = lg.ne(g)
	b = Vector{VT}(undef,nv)
	AI = Vector{T}(undef,2*ne)
	AJ = Vector{T}(undef,2*ne)
	AV = Vector{VT}(undef,2*ne)
	edge_costs = Vector{VT}(undef,ne)
	upperbounds = Vector{VT}(undef,ne)
	lowerbounds = Vector{VT}(undef,ne)
	sense = ['=' for _ in 1:nv]

	has_edge_demand= edge_demand !== nothing

	inds = [-1,0]
	for (n,e) in enumerate(lg.edges(g))
		inds .+= 2
		s = lg.src(e)
		t = lg.dst(e)
		AI[inds] = [s,t] #the node index
		AJ[inds] = [n,n]  #the edge index
		AV[inds] = [1,-1]
		upperbounds[n] = edge_capacity[s,t]
		lowerbounds[n] = has_edge_demand ? edge_demand[s,t] : zero(VT)
		if edge_demand_exact && has_edge_demand && lowerbounds[n] != 0
			upperbounds[n] = lowerbounds[n]
		end

		edge_costs[n] = edge_cost[s,t]
	end
	for n = 1:nv
		b[n] = node_demand[n]
	end

	for n in source_nodes
		sense[n] = '>'
	end
	for n in sink_nodes
		sense[n] = '<'
	end

	A = sparse(AI,AJ,AV,nv,ne)
	# return (edge_costs, A, sense, b, lowerbounds, upperbounds, solver)
	sol = linprog(edge_costs, A, sense, b, lowerbounds, upperbounds, solver)
	sol_sparse = sparse(view(AI,1:2:2*ne),view(AI,2:2:2*ne),sol.sol,nv,nv)
	return sol_sparse
end

@deprecate mincost_flow(g::lg.DiGraph,
		capacity::AbstractMatrix,
		demand::AbstractMatrix,
		cost::AbstractMatrix,
		solver::AbstractMathProgSolver,
		source::Int, # if source and/or sink omitted or not in nodes, circulation problem
		sink::Int)  mincost_flow(g,spzeros(lg.nv(g)),capacity,cost,solver,edge_demand=demand,source_nodes=[source],sink_nodes=[sink])

@deprecate mincost_flow(g::lg.DiGraph,
		capacity::AbstractMatrix,
		demand::AbstractMatrix,
		cost::AbstractMatrix,
		solver::AbstractMathProgSolver
		)  mincost_flow(g,spzeros(lg.nv(g)),capacity,cost,solver,edge_demand=demand)


text:
"""
`@map(d, x)`

Apply the expression `x` row by row in `d`: return the result as an array or as a table
(if the elements are `Tuples` or `NamedTuples`). Use `{}` syntax for automatically named `NamedTuples`.
Symbols refer to fields of the row.
In this context, `_` refers to the whole row. To use actual symbols, escape them with `^`, as in `^(:a)`.
Use `cols(c)` to refer to field `c` where `c` is a variable that evaluates to a symbol. `c` must be available in
the scope where the macro is called.

## Examples

```jldoctest map
julia> t = table((a = [1,2,3], b = ["x","y","z"]));

julia> @map t :b*string(:a)
3-element Array{String,1}:
 "x1"
 "y2"
 "z3"

julia> @map t {:a, copy = :a, :b}
Table with 3 rows, 3 columns:
a  copy  b
────────────
1  1     "x"
2  2     "y"
3  3     "z"
```
"""
macro map(args...)
    esc(map_helper(args...))
end

function map_helper(d, x)
    anon_func, syms = extract_anonymous_function(x, replace_field)
    if !isempty(syms) && !(:(_) in syms)
        fields = Expr(:call, :(JuliaDBMeta.distinct_tuple), syms...)
        :(map($anon_func, (JuliaDBMeta._table)($d), select = $fields))
    else
        :(map($anon_func, (JuliaDBMeta._table)($d)))
    end
end

function map_helper(x)
    i = gensym()
    Expr(:(->), i, map_helper(i, x))
end

replace_field(iter, x) =  Expr(:call, :getfield, iter, x)
replace_field(iter) = iter


text:
export stablepartition!, isstableset

function isstableset(S::Set{T}, E::Set{Tuple{T,T}}) where T
  for (i,j) in E
    if i in S && j in S
      return false
    end
  end
  return true
end

function stablepartition!(conflicts::Set{Tuple{T, T}},
                             stablesets::Array{Set{T}, 1} = Set{T}[])  where T

  # elements already in stablesets are already taken care of
  prepartitioned = union(stablesets...)

  # convert the set of conflicts into a dictionary
  # mapping each id to the ids of its neighbors in the graph
  id2neighborid = Dict{T, Set{T}}() # map from variable id to neighbor ids
  for e in conflicts
    for i in 1:2
      v = e[i]
      n = e[i%2+1]
      if v in prepartitioned
        continue
      elseif !(v in keys(id2neighborid))
        # initialize v and add neighbor n
        id2neighborid[v] = Set{T}([n])
      else
        # just add neighbor n
        union!(id2neighborid[v], n)
      end
    end
  end

  # apportion variables not prepartitioned into the stable sets,
  # creating additional stable sets if necessary
  # might be better to produce a set cover with maximal stable sets,
  # but that's (NP) hard
  stablepartition!(id2neighborid, stablesets)

  return stablesets
end

# id2neighborhoodid maps node ids to the list of the node's neighbors in the graph
# returns a list of stable sets in the graph whose union is the full set of nodes
# uses a greedy algorithm
function stablepartition!(id2neighborid::Dict{T, Set{T}},
                             stablesets::Array{Set{T}, 1} = Set{T}[])  where T
  for varid in keys(id2neighborid)
    foundhome = false
    varid in id2neighborid[varid] && error("cannot partition graph with self loop")
    for s in stablesets
      if length(intersect(s, id2neighborid[varid])) == 0
        # var has no neighbors in s, so add it to s
        union!(s, varid)
        foundhome = true
        break
      end
    end
    if !foundhome
      # var conflicts with all previous subsets, so make a new subset
      push!(stablesets, Set{T}([varid]))
    end
  end
  return stablesets
end


text:
using ITensors
using ITensorMPOCompression
using Test
using Revise

@testset "Blocking functions" begin

    is=Index(2,"Site,SpinHalf,n=1")
    V=ITensor(1.0,Index(3,"Link,l=0"),Index(3,"Link,qx"),is,is')

    W=ITensor(0.0,Index(4,"Link,l=0"),Index(4,"Link,l=1"),is,is')
    W=setV(W,V,matrix_state(lower,left))
    @test matrix(slice(W,is=>1,is'=>1)) ==
    [0.0 0.0 0.0 0.0;
     0.0 1.0 1.0 1.0;
     0.0 1.0 1.0 1.0;
     0.0 1.0 1.0 1.0]

    W=ITensor(0.0,Index(4,"Link,l=0"),Index(4,"Link,l=1"),is,is')
    W=setV(W,V,matrix_state(lower,right))
    @test matrix(slice(W,is=>1,is'=>1)) ==
    [1.0 1.0 1.0 0.0;
     1.0 1.0 1.0 0.0;
     1.0 1.0 1.0 0.0;
     0.0 0.0 0.0 0.0]

    lWlink=Index(4,"Link,l=1")
    L=ITensor(2.0,Index(3,"Link,ql"),Index(3,"Link,l=1"))
    Lplus,il=growRL(L,lWlink,V_offsets(1,1))
    @test matrix(Lplus) ==
    [1.0 0.0 0.0 0.0;
     0.0 2.0 2.0 2.0;
     0.0 2.0 2.0 2.0;
     0.0 2.0 2.0 2.0]

    #
    # these test are potentially flaky because they depend on the
    # index ordering of W.  In other words matrix(W) is ill defined, the transpose
    # of what you expect could be returned.
    #
    A=[i*1.0 for i in 1:4*4*2*2]
    V=ITensor(0.5,Index(3,"Link,l=0"),Index(3,"Link,qx"),is,is')
    W=ITensor(eltype(A),A,Index(4,"Link,l=0"),Index(4,"Link,l=1"),is,is')
    W=setV(W,V,matrix_state(lower,left))
    @test matrix(slice(W,is=>1,is'=>1)) ==
    [1.0 5.0 9.0 13.0 ;
     2.0 0.5 0.5  0.5 ;
     3.0 0.5 0.5  0.5 ;
     4.0 0.5 0.5  0.5 ]

    # Now force resizing.  It shoud preserve the last row and col of W
    V=ITensor(0.5,Index(3,"Link,l=0"),Index(2,"Link,qx"),is,is')
    W=ITensor(eltype(A),A,Index(4,"Link,l=0"),Index(4,"Link,l=1"),is,is')
    W=setV(W,V,matrix_state(lower,left))
    @test transpose(matrix(slice(W,is=>1,is'=>1))) ==
    [1.0 0.0  0.0 ; #top row gets zeroed out on resizing
     2.0 0.5  0.5 ;
     3.0 0.5  0.5 ;
     4.0 0.5  0.5 ]

    V=ITensor(0.5,Index(3,"Link,l=0"),Index(2,"Link,qx"),is,is')
    W=ITensor(eltype(A),A,Index(4,"Link,l=0"),Index(4,"Link,l=1"),is,is')

    W=setV(W,V,matrix_state(lower,right))

    @test matrix(slice(W,is=>1,is'=>1)) ==
    [0.5   0.5  0.5  0.0 ;
     0.5   0.5  0.5  0.0 ;
     13.0 14.0 15.0 16.0 ]

    V=ITensor(0.5,Index(2,"Link,qx"),Index(3,"Link,l=1"),is,is')

    W=ITensor(eltype(A),A,Index(4,"Link,l=0"),Index(4,"Link,l=1"),is,is')
    W=setV(W,V,matrix_state(lower,left))
    @test matrix(slice(W,is=>1,is'=>1)) ==
    [1.0 5.0 9.0 13.0;
     0.0 0.5 0.5  0.5;
     0.0 0.5 0.5  0.5]

    W=ITensor(eltype(A),A,Index(4,"Link,l=0"),Index(4,"Link,l=1"),is,is')

    W=setV(W,V,matrix_state(lower,right))
    @test matrix(slice(W,is=>1,is'=>1)) ==
    [0.5 0.5  0.5 0.0;
     0.5 0.5  0.5 0.0;
     4.0 8.0 12.0 16.0]

    #
    #  GrowRL tests
    #
    Lplus,il=growRL(L,lWlink,V_offsets(1,0))
    @test matrix(Lplus) ==
    [1.0 0.0 0.0 0.0;
     2.0 2.0 2.0 0.0;
     2.0 2.0 2.0 0.0;
     2.0 2.0 2.0 1.0]
    #
    Lplus,il=growRL(L,lWlink,V_offsets(0,1))
    @test matrix(Lplus) ==
    [1.0 2.0 2.0 2.0;
     0.0 2.0 2.0 2.0;
     0.0 2.0 2.0 2.0;
     0.0 0.0 0.0 1.0]
    #
    Lplus,il=growRL(L,lWlink,V_offsets(0,0))
    @test matrix(Lplus) ==
    [2.0 2.0 2.0 0.0;
     2.0 2.0 2.0 0.0;
     2.0 2.0 2.0 0.0;
     0.0 0.0 0.0 1.0]

end


text:
@defcomp climateco2cycle begin
    # Anthropogenic CO2 emissions in Mt of C
    mco2 = Parameter(index=[time])

    # Terrestrial biosphere CO2 emissions in Mt of C
    terrestrialco2 = Variable(index=[time])

    # Net CO2 emissions in Mt of C
    globc = Variable(index=[time])

    # Carbon boxes
    cbox = Variable(index=[time,5])

    # Initial carbon box 1
    cbox10 = Parameter()

    # Initial carbon box 2
    cbox20 = Parameter()

    # Initial carbon box 3
    cbox30 = Parameter()

    # Initial carbon box 4
    cbox40 = Parameter()

    # Initial carbon box 5
    cbox50 = Parameter()

    # Carbon decay in box 1
    co2decay1 = Variable()

    # Carbon decay in box 2
    co2decay2 = Variable()

    # Carbon decay in box 3
    co2decay3 = Variable()

    # Carbon decay in box 4
    co2decay4 = Variable()

    # Carbon decay in box 5
    co2decay5 = Variable()

    # Carbon decay in box 1
    lifeco1 = Parameter(default = 1.0)

    # Carbon decay in box 2
    lifeco2 = Parameter(default = 363.0)

    # Carbon decay in box 3
    lifeco3 = Parameter(default = 74.0)

    # Carbon decay in box 4
    lifeco4 = Parameter(default = 17.0)

    # Carbon decay in box 5
    lifeco5 = Parameter(default = 2.0)

    # Fraction of carbon emission in box 1
    co2frac1 = Parameter(default = 0.13)

    # Fraction of carbon emission in box 2
    co2frac2 = Parameter(default = 0.2)

    # Fraction of carbon emission in box 3
    co2frac3 = Parameter(default = 0.32)

    # Fraction of carbon emission in box 4
    co2frac4 = Parameter(default = 0.25)

    # Fraction of carbon emission in box 5
    co2frac5 = Parameter(default = 0.1)

    # Atmospheric CO2 concentration
    acco2 = Variable(index=[time])

    # Stock of CO2 in the terrestrial biosphere
    terrco2stock = Variable(index=[time])

    # Temperature
    temp = Parameter(index=[time])

    terrco2sens = Parameter(default = 2600.000000000002)
    terrco2stock0 = Parameter(default = 1.9e6)

    tempin2010 = Variable()

    function run_timestep(p, v, d, t)

        if is_first(t)
            v.co2decay1 = p.lifeco1
            v.co2decay2 = exp(-1.0 / p.lifeco2)
            v.co2decay3 = exp(-1.0 / p.lifeco3)
            v.co2decay4 = exp(-1.0 / p.lifeco4)
            v.co2decay5 = exp(-1.0 / p.lifeco5)

            v.terrco2stock[t] = p.terrco2stock0

            v.cbox[t,1] = p.cbox10
            v.cbox[t,2] = p.cbox20
            v.cbox[t,3] = p.cbox30
            v.cbox[t,4] = p.cbox40
            v.cbox[t,5] = p.cbox50
            v.acco2[t] = v.cbox[t,1] + v.cbox[t,2] + v.cbox[t,3] + v.cbox[t,4] + v.cbox[t,5]
        else

            if gettime(t) == 2011
                v.tempin2010 = p.temp[TimestepValue(2010)]
            end

            if gettime(t) > 2010
                v.terrestrialco2[t] = (p.temp[t - 1] - v.tempin2010) * p.terrco2sens * v.terrco2stock[t - 1] / p.terrco2stock0
            else
                v.terrestrialco2[t] = 0
            end

            v.terrco2stock[t] = max(v.terrco2stock[t - 1] - v.terrestrialco2[t], 0.0)

            v.globc[t] = p.mco2[t] + v.terrestrialco2[t]

            # Calculate CO2 concentrations
            v.cbox[t,1] = v.cbox[t - 1,1] * v.co2decay1 + 0.000471 * p.co2frac1 * (v.globc[t])
            v.cbox[t,2] = v.cbox[t - 1,2] * v.co2decay2 + 0.000471 * p.co2frac2 * (v.globc[t])
            v.cbox[t,3] = v.cbox[t - 1,3] * v.co2decay3 + 0.000471 * p.co2frac3 * (v.globc[t])
            v.cbox[t,4] = v.cbox[t - 1,4] * v.co2decay4 + 0.000471 * p.co2frac4 * (v.globc[t])
            v.cbox[t,5] = v.cbox[t - 1,5] * v.co2decay5 + 0.000471 * p.co2frac5 * (v.globc[t])

            v.acco2[t] = v.cbox[t,1] + v.cbox[t,2] + v.cbox[t,3] + v.cbox[t,4] + v.cbox[t,5]
        end
    end
end


text:
# ---
# title: 598. Range Addition II
# id: problem598
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/range-addition-ii/description/>
# hidden: true
# ---
#
# You are given an `m x n` matrix `M` initialized with all `0`'s and an array of
# operations `ops`, where `ops[i] = [ai, bi]` means `M[x][y]` should be
# incremented by one for all `0 <= x < ai` and `0 <= y < bi`.
#
# Count and return _the number of maximum integers in the matrix after
# performing all the operations_.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/02/ex1.jpg)
#
#
#
#     Input: m = 3, n = 3, ops = [[2,2],[3,3]]
#     Output: 4
#     Explanation: The maximum integer in M is 2, and there are four of it in M. So return 4.
#
#
# **Example 2:**
#
#
#
#     Input: m = 3, n = 3, ops = [[2,2],[3,3],[3,3],[3,3],[2,2],[3,3],[3,3],[3,3],[2,2],[3,3],[3,3],[3,3]]
#     Output: 4
#
#
# **Example 3:**
#
#
#
#     Input: m = 3, n = 3, ops = []
#     Output: 9
#
#
#
#
# **Constraints:**
#
#   * `1 <= m, n <= 4 * 104`
#   * `1 <= ops.length <= 104`
#   * `ops[i].length == 2`
#   * `1 <= ai <= m`
#   * `1 <= bi <= n`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
export update_H
"""
Update the value of the (approximation) of the second derivative depending on
the approximation we want.
If we want Newton it will compute the exact second derivative.
Secant or Improved Secant will compute approximation of the second derivative.
Inputs:
  - direction, string giving the second derivative to compute (either "Nwt",
  "Sec", "SecA")
  - h, a line model (our function)
  - t, our current point
  - φt, value of φ(t). For more information on φ see the doc for phi_dphi
  - dφt, value of φ'(t).
  - φtprec, previous value of φ(t).
  - dφtprec, previous value of φ'(t).
"""
function update_H(direction :: Symbol, h :: LineModel, t :: Float64,
                  tprec :: Float64, φt :: Float64, dφt :: Float64,
                  φtprec :: Float64, dφtprec :: Float64; kwargs...)

  if direction == :Nwt
    return hess(h, t)
  elseif direction == :Sec
    s    = t-tprec
    y    = dφt - dφtprec
    seck = y / s

    return seck
  elseif direction == :SecA

    s = t   - tprec
    y = dφt - dφtprec

    Γ= 3 * (dφt + dφtprec) * s - 6 * (φt - φtprec)
    if (y * s + Γ) < eps(Float64) * s ^ 2
      yt = y
    else
      yt =  y + Γ / s
    end

    seck = yt / s

    return seck
  end
end


text:
#xor-ml.jl - set up to test the xor function in julia.

using SigmoidNumbers
using GenML
using Base.Test

PType = Float32#Posit{12,0}

#assign sigmoid to the pseudologistic
GenML.TF.sigmoid(x::Posit) = SigmoidNumbers.pseudologistic(x)
GenML.TF.dasysigmoid(y::Posit) = SigmoidNumbers.delta_psl(y)
GenML.CF.logloss(x::Posit) = SigmoidNumbers.pseudologcost(x)

xornet = GenML.MLP.MultilayerPerceptron{PType, (2,2,1)}()

#hand-written transition matrices.
xornet.layers[1].transition[:] = PType[5.0 5.0; -10.0 -10.0]
xornet.layers[1].bias[:] = PType[-7.5, 7.5]
xornet.layers[2].transition[:] = PType[-10.0 -10.0]
xornet.layers[2].bias[:] = PType[5.0]

#test the possible input/output pairs.

println("attempting premade xor neural net, with posit $PType")

@test xornet([true, true])[1] < PType(0.5)
@test xornet([true, false])[1] > PType(0.5)
@test xornet([false, true])[1] > PType(0.5)
@test xornet([false, false])[1] < PType(0.5)

println("execution of premade xor neural net successful with $PType")

println("attempting to train xor neural net, with posit $PType")

Base.randn(::Type{PType}, dims::Integer...) = map(PType, randn(dims...))
Base.rand(::Type{PType}, dims::Integer...) = map(PType, rand(dims...))

function xortrain()

  println("working on unreliable xor data set with backpropagation")

  input_matrix = rand(Bool, 10, 500)
  training_results = Array{Bool,2}(1,500)
  training_results[:] = [input_matrix[1, col] $ input_matrix[2, col] $ (rand() < 0.05) for col in 1:size(input_matrix,2)]

  xornet = GenML.MLP.MultilayerPerceptron{PType,(10,2,1)}(randn)

  for rounds = 1:40
  for idx = 1:500
    input_column = input_matrix[:, idx]
    training_column = training_results[:, idx]

    GenML.Optimizers.backpropagationoptimize(xornet, input_column, training_column, GenML.CF.crossentropy)
  end
  end

  #verify that the optimization has resulted in a good data set.
  wrongcount = 0
  for x = 1:50
    input_vector = rand(Bool, 10)
    wrongcount += (xornet(input_vector)[1] > PType(0.5)) != (input_vector[1] $ input_vector[2])
  end
  println("incorrect responses for noisy xor: $wrongcount / 50 => $(2 * (50 - wrongcount)) % correct")

  return wrongcount
end

xortrain() == 0 || xortrain() == 0 || xortrain() == 0 || xortrain() == 0 || xortrain() == 0 || xortrain() == 0 || xortrain() == 0 || xortrain() == 0

max_issues_repo_name: MohHizzani/SigmoidNumbers

max_forks_repo_name: MohHizzani/SigmoidNumbers


text:
using churrosJulia

churro("{o}====} {======={o} {{==o}")

res = churro(str,diagnostics=false)

print("done: $res")