text:
using LinearAlgebra
using GLM
using LsqFit
using Distributions
using StatsBase
using Printf
using Plots

#linear fit wrapper

function lfit(ndf::DataFrame)
    lr = lm(@formula(y_mean ~ x_mean), ndf)
    c = coef(lr)
    return x -> c[1] + c[2]*x, predict(lr), c
end

struct RFit
	fitpar::Vector{Number}
	fitstd::Vector{Number}
	ci::Vector{Tuple{Number, Number}}
	g::Function
end

function gpol1(ct)
    function f1(z)
        return ct[1] + ct[2] * z
    end
    return f1
end

function gpol2(ct)
    function f2(z)
        return ct[1] + ct[2] * z + ct[3] * z^2
    end
    return f2
end

function gpol3(ct)
    function f3(z)
        return ct[1] + ct[2] * z + ct[3] * z^2 + ct[4] * z^3
    end
    return f3
end

function polfit(pol, x::Vector{<:Real},y::Vector{<:Real}, p0::Vector{<:Real})
    fq = curve_fit(pol, x, y, p0)
    cfq = coef(fq)
    @info "coef(fq)" cfq
    sfq = stderror(fq)
    @info "std(fq)" sfq
    @info "margin_of_error (90%)" margin_error(fq, 0.1)
    @info " confidence_interval (90%)" confidence_interval(fq, 0.1)
    return cfq, sfq
end

function fit_pol1(x::Vector{<:Real},y::Vector{<:Real}, ci=0.1)
    @. pol(x, p) = p[1] + p[2] * x
    p0 = [1.0, 1.0]
    fq = curve_fit(pol, x, y, p0)

    RFit(coef(fq), stderror(fq), confidence_interval(fq, ci),
         gpol1(coef(fq)))
end

function fit_pol2(x::Vector{<:Real},y::Vector{<:Real}, ci=0.1)
    @. pol(x, p) = p[1] + p[2] * x + p[3] * x^2
    p0 = [1.0, 1.0, 1.0]
    fq = curve_fit(pol, x, y, p0)
    RFit(coef(fq), stderror(fq), confidence_interval(fq, ci),
         gpol2(coef(fq)))
end

function fit_pol3(x::Vector{<:Real},y::Vector{<:Real}, ci=0.1)
    @. pol(x, p) = p[1] + p[2] * x + p[3] * x^2 + p[4] * x^3
    p0 = [1.0, 1.0, 1.0, 1.0]
    fq = curve_fit(pol, x, y, p0)
    RFit(coef(fq), stderror(fq), confidence_interval(fq, ci),
         gpol3(coef(fq)))
end

"""
	fit_gauss(y, xmin, xmax, bins=25)

Fit a normal distribution to data
"""

struct FGauss
	mu::Vector{Number}
	std::Vector{Number}
	C::Vector{Number}
	h::Histogram
	X::Vector{Number}
	Y::Vector{Number}
	g::Vector{Function}
end

function gausg(μ::Real, σ::Real, C::Real)
	function gausx(x)
		return C * pdf(Normal(μ, σ,), x)
	end
	return gausx
end

function gausg2(μ1::Real, σ1::Real, C1::Real, μ2::Real, σ2::Real, C2::Real)
	function gausx(x)
		return C1 * pdf(Normal(μ1, σ1,), x) + C2 * pdf(Normal(μ2, σ2,), x)
	end
	return gausx
end

@. gauss1fm(x, p) = p[1]* pdf(Normal(p[2], p[3]), x)
@. gauss1(x, p) = p[1]* pdf(Normal(p[2], p[3]), x)
@. gauss2(x, p) = p[1]* pdf(Normal(p[2], p[3]), x) +  p[4]* pdf(Normal(p[5], p[6]), x)
@. gausscm(x, p) = p[1] * pdf(Normal(p[2], p[3]), x) + p[4] * pdf(Normal(p[2], p[5]), x)

function gaussfm(mu::Real)
	function gauss(x::Vector{<:Real}, p::Vector{<:Real})
		return @. p[1]* pdf(Normal(mu, p[2]), x)
	end
	return gauss
end

function gauss2fm(mu::Real)
	function gauss2(x::Vector{<:Real}, p::Vector{<:Real})
		return @. p[1]* pdf(Normal(mu, p[2]), x) + p[3]* pdf(Normal(mu, p[4]), x)
	end
	return gauss2
end

function gauss1fm(mu)
    function gauss1(x,p)
        return @. p[1]* pdf(Normal(mu, p[2]), x)
    end
    return gauss1
end

function cfit(ffit::Function, x::Vector{Float64}, y::Vector{Float64},
	                         p0::Vector{Float64}, lb::Vector{Float64}, ub::Vector{Float64})
	fq = curve_fit(ffit, x, y, p0, lower=lb, upper=ub)
    cfq = coef(fq)
    @debug "coef(fq)" cfq
    return cfq
end

"""
    cfit

Take into account the errors on the measured variables. Assuming weights
of 1/sigma^2 as in standard least squares.
"""
function cfit(ffit::Function, x::Vector{<:Real}, y::Vector{<:Real},
	          yerr::Vector{<:Real}, p0::Vector{<:Real},
              lb::Vector{<:Real}, ub::Vector{<:Real})
	fq = curve_fit(ffit, x, y, yerr.^-2, p0, lower=lb, upper=ub)
    return fq
end

function fit_gauss(h::Histogram)
	c = centers(h)
	w = h.weights * 1.0
	@debug "histo"  w c
	mu, sigma = mean_and_std(c, Weights(h.weights); corrected = false)
	@debug "mu, std" mu, sigma

	# fit parameters lb, ub, po are lower, upper bounds and pars
    lb = [0., mu - 100.0*sigma, sigma/100.0]
    ub = [100*sum(w), mu + 100.0*sigma, 100.0*sigma]
    p0_bounds = [sum(w), mu, sigma]
	CC, μ, σ  = cfit(gauss1, c, w, p0_bounds, lb, ub)
	gx = gausg(μ, σ, CC)
	return FGauss([μ], [σ], [CC], h, c, gx.(c), [gx])

end

fit_gauss(x::Vector{Float64}, xmin::Float64, xmax::Float64;
	      bins::Integer=50, norm=false) =fit_gauss(hist1d(x, bins, xmin, xmax, norm))

"""
	fit_gauss_fm(y::Vector{Float64}, xmin::Float64, xmax::Float64, bins=50, fm=0.0)

Fit a gaussian with a fixed mean fm
"""
function fit_gauss_fm(y::Vector{Float64}, xmin::Float64, xmax::Float64;
	                  bins=50, norm=false, fm=0.0)

	# fit the unbinned distribution
	x  =  in_range(y, xmin, xmax)
	σ = std(x)
	@debug "gfit_gauss_fm: σ = $σ"

	# bin distribution
    h = hist1d(x, bins, xmin, xmax, norm)
    c = centers(h)
    w = h.weights *1.0
    @debug "histo w and c"  w c

	# fit parameters lb, ub, po are lower, upper bounds and pars
    lb = [0.0, σ/10.0]
    ub = [1.0e+6*sum(w),  σ * 10.0]
    p0_bounds = [0.0, σ]

	g1 = gaussfm(fm)
	CC, sigma  = cfit(g1, c, w, p0_bounds, lb, ub)

	@debug "CC,  sigma"  CC  sigma
	gx = gausg(fm, sigma, CC)

	return FGauss([fm], [sigma], [CC], h, c, gx.(c), [gx])
end

function plot_fit_gauss(x::Vector{Float64}, xs::String, ys::String,
    bins::Integer, xmin::Float64, xmax::Float64;
    xgmin::Float64, xgmax::Float64, gbins::Integer=50)

    h,p = hist1d(x, xs, bins, xmin, xmax, norm=true)
    fg  = fit_gauss(x, xgmin, xgmax, bins=gbins, norm=true)
    gx  = fg.g[1]
    X   = centers(h)
    Y   = h.weights
    σY  = sqrt.(Y)
    lbl = @sprintf " μ=%5.1f, σ =%5.1f " fg.mu[1] fg.std[1]
    lbl =string("gaussian fit:\n", lbl)
    p = scatter(X, Y, yerror=σY,fmt = :png,
    shape = :circle, color = :black, label="data", legend=true)
    p = plot!(p, X, gx.(X), lw=2, label=lbl, legend=true, fmt = :png)
    xlabel!(xs)
    xlabel!(ys)
    return fg, p
end

"""
	fitg1(x, xs, xmin, xmax, xgmin, xgmax; bins=100)

returns the fit and the plot
"""
function fitg1(x, xs, bins, xmin, xmax;
	           xgmin, xgmax, fbins=100, norm=true, fm=0.0, flex_mean=false)

	h, p = hist1d(x, xs, bins, xmin, xmax, norm=norm, legend=true)
    if flex_mean
        fg = fit_gauss(x, xgmin, xgmax, bins=fbins, norm=norm)
    else
        fg = fit_gauss_fm(x, xgmin, xgmax, bins=fbins, norm=norm, fm=fm)
    end
	gx = fg.g[1]
	X = centers(h)
	lbl = @sprintf "σ =%4.1f " fg.std[1]
    p = plot!(p, X, gx.(X), lw=2, label=lbl, legend=true, fmt = :png)
	xlabel!(xs)
	return fg, p
end

"""
	gfit_gauss2_cmean(y, xmin, xmax, bins, sigmas, cs, cmean=0.0)

Fit a double gaussian (with sigmas -->[sigma1, sigma2] cs -->[c1, c2] )
and a cmean to data.
"""
function gfit_gauss2_cmean(y::Vector{Float64}, xmin::Float64, xmax::Float64,
	                       bins::Integer, sigmas::Vector{Float64}, cs::Vector{Float64},
						   norm=false, cmean=0.0)

    x =  in_range(y, xmin, xmax)
    h = hist1d(x, bins, xmin, xmax, norm)
    c = centers(h)
    w = h.weights
    @debug "histo centers and weights in full region"  w c

    g2 = gauss2fm(cmean)
    # fit parameters lb, ub, po are lower, upper bounds and pars

    lb = [cs[1]/100.0, sigmas[1]/5.0, cs[2]/100.0, sigmas[2]/5.0]
    ub = [cs[1]*100.0, sigmas[1]*5.0, cs[2]*100.0, sigmas[2]*5.0]
    p0_bounds = [cs[1], sigmas[1], cs[2], sigmas[2]]

    @debug "pars" p0 lb ub
    # fit double gaussian
    fq = curve_fit(g2, c, w, p0_bounds, lower=lb, upper=ub)
    C1, sigma1, C2,  sigma2   = coef(fq)
    @debug "C1 sigma1 C2 sigma2" C1 sigma1 C2  sigma2

    #
    gx = gausg2(cmean, sigma1, C1, cmean, sigma2, C2)
    gx1 = gausg(cmean, sigma1, C1)
    gx2 = gausg(cmean, sigma2, C2)
    #return (sigma1 = sigma1, sigma2 = sigma2, C1 = C1, C2=C2,
    #        h = h, xg = c, yg = gx.(c), gx = gx, gx1=gx1, gx2=gx2)

	return FGauss([cmean, cmean], [sigma1, sigma2], [C1, C2],
			       h, c, gx.(c), [gx, gx1, gx2])
    end

"""
	fit_2gauss_cmean(data, gp, g1p, g2p, cm)

Fit two gaussian with common mean
"""
function fit_2gauss_cmean(data, gp, g1p, g2p, cm, norm=false)
    gf1 = fit_gauss_fm(data, g1p.xmin,g1p.xmax, bins=g1p.nbin, norm=norm, fm=cm)
    @debug gf1
    gf2 = fit_gauss_fm(data, g2p.xmin,g2p.xmax,bins=g2p.nbin, norm= norm, fm=cm)
    @debug gf2
    gf = gfit_gauss2_cmean(data, gp.xmin,gp.xmax,gp.nbin,
	                       [gf1.std[1], gf2.std[1]], [gf1.C[1], gf2.C[1]], norm)
    @debug gf
    return gf
end

# fit_gauss2(x::Vector{Float64},
#            xmin::Vector{Float64},
# 		   xmax::Vector{Float64},
# 		   bins::Vector{Int64}) = gfit_gauss2(x, xmin, xmax, bins)

"""
    fitg2(x, xs, xmin, xmax, xg1min, xg1max, xg2min, xg2max, xgmin, xgmax; bins=100)
Fits 2 gaussians with common mean (0 by default) to vector x.
"""
function fitg2(x, xs, bins, xmin, xmax;
	           xg1min, xg1max, xg2min, xg2max, xgmin, xgmax, cm=0.0,
         	   g1bins=100, g2bins=100, gbins=100, norm=true)

	hp, p = hist1d(x, xs, bins, xmin, xmax, norm=norm, legend=true)

    g1p = (xmin = xg1min, xmax = xg1max, nbin=g1bins)
    g2p = (xmin= xg2min, xmax=  xg2max, nbin=g2bins)
    gp  = (xmin= xgmin, xmax=  xgmax, nbin=gbins)

    fg  = fit_2gauss_cmean(x, gp, g1p, g2p, cm, norm)
	gx  = fg.g[1]
	gx1 = fg.g[2]
	gx2 = fg.g[3]

	lbl = @sprintf "σt =%4.1f mm, σ =%4.1f mm" fg.std[1] fg.std[2]
	st  = @sprintf "σt =%4.1f mm " fg.std[1]
	sf  = @sprintf "σ  =%4.1f mm" fg.std[2]
    p = plot!(p,fg.X, fg.Y, label=lbl, lw=2, fmt = :png)
    p = plot!(p,fg.X, gx1.(fg.X), label=st, lw=1, fmt = :png)
    p = plot!(p,fg.X, gx2.(fg.X), label=sf, lw=1, fmt = :png)
    #xlabel!(xs)
    return fg, p
