text:
using KitBase, FluxReconstruction, OrdinaryDiffEq, Langevin, LinearAlgebra, Plots, JLD2
using KitBase.ProgressMeter: @showprogress
using Base.Threads: @threads

begin
    x0 = 0
    x1 = 3
    ncell = 100#300
    nface = ncell + 1
    dx = (x1 - x0) / ncell
    deg = 3 # polynomial degree
    nsp = deg + 1
    cfl = 0.05
    dt = cfl * dx / 10
    t = 0.0

    uqMethod = "galerkin"
    nr = 9
    nRec = 18
    opType = "uniform"
    parameter1 = -1.0
    parameter2 = 1.0
end

ps = FRPSpace1D(x0, x1, ncell, deg)
uq = UQ1D(nr, nRec, parameter1, parameter2, opType, uqMethod)
VInv = inv(Array(ps.V))
l2 = [uq.t2Product[j-1, j-1] for j = 1:uq.nm+1]

cd(@__DIR__)
include("../../filter.jl")

function l2_strength(a1, a2=a1)
    (1/a1 - 1)/2/ps.deg^2 / (ps.deg+1)^2, (1/a2-1)/2/uq.nm^2 / (uq.nm+1)^2
end

u0 = zeros(ncell, nsp, uq.nm+1)
u0q = zeros(ncell, nsp, uq.nq)
let p0 = 0.5, p1 = 1.5, σ = 0.2, val0 = 11.0, val1 = 1.0
    for i = 1:ncell
        for j = 1:uq.nq
            if ps.x[i] < p0 + σ * uq.pceSample[j]
                u0q[i, :, j] .= val0
            elseif p0 + σ * uq.pceSample[j] <= ps.x[i] <= p1 + σ * uq.pceSample[j]
                u0q[i, :, j] .= val0 + (val1 - val0) / (p0 - p1) * (p0 + σ * uq.pceSample[j] - ps.x[i])
            else
                u0q[i, :, j] .= val1
            end
        end

        for j = 1:nsp
            u0[i, j, :] .= ran_chaos(u0q[i, j, :], uq)
        end
    end
end

function dudt!(du, u, p, t)
    du .= 0.0
    J, ll, lr, lpdm, dgl, dgr, Δx, uq = p

    nm = uq.nm
    nq = uq.nq

    ncell = size(u, 1)
    nsp = size(u, 2)

    u_ran = zeros(ncell, nsp, nq)
    @threads for j = 1:nsp
        for i = 1:ncell
            @inbounds u_ran[i, j, :] .= chaos_ran(u[i, j, :], uq)
        end
    end

    f = zeros(ncell, nsp, nm+1)
    @threads for j = 1:nsp
        for i = 1:ncell
            _f = zeros(nq)
            for k = 1:nq
                _f[k] = 0.5 * u_ran[i, j, k]^2 / J[i]
            end

            f[i, j, :] .= ran_chaos(_f, uq)
        end
    end

    u_face = zeros(ncell, nm+1, 2)
    f_face = zeros(ncell, nm+1, 2)
    @inbounds @threads for k = 1:nm+1
        for i = 1:ncell
            u_face[i, k, 1] = dot(u[i, :, k], lr)
            f_face[i, k, 1] = dot(f[i, :, k], lr)
            u_face[i, k, 2] = dot(u[i, :, k], ll)
            f_face[i, k, 2] = dot(f[i, :, k], ll)
        end
    end

    f_interaction = zeros(ncell + 1, nm+1)
    @threads for i = 2:ncell
        @inbounds f_interaction[i, :] .= (f_face[i-1, :, 1] .+ f_face[i, :, 2]) ./ 2 -
            (Δx[i-1] + Δx[i]) / 2 * (u_face[i, :, 2] .- u_face[i-1, :, 1])
    end
    #f_interaction[1, :] .= (f_face[ncell, :, 1] .+ f_face[1, :, 2]) ./ 2 -
    #    (Δx[ncell] + Δx[1]) / 2 * (u_face[1, :, 2] .- u_face[ncell, :, 1])
    #f_interaction[ncell+1, :] .= (f_face[ncell, :, 1] .+ f_face[1, :, 2]) ./ 2 -
    #    (Δx[ncell] + Δx[1]) / 2 * (u_face[1, :, 2] .- u_face[ncell, :, 1])

    rhs1 = zero(u)
    for i = 1:ncell, ppp1 = 1:nsp, l = 1:nm+1
        @inbounds rhs1[i, ppp1, l] = dot(f[i, :, l], lpdm[ppp1, :])
    end

    idx = 2:ncell-1
    @threads for l = 1:nm+1
        for i in idx, ppp1 = 1:nsp
            @inbounds du[i, ppp1, l] =
                -(
                    rhs1[i, ppp1, l] +
                    (f_interaction[i, l] - f_face[i, l, 2]) * dgl[ppp1] +
                    (f_interaction[i+1, l] - f_face[i, l, 1]) * dgr[ppp1]
                )
        end
    end
end

tspan = (0.0, 0.1)
p = (ps.J, ps.ll, ps.lr, ps.dl, ps.dhl, ps.dhr, ps.dx, uq)
prob = ODEProblem(dudt!, u0, tspan, p)
nt = tspan[2] ÷ dt |> Int
itg = init(prob, Tsit5(), saveat = tspan[2], adaptive = false, dt = dt)

@showprogress for iter = 1:nt
    step!(itg)

    # filter
    for i = 1:size(itg.u, 1)
        ũ = VInv * itg.u[i, :, :]
        #su = maximum([ũ[end, j]^2 / sum(ũ[:, j].^2) for j = 1:uq.nm+1])
        #sv = maximum([ũ[j, end]^2 / sum(ũ[j, :].^2) for j = 1:nsp])
        #isShock = max(detector(log10(su), ps.deg), detector(log10(sv), ps.deg))
        #=if isShock
            λ1 = dt * (su)
            λ2 = dt * (sv)

            #FR.filter_exp!(uModal, 10, 100)
            FR.modal_filter!(ũ, λ1, λ2; filter = :l2)
            #FR.modal_filter!(ũ, λ1, λ2; filter = :l2opt)
            #FR.modal_filter!(ũ; filter = :lasso)

            itg.u[i, :, :] .= ps.V * ũ
        end=#

        su = maximum([ũ[end, j]^2 / (sum(ũ[:, j].^2) + 1e-4) for j = 1:uq.nm+1])
        sv = maximum([ũ[j, end]^2 * uq.t2Product[uq.nm, uq.nm] / (sum(ũ[j, :].^2 .* l2) + 1e-4) for j = 1:nsp])
        isShock = max(
            #shock_detector(log10(su), ps.deg, -3 * log10(ps.deg), 1.0),
            shock_detector(log10(sv), ps.deg, -2 * log10(ps.deg), 1.0),
        )

        if isShock
            #FR.modal_filter!(ũ, 16; filter = :exp)

            #FR.modal_filter!(ũ, 8, 8; filter = :houli)

            #λ1 = dt * sqrt(su) * 2.0
            #λ2 = dt * sqrt(sv) * 2.0

            #λ1 = 3.5e-2
            #λ2 = λ1 * ps.deg^2 * (ps.deg+1)^2 / uq.nm^2 / (uq.nm+1)^2

            λ1, λ2 = l2_strength(0.6)

            #FR.modal_filter!(ũ, λ1, λ2; filter = :l2)
            #FR.modal_filter!(ũ, λ1, λ2; filter = :l2opt)
            #FR.modal_filter!(ũ, PhiL1; filter = :lasso)
            FR.modal_filter!(ũ, 3, 3, 36*dt; filter = :exp)

            itg.u[i, :, :] .= ps.V * ũ
        end

        #=tmp = [chaos_ran(itg.u[i, j, :], uq) for j = 1:nsp]
        uquad = [tmp[i][j] for i=1:nsp, j=1:uq.nq]
        δu = [maximum(uquad[:, j]) - minimum(uquad[:, j]) for j = 1:uq.nm+1] |> maximum
        δv = [maximum(uquad[j, :]) - minimum(uquad[j, :]) for j = 1:nsp] |> maximum
        λ1 = dt * relu(δu - 7.0)
        λ2 = dt * relu(δv - 7.0)
        uModal = VInv * itg.u[i, :, :]
        FR.modal_filter!(uModal, λ1, λ2; filter = :l2)
        itg.u[i, :, :] .= ps.V * uModal=#
    end
end

begin
    x = zeros(ncell * nsp)
    w = zeros(ncell * nsp, uq.nm+1)
    for i = 1:ncell
        idx0 = (i - 1) * nsp

        for j = 1:nsp
            idx = idx0 + j
            x[idx] = ps.xpg[i, j]

            w[idx, :] .= itg.u[i, j, :]
        end
    end

    sol = zeros(ncell*nsp, 2)
    for i in axes(sol, 1)
        sol[i, 1] = mean(w[i, :], uq.op)
        sol[i, 2] = std(w[i, :], uq.op)
    end

    pic1 = plot(x, sol[:, 1], label="u", xlabel="x", ylabel="mean")
    pic2 = plot(x, sol[:, 2], label="u", xlabel="x", ylabel="std")
    plot(pic1, pic2)
end

plot(x, sol0[:, 2], label="Lasso", line=:dash)
plot!(x, sol[:, 2], label="L2", xlabel="x", ylabel="std")

uξ = chaos_ran(itg.u[53, 2, :], uq)
plot(uq.op.quad.nodes, uξ)

sol0 = deepcopy(sol)
#@save "nofilter.jld2" x sol
#@save "l2_apt.jld2" x sol
#@save "lasso_apt.jld2" x sol
#@save "l2opt_apt.jld2" x sol
#@save "exp.jld2" x sol
#@save "houli.jld2" x sol
#@save "lasso.jld2" x sol
#@save "l2.jld2" x sol
#@save "l2_apt.jld2" x sol

max_stars_repo_name: CSMMLab/FRSG

max_issues_repo_name: CSMMLab/FRSG

max_forks_repo_name: CSMMLab/FRSG


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible ideal GLM-MHD equations

equations = IdealGlmMhdEquations2D(1.4)

initial_condition = initial_condition_weak_blast_wave

surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell)
volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
dg = DGMulti(polydeg=3, element_type = Quad(), approximation_type = Polynomial(),
             surface_integral = SurfaceIntegralWeakForm(surface_flux),
             volume_integral = VolumeIntegralFluxDifferencing(volume_flux))

mesh = DGMultiMesh(dg, cells_per_dimension=(16, 16),
                   coordinates_min=(-2.0, -2.0), coordinates_max=(2.0, 2.0),
                   periodicity=true)
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 0.4)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
alive_callback = AliveCallback(analysis_interval=analysis_interval)

cfl = 1.0
stepsize_callback = StepsizeCallback(cfl=cfl)

glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl)

callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        stepsize_callback,
                        alive_callback,
                        glm_speed_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);

summary_callback() # print the timer summary

ext: jl


text:
using SafeTestsets

@safetestset "Tool tests" begin include("Test_costs.jl") end


text:
using FactCheck

include("lasso.jl")
include("fusedlasso.jl")
include("trendfiltering.jl")

FactCheck.exitstatus()


text:
# SPDX-License-Identifier: X11
# 2020-11-19
# Vacation (100pt)

function main()
  N = parse(Int, readline())
  vals = zeros(Int, N + 1, 3)

  for i ∈ 2:N + 1
    A,B,C = parse.(Int, (split ∘ readline)())

    vals[i, 1] = A + max(vals[i - 1, 2], vals[i - 1, 3])
    vals[i, 2] = B + max(vals[i - 1, 1], vals[i - 1, 3])
    vals[i, 3] = C + max(vals[i - 1, 1], vals[i - 1, 2])
  end

  println(max(vals[N + 1, 1], vals[N + 1, 2], vals[N + 1, 3]))
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
# heavily influenced by and mostly of the Arb C interface used in Nemo.jl
#=
    There is no direct support for matricies of ArbFloats.
    We partially overload ArbRealMatrix to zero the radii.
    And provide conversion from Array{ArbFloat{P},2}.
=#

mutable struct ArbFloatMatrix{P} <: AbstractArbMatrix{P, ArbFloat}
    arbrealmatrix::ArbRealMatrix{P}

    function ArbFloatMatrix{P}(rowcount::Int, colcount::Int) where {P}
        z = ArbRealMatrix{P}(rowcount, colcount)
        return new{P}(z)
    end
end

ArbFloatMatrix(x::ArbFloatMatrix{P}) where {P} = x
ArbFloatMatrix{P}(x::ArbFloatMatrix{P}) where {P} = x

@inline function ArbFloatMatrix(nrows::Int, ncols::Int)
    P = workingprecision(ArbFloat)
	return ArbFloatMatrix{P}(nrows, ncols)
end

function ArbFloatMatrix(fpm::Array{ArbFloat, 2})
    P = workingprecision(ArbFloat)
    nrows, ncols = size(fpm)
    arm = ArbFloatMatrix{P}(nrows, ncols)
    for r = 1:nrows
       for c = 1:ncols
	       arm[r,c] = fpm[r,c]
	   end
    end
    return arm
end

function ArbFloatMatrix(fpm::Array{ArbFloat{P},2}) where {P}
    nrows, ncols = size(fpm)
    arm = ArbFloatMatrix{P}(nrows, ncols)
    for r = 1:nrows
	    for c = 1:ncols
			arm.arbrealmatrix[r,c] = ArbReal{P}(fpm[r,c])
	    end
	end
	return arm
end

function ArbFloatMatrix{P}(fpm::Array{ArbFloat{P},2}) where {P}
    nrows, ncols = size(fpm)
    arm = ArbFloatMatrix{P}(nrows, ncols)
    for r = 1:nrows
		for c = 1:ncols
			arm.arbrealmatrix[r,c] = ArbReal{P}(fpm[r,c])
	    end
	end
	return arm
end

function ArbFloatMatrix{P}(fpm::Array{ArbFloat,2}) where {P}
    nrows, ncols = size(fpm)
    arm = ArbFloatMatrix{P}(nrows, ncols)
    for r = 1:nrows
		for c = 1:ncols
			arm.arbrealmatrix[r,c] = ArbReal{P}(fpm[r,c])
	    end
	end
	return arm
end

function ArbFloatMatrix(fpm::Array{F,2}) where {F<:Union{Integer,AbstractFloat}}
    P = workingprecision(ArbFloat)
    nrows, ncols = size(fpm)
    arm = ArbFloatMatrix{P}(nrows, ncols)
    for r = 1:nrows
		for c = 1:ncols
			arm.arbrealmatrix[r,c] = ArbReal{P}(fpm[r,c])
	    end
	end
	return arm
end

function Matrix(arm::ArbFloatMatrix{P}) where {P}
    nrows, ncols = rowcount(arm), colcount(arm)
    m = zeros(ArbFloat{P}, nrows, ncols)
    for r = 1:nrows
        for c = 1:ncols
            m[r,c] = ArbFloat{P}(arm.arbrealmatrix[r,c])
        end
    end
    return m
end

