text:
using CairoMakie, AbstractPlotting, LinearAlgebra
# %%

const SIZE = 1 << 8
const MID = SIZE ÷ 2
const RAND_BORDER = SIZE ÷ 6
const RAND_RANGE = range(RAND_BORDER, SIZE-RAND_BORDER, step=1)
const RAND_DEN = 0.5

# %%
mutable struct Lenia
    R::Int64
    peaks::Array{Float64, 1}
    μ::Float64
    σ::Float64
    dt::Float64
    kernel_type::Int64
    delta_type::Int64
    kernel_core::Function
    kernel_shell::Function
    delta_func::Function
end

# %%
"""
I collapsed the `set_params` function into the struct constructor
"""
function Lenia(R::Int64, peaks::Array{Float64, 1}, μ::Float64, σ::Float64, dt::Float64, kernel_type::Int64=0, delta_type::Int64=0)

    function kernel_core(s::Lenia, r::Float64)
        rm = min(r, 1)
        if s.kernel_type == 0
            return (4 * rm * (1-rm))^4
        else
            return exp( 4 - 1 / (rm * (1-rm)) )
        end
    end

    function kernel_shell(s::Lenia, r::Float64)
        k = length(s.peaks)
        kr = k * r
        peak = s.peaks[min(floor(Int64, kr)+1, k)] ## Julia's indexing starts from 1
        println(peak)
        return (r<1) * s.kernel_core(s, kr % 1) * peak
    end

    function delta_func(s::Lenia, n::Float64)
        if s.delta_type == 0
            return max(0, 1 - (n - s.μ)^2 / (s.σ^2 * 9) )^4 * 2 - 1
        else
            return exp( - (n - s.μ)^2 / (s.σ^2 * 2) ) * 2 - 1
        end
    end

    Lenia(R, peaks, μ, σ, dt, kernel_type, delta_type, kernel_core, kernel_shell, delta_func)
end
# %%

lenia = Lenia(52, [1/2, 2/3, 1], 0.3, 0.03, 0.1)
# %%
k_core(x) = lenia.kernel_core(lenia, x)
k_shell(x) = lenia.kernel_shell(lenia, x)
δ(n) = lenia.delta_func(lenia, n)
# %%

function plot_funcs()
    x = LinRange(0, 1.2, 120)
    fig = Figure(resolution = (1024, 512))
    ax1 = Axis(fig[1, 1], title = "Kernel's core and shell", xlims=(0, 1.2))
    lines!(ax1, x, k_core.(x), color = :cyan)
    lines!(ax1, x, k_shell.(x), color=:blue)
    ax2 = Axis(fig[1, 2], title="δ function")
    lines!(ax2, x, δ.(x), color =:green)
    fig
end

# %%
fig = plot_funcs()
# %%
#-----## Here ends the Julia code ##-----#
# %%

##

# %%

def show_world(s, A, vmin=0, vmax=1, is_display=True):
    dpi = 80
    fig = plt.figure(figsize=(np.shape(A)[1]/dpi, np.shape(A)[0]/dpi), dpi=dpi)
    ax = fig.add_axes([0, 0, 1, 1])
    ax.axis('off')
    img = ax.imshow(A, cmap='jet', interpolation='none', aspect=1, vmin=vmin, vmax=vmax)
    #fig.figimage(A)
    #plt.subplot(1, 2, 1)
    #plt.imshow(A, cmap='jet', interpolation='none') #nipy_spectral gist_ncar
    if is_display: plt.show()
    return fig, img
Lenia.show_world = show_world

def calc_kernel(s):
    I = np.array([np.arange(SIZE),]*SIZE)
    X = (I-MID) / s.R
    Y = X.T
    D = np.sqrt(X**2 + Y**2)

    s.kernel = s.kernel_shell(D)
    s.kernel_sum = np.sum(s.kernel)
    kernel_norm = s.kernel / s.kernel_sum
    s.kernel_FFT = np.fft.fft2(kernel_norm)
Lenia.calc_kernel = calc_kernel

lenia.calc_kernel()
#lenia.show_world(lenia.D, 0, None)
fig = lenia.show_world(lenia.kernel)
#lenia.show_world(np.fft.fftshift(abs(lenia.kernel_FFT)))

def load_cells(s, id):
    if id==0: s.name='Orbium bicaudatus'; s.set_params(R=13,peaks=np.array([1]),mu=0.15,sigma=0.014,dt=0.1);s.cells=np.array([[0,0,0,0,0,0,0.1,0.14,0.1,0,0,0.03,0.03,0,0,0.3,0,0,0,0],[0,0,0,0,0,0.08,0.24,0.3,0.3,0.18,0.14,0.15,0.16,0.15,0.09,0.2,0,0,0,0],[0,0,0,0,0,0.15,0.34,0.44,0.46,0.38,0.18,0.14,0.11,0.13,0.19,0.18,0.45,0,0,0],[0,0,0,0,0.06,0.13,0.39,0.5,0.5,0.37,0.06,0,0,0,0.02,0.16,0.68,0,0,0],[0,0,0,0.11,0.17,0.17,0.33,0.4,0.38,0.28,0.14,0,0,0,0,0,0.18,0.42,0,0],[0,0,0.09,0.18,0.13,0.06,0.08,0.26,0.32,0.32,0.27,0,0,0,0,0,0,0.82,0,0],[0.27,0,0.16,0.12,0,0,0,0.25,0.38,0.44,0.45,0.34,0,0,0,0,0,0.22,0.17,0],[0,0.07,0.2,0.02,0,0,0,0.31,0.48,0.57,0.6,0.57,0,0,0,0,0,0,0.49,0],[0,0.59,0.19,0,0,0,0,0.2,0.57,0.69,0.76,0.76,0.49,0,0,0,0,0,0.36,0],[0,0.58,0.19,0,0,0,0,0,0.67,0.83,0.9,0.92,0.87,0.12,0,0,0,0,0.22,0.07],[0,0,0.46,0,0,0,0,0,0.7,0.93,1,1,1,0.61,0,0,0,0,0.18,0.11],[0,0,0.82,0,0,0,0,0,0.47,1,1,0.98,1,0.96,0.27,0,0,0,0.19,0.1],[0,0,0.46,0,0,0,0,0,0.25,1,1,0.84,0.92,0.97,0.54,0.14,0.04,0.1,0.21,0.05],[0,0,0,0.4,0,0,0,0,0.09,0.8,1,0.82,0.8,0.85,0.63,0.31,0.18,0.19,0.2,0.01],[0,0,0,0.36,0.1,0,0,0,0.05,0.54,0.86,0.79,0.74,0.72,0.6,0.39,0.28,0.24,0.13,0],[0,0,0,0.01,0.3,0.07,0,0,0.08,0.36,0.64,0.7,0.64,0.6,0.51,0.39,0.29,0.19,0.04,0],[0,0,0,0,0.1,0.24,0.14,0.1,0.15,0.29,0.45,0.53,0.52,0.46,0.4,0.31,0.21,0.08,0,0],[0,0,0,0,0,0.08,0.21,0.21,0.22,0.29,0.36,0.39,0.37,0.33,0.26,0.18,0.09,0,0,0],[0,0,0,0,0,0,0.03,0.13,0.19,0.22,0.24,0.24,0.23,0.18,0.13,0.05,0,0,0,0],[0,0,0,0,0,0,0,0,0.02,0.06,0.08,0.09,0.07,0.05,0.01,0,0,0,0,0]])
    elif id==1: s.name='Rotorbium'; s.set_params(R=13,peaks=np.array([1]),mu=0.156,sigma=0.0224,dt=0.1);s.cells=np.array([[0,0,0,0,0,0,0,0,0.003978,0.016492,0.004714,0,0,0,0,0,0,0,0,0],[0,0,0,0,0,0.045386,0.351517,0.417829,0.367137,0.37766,0.426948,0.431058,0.282864,0.081247,0,0,0,0,0,0],[0,0,0,0,0.325473,0.450995,0.121737,0,0,0,0.003113,0.224278,0.47101,0.456459,0.247231,0.071609,0.013126,0,0,0],[0,0,0,0.386337,0.454077,0,0,0,0,0,0,0,0.27848,0.524466,0.464281,0.242651,0.096721,0.038476,0,0],[0,0,0.258817,0.583802,0.150994,0,0,0,0,0,0,0,0.226639,0.548329,0.550422,0.334764,0.153108,0.087049,0.042872,0],[0,0.008021,0.502406,0.524042,0.059531,0,0,0,0,0,0,0.033946,0.378866,0.615467,0.577527,0.357306,0.152872,0.090425,0.058275,0.023345],[0,0.179756,0.596317,0.533619,0.162612,0,0,0,0,0.015021,0.107673,0.325125,0.594765,0.682434,0.594688,0.381172,0.152078,0.073544,0.054424,0.030592],[0,0.266078,0.614339,0.605474,0.379255,0.195176,0.16516,0.179148,0.204498,0.299535,0.760743,1,1,1,1,0.490799,0.237826,0.069989,0.043549,0.022165],[0,0.333031,0.64057,0.686886,0.60698,0.509866,0.450525,0.389552,0.434978,0.859115,0.94097,1,1,1,1,1,0.747866,0.118317,0.037712,0.006271],[0,0.417887,0.6856,0.805342,0.824229,0.771553,0.69251,0.614328,0.651704,0.843665,0.910114,1,1,0.81765,0.703404,0.858469,1,0.613961,0.035691,0],[0.04674,0.526827,0.787644,0.895984,0.734214,0.661746,0.670024,0.646184,0.69904,0.723163,0.682438,0.618645,0.589858,0.374017,0.30658,0.404027,0.746403,0.852551,0.031459,0],[0.130727,0.658494,0.899652,0.508352,0.065875,0.009245,0.232702,0.419661,0.461988,0.470213,0.390198,0.007773,0,0.010182,0.080666,0.17231,0.44588,0.819878,0.034815,0],[0.198532,0.810417,0.63725,0.031385,0,0,0,0,0.315842,0.319248,0.321024,0,0,0,0,0.021482,0.27315,0.747039,0,0],[0.217619,0.968727,0.104843,0,0,0,0,0,0.152033,0.158413,0.114036,0,0,0,0,0,0.224751,0.647423,0,0],[0.138866,1,0.093672,0,0,0,0,0,0.000052,0.015966,0,0,0,0,0,0,0.281471,0.455713,0,0],[0,1,0.145606,0.005319,0,0,0,0,0,0,0,0,0,0,0,0.016878,0.381439,0.173336,0,0],[0,0.97421,0.262735,0.096478,0,0,0,0,0,0,0,0,0,0,0.013827,0.217967,0.287352,0,0,0],[0,0.593133,0.2981,0.251901,0.167326,0.088798,0.041468,0.013086,0.002207,0.009404,0.032743,0.061718,0.102995,0.1595,0.24721,0.233961,0.002389,0,0,0],[0,0,0.610166,0.15545,0.200204,0.228209,0.241863,0.243451,0.270572,0.446258,0.376504,0.174319,0.154149,0.12061,0.074709,0,0,0,0,0],[0,0,0.354313,0.32245,0,0,0,0.151173,0.479517,0.650744,0.392183,0,0,0,0,0,0,0,0,0],[0,0,0,0.329339,0.328926,0.176186,0.198788,0.335721,0.534118,0.549606,0.361315,0,0,0,0,0,0,0,0,0],[0,0,0,0,0.090407,0.217992,0.190592,0.174636,0.222482,0.375871,0.265924,0,0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0,0,0.050256,0.235176,0,0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0,0,0,0,0.180145,0.132616,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0,0,0,0,0,0,0.092581,0.188519,0.118256,0,0,0,0]])
Lenia.load_cells = load_cells

def clear_world(s):
    s.world = np.zeros((SIZE, SIZE))
Lenia.clear_world = clear_world

def add_cells(s):
    w = np.shape(s.cells)[1]
    h = np.shape(s.cells)[0]
    i = MID - int(w / 2)
    j = MID - int(h / 2)
    s.world[j:j+h, i:i+w] = s.cells
Lenia.add_cells = add_cells

def multiply_cells(s, n=2):
    w = np.shape(s.cells)[1] * n
    h = np.shape(s.cells)[0] * n
    cells2 = np.zeros((h, w))
    for i in range(n):
        for j in range(n):
            cells2[i:h:n, j:w:n] = s.cells
    s.cells = cells2
    s.R = s.R * n
Lenia.multiply_cells = multiply_cells

lenia.clear_world()
lenia.load_cells(0)
lenia.multiply_cells(4)
lenia.calc_kernel()
lenia.add_cells()
fig = lenia.show_world(lenia.world)
#fig = lenia.show_world(lenia.kernel)

def calc_once(s, is_update=True):
    world_FFT = np.fft.fft2(s.world)
    s.potential = np.roll(np.real(np.fft.ifft2(s.kernel_FFT * world_FFT)), MID, (0, 1))
    s.delta = s.delta_func(s.potential)
    if is_update: s.world = np.maximum(0, np.minimum(1, s.world + s.delta * s.dt))
Lenia.calc_once = calc_once

lenia.calc_once()
fig, img = lenia.show_world(lenia.world, is_display=False)
#img = lenia.show_world(lenia.potential, 0, 0.3)
#img = lenia.show_world(lenia.delta, -1, 1)

def animate(i):
    lenia.calc_once()
    img.set_array(lenia.world)
    return img,

plt.rcParams["animation.html"] = "jshtml"

anim = animation.FuncAnimation(fig, animate, frames=100, interval=50, blit=True)
HTML(anim.to_jshtml())

def run_lenia(animal_id):
    def run_animate(i):
        lenia.calc_once()
        img.set_array(lenia.world)
        return img,

    lenia = Lenia()
    lenia.clear_world()
    lenia.load_cells(animal_id)
    lenia.multiply_cells(4)
    lenia.calc_kernel()
    lenia.add_cells()
    fig, img = lenia.show_world(lenia.world, is_display=False)
    return animation.FuncAnimation(fig, run_animate, frames=200, interval=50, blit=True)

i = interact(run_lenia, animal_id={
    'Orbium bicaudatus':0,
    'Rotorbium':1})

max_stars_repo_name: izzortsi/Lenia

max_issues_repo_name: izzortsi/Lenia

max_forks_repo_name: izzortsi/Lenia


text:
module ExponentialTest

using Test
using ForneyLab
using ForneyLab: outboundType, isApplicable
using ForneyLab: SPExponentialOutNG, SPExponentialOutNP, SPExponentialIn1LN, SPExponentialIn1PN