end

"""
	fit_profile(x1, x2, tx1, ty1, fit="pol1", bins=25)
    fit_profile(df1, c1, c2, tx1, ty1, fit="pol1", bins=25)
Create and fit a profile with pol1 or poli2 functions.
Return fit parameters, fit function and plot
"""
fit_profile(df1::DataFrame, c1::String, c2::String,
            tx1::String, ty1::String, fit="pol1", bins=25;
            ybin_width::Real=0.1, ymin::Real=352.4, ymax::Real=392.4,
            min_proportion::Real=0.0) =
	fit_profile(df1[!,c1], df1[!,c2], tx1, ty1, fit, bins,
                ybin_width=ybin_width, ymin=ymin, ymax=ymax, min_proportion=min_proportion)

function fit_profile(x1::Vector{<:Real}, x2::Vector{<:Real},
	                 tx1::String, ty1::String, fit="pol1", bins=25;
                     ybin_width::Real=0.1, ymin::Real=minimum(x2),
                     ymax::Real=maximum(x2), min_proportion::Real=0.0)

    pdf1, _ = p1df(x1,x2, bins, ybin_width=ybin_width, ymin=ymin, ymax=ymax, min_proportion=min_proportion)

    if fit == "pol1"
        fr = fit_pol1(pdf1.x_mean, pdf1.y_mean)
    elseif fit == "pol2"
        fr = fit_pol2(pdf1.x_mean, pdf1.y_mean)
	elseif fit == "pol3"
		fr = fit_pol3(pdf1.x_mean, pdf1.y_mean)
	else
		println("option not implemented")
		return nothing
    end

    p1 = scatter(pdf1.x_mean,pdf1.y_mean, yerror=pdf1.y_std,fmt = :png,
	          shape = :circle, color = :black, legend=false)
    p1 = plot!(p1, pdf1.x_mean, fr.g.(pdf1.x_mean), fmt = :png)
    xlabel!(tx1)
    ylabel!(ty1)

   return fr, p1
end

max_stars_repo_name: jjgomezcadenas/ATools

max_issues_repo_name: jjgomezcadenas/ATools

max_forks_repo_name: jjgomezcadenas/ATools


text:
using JuMP, FactCheck, AmplNLWriter

## Solve test problem with sind and cosd functions
 #
 #  min   (7 - (3*cosd(x1) + 5*cosd(x2)))^2 + (0 - (3*sind(x1) + 5*sind(x2)))^2
 #  s.t.  x1, x2 continuous
 #
 #  The optimal objective value is 0
 ##

# Allow resolving the model from multiple starts after NLP changes in JuMP 0.12
EnableNLPResolve()

if !isdefined(:solver); solver = CouenneNLSolver(); end

m = Model(solver=solver)
@variable(m, x[1:2])

@NLobjective(m, Min, (7 - (3*cosd(x[1]) + 5*cosd(x[2])))^2 + (0 - (3*sind(x[1]) + 5*sind(x[2])))^2)

context("example: jump_nltrig") do
    setvalue(x[1], 30)
    setvalue(x[2], -50)
    @fact solve(m) --> :Optimal
    @fact getvalue(x)[:] --> roughly([38.21321, -21.78678], 1e-5)
    @fact getobjectivevalue(m) --> roughly(0.0, 1e-5)
    # Now try from the other side
    setvalue(x[1], -30)
    setvalue(x[2], 50)
    @fact solve(m) --> :Optimal
    @fact getvalue(x)[:] --> roughly([-38.21321, 21.78678], 1e-5)
    @fact getobjectivevalue(m) --> roughly(0.0, 1e-5)
end

ext: jl


text:
module test_octo_raw

using Test
using Octo.Adapters.SQL # Raw to_sql

text = """
SELECT a
FROM   b
WHERE  c"""

@test to_sql([Raw(text)]) == text

end # module test_octo_raw


text:
#!/usr/bin/env julia

const DATA_DIR   = "../18_data"
const INPUT_FILE = "input.txt"

const LPAREN = '('
const RPAREN = ')'

const OP_PLUS  = '+'
const OP_TIMES = '*'

# ----------------------------------------
macro assert(bool_expr)
    message = string("Assertion: ", bool_expr, " failed")
    :($(esc(bool_expr)) || error($message))
end
# ----------------------------------------
@enum Op op_plus op_times

mutable struct PartialExpr
    is_subexpr::Bool
    lval::Union{Int, Nothing}
    op::Union{Op, Nothing}
    # --------------------
    function PartialExpr(is_subexpr=false)
        return new(is_subexpr, nothing, nothing)
    end
end
# ----------------------------------------

function pexpr_addval!(pexpr, val)
    if isnothing(pexpr.lval)
        pexpr.lval = val
    else
        if isnothing(pexpr.op)
            error("Encountered digit where operator expected")
        else
            new_lval = eval_op(pexpr.lval, pexpr.op, val)
            pexpr.lval = new_lval
            pexpr.op = nothing
        end
    end
end

function pexpr_print(pexpr::PartialExpr)
    op_str = isnothing(pexpr.op) ? "" : (pexpr.op == op_plus ? "+" : "*")
    subexpr_str = pexpr.is_subexpr ? "(subexpr)" : ""
    output = "PExpr: $(pexpr.lval) $(op_str) $(subexpr_str)"
    println(strip(output))
end

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

function eval_op(arg_left::Int, op::Op, arg_right::Int)
    if op == op_plus
        return arg_left + arg_right
    elseif op == op_times
        return arg_left * arg_right
    elseif isnothing(op)
        error("operator is not set")
    else
        error("Unrecognized operator value")
    end
end

function eval_expr(expr::String)
    partials = Array{PartialExpr,1}()
    push!(partials, PartialExpr())
    for c in expr
        if isspace(c)
            continue
        end
        top = last(partials)
        if isdigit(c)
            val = parse(Int, c)
            pexpr_addval!(top, val)
        elseif c == LPAREN
            push!(partials, PartialExpr(true))
        elseif c == RPAREN
            @assert top.is_subexpr
            @assert !isnothing(top.lval)
            @assert isnothing(top.op)
            val = top.lval
            _ = pop!(partials)
            top = last(partials)
            pexpr_addval!(top, val)
        elseif c == OP_PLUS
            @assert !isnothing(top.lval)
            @assert isnothing(top.op)
            top.op = op_plus
        elseif c == OP_TIMES
            @assert !isnothing(top.lval)
            @assert isnothing(top.op)
            top.op = op_times
        else
            error("Unrecognized character: $(c)")
        end
    end
    @assert length(partials) == 1
    only = partials[1]
    @assert !isnothing(only.lval)
    @assert isnothing(only.op)
    return only.lval
end

function main(filename)
    lines::Array{String} = readlines(filename)
    result = 0
    for line in lines
        # println("line=$(line)")
        val = eval_expr(line)
        # println("\tvalue=$(val)")
        result += val
    end
    println("Total of all exprs=$(result)")
end

main("$DATA_DIR/$INPUT_FILE")

max_stars_repo_name: jaycoskey/AdventOfCode

max_issues_repo_name: jaycoskey/AdventOfCode

max_forks_repo_name: jaycoskey/AdventOfCode


text:
using VulkanGen

const destdir = joinpath(dirname(dirname(@__DIR__)), "generated")

wrapper_config(platform, destfile) = WrapperConfig(platform, joinpath(destdir, destfile))

configs = [
    wrapper_config(Linux(), "linux.jl"),
    wrapper_config(MacOS(), "macos.jl"),
    wrapper_config(BSD(), "bsd.jl"),
    wrapper_config(Windows(), "windows.jl"),
]

for config in configs
    @time vw = VulkanWrapper(config)
    @time write(vw, config)
    @info "Vulkan successfully wrapped at $(config.destfile)"
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
function movetostring(move::Move)
    # extract move details
    from_sqr = from(move)
    to_sqr = to(move)
    flags = flag(move)
    # convert to string
    ucistring = LABELLED_SQUARES[from_sqr]
    ucistring *= LABELLED_SQUARES[to_sqr]
    if flags == __KNIGHT_PROMO
        ucistring *= "n"
    elseif flags == __BISHOP_PROMO
        ucistring *= "b"
    elseif flags == __ROOK_PROMO
        ucistring *= "r"
    elseif flags == __QUEEN_PROMO
        ucistring *= "q"
    end
    return ucistring
end

function uci_main()
    if iszero(length(ARGS))
        # by default, load starting position
        board = importfen(START_FEN)
        threads = ThreadPool(1)
    else
        # else, a FEN was supplied in command line (assume)
        board = importfen(ARGS[1])
        println(io, ARGS[2])
        threads = ThreadPool(Threads.nthreads())
    end
    # init thread boards
    setthreadpoolboard!(threads, board)
    ttable = TT_Table()

    # main UCI loop
    while true
        line = readline()

        if line == "uci"
            uci_engine()
            continue

        elseif line == "isready"
            uci_isready()
            continue

        elseif line == "ucinewgame"
            uci_newgame!(threads, ttable)
        elseif line == "quit"
            break
        end

        splitlines = split(line)

        if splitlines[1] == "go"
            uci_go(threads, ttable, splitlines)
            continue

        elseif splitlines[1] == "position"
            uci_position!(threads, splitlines)

        elseif splitlines[1] == "perft"
            uci_perft(threads, splitlines)

        elseif splitlines[1] == "setoption"
            uci_setoptions(threads, splitlines, ttable)
        end
        # additional options currently unsupported
    end
    return
end

function uci_engine()
    print("id name Astellarn ", ASTELLARN_VERSION, "\nid author Jeremy Revell\nuciok\n")
end

function uci_isready()
    print("readyok\n")
end

function uci_newgame!(threads::ThreadPool, ttable::TT_Table)
    board = importfen(START_FEN)
    threads[1] = Thread()
    setthreadpoolboard!(threads, board)
    ttable.table = Dict{UInt64, TT_Entry}()
    return
end

function uci_perft(threads::ThreadPool, splitlines::Vector{SubString{String}})
    depth = parse(Int, splitlines[2])
    start_time = time()
    nodes = perft(threads[1].board, depth)
    stop_time = time()
    elapsed = stop_time - start_time
    nps = nodes/elapsed
    @printf("Total time (ms) : %d\n", elapsed*1000)
    @printf("Nodes searched : %d\n", nodes)
    @printf("Nodes/second :  %d\n", nps)
    return
end

function uci_go(threads::ThreadPool, ttable::TT_Table, splitlines::Vector{SubString{String}})
    # First, take note of the starting time.
    start_time = time()

    # Initialise the default search + time management values.
    infinite = false
    depth = wtime = btime = movetime = winc = binc = 0
    movestogo = -1

    # Extract search and time management parameters from the UCI command.
    for i in eachindex(splitlines)
        if splitlines[i] == "infinite"
            infinite = true
            break # We can break as we expect no more inputs
        end
        # The explicit depth.
        if splitlines[i] == "depth"
            depth = parse(Int, splitlines[i + 1])
            infinite = true # So that we know we are not time limited
            break # We can assume no more inputs?
        end
        # White's clock time.
        if splitlines[i] == "wtime"
            wtime = parse(Int, splitlines[i + 1])
            continue
        end
        # Black's clock time.
        if splitlines[i] == "btime"
            btime = parse(Int, splitlines[i + 1])
            continue
        end
        # White increment amount.
        if splitlines[i] == "winc"
            winc = parse(Int, splitlines[i + 1])
            continue
        end
        # Black increment amount.
        if splitlines[i] == "binc"
            binc = parse(Int, splitlines[i + 1])
            continue
        end
        # The movetime we have available.
        if splitlines[i] == "movetime"
            movetime = parse(Int, splitlines[i + 1])
            continue
        end
        # Find out the movestogo.
        if splitlines[i] == "movestogo"
            movestogo = parse(Int, splitlines[i + 1])
            continue
        end
    end

    # Which colour's turn is it?
    if threads[1].board.turn == WHITE
        clock_time = wtime
        inc_time = winc
    else
        clock_time = btime
        inc_time = binc
    end

    # Init the time manager and set to the thread
    threads[1].timeman.isinfinite = infinite
    threads[1].timeman.movetime   = movetime
    threads[1].timeman.depth      = depth
    threads[1].timeman.start_time = start_time
    threads[1].timeman.clock_time = clock_time
    threads[1].timeman.inc_time   = inc_time
    threads[1].timeman.movestogo  = movestogo
    threads[1].timeman.ideal_time = 0
    threads[1].timeman.max_time   = 0

    threads[1].stop = false

    # following section is a hack, while multithreadding is not supported.
    threads[1].ss.nodes = 0
    threads[1].ss.depth = 0
    threads[1].ss.seldepth = 0
    threads[1].ss.tbhits = 0

    # Calling the main search function.
    if depth == 0
        depth = 50
    end
    eval = find_best_move(threads[1], ttable, depth)

    # Extract the search info and stats.
    move = threads[1].pv[1][1]
    nodes = threads[1].ss.nodes
    elapsed = elapsedtime(threads[1].timeman)
    nps = nodes*1000/elapsed
    hf = hashfull(ttable)

    # Report back to the UCI
    ucistring = movetostring(move)
    if !threads[1].stop
        @printf("info depth %d seldepth %d nodes %d nps %d tbhits %d score cp %d hashfull %d pv ", threads[1].ss.depth, threads[1].ss.seldepth, nodes, nps, threads[1].ss.tbhits, eval, hf)
        print(join(movetostring.(threads[1].pv[1]), " "))
        print("\n")
    end
    print("bestmove ", ucistring, "\n")
    return