Matrix(x::Array{ArbFloat, 2}) = x
Matrix(x::Array{ArbFloat{P}, 2}) where {P} = x

@inline rowcount(x::ArbFloatMatrix{P}) where {P} = getfield(x.arbrealmatrix, :rowcount)
@inline colcount(x::ArbFloatMatrix{P}) where {P} = getfield(x.arbrealmatrix, :colcount)
@inline eachcell(x::ArbFloatMatrix{P}) where {P} = getfield(x.arbrealmatrix, :eachcell)
@inline eachrow(x::ArbFloatMatrix{P}) where {P} = getfield(x.arbrealmatrix, :eachrow)

@inline cellvalue(x::ArbFloatMatrix{P}, row::Int, col::Int) where {P} = eachrow(x)[row][col-1]

arbzeros(::Type{ArbFloat{P}},rowcount::SI, colcount::SI) where {P, SI<:Signed} =
    ArbFloatMatrix{P}(rowcount, colcount)

arbzeros(::Type{ArbFloat},rowcount::SI, colcount::SI) where {P, SI<:Signed} =
    ArbFloatMatrix(rowcount, colcount)

@inline Base.isempty(x::ArbFloatMatrix{P}) where {P} =
    rowcount(x) === ArbMatIdx0 || colcount(x) === ArbMatIdx0

@inline Base.size(x::ArbFloatMatrix{P}) where {P} = (rowcount(x), colcount(x))

@inline function issquare(x::ArbFloatMatrix{P}) where {P}
	return rowcount(x) === colcount(x)
end

@inline function Base.getindex(x::ArbFloatMatrix{P}, linearidx::Int) where {P}
    rowidx, colidx = linear_to_cartesian(rowcount(x), linearidx)
    return getindex(x, rowidx, colidx)
end

@inline function Base.getindex(x::ArbFloatMatrix{P}, rowidx::Int, colidx::Int) where {P}
    checkbounds(x.arbrealmatrix, rowidx, colidx)
    return getindexˌ(x.arbrealmatrix, rowidx, colidx)
end

function Base.setindex!(x::ArbFloatMatrix{P}, z::ArbFloat{P}, linearidx::Int) where {P}
    rowidx, colidx = linear_to_cartesian(rowcount(x), linearidx)
    return setindex!(x.arbrealmatrix, ArbReal{P}(z), rowidx, colidx)
end

Base.setindex!(x::ArbFloatMatrix{P1}, z::ArbFloat{P2}, linearidx::Int) where {P1,P2} =
    setindex!(x.arbrealmatrix, ArbReal{P1}(z), linearidx)

Base.setindex!(x::ArbFloatMatrix{P}, z::F, linearidx::Int) where {P, F<:Union{Signed,IEEEFloat}} =
    setindex!(x.arbrealmatrix, ArbReal{P}(z), linearidx)

Base.setindex!(x::ArbFloatMatrix{P}, z::R, linearidx::Int) where {P, R<:Real} =
    setindex!(x.arbrealmatrix, ArbReal{P}(BigFloat(z)), linearidx)

function Base.setindex!(x::ArbFloatMatrix{P}, z::ArbFloat{P},
                        rowidx::Int, colidx::Int) where {P}
    checkbounds(x, rowidx, colidx)
    return setindexˌ!(x.arbrealmatrix, ArbReal{P}(z), rowidx, colidx)
end

Base.setindex!(x::ArbFloatMatrix{P1}, z::ArbFloat{P2}, rowidx::Int, colidx::Int) where {P1,P2} =
    setindex!(x.arbrealmatrix, ArbReal{P1}(z), rowidx, colidx)

Base.setindex!(x::ArbFloatMatrix{P}, z::F, rowidx::Int, colidx::Int) where {P, F<:Union{Signed,IEEEFloat}} =
    setindex!(x.arbrealmatrix, ArbReal{P}(z), rowidx, colidx)

Base.setindex!(x::ArbFloatMatrix{P}, z::R, rowidx::Int, colidx::Int) where {P, R<:Real} =
    setindex!(x.arbrealmatrix, ArbReal{P}(BigFloat(z)), rowidx, colidx)

# transpose

#=
function transpose(x::Array{ArbFloat,2})
    y = transpose(Array{ArbReal,2}(x))
    z = Array{ArbFloat,2}(undef, size(y))
    z[:,:] = y[:,:]
    return (Transpose{ArbFloat,Array{ArbFloat,2}}).(z)
end
=#

# matrix multiply

if Sys.CPU_THREADS <= 2

    function mul!(z::ArbFloatMatrix{P}, x::ArbFloatMatrix{P}, y::ArbFloatMatrix{P}) where {P}
        ccall(@libarb(arb_mat_approx_mul), Cvoid, (Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Cint),
               z.arbrealmatrix, x.arbrealmatrix, y.arbrealmatrix, P)
        return nothing
    end

    function matmul(x::ArbFloatMatrix{P}, y::ArbFloatMatrix{P}) where {P}
        z = ArbFloatMatrix{P}(rowcount(x), colcount(y))
        ccall(@libarb(arb_mat_approx_mul), Cvoid, (Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Cint),
              z.arbrealmatrix, x.arbrealmatrix, y.arbrealmatrix, P)
        return ArbFloat{P}.(Matrix(z))
    end

else

    function mul!(z::ArbFloatMatrix{P}, x::ArbFloatMatrix{P}, y::ArbFloatMatrix{P}) where {P}
        ccall(@libarb(arb_mat_mul_threaded), Cvoid, (Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Cint),
               z.arbrealmatrix, x.arbrealmatrix, y.arbrealmatrix, P)
        return nothing
    end

    function matmul(x::ArbFloatMatrix{P}, y::ArbFloatMatrix{P}) where {P}
        z = ArbFloatMatrix{P}(rowcount(x), colcount(y))
        ccall(@libarb(arb_mat_mul_threaded), Cvoid, (Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Cint),
              z.arbrealmatrix, x.arbrealmatrix, y.arbrealmatrix, P)
        return ArbFloat{P}.(Matrix(z))
    end

end

function matmul(x::Array{ArbFloat, 2}, y::Array{ArbFloat, 2})
    xx = ArbFloatMatrix(x)
    yy = ArbFloatMatrix(x)
    xy = *(xx, yy)
    return Matrix(xy)
end

function matmul(x::Array{ArbFloat{P}, 2}, y::Array{ArbFloat{P},2}) where {P}
    xx = ArbFloatMatrix{P}(x)
    yy = ArbFloatMatrix{P}(x)
    xy = *(xx, yy)
    return Matrix(xy)
end

function Base.:(*)(x::ArbFloatMatrix{P}, y::ArbFloatMatrix{P}) where {P}
    checkmulable(x, y)
    return matmul(x, y)
end

@inline function Base.:(*)(x::Array{ArbFloat{P},2}, y::Array{ArbFloat{P},2}) where {P}
    checkmulable(x, y)
    return matmul(ArbFloatMatrix{P}(x), ArbFloatMatrix{P}(y))
end

@inline function Base.:(*)(x::Array{ArbFloat,2}, y::Array{ArbFloat,2})
    checkmulable(x, y)
    P = workingprecision(ArbFloat)
    return matmul(ArbFloatMatrix{P}(x), ArbFloatMatrix{P}(y))
end

# checks for validity

@inline function checkbounds(x::ArbFloatMatrix{P}, r::Int, c::Int) where {P}
    withinbounds = 0 < r <= rowcount(x) && 0 < c <= colcount(x)
    withinbounds && return nothing
    throw(BoundsError("($r, $c) not in 1:$(rowcount(x)), 1:$(colcount(x))"))
end

@inline function checkbounds(x::ArbFloatMatrix{P}, linear_rc::Int) where {P}
    withinbounds = 0 < linear_rc <= rowcount(x) * colcount(x)
    withinbounds && return nothing
    throw(BoundsError("($rc) not in 1:$(rowcount(x) * colcount(x))"))
end

@inline function checksquare(x::ArbFloatMatrix{P}) where {P}
    issquare(x) && return nothing
    throw(DomainError("matrix is not square ($rowcount(x), $colcount(x))"))
end

@inline function checkcompat(x::ArbFloatMatrix{P}, y::ArbFloatMatrix{P}) where {P}
    compat = (rowcount(x) === colcount(y) && rowcount(y) === colcount(x))
    compat && return nothing
    throw(DomainError("incompatible matrices: ($rowcount(x), $colcount(x)), ($rowcount(y), $colcount(y))"))
end

@inline function checkmulable(x::ArbFloatMatrix{P}, y::ArbFloatMatrix{P}) where {P}
    mulable = colcount(x) === rowcount(y)
    mulable && return nothing
    throw(ErrorException("Dimension Mismatach: ($rowcount(x), $colcount(x)), ($rowcount(y), $colcount(y))"))
end

# matrix functions: exp

function exp(x::Array{ArbFloat{P},2}) where {P}
    checksquare(x)
    y = ArbFloatMatrix{P}(x)
    z = ArbFloatMatrix{P}(rowcount(x), colcount(x))
    ccall(@libarb(arb_mat_exp), Cvoid, (Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Cint),
      z.arbrealmatrix, y.arbrealmatrix, P)
    return ArbFloat{P}.(Matrix(z))
end

function exp(x::Array{ArbFloat,2})
    checksquare(x)
    P = workingprecision(ArbFloat)
    y = ArbFloatMatrix{P}(x)
    z = ArbFloatMatrix{P}(rowcount(x), colcount(x))
    ccall(@libarb(arb_mat_exp), Cvoid, (Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Cint),
      z.arbrealmatrix, y.arbrealmatrix, P)
    return ArbFloat{P}.(Matrix(z))
end

#=
function exp(x::ArbFloatMatrix{P}) where {P}
    y = ArbFloatMatrix{P}(x)
    z = ArbFloatMatrix{P}(rowcount(x), colcount(x))
    ccall(@libarb(arb_mat_exp), Cvoid, (Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Cint),
      z.arbrealmatrix, y.arbrealmatrix, P)
    return ArbFloat{P}.(Matrix(z))
end
=#


text:
# These routines should be added to CartesianGrids

import CartesianGrids: cross!

"""
    cross!(C::VectorData,A::ScalarData/VectorData,B::VectorData/ScalarData) -> VectorData

Compute the cross product between the point data `A` and `B`, one of which
is scalar data and treated as an out-of-plane component of a vector, while
the other is in-plane vector data, and return the result as vector data `C`.
"""
function cross!(C::VectorData{N},A::ScalarData{N},B::VectorData{N}) where {N}
    @. C.u = -A*B.v
    @. C.v = A*B.u
    return C
end

function cross!(C::VectorData{N},A::VectorData{N},B::ScalarData{N}) where {N}
    @. C.u = A.v*B
    @. C.v = -A.u*B
    return C
end


text:
""" Limits and colimits in a category.
"""
module Limits
export AbstractLimit, AbstractColimit, Limit, Colimit,
  LimitAlgorithm, ColimitAlgorithm,
  ob, cone, cocone, apex, legs, limit, colimit, universal,
  Terminal, Initial, terminal, initial, delete, create, factorize,
  BinaryProduct, Product, product, proj1, proj2, pair,
  BinaryPullback, Pullback, pullback,
  BinaryEqualizer, Equalizer, equalizer, incl,
  BinaryCoproduct, Coproduct, coproduct, coproj1, coproj2, copair,
  BinaryPushout, Pushout, pushout, pushout_complement, can_pushout_complement,
  BinaryCoequalizer, Coequalizer, coequalizer, proj,
  @cartesian_monoidal_instance, @cocartesian_monoidal_instance,
  ComposeProductEqualizer, ComposeCoproductCoequalizer,
  SpecializeLimit, SpecializeColimit, ToBipartiteLimit, ToBipartiteColimit

using AutoHashEquals
using StaticArrays: StaticVector, SVector

using ...GAT, ...Theories
import ...Theories: ob, terminal, product, proj1, proj2, equalizer, incl,
  initial, coproduct, coproj1, coproj2, coequalizer, proj,
  delete, create, pair, copair, factorize
using ...CSetDataStructures, ..FinCats, ..FreeDiagrams
import ..FreeDiagrams: apex, legs

# Data types for limits
#######################

""" Abstract type for limit in a category.

The standard concrete subtype is [`Limit`](@ref), although for computational
reasons certain categories may use different subtypes to include extra data.
"""
abstract type AbstractLimit{Ob,Diagram} end

ob(lim::AbstractLimit) = apex(lim)
cone(lim::AbstractLimit) = lim.cone
apex(lim::AbstractLimit) = apex(cone(lim))
legs(lim::AbstractLimit) = legs(cone(lim))

Base.iterate(lim::AbstractLimit, args...) = iterate(cone(lim), args...)
Base.eltype(lim::AbstractLimit) = eltype(cone(lim))
Base.length(lim::AbstractLimit) = length(cone(lim))

""" Limit in a category.
"""
@auto_hash_equals struct Limit{Ob,Diagram,Cone<:Multispan{Ob}} <:
    AbstractLimit{Ob,Diagram}
  diagram::Diagram
  cone::Cone
end

const Terminal{Ob} = AbstractLimit{Ob,<:EmptyDiagram}
const BinaryProduct{Ob} = AbstractLimit{Ob,<:ObjectPair}
const Product{Ob} = AbstractLimit{Ob,<:DiscreteDiagram}
const BinaryPullback{Ob} = AbstractLimit{Ob,<:Cospan}
const Pullback{Ob} = AbstractLimit{Ob,<:Multicospan}
const BinaryEqualizer{Ob} = AbstractLimit{Ob,<:ParallelPair}
const Equalizer{Ob} = AbstractLimit{Ob,<:ParallelMorphisms}

proj1(lim::Union{BinaryProduct,BinaryPullback}) = first(legs(lim))
proj2(lim::Union{BinaryProduct,BinaryPullback}) = last(legs(lim))
incl(eq::Equalizer) = only(legs(eq))

""" Algorithm for computing limits.
"""
abstract type LimitAlgorithm end

# Data types for colimits
#########################

""" Abstract type for colimit in a category.

The standard concrete subtype is [`Colimit`](@ref), although for computational
reasons certain categories may use different subtypes to include extra data.
"""
abstract type AbstractColimit{Ob,Diagram} end

ob(colim::AbstractColimit) = apex(colim)
cocone(colim::AbstractColimit) = colim.cocone
apex(colim::AbstractColimit) = apex(cocone(colim))
legs(colim::AbstractColimit) = legs(cocone(colim))

Base.iterate(colim::AbstractColimit, args...) = iterate(cocone(colim), args...)
Base.eltype(colim::AbstractColimit) = eltype(cocone(colim))
Base.length(colim::AbstractColimit) = length(cocone(colim))

""" Colimit in a category.
"""
@auto_hash_equals struct Colimit{Ob,Diagram,Cocone<:Multicospan{Ob}} <:
    AbstractColimit{Ob,Diagram}
  diagram::Diagram
  cocone::Cocone
end