@testset "Exponential node construction through exp() syntax" begin
    g = FactorGraph()

    @RV x ~ GaussianMeanVariance(constant(0.0), constant(1.0))
    @RV z = exp(x)

    @test isa(z, Variable)
    @test isa(g.nodes[:exponential_1], Exponential)
end

#-------------
# Update rules
#-------------

@testset "SPExponentialOutNG" begin
    @test SPExponentialOutNG <: SumProductRule{Exponential}
    @test outboundType(SPExponentialOutNG) == Message{LogNormal}
    @test isApplicable(SPExponentialOutNG, [Nothing, Message{Gaussian}])

    @test ruleSPExponentialOutNG(nothing, Message(Univariate, GaussianMeanVariance, m=1.0, v=2.0)) == Message(Univariate, LogNormal, m=1.0, s=2.0)
end

@testset "SPExponentialOutNP" begin
    @test SPExponentialOutNP <: SumProductRule{Exponential}
    @test outboundType(SPExponentialOutNP) == Message{PointMass}
    @test isApplicable(SPExponentialOutNP, [Nothing, Message{PointMass}])

    @test ruleSPExponentialOutNP(nothing, Message(Univariate, PointMass, m=2.0)) == Message(Univariate, PointMass, m=exp(2.0))
end

@testset "SPExponentialIn1LN" begin
    @test SPExponentialIn1LN <: SumProductRule{Exponential}
    @test outboundType(SPExponentialIn1LN) == Message{GaussianMeanVariance}
    @test isApplicable(SPExponentialIn1LN, [Message{LogNormal}, Nothing])

    @test ruleSPExponentialIn1LN(Message(Univariate, LogNormal, m=1.0, s=2.0), nothing) == Message(Univariate, GaussianMeanVariance, m=1.0, v=2.0)
end

@testset "SPExponentialIn1PN" begin
    @test SPExponentialIn1PN <: SumProductRule{Exponential}
    @test outboundType(SPExponentialIn1PN) == Message{PointMass}
    @test isApplicable(SPExponentialIn1PN, [Message{PointMass}, Nothing])

    @test ruleSPExponentialIn1PN(Message(Univariate, PointMass, m=2.0), nothing) == Message(Univariate, PointMass, m=log(2.0))
end

end # module

ext: jl


text:
module ArrayBenchmarks

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

using .RandUtils
using BenchmarkTools
using Compat

import Compat: UTF8String, view

const SUITE = BenchmarkGroup()

#############################################################################
# basic array-math reduction-like functions

afloat = samerand(10^3)
aint = samerand(Int, 10^3)
acomplex = samerand(Complex{Float64}, 10^3)
g = addgroup!(SUITE, "reductions", ["sum", "array", "reduce"])
norm1(x) = norm(x, 1)
norminf(x) = norm(x, Inf)
perf_reduce(x) = reduce((x,y) -> x + 2y, real(zero(eltype(x))), x)
perf_mapreduce(x) = mapreduce(x -> real(x)+imag(x), (x,y) -> x + 2y, real(zero(eltype(x))), x)
for a in (afloat, aint)
    for fun in (sum, norm, norm1, norminf, mean, var, perf_reduce, perf_mapreduce)
        g[string(fun), string(eltype(a))] = @benchmarkable $fun($a)
    end
    g["sumabs2", string(eltype(a))] = @benchmarkable sum(abs2, $a)
    g["sumabs", string(eltype(a))] = @benchmarkable sum(abs, $a)
    g["maxabs", string(eltype(a))] = @benchmarkable maximum(abs, $a)
end

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

############
# indexing #
############

# #10525 #
#--------#

include("sumindex.jl")

s = 500
A3d = samerand(11,11,11)
S3d = view(A3d, 1:10, 1:10, 1:10)
arrays = (makearrays(Int32, s, s)..., makearrays(Float32, s, s)..., trues(s, s), A3d, S3d)
ranges = (1:10^5, 10^5:-1:1, 1.0:1e5, linspace(1,2,10^4))
arrays_iter = map(x -> (x, string(typeof(x))), arrays)
ranges_iter = map(x -> (x, repr(x)), ranges)
g = addgroup!(SUITE, "index", ["sum", "simd"])

for (A, str) in (arrays_iter..., ranges_iter...)
    g["sumelt", str]             = @benchmarkable perf_sumelt($A)
    g["sumelt_boundscheck", str] = @benchmarkable perf_sumelt_boundscheck($A)
    g["sumeach", str]            = @benchmarkable perf_sumeach($A)
    g["sumlinear", str]          = @benchmarkable perf_sumlinear($A)
    g["sumcartesian", str]       = @benchmarkable perf_sumcartesian($A)
    g["sumeach_view", str]       = @benchmarkable perf_sumeach_view($A)
    g["sumlinear_view", str]     = @benchmarkable perf_sumlinear_view($A)
    g["sumcartesian_view", str]  = @benchmarkable perf_sumcartesian_view($A)
    if ndims(A) == 2
        g["mapr_access", str]    = @benchmarkable perf_mapr_access($A, B, zz, n) setup = begin B, zz, n = setup_mapr_access($A) end #20517
    end
    if ndims(A) <= 2
        g["sumcolon", str]       = @benchmarkable perf_sumcolon($A)
        g["sumrange", str]       = @benchmarkable perf_sumrange($A)
        g["sumlogical", str]     = @benchmarkable perf_sumlogical($A)
        g["sumvector", str]      = @benchmarkable perf_sumvector($A)
        g["sumcolon_view", str]  = @benchmarkable perf_sumcolon_view($A)
        g["sumrange_view", str]  = @benchmarkable perf_sumrange_view($A)
        g["sumlogical_view", str]= @benchmarkable perf_sumlogical_view($A)
        g["sumvector_view", str] = @benchmarkable perf_sumvector_view($A)
    end
end
g["sub2ind"] = @benchmarkable perf_sub2ind((1000,1000,1000), 1:1000, 1:1000, 1:1000)
g["ind2sub"] = @benchmarkable perf_ind2sub((100,100,10), 1:10^5)
g["sum", "3darray"] = @benchmarkable sum($A3d)
g["sum", "3dsubarray"] = @benchmarkable sum($S3d)

for b in values(g)
    b.params.time_tolerance = 0.50
end

# #18774 #
#--------#
include("generate_kernel.jl")

nmax = 6  # maximum dimensionality is nmax + 1
# get path to current directory and append file name to it
fname = joinpath(dirname(@__FILE__), "hdindexing.jl")
make_stencil(fname, nmax)  # generate source file
include("hdindexing.jl")

npts_dir = [10000, 80, 20, 12, 9, 6]  # number of points in each direction
for i=1:nmax
  dims_vec = zeros(Int, i+1)
  fill!(dims_vec, npts_dir[i])
  dims_vec[end] = 2  # last dimension must be 2
  u_i = samerand(dims_vec...)
  u_ip1 = zeros(dims_vec...)

  str = string(i+1, "d")

  # perf_hdindexing is defined in the generated source file hdindexing.jl
  g[str] = @benchmarkable perf_hdindexing5($u_i, $u_ip1)
end

# #10301 #
#--------#

include("revloadindex.jl")

v = samerand(10^6)
n = samerand()

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

g["rev_load_slow!"]    = @benchmarkable perf_rev_load_slow!(fill!($v, $n))
g["rev_load_fast!"]    = @benchmarkable perf_rev_load_fast!(fill!($v, $n))
g["rev_loadmul_slow!"] = @benchmarkable perf_rev_loadmul_slow!(fill!($v, $n), $v)
g["rev_loadmul_fast!"] = @benchmarkable perf_rev_loadmul_fast!(fill!($v, $n), $v)

# #9622 #
#-------#

perf_setindex!(A, val, inds) = setindex!(A, val, inds...)

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

for s in (1, 2, 3, 4, 5)
    A = samerand(Float64, ntuple(one, s)...)
    y = one(eltype(A))
    i = length(A)
    g["setindex!", ndims(A)] = @benchmarkable perf_setindex!(fill!($A, $y), $y, $i)
end

###############################
# SubArray (views vs. copies) #
###############################

# LU factorization with complete pivoting. These functions deliberately allocate
# a lot of temprorary arrays by working on vectors instead of looping through
# the elements of the matrix. Both a view (SubArray) version and a copy version
# are provided.

include("subarray.jl")

n = samerand()

g = addgroup!(SUITE, "subarray", ["lucompletepiv"])

for s in (100, 250, 500, 1000)
    m = samerand(s, s)
    g["lucompletepivCopy!", s] = @benchmarkable perf_lucompletepivCopy!(fill!($m, $n))
    g["lucompletepivSub!", s]  = @benchmarkable perf_lucompletepivSub!(fill!($m, $n))
end

#################
# concatenation #
#################

include("cat.jl")

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

for s in (5, 500)
    A = samerand(s, s)
    g["hvcat", s]        = @benchmarkable perf_hvcat($A, $A)
    g["hcat", s]         = @benchmarkable perf_hcat($A, $A)
    g["vcat", s]         = @benchmarkable perf_vcat($A, $A)
    g["catnd", s]        = @benchmarkable perf_catnd($s)
    g["hvcat_setind", s] = @benchmarkable perf_hvcat_setind($A, $A)
    g["hcat_setind", s]  = @benchmarkable perf_hcat_setind($A, $A)
    g["vcat_setind", s]  = @benchmarkable perf_vcat_setind($A, $A)
    g["catnd_setind", s] = @benchmarkable perf_catnd_setind($s)
end

############################
# in-place growth (#13977) #
############################

function perf_push_multiple!(collection, items)
    for item in items
        push!(collection, item)
    end
    return collection
end

g = addgroup!(SUITE, "growth", ["push!", "append!", "prepend!"])

for s in (8, 256, 2048)
    v = samerand(s)
    g["push_single!", s]   = @benchmarkable push!(x, samerand())       setup=(x = copy($v))
    g["push_multiple!", s] = @benchmarkable perf_push_multiple!(x, $v) setup=(x = copy($v))
    g["append!", s]        = @benchmarkable append!(x, $v)             setup=(x = copy($v))
    g["prerend!", s]       = @benchmarkable prepend!(x, $v)            setup=(x = copy($v))
end

##########################
# comprehension (#13401) #
##########################

perf_compr_collect(X) = [x for x in X]
perf_compr_iter(X) = [sin(x) + x^2 - 3 for x in X]
perf_compr_index(X) = [sin(X[i]) + (X[i])^2 - 3 for i in eachindex(X)]

ls = linspace(0,1,10^7)
rg = 0.0:(10.0^(-7)):1.0
arr = collect(ls)

g = addgroup!(SUITE, "comprehension", ["iteration", "index", "linspace", "collect", "range"])

for X in (ls, rg, arr)
    T = string(typeof(X))
    g["collect", T] = @benchmarkable collect($X)
    g["comprehension_collect", T]   = @benchmarkable perf_compr_collect($X)
    g["comprehension_iteration", T] = @benchmarkable perf_compr_iter($X)
    g["comprehension_indexing", T]  = @benchmarkable perf_compr_index($X) time_tolerance=0.30
end

###############################
# BoolArray/BitArray (#13946) #
###############################

function perf_bool_load!(result, a, b)
    for i in eachindex(result)
        result[i] = a[i] != b
    end
    return result
end

function perf_true_load!(result)
    for i in eachindex(result)
        result[i] = true
    end
    return result
end

n, range = 10^6, -3:3
a, b = samerand(range, n), samerand(range)
boolarr, bitarr = Vector{Bool}(n), BitArray(n)

g = addgroup!(SUITE, "bool", ["index", "bitarray", "fill!"])

g["bitarray_bool_load!"]  = @benchmarkable perf_bool_load!($bitarr, $a, $b)
g["boolarray_bool_load!"] = @benchmarkable perf_bool_load!($boolarr, $a, $b)
g["bitarray_true_load!"]  = @benchmarkable perf_true_load!($bitarr)
g["boolarray_true_load!"] = @benchmarkable perf_true_load!($boolarr)
g["bitarray_true_fill!"]  = @benchmarkable fill!($bitarr, true)
g["boolarray_true_fill!"] = @benchmarkable fill!($boolarr, true)

####################################
# Float to Int conversion (#18954) #
####################################

function perf_convert!(a, x)
    for i=1:length(x)
        a[i] = x[i]
    end
    return a
end

x_int = rand(1:1000000,100,100)
x_float = 1.0 * x_int
x_complex = x_float .+ 0.0im

g = addgroup!(SUITE, "convert", ["Int"])
g["Int", "Float64"] = @benchmarkable  perf_convert!($x_int, $x_float)
g["Float64", "Int"] = @benchmarkable  perf_convert!($x_float, $x_int)
g["Complex{Float64}", "Int"] = @benchmarkable  perf_convert!($x_complex, $x_int)
g["Int", "Complex{Float64}"] = @benchmarkable  perf_convert!($x_int, $x_complex)

end # module


text:
module AdventOfCode

using Combinatorics
mutable struct Moon pos; vel end

gravity!(m1::Moon, m2::Moon) = (dv = (m1.pos .< m2.pos) - (m1.pos .> m2.pos); m1.vel += dv; m2.vel -= dv)
gravity!(moons) = map(ms -> gravity!(ms...), combinations(moons, 2))
velocity!(m::Moon) = m.pos += m.vel
energy(m::Moon) = sum(abs.(m.pos)) * sum(abs.(m.vel))

# Part 1: compute total energy after 1000 simulation steps
initial_moons() = [
    Moon([-7 -8 9], [0 0 0]),
    Moon([-12 -3 -4], [0 0 0]),
    Moon([6 -17 -9], [0 0 0]),
    Moon([4 -10 -6], [0 0 0])
]
jupiter_moons = initial_moons()
for step in 1:1000
    gravity!(jupiter_moons)
    velocity!.(jupiter_moons)
end
println("Total energy: $(sum(energy.(jupiter_moons)))")

# Part 2: find how long it takes before first repeated step
# Approach: x, y and z dimension equations are independent, maybe each one has a shorter period?
dim_history = [Dict() Dict() Dict()]
cycle_ends = Array{Any}([nothing nothing nothing])
jupiter_moons = initial_moons()
for step in 0:10000000
    gravity!(jupiter_moons)
    velocity!.(jupiter_moons)
    for dim in 1:3
        dim_state = map(m -> (m.pos[dim], m.vel[dim]), jupiter_moons)
        if isnothing(cycle_ends[dim])
            if dim_state ∈ keys(dim_history[dim])
                cycle_ends[dim] = (dim_history[dim][dim_state], step)
            else
                dim_history[dim][dim_state] = step
            end
        end
    end
    if !any(isnothing.(cycle_ends)) break end