end

function uci_position!(threads::ThreadPool, splitlines::Vector{SubString{String}})
    # import a given FEN string
    if splitlines[2] == "startpos"
        board = importfen(START_FEN)
    elseif splitlines[2] == "fen"
        fen = join(splitlines[3:8], " ")
        board = importfen(fen)
    end
    # make the given moves, if any
    makemoves = false
    movestack = MoveStack(MAX_MOVES)
    for i in eachindex(splitlines)
        if splitlines[i] == "moves"
            makemoves = true
            continue
        end
        if makemoves == true
            gen_moves!(movestack, board)
            for move in movestack
                if movetostring(move) == splitlines[i]
                    apply_move!(board, move)
                    break
                end
            end
            clear!(movestack)
        end
    end
    setthreadpoolboard!(threads, board)
    return
end

function uci_setoptions(threads::ThreadPool, splitlines::Vector{SubString{String}}, ttable::TT_Table)
    if splitlines[3] == "Hash"
        newtable = TT_Table(parse(Int, splitlines[5]))
        ttable.table = newtable.table
        ttable.hashmask = newtable.hashmask
        println("info string set Hash to ", splitlines[5])
    end

    if splitlines[3] == "Threads"
        num_threads = parse(Int, splitlines[5])
        createthreadpool(threads[1].board, num_threads)
        println("info string set Threads to ", num_threads)
    end

    if splitlines[3] == "SyzygyPath"
        tb_init(splitlines[5])
        println("info string set SyzygyPath to ", splitlines[5])
    end

    if splitlines[3] == "MoveOverhead"
        threads[1].timeman.move_overhead = parse(Int, splitlines[5])
        println("info string set MoveOverhead to ", splitlines[5])
    end
    return
end

# This function is called when the thread (or search) wishes to output stats mid-execution.
function uci_report(thread::Thread, ttable::TT_Table, α::Int, β::Int, value::Int)
    score = max(α, min(β, value))
    elapsed = elapsedtime(thread.timeman)
    nps = thread.ss.nodes*1000 / elapsed
    hf = hashfull(ttable)
    @printf("info depth %d seldepth %d nodes %d nps %d tbhits %d score cp %d hashfull %d pv ", thread.ss.depth, thread.ss.seldepth, thread.ss.nodes, nps, thread.ss.tbhits, score, hf)
    print(join(movetostring.(thread.pv[1]), " "))
    print("\n")
    return
end

max_stars_repo_name: J-Revell/Astellarn

max_issues_repo_name: J-Revell/Astellarn

max_forks_repo_name: J-Revell/Astellarn


text:
using SymbolicRegression
using Test
using Random

begin
    X = 2 .* randn(MersenneTwister(0), Float32, 2, 1000)
    y = 3 * cos.(X[2, :]) + X[1, :] .^ 2 .- 2

    options = SymbolicRegression.Options(;
        binary_operators=(+, *, /, -),
        unary_operators=(cos,),
        crossoverProbability=0.0,  # required for recording, as not set up to track crossovers.
        max_evals=10000,
        deterministic=true,
        seed=0,
        verbosity=0,
        progress=false,
    )

    all_outputs = []
    for i in 1:2
        local hall_of_fame = EquationSearch(
            X, y; niterations=5, options=options, multithreading=false, numprocs=0
        )
        local dominating = calculate_pareto_frontier(X, y, hall_of_fame, options)
        push!(all_outputs, dominating[end].tree)
    end

    @test string(all_outputs[1]) == string(all_outputs[2])
end


text:
using Test
using NQCDynamics
using ComponentArrays

atoms = Atoms([:C])
model = NQCModels.Harmonic()
cell = PeriodicCell(hcat(1))
sim = Simulation(atoms, model; cell=cell)

z = ComponentVector(v=fill(1.0, 1, length(atoms)), r=zeros(1, length(atoms)))

solution = run_trajectory(z, (0.0, 300), sim; dt=0.1)
Rs = vcat(get_positions.(solution.u)...)
@test !isapprox(minimum(Rs), 0, atol=1e-4) | (minimum(Rs) >= 0) # check atom leaves cell

solution = run_trajectory(z, (0.0, 300), sim; callback=DynamicsUtils.CellBoundaryCallback(), dt=0.1)
Rs = vcat(get_positions.(solution.u)...)
@test isapprox(minimum(Rs), 0, atol=1e-4) | (minimum(Rs) >= 0) # Check atom remains inside


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule qwt_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("qwt")
JLLWrappers.@generate_main_file("qwt", UUID("ed0789fa-10db-50b3-94da-03266d70be0f"))
end  # module qwt_jll


text:
#
# Univariate Calculus
#

macro test_simplified_isequal(lhs, rhs)
    quote
        @test isequal(simplify($(esc(lhs))), $(esc(rhs)))
    end
end

@test isequal(differentiate(:(2), :x), 0)
@test isequal(differentiate(:(x), :x), 1)
@test_simplified_isequal differentiate(:(x + x), :x) 2
@test_simplified_isequal differentiate(:(x - x), :x) 0
@test_simplified_isequal differentiate(:(2 * x), :x) 2
@test_simplified_isequal differentiate(:(2 / x), :x) :(-2 / x^2)
@test_simplified_isequal differentiate(:(x / 2), :x) 0.5
@test_simplified_isequal differentiate(:(sin(x) / x), :x) :((cos(x) * x - sin(x)) / x^2)
@test_simplified_isequal differentiate(:(x * 2), :x) 2
@test_simplified_isequal differentiate(:(a * x), :x) :a
@test_simplified_isequal differentiate(:(x * a), :x) :a
@test_simplified_isequal differentiate(:(x ^ 2), :x) :(2 * x)
@test_simplified_isequal differentiate(:(a * x ^ 2), :x) :(a * (2 * x))
@test_simplified_isequal differentiate(:(2 ^ x), :x) :(*(0.6931471805599453, ^(2, x)))
@test_simplified_isequal differentiate(:(sin(x)), :x) :(cos(x))
@test_simplified_isequal differentiate(:(cos(x)), :x) :(-sin(x))
@test_simplified_isequal differentiate(:(tan(x)), :x) :(1 + tan(x)^2)
@test_simplified_isequal differentiate(:(exp(x)), :x) :(exp(x))
@test_simplified_isequal differentiate(:(log(x)), :x) :(1 / x)
@test_simplified_isequal differentiate(:(sin(x) + sin(x)), :x) :(cos(x) + cos(x))
@test_simplified_isequal differentiate(:(sin(x) - cos(x)), :x) :(cos(x) + sin(x))
@test_simplified_isequal differentiate(:(x * sin(x)), :x) :(sin(x) + x * cos(x))
@test_simplified_isequal differentiate(:(x / sin(x)), :x) :((sin(x) - x * cos(x)) / (sin(x)^2))
@test_simplified_isequal differentiate(:(sin(sin(x))), :x) :(*(cos(x),cos(sin(x))))
@test_simplified_isequal differentiate(:(sin(cos(x) + sin(x))), :x) :(*(+(-sin(x),cos(x)),cos(+(cos(x),sin(x)))))
@test_simplified_isequal differentiate(:(exp(-x)), :x) :(-exp(-x))
@test_simplified_isequal differentiate(:(log(x^2)), :x) :((2x) * (1 / x^2))
@test_simplified_isequal differentiate(:(abs2(x)), :x) :(2x)
@test_simplified_isequal differentiate(:(inv(x)), :x) :(-abs2(inv(x)))
@test_simplified_isequal differentiate(:(x^n), :x) :(*(n, ^(x, -(n, 1))))
@test_simplified_isequal differentiate(:(n^x), :x) :(*(^(n, x), log(n)))
@test_simplified_isequal differentiate(:(n^n), :x) 0

#
# Multivariate Calculus
#

@test isequal(simplify.(differentiate(:(sin(x) + sin(y)), [:x, :y])), [:(cos(x)), :(cos(y))])
@test isequal(simplify.(differentiate(:(x^2), [:x, :y])), Any[:(2*x), 0])

# TODO: Get the generalized power rule right.
# @test isequal(differentiate(:(sin(x)^2), :x), :(2 * sin(x) * cos(x)))

#
# Strings instead of symbols
#

# @test isequal(differentiate("sin(x) + cos(x)^2"), :(+(cos(x),*(2,cos(x)))))
@test_simplified_isequal differentiate("x + exp(-x) + sin(exp(x))", :x) :(+(1,-exp(-x),*(exp(x),cos(exp(x)))))

# TODO: Make these work
# differentiate(:(sin(x)), :x)(0.0)
# differentiate(:(sin(x)), :x)(1.0)
# differentiate(:(sin(x)), :x)(pi)

#
# SymbolicVariable use
#

x = BasicVariable(:x)
y = BasicVariable(:y)

@test isequal(@sexpr(x + y), :($x + $y))
@test_simplified_isequal differentiate(@sexpr(3 * x), x) 3
@test_simplified_isequal differentiate(:(sin(sin(x))), :x) :(*(cos(x),cos(sin(x))))
@test_simplified_isequal differentiate(@sexpr(sin(sin(x))), x) :(*(cos($x),cos(sin($x))))

function testfun(x)
    z = BasicVariable(:z)
    simplify(differentiate(@sexpr(3*x + x^2*z), z))
end

@test isequal(testfun(x), :(^($(x),2)))
@test isequal(testfun(3), 9)
@test isequal(testfun(@sexpr(x+y)), :(^(+($x,$y),2)))

#
# Simplify tests
#

@test isequal(simplify(:(x-0)), :x)
@test isequal(simplify(:(0-x)), :(-x))

@test isequal(simplify(:(x+y)), :(+(x,y)))
@test isequal(simplify(:(x+3)), :(+(3,x)))
@test isequal(simplify(:(x+3+4)), :(+(7,x)))
@test isequal(simplify(:(2+y+x+3)), :(+(5,y,x)))

@test isequal(simplify(:(x*y)), :(*(x,y)))
@test isequal(simplify(:(x*3)), :(*(3,x)))
@test isequal(simplify(:(x*3*4)), :(*(12,x)))
@test isequal(simplify(:(2*y*x*3)), :(*(6,y,x)))

@test isequal(simplify(:(sin((1*(sin(0/0)))))), NaN)

#
# Tests with ifelse
#
@test_simplified_isequal differentiate(:(ifelse(x < 1, exp(x^2), 1/x)), :x) :(ifelse(x < 1,2x * exp(x^2), -1/x^2))
@test_simplified_isequal differentiate(:(ifelse(x <= 0, 0, ifelse(x > 1, 1, x))), :x) :(ifelse(x <= 0, 0, ifelse(x > 1, 0, 1)))


text:
using Newtonsmethod
using Test
using LinearAlgebra

# Testing the first function (unit test)
f(x) = x^2-16.0
f′(x) = 2.0*x
@test newtonroot(f,f′;x₀=10.0)[1] ≈ 4.0 atol=0.00001

f(x) = log(x)-1.0
f′(x) = 1/x
@test newtonroot(f,f′;x₀=3.0)[1] ≈ 2.71828182 atol=0.00001

f(x) = exp(x) - 2.5
f′(x) = exp(x)
@test newtonroot(f,f′;x₀=10.0)[1] ≈ 0.9162907318741551 atol=0.00001

# Testing the second function
f(x) = x^2-16.0
@test newtonroot(f;x₀=10.0)[1] ≈ 4.0 atol=0.00001

f(x) = log(x)-1.0
@test newtonroot(f;x₀=3.0)[1] ≈ 2.71828182 atol=0.00001

f(x) = exp(x) - 2.5
@test newtonroot(f;x₀=10.0)[1] ≈ 0.9162907318741551 atol=0.00001

# Testing the BigFloat
f(x) = exp(x) - 2.5
@test newtonroot(f;x₀ = BigFloat(10.0), tolerance = BigFloat(1E-10), maxiter = BigFloat(1000))[1] ≈ 0.9162907318741551 atol=0.00001

# Testing non-convergence
f(x) = x^2+2.0
f′(x) = 2.0*x
@test newtonroot(f,f′;x₀=1.0) == nothing

# Testing maxiter
f(x) = exp(x) - 2.5
a=newtonroot(f;x₀=10.0)[1]
b=newtonroot(f;x₀=10.0,maxiter=5)
@testset "maxiter" begin
    @test a ≈ 0.9162907318741551 atol=0.00001
    @test b == nothing
end;

# Testing tolerance level
g(x) = exp(x) - 2.5
x₁ = newtonroot(g;x₀=10.0)[1]
x₂ = newtonroot(g;x₀=10.0,tolerance=1E-3)[1]
x₃ = newtonroot(g;x₀=10.0,tolerance=1)[1]
@testset "tolerance" begin
    @test norm(g(x₁)) < norm(g(x₂))
    @test norm(g(x₁)) < norm(g(x₃))
    @test norm(g(x₂)) < norm(g(x₃))
end;

# Give a fail test
# g(x) = exp(x) - 2.5
# @test newtonroot(g;x₀=1.0)[1] == 1.0