const Initial{Ob} = AbstractColimit{Ob,<:EmptyDiagram}
const BinaryCoproduct{Ob} = AbstractColimit{Ob,<:ObjectPair}
const Coproduct{Ob} = AbstractColimit{Ob,<:DiscreteDiagram}
const BinaryPushout{Ob} = AbstractColimit{Ob,<:Span}
const Pushout{Ob} = AbstractColimit{Ob,<:Multispan}
const BinaryCoequalizer{Ob} = AbstractColimit{Ob,<:ParallelPair}
const Coequalizer{Ob} = AbstractColimit{Ob,<:ParallelMorphisms}

coproj1(colim::Union{BinaryCoproduct,BinaryPushout}) = first(legs(colim))
coproj2(colim::Union{BinaryCoproduct,BinaryPushout}) = last(legs(colim))
proj(coeq::Coequalizer) = only(legs(coeq))

""" Algorithm for computing colimits.
"""
abstract type ColimitAlgorithm end

# Generic (co)limits
####################

""" Limit of a diagram.

To define limits in a category with objects `Ob`, override the method
`limit(::FreeDiagram{Ob})` for general limits or `limit(::D)` with suitable type
`D <: FixedShapeFreeDiagram{Ob}` for limits of specific shape, such as products
or equalizers.

See also: [`colimit`](@ref)
"""
limit(diagram; kw...) = limit(diagram_type(diagram), diagram; kw...)
limit(diagram, ::Nothing; kw...) = limit(diagram; kw...) # alg == nothing

""" Colimit of a diagram.

To define colimits in a category with objects `Ob`, override the method
`colimit(::FreeDiagram{Ob})` for general colimits or `colimit(::D)` with
suitable type `D <: FixedShapeFreeDiagram{Ob}` for colimits of specific shape,
such as coproducts or coequalizers.

See also: [`limit`](@ref)
"""
colimit(diagram; kw...) = colimit(diagram_type(diagram), diagram; kw...)
colimit(diagram, ::Nothing; kw...) = colimit(diagram; kw...) # alg == nothing

""" Universal property of (co)limits.

Compute the morphism whose existence and uniqueness is guaranteed by the
universal property of (co)limits.

See also: [`limit`](@ref), [`colimit`](@ref).
"""
function universal end

# Specific (co)limits
#####################

""" Terminal object.

To implement for a type `T`, define the method `limit(::EmptyDiagram{T})`.
"""
terminal(T::Type; kw...) = limit(EmptyDiagram{T}(); kw...)

""" Initial object.

To implement for a type `T`, define the method `colimit(::EmptyDiagram{T})`.
"""
initial(T::Type; kw...) = colimit(EmptyDiagram{T}(); kw...)

""" Unique morphism into a terminal object.

To implement for a type `T`, define the method
`universal(::Terminal{T}, ::SMultispan{0,T})`.
"""
delete(A::T) where T = delete(terminal(T), A)
delete(lim::Terminal, A) = universal(lim, SMultispan{0}(A))

""" Unique morphism out of an initial object.

To implement for a type `T`, define the method
`universal(::Initial{T}, ::SMulticospan{0,T})`.
"""
create(A::T) where T = create(initial(T), A)
create(colim::Initial, A) = universal(colim, SMulticospan{0}(A))

""" Product of objects.

To implement for a type `T`, define the method `limit(::ObjectPair{T})` and/or
`limit(::DiscreteDiagram{T})`.
"""
product(A, B; kw...) = limit(ObjectPair(A, B); kw...)
product(As::AbstractVector; kw...) = limit(DiscreteDiagram(As); kw...)

""" Coproduct of objects.

To implement for a type `T`, define the method `colimit(::ObjectPair{T})` and/or
`colimit(::DiscreteDiagram{T})`.
"""
coproduct(A, B; kw...) = colimit(ObjectPair(A, B); kw...)
coproduct(As::AbstractVector; kw...) = colimit(DiscreteDiagram(As); kw...)

""" Equalizer of morphisms with common domain and codomain.

To implement for a type `T`, define the method `limit(::ParallelPair{T})` and/or
`limit(::ParallelMorphisms{T})`.
"""
equalizer(f, g; kw...) = limit(ParallelPair(f, g); kw...)
equalizer(fs::AbstractVector; kw...) = limit(ParallelMorphisms(fs); kw...)

""" Coequalizer of morphisms with common domain and codomain.

To implement for a type `T`, define the method `colimit(::ParallelPair{T})` or
`colimit(::ParallelMorphisms{T})`.
"""
coequalizer(f, g; kw...) = colimit(ParallelPair(f, g); kw...)
coequalizer(fs::AbstractVector; kw...) = colimit(ParallelMorphisms(fs); kw...)

""" Pullback of a pair of morphisms with common codomain.

To implement for a type `T`, define the method `limit(::Cospan{T})` and/or
`limit(::Multicospan{T})` or, if you have already implemented products and
equalizers, rely on the default implementation.
"""
pullback(f, g; kw...) = limit(Cospan(f, g); kw...)
pullback(fs::AbstractVector; kw...) = limit(Multicospan(fs); kw...)

""" Pushout of a pair of morphisms with common domain.

To implement for a type `T`, define the method `colimit(::Span{T})` and/or
`colimit(::Multispan{T})` or, if you have already implemented coproducts and
coequalizers, rely on the default implementation.
"""
pushout(f, g; kw...) = colimit(Span(f, g); kw...)
pushout(fs::AbstractVector; kw...) = colimit(Multispan(fs); kw...)

""" Pairing of morphisms: universal property of products/pullbacks.

To implement for products of type `T`, define the method
`universal(::BinaryProduct{T}, ::Span{T})` and/or
`universal(::Product{T}, ::Multispan{T})` and similarly for pullbacks.
"""
pair(f, g) = pair(product(codom(f), codom(g)), f, g)
pair(fs::AbstractVector) = pair(product(map(codom, fs)), fs)
pair(lim::Union{BinaryProduct,BinaryPullback}, f, g) =
  universal(lim, Span(f, g))
pair(lim::Union{Product,Pullback}, fs::AbstractVector) =
  universal(lim, Multispan(fs))

""" Copairing of morphisms: universal property of coproducts/pushouts.

To implement for coproducts of type `T`, define the method
`universal(::BinaryCoproduct{T}, ::Cospan{T})` and/or
`universal(::Coproduct{T}, ::Multicospan{T})` and similarly for pushouts.
"""
copair(f, g) = copair(coproduct(dom(f), dom(g)), f, g)
copair(fs::AbstractVector) = copair(coproduct(map(dom, fs)), fs)
copair(colim::Union{BinaryCoproduct,BinaryPushout}, f, g) =
  universal(colim, Cospan(f, g))
copair(colim::Union{Coproduct,Pushout}, fs::AbstractVector) =
  universal(colim, Multicospan(fs))

""" Factor morphism through (co)equalizer, via the universal property.

To implement for equalizers of type `T`, define the method
`universal(::Equalizer{T}, ::SMultispan{1,T})`. For coequalizers of type `T`,
define the method `universal(::Coequalizer{T}, ::SMulticospan{1,T})`.
"""
factorize(lim::Equalizer, h) = universal(lim, SMultispan{1}(h))
factorize(colim::Coequalizer, h) = universal(colim, SMulticospan{1}(h))

""" Pushout complement: extend composable pair to a pushout square.

[Pushout complements](https://ncatlab.org/nlab/show/pushout+complement) are the
essential ingredient for double pushout (DPO) rewriting.
"""
pushout_complement(f, g) = pushout_complement(ComposablePair(f, g))

""" Can a pushout complement be constructed for a composable pair?

Even in nice categories, this is not generally possible.
"""
can_pushout_complement(f, g) = can_pushout_complement(ComposablePair(f, g))

# (Co)cartesian monoidal categories
###################################

""" Define cartesian monoidal structure using limits.

Implements an instance of [`CartesianCategory`](@ref) assuming that finite
products have been implemented following the limits interface.
"""
macro cartesian_monoidal_instance(Ob, Hom)
  esc(quote
    import Catlab.Theories: CartesianCategory, otimes, ⊗, munit, braid, σ,
      mcopy, delete, pair, proj1, proj2, Δ, ◊

    @instance CartesianCategory{$Ob, $Hom} begin
      @import dom, codom, compose, ⋅, id, munit, delete, pair

      otimes(A::$Ob, B::$Ob) = ob(product(A, B))

      function otimes(f::$Hom, g::$Hom)
        π1, π2 = product(dom(f), dom(g))
        pair(product(codom(f), codom(g)), π1⋅f, π2⋅g)
      end

      function braid(A::$Ob, B::$Ob)
        AB, BA = product(A, B), product(B, A)
        pair(BA, proj2(AB), proj1(AB))
      end

      mcopy(A::$Ob) = pair(id(A),id(A))
      proj1(A::$Ob, B::$Ob) = proj1(product(A, B))
      proj2(A::$Ob, B::$Ob) = proj2(product(A, B))
    end

    otimes(As::AbstractVector{<:$Ob}) = ob(product(As))

    function otimes(fs::AbstractVector{<:$Hom})
      Π, Π′ = product(map(dom, fs)), product(map(codom, fs))
      pair(Π′, map(compose, legs(Π), fs))
    end

    munit(::Type{T}) where T <: $Ob = ob(terminal(T))
  end)
end

""" Define cocartesian monoidal structure using colimits.

Implements an instance of [`CocartesianCategory`](@ref) assuming that finite
coproducts have been implemented following the colimits interface.
"""
macro cocartesian_monoidal_instance(Ob, Hom)
  esc(quote
    import Catlab.Theories: CocartesianCategory, oplus, ⊕, mzero, swap,
      plus, zero, copair, coproj1, coproj2

    @instance CocartesianCategory{$Ob, $Hom} begin
      @import dom, codom, compose, ⋅, id, mzero, copair

      oplus(A::$Ob, B::$Ob) = ob(coproduct(A, B))

      function oplus(f::$Hom, g::$Hom)
        ι1, ι2 = coproduct(codom(f), codom(g))
        copair(coproduct(dom(f), dom(g)), f⋅ι1, g⋅ι2)
      end

      function swap(A::$Ob, B::$Ob)
        AB, BA = coproduct(A, B), coproduct(B, A)
        copair(AB, coproj2(BA), coproj1(BA))
      end

      plus(A::$Ob) = copair(id(A),id(A))
      zero(A::$Ob) = create(A)
      coproj1(A::$Ob, B::$Ob) = coproj1(coproduct(A, B))
      coproj2(A::$Ob, B::$Ob) = coproj2(coproduct(A, B))
    end

    oplus(As::AbstractVector{<:$Ob}) = ob(coproduct(As))

    function oplus(fs::AbstractVector{<:$Hom})
      ⊔, ⊔′ = coproduct(map(dom, fs)), coproduct(map(codom, fs))
      copair(⊔, map(compose, fs, legs(⊔′)))
    end

    mzero(::Type{T}) where T <: $Ob = ob(initial(T))
  end)
end

# (Co)limit algorithms
######################

# Specialize (co)limits
#----------------------

""" Meta-algorithm attempting to reduce general limits to common special cases.

Specifically, reduce limits of free diagrams that happen to be discrete to
products. TODO: Handle pullbacks similarly.
"""
@Base.kwdef struct SpecializeLimit <: LimitAlgorithm
  fallback::Union{LimitAlgorithm,Nothing} = nothing
end

limit(diagram, alg::SpecializeLimit) = limit(diagram, alg.fallback)

function limit(F::FinDomFunctor{<:FinCat{Int},<:TypeCat{Ob,Hom}},
               alg::SpecializeLimit) where {Ob,Hom}
  if is_discrete(dom(F))
    return limit(DiscreteDiagram(collect_ob(F), Hom), alg)
  end
  limit(F, alg.fallback)
end
limit(diagram::FreeDiagram, alg::SpecializeLimit) =
  limit(FinDomFunctor(diagram), alg)

function limit(Xs::DiscreteDiagram{Ob,Hom,Obs},
               alg::SpecializeLimit) where {Ob,Hom,Obs}
  if !(Obs <: StaticVector) && length(Xs) <= 3
    return limit(DiscreteDiagram(SVector{length(Xs)}(ob(Xs)), Hom), alg)
  end
  limit(Xs, alg.fallback)
end

""" Meta-algorithm attempting to reduce general colimits to common special cases.

Dual to [`SpecializeLimit`](@ref).
"""
@Base.kwdef struct SpecializeColimit <: ColimitAlgorithm
  fallback::Union{ColimitAlgorithm,Nothing} = nothing
end

colimit(diagram, alg::SpecializeColimit) = colimit(diagram, alg.fallback)

function colimit(F::FinDomFunctor{<:FinCat{Int},<:TypeCat{Ob,Hom}},
                 alg::SpecializeColimit) where {Ob,Hom}
  if is_discrete(dom(F))
    return colimit(DiscreteDiagram(collect_ob(F), Hom), alg)
  end
  colimit(F, alg.fallback)
end
colimit(diagram::FreeDiagram, alg::SpecializeColimit) =
  colimit(FinDomFunctor(diagram), alg)

function colimit(Xs::DiscreteDiagram{Ob,Hom,Obs},
                 alg::SpecializeColimit) where {Ob,Hom,Obs}
  if !(Obs <: StaticVector) && length(Xs) <= 3
    return colimit(DiscreteDiagram(SVector{length(Xs)}(ob(Xs)), Hom), alg)
  end
  colimit(Xs, alg.fallback)
end

""" Limit of a singleton diagram.
"""
struct SingletonLimit{Ob,Diagram<:SingletonDiagram{Ob}} <: AbstractLimit{Ob,Diagram}
  diagram::Diagram
end

cone(lim::SingletonLimit) = let x = only(lim.diagram)
  SMultispan{1}(x, id(x))
end
universal(::SingletonLimit, cone::Multispan) = only(cone)

limit(diagram::SingletonDiagram, ::SpecializeLimit) = SingletonLimit(diagram)

""" Colimit of a singleton diagram.
"""
struct SingletonColimit{Ob,Diagram<:SingletonDiagram{Ob}} <: AbstractColimit{Ob,Diagram}
  diagram::Diagram
end

cocone(colim::SingletonColimit) = let x = only(colim.diagram)
  SMulticospan{1}(x, id(x))
end
universal(::SingletonColimit, cocone::Multicospan) = only(cocone)

colimit(diagram::SingletonDiagram, ::SpecializeColimit) =
  SingletonColimit(diagram)

# Composite (co)limits
#---------------------

""" Compute pullback by composing a product with an equalizer.

See also: [`ComposeCoproductCoequalizer`](@ref).
"""
struct ComposeProductEqualizer <: LimitAlgorithm end

""" Pullback formed as composite of product and equalizer.

The fields of this struct are an implementation detail; accessing them directly
violates the abstraction. Everything that you can do with a pushout, including
invoking its universal property, should be done through the generic interface
for limits.

See also: [`CompositePushout`](@ref).
"""
struct CompositePullback{Ob, Diagram<:Multicospan, Cone<:Multispan{Ob},
    Prod<:Product, Eq<:Equalizer} <: AbstractLimit{Ob,Diagram}
  diagram::Diagram
  cone::Cone
  prod::Prod
  eq::Eq