end

if all(0 .== first.(cycle_ends))
    period = lcm(last.(cycle_ends)...)
    println("Loop start/end: $cycle_ends, repeat period: $period")
else
    println("Ugh, the maths is harder, I have to think more to implement this")
end

end

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
@doc """
    postproc(model::Model, param::Parameter, obs::MCObservableSet)

Post process of observables. For example, Specific heat will be calculated from energy, energy^2, and temperature.
"""
function postproc end


text:
dv = dataones(3)
push!(dv, 3.0)
push!(dv, NA)

@assert isequal(dv, DataVector[1.0, 1.0, 1.0, 3.0, NA])

a, b = pop!(dv), pop!(dv)
@assert isna(a)
@assert b == 3.0

enqueue!(dv, 3.0)
enqueue!(dv, NA)

@assert isequal(dv, DataVector[NA, 3.0, 1.0, 1.0, 1.0])

a, b = shift(dv), shift(dv)
@assert isna(a)
@assert b == 3.0


text:
export
FactorFunction,
ProbabilityDistribution,
Univariate,
Multivariate,
MatrixVariate,
PointMass,
mean,
mode,
var,
cov,
differentialEntropy,
conditionalDifferentialEntropy,
averageEnergy,
==,
vague,
sample,
dims,
logPdf

abstract type VariateType end
abstract type Univariate <: VariateType end
abstract type Multivariate <: VariateType end
abstract type MatrixVariate <: VariateType end

"""Types through which a probability distribution may be defined"""
const FactorFunction = Union{FactorNode, Function}

"""Encodes a probability distribution as a `FactorFunction` of type `family` with fixed interfaces"""
struct ProbabilityDistribution{var_type<:VariateType, family<:FactorFunction}
    params::Dict
end

"""Sample multiple realizations from a probability distribution"""
sample(dist::ProbabilityDistribution, n_samples::Int64) = [sample(dist) for i in 1:n_samples] # TODO: individual samples can be optimized

"""Extract VariateType from dist"""
variateType(dist::ProbabilityDistribution{V, F}) where {V<:VariateType, F<:FactorFunction} = V

show(io::IO, dist::ProbabilityDistribution) = println(io, format(dist))

matches(Ta::Type{Pa}, Tb::Type{Pb}) where {Pa<:ProbabilityDistribution, Pb<:ProbabilityDistribution} = (Pa<:Pb)
matches(::Type{Nothing}, ::Type{T}) where T<:ProbabilityDistribution = false

mean(dist::ProbabilityDistribution) = isProper(dist) ? unsafeMean(dist) : error("mean($(dist)) is undefined because the distribution is improper.")
mode(dist::ProbabilityDistribution) = isProper(dist) ? unsafeMode(dist) : error("mode($(dist)) is undefined because the distribution is improper.")
var(dist::ProbabilityDistribution) = isProper(dist) ? unsafeVar(dist) : error("var($(dist)) is undefined because the distribution is improper.")
cov(dist::ProbabilityDistribution) = isProper(dist) ? unsafeCov(dist) : error("cov($(dist)) is undefined because the distribution is improper.")
logPdf(dist::ProbabilityDistribution) = isProper(dist) ? logPdf(dist) : error("logPdf($(dist)) is undefined.")

"""
`PointMass` is an abstract type used to describe point mass distributions.
It never occurs in a `FactorGraph`, but it is used as a probability distribution
type.
"""
abstract type PointMass <: DeltaFactor end

slug(::Type{PointMass}) = "δ"

format(dist::ProbabilityDistribution{V, PointMass}) where V<:VariateType = "$(slug(PointMass))(m=$(format(dist.params[:m])))"

dims(dist::ProbabilityDistribution{Univariate, PointMass}) = 1
dims(dist::ProbabilityDistribution{Multivariate, PointMass}) = length(dist.params[:m])
dims(dist::ProbabilityDistribution{MatrixVariate, PointMass}) = size(dist.params[:m])

# PointMass distribution constructors
ProbabilityDistribution(::Type{Univariate}, ::Type{PointMass}; m::Number=1.0) = ProbabilityDistribution{Univariate, PointMass}(Dict(:m=>m))
ProbabilityDistribution(::Type{PointMass}; m::Number=1.0) = ProbabilityDistribution{Univariate, PointMass}(Dict(:m=>m))
ProbabilityDistribution(::Type{Multivariate}, ::Type{PointMass}; m::Vector=[1.0]) = ProbabilityDistribution{Multivariate, PointMass}(Dict(:m=>m))
ProbabilityDistribution(::Type{MatrixVariate}, ::Type{PointMass}; m::AbstractMatrix=transpose([1.0])) = ProbabilityDistribution{MatrixVariate, PointMass}(Dict(:m=>m))

unsafeMean(dist::ProbabilityDistribution{T, PointMass}) where T<:VariateType = deepcopy(dist.params[:m])

unsafeMode(dist::ProbabilityDistribution{T, PointMass}) where T<:VariateType = deepcopy(dist.params[:m])

unsafeMeanVector(dist::ProbabilityDistribution{Univariate, PointMass}) = [dist.params[:m]]
unsafeMeanVector(dist::ProbabilityDistribution{Multivariate, PointMass}) = deepcopy(dist.params[:m])

unsafeInverseMean(dist::ProbabilityDistribution{Univariate, PointMass}) = 1.0/dist.params[:m]
unsafeInverseMean(dist::ProbabilityDistribution{MatrixVariate, PointMass}) = cholinv(dist.params[:m])

unsafeLogMean(dist::ProbabilityDistribution{Univariate, PointMass}) = log(clamp(dist.params[:m], tiny, Inf))
unsafeLogMean(dist::ProbabilityDistribution{Multivariate, PointMass}) = log.(clamp.(dist.params[:m], tiny, Inf))
unsafeLogMean(dist::ProbabilityDistribution{MatrixVariate, PointMass}) = log.(clamp.(dist.params[:m], tiny, Inf))

unsafeDetLogMean(dist::ProbabilityDistribution{Univariate, PointMass}) = log(dist.params[:m])
unsafeDetLogMean(dist::ProbabilityDistribution{MatrixVariate, PointMass}) = logdet(dist.params[:m])

unsafeMirroredLogMean(dist::ProbabilityDistribution{Univariate, PointMass}) = log(1.0 - dist.params[:m])

unsafeVar(dist::ProbabilityDistribution{Univariate, PointMass}) = 0.0
unsafeVar(dist::ProbabilityDistribution{Multivariate, PointMass}) = zeros(dims(dist))

unsafeCov(dist::ProbabilityDistribution{Univariate, PointMass}) = 0.0
unsafeCov(dist::ProbabilityDistribution{Multivariate, PointMass}) = zeros(dims(dist), dims(dist))

unsafeMeanCov(dist::ProbabilityDistribution) = (unsafeMean(dist), unsafeCov(dist)) # Can be overloaded for efficiency

unsafeWeightedMeanPrecision(dist::ProbabilityDistribution) = (unsafeWeightedMean(dist), unsafePrecision(dist)) # Can be overloaded for efficiency

isProper(::ProbabilityDistribution{T, PointMass}) where T<:VariateType = true

# Probability distribution parametrized by function
slug(::Type{Function}) = "f"

format(dist::ProbabilityDistribution{V, Function}) where V<:VariateType = "$(dist.params)"

# Distribution constructors
ProbabilityDistribution(::Type{V}, ::Type{Function}; kwargs...) where V<:VariateType = ProbabilityDistribution{V, Function}(Dict{Symbol,Any}(kwargs))
ProbabilityDistribution(::Type{Function}; kwargs...) = ProbabilityDistribution{Univariate, Function}(Dict{Symbol,Any}(kwargs))

unsafeMode(dist::ProbabilityDistribution{T, Function}) where T<:VariateType = deepcopy(dist.params[:mode])

vague(::Type{Function}) = ProbabilityDistribution(Univariate, Function)

isProper(dist::ProbabilityDistribution{<:VariateType, Function}) = haskey(dist.params, :log_pdf)

logPdf(dist::ProbabilityDistribution{<:VariateType, Function}, x) = dist.params[:log_pdf](x)

# Convert VariateTypes
convert(::Type{ProbabilityDistribution{Multivariate, PointMass}}, dist::ProbabilityDistribution{Univariate, PointMass}) =
    ProbabilityDistribution(Multivariate, PointMass, m=[dist.params[:m]])
convert(::Type{ProbabilityDistribution{MatrixVariate, PointMass}}, dist::ProbabilityDistribution{Univariate, PointMass}) =
    ProbabilityDistribution(MatrixVariate, PointMass, m=mat(dist.params[:m]))
convert(::Type{ProbabilityDistribution{MatrixVariate, PointMass}}, dist::ProbabilityDistribution{Multivariate, PointMass}) =
    ProbabilityDistribution(MatrixVariate, PointMass, m=reshape(dist.params[:m], dims(dist), 1))

sample(dist::ProbabilityDistribution{T, PointMass}) where T<:VariateType = deepcopy(dist.params[:m])

"""
Compute conditional differential entropy: H(Y|X) = H(Y, X) - H(X)
"""
conditionalDifferentialEntropy(marg_joint::ProbabilityDistribution{Multivariate}, marg_condition::Vararg{ProbabilityDistribution}) = differentialEntropy(marg_joint) - sum([differentialEntropy(marg) for marg in marg_condition])

"""
This function ensures the argument expression is evaluated at runtime, allowing access to local variables
"""
pack(expr) = expr

function ==(t::ProbabilityDistribution{var_t, fam_t}, u::ProbabilityDistribution{var_u, fam_u}) where {fam_t, fam_u, var_t, var_u}
    (fam_t == fam_u) || return false
    (var_t == var_u) || return false

    # Parameter values are checked for approximate equality
    t.params === u.params && return true
    if length(t.params) != length(u.params) return false end
    for pair in t.params
        if !in(pair, u.params, (a,b)->isapprox(a, b, atol=tiny))
            return false
        end
    end

    return true
end

function gaussianQuadrature(h::Function; m::Float64=0.0, v::Float64=1.0)
    abscissas = [-7.12581, -6.4095, -5.81223, -5.27555, -4.77716, -4.30555, -3.85376, -3.41717, -2.99249, -2.57725, -2.1695, -1.76765, -1.37038, -0.9765, -0.584979, -0.194841, 0.194841, 0.584979, 0.9765, 1.37038, 1.76765, 2.1695, 2.57725, 2.99249, 3.41717, 3.85376, 4.30555, 4.77716, 5.27555, 5.81223, 6.4095, 7.12581]
    weights = [7.31068e-23, 9.23174e-19, 1.19734e-15, 4.21501e-13, 5.93329e-11, 4.09883e-9, 1.57417e-7, 3.65059e-6, 5.41658e-5, 0.000536268, 0.00365489, 0.0175534, 0.0604581, 0.15127, 0.277458, 0.375238, 0.375238, 0.277458, 0.15127, 0.0604581, 0.0175534, 0.00365489, 0.000536268, 5.41658e-5, 3.65059e-6, 1.57417e-7, 4.09883e-9, 5.93329e-11, 4.21501e-13, 1.19734e-15, 9.23174e-19, 7.31068e-23]

    1/sqrt(pi) * sum([weights[i]*h(sqrt(2*v)*abscissas[i] + m) for i=1:32])
end


text:
using ToyFHE
using GaloisFields
using OffsetArrays
using Test

params = BGVParams(
    ToyFHE.CryptParameters.palisade_parameters[4096],
    256,
    8/√(2π)
)
kp = keygen(params)

plain = zero(plaintext_space(params))
plain[0] = 6

c = encrypt(kp, plain)
@test decrypt(kp, c)[0] == 6

let y = c*c
    @test decrypt(kp, y)[0] == 0x24
end


text:
module FEMSolids

using Reexport
@reexport using Ferrite

include("linear_elasticity.jl")
include("element_routine.jl")
include("adaptivity.jl")

export LinearElasticity
export bulk_modulus, shear_modulus
export element_routine!, Primal
export refine, linear_to_quadratic

end


text:
include("header.jl")

# info("Test neural net...")
n0(w,x,y)=sum(abs2,((w[3]*max.(0,w[1]*x.+w[2]).+w[4])-y))
n1 = grad(n0)
n1sum(w,x,y)=sum(map(sum,n1(w,x,y)))
n1sumd(w,x,y)=sum(map(sum,values(n1(w,x,y))))
wa = Any[rand(2,3),rand(2),rand(2,2),rand(2)]
wt = (wa...,)
wd = Dict(); for i=1:length(wa); wd[i]=wa[i]; end

@testset "neuralnet" begin
    @test gradcheck(n0, wa, rand(3,10), rand(2,10))
    @test gradcheck(n0, wt, rand(3,10), rand(2,10))
    @test gradcheck(n0, wd, rand(3,10), rand(2,10))
    @test gradcheck(n1sum, wa, rand(3,10), rand(2,10))
    @test gradcheck(n1sum, wt, rand(3,10), rand(2,10))
    @test gradcheck(n1sumd, wd, rand(3,10), rand(2,10))
end

nothing


text:
module ScalingDimspecTests

using Interpolations, DualNumbers, Base.Test

xs = -pi:(2pi/10):pi-2pi/10
ys = -2:.1:2
f(x,y) = sin(x) * y^2

itp = interpolate(Float64[f(x,y) for x in xs, y in ys], (BSpline(Quadratic(Periodic())), BSpline(Linear())), OnGrid())
sitp = scale(itp, xs, ys)

for (ix,x) in enumerate(xs), (iy,y) in enumerate(ys)
    @test_approx_eq_eps sitp[x,y] f(x,y) sqrt(eps(1.0))

    g = gradient(sitp, x, y)
    fx = epsilon(sitp[dual(x,1), dual(y,0)])
    fy = epsilon(sitp[dual(x,0), dual(y,1)])

    @test_approx_eq_eps g[1] fx sqrt(eps(1.0))
    @test_approx_eq_eps g[2] fy sqrt(eps(1.0))
end

end


text:
type MethodSig
    name::Symbol
    params::Array{Any, 1}
end