text:
function breaksymmetries(tp::Union{NodeTreeProblem,TreeProblem};
    rules::Vector{Symbol} = [:leftfirst])
    bt = tp.bt
    pd = tp.pd
    levels = 1:tp.nlevels
    function addleftfirstrule(cb)
        xval = JuMP.getvalue(tp.assign)
        for u in getnodes(bt, bt.depth-1)
            children = getchildren(bt, u)
            left = children[1]; right = children[2]
            if round(sum(xval[:,left,levels])) < round(sum(xval[:,right,levels]))
                JuMP.@lazyconstraint(cb,
                    sum(tp.assign[:,right,levels]) <= sum(tp.assign[:,left,levels]))
            end
        end
    end
    function addalphabetizerule(cb)
        xval = JuMP.getvalue(tp.assign)
        for u in getnodes(bt, bt.depth-1)
            children = getchildren(bt, u)
            in(tp.outgroupnode, children) && continue
            left = children[1]; right = children[2]
            if round(sum(xval[:,left,1])) + round(sum(xval[:,right,1])) == 2
                a = findfirst(xval[:,left,1] .> 0.1)
                b = findfirst(xval[:,right,1] .> 0.1)
                if a == nothing
                    println(xval[:,left,1])
                    println(xval[:,right,1])
                end
                if b == nothing
                    println(xval[:,right,1])
                    println(xval[:,left,1])
                end
                a <= b && continue
                JuMP.@lazyconstraint(cb,
                    sum(tp.assign[1:b,right,1]) + sum(tp.assign[(b+1):pd.npop,left,1]) <= 1)
            end
        end
    end
    if in(:leftfirst, rules)
        JuMP.addlazycallback(tp.model, addleftfirstrule)
    end
    if in(:alphabetize, rules)
        JuMP.addlazycallback(tp.model, addalphabetizerule)
    end
end


text:
using BinaryBuilder

name = "CeresSolver"
version = v"1.14.0"

sources = [
    ArchiveSource("http://ceres-solver.org/ceres-solver-$(version).tar.gz",
                  "4744005fc3b902fed886ea418df70690caa8e2ff6b5a90f3dd88a3d291ef8e8e")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/ceres-solver-1.14.0

mkdir cmake_build
cd cmake_build/

if [[ "$target" == *-freebsd* || "$target" == *-apple-* ]]; then
  # Clang doesn't play nicely with OpenMP and
  # compilation fails with glog due to a c++11 error
  CMAKE_FLAGS=(-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN%.*}_gcc.cmake
               -DMINIGLOG=ON)
else
  CMAKE_FLAGS=(-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN}
               -DGLOG_INCLUDE_DIR_HINTS=${prefix}/include)
fi

CMAKE_FLAGS+=(-DCMAKE_INSTALL_PREFIX=${prefix}
              -DBUILD_SHARED_LIBS=ON
              -DBUILD_EXAMPLES=OFF
              -DBUILD_TESTING=OFF
              -DOPENMP=ON
              -DTBB=OFF
              -DBLAS_LIBRARIES=${libdir}/libopenblas.${dlext}
              -DLAPACK_LIBRARIES=${libdir}/libopenblas.${dlext}
              -DMETIS_LIBRARY=${libdir}/libmetis.${dlext}
              -DSUITESPARSE_INCLUDE_DIR_HINTS=${prefix}/include
              -DAMD_LIBRARY="${libdir}/libamd.${dlext} ${libdir}/libsuitesparseconfig.${dlext}"
              -DCAMD_LIBRARY="${libdir}/libcamd.${dlext} ${libdir}/libsuitesparseconfig.${dlext}"
              -DCCOLAMD_LIBRARY="${libdir}/libccolamd.${dlext} ${libdir}/libsuitesparseconfig.${dlext}"
              -DCHOLMOD_LIBRARY="${libdir}/libcholmod.${dlext} ${libdir}/libsuitesparseconfig.${dlext}"
              -DCOLAMD_LIBRARY="${libdir}/libcolamd.${dlext} ${libdir}/libsuitepsarseconfig.${dlext}"
              -DSUITESPARSEQR_LIBRARY="${libdir}/libspqr.${dlext} ${libdir}/libsuitesparseconfig.${dlext}"
              -DSUITESPARSE_CONFIG_LIBRARY="${libdir}/libsuitesparseconfig.${dlext}"
              )

cmake ${CMAKE_FLAGS[@]} ..

make -j${nproc}
make install