end

function limit(cospan::Multicospan, ::ComposeProductEqualizer)
  prod = product(feet(cospan))
  (ι,) = eq = equalizer(map(compose, legs(prod), legs(cospan)))
  cone = Multispan(map(π -> ι⋅π, legs(prod)))
  CompositePullback(cospan, cone, prod, eq)
end

function universal(lim::CompositePullback, cone::Multispan)
  factorize(lim.eq, universal(lim.prod, cone))
end

""" Compute pushout by composing a coproduct with a coequalizer.

See also: [`ComposeProductEqualizer`](@ref).
"""
struct ComposeCoproductCoequalizer <: ColimitAlgorithm end

""" Pushout formed as composite of coproduct and equalizer.

See also: [`CompositePullback`](@ref).
"""
struct CompositePushout{Ob, Diagram<:Multispan, Cocone<:Multicospan{Ob},
    Coprod<:Coproduct, Coeq<:Coequalizer} <: AbstractColimit{Ob,Diagram}
  diagram::Diagram
  cocone::Cocone
  coprod::Coprod
  coeq::Coeq
end

function colimit(span::Multispan, ::ComposeCoproductCoequalizer)
  coprod = coproduct(feet(span))
  (π,) = coeq = coequalizer(map(compose, legs(span), legs(coprod)))
  cocone = Multicospan(map(ι -> ι⋅π, legs(coprod)))
  CompositePushout(span, cocone, coprod, coeq)
end

function universal(lim::CompositePushout, cone::Multicospan)
  factorize(lim.coeq, universal(lim.coprod, cone))
end

# Bipartite (co)limits
#---------------------

""" Compute a limit by reducing the diagram to a free bipartite diagram.
"""
struct ToBipartiteLimit <: LimitAlgorithm end

""" Limit computed by reduction to the limit of a free bipartite diagram.
"""
struct BipartiteLimit{Ob, Diagram, Cone<:Multispan{Ob},
                      Lim<:AbstractLimit} <: AbstractLimit{Ob,Diagram}
  diagram::Diagram
  cone::Cone
  limit::Lim
end

function limit(F::Union{Functor,FreeDiagram}, ::ToBipartiteLimit)
  d = BipartiteFreeDiagram(F)
  lim = limit(d)
  cone = Multispan(apex(lim), map(incident(d, :, :orig_vert₁),
                                  incident(d, :, :orig_vert₂)) do v₁, v₂
    if v₁ == 0
      e = first(incident(d, v₂, :tgt))
      compose(legs(lim)[src(d, e)], hom(d, e))
    else
      legs(lim)[v₁]
    end
  end)
  BipartiteLimit(F, cone, lim)
end

function universal(lim::BipartiteLimit, cone::Multispan)
  lim = lim.limit
  cone = Multispan(apex(cone), legs(cone)[lim.diagram[:orig_vert₁]])
  universal(lim, cone)
end

""" Compute a colimit by reducing the diagram to a free bipartite diagram.
"""
struct ToBipartiteColimit <: ColimitAlgorithm end

""" Limit computed by reduction to the limit of a free bipartite diagram.
"""
struct BipartiteColimit{Ob, Diagram, Cocone<:Multicospan{Ob},
                        Colim<:AbstractColimit} <: AbstractColimit{Ob,Diagram}
  diagram::Diagram
  cocone::Cocone
  colimit::Colim
end

function colimit(F::Union{Functor,FreeDiagram}, ::ToBipartiteColimit)
  d = BipartiteFreeDiagram(F, colimit=true)
  colim = colimit(d)
  cocone = Multicospan(apex(colim), map(incident(d, :, :orig_vert₁),
                                        incident(d, :, :orig_vert₂)) do v₁, v₂
    if v₂ == 0
      e = first(incident(d, v₁, :src))
      compose(hom(d, e), legs(colim)[tgt(d, e)])
    else
      legs(colim)[v₂]
    end
  end)
  BipartiteColimit(F, cocone, colim)
end

function universal(colim::BipartiteColimit, cocone::Multicospan)
  colim = colim.colimit
  cocone = Multicospan(apex(cocone), legs(cocone)[colim.diagram[:orig_vert₂]])
  universal(colim, cocone)
end

end


text:
"""
    komptrap(f,a,b,M)
adalah fungsi yang digunakan untuk menghitung nilai integral `f` secara numerik
pada interval `[a,b]` menggunakan komposit trapesium dengan `M` sub-interval.

# Example
```jl
julia> a = 1; b = 6;

julia> f(x) = 2+sin(2*sqrt(x));

julia> y = komptrap(f,a,b,10)
8.193854565172531
```
return solusi `y`
"""
function komptrap(f,a,b,M)
  h = (b-a)/M;
  s = 0;
  for k = 1:M-1
    x = a+k*h;
    s = s+f(x);
  end
  y = h/2*(f(a)+f(b)+2*s)
end


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

struct DummyVariableBridge <: MOIB.Variable.AbstractBridge end

@testset "AbstractBridge" begin
    model = MOIU.Model{Float64}()
    bridge = DummyVariableBridge()
    attr = MOI.VariablePrimalStart()
    @test !MOI.supports(model, attr, typeof(bridge))
    i = MOIB.Variable.IndexInVector(1)
    @test_throws MOI.UnsupportedAttribute(attr) MOI.set(
        model,
        attr,
        bridge,
        1.0,
    )
    @test_throws MOI.UnsupportedAttribute(attr) MOI.set(
        model,
        attr,
        bridge,
        1.0,
        i,
    )
    attr = MOI.VariablePrimal()
    err = MOI.SetAttributeNotAllowed(attr)
    @test_throws err MOI.set(model, attr, bridge, 1.0)
    @test_throws err MOI.set(model, attr, bridge, 1.0, i)
end


text:
module TupleBenchmarks

include(joinpath(dirname(@__FILE__), "..", "utils", "RandUtils.jl"))

using .RandUtils
using BenchmarkTools
using Compat

import Compat: UTF8String, view

const SUITE = BenchmarkGroup()

###############
# issue #5274 #
###############

immutable TupleWrapper{N, T}
    data::NTuple{N, T}
end

Base.eltype{N,T}(::TupleWrapper{N,T}) = T
Base.length{N,T}(::TupleWrapper{N,T}) = N

function get_index(n::NTuple, i::Int)
    @inbounds v = n[i]
    return v
end

function get_index(n::TupleWrapper, i::Int)
    @inbounds v = n.data[i]
    return v
end

function sum_tuple{N, T}(n::Union{NTuple{N, T}, TupleWrapper{N, T}})
    s = zero(T)
    for i in 1:N
        s += get_index(n, i)
    end
    return s
end

const TUPLE_SUM_SIZES = (3, 8, 30, 60)
const TUPLE_SUM_TYPES = (Float32, Float64)

g = addgroup!(SUITE, "index", ["sum"])

for s in TUPLE_SUM_SIZES, T in TUPLE_SUM_TYPES
    tup = tuple(samerand(T, s)...)
    tupwrap = TupleWrapper(tup)
    g["sumelt", "NTuple", s, T] = @benchmarkable sum_tuple($tup) time_tolerance=0.40
    g["sumelt", "TupleWrapper", s, T] = @benchmarkable sum_tuple($tupwrap) time_tolerance=0.40
end

#####################
# Fixed Size Arrays #
#####################

# Short fixed size array implementation

@compat abstract type FixedArray{T, N} <: AbstractArray{T, N} end

@compat Base.IndexStyle(::Type{<: FixedArray}) = IndexLinear()
Base.getindex(fsa::FixedArray, i::Int) = fsa.data[i]

immutable FixedVector{L, T} <: FixedArray{T, 1}
    data::NTuple{L, T}
end

Base.size{L}(::FixedVector{L}) = (L,)
Base.size{L, T}(::Type{FixedVector{L, T}}) = (L,)
Base.length{L}(::FixedVector{L}) = L

immutable FixedMatrix{R, C, T, RC} <: FixedArray{T, 2}
    data::NTuple{RC, T}
end

Base.size{R, C}(::FixedMatrix{R, C}) = (R, C)
Base.size{R, C, T, RC}(::Type{FixedMatrix{R, C, T, RC}}) = (R, C)
Base.length{R, C, T, RC}(::FixedMatrix{R, C, T, RC}) = RC

# Reductions

@inline function perf_reduce(op, a::FixedArray)
    if length(a) == 1
        return a[1]
    else
        s = op(a[1], a[2])
        for j = 3:length(a)
            s = op(s, a[j])
        end
        return s
    end
end

perf_minimum(a::FixedArray) = perf_reduce(min, a)

@inline function perf_reduce(op, v0, a::FixedArray)
    if length(a) == 0
        return v0
    else
        s = v0
        @inbounds @simd for j = 1:length(a)
            s = op(s, a[j])
        end
        return s
    end
end

perf_sum{T}(v::FixedArray{T}) = perf_reduce(+, zero(T), v)

@inline function perf_mapreduce(f, op, v0, a1::FixedArray)
    if length(a1) == 0
        return v0
    else
        s = op(v0, f(a1[1]))
        for j = 2:length(a1)
            s = op(s, f(a1[j]))
        end
        return s
    end
end

perf_sumabs2{T}(a::FixedArray{T}) = perf_mapreduce(abs2, +, zero(T), a)

# Linear Algebra

@generated function perf_matvec{R, C, T}(A::FixedMatrix{R, C, T}, b::FixedVector{C, T})
    sA = size(A)
    sB = size(b)
    exprs = Expr(:tuple, [reduce((ex1,ex2) -> :(+($ex1,$ex2)),
                [:(A[$(sub2ind(sA, k, j))]*b[$j]) for j = 1:sA[2]]) for k = 1:sA[1]]...)
    return quote
        @inbounds return FixedVector{R, T}($exprs)
    end
end

@generated function perf_matmat{R1, R2, C, T}(A::FixedMatrix{R1, C, T}, B::FixedMatrix{C, R2, T})
    sA = size(A)
    sB = size(B)
    exprs =  Expr(:tuple, [reduce((ex1,ex2) -> :(+($ex1,$ex2)),
                [:(A[$(sub2ind(sA, k1, j))] * B[$(sub2ind(sB, j, k2))]) for j = 1:sA[2]]) for k1 = 1:sA[1], k2 = 1:sB[2]]...)
    result_type = FixedMatrix{R1, R2, T, (R1 * R2)}
    return quote
        @inbounds return $result_type($exprs)
    end
end

if VERSION >= v"0.5"

# Benchmarks #
##############
v2, v4, v8, v16 = [FixedVector((rand(i)...)) for i in (2, 4, 8, 16)]
m2x2, m4x4, m8x8, m16x16 = [FixedMatrix{i,i, Float64, i*i}((rand(i*i)...)) for i in (2, 4, 8, 16)]

# Reductions
g = addgroup!(SUITE, "reduction", ["tuple"])

for mv in (v2, v4, v8, v16, m2x2, m4x4, m8x8, m16x16)
    g["sum", size(mv)] = @benchmarkable perf_sum($mv)
    g["sumabs", size(mv)] = @benchmarkable perf_sumabs2($mv)
    g["minimum", size(mv)] = @benchmarkable perf_minimum($mv)
end

# Linear algebra

g = addgroup!(SUITE, "linear algebra", ["tuple"])

for (m, v) in zip((m2x2, m4x4, m8x8, m16x16), (v2, v4, v8, v16 ))
    g["matvec", size(m), size(v)] = @benchmarkable perf_matvec($m, $v)
    g["matmat", size(m), size(m)] = @benchmarkable perf_matmat($m, $m)
end

end # version

end # module


text:
using DistanceTransforms
using Test
using TestSetExtensions


text:
__precompile__()

module Mesh

export squaremesh, cubemesh, gmshcall
export readmesh, writemesh

include("squaremesh.jl");
include("cubemesh.jl");
include("gmshcall.jl");
include("readmesh.jl");
include("writemesh.jl");

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
using JUDI, SegyIO, LinearAlgebra, PyPlot

# Set up model structure
n = (121, 101)   # (x,y,z) or (x,z)
d = (2.5f0, 2.5f0) # in mm
o = (0., 0.)

# Velocity [km/s]
v = ones(Float32,n) .+ 0.5f0
v0 = ones(Float32,n) .+ 0.5f0
v[:,Int(round(end/2)):end] .= 4f0

# Slowness squared [s^2/km^2]
m = (1f0 ./ v).^2
m0 = (1f0 ./ v0).^2
dm = vec(m - m0)

# Setup info and model structure
nsrc = 1	# number of sources
model = Model(n, d, o, m)
model0 = Model(n, d, o, m0)

## Set up receiver geometry
nxrec = 120
xrec = range(d[1], stop=d[2]*(n[1]-1), length=nxrec)
yrec = 0f0
zrec = range(d[1], stop=d[1], length=nxrec)

# receiver sampling and recording time
timeR = 250f0   # receiver recording time [ms]
dtR = 0.2f0    # receiver sampling interval [ms]

# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)

## Set up source geometry (cell array with source locations for each shot)
xsrc = d[1]*61
ysrc = 0f0
zsrc = d[1]

# source sampling and number of time steps
timeS = 250f0     # ms
dtS =  0.2f0    # ms

# Set up source structure
srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dtS, t=timeS)

# setup wavelet
f0 = .05  # kHz
wavelet = ricker_wavelet(timeS, dtS, f0)
q = judiVector(srcGeometry, wavelet)

# Set up info structure for linear operators
ntComp = get_computational_nt(srcGeometry, recGeometry, model)
info = Info(prod(n), nsrc, ntComp)

###################################################################################################

opt = Options(isic=true)
# Setup operators
Pr = judiProjection(info, recGeometry)
F = judiModeling(info, model)
F0 = judiModeling(info, model0; options=opt)
Ps = judiProjection(info, srcGeometry)
J = judiJacobian(Pr*F0*adjoint(Ps), q)

# Nonlinear modeling
dobs = Pr*F*adjoint(Ps)*q

# With a transducer source pointing down so pi/2 angle and radius 5mm (1cm diameter)
q2 = transducer(q, model.d, 5, pi/2 .* ones(q.nsrc))
Ps2 = judiProjection(info, q2.geometry)

dobs2 = Pr*F*adjoint(Ps2)*q2

a = 1e-1
figure()
subplot(121)
imshow(dobs.data[1], vmin=-a, vmax=a, cmap="seismic", aspect=.25)
title("Point source")
subplot(122)
imshow(dobs2.data[1], vmin=-a, vmax=a, cmap="seismic", aspect=.25)
title("Transducer source")

dm = J'*dobs2


text:
using Documenter, Makie

makedocs(
    modules = [Makie],
    format = :html,
    sitename = "Plotting in pure Julia",
    pages = [
        "Home" => "index.md",
        "Basics" => [
            "scene.md",
            "conversions.md",
            "functions.md",
            "documentation.md",
            "backends.md",
            "extending.md",
            "themes.md",
            "interaction.md",
            "axis.md",
            "legends.md",
            "output.md",
            "reflection.md",
            "layout.md"
        ],
        "Developper Documentation" => [
            "devdocs.md",
        ],
    ]
)