@doc doc"""
    A simple macro to automate creating interface code.

    Usage:
        create an interface with
        ```
        @interface MyInterfaceType begin
            func1(self::MyInterfaceType, x::Int)
            func2(self::MyInterfaceType, y)
        end
        ```
""" ->
macro interface(name, prototypes)
    @assert isa(name, Symbol)
    @assert prototypes.head == :block

    typename = esc(name)
    # hidden_typename is just used as a hack to make sure updateing
    # the mutable unions modifies what is excepted by functions that use it
    hidden_typename = esc("_$(name)")
    strname = string(name)
    code = quote
        typealias $(typename) ccall(Libdl.dlsym(MUTABLE_UNION_LIB, :jl_type_mutable_union), Any, ())
        eval(parse(string("type ", $(hidden_typename), " end")))
        ccall(Libdl.dlsym(MUTABLE_UNION_LIB, :jl_type_mutable_union_append),
                Void, (Any, Any), $(typename), $(hidden_typename)
            )
    end

    methods = []
    for line in prototypes.args
        if isa(line, Symbol)
            println("Providing only function names isn't supported yet.")
            continue
        elseif line.head == :line
            continue
        elseif isa(line, Expr)
            fname = line.args[1]
            params = []
            for arg in line.args[2:end]
                push!(params, arg.args[2])
            end
            push!(methods, MethodSig(fname, params))
        end
    end

    func = esc(:(get_methods(interface::Type{$(name)}) = return $(methods)))

    code = quote
        $code
        $func
        Base.show(io::IO, obj::Type{$(typename)}) = print(io, $(strname))
    end

    return code
end

@doc doc"""
    A second macro for making a type implement a declared interface. This is really
    just syntactic sugar for the implements function. The macro takes an
    expression of the forms provided below and calls the implements function

    Usage:
        *Assuming MyType and MyInterface exist and MyType has
        all the required methods defined
        ```
        @implements MyType <: MyInterface
        ```

    or to implement multiple interfaces at once
        ```
        @implements MyType <: (MyInterface, OtherInterface, ...)
        ```
""" ->
macro implements(expr)
   @assert isa(expr, Expr)
   @assert length(expr.args) == 3
   @assert expr.args[2] == :<:

   typename = esc(expr.args[1])
   interfaces = esc(expr.args[3])
   code = quote
       Interfaces.implements($(typename), $(interfaces))
   end
   return code
end

@doc doc"""
    The implements function called from the @implements macro
    takes a single type and either a single Interface or a Tuple of Interfaces.
    The implements function uses methods with to ensure that Type provided implements
    each of the required methods for each Interfaces and then updates the typealias
    union for that Interface with the Type.

    Usage:
        *Assuming MyType and MyInterface exist and MyType has
        all the required methods defined
        ```
        implements(MyType, MyInterface)
        ```

        or for a batch implements
        ```
        implements(MyType, (MyInterface, MyOtherInterface))
        ```
""" ->
function implements(name::Type, interfaces::Union{Type,Tuple})
    interface_names = None
    if isa(interfaces, Type)
        interface_names = (interfaces,)
    elseif isa(interfaces, Tuple)
        interface_names = interfaces
    end

    for i in interface_names
        # Only bother checking the methods and updating the typealias
        # If the interface != type, since methodswith is slow
        # (Note: Foo == Union{Foo})
        if i != name
            methods_exist(i, name, current_module())

            ccall(Libdl.dlsym(MUTABLE_UNION_LIB, :jl_type_mutable_union_append),
                Void, (Any, Any), i, name
            )
        end
    end
end

@doc doc"""
    Determines whether type atype supports all the same methods
    as the interface.
""" ->
function methods_exist(interface::Type, atype::Type, mod)
    get_methods_func = eval(mod, :get_methods)
    interface_type = eval(mod, interface)
    for interface_method in get_methods_func(interface_type)
        fname = interface_method.name
        params = []

        for p in interface_method.params
            param = eval(mod, p)
            if param == interface
                push!(params, atype)
            else
                push!(params, param)
            end
        end

        func = eval(mod, fname)
        if !method_exists(func, tuple(params...))
            error("Required method $(fname)($(params)) not implemented for $(atype)")
        end
    end
end


text:
module Compiler

using GPUCompiler
using LLVM
using LLVM.Interop

import GPUCompiler: CompilerJob, FunctionSpec, codegen

import Enzyme_jll
import Libdl

function __init__()
    Libdl.dlopen(Enzyme_jll.libEnzyme_path, Libdl.RTLD_GLOBAL)
    LLVM.clopts("-enzyme_preopt=0")
end

# Define EnzymeTarget
using LLVM: Target, TargetMachine

Base.@kwdef struct EnzymeTarget <: AbstractCompilerTarget
end
GPUCompiler.llvm_triple(::EnzymeTarget) = Sys.MACHINE

# GPUCompiler.llvm_datalayout(::EnzymeTarget) =  nothing

function GPUCompiler.llvm_machine(target::EnzymeTarget)
    triple = GPUCompiler.llvm_triple(target)
    t = Target(; triple = triple)
    tm = TargetMachine(t, triple)
    LLVM.asm_verbosity!(tm, true)

    return tm
end

module Runtime
    # the runtime library
    signal_exception() = return
    malloc(sz) = Base.Libc.malloc(sz)
    report_oom(sz) = return
    report_exception(ex) = return
    report_exception_name(ex) = return
    report_exception_frame(idx, func, file, line) = return
end

struct EnzymeCompilerParams <: AbstractCompilerParams end

## job

# TODO: We shouldn't blancket opt-out
GPUCompiler.check_invocation(job::CompilerJob{EnzymeTarget}, entry::LLVM.Function) = nothing

GPUCompiler.runtime_module(target::CompilerJob{EnzymeTarget}) = Runtime
GPUCompiler.isintrinsic(::CompilerJob{EnzymeTarget}, fn::String) = true
GPUCompiler.can_throw(::CompilerJob{EnzymeTarget}) = true

# TODO: encode debug build or not in the compiler job
#       https://github.com/JuliaGPU/CUDAnative.jl/issues/368
GPUCompiler.runtime_slug(job::CompilerJob{EnzymeTarget}) = "enzyme"

include("compiler/optimize.jl")
include("compiler/cassette.jl")
include("compiler/validation.jl")

end


text:
using Test
using DataTypesBasic

@test isidentity(Option(5))
@test !isidentity(Option(nothing))
@test isconst(Option(nothing))
@test !isconst(Option(5))

@test iftrue(false) do
  54
end == Option()
@test iftrue(true) do
  54
end == Option(54)
@test iftrue(false, :Something) == Option()
@test iftrue(true, :Something) == Option(:Something)

@test iffalse(false, 54) == Option(54)
@test iffalse(true, 54) == Option()
@test iffalse(() -> "hi", false) == Option("hi")
@test iffalse(() -> "hi", true) == Option()

@test eltype(Option{Int}) == Int
@test eltype(Option) == Any

@test isoption(Identity(3))
@test isoption(Option())
@test !isoption(Const("anything but nothing"))

@test issome(Option(3))
@test !issome(Option())
@test_throws MethodError issome(Const("other"))

@test !isnone(Option(3))
@test isnone(Option())
@test_throws MethodError isnone(Const("other"))


text:
module GuaranteedEllipseFitting

using LinearAlgebra
using Parameters
using Random
using StaticArrays
using Setfield

include("util.jl")
include("types.jl")
include("observations.jl")
include("conversion.jl")
include("sampling.jl")
include("normalise_data_context.jl")
include("direct_ellipse_fit.jl")
include("guaranteed_ellipse_fit.jl")
include("fit_ellipse.jl")
include("covariance.jl")

export Observations,
       UncertainObservations,
       NormaliseDataContext,
       ToNormalisedSpace,
       FromNormalisedSpace,
       DirectEllipseFit,
       GuaranteedEllipseFit,
       AlgebraicToGeometric,
       GeometricToAlgebraic,
       LevenbergMarquardt,
       sample_points_on_ellipse,
       fit_ellipse,
       determine_algebraic_covariance,
       determine_geometric_covariance
end


text:
# MinHash
# ========
#
# Functions to generate MinHash sketches of biological sequences
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md
using DataStructures: SortedSet

"""
MinHash Sketch type

MinHash sketches are a sorted set of kmer hashes, containing the smallest `s`
hashes for kmers of a given length. The type contains two parameters:

* .sketch: a sorted set of hashes
* .kmersize: the length of kmers used to generate the sketch
"""
struct MinHashSketch
    sketch::Vector{UInt64}
    kmersize::Int

    function MinHashSketch(sketch::Vector, kmersize::Integer)
        length(sketch) > 0 || error("Sketch cannot be empty")
        kmersize > 0 || error("Kmersize must be greater than 0")
        return new(sketch, kmersize)
    end
end

Base.getindex(s::MinHashSketch, part) = getindex(s.sketch, part)
Base.length(s::MinHashSketch) = length(s.sketch)

Base.iterate(s::MinHashSketch) = iterate(s.sketch)
Base.iterate(s::MinHashSketch, state) = iterate(s.sketch, state)

function Base.:(==)(a::MinHashSketch, b::MinHashSketch)
    return a.kmersize == b.kmersize && a.sketch == b.sketch
end

function kmerminhash!(::Type{DNAKmer{k}}, seq::BioSequence, s::Integer, kmerhashes::Vector{UInt64}) where {k}
    # generate first `s` kmers
    iter = each(DNAKmer{k}, seq)
    iter_value = iterate(iter)
    while length(kmerhashes) < s && iter_value !== nothing
        (_, kmer), state = iter_value
        iter_value = iterate(iter, state)
        h = hash(canonical(kmer)) # hash lexigraphic minimum of kmer and reverse compliment of kmer
        if h ∉ kmerhashes
            push!(kmerhashes, h)
        end
    end

    sort!(kmerhashes)

    # scan `seq` to make a minhash
    while iter_value !== nothing
        (_, kmer), state = iter_value
        iter_value = iterate(iter, state)
        h = hash(canonical(kmer)) # hash lexigraphic minimum of kmer and reverse compliment of kmer
        if h < kmerhashes[end]
            i = searchsortedlast(kmerhashes, h)
            if i == 0 && h != kmerhashes[1]
                pop!(kmerhashes)
                pushfirst!(kmerhashes, h)
            elseif h != kmerhashes[i]
                pop!(kmerhashes)
                insert!(kmerhashes, i+1, h)
            end
        end
    end

    return kmerhashes
end

"""
    minhash(seq, k::Integer, s::Integer)

Generate a MinHash sketch of size `s` for kmers of length `k`.
"""
function minhash(seq::BioSequence, k::Integer, s::Integer)
    kmerhashes = kmerminhash!(DNAKmer{k}, seq, s, UInt64[])
    length(kmerhashes) < s && error("failed to generate enough hashes")

    return MinHashSketch(kmerhashes, k)
end

function minhash(seqs::Vector{T}, k::Integer, s::Integer) where {T<:BioSequence}
    kmerhashes = UInt64[]
    for seq in seqs
        kmerminhash!(DNAKmer{k}, seq, s, kmerhashes)
    end

    length(kmerhashes) < s && error("failed to generate enough hashes")
    return MinHashSketch(kmerhashes, k)
end

function minhash(seqs::BioCore.IO.AbstractReader, k::Integer, s::Integer)
    kmerhashes = UInt64[]
    for seq in seqs
        kmerminhash!(DNAKmer{k}, sequence(seq), s, kmerhashes)
    end

    length(kmerhashes) < s && error("failed to generate enough hashes")
    return MinHashSketch(kmerhashes, k)
end


text:
# ================================
#   CAKE WAVELETS LIFT
# ================================

# Cut-off function in the frequency plane
mn(ρ::Real; t=1, N = 12) = exp(-ρ^2/t)*sum([(ρ^2/t)^k/factorial(k) for k=0:N])

# B-spline for the radial part
function b_spline(x::T,t1::T,t2::T,t3::T,t4::T) where T<:Real
    if x >= t1 && x<t2
        return (x-t1)^2/(t3-t1)/(t2-t1)
    elseif x >= t2 && x<t3
        return (x-t1)*(t3-x)/(t3-t1)/(t3-t2)+(t4-x)*(x-t2)/(t4-t2)/(t3-t2)
    elseif x >= t3 && x<t4
        return (t4-x)^2/(t4-t2)/(t4-t3)
    else
        return 0
    end
end
b_spline(x) = b_spline(x,-1.,-.5,.5,1.)

# Fourier part of the cake
function cake(ϕ::Real, θs::Int; m::Int = 100, bw::Real = 5, oriented::Bool = false, args...)
    # ϕ = orientation of the cake
    # m = the cake wavelet will be an m x m matrix
    # θs = angular resolution
    # bw = bandwith at which to cut
    @assert m%2 == 0

    x(i) = 2bw*(i-1)/m - bw
    ρ(i,j) = sqrt(x(i)^2 +x(j)^2)
    θ(i,j) = atan(x(j),x(i))
    cake(x) = [mn(ρ(i,j))*b_spline((mod(θ(i,j)-x ,2π)- π/2)*θs/(2*pi)) for i in 1:m, j in 1:m]

    if oriented
        cake(ϕ)
    else
        cake(ϕ) + cake(ϕ+π)
    end
end

function lift(img::Array{T,2}, θs::Int = 30 ; m::Int = 100, args...) where T<:Real
    # θs = angular resolution of the cake wvlts
    # m = resolution of the kernel, if 0 it is chosen as the size of the image
    # (to have perfect reconstruction)
    @assert(m%2 == 0)
    @assert(size(img, 1)%2 == size(img, 2) %2 == 0)

    m == 0 ? m = minimum(size(img)) : nothing

    lift = zeros(Float64, (size(img)...,θs))
    fimg = img |> fft |> fftshift
    bnd1 = Int(ceil((size(img,1)-m)/2))
    bnd2 = Int(ceil((size(img,2)-m)/2))
    NORM = 1.3333333333333384 # NORM = maximum(cake)
    for k in 1:θs
        flift = zeros(Complex{Float64},size(img))
        cake_wavelet = cake((k-1)*π/θs-π/2, θs; m = m, args...)/NORM
        for r in (bnd1+1):(bnd1 + m), h in (bnd2+1):(bnd2 + m)
            @inbounds  flift[r,h] = cake_wavelet[r-bnd1,h-bnd2]*fimg[r,h]
        end
        @inbounds lift[:,:,k] = flift |> ifftshift |> ifft |> real
    end
    return lift
end

max_stars_repo_name: dprn/WCvsLHE

max_issues_repo_name: dprn/WCvsLHE