cd ..
"""

platforms = expand_cxxstring_abis(supported_platforms())

# The products that we will ensure are always built
products = Product[
    LibraryProduct("libceres", :libceres),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Eigen_jll"),
    Dependency("glog_jll"),
    Dependency("METIS_jll"),
    Dependency("OpenBLAS32_jll"),
    Dependency("SuiteSparse_jll"),
    Dependency("CompilerSupportLibraries_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module TestPrint

using Test

using VecMathOptInterface
const MOI = VecMathOptInterface
const MOIU = MOI.Utilities

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

const LATEX = MIME("text/latex")
const PLAIN = MIME("text/plain")
const IN = Sys.iswindows() ? "in" : "∈"

# Windows fun...
function _string_compare(a, b)
    @test replace(a, "\r\n" => "\n") == replace(b, "\r\n" => "\n")
    return
end

function test_nonname_variable()
    model = MOIU.Model{Float64}()
    x = MOI.add_variable(model)
    @test MOIU._to_string(PLAIN, model, x) == "v[1]"
    @test MOIU._to_string(LATEX, model, x) == "v_{1}"
end

function test_numbers()
    options =
        MOIU._PrintOptions(MIME("text/plain"); simplify_coefficients = true)
    @test MOIU._to_string(options, 1.0, "x"; is_first = true) == "x"
    @test MOIU._to_string(options, 1.0, "x"; is_first = false) == " + x"
    @test MOIU._to_string(options, -1.0, "x"; is_first = true) == "-x"
    @test MOIU._to_string(options, -1.0, "x"; is_first = false) == " - x"
    @test MOIU._to_string(options, 1.2, "x"; is_first = true) == "1.2 x"
    @test MOIU._to_string(options, 1.2, "x"; is_first = false) == " + 1.2 x"
    @test MOIU._to_string(options, -1.2, "x"; is_first = true) == "-1.2 x"
    @test MOIU._to_string(options, -1.2, "x"; is_first = false) == " - 1.2 x"

    options =
        MOIU._PrintOptions(MIME("text/plain"); simplify_coefficients = false)
    @test MOIU._to_string(options, 1.0, "x"; is_first = true) == "1.0 x"
    @test MOIU._to_string(options, 1.0, "x"; is_first = false) == " + 1.0 x"
    @test MOIU._to_string(options, -1.0, "x"; is_first = true) == "-1.0 x"
    @test MOIU._to_string(options, -1.0, "x"; is_first = false) == " - 1.0 x"
    @test MOIU._to_string(options, 1.2, "x"; is_first = true) == "1.2 x"
    @test MOIU._to_string(options, 1.2, "x"; is_first = false) == " + 1.2 x"
    @test MOIU._to_string(options, -1.2, "x"; is_first = true) == "-1.2 x"
    @test MOIU._to_string(options, -1.2, "x"; is_first = false) == " - 1.2 x"
end

function test_variable()
    model = MOIU.Model{Float64}()
    x = MOI.add_variable(model)
    for (name, latex_name) in
        [("x", "x"), ("x_y", "x\\_y"), ("x^2", "x\\^2"), ("x[a,b]", "x_{a,b}")]
        MOI.set(model, MOI.VariableName(), x, name)
        @test MOIU._to_string(PLAIN, model, x) == name
        @test MOIU._to_string(LATEX, model, x) == latex_name
    end
end

function test_ScalarAffineTerm()
    model = MOIU.Model{Float64}()
    x = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, "x")
    @test MOIU._to_string(
        PLAIN,
        model,
        MOI.ScalarAffineTerm(-1.2, x);
        is_first = false,
    ) == " - 1.2 x"
    @test MOIU._to_string(
        LATEX,
        model,
        MOI.ScalarAffineTerm(1.2, x);
        is_first = false,
    ) == " + 1.2 x"
end

function test_ScalarAffineFunction()
    model = MOIU.Model{Float64}()
    x = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, "x")
    f = MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.([-1.2, 1.3], [x, x]),
        1.4,
    )
    @test MOIU._to_string(PLAIN, model, f) == "1.4 - 1.2 x + 1.3 x"
    @test MOIU._to_string(LATEX, model, f) == "1.4 - 1.2 x + 1.3 x"
end

function test_ScalarQuadraticTerm()
    model = MOIU.Model{Float64}()
    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, "x")
    MOI.set(model, MOI.VariableName(), y, "y")
    term = MOI.ScalarQuadraticTerm(-1.2, x, x)
    @test MOIU._to_string(PLAIN, model, term; is_first = false) == " - 0.6 x²"
    @test MOIU._to_string(LATEX, model, term; is_first = false) == " - 0.6 x^2"
    term = MOI.ScalarQuadraticTerm(1.2, x, y)
    @test MOIU._to_string(PLAIN, model, term; is_first = false) == " + 1.2 x*y"
    @test MOIU._to_string(LATEX, model, term; is_first = false) ==
          " + 1.2 x\\times y"
end

function test_ScalarQuadraticFunction()
    model = MOIU.Model{Float64}()
    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, "x")
    MOI.set(model, MOI.VariableName(), y, "y")
    f = MOI.ScalarQuadraticFunction(
        MOI.ScalarQuadraticTerm.([0.5, 0.6], [x, x], [x, y]),
        MOI.ScalarAffineTerm.([-1.2, 1.3], [x, x]),
        1.4,
    )
    @test MOIU._to_string(PLAIN, model, f) ==
          "1.4 - 1.2 x + 1.3 x + 0.25 x² + 0.6 x*y"
    @test MOIU._to_string(LATEX, model, f) ==
          "1.4 - 1.2 x + 1.3 x + 0.25 x^2 + 0.6 x\\times y"
end

function test_VectorOfVariables()
    model = MOIU.Model{Float64}()
    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, "x")
    MOI.set(model, MOI.VariableName(), y, "y")
    f = MOI.VectorOfVariables([x, y])
    @test MOIU._to_string(PLAIN, model, f) == "┌ ┐\n│x│\n│y│\n└ ┘"
    @test MOIU._to_string(LATEX, model, f) ==
          "\\begin{bmatrix}\nx\\\\\ny\\end{bmatrix}"
end

function test_LessThan()
    s = MOI.LessThan(1.2)
    @test MOIU._to_string(PLAIN, s) == "<= 1.2"
    @test MOIU._to_string(LATEX, s) == "\\le 1.2"
end

function test_GreaterThan()
    s = MOI.GreaterThan(1.2)
    @test MOIU._to_string(PLAIN, s) == ">= 1.2"
    @test MOIU._to_string(LATEX, s) == "\\ge 1.2"
end

function test_EqualTo()
    s = MOI.EqualTo(1.2)
    @test MOIU._to_string(PLAIN, s) == "== 1.2"
    @test MOIU._to_string(LATEX, s) == "= 1.2"
end

function test_Interval()
    s = MOI.Interval(1.2, 1.3)
    @test MOIU._to_string(PLAIN, s) == "$(IN) [1.2, 1.3]"
    @test MOIU._to_string(LATEX, s) == "\\in \\[1.2, 1.3\\]"
end

function test_ZeroOne()
    s = MOI.ZeroOne()
    @test MOIU._to_string(PLAIN, s) == "$(IN) {0, 1}"
    @test MOIU._to_string(LATEX, s) == "\\in \\{0, 1\\}"
end

function test_Integer()
    s = MOI.Integer()
    @test MOIU._to_string(PLAIN, s) == "$(IN) ℤ"
    @test MOIU._to_string(LATEX, s) == "\\in \\mathbb{Z}"
end

function test_ExponentialCone()
    s = MOI.ExponentialCone()
    @test MOIU._to_string(PLAIN, s) == "$(IN) ExponentialCone()"
    @test MOIU._to_string(LATEX, s) == "\\in \\text{ExponentialCone()}"
end

function test_feasibility()
    model = MOIU.Model{Float64}()
    @test sprint(print, model) == "Feasibility\n\nSubject to:\n"
    _string_compare(
        sprint(print, MOIU.latex_formulation(model)),
        raw"""
        $$ \begin{aligned}
        \text{feasibility}\\
        \text{Subject to}\\
        \end{aligned} $$""",
    )
    return
end

function test_min()
    model = MOIU.Model{Float64}()
    MOIU.loadfromstring!(model, "variables: x\nminobjective: x")
    @test sprint(print, model) == """
    Minimize VariableIndex:
     x

    Subject to:
    """
    _string_compare(
        sprint(print, MOIU.latex_formulation(model)),
        raw"""
        $$ \begin{aligned}
        \min\quad & x \\
        \text{Subject to}\\
        \end{aligned} $$""",
    )
    return
end

function test_max()
    model = MOIU.Model{Float64}()
    MOIU.loadfromstring!(model, "variables: x\nmaxobjective: x")
    @test sprint(print, model) == """
    Maximize VariableIndex:
     x

    Subject to:
    """
    _string_compare(
        sprint(print, MOIU.latex_formulation(model)),
        raw"""
        $$ \begin{aligned}
        \max\quad & x \\
        \text{Subject to}\\
        \end{aligned} $$""",
    )
    return
end

function test_model()
    model = MOIU.Model{Float64}()
    MOIU.loadfromstring!(
        model,
        """
        variables: x, y, z
        minobjective: x + 2 + 3.1*y + -1.2*z
        c1: x >= 0.1
        c2: y in ZeroOne()
        c2: z in Integer()
        c3: [x, y] in SecondOrderCone(2)
        c4: [1, x, y] in SecondOrderCone(2)
        c4: [1.0 * x * x, y, 1] in ExponentialCone()
        c4: [1, 1.0 * x * x, y] in ExponentialCone()
        c2: x in ZeroOne()
        c5: 2.0 * x * x + y + -1 * z <= 1.0
        c5: x + x >= 1.0
        c5: x + x in Interval(1.0, 2.0)
        c5: x + -1 * y == 0.0
        """,
    )
    @test sprint(print, model) == """
    Minimize ScalarAffineFunction{Float64}:
     2.0 + 1.0 x + 3.1 y - 1.2 z

    Subject to:

    ScalarAffineFunction{Float64}-in-EqualTo{Float64}
     0.0 + 1.0 x - 1.0 y == 0.0

    ScalarAffineFunction{Float64}-in-GreaterThan{Float64}
     0.0 + 2.0 x >= 1.0

    ScalarAffineFunction{Float64}-in-Interval{Float64}
     0.0 + 2.0 x $(IN) [1.0, 2.0]

    ScalarQuadraticFunction{Float64}-in-LessThan{Float64}
     0.0 + 1.0 y - 1.0 z + 2.0 x² <= 1.0

    VectorOfVariables-in-SecondOrderCone
     ┌ ┐
     │x│
     │y│
     └ ┘ $(IN) SecondOrderCone(2)

    VectorAffineFunction{Float64}-in-SecondOrderCone
     ┌           ┐
     │1.0        │
     │0.0 + 1.0 x│
     │0.0 + 1.0 y│
     └           ┘ $(IN) SecondOrderCone(2)

    VectorQuadraticFunction{Float64}-in-ExponentialCone
     ┌            ┐
     │0.0 + 1.0 x²│
     │0.0 + 1.0 y │
     │1.0         │
     └            ┘ $(IN) ExponentialCone()
     ┌            ┐
     │1.0         │
     │0.0 + 1.0 x²│
     │0.0 + 1.0 y │
     └            ┘ $(IN) ExponentialCone()

    VariableIndex-in-GreaterThan{Float64}
     x >= 0.1

    VariableIndex-in-Integer
     z $(IN) ℤ

    VariableIndex-in-ZeroOne
     x $(IN) {0, 1}
     y $(IN) {0, 1}
    """
end

function test_latex()
    model = MOIU.Model{Float64}()
    MOIU.loadfromstring!(
        model,
        """
        variables: x, y, z
        minobjective: x + 2 + 3.1*y + -1.2*z
        c1: x >= 0.1
        c2: y in ZeroOne()
        c2: z in Integer()
        c3: [x, y] in SecondOrderCone(2)
        c4: [1, x, y] in SecondOrderCone(2)
        c4: [1.0 * x * x, y, 1] in ExponentialCone()
        c4: [1, 1.0 * x * x, y] in ExponentialCone()
        c2: x in ZeroOne()
        c5: 2.0 * x * x + y + -1 * z <= 1.0
        c5: x + x >= 1.0
        c5: x + x in Interval(1.0, 2.0)
        c5: x + -1 * y == 0.0
        """,
    )
    _string_compare(
        sprint(
            io -> show(io, MIME("text/latex"), MOIU.latex_formulation(model)),
        ),
        raw"""
        $$ \begin{aligned}
        \min\quad & 2.0 + 1.0 x + 3.1 y - 1.2 z \\
        \text{Subject to}\\
         & \text{ScalarAffineFunction{Float64}-in-EqualTo{Float64}} \\
         & 0.0 + 1.0 x - 1.0 y = 0.0 \\
         & \text{ScalarAffineFunction{Float64}-in-GreaterThan{Float64}} \\
         & 0.0 + 2.0 x \ge 1.0 \\
         & \text{ScalarAffineFunction{Float64}-in-Interval{Float64}} \\
         & 0.0 + 2.0 x \in \[1.0, 2.0\] \\
         & \text{ScalarQuadraticFunction{Float64}-in-LessThan{Float64}} \\
         & 0.0 + 1.0 y - 1.0 z + 2.0 x^2 \le 1.0 \\
         & \text{VectorOfVariables-in-SecondOrderCone} \\
         & \begin{bmatrix}
        x\\
        y\end{bmatrix} \in \text{SecondOrderCone(2)} \\
         & \text{VectorAffineFunction{Float64}-in-SecondOrderCone} \\
         & \begin{bmatrix}
        1.0\\
        0.0 + 1.0 x\\
        0.0 + 1.0 y\end{bmatrix} \in \text{SecondOrderCone(2)} \\
         & \text{VectorQuadraticFunction{Float64}-in-ExponentialCone} \\
         & \begin{bmatrix}
        0.0 + 1.0 x^2\\
        0.0 + 1.0 y\\
        1.0\end{bmatrix} \in \text{ExponentialCone()} \\
         & \begin{bmatrix}
        1.0\\
        0.0 + 1.0 x^2\\
        0.0 + 1.0 y\end{bmatrix} \in \text{ExponentialCone()} \\
         & \text{VariableIndex-in-GreaterThan{Float64}} \\
         & x \ge 0.1 \\
         & \text{VariableIndex-in-Integer} \\
         & z \in \mathbb{Z} \\
         & \text{VariableIndex-in-ZeroOne} \\
         & x \in \{0, 1\} \\
         & y \in \{0, 1\} \\
        \end{aligned} $$""",
    )
    return
end

function test_latex_simplified()
    model = MOIU.Model{Float64}()
    MOIU.loadfromstring!(
        model,
        """
        variables: x, y, z
        minobjective: x + 2 + 3.1*y + -1.2*z
        c1: x >= 0.1
        c2: y in ZeroOne()
        c2: z in Integer()
        c3: [x, y] in SecondOrderCone(2)
        c4: [1, x, y] in SecondOrderCone(2)
        c4: [1.0 * x * x, y, 1] in ExponentialCone()
        c4: [1, 1.0 * x * x, y] in ExponentialCone()
        c2: x in ZeroOne()
        c5: 2.0 * x * x + y + -1 * z <= 1.0
        c5: x + x >= 1.0
        c5: x + x in Interval(1.0, 2.0)
        c5: x + -1 * y == 0.0
        """,
    )
    model_string = sprint() do io
        return MOIU._print_model(
            io,
            MOIU._PrintOptions(
                MIME("text/latex");
                simplify_coefficients = true,
                print_types = false,
            ),
            model,
        )
    end
    _string_compare(
        model_string,
        raw"""
        $$ \begin{aligned}
        \min\quad & 2 + x + 3.1 y - 1.2 z \\
        \text{Subject to}\\
         & x - y = 0 \\
         & 2 x \ge 1 \\
         & 2 x \in \[1, 2\] \\
         & y - z + 2 x^2 \le 1 \\
         & \begin{bmatrix}
        x\\
        y\end{bmatrix} \in \text{SecondOrderCone(2)} \\
         & \begin{bmatrix}
        1\\
        x\\
        y\end{bmatrix} \in \text{SecondOrderCone(2)} \\
         & \begin{bmatrix}
        x^2\\
        y\\
        1\end{bmatrix} \in \text{ExponentialCone()} \\
         & \begin{bmatrix}
        1\\
        x^2\\
        y\end{bmatrix} \in \text{ExponentialCone()} \\
         & x \ge 0.1 \\
         & z \in \mathbb{Z} \\
         & x \in \{0, 1\} \\
         & y \in \{0, 1\} \\
        \end{aligned} $$""",
    )
    return
end

function test_plain_simplified()
    model = MOIU.Model{Float64}()
    MOIU.loadfromstring!(
        model,
        """
        variables: x, y, z
        minobjective: x + -2 + 3.1*y + -1.2*z
        c1: x >= 0.1
        c2: y in ZeroOne()
        c2: z in Integer()
        c3: [x, y] in SecondOrderCone(2)
        c4: [1, x, y] in SecondOrderCone(2)
        c4: [1.0 * x * x, y, 1] in ExponentialCone()
        c4: [1, 1.0 * x * x, y] in ExponentialCone()
        c2: x in ZeroOne()
        c5: 2.0 * x * x + y + -1 * z <= 1.0
        c5: x + x >= 1.0
        c5: x + x in Interval(1.0, 2.0)
        c5: x + -1 * y == 0.0
        """,
    )
    model_string = sprint() do io
        return MOIU._print_model(
            io,
            MOIU._PrintOptions(
                MIME("text/plain");
                simplify_coefficients = true,
                print_types = false,
            ),
            model,
        )
    end
    @test model_string == """
    Minimize: -2 + x + 3.1 y - 1.2 z

    Subject to:
     x - y == 0
     2 x >= 1
     2 x $(IN) [1, 2]
     y - z + 2 x² <= 1
     ┌ ┐
     │x│
     │y│
     └ ┘ $(IN) SecondOrderCone(2)
     ┌ ┐
     │1│
     │x│
     │y│
     └ ┘ $(IN) SecondOrderCone(2)
     ┌  ┐
     │x²│
     │y │
     │1 │
     └  ┘ $(IN) ExponentialCone()
     ┌  ┐
     │1 │
     │x²│
     │y │
     └  ┘ $(IN) ExponentialCone()
     x >= 0.1
     z $(IN) ℤ
     x $(IN) {0, 1}
     y $(IN) {0, 1}
    """
end

function test_nlp()
    model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
    v = MOI.add_variables(model, 4)
    l = [1.1, 1.2, 1.3, 1.4]
    u = [5.1, 5.2, 5.3, 5.4]
    MOI.add_constraint.(model, v, MOI.GreaterThan.(l))
    MOI.add_constraint.(model, v, MOI.LessThan.(u))
    for i in 1:4
        MOI.set(model, MOI.VariableName(), v[i], "x[$i]")
    end
    lb, ub = [25.0, 40.0], [Inf, 40.0]
    evaluator = MOI.Test.HS071(true)
    block_data = MOI.NLPBlockData(MOI.NLPBoundsPair.(lb, ub), evaluator, true)
    MOI.set(model, MOI.NLPBlock(), block_data)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    @test sprint(print, model) == """
    Minimize Nonlinear:
     x[1] * x[4] * (x[1] + x[2] + x[3]) + x[3]

    Subject to:

    VariableIndex-in-GreaterThan{Float64}
     x[1] >= 1.1
     x[2] >= 1.2
     x[3] >= 1.3
     x[4] >= 1.4

    VariableIndex-in-LessThan{Float64}
     x[1] <= 5.1
     x[2] <= 5.2
     x[3] <= 5.3
     x[4] <= 5.4

    Nonlinear
     x[1] * x[2] * x[3] * x[4] >= 25.0
     x[1] ^ 2 + x[2] ^ 2 + x[3] ^ 2 + x[4] ^ 2 == 40.0
    """
    _string_compare(
        sprint(print, MOIU.latex_formulation(model)),
        raw"""
        $$ \begin{aligned}
        \min\quad & x_{1} \times x_{4} \times (x_{1} + x_{2} + x_{3}) + x_{3} \\
        \text{Subject to}\\
         & \text{VariableIndex-in-GreaterThan{Float64}} \\
         & x_{1} \ge 1.1 \\
         & x_{2} \ge 1.2 \\
         & x_{3} \ge 1.3 \\
         & x_{4} \ge 1.4 \\
         & \text{VariableIndex-in-LessThan{Float64}} \\
         & x_{1} \le 5.1 \\
         & x_{2} \le 5.2 \\
         & x_{3} \le 5.3 \\
         & x_{4} \le 5.4 \\
         & \text{Nonlinear} \\
         & x_{1} \times x_{2} \times x_{3} \times x_{4} \ge 25.0 \\
         & x_{1} ^ 2 + x_{2} ^ 2 + x_{3} ^ 2 + x_{4} ^ 2 = 40.0 \\
        \end{aligned} $$""",
    )
    return
end

function test_nlp_no_objective()
    model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
    v = MOI.add_variables(model, 4)
    for i in 1:4
        MOI.set(model, MOI.VariableName(), v[i], "x[$i]")
    end
    lb, ub = [25.0, 40.0], [Inf, 40.0]
    evaluator = MOI.Test.HS071(true)
    block_data = MOI.NLPBlockData(MOI.NLPBoundsPair.(lb, ub), evaluator, false)
    MOI.set(model, MOI.NLPBlock(), block_data)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), v[1])
    @test sprint(print, model) == """
    Minimize VariableIndex:
     x[1]

    Subject to:

    Nonlinear
     x[1] * x[2] * x[3] * x[4] >= 25.0
     x[1] ^ 2 + x[2] ^ 2 + x[3] ^ 2 + x[4] ^ 2 == 40.0
    """
    _string_compare(
        sprint(print, MOIU.latex_formulation(model)),
        raw"""
        $$ \begin{aligned}
        \min\quad & x_{1} \\
        \text{Subject to}\\
         & \text{Nonlinear} \\
         & x_{1} \times x_{2} \times x_{3} \times x_{4} \ge 25.0 \\
         & x_{1} ^ 2 + x_{2} ^ 2 + x_{3} ^ 2 + x_{4} ^ 2 = 40.0 \\
        \end{aligned} $$""",
    )
    return
end

function test_print_with_acronym()
    @test sprint(MOIU.print_with_acronym, "VecMathOptInterface") == "MOI"
    @test sprint(
        MOIU.print_with_acronym,
        "VecMathOptInterface.VecMathOptInterface",
    ) == "MOI.MOI"
    @test sprint(
        MOIU.print_with_acronym,
        "VecMathOptInterface.Utilities.VecMathOptInterface",
    ) == "MOIU.MOI"
    @test sprint(MOIU.print_with_acronym, "VecMathOptInterfaceXXBridges") ==
          "MOIXXBridges"
    @test sprint(MOIU.print_with_acronym, "VecMathOptInterface.BridgesXX") ==
          "MOIBXX"
    @test sprint(MOIU.print_with_acronym, "VecMathOptInterface.Test.x") == "MOIT.x"
    @test sprint(MOIU.print_with_acronym, "VecMathOptInterface.x.Test") ==
          "MOI.x.Test"
    @test sprint(MOIU.print_with_acronym, "VecMathOptInterface.Utilities.Test") ==
          "MOIU.Test"
    @test sprint(MOIU.print_with_acronym, "VecMathOptInterface.Utilities.Test") ==
          "MOIU.Test"
    return
end

end

TestPrint.runtests()


text:
# DataValueArray is dense and allows fast linear indexing.
Base.IndexStyle(::Type{<:DataValueArray}) = Base.IndexLinear()

# resolve ambiguity created by the two definitions that follow.
function Base.getindex(X::DataValueArray{T,N}) where {T,N}
    return X[1]
end

"""
    getindex{T,N}(X::DataValueArray{T,N}, I::Int...)