ENV["TRAVIS_BRANCH"] = "latest"
ENV["TRAVIS_PULL_REQUEST"] = "false"
ENV["TRAVIS_REPO_SLUG"] = "github.com/SimonDanisch/Makie.jl.git"
ENV["TRAVIS_TAG"] = "tag"
ENV["TRAVIS_OS_NAME"] = "linux"
ENV["TRAVIS_JULIA_VERSION"] = "0.6"

deploydocs(
    deps   = Deps.pip("mkdocs", "python-markdown-math", "mkdocs-cinder"),
    repo   = "github.com/SimonDanisch/Makie.jl.git",
    julia  = "0.6",
    target = "build",
    osname = "linux",
    make = nothing
)


text:
#SSH Cluster manager...
# machines => array of machine elements
# machine => address or (address, cnt)
# address => string of form `[user@]host[:port] bind_addr[:bind_port]`
# cnt => :auto or number
# :auto launches NUM_CORES number of workers at address
# number launches the specified number of workers at address

struct SSHMachine
    address::String #string of form `[user@]host[:port] bind_addr[:bind_port]`
    core_count::Union{Symbol, Int}
end

function machine(addr::String, core_count::Symbol)
    @assert( core_count == :auto, "Only accepted core Symbol is `:auto`, which uses NUM_CORES workers." )
    return SSHMachine(addr, core_count)
end

function machine(addr::String, core_count::Int)
    @assert( core_count >= 0, "Cannot make negative workers.." )
    return SSHMachine(addr, core_count)
end

max_stars_repo_name: caseykneale/MrPhelps

max_issues_repo_name: caseykneale/MrPhelps

max_forks_repo_name: caseykneale/MrPhelps


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>5.0),
    "X0"   => [(1, Hyperrectangle([2.0, 1], [0.5, 0]))],
    "info" => Dict(
        "name" => "two_tanks",
        "dim" => 2,
        "linear" => true,
        "hybrid" => true
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
#= src/id.jl

References:

  H. Cheng, Z. Gimbutas, P.G. Martinsson, V. Rokhlin. On the compression of low
    rank matrices. SIAM J. Sci. Comput. 26 (4): 1389-1404, 2005.

  E. Liberty, F. Woolfe, P.-G. Martinsson V. Rokhlin, M. Tygert. Randomized
    algorithms for the low-rank approximation of matrices. Proc. Natl. Acad.
    Sci. U.S.A. 104 (51): 20167-20172, 2007.
=#

# IDPackedV

mutable struct IDPackedV{S} <: Factorization{S}
  sk::Vector{Int}
  rd::Vector{Int}
  T::Matrix{S}
end

conj!(A::IDPackedV) = IDPackedV(A.sk, A.rd, conj!(A.T))
conj(A::IDPackedV) = IDPackedV(A.sk, A.rd, conj(A.T))

convert(::Type{IDPackedV{T}}, A::IDPackedV) where {T} =
  IDPackedV(A.sk, A.rd, convert(Array{T}, A.T))
convert(::Type{Factorization{T}}, A::IDPackedV) where {T} = convert(IDPackedV{T}, A)
convert(::Type{Array}, A::IDPackedV) = Matrix(A)
convert(::Type{Array{T}}, A::IDPackedV) where {T} = convert(Array{T}, Matrix(A))

copy(A::IDPackedV) = IDPackedV(copy(A.sk), copy(A.rd), copy(A.T))

  function _Matrix!(trans::Symbol, A::AbstractMatrix{T}, V::IDPackedV{T}) where T
    chktrans(trans)
    k, n = size(V)
    if trans == :n
      size(A) == (k, n) || throw(DimensionMismatch)
      @inbounds for j = 1:k
        @simd for i = 1:k
          A[i,j] = i == j ? 1 : 0
        end
      end
      A[:,k+1:n] = V[:T]
      rmul!(A, V[:P]')
    else
      size(A) == (n, k) || throw(DimensionMismatch)
      @inbounds for j = 1:k
        @simd for i = 1:k
          A[i,j] = i == j ? 1 : 0
        end
      end
      adjoint!(view(A,k+1:n,:), V[:T])
      lmul!(V[:P], A)
    end
    A
  end

_Matrix!(A::AbstractMatrix{T}, V::IDPackedV{T}) where {T} = _Matrix!(:n, A, V)
function _Matrix(trans::Symbol, A::IDPackedV{T}) where T
  chktrans(trans)
  k, n = size(A)
  if trans == :n  B = Array{T}(undef, k, n)
  else            B = Array{T}(undef, n, k)
  end
  _Matrix!(trans, B, A)
end
Matrix(A::IDPackedV) = _Matrix(:n, A)

  adjoint(A::IDPackedV) = Adjoint(A)
  transpose(A::IDPackedV) = Transpose(A)

function getindex(A::IDPackedV, d::Symbol)
  if     d == :P   return ColumnPermutation(A[:p])
  elseif d == :T   return A.T
  elseif d == :k   return length(A.sk)
  elseif d == :p   return [A.sk; A.rd]
  elseif d == :rd  return A.rd
  elseif d == :sk  return A.sk
  else             throw(KeyError(d))
  end
end

ishermitian(::IDPackedV) = false
issym(::IDPackedV) = false

isreal(::IDPackedV{T}) where {T} = T <: Real

ndims(::IDPackedV) = 2

size(A::IDPackedV) = (size(A.T,1), sum(size(A.T)))
size(A::IDPackedV, dim::Integer) =
  (dim == 1 ? size(A.T,1) : (dim == 2 ? sum(size(A.T)) : 1))

## BLAS/LAPACK multiplication routines
  ### left-multiplication
  function mul!!(y::AbstractVector{T}, A::IDPackedV{T}, x::AbstractVector{T}) where T<:BlasFloat
    k, n = size(A)
    lmul!(A[:P]', x)
    copyto!(y, view(x,1:k))
    BLAS.gemv!('N', one(T), A[:T], view(x,k+1:n), one(T), y)
  end  # overwrites x
  function mul!!(C::AbstractMatrix{T}, A::IDPackedV{T}, B::AbstractMatrix{T}) where T<:BlasFloat
    k, n = size(A)
    lmul!(A[:P]', B)
    copyto!(C, view(B,1:k,:))
    BLAS.gemm!('N', 'N', one(T), A[:T], view(B,k+1:n,:), one(T), C)
  end  # overwrites B
  mul!(C::AbstractVector{T}, A::IDPackedV{T}, B::AbstractVector{T}) where {T} =
    mul!!(C, A, copy(B))
  mul!(C::AbstractMatrix{T}, A::IDPackedV{T}, B::AbstractMatrix{T}) where {T} =
    mul!!(C, A, copy(B))

  for Adj in (:Transpose, :Adjoint)
    @eval begin
      function mul!(C::AbstractMatrix{T}, A::IDPackedV{T}, Bc::$Adj{T,<:AbstractMatrix{T}}) where T<:BlasFloat
        k, n = size(A)
        tmp = $Adj(A[:P]) * Bc
        copyto!(C, view(tmp,1:k,:))
        BLAS.gemm!('N', 'N', one(T), A[:T], view(tmp,k+1:n,:), one(T), C)
      end
      function mul!(C::AbstractVector{T}, Ac::$Adj{T,IDPackedV{T}}, B::AbstractVector{T}) where T
        A = parent(Ac)
        k, n = size(A)
        copyto!(view(C,1:k,:), B)
        mul!(view(C,k+1:n,:), $Adj(A[:T]), B)
        lmul!(A[:P], C)
      end
      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,IDPackedV{T}}, B::AbstractMatrix{T}) where T
        A = parent(Ac)
        k, n = size(A)
        copyto!(view(C,1:k,:), B)
        mul!(view(C,k+1:n,:), $Adj(A[:T]), B)
        lmul!(A[:P], C)
      end
    end
  end

  for (Adj, adj!) in ((:Transpose, :transpose!), (:Adjoint,:adjoint!))
    @eval begin
      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,IDPackedV{T}}, Bc::$Adj{T,<:AbstractMatrix{T}}) where T
        A = parent(Ac)
        B = parent(Bc)
        k, n = size(A)
        $adj!(view(C,1:k,:), B)
        mul!(view(C,k+1:n,:), $Adj(A[:T]), $Adj(B))
        lmul!(A[:P], C)
      end
    end
  end

  ### right-multiplication

  function mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, B::IDPackedV{T}) where T
    k, n = size(B)
    copyto!(view(C,:,1:k), A)
    mul!(view(C,:,k+1:n), A, B[:T])
    rmul!(C, B[:P]')
  end

  for (Adj, trans) in ((:Adjoint, 'C'), (:Transpose, 'T'))
    @eval begin
      function mul!!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, Bc::$Adj{T,IDPackedV{T}}) where T<:BlasFloat
        B = parent(Bc)
        k, n = size(B)
        rmul!(A, B[:P])
        copyto!(C, view(A,:,1:k))
        BLAS.gemm!('N', $trans, one(T), view(A,:,k+1:n), B[:T], one(T), C)
      end  # overwrites A
      mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, Bc::$Adj{T,IDPackedV{T}}) where {T} =
        mul!!(C, copy(A), Bc)
    end
  end

  for (Adj, adj!) in ((:Transpose, :transpose!), (:Adjoint,:adjoint!))
    @eval begin
      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,<:AbstractMatrix{T}}, B::IDPackedV{T}) where T
        A = parent(Ac)
        k, n = size(B)
        $adj!(view(C,:,1:k), A)
        mul!(view(C,:,k+1:n), $Adj(A), B[:T])
        rmul!(C, $Adj(B[:P]))
      end
    end
  end

  for (Adj, trans) in ((:Adjoint, 'C'), (:Transpose, 'T'))
    @eval begin
      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,<:AbstractMatrix{T}}, Bc::$Adj{T,IDPackedV{T}}) where T<:BlasFloat
        B = parent(Bc)
        k, n = size(B)
        tmp = Ac * B[:P]
        copyto!(C, view(tmp,:,1:k))
        BLAS.gemm!('N', $trans, one(T), view(tmp,:,k+1:n), B[:T], one(T), C)
      end
    end
  end

# ID

mutable struct ID{S} <: Factorization{S}
  sk::Vector{Int}
  rd::Vector{Int}
  C::Matrix{S}
  T::Matrix{S}
  ID{S}(sk::Vector{Int}, rd::Vector{Int}, C::Matrix{S}, T::Matrix{S}) where S =
    new{S}(sk, rd, C, T)
end
ID{S}(sk::Vector{Int}, rd::Vector{Int}, C::AbstractMatrix{S}, T::AbstractMatrix{S}) where S =
  ID{S}(sk, rd, convert(Matrix{S}, C), convert(Matrix{S}, T))

function ID(trans::Symbol, A::AbstractMatOrLinOp{T}, V::IDPackedV{T}) where T
  chktrans(trans)
  ID{T}(V.sk, V.rd, getcols(trans, A, V.sk), V.T)
end
ID(trans::Symbol, A::AbstractMatOrLinOp, sk, rd, T) =
  ID(trans, A, IDPackedV(sk, rd, T))
ID(A::AbstractMatOrLinOp, args...) = ID(:n, A, args...)
ID(sk::Vector{Int}, rd::Vector{Int}, C::Matrix{S}, T::Matrix{S}) where S =
  ID{S}(sk, rd, C, T)
ID(A, args...) = ID(LinOp(A), args...)

conj!(A::ID) = ID(A.sk, A.rd, conj!(A.C), conj!(A.T))
conj(A::ID) = ID(A.sk, A.rd, conj(A.C), conj(A.T))

convert(::Type{ID{T}}, A::ID) where {T} =
  ID(A.sk, A.rd, convert(Array{T}, A.C), convert(Array{T}, A.T))
convert(::Factorization{T}, A::ID) where {T} = convert(ID{T}, A)
convert(::Type{Array}, A::ID) = Matrix(A)
convert(::Type{Array{T}}, A::ID) where {T} = convert(Array{T}, Matrix(A))

copy(A::ID) = ID(copy(A.sk), copy(A.rd), copy(A.C), copy(A.T))

Matrix(A::ID) = A[:C]*A[:V]

  adjoint(A::ID) = Adjoint(A)
  transpose(A::ID) = Transpose(A)

function getindex(A::ID{T}, d::Symbol) where T
  if     d == :C   return A.C
  elseif d == :P   return ColumnPermutation(A[:p])
  elseif d == :T   return A.T
  elseif d == :V   return IDPackedV(A.sk, A.rd, A.T)
  elseif d == :k   return length(A.sk)
  elseif d == :p   return [A.sk; A.rd]
  elseif d == :rd  return A.rd
  elseif d == :sk  return A.sk
  else             throw(KeyError(d))
  end
end

ishermitian(::ID) = false
issym(::ID) = false

isreal(::ID{T}) where {T} = T <: Real

ndims(::ID) = 2

size(A::ID) = (size(A.C,1), sum(size(A.T)))
size(A::ID, dim::Integer) =
  (dim == 1 ? size(A.C,1) : (dim == 2 ? sum(size(A.T)) : 1))

## BLAS/LAPACK multiplication routines

### left-multiplication

function mul!!(y::AbstractVector{T}, A::ID{T}, x::AbstractVector{T}) where T
  tmp = Array{T}(undef, A[:k])
  mul!!(tmp, A[:V], x)
  mul!(y, A[:C], tmp)
end  # overwrites x
function mul!!(C::AbstractMatrix{T}, A::ID{T}, B::AbstractMatrix{T}) where T
  tmp = Array{T}(undef, A[:k], size(B,2))
  mul!!(tmp, A[:V], B)
  mul!(C, A[:C], tmp)
end  # overwrites B
mul!(C::AbstractVector{T}, A::ID{T}, B::AbstractVector{T}) where {T} =
  mul!!(C, A, copy(B))