max_forks_repo_name: dprn/WCvsLHE


text:
using LinearAlgebra
epsilon = 0.00000001 # default, be small anyway

function simplex_method_loop(A, b, c, initial_point, dimension, constraint)
    inter_set = zeros(Float64, (dimension, dimension))
    inter_set_b = zeros(Float64, dimension)
    other_set = zeros(Float64, (constraint - dimension, dimension))
    other_set_b = zeros(Float64, constraint - dimension)
    new_point = zeros(Float64, (dimension, constraint - dimension, dimension))
    new_data = zeros(Float64, (dimension, constraint - dimension))
    new_data = zeros(Float64, (dimension, constraint - dimension))
    new_data[1,1] = 100
    inter = A * initial_point
    j = 0
    l = 0
    instance = 0
    index = 0
    for i = 1:constraint
        if j != dimension
            if abs(inter[i] - b[i]) < epsilon
                j += 1
                inter_set[j,:] = A[i, :]
                inter_set_b[j] = b[i]
            else
                l += 1
                other_set[l,:] = A[i,:]
                other_set_b[l] = b[i]
            end
        else
            index = 1
            instance = i
            break
        end
    end
    if index == 1
        for i = instance:constraint
            l += 1
            other_set[l,:] = A[i,:]
            other_set_b[l] = b[i]
        end
    end
    while maximum(new_data) > epsilon
        for i = 1 : dimension
            for k = 1 : constraint - dimension
                temp = inter_set[i,:]
                temp_b = inter_set_b[i]
                inter_set[i,:] = other_set[k,:]
                inter_set_b[i] = other_set_b[k]
                if det(inter_set) != 0
                    new_point[i,k,:] = inv(inter_set) * inter_set_b
                    new_data[i,k] = c' * (new_point[i,k,:] - initial_point)
                end
                inter_set[i,:] = temp
                inter_set_b[i] = temp_b
            end
        end
        index = argmax(new_data)
        initial_point = new_point[index[1],index[2],:]
        temp = other_set[index[2],:]
        temp_b = other_set_b[index[2]]
        other_set[index[2],:] = inter_set[index[1],:]
        other_set_b[index[2]] = inter_set_b[index[1]]
        inter_set[index[1],:] = temp
        inter_set_b[index[1]] = temp_b
    end
    return initial_point
end

function geometry_find(A, b, c, dimension, constraint)
    initial_point = inv(A[1:dimension,:]) * b[1:dimension]
    abandoned_set = zeros(Float64, (constraint - dimension, dimension))
    abandoned_set_b = zeros(Float64, constraint - dimension)
    for i = dimension + 1:constraint
        if sum(A[i,:] .* initial_point) - b[i] > epsilon
            for j = 1:dimension
                abandoned_set[i - dimension,:] = A[j,:]
                abandoned_set_b[i - dimension] = b[j]
                A[j,:] = A[i,:]
                b[j] = b[i]
                abandon = abandoned_set * inv(A[1:dimension,:]) *  b[1:dimension] - abandoned_set_b
                if maximum(abandon) < epsilon
                    initial_point = inv(A[1:dimension,:]) *  b[1:dimension]
                    break
                else
                    A[j,:] = abandoned_set[i - dimension,:]
                    b[j] = abandoned_set_b[i - dimension]
                    abandoned_set[i - dimension,:] = zeros(Float64, dimension)
                    abandoned_set_b[i - dimension] = 0
                end
            end
        end
    end
    return initial_point
end

dimension = 20 # default
constraint = 100 # default
A = 100 * rand(Float64, (constraint, dimension)) # default
constraint_matrix = A
b = 100 * rand(Float64, constraint) # default
constraint_matrix_b = b
c = 100 * rand(dimension) # default
# objective(x) = c' * x  default
# CAUTION: Ax <= b is the constraint type
# CAUTION: the program is intended to solve Maximium
# CAUTION: you should ensure the affine independence of your constraint
initial_point = geometry_find(A, b, c, dimension, constraint)

# print(initial_point)
optimal_point = simplex_method_loop(constraint_matrix, constraint_matrix_b, c, initial_point, dimension, constraint)
# print(optimal_point)|

max_stars_repo_name: zhuhs-q/basic-optimal-method

max_issues_repo_name: zhuhs-q/basic-optimal-method

max_forks_repo_name: zhuhs-q/basic-optimal-method


text:
## File that runs the entire cross-val analysis and dumps it for the traffic stuff

#include("fitTraffic.jl")
#using PyPlot

##########
#Generate the simulated data
##########
numData = 1; sigma = .05
flow_data = Array(Float64, length(arcs), numData)
demand_data = Dict{(Int, Int), Array{Float64, 1}}()

numNodes = maximum(map(pair->pair[1], keys(demands)))
g = simple_inclist(numNodes, is_directed=true)
vArcs = Arc[]
for arc in values(arcs)
    arc.obsflow = arc.trueflow
    add_edge!(g, arc.initNode, arc.termNode)
    push!(vArcs, arc)
end

#srand(8675309)
for iRun = 1:numData
    #perturb the demand_data
    for odpair in keys(demands)
        if ! haskey(demand_data, odpair)
            demand_data[odpair] = [demands[odpair] * (1 + sigma * rand()), ]
        else
            push!(demand_data[odpair], demands[odpair] * (1  + sigma * rand()))
        end
    end
    #solve using FW and record
    conv_tol, flow_data[:, iRun] = frank_wolfe(g, vArcs, demand_data, iRun)
    println(conv_tol)
end

#Randomzie the flow data a little bit too
for i = 1:size(flow_data, 1)
    for j = 1:size(flow_data, 2)
        flow_data[i, j] *= (1 + sigma * rand() )
    end
end

max_stars_repo_name: jingzbu/InverseVITraffic

max_issues_repo_name: jingzbu/InverseVITraffic

max_forks_repo_name: jingzbu/InverseVITraffic


text:
"""
Generate a new action when the set of actions is widened.
"""
function next_action end

mutable struct RandomActionGenerator{RNG<:AbstractRNG}
    rng::RNG
end
RandomActionGenerator() = RandomActionGenerator(Random.GLOBAL_RNG)

function next_action(gen::RandomActionGenerator, mdp::Union{POMDP,MDP}, s, snode::AbstractStateNode)
    rand(gen.rng, actions(mdp, s))
end
next_action(f::Function, mdp::Union{POMDP,MDP}, s, snode::AbstractStateNode) = f(mdp, s, snode)


text:
module Example

Base.@ccallable empty()::Cint = 10
Base.@ccallable foo(x::Cint)::Cint = 3
Base.@ccallable function bar(x::Cint)::Cint
    return 3
end

Base.@ccallable function just_print()::Cvoid
    print("Example")
end

Base.@ccallable function bar_segfault(x::Cint)::Cint
    println("Example")
    return 3
end

Base.@ccallable function example(x::Cint)::Cint
    return x
end
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# transpose
function ad_transpose(x::AD)
  x.value = (x.value)'
  x.grad = (x.grad)'
  return x
end
transpose(x::AD) = ad_transpose(x)

# addition
function ad_add(x::AD, y::AD)
  gradSize = (size(x.value) > size(y.value))? size(x.value):size(y.value)
  result = AD("+", x.value + y.value, true, zeros(gradSize), ad_addD, (x.level > y.level ? x.level:y.level)+1)
  push!(result.parents, x)
  push!(result.parents, y)
  return result
end
function ad_addD(this::AD)
  prevGradSize = size(this.grad)
  adNode1Size = size(this.parents[1].value)
  adNode2Size = size(this.parents[2].value)

  if adNode1Size == prevGradSize
    this.parents[1].grad += 1 * this.grad
  elseif adNode1Size[1] == prevGradSize[1]
    this.parents[1].grad += 1 * sum(this.grad,2)
  else
    this.parents[1].grad += 1 * sum(this.grad,1)
  end

  if adNode2Size == prevGradSize
    this.parents[2].grad += 1 * this.grad
  elseif adNode2Size[1] == prevGradSize[1]
    this.parents[2].grad += 1 * sum(this.grad,2)
  else
    this.parents[2].grad += 1 * sum(this.grad,1)
  end
  return
end
+(x::AD, y::AD) = ad_add(x, y)

function ad_eladd(x::AD, y::AD)
  gradSize = (size(x.value) > size(y.value))? size(x.value):size(y.value)
  result = AD(".+", x.value .+ y.value, true, zeros(gradSize), ad_eladdD, (x.level > y.level ? x.level:y.level)+1)
  push!(result.parents, x)
  push!(result.parents, y)
  return result
end
function ad_eladdD(this::AD)
  prevGradSize = size(this.grad)
  adNode1Size = size(this.parents[1].value)
  adNode2Size = size(this.parents[2].value)

  if adNode1Size == prevGradSize
    this.parents[1].grad += 1 * this.grad
  elseif adNode1Size[1] == prevGradSize[1]
    this.parents[1].grad += 1 * sum(this.grad,2)
  else
    this.parents[1].grad += 1 * sum(this.grad,1)
  end

  if adNode2Size == prevGradSize
    this.parents[2].grad += 1 * this.grad
  elseif adNode2Size[1] == prevGradSize[1]
    this.parents[2].grad += 1 * sum(this.grad,2)
  else
    this.parents[2].grad += 1 * sum(this.grad,1)
  end

  return
end
.+(x::AD, y::AD) = ad_eladd(x, y)

# multiplication
function ad_mul(x::AD, y::AD)
  if size(x.value) == (1,1) || size(y.value) == (1,1)
    return ad_elmul(x, y)
  else
    resultValue = x.value * y.value
    result = AD("*", resultValue, true, zeros(size(resultValue)), ad_mulD, (x.level > y.level ? x.level:y.level)+1)
    push!(result.parents, x)
    push!(result.parents, y)
    return result
  end
end
function ad_mulD(this::AD)
  rowNode1, colNode1 = size(this.parents[1].value)
  rowNode2, colNode2 = size(this.parents[2].value)

  for i = 1:rowNode1, j = 1:colNode2
    pd = this.grad[i,j]
    for k = 1:colNode1
      this.parents[1].grad[i,k] += this.parents[2].value[k,j] * pd
      this.parents[2].grad[k,j] += this.parents[1].value[i,k] * pd
    end
  end
end
*(x::AD, y::AD) = ad_mul(x, y)

# element-wise multiplication
function ad_elmul(x::AD, y::AD)
  gradSize = (size(x.value) > size(y.value))? size(x.value):size(y.value)
  result = AD(".*", x.value .* y.value, true, zeros(gradSize), ad_elmulD, (x.level > y.level ? x.level:y.level)+1)
  push!(result.parents, x)
  push!(result.parents, y)
  return result
end
function ad_elmulD(this::AD)
  this.parents[1].grad += this.parents[2].value .* this.grad
  this.parents[2].grad += this.parents[1].value .* this.grad
  return
end
.*(x::AD, y::AD) = ad_elmul(x, y)

# relu
function reluMat(z::Matrix{Float64})
    return max(0.0, z);
end
function relu(x::AD)
  result = AD("relu", reluMat(x.value), true, zeros(size(x.value)), ad_reluD, x.level + 1)
  push!(result.parents, x)
  return result
end
function reluGradientMat(z::Matrix{Float64})
    grad = ones(z);
    grad[z.<=0] = 0;
    return grad;
end
function ad_reluD(this::AD)
  this.parents[1].grad += reluGradientMat(this.parents[1].value) .* this.grad
  return
end

# sigmoid
function sigmoidMat(z::Matrix{Float64})
    g = 1.0 ./ (1.0 + exp(-z));
    return g;
end
function sigmoid(x::AD)
  result = AD("sigmoid", sigmoidMat(x.value), true, zeros(size(x.value)), ad_sigmoidD, x.level + 1)
  push!(result.parents, x)
  return result
end
function sigmoidGradientMat(z::Matrix{Float64})
  return z.*(1.0 - z);
end
function ad_sigmoidD(this::AD)
  this.parents[1].grad += sigmoidGradientMat(this.parents[1].value) .* this.grad
  return
end

# tanh
function tanh(x::AD)
  result = AD("tanh", tanh(x.value), true, zeros(size(x.value)), ad_tanhD, x.level + 1)
  push!(result.parents, x)
  return result
end
function tanhGradientMat(z::Matrix{Float64})
  return (1.0 - z^2)
end
function ad_tanhD(this::AD)
  this.parents[1].grad += tanhGradientMat(this.parents[1].value) .* this.grad
  return
end

# softmax
function softmaxMat(z::Matrix{Float64}, dimOrder::Int64)
    g = exp(z)
    g ./=sum(g, dimOrder)
    return g;
end

# Using cross-entropy loss function.
# evaluating d(loss)/d(loss)*d(loss)/d(softmax)*d(softmax)/d(z)
# because calculating d(softmax)/d(z) individually requires too much memory.
function softmaxLoss(y::Matrix{Float64}, x::AD, dimOrder::Int64, computeLoss::Bool, reg::Float64)
   # dimOrder: along which the dimentions/classes are present
   probs = softmaxMat(x.value, dimOrder)
   # compute the loss: average cross-entropy loss and regularization
   loss = 0.0
   if computeLoss
     dataInstanceOrder = (dimOrder == 1? 2 : 1)
     numInstance = size(y, dataInstanceOrder)
     correctLogProbs = zeros(size(y))
     for j in 1:numInstance
       correctLogProbs[j, :] = -log(probs[j,Int64(dataInstanceOrder==1?y[j, 1]:y[1,j])]);
     end
     loss = sum(correctLogProbs)/numInstance + reg
   end

   # return node
   result = AD("softmaxLoss", probs, true, zeros(size(x.value)), ad_softmaxLossD, x.level + 1, y, dataInstanceOrder, loss)
   push!(result.parents, x)
   return result
end

softmaxLoss(y::Matrix{Float64}, x::AD, dimOrder::Int64, computeLoss::Bool) = softmaxLoss(y, x, dimOrder, computeLoss, 0.0)
softmaxLoss(y::Matrix{Float64}, x::AD, dimOrder::Int64) = softmaxLoss(y, x, dimOrder, false, 0.0)

function ad_softmaxLossD(this::AD)
  this.parents[1].grad += this.value .* this.grad
  if this.outputOrder == 1
    for j in 1:size(this.parents[1].grad, this.outputOrder)
      this.parents[1].grad[j,Int64(this.output[j,1])] -= 1;
    end
  elseif this.outputOrder == 2
    for j in 1:size(this.parents[1].grad, this.outputOrder)
      this.parents[1].grad[Int64(this.output[1,j]), j] -= 1;
    end
  end
  this.parents[1].grad /= size(this.output, this.outputOrder);
  return