Retrieve a single entry from a `DataValueArray`. If the value in the entry
designated by `I` is present, then it will be returned wrapped in a
`DataValue{T}` container. If the value is missing, then this method returns
`DataValue{T}()`.
"""
@inline function Base.getindex(X::DataValueArray{T,N}, i::Int) where {T,N}
    @boundscheck checkbounds(X, i)
    if isbitstype(T)
        @inbounds r = DataValue{T}(X.values[i], !X.isna[i])
        return r
    else
        @inbounds if X.isna[i]
            return DataValue{T}()
        else
            return DataValue{T}(X.values[i])
        end
    end
end

"""
    getindex{T,N}(X::DataValueArray{T,N}, I::DataValue{Int}...)

Just as above, with the additional behavior that this method throws an error if
any component of the index `I` is missing.
"""
@inline function Base.getindex(X::DataValueArray{T,N}, I::DataValue{Int}...) where {T,N}
    any(isna, I) && throw(DataValueException())
    values = [ get(i) for i in I ]
    return getindex(X, values...)
end

"""
    setindex!(X::DataValueArray, v::DataValue, I::Int...)

Set the entry of `X` at position `I` equal to a `DataValue` value `v`. If
`v` is missing, then only `X.isna` is updated to indicate that the entry at
index `I` is null. If `v` is not null, then `X.isna` is updated to indicate
that the entry at index `I` is present and `X.values` is updated to store the
value wrapped in `v`.
"""
@inline function Base.setindex!(X::DataValueArray, v::DataValue, I::Int...)
    if isna(v)
        X.isna[I...] = true
    else
        X.isna[I...] = false
        X.values[I...] = get(v)
    end
    return v
end

"""
    setindex!(X::DataValueArray, v::Any, I::Int...)

Set the entry of `X` at position `I` equal to `v`. This method always updates
`X.isna` to indicate that the entry at index `I` is present and `X.values`
to store `v` at `I`.
"""
@inline function Base.setindex!(X::DataValueArray, v::Any, I::Int...)
    X.values[I...] = v
    X.isna[I...] = false
    return v
end

@inline function Base.setindex!(X::DataValueArray, v::DataValue{Union{}}, I::Int...)
    X.isna[I...] = true
    return v
end

"""
    unsafe_setindex_value!(X::DataValueArray, v, I::Int...)

Set the value entry of `X` at position `I` equal to `v`. This method
does not update X.isna.
"""
@inline function unsafe_setindex_value!(x::DataValueArray, v, I::Int...)
    x.values[I...] = v
    return
end

"""
    unsafe_setindex_isna!(X::DataValueArray, v, I::Int...)

Set the isna entry of `X` at position `I` equal to `v`. This method
does not update the actual value at position `I`.
"""
@inline function unsafe_setindex_isna!(x::DataValueArray, v, I::Int...)
    x.isna[I...] = v
    return
end

# DA TODO disabled
# function unsafe_getindex_notnull(X::DataValueArray, I::Int...)
#     return DataValue(getindex(X.values, I...))
# end

# DA TODO disabled
# function unsafe_getvalue_notnull(X::DataValueArray, I::Int...)
#     return getindex(X.values, I...)
# end

# DA TODO disabled
# function Base.checkindex(::Type{Bool}, inds::AbstractUnitRange, i::DataValue)
#     isna(i) ? throw(DataValueException()) : checkindex(Bool, inds, get(i))
# end

# DA TODO disabled
# function Base.checkindex{N}(::Type{Bool}, inds::AbstractUnitRange, I::DataValueArray{Bool, N})
#     any(isna, I) && throw(DataValueException())
#     checkindex(Bool, inds, I.values)
# end

# DA TODO disabled
# function Base.checkindex{T<:Real}(::Type{Bool}, inds::AbstractUnitRange, I::DataValueArray{T})
#     any(isna, I) && throw(DataValueException())
#     b = true
#     for i in 1:length(I)
#         @inbounds v = unsafe_getvalue_notnull(I, i)
#         b &= checkindex(Bool, inds, v)
#     end
#     return b
# end

# DA TODO disabled
# function Base.to_index(X::DataValueArray)
#     any(isna, X) && throw(DataValueException())
#     Base.to_index(X.values)
# end


text:
# Simultaneous function and Jacobian evaluation

"""
    value_jacobian(f::Function, x::AbstractVector)