mul!(C::AbstractMatrix{T}, A::ID{T}, B::AbstractMatrix{T}) where {T} =
  mul!!(C, A, copy(B))

  for Adj in (:Transpose, :Adjoint)
    @eval begin
      function mul!(C::AbstractMatrix{T}, A::ID{T}, Bc::$Adj{T,<:AbstractMatrix{T}}) where T
        tmp = A[:V] * Bc
        mul!(C, A[:C], tmp)
      end
      function mul!(C::AbstractVector{T}, Ac::$Adj{T,ID{T}}, B::AbstractVector{T}) where T
        A = parent(Ac)
        tmp = $Adj(A[:C]) * B
        mul!(C, $Adj(A[:V]), tmp)
      end
      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,ID{T}}, B::AbstractMatrix{T}) where T
        A = parent(Ac)
        tmp = $Adj(A[:C]) * B
        mul!(C, $Adj(A[:V]), tmp)
      end
      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,ID{T}}, Bc::$Adj{T,<:AbstractMatrix{T}}) where T
        A = parent(Ac)
        tmp = $Adj(A[:C]) * Bc
        mul!(C, $Adj(A[:V]), tmp)
      end
        ### right-multiplication
      function mul!!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, Bc::$Adj{T,ID{T}}) where T
        B = parent(Bc)
        tmp = Array{T}(undef, size(A,1), B[:k])
        mul!!(tmp, A, $Adj(B[:V]))
        mul!(C, tmp, $Adj(B[:C]))
      end  # overwrites A
      mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, Bc::$Adj{T,ID{T}}) where {T} =
        mul!!(C, copy(A), Bc)
      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,<:AbstractMatrix{T}}, B::ID{T}) where T
        tmp = Ac * B[:C]
        mul!(C, tmp, B[:V])
      end
      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,<:AbstractMatrix{T}}, Bc::$Adj{T,ID{T}}) where T
        B = parent(Bc)
        tmp = Ac * $Adj(B[:V])
        mul!(C, tmp, $Adj(B[:C]))
      end
    end
  end

  ### right-multiplication

  mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, B::ID{T}) where {T} =
    mul!(C, A*B[:C], B[:V])

  # standard operations

  ## left-multiplication
  for t in (:IDPackedV, :ID)
    @eval begin
      function *(A::$t{TA}, B::AbstractVector{TB}) where {TA,TB}
        T = promote_type(TA, TB)
        AT = convert($t{T}, A)
        BT = (T == TB ? B : convert(Array{T}, B))
        CT = Array{T}(undef, size(A,1))
        mul!(CT, AT, BT)
      end
      function *(A::$t{TA}, B::AbstractMatrix{TB}) where {TA,TB}
        T = promote_type(TA, TB)
        AT = convert($t{T}, A)
        BT = (T == TB ? B : convert(Array{T}, B))
        CT = Array{T}(undef, size(A,1), size(B,2))
        mul!(CT, AT, BT)
      end
      function *(A::AbstractMatrix{TA}, B::$t{TB}) where {TA,TB}
        T = promote_type(TA, TB)
        AT = (T == TA ? A : convert(Array{T}, A))
        BT = convert($t{T}, B)
        CT = Array{T}(undef, size(A,1), size(B,2))
        mul!(CT, AT, BT)
      end
    end
    for Adj in (:Transpose, :Adjoint)
      @eval begin
        function *(Ac::$Adj{TA,$t{TA}}, B::AbstractVector{TB}) where {TA,TB}
          A = parent(Ac)
          T = promote_type(TA, TB)
          AT = convert($t{T}, A)
          BT = (T == TB ? B : convert(Array{T}, B))
          CT = Array{T}(undef, size(Ac,1))
          mul!(CT, $Adj(AT), BT)
        end
        function *(A::$t{TA}, Bc::$Adj{TB,<:AbstractMatrix{TB}}) where {TA,TB}
          B = parent(Bc)
          T = promote_type(TA, TB)
          AT = convert($t{T}, A)
          BT = (T == TB ? B : convert(Array{T}, B))
          CT = Array{T}(undef, size(A,1), size(Bc,2))
          mul!(CT, AT, $Adj(BT))
        end
        function *(Ac::$Adj{TA,$t{TA}}, B::AbstractMatrix{TB}) where {TA,TB}
          A = parent(Ac)
          T = promote_type(TA, TB)
          AT = convert($t{T}, A)
          BT = (T == TB ? B : convert(Array{T}, B))
          CT = Array{T}(undef, size(Ac,1), size(B,2))
          mul!(CT, $Adj(AT), BT)
        end
        function *(Ac::$Adj{TA,$t{TA}}, Bc::$Adj{TB,<:AbstractMatrix{TB}}) where {TA,TB}
          A = parent(Ac)
          B = parent(Bc)
          T = promote_type(TA, TB)
          AT = convert($t{T}, A)
          BT = (T == TB ? B : convert(Array{T}, B))
          CT = Array{T}(undef, size(Ac,1), size(Bc,2))
          mul!(CT, $Adj(AT), $Adj(BT))
        end
        function *(A::AbstractMatrix{TA}, Bc::$Adj{TB,$t{TB}}) where {TA,TB}
          B = parent(Bc)
          T = promote_type(TA, TB)
          AT = (T == TA ? A : convert(Array{T}, A))
          BT = convert($t{T}, B)
          CT = Array{T}(undef, size(A,1), size(Bc,2))
          mul!(CT, AT, $Adj(BT))
        end
        function *(Ac::$Adj{TA,<:AbstractMatrix{TA}}, B::$t{TB}) where {TA,TB}
          A = parent(Ac)
          T = promote_type(TA, TB)
          AT = (T == TA ? A : convert(Array{T}, A))
          BT = convert($t{T}, B)
          CT = Array{T}(undef, size(Ac,1), size(B,2))
          mul!(CT, $Adj(AT), BT)
        end
        function *(Ac::$Adj{TA,<:AbstractMatrix{TA}}, Bc::$Adj{TB,$t{TB}}) where {TA,TB}
          A = parent(Ac)
          B = parent(Bc)
          T = promote_type(TA, TB)
          AT = (T == TA ? A : convert(Array{T}, A))
          BT = convert($t{T}, B)
          CT = Array{T}(undef, size(Ac,1), size(Bc,2))
          mul!(CT, $Adj(AT), $Adj(BT))
        end
      end
    end
  end

# factorization routines

for sfx in ("", "!")
  f = Symbol("idfact", sfx)
  g = Symbol("pqrfact", sfx)
  h = Symbol("id", sfx)
  @eval begin
    function $f(
        trans::Symbol, A::AbstractMatOrLinOp{T}, opts::LRAOptions=LRAOptions(T);
        args...) where T
      opts = copy(opts; args...)
      opts.pqrfact_retval = "t"
      chkopts!(opts, A)
      if opts.sketch == :none
        F = $g(trans, A, opts)
      else
        F = sketchfact(:left, trans, A, opts)
      end
      k = F[:k]
      IDPackedV(F.p[1:k], F.p[k+1:end], get(F.T))
    end
    $f(trans::Symbol, A, args...; kwargs...) =
      $f(trans, LinOp(A), args...; kwargs...)
    $f(A, args...; kwargs...) = $f(:n, A, args...; kwargs...)

    function $h(trans::Symbol, A, args...; kwargs...)
      V = $f(trans, A, args...; kwargs...)
      V.sk, V.rd, V.T
    end
    $h(A, args...; kwargs...) = $h(:n, A, args...; kwargs...)
  end
end


text:
import Base.length

using Distributions:Categorical
using StatsFuns: logsumexp
using Scruff.SFuncs: Cat

export
    Particles,
    effective_sample_size,
    normalize_weights,
    log_prob_evidence,
    resample,
    probability,
    marginal

"Sample = Dict{Symbol, Any}"
Sample = Dict{Symbol, Any}

"""
    struct Particles

A structure of particles containing a vector of `Samples` and of log_weights.
"""
struct Particles
    samples::Vector{Sample}
    log_weights::Vector{Float64}
end

"""
    probability(parts::Particles, predicate::Sample -> Bool)::Float64

    Returns the probability that the predicate is satisfied
"""
function probability(parts::Particles, predicate::Function)::Float64
    sum = 0.0
    tot = 0.0
    for i in 1:length(parts.samples)
        w = exp(parts.log_weights[i])
        s = parts.samples[i]
        if predicate(s)
            sum += w
        end
        tot += w
    end
    return sum / tot
end

"""
    marginal(parts::Particles, x::Symbol)::Cat

Returns a Cat representing the marginal distribution over the given symbol according to parts
"""
function marginal(parts::Particles, x::Symbol)::Cat
    d = Dict{Any, Float64}()
    xs = [s[x] for s in parts.samples]
    lws = normalize_weights(parts.log_weights)
    for i in 1:length(parts.samples)
        d[xs[i]] = get(d, xs[i], 0.0) + exp(lws[i])
    end
    ks = Vector{Any}()
    ps = Vector{Float64}()
    for (k,p) in d
        push!(ks, k)
        push!(ps, p)
    end
    return Cat(ks, ps)
end

"""
    expected_sample_size(log_weights::Vector{Float64})

Effective sample size
"""
function effective_sample_size(lws::Vector{Float64})
    log_ess = (2 * logsumexp(lws)) - logsumexp(2 .* lws)
    return exp(log_ess)
end

"""
    normalize_weights(log_weights::Vector{Float64})

Normalize weights
"""
function normalize_weights(lws::Vector{Float64})
    return lws .- logsumexp(lws)
end

function log_prob_evidence(lws::Vector{Float64})
    return logsumexp(lws) - log(length(lws))
end

function resample(ps::Particles, target_num_particles::Int = length(ps.samples))
    lnws = normalize_weights(ps.log_weights)
    weights = exp.(lnws)
    selections = rand(Categorical(weights/sum(weights)), target_num_particles)
    samples = map(selections) do ind
        ps.samples[ind]
    end
    return Particles(samples, zeros(target_num_particles))
end


text:
using Documenter
makedocs(;
    sitename="Mizuto Kadowaki"
)

deploydocs(;
    devbranch="main",
    target="build",
    repo="github.com/MizutoKadowaki0312/Documenter_website",
    versions=nothing
)


text:
let
    basepath = joinpath(pwd(),"foo")
    subdirs = ["bar","baz"]
    files = ["hello.a","goodbye.a","hello.b","goodbye.b"]
    for d in subdirs
        for f in files
            mkpath(joinpath(basepath,d,f))
        end
    end
    for ext in [".a",".b"]
        for kw in ["hello","goodbye"]
            @test all(p->findfirst(ext,p) != nothing, get_files_matching(basepath,ext))
            @test all(p->findfirst(kw,p) != nothing, get_files_matching(basepath,ext,[kw]))
        end
    end
    rm(basepath;recursive=true)
end
let
    logfile = "tmp.log"
    errfile = "tmp.err"
    logmsg = "hello world"
    errmsg = "error world"
    redirect_to_files(logfile,errfile) do
        println(logmsg)
        println(stderr,errmsg)
    end
    @test readline(logfile) == logmsg
    @test readline(errfile) == errmsg
    rm(logfile)
    rm(errfile)
end


text:
export BERT
using Knet.Ops21
using Knet.Layer21

struct BERT

end

function BERT(;
              vocab_size = 30522,
              hidden_size = 768,
              num_hidden_layers = 12,
              num_attention_layers = 12,
              intermediate_size = 3072,
              hidden_act = gelu,
              hidden_dropout_prob = 0.1,
              attention_probs_dropout_prob = 0.1,
              max_position_embeddings = 512,
              type_vocab_size = 2,
              initializer_range = 0.02,
              layer_norm_eps = 1e-12,
              gradient_checkpointing = false,
              )

end


text:
using SnowyOwl
using DataFrames
using LinearAlgebra
using SparseArrays
using Test

tests = [
    "object",
    # "io"
    "filter",
    "moments",
]

@testset "SnowyOwl.jl" begin
    for t in tests
        include("$(t).jl")
    end
end