end


text:
using ReinforcementLearningCore

using AbstractTrees
using TensorBoardLogger: TBLogger

RLCore.is_expand(::TBLogger) = false

"""
    EnrichedAction(action;kwargs...)

Inject some runtime info into the action
"""
struct EnrichedAction{A,M}
    action::A
    meta::M
end

EnrichedAction(action;kwargs...) = EnrichedAction(action, kwargs.data)

(env::AbstractEnv)(action::EnrichedAction) = env(action.action)


text:
# Can we figure out how to get rid of this?
(p::DiffEqBase.DefaultLinSolve)(x, A::ComponentMatrix, b, update_matrix=false; tol=nothing, kwargs...) = p(x, getdata(A), b, update_matrix; tol, kwargs...)

DiffEqBase.diffeqbc(x::ComponentArray) = DiffEqBase.DiffEqBC(x)


text:
module Var

using Random: AbstractRNG
using ..Basis

"""
`recurse(f, ω)`
Recursively apply contextual execution to internals of `f`"""
function recurse end

include("variable.jl")          # Random / Parametric Variables
include("multivariate.jl")      # Multivariate Distributions
# include("typevar.jl")           # Type Variables
include("member.jl")            # Families
include("primparam.jl")         # Primitive Parameters
include("primdist.jl")          # Primitive Distributions
include("pointwise.jl")         # Point wise variable application
include("dispatch.jl")          # Contextual application