Returns the value and Jacobian of the function `f` at the point `x`. The automatic differentiation package, Zygote, is used by default. To define a custom adjoint rule for a function to be used in a constraint or objective, use [`ChainRulesCore.jl`](https://github.com/JuliaDiff/ChainRulesCore.jl). The Jacobian is returned as an instance of `Adjoint` for cache efficiency in the optimizer.
"""
function value_jacobian(f::Function, x::AbstractVector)
    out, pullback = Zygote.pullback(f, x)
    if out isa Number
        grad = pullback(1.0)[1]
        return out, grad'
    else
        dim = getdim(f)
        jact = mapreduce(hcat, 1:dim; init = zeros(length(x), 0)) do i
            pullback(I(dim)[:,i])[1]
        end
        return out, jact'
    end
end
function value_jacobian(c::VectorOfFunctions, x::AbstractVector)
    vals_jacs = map(c.fs) do f
        value_jacobian(f, x)
    end
    vals = mapreduce(vcat, vals_jacs; init = zeros(0)) do val_jac
        val_jac[1]
    end
    jact = mapreduce(hcat, vals_jacs; init = zeros(length(x), 0)) do val_jac
        val_jac[2]'
    end
    return vals, jact'
end

"""
    value_jacobian_transpose(f::Function, x::AbstractVector)

Returns the value and transpose of the Jacobian of `f` at the point `x`. This calls [`value_jacobian`](@ref) and transposes the Jacobian.
"""
function value_jacobian_transpose(f::Function, x::AbstractVector)
    val, jac = value_jacobian(f, x)
    return val, jac'
end

"""
    value_gradient(f::Function, x::AbstractVector)

Returns the value and gradient of the scalar-valued function `f` at the point `x`. This is a convenience function for scalar-valued functions that simply calls [`value_jacobian_transpose`](@ref) on `f` and `x`.
"""
function value_gradient(f::Function, x::AbstractVector)
    return value_jacobian_transpose(f, x)
end


text:
# Load package
using Pkg;
using Distributions, Statistics, LinearAlgebra
using Plots

# Define population, N()
n_pop = 2_500; # Population size
n_iter = 5000; # Iteration of collision

# Distribution of income
μ_income = 50_000; # Mean income
σ_income = 12_000; # Standard deviation of income
𝐃 = Normal(μ_income, σ_income)

# Decider for selecting two agents to collide
𝐔 = Uniform(1, n_pop)

# vector of initial income of individual
val_inc_ini = deepcopy(rand(𝐃, n_pop));
# Equal splits
val_inc_trc = zeros(n_pop, n_iter); # Trace the changes in distribuiton
val_inc_trc[:,1] = val_inc_ini; # the very first distributional look
# Uniform random split
val_inc_trc_rand = zeros(n_pop, n_iter); # Trace the changes in distribuiton
val_inc_trc_rand[:,1] = val_inc_ini; # the very first distributional look

# Agent collide-split economy: two agents meet at a time, their wealth is equally distributed

for iter ∈ 2:n_iter
    # initialize the starting wealth distributio of the period
    val_inc_trc[:, iter] = val_inc_trc[:, (iter - 1)];
    # For each iteration, we select two agents who will collide
    agt_1, agt_2 = sample(1:n_pop, 2; replace = false)
    # collect their financial values
    val_inc_agt_1, val_inc_agt_2 = val_inc_trc[agt_1, iter], val_inc_trc[agt_2, iter];
    # divide up their wealth equally, and save
    val_inc_distributed = mean([val_inc_agt_1, val_inc_agt_2]);
    val_inc_trc[agt_1, iter] = val_inc_distributed;
    val_inc_trc[agt_2, iter] = val_inc_distributed;

    # Iteration 2 - uniform wealth redistribution, ranom split =========================
    # initialize the starting wealth distributio of the period
    val_inc_trc_rand[:, iter] = val_inc_trc_rand[:, (iter - 1)];
    # For each iteration, we select two agents who will collide
    agt_1, agt_2 = sample(1:n_pop, 2; replace = false)
    # collect their financial values
    total_wealth = val_inc_trc_rand[agt_1, iter] + val_inc_trc_rand[agt_2, iter];
    # How to split the wealth ..........................................................
    ρ = rand(1)[1]; #proportion to be split`
    # Redistribute wealth
    #   because of randomization, no need for random in order,
    val_inc_agt_1, val_inc_agt_2 = total_wealth * ρ, total_wealth * (1 - ρ)
    val_inc_trc_rand[agt_1, iter] = val_inc_agt_1;
    val_inc_trc_rand[agt_2, iter] = val_inc_agt_2;
end

println("Mean - true: $(μ_income); std - true: $(σ_income)")
println("Mean: $(mean(val_inc_trc[:,end])); std: $(std(val_inc_trc[:,end]))")
@gif for i ∈ 1:10:n_iter
    histogram(val_inc_trc[:,i], title = "Iteration: $(i)", xlims = (0,1e5))
end every 5

println("Mean - true: $(μ_income); std - true: $(σ_income)")
println("Mean: $(mean(val_inc_trc_rand[:,end])); std: $(std(val_inc_trc_rand[:,end]))")
@gif for i ∈ 1:10:n_iter
    histogram(val_inc_trc_rand[:,i], title = "Iteration: $(i)", xlims = (0,2e5))
end every 5

max_stars_repo_name: justinjoliver/simulation-discrete

max_issues_repo_name: justinjoliver/simulation-discrete

max_forks_repo_name: justinjoliver/simulation-discrete


text:
using Colors

import Ls
import Ls.LineSearch
import GdOpt
import GdOpt.GdOptimizer
import Obj
import Plotting

# Quick test functions for GD.
mnist_gd() = mnist(GdOpt.gd(), Ls.bt())
mnist_gd_boost() = mnist_boost(GdOpt.gd(), Ls.bt())

# Quick test functions for SGD.
mnist_sgd(batchsize=1) = mnist(GdOpt.sgd(), Ls.sbt(), ϵ=0.0, batchsize=batchsize)
mnist_sgd_boost(batchsize=1) = mnist_boost(GdOpt.sgd(), Ls.sbt(), ϵ=0.0, batchsize=batchsize)

# Quick test functions for SVRG.
mnist_svrg(; batchsize=10, estiter=10, strategy=:last, ϵ=0.0) =
    mnist(GdOpt.svrg(estiter, strategy), Ls.sbt(), batchsize=batchsize, ϵ=ϵ, assumedgrad=true)
mnist_svrg_boost(; batchsize=10, estiter=10, strategy=:last, ϵ=0.0) =
    mnist_boost(GdOpt.svrg(estiter, strategy), Ls.sbt(), batchsize=batchsize, ϵ=ϵ, assumedgrad=true)

"""
    mnist_massive([basefolder; seeds, gdopts, batchsizes])

    Massive evaluation using all of the optimizers with different parameters and seeds.
    In default configuration, this may run up to 11hrs.
"""
function mnist_massive(basefolder::ASCIIString="seven_vs_all";
                       seeds::UnitRange{Int}      = 1337:(1337+9),
                       gdopts::Array{GdOptimizer} = [GdOpt.gd(), GdOpt.sgd(),
                                                     GdOpt.svrg(5, :last),  GdOpt.svrg(5, :avg),
                                                     GdOpt.svrg(10, :last), GdOpt.svrg(10, :avg),
                                                     GdOpt.svrg(25, :last), GdOpt.svrg(25, :avg)],
                       batchsizes::Array{Int,1}   = [1, 10, 25, 100])

    println("\nSeeds:       $seeds")
    println("Optimizers:  ", join([gdopt.name for gdopt in gdopts], ", "))
    println("Batch sizes: ", join(batchsizes, ", "))

    gdcontained = false
    for gdopt in gdopts
        if gdopt.name == GdOpt.gd().name
            gdcontained = true
            break
        end
    end

    # let user confirm his madness
    numexperiments = length(seeds) * length(gdopts) * length(batchsizes) * 4 + gdcontained * length(seeds) * 4
    print("\nThis results in $numexperiments experiments, ",
          "that may last up to $(numexperiments/120)hrs and use up to $(numexperiments*6)MB of storage.\nPlease confirm: y/N? ")
    if !startswith(readline(STDIN), "y")
        return nothing
    end

    # iteration printer
    expiter = 0
    printiter(add::Int) = println("\n", repeat("-", 25),
                                  " Finished iteration $(expiter += add), ",
                                  "still $((numexperiments-expiter)) ",
                                  "($(floor((numexperiments-expiter)/6)/20)hrs) to go ",
                                  repeat("-", 25), "\n")

    # go for it
    for seed in seeds

        X_train, y_train, X_test, y_test = _mnist_data(basefolder*"/"*string(seed))  # read data

        df = LsoBase.new_acc()
        for gdopt in gdopts

            if gdopt.name == "gd"
                try
                    df = vcat(df, mnist(gdopt, Ls.bt(), basefolder, seed,
                                        X_train=X_train, y_train=y_train,
                                        X_test=X_test, y_test=y_test))
                catch e
                    warn(e)
                    writedlm("./results/$basefolder/$seed/$(gdopt.name)$((batchsize>0)?batchsize:"")_$(Ls.bt().name).dlm", string(e))
                end
                printiter(1)

                try
                    df = vcat(df, mnist_boost(gdopt, Ls.bt(), basefolder, seed,
                                              X_train=X_train, y_train=y_train,
                                              X_test=X_test, y_test=y_test))
                catch e
                    warn(e)
                    writedlm("./results/$basefolder/$seed/boost_$(gdopt.name)$((batchsize>0)?batchsize:"")_$(Ls.bt().name).dlm", string(e))
                end
                printiter(3) # boosting has 3 iterations
            end

            for batchsize in batchsizes
                ϵ=1e-3
                if (gdopt.name == "sgd" && batchsize != 100) || (gdopt.name == "svrg" && batchsize == 1)
                    ϵ=0.0
                end

                try
                    df = vcat(df, mnist(gdopt, Ls.sbt(), basefolder, seed, batchsize=batchsize, ϵ=ϵ,
                                        X_train=X_train, y_train=y_train,
                                        X_test=X_test, y_test=y_test))
                catch e
                    warn(e)
                    writedlm("./results/$basefolder/$seed/$(gdopt.name)$((batchsize>0)?batchsize:"")_$(Ls.sbt().name).dlm", string(e))
                end
                printiter(1)

                try
                    df = vcat(df, mnist_boost(gdopt, Ls.sbt(), basefolder, seed, batchsize=batchsize, ϵ=ϵ,
                                              X_train=X_train, y_train=y_train,
                                              X_test=X_test, y_test=y_test))
                catch e
                    warn(e)
                    writedlm("./results/$basefolder/$seed/boost_$(gdopt.name)$((batchsize>0)?batchsize:"")_$(Ls.sbt().name).dlm", string(e))
                end
                printiter(3) # boosting has 3 iterations
            end

        end
        LsoBase.store(df, "./results/$basefolder/$seed/acc.csv")
    end

    println("\nDone!\n")

end

"""
    mnist(gdopt, ls, [, basefolder, seed; ϵ, maxtime, batchsize, assumedgrad, storage])

    Evaluate optimizer on MNIST data with Logistic Regression.
"""
function mnist(gdopt::GdOptimizer, ls::LineSearch, basefolder::ASCIIString="seven_vs_all", seed::Int=1337;
               λ::Float64=0.0, ϵ::Float64=1e-3, maxtime::Float64=30.0, batchsize::Int=-1,
               assumedgrad=false, storage=true,
               X_train::Union{Void,Array{Float64,2}}=nothing, y_train::Union{Void,Array{Float64,1}}=nothing,
               X_test::Union{Void,Array{Float64,2}}=nothing,  y_test::Union{Void,Array{Float64,1}}=nothing)

    println("\nRunning Eval.mnist with seed $seed on $(gdopt.name) with $(ls.name) and batchsize $batchsize...")
    srand(seed)
    X_train, y_train, X_test, y_test = _mnist_data(basefolder*"/"*string(seed), X_train, y_train, X_test, y_test)  # read data

    # optimize
    inf, obj, w, acc_train, acc_test, time, iterrate = _mnist_opt(
            gdopt, ls, X_train, y_train, X_test, y_test,
            λ, ϵ, maxtime, batchsize
    )

    # plotting and storage
    if storage
        filename = "./results/$basefolder/$seed/$(gdopt.name)$((batchsize>0)?batchsize:"")_$(ls.name)"
        LsoBase.store(inf, filename * ".csv")
        Plotting.draw_plot(Plotting.plot_inf(inf, assumedgrad), filename * ".pdf")
    end

    println("")
    df = LsoBase.new_acc()
    LsoBase.push_acc!(df, gdopt.name, ls.name, batchsize, seed, basefolder, "LogReg",
                      size(X_train)[1], acc_train, acc_test, time, iterrate)
    return df

end

"""
    mnist(gdopt, ls, [, basefolder, seed; ϵ, maxtime, batchsize, assumedgrad, storage])

    Evaluate optimizer on MNIST data with historic boosting.
"""
function mnist_boost(gdopt::GdOptimizer, ls::LineSearch, basefolder::ASCIIString="seven_vs_all", seed::Int=1337;
                     λ::Float64=0.0, ϵ::Float64=1e-3, maxtime::Float64=30.0, batchsize::Int=-1,
                     assumedgrad=false, frac1=.5, storage=true,
                     X_train::Union{Void,Array{Float64,2}}=nothing, y_train::Union{Void,Array{Float64,1}}=nothing,
                     X_test::Union{Void,Array{Float64,2}}=nothing,  y_test::Union{Void,Array{Float64,1}}=nothing)

    println("\nRunning Eval.mnist_boost with seed $seed on $(gdopt.name) with $(ls.name) and batchsize $batchsize...")
    srand(seed)
    X_train, y_train, X_test, y_test = _mnist_data(basefolder*"/"*string(seed), X_train, y_train, X_test, y_test)  # read data

    ################
    # ITERATION 1
    ################

    # random subset of relative size frac1
    train1 = randperm(length(y_train))[1:convert(Int, floor(frac1*length(y_train)))]
    X_train1 = X_train[train1,:]
    y_train1 = y_train[train1]

    # optimize
    inf1, obj1, w1, acc_train1, acc_test1, time1, iterrate1 = _mnist_opt(
            gdopt, ls, X_train1, y_train1, X_test, y_test,
            λ, ϵ, maxtime, batchsize
    )

    ################
    # ITERATION 2
    ################

    # find subset that is classified with 50% acc
    y_pred1    = Obj.logreg_predict(w1, X_train)
    correct1   = shuffle(find(y_train .== y_pred1))
    false1     = shuffle(find(y_train .!= y_pred1))
    train2size = min(length(correct1), length(false1))
    X_train2   = getindex(X_train, vcat(correct1[1:train2size], false1[1:train2size]), :)
    y_train2   = getindex(y_train, vcat(correct1[1:train2size], false1[1:train2size]))

    # optimize
    inf2, obj2, w2, acc_train2, acc_test2, time2, iterrate2 = _mnist_opt(
            gdopt, ls, X_train2, y_train2, X_test, y_test,
            λ, ϵ, maxtime, batchsize
    )

    ################
    # ITERATION 3
    ################

    # find subset for which predictors decide differently
    y_pred2 = Obj.logreg_predict(w2, X_train)
    X_train3 = X_train[(y_pred1 .!= y_pred2), :]
    y_train3 = y_train[(y_pred1 .!= y_pred2)]

    # optimize
    inf3, obj3, w3, acc_train3, acc_test3, time3, iterrate3 = _mnist_opt(
            gdopt, ls, X_train3, y_train3, X_test, y_test,
            λ, ϵ, maxtime, batchsize
    )

    ################
    # MAJORITY VOTE
    ################

    # training prediction
    y_trainpred1 = Obj.logreg_predict(w1, X_train)
    y_trainpred2 = Obj.logreg_predict(w2, X_train)
    y_trainpred3 = Obj.logreg_predict(w3, X_train)
    y_trainpred  = vec(sign(sum([y_trainpred1 y_trainpred2 y_trainpred3], 2)))

    # test prediction
    y_testpred1  = Obj.logreg_predict(w1, X_test)
    y_testpred2  = Obj.logreg_predict(w2, X_test)
    y_testpred3  = Obj.logreg_predict(w3, X_test)
    y_testpred   = vec(sign(sum([y_testpred1 y_testpred2 y_testpred3], 2)))

    # acc
    acc_train = LsoBase.acc(y_train, y_trainpred)
    acc_test  = LsoBase.acc(y_test, y_testpred)
    timesum   = sum([time1, time2, time3])
    numex     = size(unique(vcat(X_train1, X_train2, X_train3), 1))[1]
    iterrate  = mean([iterrate1, iterrate2, iterrate3])

    # conclusion
    df = LsoBase.new_acc()
    LsoBase.push_acc!(df, gdopt.name, ls.name, batchsize, seed, basefolder, "LogReg B1",
                      size(X_train1)[1], acc_train1, acc_test1, time1,   iterrate1)
    LsoBase.push_acc!(df, gdopt.name, ls.name, batchsize, seed, basefolder, "LogReg B2",
                      size(X_train2)[1], acc_train2, acc_test2, time2,   iterrate2)
    LsoBase.push_acc!(df, gdopt.name, ls.name, batchsize, seed, basefolder, "LogReg B3",
                      size(X_train3)[1], acc_train3, acc_test3, time3,   iterrate3)
    LsoBase.push_acc!(df, gdopt.name, ls.name, batchsize, seed, basefolder, "Boosting B1-3",
                      numex,             acc_train,  acc_test,  timesum, iterrate)
    println(string(df), "\n")

    # concatenate inf
    inf2[:time] += time1
    inf2[:iter] += inf1[end, :iter]
    inf3[:time] += time2 + time1
    inf3[:iter] += inf2[end, :iter]
    inf = vcat(inf1, inf2, inf3)

    # plotting and storage
    if storage
        filename = "./results/$basefolder/$seed/boost_$(gdopt.name)$((batchsize>0)?batchsize:"")_$(ls.name)"
        LsoBase.store(inf, filename * ".csv")
        Plotting.draw_plot(Plotting.plot_inf(inf, assumedgrad, maxtime=3*maxtime), filename * ".pdf")
    end

    println("")
    return df

end

function _mnist_data(folder::ASCIIString,
                     X_train::Union{Void,Array{Float64,2}}=nothing, y_train::Union{Void,Array{Float64,1}}=nothing,
                     X_test::Union{Void,Array{Float64,2}}=nothing,  y_test::Union{Void,Array{Float64,1}}=nothing)

    if X_train == nothing || y_train == nothing || X_test == nothing || y_test == nothing
        println("Reading data from ./data/$folder...")
    else
        println("Data already present - nothing to read.")
    end

    if X_train == nothing
        X_train = readdlm("./data/$folder/X_train.dlm")
    end
    if y_train == nothing
        y_train = vec(readdlm("./data/$folder/y_train.dlm"))
    end
    if X_test == nothing
        X_test  = readdlm("./data/$folder/X_test.dlm")
    end
    if y_test == nothing
        y_test  = vec(readdlm("./data/$folder/y_test.dlm"))
    end

    return X_train, y_train, X_test, y_test

end

function _mnist_opt(gdopt::GdOptimizer, ls::LineSearch, X_train, y_train, X_test, y_test,
                    λ::Float64, ϵ::Float64, maxtime::Float64, batchsize::Int)

    w_0 = zeros(784)
    obj = Obj.logreg(X_train, y_train, λ)

    if size(X_train)[1] == 0
        println("\nSkipping experiment because of empty training set...")
        inf = LsoBase.new_inf()
        LsoBase.push_inf!(inf, w_0, 1e4, 1e4, 0, 0, 0.0)
        LsoBase.push_truth!(inf, obj.f, obj.g)  # true f and opt values
        return inf, obj, w_0, 0.0, 0.0, 0.0, 0.0
    else

        println("\nNow considering $(size(X_train)[1]) training examples...")

        # optimize

        @time inf = GdOpt.opt(gdopt, ls, obj, w_0,
                              ϵ=ϵ, maxtime=maxtime, batchsize=batchsize)
        w = inf[end, :w]
        time = inf[end, :time]
        iterrate = inf[end, :iter] / inf[end, :time]

        LsoBase.push_truth!(inf, obj.f, obj.g)  # true f and opt values

        # acc
        acc_train = LsoBase.acc(y_train, Obj.logreg_predict(w, X_train))
        acc_test  = LsoBase.acc(y_test, Obj.logreg_predict(w, X_test))
        println("\n  Training set acc: $acc_train")
        println(  "      Test set acc: $acc_test")
        println(  "  Iterations / sec: $iterrate\n")

        return inf, obj, w, acc_train, acc_test, time, iterrate

    end

end

max_stars_repo_name: mirkobunse/lso

max_issues_repo_name: mirkobunse/lso

max_forks_repo_name: mirkobunse/lso


text:
using Javis

function object(r,p=O, color="black")
    sethue(color)
    circle(p, r, :fill)
    return p
end

function connector!(connection,p1, p2, color)
    sethue(color)
    push!(connection,[p1,p2])
    map(x->line(x[1],x[2], :stroke),connection)
end

function ground(args...)
    background("black") # canvas background
    sethue("white") # pen color
    # circle(O,50,:fill)
end

function circ(p = O, color = "black", action = :fill, radius = 25, edge = "solid")
    sethue(color)
    setdash(edge)
    circle(p, radius, action)
end

function make_animation()

    # to store the connectors
    connection = []

    frames = 1000

    # setup the video
    myvideo = Video(500,500)
    Background(1:frames,ground)

    # add the objects
    earth = Object(1:frames,(args...)->object(5,O,"blue"),Point(200,0))
    venus = Object(1:frames,(args...)->object(4,O,"red"),Point(144,0))
    #sun = Object(1:frames,(aegs...)->object(50,O,"yellow"),Point(0,0))

    # draw the orbits
    earth_orbit = Object((args...) -> circ(O, "white", :stroke, 200))
    venus_orbit = Object((args...) -> circ(O, "white", :stroke, 144))

    # move the planets
    act!(earth, Action(anim_rotate_around(12.5*2π, O)))
    act!(venus, Action(anim_rotate_around(12.5*2π*(224.7/365), O)))

    # draw the connectors
    Object(1:frames, (args...)->connector!(connection,pos(earth), pos(venus), "#f05a4f"))

    # render
    render(myvideo,pathname = "cosmic_dance.gif")
end

make_animation()

max_stars_repo_name: Ved-Mahajan/JavisIsAwsome

max_issues_repo_name: Ved-Mahajan/JavisIsAwsome

max_forks_repo_name: Ved-Mahajan/JavisIsAwsome


text:
# ---
# title: 15. 3Sum
# id: problem15
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Two Pointers
# link: <https://leetcode.com/problems/3sum/description/>
# hidden: true
# ---
#
# Given an array `nums` of _n_ integers, are there elements _a_ , _b_ , _c_ in
# `nums` such that _a_ \+ _b_ \+ _c_ = 0? Find all unique triplets in the array
# which gives the sum of zero.
#
# Notice that the solution set must not contain duplicate triplets.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [-1,0,1,2,-1,-4]
#     Output: [[-1,-1,2],[-1,0,1]]
#
#
# **Example 2:**
#
#
#
#     Input: nums = []
#     Output: []
#
#
# **Example 3:**
#
#
#
#     Input: nums = [0]
#     Output: []
#
#
#
#
# **Constraints:**
#
#   * `0 <= nums.length <= 3000`
#   * `-105 <= nums[i] <= 105`
#
#
## @lc code=start
using LeetCode

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


text:
# Autogenerated wrapper script for Luna_jll for x86_64-linux-musl-cxx11
export libluna, luna

using FFTW_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("Luna")
JLLWrappers.@declare_library_product(libluna, "libluna.so")
JLLWrappers.@declare_executable_product(luna)
function __init__()
    JLLWrappers.@generate_init_header(FFTW_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libluna,
        "lib/libluna.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

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

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"Evaluate the Laplace 2D single layer potential kernel"
function laplace_slp(x, y)
    T = eltype(x)
    if norm(x-y) < eps(T)
         return zero(Complex{T})
    else
        return complex(1/2/pi*log(norm(x-y)))
    end
end

export Laplace_SLP_2D
"The single layer potential kernel of the Laplace equation."
struct Laplace_SLP_2D{T} <: BoundaryKernel
end

Laplace_SLP_2D() = Laplace_SLP_2D{Float64}()

(kernel::Laplace_SLP_2D)(t, tau, param, x, y) = laplace_slp(x, y)
(kernel::Laplace_SLP_2D)(x, tau, param, y) = laplace_slp(x, y)

is_symmetric(::Laplace_SLP_2D) = true

singularity(::Laplace_SLP_2D) = LogSingularDiagonal()

"Evaluate the Laplace 2D double layer potential kernel"
laplace_dlp(x, y, normal_y, z = norm(x-y)) =
    (normal_y' * (x-y)) / z^2/2/pi

"Evaluate the Laplace 2D double layer potential kernel"

function laplace_dlp_kernel(x, y, tau, param)
    z = norm(x-y)
    if abs(z) > eps(eltype(x))
        return Complex{Float64}(laplace_dlp(x, y, normal(param, tau), z))
    else
        return Complex{Float64}(0)
    end
end

export Laplace_DLP_2D

"The double layer potential kernel of the Laplace equation in 2D."
struct Laplace_DLP_2D{T} <: BoundaryKernel
end

Laplace_DLP_2D() = Laplace_DLP_2D{Float64}()

BasisFunctions.name(kernel::Laplace_DLP_2D) = "2D Laplace double layer potential kernel"

(kernel::Laplace_DLP_2D)(t, tau, param, x, y) =
    laplace_dlp_kernel(x, y, tau, param)
(kernel::Laplace_DLP_2D)(x, tau, param, y) =
    laplace_dlp_kernel(x, y, tau, param)

is_symmetric(::Laplace_DLP_2D) = false

# The kernel is continuous but its derivative is singular
singularity(::Laplace_DLP_2D) = LogSingularDiagonal()

"Adjoint Double layer potential for Laplace equation"

export Laplace_adjDLP_2D
"The adjoint double layer potential kernel of the Laplace equation in 2D."
struct Laplace_adjDLP_2D{T} <: BoundaryKernel
end

laplace_adjdlp(x, y, normal_x, z = norm(x-y)) =
    (normal_x' * (y-x)) / z^2/2/pi

"Evaluate the adjoint of the Laplace 2D double layer potential kernel"

function laplace_adjdlp_kernel(x, y, t, param)
    z = norm(x-y)
    if abs(z) > 10*eps(eltype(x))
        return Complex{Float64}(laplace_adjdlp(x, y, normal(param, t), z))
    else
        return Complex{Float64}(0)
    end
end

Laplace_adjDLP_2D() = Laplace_adjDLP_2D{Float64}()

BasisFunctions.name(kernel::Laplace_adjDLP_2D) = "Adjoint of 2D Laplace double layer potential kernel"

(kernel::Laplace_adjDLP_2D)(t, tau, param, x, y) =
    laplace_adjdlp_kernel(x, y, t, param)

is_symmetric(::Laplace_adjDLP_2D) = false

# The kernel is continuous but its derivative is singular
singularity(::Laplace_adjDLP_2D) = LogSingularDiagonal()
## Some boundary conditions

export laplace_bcond_harmonic, laplace_param_bcond_harmonic, laplace_parboundary_condition_potential_flow

harmonic_poly(z) = z^2-z+2

boundary_condition_harmonic(x) = real(harmonic_poly(x[1]+im*x[2]))

"Return a field function that evaluates to a plane wave."
laplace_bcond_harmonic() = (x,y) -> boundary_condition_harmonic(SVector(x,y))

"Return a parametric boundary function that evaluates to a plane wave."
laplace_param_bcond_harmonic(param) =
    t -> boundary_condition_harmonic(applymap(param, t))

"Return a parametric boundary function corresponding to the unperturbed flow at r=∞"
laplace_parboundary_condition_potential_flow(param, velocity) =
    t -> velocity'*normal(param, t)


text:
function EventBasedScheduler(arg0::MeasurementBuilder, arg1::DatesSelector, arg2::Propagator, arg3::EventDetector, arg4::SignSemantic)
    return EventBasedScheduler((MeasurementBuilder, DatesSelector, Propagator, EventDetector, SignSemantic), arg0, arg1, arg2, arg3, arg4)
end

function generate(obj::EventBasedScheduler, arg0::List)
    return jcall(obj, "generate", SortedSet, (List,), arg0)
end


text:
using DifferentialEquations
using Plots

n_parasites = 100;

c = 4.0;
ux = 0.2;
ux1 = fill(0.2, n_parasites); #le ux et uy 1000 à cause de la β # une autre façon :[0.2 for x in 1:1000]
Random.seed!(1234);
uy = rand(200:1000, n_parasites)/1000;
Random.seed!(1235);
r1 = rand(Float64, n_parasites)
Random.seed!(1236);
r2 = rand(Float64, n_parasites)
Random.seed!(1237);
r3 = rand(0.0:1000, n_parasites)/1000;
bx = 1.0;
by = bx .* r1 .* (1 .- r1 .* r2);
V0 = (by .* ux) ./ (bx .* uy);
α = 1 .- V0;
βy = r1 .-(α.* by)/bx;
ey = bx .* (1 .- r3) .* (1 .- (r1 .* r2));

Y = zeros(Float64, length(ey));
Y[1] = 1.0;
X0 = 10.0;

# differential equations formula for finding densities of each pop at next time step
function fonction(u, p, t)
    x = u[1]
    y = u[2:end]
    regul = (1-(sum(u))/(p.K))
    dx = (p.bx*x + sum(p.ey.*y))*regul - p.ux*x-p.c*sum(p.βy.*y)*x
    dy = p.by .* y .* regul .- p.uy .* y .+ p.c .* p.βy .* x .* y
    return vcat(dx, dy)
end

debut = 0.0;
duree = 1000.0;
fin = debut + duree;
N = zeros(Float64, (n_parasites+1, (n_parasites-1)*Int(duree)+1));
new_U = vcat(X0, Y);
parameters = (bx = bx, βy = βy, ey = ey, c = c, K = 80.0, ux = ux, by = by, uy = uy);

# each strain introduction (1000x)
@progress "Simulation" for i in 2:length(Y)
    # initial conditions
    prob = ODEProblem(fonction, new_U, (debut,fin), parameters)
    solution = solve(prob, saveat=debut:1.0:fin)
    for t in eachindex(solution.t)
        pop = solution.u[t]
        for i in 1:n_parasites
            if (pop.<0)[i]
                pop[i] = 0
            end
        end
        N[:,Int(solution.t[t]+1)] = pop
    end

    # set conditions & new parasite for next loop
    global new_U = solution[end]
    new_y = findfirst(x -> x == 0.0, new_U)
    new_U[new_y] = 1.0

    # set limits for next loop
    global debut = fin
    #global fin = Float64((i+1).* 1000.0)
    global fin = debut + duree
end

Np = N';

#strains that survive at each time step
survival = (Np.>0.0)[:,2:end];

#mean β
survived_βy = survival.*βy';
avg_survived_βy = sum(survived_βy; dims=2)./sum(survival; dims=2);
βy_avg = avg_survived_βy;

#weighted β
avg_w_survived_βy = sum(Np[:,2:end].*βy'; dims=2)./sum(Np[:,2:end]; dims=2);
βi_avg = avg_w_survived_βy;

#mean virulence
Vir = (1 .- (by .+ ey).*ux./(bx.*uy));

survived_Vir = survival.*Vir';
avg_survived_Vir = sum(survived_Vir; dims=2)./sum(survival; dims=2);
Vir_avg = avg_survived_Vir;

#weighted virulence
avg_w_survived_Vir = sum(Np[:,2:end].*Vir'; dims=2)./sum(Np[:,2:end]; dims=2);
Vir_avg_w = avg_w_survived_Vir;

#plot(βy_avg, title = "Mean virulence and beta", label = "beta")
#plot!(Vir_avg, label = "virulence")
plot(βi_avg,c=:black, title = "Virulence and beta", label = "Beta", xlabel = "Time", ylabel = "Mean virulence & \n Mean Beta", ylims = (0,1))
plot!(Vir_avg_w,c=:blue, label = "Virulence")

# png("Figure 3/graph_3h.png")


text:
type SDF_PID_VD <: SDFDynamicSystem
    omg::Float64
    zta::Float64
    mu::Float64
    xi::Float64
    kpa::Float64
    alpha::Float64
    u0::Float64
    v0::Float64
    ud0::Float64
    vd0::Float64

    SDF_PID_VD(omg::Float64, zta::Float64, mu::Float64, xi::Float64, kpa::Float64, alpha::Float64) = new(omg, zta, mu, xi, kpa, alpha, 0.0, 0.0, 0.0, 0.0)
    SDF_PID_VD(omg::Float64, zta::Float64, mu::Float64, xi::Float64, kpa::Float64) = new(omg, zta, mu, xi, kpa, 1.0, 0.0, 0.0, 0.0, 0.0)
end

function systemEquation(s::SDF_PID_VD, p::Function)
    function dydt(y::AbstractArray{Float64,1}, t::Float64)
        u, v, ud, vd = y
        omg2 = s.omg*s.omg
        dvdt = p(t)-2.0*s.zta*s.omg*v-omg2*u-s.kpa*omg2*(u-ud)

        if s.alpha == 1.0
            dvddt = (s.kpa*omg2*(u-ud) - 2.0*s.xi*s.omg*vd)/s.mu
        else
            dvddt = (s.kpa*omg2*(u-ud) - 2.0*s.xi*s.omg*sign(vd)*abs(vd)^s.alpha)/s.mu
        end
        [v; dvdt; vd; dvddt]
    end
    dydt
end

function response(s::SDF_PID_VD, p::Function, t::AbstractArray{Float64,1})
    y0 = [s.u0; s.v0; s.ud0; s.vd0]

    f = systemEquation(s, p)

    y = ode4(f, y0, t)

    u = y[:,1]
    v = y[:,2]
    ud = y[:,3]
    vd = y[:,4]
    fdz = 2.0*s.zta*s.omg*v
    fdx = s.kpa*s.omg*s.omg*(u-ud)
    fdd = 2.0*s.xi*s.omg*sign(vd).*abs(vd).^s.alpha
    a = map(p,t)-fdz-fdx-s.omg*s.omg*u

    [u v a fdz fdx ud fdd]
end

max_stars_repo_name: Panchatantra/DSR

max_issues_repo_name: Panchatantra/DSR

max_forks_repo_name: Panchatantra/DSR


text:
function mark_all!(r::AbstractChessRobot)
    side1 = Sud
    side2 = West
    h, w = go_to_corner!(r, side1, side2)
    s = Ost
    mark_snake!(r, s)
    go_to_corner!(r, side1, side2)
    move_back(r, inverse_side(side1), inverse_side(side2), h, w)
end

max_stars_repo_name: dmitrii-a-piskunov/mirea-progs

max_issues_repo_name: dmitrii-a-piskunov/mirea-progs

max_forks_repo_name: dmitrii-a-piskunov/mirea-progs