text:
@turing_testset "gibbs conditionals" begin
    Random.seed!(100)

    @turing_testset "gdemo" begin
        # We consider the model
        # ```math
        # s ~ InverseGamma(2, 3)
        # m ~ Normal(0, √s)
        # xᵢ ~ Normal(m, √s), i = 1, …, N,
        # ```
        # with ``N = 2`` observations ``x₁ = 1.5`` and ``x₂ = 2``.

        # The conditionals and posterior can be formulated in terms of the following statistics:
        N = 2
        x_mean = 1.75 # sample mean ``∑ xᵢ / N``
        x_var = 0.0625 # sample variance ``∑ (xᵢ - x_bar)^2 / N``
        m_n = 3.5 / 3 # ``∑ xᵢ / (N + 1)``

        # Conditional distribution
        # ```math
        # m | s, x ~ Normal(m_n, sqrt(s / (N + 1)))
        # ```
        cond_m = let N=N, m_n=m_n
            c -> Normal(m_n, sqrt(c.s / (N + 1)))
        end

        # Conditional distribution
        # ```math
        # s | m, x ~ InverseGamma(2 + (N + 1) / 2, 3 + (m^2 + ∑ (xᵢ - m)^2) / 2) =
        #            InverseGamma(2 + (N + 1) / 2, 3 + m^2 / 2 + N / 2 * (x_var + (x_mean - m)^2))
        # ```
        cond_s = let N=N, x_mean=x_mean, x_var=x_var
            c -> InverseGamma(2 + (N + 1) / 2, 3 + c.m^2 / 2 + N / 2 * (x_var + (x_mean - c.m)^2))
        end

        # Three Gibbs samplers:
        # one for each variable fixed to the posterior mean
        s_posterior_mean = 49/24
        sampler1 = Gibbs(
            GibbsConditional(:m, cond_m),
            GibbsConditional(:s, _ -> Normal(s_posterior_mean, 0)),
        )
        chain = sample(gdemo_default, sampler1, 10_000)
        cond_m_mean = mean(cond_m((s = s_posterior_mean,)))
        check_numerical(chain, [:m, :s], [cond_m_mean, s_posterior_mean])
        @test all(==(s_posterior_mean), chain[:s][2:end])

        m_posterior_mean = 7/6
        sampler2 = Gibbs(
            GibbsConditional(:m, _ -> Normal(m_posterior_mean, 0)),
            GibbsConditional(:s, cond_s),
        )
        chain = sample(gdemo_default, sampler2, 10_000)
        cond_s_mean = mean(cond_s((m = m_posterior_mean,)))
        check_numerical(chain, [:m, :s], [m_posterior_mean, cond_s_mean])
        @test all(==(m_posterior_mean), chain[:m][2:end])

        # and one for both using the conditional
        sampler3 = Gibbs(GibbsConditional(:m, cond_m), GibbsConditional(:s, cond_s))
        chain = sample(gdemo_default, sampler3, 10_000)
        check_gdemo(chain)
    end

    @turing_testset "GMM" begin
        # We consider the model
        # ```math
        # μₖ ~ Normal(m, σ_μ), k = 1, …, K,
        # zᵢ ~ Categorical(π), i = 1, …, N,
        # xᵢ ~ Normal(μ_{zᵢ}, σₓ), i = 1, …, N,
        # ```
        # with ``K = 2`` clusters, ``N = 20`` observations, and the following parameters:
        K = 2 # number of clusters
        π = fill(1/K, K) # uniform cluster weights
        m = 0.5 # prior mean of μₖ
        σ²_μ = 4.0 # prior variance of μₖ
        σ²_x = 0.01 # observation variance
        N = 20  # number of observations

        # We generate data
        μ_data = rand(Normal(m, sqrt(σ²_μ)), K)
        z_data = rand(Categorical(π), N)
        x_data = rand(MvNormal(μ_data[z_data], σ²_x * I))

        @model function mixture(x)
            μ ~ $(MvNormal(fill(m, K), σ²_μ * I))
            z ~ $(filldist(Categorical(π), N))
            x ~ MvNormal(μ[z], $(σ²_x * I))
            return x
        end
        model = mixture(x_data)

        # Conditional distribution ``z | μ, x``
        # see http://www.cs.columbia.edu/~blei/fogm/2015F/notes/mixtures-and-gibbs.pdf
        cond_z = let x=x_data, log_π=log.(π), σ_x=sqrt(σ²_x)
            c -> begin
                dists = map(x) do xi
                    logp = log_π .+ logpdf.(Normal.(c.μ, σ_x), xi)
                    return Categorical(StatsFuns.softmax!(logp))
                end
                return arraydist(dists)
            end
        end

        # Conditional distribution ``μ | z, x``
        # see http://www.cs.columbia.edu/~blei/fogm/2015F/notes/mixtures-and-gibbs.pdf
        cond_μ = let K=K, x_data=x_data, inv_σ²_μ=inv(σ²_μ), inv_σ²_x=inv(σ²_x)
            c -> begin
                # Convert cluster assignments to one-hot encodings
                z_onehot = c.z .== (1:K)'

                # Count number of observations in each cluster
                n = vec(sum(z_onehot; dims=1))

                # Compute mean and variance of the conditional distribution
                μ_var = @. inv(inv_σ²_x * n + inv_σ²_μ)
                μ_mean = (z_onehot' * x_data) .* inv_σ²_x .* μ_var

                return MvNormal(μ_mean, Diagonal(μ_var))
            end
        end

        estimate(chain, var) = dropdims(mean(Array(group(chain, var)), dims=1), dims=1)
        function estimatez(chain, var, range)
            z = Int.(Array(group(chain, var)))
            return map(i -> findmax(counts(z[:,i], range))[2], 1:size(z,2))
        end

        lμ_data, uμ_data = extrema(μ_data)

        # Compare three Gibbs samplers
        sampler1 = Gibbs(GibbsConditional(:z, cond_z), GibbsConditional(:μ, cond_μ))
        sampler2 = Gibbs(GibbsConditional(:z, cond_z), MH(:μ))
        sampler3 = Gibbs(GibbsConditional(:z, cond_z), HMC(0.01, 7, :μ))
        for sampler in (sampler1, sampler2, sampler3)
            chain = sample(model, sampler, 10_000)

            μ_hat = estimate(chain, :μ)
            lμ_hat, uμ_hat = extrema(μ_hat)
            @test isapprox([lμ_data, uμ_data], [lμ_hat, uμ_hat], atol=0.1)

            z_hat = estimatez(chain, :z, 1:2)
            ari, _, _, _ = randindex(z_data, Int.(z_hat))
            @test isapprox(ari, 1, atol=0.1)
        end
    end
end


text:
# This program solves Problem in Section 5.1.1 in the following paper
# Rognes, M. E., Ham, D. A., Cotter, C. J., and McRae, A. T. T.
# Automating the solution of PDEs on the sphere and other manifolds in FEniCS 1.2,
# Geosci. Model Dev., 6, 2099–2119, https://doi.org/10.5194/gmd-6-2099-2013, 2013.

function solve_darcy(model,order,degree,ls=BackslashSolver())
  g(x) = sin(0.5*π*x[2]) # x[1]*x[2]*x[3]
  RT=ReferenceFE(raviart_thomas,Float64,order)
  DG=ReferenceFE(lagrangian,Float64,order)
  V = FESpace(model,RT; conformity=:Hdiv)
  Q = FESpace(model,DG; conformity=:L2)
  T = FESpace(model,DG; conformity=:L2)

  U = TrialFESpace(V)
  P = TrialFESpace(Q)
  R = TrialFESpace(T)

  Y = MultiFieldFESpace([V, Q, T])
  X = MultiFieldFESpace([U, P, R])

  trian = Triangulation(model)
  dΩ = Measure(trian,degree)
  dω = Measure(trian,degree,ReferenceDomain())

  aref((u,p,r),(v,q,t)) = ∫( v⋅u + q*r + t*p )*dΩ +  ∫(q*DIV(u)+DIV(v)*p)*dω
  aphy((u,p,r),(v,q,t)) = ∫( v⋅u + q*r + t*p + q*(∇⋅u) + (∇⋅v)*p)*dΩ
  b((v,q,t)) = ∫(q*g)*dΩ

  op = AffineFEOperator(aref,b,X,Y)
  fels = LinearFESolver(ls)
  xh = solve(fels,op)
  uh, ph, rh = xh
  e = g-rh
  sqrt(sum(∫(e*e)dΩ))
end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
using FunctionalPipes
using Test

@testset "FunctionalPipes.jl" begin

    add(n) = x -> x + n
    mul(n) = x -> x * n
    zero() = 0

    # init value provided, single function
    @test pipe(1, add(1)) == 2

    # init value provided, multiple functions
    @test pipe(1, add(1), mul(2), add(3)) == 7

    # no init value, single function
    @test pipe(zero, add(1)) == 1

    # no init value, multiple functions
    @test pipe(zero, add(1), mul(2)) == 2

    # nothing as init value is the same as no init value
    @test pipe(nothing, zero, add(1)) == 1

    # dispatch with side effect
    @test pipe(
        zero,
        add(1),
        dispatch(println, "good job: ", 🔥),
        mul(2)
    ) == 2

    @test pipe(
        zero,
        add(3),
        dispatch((x,y,z)->x*y*z, 2, 🔥, 4),
        dispatch(println, "variation 1: ", 🔥)
    ) == 24

    @test pipe(
        zero,
        add(3),
        dispatch((x,y,z)->x*y*z, 🔥, 2, 4),
        dispatch(println, "variation 2: ", 🔥)
    ) == 24

    @test pipe(
        zero,
        add(3),
        dispatch((x,y,z)->x*y*z, 2, 4, 🔥),
        dispatch(println, "variation 3: ", 🔥)
    ) == 24

    @test pipe(
        1:10,
        collect,
        dispatch(filter, isodd, 🔥),
        dispatch(println, "values: ", 🔥),
        sum,
        dispatch(println, "total: ", 🔥)
    ) == 25

end


text:
#!/usr/local/julia/bin/julia

# @license BSD-3 https://opensource.org/licenses/BSD-3-Clause
# Copyright (c) 2020, Tim Redick
# All rights reserved.

using Pkg

function has_manifest_file()
    return isfile(joinpath(pwd(),"Manifest.toml")) || isfile(joinpath(pwd(),"JuliaManifest.toml"))
end

function has_project_file()
    return isfile(joinpath(pwd(),"Project.toml")) || isfile(joinpath(pwd(),"JuliaProject.toml"))
end

function is_package()
    return has_project_file() && !has_manifest_file()
end

function is_app()
    return has_project_file() && has_manifest_file()
end

if is_app()
    Pkg.activate(pwd())
    Pkg.instantiate()
elseif is_package()
    Pkg.develop(path=pwd())
end

max_stars_repo_name: Tuebel/BayesionMethodsForHackers

max_issues_repo_name: Tuebel/BayesionMethodsForHackers

max_forks_repo_name: Tuebel/BayesianMethodsForHackers


text:
module ReducedNetwork
using DifferentialEquations, LightGraphs, SimpleWeightedGraphs, Statistics,
    Plots, GraphPlot, Compose, Cairo, Fontconfig, Printf # for plotting
using SpikingNeuralNets.GraphUtils: mesh_graph

# Effective transfer function
x(S1, S2, Inet, I, Inoise, t, p) = p.JNex*S1 + p.JNin*S2 + Inet + p.Iext + I*(t >= 0.0) + Inoise
function H(x1, x2, p)
    v = p.a*x1 - p.fA(x2) - p.b
    return v / (1 - exp(p.d * v))
end

"""
    Inetwork(S::AbstractVector{T}, g::AbstractGraph)::Vector{T} where {T <: Number}

For each vertex in `g` calculate a weighted average of S from all inneighbors of `v`.
"""
function Inetwork(S::AbstractVector{T}, g::AbstractGraph)::Vector{T} where {T <: Number}
    I = zeros(T, size(S))
    @inbounds for v in vertices(g)
        in = inneighbors(g, v)
        if !iszero(length(in))
            I[v] = sum(weights(g)[in,v] .* S[in])
        end
    end
    return I
end

"""
    ReducedModel(du, u, p, t, W)

Calculate the change mean gating variables `du` given the current gating
variables `u`, named parameters `p`, time `t`, and noise currents `W`.
"""
function ReducedModel(du, u, p, t, W)
    Sa, Sb = eachcol(u)
    Inoisea, Inoiseb = eachcol(W)

    Ineta = Inetwork(p.Jex*Sa + p.Jin*Sb, p.g)          # network inputs
    Inetb = Inetwork(p.Jex*Sb + p.Jin*Sa, p.g)
    xa = x.(Sa, Sb, Ineta, p.Ia, Inoisea, t, Ref(p))    # intermediate variables
    xb = x.(Sb, Sa, Inetb, p.Ib, Inoiseb, t, Ref(p))
    du[:,1] .= .-(Sa/p.τs) .+ (1 .- Sa) .* p.γ .* H.(xa, xb, Ref(p))
    du[:,2] .= .-(Sb/p.τs) .+ (1 .- Sb) .* p.γ .* H.(xb, xa, Ref(p))
end

"""
    run(g::AbstractGraph; c=fill(0, nv(g)), μ0=fill(30.0, nv(g)), S0=zeros(nv(g), 2),
        Inoise0=zeros(nv(g), 2), tstart=-0.5, tstop=1.0, reltol=1e-15,
        abstol=1e-15, dt=0.1e-3, Jex=1.0, Jin=0.0, Jext=0.2243e-3, JNex=0.1561,
        JNin=-0.0264, JAex=9.9026e-4, JAin=-6.5177e-5, τs=100e-3, τAMPA=2e-3,
        γ=0.641, Iext=0.35, σnoise=0.007, display=true, color1=:red, color2=:blue,
        kwargs...)

Run a reduced 2-variable model of neural populations during a random dot motion task.
Return a tuple (g, t, Sa, Sb, Inoisea, Inoiseb, Ineta, Inetb, xa, xb, ra, rb).

`g`: a graph of neural populations, where each vertex represents a population of
     neurons, and each edge represents a connection between neural populations.
     Weighted graphs (e.g. `SimpleWeightedDiGraph`) a.
`c`: the coherence of external input to each node.
`μ0`: the firing rate (in Hz) of external input to each node,
`S0`: the initial average NMDA gating variables for each node,
`Inoise0`: are the initial noise currents for each node,
`tstart`: the starting time of the simulation (negative times having no stimulus)
`tstop`: is the stopping time of the simulation
`dt`: is the timestep of the simulation.
`Jex`: the weight of recurrent input from excitatory sub-populations
`Jin`: the weight of recurrent input from inhibitory sub-populations
`Jext`: the weight of external inputs
`JNex`: the excitatory (> 0) NMDA weight
`JNin`: the inhibitory (< 0) NMDA weight
`JAex`: the excitatory (> 0) AMPA weight
`JAin`: the excitatory (> 0) AMPA weight
`τs`: the time constant of NMDA channels
`τAMPA`: the time constant of AMPA channels
`γ`: a factor applied to firing rate to update gating variables
`Iext`: baseline external input current (nV)
`σnoise`: standard deviation of noise current (nV).
`display`: if true, draw a plot of gating variables and firing rates over time
`color1`: the line color of the correct direction in the output plot
`color2`: the line color of the incorrect direction in the output plot
Additional keyword arguments are passed to DifferentialEquations.solve.
"""
function run(g::AbstractGraph; c=fill(0, nv(g)), μ0=fill(30.0, nv(g)), S0=zeros(nv(g), 2),
             Inoise0=zeros(nv(g), 2), tstart=-0.5, tstop=1.0, reltol=1e-15,
             abstol=1e-15, dt=0.1e-3, Jex=1.0, Jin=0.0, Jext=0.2243e-3, JNex=0.1561,
             JNin=-0.0264, JAex=9.9026e-4, JAin=-6.5177e-5, τs=100e-3, τAMPA=2e-3,
             γ=0.641, Iext=0.35, σnoise=0.007, display=true, color1=:red, color2=:blue,
             kwargs...)
    # set up parameters
    Ia = @. Jext*μ0*(1 + c/100.0)
    Ib = @. Jext*μ0*(1 - c/100.0)
    a = 239400JAex + 270
    b = 97000JAex + 108
    d = 30JAex - 0.154
    fA(x) = JAin*(276x - 106) * (x >= 0.4)
    params = (g=g, a=a, b=b, c=c, d=d, fA=fA, μ0=μ0, Ia=Ia, Ib=Ib, Jex=Jex, Jin=Jin,
              Jext=Jext, JNex=JNex, JNin=JNin, JAex=JAex, JAin=JAin,
              τs=τs, τAMPA=τAMPA, γ=γ, Iext=Iext)

    # noise process: τAMPA dI(t)/dt = -I(t) + W(t)*σnoise*√τAMPA
    W = OrnsteinUhlenbeckProcess(τAMPA, 0.0, σnoise*√τAMPA, tstart, Inoise0)
    prob = RODEProblem(ReducedModel, S0, (tstart, tstop), params, noise=W)
    sol = solve(prob, reltol=reltol, abstol=abstol, dt=dt, kwargs...)

    if display
        # Reconstruct the variables
        Sa, Sb = sol[:,1,:], sol[:,2,:]
        Inoisea, Inoiseb = sol.W[:,1,1:length(sol.t)], sol.W[:,2,1:length(sol.t)]
        Ineta = mapslices((s)->Inetwork(s, g), Jex*Sa + Jin*Sb; dims=1)
        Inetb = mapslices((s)->Inetwork(s, g), Jex*Sb + Jin*Sa; dims=1)

        # reconstruct the firing rates
        xa, xb = similar(Sa), similar(Sb)
        @views for v in vertices(g)
            xa[v,:] .= x.(Sa[v,:], Sb[v,:], Ineta[v,:], Ia[v], Inoisea[v,:], sol.t, Ref(params))
            xb[v,:] .= x.(Sb[v,:], Sa[v,:], Inetb[v,:], Ib[v], Inoiseb[v,:], sol.t, Ref(params))
        end
        ra = H.(xa, xb, Ref(params))
        rb = H.(xb, xa, Ref(params))

        Base.display(plot(plot(sol.t, [Sa' Sb'], title="Gating Variables"),
                     plot(sol.t, [ra' rb'], title="Firing Rates (Hz)"),
                     leg=false, linecolor=[fill(:red, 1, nv(g)) fill(:blue, 1, nv(g))], layout=(2, 1)))
        return (g, sol.t, Sa, Sb, Inoisea, Inoiseb, Ineta, Inetb, xa, xb, ra, rb)
    else
        return g, sol
    end
end

"""
    run(g::AbstractGraph, inputs; c=0.0, μ0=30.0, kwargs...)

Run a reduced 2-variable model of neural populations during a random dot motion task,
where `inputs` contains all of the vertices of `g` receiving external input, and
`c` and `μ0` are the coherence and strength of external input. Return a tuple
(g, t, Sa, Sb, Inoisea, Inoiseb, Ineta, Inetb, xa, xb, ra, rb).
"""
function run(g, inputs; c=0.0, μ0=30.0, kwargs...)
    C = fill(c, nv(g))
    μ = zeros(nv(g))
    μ[inputs] .= μ0
    return run(g; c=C, μ0=μ, kwargs...)
end

"""
    run(dims::Tuple{Vararg{<:Integer, N}}; inputs=ifelse(N==1, 1, 1:first(dims)), weight=0.015, c=0, μ0=30.0, kwargs...) where N

Run a reduced 2-variable model of a mesh_graph of neural populations during a
random dot motion task, where `dims` is the size of the mesh graph, `inputs`
contains all of the vertices of `g` receiving external input, and `c` and `μ0`
are the coherence strength of external input. Return a tuple
(g, t, Sa, Sb, Inoisea, Inoiseb, Ineta, Inetb, xa, xb, ra, rb).
"""
function run(dims::Tuple{Vararg{<:Integer, N}}; inputs=ifelse(N==1, 1, 1:first(dims)), weight=0.015, c=0, μ0=30.0, kwargs...) where N
    return run(mesh_graph(dims; self=false, weight=weight), inputs; kwargs...)
end

"""
    run(dims::Vararg{<:Integer, N}; kwargs...) where N

Run a reduced 2-variable model of a mesh_graph of neural populations during a
random dot motion task, where `dims` is the size of the mesh graph, `inputs`
contains all of the vertices of `g` receiving external input, and `c` and `μ0`
are the coherence strength of external input. Return a tuple
(g, t, Sa, Sb, Inoisea, Inoiseb, Ineta, Inetb, xa, xb, ra, rb).
"""
function run(dims::Vararg{<:Integer, N}; kwargs...) where N
    return run(dims; kwargs...)
end

"""
    RBplot(g::AbstractGraph, Sa::Vector{<:Real}, Sb::Vector{<:Real};
           color1=colorant"red", color2=colorant"yellow",
           locs_x=nothing, locs_y=nothing, kwargs...)

Draw a plot of the graph g, where each vertex `v` is colored by the values
`Sa[v]` and `Sb[v]` on a scale from `color1` to `color2`.
"""
function RBplot(g::AbstractGraph, Sa::Vector{<:Real}, Sb::Vector{<:Real};
                color1=colorant"red", color2=colorant"yellow",
                locs_x=nothing, locs_y=nothing, kwargs...)
    nodefillc = @. weighted_color_mean((Sa-Sb+1)/2, color1, color2)
    if locs_x != nothing && locs_y != nothing
        return gplot(DiGraph(weights(g)), locs_x, locs_y; nodefillc=nodefillc, kwargs...)
    end
    return gplot(DiGraph(weights(g)); nodefillc=nodefillc, kwargs...)
end

"""
    RBgif(g::AbstractGraph, Sa::Matrix{<:Real}, Sb::Matrix{<:Real};
          layout::Function=spring_layout,fps:Integer=100, skip::Integer=1,
          width::Integer=800, height::Integer=600, kwargs...)

Draw a gif of the graph g, where each vertex `v` is colored over time by the
values `Sa[v, t]` and `Sb[v, t]` on a scale from `color1` to `color2`.
"""
function RBgif(g::AbstractGraph, Sa::Matrix{<:Real}, Sb::Matrix{<:Real};
               layout::Function=spring_layout,
               filename=nothing, fps::Integer=100, skip::Integer=1,
               width::Integer=800, height::Integer=600, kwargs...)
    # use the Animation object manually since frame() doesn't work on Compose objects
    anim = Animation()
    locs_x, locs_y = layout(g)
    for t in 1:skip:size(Sa, 2)
        fname = @sprintf("%06d.png", length(anim.frames)+1)
        draw(PNG(joinpath(anim.dir, fname), width, height),
             RBplot(g, Sa[:,t], Sb[:,t]; locs_x=locs_x, locs_y=locs_y, kwargs...))
        push!(anim.frames, fname)
    end
    if isnothing(filename)
        return gif(anim; fps=fps)
    end
    return gif(anim, filename; fps=fps)
end

function phaseplot(g::AbstractGraph, inputs; iterations=1, threaded=true,
                   Jex=-1.0:0.01:1.0, Jin=-1.0:0.01:1.0, kwargs...)
    diff = Array{Float64}(undef, length(Jex), length(Jin), nv(g))

    function difference(Jex, Jin)
        g, sol = run(g, inputs; Jex=Jex, Jin=Jin, display=false, kwargs...)
        return sol[:,1,end] .- sol[:,2,end]
    end

    if threaded
        @inbounds @Threads.threads for i in eachindex(Jex)
            for j in eachindex(Jin)
                D = [difference(Jex[i], Jin[j]) for itr in 1:iterations]
                diff[i, j, :] .= mean(D)
            end
        end
    else
        @inbounds for i in eachindex(Jex), j in eachindex(Jin)
            diff[i, j, :] .= mean(difference(Jex[i], Jin[j]) for itr in 1:iterations)
        end
    end
    return diff
end

end


text:
module test_julia_generated

using Test

generated_arg_types = []
non_generated_args = []

function f(x)
    if @generated
        push!(generated_arg_types, x)
        :x
    else
        push!(non_generated_args, x)
        0
    end
end

f(1)
@test generated_arg_types == [Int]
f("")
@test generated_arg_types == [Int, String]
f(2)
@test generated_arg_types == [Int, String]

# https://discourse.julialang.org/t/testing-function-with-if-generated/17823/8
function _invoke(f, x...; generated)
    ci = code_lowered(f, Core.Typeof(x); generated)[1]
    g = @eval @generated function g($([gensym() for _ in x]...))
        return $(QuoteNode(ci))
    end
    Base.invokelatest(g, x...)
end

_invoke(f, 1; generated=false)
@test non_generated_args == [1]

end # module test_julia_generated


text:
@params mutable struct Compliance{T} <: AbstractFunction{T}
	problem::StiffnessTopOptProblem
    solver::AbstractDisplacementSolver
    comp::T
    cell_comp::AbstractVector
    grad::AbstractVector
    tracing::Bool
    topopt_trace::TopOptTrace{T}
    reuse::Bool
    fevals::Integer
    logarithm::Bool
    maxfevals::Int
end
Utilities.getpenalty(c::Compliance) = c |> getsolver |> getpenalty
Utilities.setpenalty!(c::Compliance, p) = setpenalty!(getsolver(c), p)
Nonconvex.getdim(::Compliance) = 1

function Compliance(problem, solver::AbstractDisplacementSolver, args...; kwargs...)
    Compliance(whichdevice(solver), problem, solver, args...; kwargs...)
end
function Compliance(
    ::CPU,
    problem::StiffnessTopOptProblem{dim, T},
    solver::AbstractDisplacementSolver,
    ::Type{TI} = Int;
    tracing = false,
    logarithm = false,
    maxfevals = 10^8,
) where {dim, T, TI}
    comp = T(0)
    cell_comp = zeros(T, getncells(problem.ch.dh.grid))
    grad = fill(T(NaN), length(cell_comp) - sum(problem.black) - sum(problem.white))
    topopt_trace = TopOptTrace{T,TI}()
    reuse = false
    fevals = TI(0)
    return Compliance(problem, solver, comp, cell_comp, grad, tracing, topopt_trace, reuse, fevals, logarithm, maxfevals)
end

function (o::Compliance{T})(x, grad = o.grad) where {T}
    @unpack cell_comp, solver, tracing, topopt_trace = o
    @unpack elementinfo, u, xmin = solver
    @unpack metadata, Kes, black, white, varind = elementinfo
    @unpack cell_dofs = metadata

    @timeit to "Eval obj and grad" begin
        #if o.solver.vars ≈ x && getpenalty(o).p ≈ getprevpenalty(o).p
        #    grad .= o.grad
        #    return o.comp
        #end
        penalty = getpenalty(o)
        if o.reuse
            if !tracing
                o.reuse = false
            end
        else
            o.fevals += 1
            setpenalty!(solver, penalty.p)
            solver.vars .= x
            solver()
        end
        obj = compute_compliance(cell_comp, grad, cell_dofs, Kes, u,
                                    black, white, varind, solver.vars, penalty, xmin)

        if o.logarithm
            o.comp = log(obj)
            grad ./= obj
        else
            o.comp = obj
            #scale!(grad, 1/length(cell_comp))
            #o.comp = obj
        end
        if o.grad !== grad
            copyto!(o.grad, grad)
        end

        if o.tracing
            if o.reuse
                o.reuse = false
            else
                push!(topopt_trace.c_hist, obj)
                push!(topopt_trace.x_hist, copy(x))
                if length(topopt_trace.x_hist) == 1
                    push!(topopt_trace.add_hist, 0)
                    push!(topopt_trace.rem_hist, 0)
                else
                    push!(topopt_trace.add_hist, sum(topopt_trace.x_hist[end] .> topopt_trace.x_hist[end-1]))
                    push!(topopt_trace.rem_hist, sum(topopt_trace.x_hist[end] .< topopt_trace.x_hist[end-1]))
                end
            end
        end
    end
    return o.comp::T
end

function ChainRulesCore.rrule(comp::Compliance, x)
    out = comp(x, comp.grad)
    out_grad = copy(comp.grad)
    return out, Δ -> (nothing, out_grad * Δ)
end

function compute_compliance(cell_comp::Vector{T}, grad, cell_dofs, Kes, u,
                            black, white, varind, x, penalty, xmin) where {T}
    obj = zero(T)
    grad .= 0
    @inbounds for i in 1:size(cell_dofs, 2)
        cell_comp[i] = zero(T)
        Ke = rawmatrix(Kes[i])
        for w in 1:size(Ke,2)
            for v in 1:size(Ke, 1)
                cell_comp[i] += u[cell_dofs[v,i]]*Ke[v,w]*u[cell_dofs[w,i]]
            end
        end

        if black[i]
            obj += cell_comp[i]
        elseif white[i]
            if PENALTY_BEFORE_INTERPOLATION
                obj += xmin * cell_comp[i]
            else
                p = penalty(xmin) * cell_comp[i]
            end
        else
            d = ForwardDiff.Dual{T}(x[varind[i]], one(T))
            if PENALTY_BEFORE_INTERPOLATION
                p = density(penalty(d), xmin)
            else
                p = penalty(density(d, xmin))
            end
            grad[varind[i]] = -p.partials[1] * cell_comp[i]
            obj += p.value * cell_comp[i]
        end
    end

    return obj
end

function compute_inner(inner, u1, u2, solver)
    @unpack elementinfo, u, xmin = solver
    @unpack metadata, Kes, black, white, varind = elementinfo
    @unpack cell_dofs = metadata
    penalty = getpenalty(solver)
    return compute_inner(inner, u1, u2, cell_dofs, Kes,
        black, white, varind, solver.vars, penalty, xmin)
end
function compute_inner(inner::AbstractVector{T}, u1, u2, cell_dofs, Kes,
                            black, white, varind, x, penalty, xmin) where {T}
    obj = zero(T)
    @inbounds for i in 1:size(cell_dofs, 2)
        inner[i] = zero(T)
        cell_comp = zero(T)
        Ke = rawmatrix(Kes[i])
        if !black[i] && !white[i]
            for w in 1:size(Ke,2)
                for v in 1:size(Ke, 1)
                    cell_comp += u1[cell_dofs[v,i]]*Ke[v,w]*u2[cell_dofs[w,i]]
                end
            end
            d = ForwardDiff.Dual{T}(x[varind[i]], one(T))
            if PENALTY_BEFORE_INTERPOLATION
                p = density(penalty(d), xmin)
            else
                p = penalty(density(d, xmin))
            end
            inner[varind[i]] = -p.partials[1] * cell_comp
        end
    end

    return inner
end


text:
using PkgVersionHelper
using Test

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


text:
################################
######     Fixations      ######
################################

# E[Dn]  = LT(E[Dn+] + E[Dn-] + E[Dns])

# Neutrshape fixation corrected by background selection
"""

 - fix_neut()

Expected neutrshape fixations rate reduce by B vshapeue.

```math
\\mathbb{E}[D_{s}] = (1 - p_{-} - p_{+}) B \\frac{1}{2N}
```
# Returns
 - `Float64`: expected rate of neutrshape fixations.

"""
function fix_neut(param::parameters)
	@unpack B,NN = param;
	# Synonymous probabilty * (fixation probabilty corrected by BGS vshapeue)
	out::Float64 = 0.25*(1.0/(B*NN))
	return out
end

# Negative fixations corrected by background selection
"""

	fix_neg_b(p)

Expected fixation rate from negative DFE.

```math
\\mathbb{E}[D_{n-}] =  p_{-}\\left(2^-\\shapepha\\scaleta^\\shapepha\\left(-\\zeta[\\shapepha,\\frac{2+\\scaleta}{2}] + \\zeta[\\shapepha,1/2(2-\\frac{1}{N+\\scaleta})]\\right)\\right)
```

# Arguments
 - `p::Float64`: positive selected alleles probabilty.

# Returns
 - `Float64`: expected rate of fixations from negative DFE.

"""
function fix_neg_b(param::parameters,p::Float64)
	@unpack shape, B, scale, N = param;
	# Non-synonymous proportion * negative alleles probability * fixation probability from gamma distribution
	out::Float64 = 0.75*(1-p)*(2^(-shape))*(B^(-shape))*(scale^shape)*(-zeta(shape,1.0+scale/(2.0*B))+zeta(shape,0.5*(2-1.0/(N*B)+scale/B)))
	return out
end

# Positive fixations
"""
	p_fix()

Expected positive fixation rate.

```math
\\mathbb{E}[D_{n+}] =  p_{+} \\cdot B \\cdot (1 - e^{(-2s)})
```

# Arguments
 - `p::Float64`: positive selected alleles probabilty.

# Returns
 - `Float64`: expected rate of positive fixation.

"""
function p_fix(param::parameters,s::Int64)

	# Fixation probability
	# @unpack NN   = param;
	s::Float64     = s/(param.NN);
	p_fix::Float64 = (1.0-ℯ^(-2.0*s))/(1.0-ℯ^(-2.0*s));

	# Correcting p_fix for large s following Uricchio et shape. 2014
	if s >= 0.1
		p_fix = ℯ^(-(1.0+s));
		lim   = 0;
		while(lim < 200)
			p_fix = ℯ^((1.0+s)*(p_fix-1.0));
			lim   = lim + 1;
		end
		p_fix = 1 - p_fix;
	end

	return p_fix
end

# Positive fixations after apply Φ. reduction of positive fixations due deleterious linkage given a vshapeue B of background selection
"""

	fixPosSim(gamma,p)

Expected positive fixations rate reduced due to the impact of background selection and linkage. The probabilty of fixation of positively selected alleles is reduced by a factor Φ across shapel deleterious linked sites [`Anshapeyticshape.phiReduction`](@ref).

```math
\\mathbb{E}[D_{n+}] =  \\Phi \\cdot \\mathbb{E}[D_{n+}]
```

# Arguments
 - `p::Float64`: positive selected alleles probabilty.

# Returns
 - `Float64`: expected rate of positive fixations under background selection

"""
function fix_pos_sim(param::parameters,s::Int64,p::Float64)

	# Non-synonymous * positive alleles probability * B reduction * fixation probility
	red_plus     = Φ(param,s)
	out::Float64 = 0.75*p*red_plus*p_fix(param,s)
	return out
end