end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
@testset "LaguerreFunction" begin
    @testset "LaguerreFunction construction" begin
        # test construction of LaguerreFunction
        orders = [0,2,3]
        coeffs = [1.0, 2, 0.5]
        L = LaguerreFunction(orders, coeffs)
        @test L.orders == orders
        @test all(isapprox.(coeffs, L.coeffs))
        # test single Laguerre Function
        L = LaguerreFunction(2, 1 // 2)
        @test L.orders == [2]
        @test L.coeffs == [1 // 2]

        # Laguerre Function L3 with coefficient 1.
        L = LaguerreFunction(3)
        @test L.orders == [3]
        @test L.coeffs == [1]
        # test throw
        @test_throws DimensionMismatch LaguerreFunction([1,2], [2,3,2])
        @test_throws DomainError LaguerreFunction(-2, 1)
        # test string representation
        L = LaguerreFunction([2,3], [1,2])
        st = Laguerre.lashow(L)
        @test st == "(1)Lf2 + (2)Lf3"
    end

    @testset "LaguerreFunction addition" begin
        L1 = LaguerreFunction(2, 1 // 2)
        L2 = LaguerreFunction(3, 1.0)
        L = L1 + L2
        @test typeof(L) == LaguerreFunction{typeof(1.0)}
        @test L.orders == [2, 3]
        @test all(isapprox.(L.coeffs, [1 / 2, 1.0]))
        #
        L1 = LaguerreFunction(2, 1)
        L2 = LaguerreFunction(2, 1 // 3)
        L = L1 + L2
        @test typeof(L) == LaguerreFunction{typeof(1 // 3)}
        @test L.orders == [2]
        @test L.coeffs == [1 + 1 // 3]

        # test subtraction
        L = L1 - L2
        @test L.orders == [2]
        @test L.coeffs == [1 - 1 // 3]

    end

    @testset "LaguerreFunction constant operation" begin
        orders = [0,2,3]
        coeffs = [1, 2, 2]
        L1 = LaguerreFunction(orders, coeffs)
        L2 = L1 * 3
        # multiplication by a number
        @test L2.orders == [0,2,3]
        @test L2.coeffs == 3 * coeffs
        # division by a number
        @test_throws ErrorException L1 / 0
        L3 = L1 / 2
        @test L3.orders == L1.orders
        @test L3.coeffs == L1.coeffs / 2
    end

    @testset "Laguerre simplify" begin
        L1 = LaguerrePolynomial([2, 3], [1,1])
        L2 = LaguerrePolynomial(2, 1)
        res = LaguerrePolynomial(3, 1)
        @test simplify(L1 - L2) == res
        @test simplify(L1 - L1) == LaguerrePolynomial(0, 0)

        L1 = LaguerreFunction([2, 3], [1.0,1.0])
        L2 = LaguerreFunction(2, 1.0)
        res = LaguerreFunction(3, 1.0)
        @test simplify(L1 - L2) == res
        @test simplify(L1 - L1) == LaguerreFunction(0, 0)
    end
end


text:
__precompile__()
module MathOptInterfaceGLPK

import Base.copy

export GLPKOptimizerLP, GLPKOptimizerMIP

using GLPK
using MathOptInterface
const MOI = MathOptInterface
using LinQuadOptInterface
const LQOI = LinQuadOptInterface

# Many functions in this module are adapted from GLPKMathProgInterface.jl. This is the copyright notice:
## Copyright (c) 2013: Carlo Baldassi
## Permission is hereby granted, free of charge, to any person obtaining a copy
## of this software and associated documentation files (the "Software"), to deal
## in the Software without restriction, including without limitation the rights
## to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
## copies of the Software, and to permit persons to whom the Software is
## furnished to do so, subject to the following conditions:
## The above copyright notice and this permission notice shall be included in
## all copies or substantial portions of the Software.
## THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
## IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
## FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
## AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
## LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
## OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
## SOFTWARE.

const SUPPORTED_OBJECTIVES = [
    LQOI.Linear
]
const SUPPORTED_CONSTRAINTS_LP = [
    (LQOI.Linear, LQOI.EQ),
    (LQOI.Linear, LQOI.LE),
    (LQOI.Linear, LQOI.GE),
    (LQOI.Linear, LQOI.IV),
    (LQOI.SinVar, LQOI.EQ),
    (LQOI.SinVar, LQOI.LE),
    (LQOI.SinVar, LQOI.GE),
    (LQOI.SinVar, LQOI.IV),
    (LQOI.VecVar, MOI.Nonnegatives),
    (LQOI.VecVar, MOI.Nonpositives),
    (LQOI.VecVar, MOI.Zeros),
    (LQOI.VecLin, MOI.Nonnegatives),
    (LQOI.VecLin, MOI.Nonpositives),
    (LQOI.VecLin, MOI.Zeros)
]
const SUPPORTED_CONSTRAINTS_MIP = vcat(SUPPORTED_CONSTRAINTS_LP, [
    (LQOI.SinVar, MOI.ZeroOne),
    (LQOI.SinVar, MOI.Integer),
    # (VecVar, LQOI.SOS1),
    # (VecVar, LQOI.SOS2),
])

abstract type GLPKOptimizer <: LQOI.LinQuadOptimizer end

import GLPK.Prob

const Model = GLPK.Prob
Model(env) = Model()

LQOI.LinQuadModel(::Type{M},env) where M<:GLPKOptimizer = GLPK.Prob()

mutable struct GLPKOptimizerLP <: GLPKOptimizer
    LQOI.@LinQuadOptimizerBase
    method::Symbol
    param::Union{GLPK.SimplexParam, GLPK.InteriorParam}
    GLPKOptimizerLP(::Void) = new()
end
function GLPKOptimizerLP(presolve = false, method = :Simplex;kwargs...)
    if !(method in [:Simplex,:Exact,:InteriorPoint])
        error("""
        Unknown method for GLPK LP solver: $method
            Allowed methods:
                :Simplex
                :Exact
                :InteriorPoint""")
    end
    if method == :Simplex || method == :Exact
        param = GLPK.SimplexParam()
        if presolve
            param.presolve = GLPK.ON
        end
    elseif method == :InteriorPoint
        param = GLPK.InteriorParam()
        if presolve
            warn("Ignored option: presolve")
        end
    else
        error("This is a bug")
    end
    param.msg_lev = GLPK.MSG_ERR
    for (k,v) in kwargs
        i = findfirst(x->x==k, fieldnames(typeof(param)))
        if i > 0
            t = typeof(param).types[i]
            setfield!(param, i, convert(t, v))
        else
            warn("Ignored option: $(string(k))")
        end
    end

    m = GLPKOptimizerLP(nothing)
    m.param = param
    m.method = method
    MOI.empty!(m)

    return m
end

mutable struct GLPKOptimizerMIP <: GLPKOptimizer

    LQOI.@LinQuadOptimizerBase

    param::GLPK.IntoptParam
    smplxparam::GLPK.SimplexParam
    # lazycb::Union{Function,Void}
    # cutcb::Union{Function,Void}
    # heuristiccb::Union{Function,Void}
    # infocb::Union{Function,Void}
    objbound::Float64
    # cbdata::MathProgCallbackData
    binaries::Vector{Int}
    userlimit::Bool

    GLPKOptimizerMIP(::Void) = new()
end

function GLPKOptimizerMIP(presolve = false; kwargs...)

    env = nothing
    lpm = GLPKOptimizerMIP(nothing)
    lpm.param = GLPK.IntoptParam()
    lpm.smplxparam = GLPK.SimplexParam()
    lpm.objbound = -Inf
    lpm.binaries = Int[]
    lpm.userlimit = false
    MOI.empty!(lpm)

    lpm.param.msg_lev = GLPK.MSG_ERR
    lpm.smplxparam.msg_lev = GLPK.MSG_ERR
    if presolve
        lpm.param.presolve = GLPK.ON
    end

    # lpm.param.cb_func = cfunction(_internal_callback, Void, (Ptr{Void}, Ptr{Void}))
    # lpm.param.cb_info = pointer_from_objref(lpm.cbdata)

    for (k,v) in kwargs
        if k in [:cb_func, :cb_info]
            warn("ignored option: $(string(k)); use the MathProgBase callback interface instead")
            continue
        end
        i = findfirst(x->x==k, fieldnames(typeof(lpm.param)))
        s = findfirst(x->x==k, fieldnames(typeof(lpm.smplxparam)))
        if !(i > 0 || s > 0)
            warn("Ignored option: $(string(k))")
            continue
        end
        if i > 0
            t = typeof(lpm.param).types[i]
            setfield!(lpm.param, i, convert(t, v))
        end
        if s > 0
            t = typeof(lpm.smplxparam).types[s]
            setfield!(lpm.smplxparam, s, convert(t, v))
        end
    end

    return lpm
end

function MOI.empty!(m::GLPKOptimizer)
    MOI.empty!(m,nothing)
end

LQOI.lqs_supported_objectives(s::GLPKOptimizer) = SUPPORTED_OBJECTIVES
LQOI.lqs_supported_constraints(s::GLPKOptimizerMIP) = SUPPORTED_CONSTRAINTS_MIP
LQOI.lqs_supported_constraints(s::GLPKOptimizerLP) = SUPPORTED_CONSTRAINTS_LP
#=
    inner wrapper
=#

#=
    Main
=#

# LinQuadSolver # Abstract type
# done above

# LQOI.lqs_setparam!(env, name, val)
# TODO fix this one
# LQOI.lqs_setparam!(m::GLPKOptimizer, name, val) = GLPK.setfield!(m.inner, string(name), val)

# LQOI.lqs_setlogfile!(env, path)
# TODO fix this one
# LQOI.lqs_setlogfile!(m::GLPKOptimizer, path) = GLPK.setlogfile(m.inner, path::String)

# LQOI.lqs_getprobtype(m)
# TODO - consider removing, apparently useless

#=
    Constraints
=#

# LQOI.lqs_chgbds!(m, colvec, valvec, sensevec)
function LQOI.lqs_chgbds!(instance::GLPKOptimizer, colvec, valvec, sensevec)
    m = instance.inner
    for i in eachindex(colvec)
        colub = Inf
        collb = -Inf
        bt = GLPK.DB
        if sensevec[i] == Cchar('E')
            colub = valvec[i]
            collb = valvec[i]
            bt = GLPK.FX
        elseif sensevec[i] == Cchar('L')
            collb = valvec[i]
            colub = Inf
            u = GLPK.get_col_ub(m, colvec[i])
            if u < 1e100
                bt = GLPK.DB
                colub = u
            else
                bt = GLPK.LO
            end
        elseif sensevec[i] == Cchar('U')
            colub = valvec[i]
            collb = -Inf
            l = GLPK.get_col_lb(m, colvec[i])
            if l > -1e100
                bt = GLPK.DB
                collb = l
            else
                bt = GLPK.UP
            end
        else
            error("invalid bound type")
        end
        if colub > 1e100 && collb < -1e100
            bt = GLPK.FR
            colub = Inf
            collb = -Inf
        elseif colub == collb
            bt = GLPK.FX
        end
        GLPK.set_col_bnds(m, colvec[i], bt, collb, colub)
    end

end

# LQOI.lqs_getlb(m, col)
LQOI.lqs_getlb(instance::GLPKOptimizer, col) = GLPK.get_col_lb(instance.inner, col)

# LQOI.lqs_getub(m, col)
LQOI.lqs_getub(instance::GLPKOptimizer, col) = GLPK.get_col_ub(instance.inner, col)

# LQOI.lqs_getnumrows(m)
LQOI.lqs_getnumrows(instance::GLPKOptimizer) = GLPK.get_num_rows(instance.inner)

# LQOI.lqs_addrows!(m, rowvec, colvec, coefvec, sensevec, rhsvec)
function LQOI.lqs_addrows!(instance::GLPKOptimizer, rowvec, colvec, coefvec, sensevec, rhsvec)
    m = instance.inner

    nrows = length(rhsvec)

    if nrows <= 0
        error("nor row to be added")
    elseif nrows == 1
        addrow!(m, colvec::Vector, coefvec::Vector, sensevec[1], rhsvec[1])

    else

        push!(rowvec, length(colvec)+1)

        for i in 1:(length(rowvec)-1)
            inds = colvec[rowvec[i]:rowvec[i+1]-1]
            coefs = coefvec[rowvec[i]:rowvec[i+1]-1]
            addrow!(m, inds, coefs, sensevec[i], rhsvec[i])
        end

    end
    nothing
end
function addrow!(lp::GLPK.Prob, colidx::Vector, colcoef::Vector, sense::Cchar, rhs::Real)
    if length(colidx) != length(colcoef)
        error("colidx and colcoef have different legths")
    end
    GLPK.add_rows(lp, 1)
    m = GLPK.get_num_rows(lp)
    GLPK.set_mat_row(lp, m, colidx, colcoef)

    if sense == Cchar('E')
        bt = GLPK.FX
        rowlb = rhs
        rowub = rhs
    elseif sense == Cchar('G')
        bt = GLPK.LO
        rowlb = rhs
        rowub = Inf
    elseif sense == Cchar('L')
        bt = GLPK.UP
        rowlb = -Inf
        rowub = rhs
    elseif sense == Cchar('R')
        # start with lower
        bt = GLPK.DB
        rowlb = rhs
        rowub = Inf
    else
        error("row type $(sense) not valid")
        bt = GLPK.FR
    end
    GLPK.set_row_bnds(lp, m, bt, rowlb, rowub)
    return
end

function setrhs(instance::GLPKOptimizer, idx::Integer, rhs::Real)
    lp = instance.inner

    l = GLPK.get_row_lb(lp, idx)
    u = GLPK.get_row_ub(lp, idx)

    rowub = Inf
    rowlb = -Inf

    if l == u
        bt = GLPK.FX
        rowlb = rhs
        rowub = rhs
    elseif l > -Inf && u < Inf
        bt = GLPK.FX
        rowlb = rhs
        rowub = rhs
    elseif l > -Inf
        bt = GLPK.LO
        rowlb = rhs
        rowub = Inf
    elseif u < Inf
        bt = GLPK.UP
        rowlb = -Inf
        rowub = rhs
    else
        error("not valid rhs")
    end

    GLPK.set_row_bnds(lp, idx, bt, rowlb, rowub)
end
# LQOI.lqs_getrhs(m, rowvec)
function LQOI.lqs_getrhs(instance::GLPKOptimizer, row)
    m = instance.inner
    sense = GLPK.get_row_type(m, row)
    if sense == GLPK.LO
        return GLPK.get_row_lb(m, row)
    elseif sense == GLPK.FX
        return GLPK.get_row_lb(m, row)
    elseif sense == GLPK.DB
        return GLPK.get_row_lb(m, row)
    else
        return GLPK.get_row_ub(m, row)
    end
end
# colvec, coef = LQOI.lqs_getrows(m, rowvec)
# TODO improve
function LQOI.lqs_getrows(instance::GLPKOptimizer, idx)
    lp = instance.inner
    colidx, coefs = GLPK.get_mat_row(lp, idx)
    return colidx-1, coefs
end

# LQOI.lqs_getcoef(m, row, col) #??
# TODO improve
function LQOI.lqs_getcoef(instance::GLPKOptimizer, row, col)
    lp = instance.inner
    colidx, coefs = GLPK.get_mat_row(lp, row)
    idx = findfirst(colidx, col)
    if idx > 0
        return coefs[idx]
    else
        return 0.0
    end
end

# LQOI.lqs_chgcoef!(m, row, col, coef)
function LQOI.lqs_chgcoef!(instance::GLPKOptimizer, row, col, coef)
    if row == 0
        lp = instance.inner
        GLPK.set_obj_coef(lp, col, coef)
    elseif col == 0
        setrhs(instance, row, coef)
    else
        chgmatcoef!(instance, row, col, coef)
    end
    return nothing
end

function chgmatcoef!(instance::GLPKOptimizer, row, col, coef)
    lp = instance.inner
    colidx, coefs = GLPK.get_mat_row(lp, row)
    idx = findfirst(colidx, col)
    if idx > 0
        coefs[idx] = coef
    else
        push!(colidx, col)
        push!(coefs, coef)
    end
    GLPK.set_mat_row(lp, row, colidx, coefs)
    return nothing
end
# LQOI.lqs_delrows!(m, row, row)
function LQOI.lqs_delrows!(instance::GLPKOptimizer, rowbeg, rowend)

    m = instance.inner

    idx = collect(rowbeg:rowend)

    GLPK.std_basis(m)
    GLPK.del_rows(m, length(idx), idx)

    nothing
end
# LQOI.lqs_chgctype!(m, colvec, typevec)
# TODO fix types
function LQOI.lqs_chgctype!(instance::GLPKOptimizer, colvec, vartype)

    lp = instance.inner
    coltype = GLPK.CV
    for i in eachindex(colvec)
        if vartype[i] == Cint('I')
            coltype = GLPK.IV
        elseif vartype[i] == Cint('C')
            coltype = GLPK.CV
        elseif vartype[i] == Cint('B')
            coltype = GLPK.IV
            GLPK.set_col_bnds(lp, colvec[i], GLPK.DB, 0.0, 1.0)
        else
            error("invalid variable type: $(vartype[i])")
        end
        GLPK.set_col_kind(lp, colvec[i], coltype)
    end

end
# LQOI.lqs_chgsense!(m, rowvec, sensevec)
# TODO fix types
function LQOI.lqs_chgsense!(instance::GLPKOptimizer, rowvec, sensevec)
    for i in eachindex(rowvec)
        changesense!(instance, rowvec[i], sensevec[i])
    end
    nothing
end
function changesense!(instance::GLPKOptimizer, row, sense)
    m = instance.inner
    oldsense = GLPK.get_row_type(m, row)
    newsense = translatesense(sense)
    if oldsense == newsense
        return nothing
    end

    if newsense == GLPK.FX
        rowub = rowlb
    elseif oldsense == GLPK.DB
        if newsense == GLPK.UP
            rowlb = -Inf
        else
            rowub = Inf
        end
    else
        rowlb = get_row_lb(m, row)
        rowub = get_row_ub(m, row)
        if newsense == GLPK.UP
            rowub = rowlb
            rowlb = -Inf
        else
            rowlb = rowub
            rowub = Inf
        end
    end

    GLPK.set_row_bnds(m, row, newsense, rowlb, rowub)

    nothing
end
function translatesense(sense)
    if sense == Cchar('E')
        return GLPK.FX
    elseif sense == Cchar('R')
        return GLPK.DB
    elseif sense == Cchar('L')
        return GLPK.UP
    elseif sense == Cchar('G')
        return GLPK.LO
    else
        error("invalid sense")
    end
end

const VAR_TYPE_MAP = Dict{Symbol,Cchar}(
    :CONTINUOUS => Cchar('C'),
    :INTEGER => Cchar('I'),
    :BINARY => Cchar('B')
)
LQOI.lqs_vartype_map(m::GLPKOptimizer) = VAR_TYPE_MAP

# LQOI.lqs_addsos(m, colvec, valvec, typ)
LQOI.lqs_addsos!(instance::GLPKOptimizer, colvec, valvec, typ) = GLPK.add_sos!(instance.inner, typ, colvec, valvec)
# LQOI.lqs_delsos(m, idx, idx)
LQOI.lqs_delsos!(instance::GLPKOptimizer, idx1, idx2) = error("cant del SOS")

const SOS_TYPE_MAP = Dict{Symbol,Symbol}(
    :SOS1 => :SOS1,#Cchar('1'),
    :SOS2 => :SOS2#Cchar('2')
)
LQOI.lqs_sertype_map(m::GLPKOptimizer) = SOS_TYPE_MAP

# LQOI.lqs_getsos(m, idx)
# TODO improve getting processes
function LQOI.lqs_getsos(instance::GLPKOptimizer, idx)
    indices, weights, types = GLPK.getsos(instance.inner, idx)

    return indices, weights, types == Cchar('1') ? :SOS1 : :SOS2
end
# LQOI.lqs_getnumqconstrs(m)
LQOI.lqs_getnumqconstrs(instance::GLPKOptimizer) = error("GLPK does not support quadratic ocnstraints")

# LQOI.lqs_addqconstr(m, cols,coefs,rhs,sense, I,J,V)
LQOI.lqs_addqconstr!(instance::GLPKOptimizer, cols,coefs,rhs,sense, I,J,V) = error("GLPK does not support quadratic ocnstraints")

# LQOI.lqs_chgrngval
function LQOI.lqs_chgrngval!(instance::GLPKOptimizer, rows, vals)
    lp = instance.inner
    for i in eachindex(rows)
        l = GLPK.get_row_lb(lp, rows[i])
        GLPK.set_row_bnds(lp, rows[i], GLPK.DB, l, l+vals[i])
    end
    nothing
end
const CTR_TYPE_MAP = Dict{Symbol,Cchar}(
    :RANGE => Cchar('R'),
    :LOWER => Cchar('G'),
    :UPPER => Cchar('L'),
    :EQUALITY => Cchar('E')
)
LQOI.lqs_ctrtype_map(m::GLPKOptimizer) = CTR_TYPE_MAP

#=
    Objective
=#

# LQOI.lqs_copyquad(m, intvec,intvec, floatvec) #?
LQOI.lqs_copyquad!(instance::GLPKOptimizer, I, J, V) = error("GLPK does no support quadratics")

# LQOI.lqs_chgobj(m, colvec,coefvec)
function LQOI.lqs_chgobj!(instance::GLPKOptimizer, colvec, coefvec)
    ncols = GLPK.get_num_cols(instance.inner)
    new_colvec = collect(1:ncols)
    new_coefvec = zeros(ncols)
    for (ind,val) in enumerate(colvec)
        new_coefvec[val] = coefvec[ind]
    end
    m = instance.inner
    for i in eachindex(new_colvec)
        GLPK.set_obj_coef(m, new_colvec[i], new_coefvec[i])
    end
    nothing
end

# LQOI.lqs_chgobjsen(m, symbol)
# TODO improve min max names
function LQOI.lqs_chgobjsen!(instance::GLPKOptimizer, sense)
    m = instance.inner
    if sense == :Min
        GLPK.set_obj_dir(m, GLPK.MIN)
    elseif sense == :Max
        GLPK.set_obj_dir(m, GLPK.MAX)
    else
        error("Unrecognized objective sense $sense")
    end
end
# LQOI.lqs_getobj(m)
function LQOI.lqs_getobj(instance::GLPKOptimizer)
    m = instance.inner
    n = GLPK.get_num_cols(m)
    obj = Array{Float64}(n)
    for c = 1:n
        l = GLPK.get_obj_coef(m, c)
        obj[c] = l
    end
    return obj
end

# lqs_getobjsen(m)
function LQOI.lqs_getobjsen(instance::GLPKOptimizer)

    s = GLPK.get_obj_dir(instance.inner)
    if s == GLPK.MIN
        return MOI.MinSense
    elseif s == GLPK.MAX
        return MOI.MaxSense
    else
        error("Internal library error")
    end
end

#=
    Variables
=#

# LQOI.lqs_getnumcols(m)
LQOI.lqs_getnumcols(instance::GLPKOptimizer) = GLPK.get_num_cols(instance.inner)

# LQOI.lqs_newcols!(m, int)
function LQOI.lqs_newcols!(instance::GLPKOptimizer, int)
    n = GLPK.get_num_cols(instance.inner)
    GLPK.add_cols(instance.inner, int)
    for i in 1:int
        GLPK.set_col_bnds(instance.inner, n+i, GLPK.FR, -Inf, Inf)
    end
    nothing
end

# LQOI.lqs_delcols!(m, col, col)
function LQOI.lqs_delcols!(instance::GLPKOptimizer, col, col2)
    idx = collect(col:col2)
    GLPK.std_basis(instance.inner)
    GLPK.del_cols(instance.inner, length(idx), idx)
end

# LQOI.lqs_addmipstarts(m, colvec, valvec)
LQOI.lqs_addmipstarts!(instance::GLPKOptimizer, colvec, valvec) = nothing
#=
    Solve
=#

# LQOI.lqs_mipopt!(m)
LQOI.lqs_mipopt!(instance::GLPKOptimizer) = opt!(instance)

# LQOI.lqs_qpopt!(m)
LQOI.lqs_qpopt!(instance::GLPKOptimizer) = error("Quadratic solving not supported")

# LQOI.lqs_lpopt!(m)
LQOI.lqs_lpopt!(instance::GLPKOptimizer) = opt!(instance)

const TERMINATION_STATUS_MAP = Dict(
)

# LQOI.lqs_terminationstatus(m)
function LQOI.lqs_terminationstatus(model::GLPKOptimizerMIP)

    if model.userlimit
        return MOI.OtherLimit
    end
    s = GLPK.mip_status(model.inner)
    if s == GLPK.UNDEF
        if model.param.presolve == GLPK.OFF && GLPK.get_status(model.inner) == GLPK.NOFEAS
            return MOI.InfeasibleNoResult
        else
            return MOI.OtherError
        end
    end
    if s == GLPK.OPT
        return MOI.Success
    elseif s == GLPK.INFEAS
        return MOI.InfeasibleNoResult
    elseif s == GLPK.UNBND
        return MOI.UnboundedNoResult
    elseif s == GLPK.FEAS
        return MOI.SlowProgress
    elseif s == GLPK.NOFEAS
        return MOI.OtherError
    elseif s == GLPK.UNDEF
        return MOI.OtherError
    else
        error("internal library error")
    end
end
function LQOI.lqs_terminationstatus(model::GLPKOptimizerLP)
    s = lp_status(model)
    if s == GLPK.OPT
        return MOI.Success
    elseif s == GLPK.INFEAS
        return MOI.InfeasibleNoResult
    elseif s == GLPK.UNBND
        return MOI.UnboundedNoResult
    elseif s == GLPK.FEAS
        return MOI.SlowProgress
    elseif s == GLPK.NOFEAS
        return MOI.InfeasibleOrUnbounded
    elseif s == GLPK.UNDEF
        return MOI.OtherError
    else
        error("Internal library error")
    end
end

function lp_status(lpm::GLPKOptimizerLP)
    if lpm.method == :Simplex || lpm.method == :Exact
        get_status = GLPK.get_status
    elseif lpm.method == :InteriorPoint
        get_status = GLPK.ipt_status
    else
        error("bug")
    end

    s = get_status(lpm.inner)
end

function LQOI.lqs_primalstatus(model::GLPKOptimizerMIP)
    m = model.inner

    s = GLPK.mip_status(model.inner)

    out = MOI.UnknownResultStatus

    if s in [GLPK.OPT]#, GLPK.FEAS]
        out = MOI.FeasiblePoint
    end
    return out
end
function LQOI.lqs_primalstatus(model::GLPKOptimizerLP)
    m = model.inner

    s = lp_status(model)

    out = MOI.UnknownResultStatus

    if s in [GLPK.OPT]#, GLPK.FEAS]
        out = MOI.FeasiblePoint
    end
    return out
end
function LQOI.lqs_dualstatus(model::GLPKOptimizerMIP)
    return MOI.UnknownResultStatus
end
function LQOI.lqs_dualstatus(model::GLPKOptimizerLP)
    m = model.inner

    s = lp_status(model)

    out = MOI.UnknownResultStatus

    if s in [GLPK.OPT]#, GLPK.FEAS]
        out = MOI.FeasiblePoint
    end
    return out
end

# LQOI.lqs_getx!(m, place)
function LQOI.lqs_getx!(instance::GLPKOptimizerMIP, place)
    lp = instance.inner
    for c in eachindex(place)
        place[c] = GLPK.mip_col_val(lp, c)
    end
end
function LQOI.lqs_getx!(lpm::GLPKOptimizerLP, place)
    lp = lpm.inner

    if lpm.method == :Simplex || lpm.method == :Exact
        get_col_prim = GLPK.get_col_prim
    elseif lpm.method == :InteriorPoint
        get_col_prim = GLPK.ipt_col_prim
    else
        error("bug")
    end

    for c in eachindex(place)
        place[c] = get_col_prim(lp, c)
    end
    return nothing
end

# LQOI.lqs_getax!(m, place)
function LQOI.lqs_getax!(instance::GLPKOptimizerMIP, place)
    lp = instance.inner
    for c in eachindex(place)
        place[c] = GLPK.mip_row_val(lp, c)
    end
end
function LQOI.lqs_getax!(lpm::GLPKOptimizerLP, place)
    lp = lpm.inner

    if lpm.method == :Simplex || lpm.method == :Exact
        get_row_prim = GLPK.get_row_prim
    elseif lpm.method == :InteriorPoint
        get_row_prim = GLPK.ipt_row_prim
    else
        error("bug")
    end

    for r in eachindex(place)
        place[r] = get_row_prim(lp, r)
    end
    return nothing
end

# LQOI.lqs_getdj!(m, place)
# no-op
function LQOI.lqs_getdj!(instance::GLPKOptimizerMIP, place) end

function LQOI.lqs_getdj!(lpm::GLPKOptimizerLP, place)
    lp = lpm.inner

    if lpm.method == :Simplex || lpm.method == :Exact
        get_col_dual = GLPK.get_col_dual
    elseif lpm.method == :InteriorPoint
        get_col_dual = GLPK.ipt_col_dual
    else
        error("bug")
    end

    for c in eachindex(place)
        place[c] = get_col_dual(lp, c)
    end
    return nothing
end

# LQOI.lqs_getpi!(m, place)
#no-op
function LQOI.lqs_getpi!(instance::GLPKOptimizerMIP, place) end
function LQOI.lqs_getpi!(lpm::GLPKOptimizerLP, place)
    lp = lpm.inner

    if lpm.method == :Simplex || lpm.method == :Exact
        get_row_dual = GLPK.get_row_dual
    elseif lpm.method == :InteriorPoint
        get_row_dual = GLPK.ipt_row_dual
    else
        error("bug")
    end

    for r in eachindex(place)
        place[r] = get_row_dual(lp, r)
    end
    return nothing
end

# LQOI.lqs_getobjval(m)
LQOI.lqs_getobjval(instance::GLPKOptimizerMIP) = GLPK.mip_obj_val(instance.inner)

function LQOI.lqs_getobjval(lpm::GLPKOptimizerLP)
    if lpm.method == :Simplex || lpm.method == :Exact
        get_obj_val = GLPK.get_obj_val
    elseif lpm.method == :InteriorPoint
        get_obj_val = GLPK.ipt_obj_val
    else
        error("bug")
    end
    return get_obj_val(lpm.inner)
end

# LQOI.lqs_getbestobjval(m)
LQOI.lqs_getbestobjval(instance::GLPKOptimizerMIP) = instance.objbound

# LQOI.lqs_getmiprelgap(m)
LQOI.lqs_getmiprelgap(instance::GLPKOptimizer) = abs(GLPK.mip_obj_val(instance.inner)-instance.objbound)/(1e-9+GLPK.mip_obj_val(instance.inner))

# LQOI.lqs_getitcnt(m)
LQOI.lqs_getitcnt(instance::GLPKOptimizer)  = -1

# LQOI.lqs_getbaritcnt(m)
LQOI.lqs_getbaritcnt(instance::GLPKOptimizer) = -1

# LQOI.lqs_getnodecnt(m)
LQOI.lqs_getnodecnt(instance::GLPKOptimizer) = -1

# LQOI.lqs_dualfarkas(m, place)
LQOI.lqs_dualfarkas!(instance::GLPKOptimizer, place) = getinfeasibilityray(instance, place)

# LQOI.lqs_getray(m, place)
LQOI.lqs_getray!(instance::GLPKOptimizer, place) = getunboundedray(instance, place)

#no-op
function MOI.free!(instance::GLPKOptimizer) end

"""
    writeproblem(m: :MOI.AbstractOptimizer, filename::String)
Writes the current problem data to the given file.
Supported file types are solver-dependent.
"""
writeproblem(instance::GLPKOptimizer, filename::String, flags::String="") = GLPK.write_model(instance.inner, filename)

LQOI.lqs_make_problem_type_continuous(instance::GLPKOptimizer) = GLPK._make_problem_type_continuous(instance.inner)

#=
    old helpers
=#
function opt!(lpm::GLPKOptimizerLP)
    # write_lp(lpm.inner, "model.lp")
    if lpm.method == :Simplex
        solve = GLPK.simplex
    elseif lpm.method == :Exact
        solve = GLPK.exact
    elseif lpm.method == :InteriorPoint
        solve = GLPK.interior
    else
        error("bug")
    end
    return solve(lpm.inner, lpm.param)
end

function opt!(lpm::GLPKOptimizerMIP)
    # write_lp(lpm.inner, "model.lp")
    vartype = getvartype(lpm)
    lb = getvarLB(lpm)
    ub = getvarUB(lpm)
    old_lb = copy(lb)
    old_ub = copy(ub)
    for c in 1:length(vartype)
        vartype[c] in [:Int,:Bin] && (lb[c] = ceil(lb[c]); ub[c] = floor(ub[c]))
        vartype[c] == :Bin && (lb[c] = max(lb[c],0.0); ub[c] = min(ub[c],1.0))
    end
    #lpm.cbdata.vartype = vartype
    try
        setvarLB!(lpm, lb)
        setvarUB!(lpm, ub)
        if lpm.param.presolve == GLPK.OFF
            ret_ps = GLPK.simplex(lpm.inner, lpm.smplxparam)
            ret_ps != 0 && return ret_ps
        end
        ret = GLPK.intopt(lpm.inner, lpm.param)
        if ret == GLPK.EMIPGAP || ret == GLPK.ETMLIM || ret == GLPK.ESTOP
            lpm.userlimit = true
        end
    finally
        setvarLB!(lpm, old_lb)
        setvarUB!(lpm, old_ub)
    end
end

function getvarLB(lpm::GLPKOptimizer)
    lp = lpm.inner
    n = GLPK.get_num_cols(lp)
    lb = Array{Float64}(n)
    for c = 1:n
        l = GLPK.get_col_lb(lp, c)
        if l <= -realmax(Float64)
            l = -Inf
        end
        lb[c] = l
    end
    return lb
end

function setvarLB!(lpm::GLPKOptimizer, collb)
    lp = lpm.inner
    n = GLPK.get_num_cols(lp)
    if nonnull(collb) && length(collb) != n
        error("invalid size of collb")
    end
    for c = 1:n
        u = GLPK.get_col_ub(lp, c)
        if u >= realmax(Float64)
            u = Inf
        end
        if nonnull(collb) && collb[c] != -Inf
            l = collb[c]
            if u < Inf
                if l != u
                    GLPK.set_col_bnds(lp, c, GLPK.DB, l, u)
                else
                    GLPK.set_col_bnds(lp, c, GLPK.FX, l, u)
                end
            else
                GLPK.set_col_bnds(lp, c, GLPK.LO, l, 0.0)
            end
        else
            if u < Inf
                GLPK.set_col_bnds(lp, c, GLPK.UP, 0.0, u)
            else
                GLPK.set_col_bnds(lp, c, GLPK.FR, 0.0, 0.0)
            end
        end
    end
end

function getvarUB(lpm::GLPKOptimizer)
    lp = lpm.inner
    n = GLPK.get_num_cols(lp)
    ub = Array{Float64}(n)
    for c = 1:n
        u = GLPK.get_col_ub(lp, c)
        if u >= realmax(Float64)
            u = Inf
        end
        ub[c] = u
    end
    return ub
end

function setvarUB!(lpm::GLPKOptimizer, colub)
    lp = lpm.inner
    n = GLPK.get_num_cols(lp)
    if nonnull(colub) && length(colub) != n
        error("invalid size of colub")
    end
    for c = 1:n
        l = GLPK.get_col_lb(lp, c)
        if l <= -realmax(Float64)
            l = -Inf
        end
        if nonnull(colub) && colub[c] != Inf
            u = colub[c]
            if l > -Inf
                if l != u
                    GLPK.set_col_bnds(lp, c, GLPK.DB, l, u)
                else
                    GLPK.set_col_bnds(lp, c, GLPK.FX, l, u)
                end
            else
                GLPK.set_col_bnds(lp, c, GLPK.UP, 0.0, u)
            end
        else
            if l > -Inf
                GLPK.set_col_bnds(lp, c, GLPK.LO, l, 0.0)
            else
                GLPK.set_col_bnds(lp, c, GLPK.FR, 0.0, 0.0)
            end
        end
    end
end

const vartype_map = Dict(
    GLPK.CV => :Cont,
    GLPK.IV => :Int,
    GLPK.BV => :Bin
)

function getvartype(lpm::GLPKOptimizer)
    lp = lpm.inner
    ncol = GLPK.get_num_cols(lp)
    coltype = Array{Symbol}(ncol)
    for i in 1:ncol
        ct = GLPK.get_col_kind(lp, i)
        coltype[i] = vartype_map[ct]
        if i in lpm.binaries
            coltype[i] = :Bin
        elseif coltype[i] == :Bin # GLPK said it was binary, but we didn't tell it
            coltype[i] = :Int
        end
    end
    return coltype
end
nonnull(x) = (x != nothing && !isempty(x))

# The functions getinfeasibilityray and getunboundedray are adapted from code
# taken from the LEMON C++ optimization library. This is the copyright notice:
#
### Copyright (C) 2003-2010
### Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
### (Egervary Research Group on Combinatorial Optimization, EGRES).
###
### Permission to use, modify and distribute this software is granted
### provided that this copyright notice appears in all copies. For
### precise terms see the accompanying LICENSE file.
###
### This software is provided "AS IS" with no warranty of any kind,
### express or implied, and with no claim as to its suitability for any
### purpose.

function getinfeasibilityray(lpm::GLPKOptimizerLP, ray)
    lp = lpm.inner

    if lpm.method == :Simplex || lpm.method == :Exact
    elseif lpm.method == :InteriorPoint
        error("getinfeasibilityray is not available when using the InteriorPoint method")
    else
        error("bug")
    end

    m = GLPK.get_num_rows(lp)

    # ray = zeros(m)
    @assert length(ray) == m

    ur = GLPK.get_unbnd_ray(lp)
    if ur != 0
        if ur <= m
            k = ur
            get_stat = GLPK.get_row_stat
            get_bind = GLPK.get_row_bind
            get_prim = GLPK.get_row_prim
            get_ub = GLPK.get_row_ub
        else
            k = ur - m
            get_stat = GLPK.get_col_stat
            get_bind = GLPK.get_col_bind
            get_prim = GLPK.get_col_prim
            get_ub = GLPK.get_col_ub
        end

        get_stat(lp, k) == GLPK.BS || error("unbounded ray is primal (use getunboundedray)")

        ray[get_bind(lp, k)] = (get_prim(lp, k) > get_ub(lp, k)) ? -1 : 1

        GLPK.btran(lp, ray)
    else
        eps = 1e-7
        for i = 1:m
            idx = GLPK.get_bhead(lp, i)
            if idx <= m
                k = idx
                get_prim = GLPK.get_row_prim
                get_ub = GLPK.get_row_ub
                get_lb = GLPK.get_row_lb
            else
                k = idx - m
                get_prim = GLPK.get_col_prim
                get_ub = GLPK.get_col_ub
                get_lb = GLPK.get_col_lb
            end

            res = get_prim(lp, k)
            if res > get_ub(lp, k) + eps
                ray[i] = -1
            elseif res < get_lb(lp, k) - eps
                ray[i] = 1
            else
                continue # ray[i] == 0
            end

            if idx <= m
                ray[i] *= GLPK.get_rii(lp, k)
            else
                ray[i] /= GLPK.get_sjj(lp, k)
            end
        end

        GLPK.btran(lp, ray)

        for i = 1:m
            ray[i] /= GLPK.get_rii(lp, i)
        end
    end

    return nothing
end

function getunboundedray(lpm::GLPKOptimizerLP, ray)
    lp = lpm.inner

    if lpm.method == :Simplex || lpm.method == :Exact
    elseif lpm.method == :InteriorPoint
        error("getunboundedray is not available when using the InteriorPoint method")
    else
        error("bug")
    end

    m = GLPK.get_num_rows(lp)
    n = GLPK.get_num_cols(lp)

    # ray = zeros(n)
    @assert length(ray) == n

    ur = GLPK.get_unbnd_ray(lp)
    if ur != 0
        if ur <= m
            k = ur
            get_stat = GLPK.get_row_stat
            get_dual = GLPK.get_row_dual
        else
            k = ur - m
            get_stat = GLPK.get_col_stat
            get_dual = GLPK.get_col_dual
            ray[k] = 1
        end

        get_stat(lp, k) != GLPK.BS || error("unbounded ray is dual (use getinfeasibilityray)")

        for (ri, rv) in zip(GLPK.eval_tab_col(lp, ur)...)
            ri > m && (ray[ri - m] = rv)
        end

        if (GLPK.get_obj_dir(lp) == GLPK.MAX) $ (get_dual(lp, k) > 0)
            scale!(ray, -1.0)
        end
    else
        for i = 1:n
            ray[i] = GLPK.get_col_prim(lp, i)
        end
    end

    return nothing
end

end # module


text:
# Autogenerated wrapper script for ReadStat_jll for armv7l-linux-gnueabihf
export libreadstat, readstat

using Libiconv_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("ReadStat")
JLLWrappers.@declare_library_product(libreadstat, "libreadstat.so.1")
JLLWrappers.@declare_executable_product(readstat)
function __init__()
    JLLWrappers.@generate_init_header(Libiconv_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libreadstat,
        "lib/libreadstat.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        readstat,
        "bin/readstat",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()