text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

function _add_production_piecewise_linear_eqs!(
    model::JuMP.Model,
    g::Unit,
    formulation_prod_vars::Gar1962.ProdVars,
    formulation_pwl_costs::CarArr2006.PwlCosts,
    formulation_status_vars::StatusVarsFormulation,
)::Nothing
    eq_prod_above_def = _init(model, :eq_prod_above_def)
    eq_segprod_limit = _init(model, :eq_segprod_limit)
    segprod = model[:segprod]
    gn = g.name

    # Gar1962.ProdVars
    prod_above = model[:prod_above]

    K = length(g.cost_segments)
    for t in 1:model[:instance].time
        gn = g.name
        for k in 1:K
            # Equation (45) in Kneuven et al. (2020)
            # NB: when reading instance, UnitCommitment.jl already calculates
            #     difference between max power for segments k and k-1 so the
            #     value of cost_segments[k].mw[t] is the max production *for
            #     that segment*
            eq_segprod_limit[gn, t, k] = @constraint(
                model,
                segprod[gn, t, k] <= g.cost_segments[k].mw[t]
            )

            # Also add this as an explicit upper bound on segprod to make the
            # solver's work a bit easier
            set_upper_bound(segprod[gn, t, k], g.cost_segments[k].mw[t])

            # Definition of production
            # Equation (43) in Kneuven et al. (2020)
            eq_prod_above_def[gn, t] = @constraint(
                model,
                prod_above[gn, t] == sum(segprod[gn, t, k] for k in 1:K)
            )

            # Objective function
            # Equation (44) in Kneuven et al. (2020)
            add_to_expression!(
                model[:obj],
                segprod[gn, t, k],
                g.cost_segments[k].cost[t],
            )
        end
    end
end


text:
module Clustering
    using Devectorize
    using Distance
    using MLBase

    import Base.show

    export k_means
    export kmeans, kmeans!, kmeans_opts, update!

    include("seeding.jl")
    include("kmeans.jl")

    export AffinityPropagationOpts
    export affinity_propagation
    include("affprop.jl")
end

ext: jl


text:
include("..\\libs\\FunctionsRobot.jl")

function krest!(r::Robot)
    for side in instances(HorizonSide)
        if isborder(r, side)
            continue
        end
        move!(r, side)
        marksLine!(r, side)
        moveWhileMarker!(r, reversSide(side))
    end
    putmarker!(r)
end


text:
using Test

@testset "Material Parser" begin
    @test NeXLCore._mp_level5("SiO2")  == Dict{Element, Int64}(n"Si" => 1, n"O" => 2)
    @test NeXLCore._mp_level5("Al2O3") == Dict{Element, Int64}(n"Al" => 2, n"O" => 3)
    @test NeXLCore._mp_level5("NaAlSi3O8") == Dict{Element, Int64}(n"Na" => 1, n"Al" => 1, n"Si" => 3, n"O" => 8)
    @test NeXLCore._mp_level5("") == Dict{Element,Int64}()
    @test_throws ErrorException NeXLCore._mp_level5("SiLO2")
    @test_throws ErrorException NeXLCore._mp_level5("SiO2_")
    @test_throws ErrorException NeXLCore._mp_level5("Si_O2")
    @test_throws ErrorException NeXLCore._mp_level5("9SiO2")

    @test NeXLCore._mp_level4("Ca5(PO4)3F") == Dict{Element, Int}(n"P" => 3, n"Ca" => 5, n"F" => 1, n"O" => 12)
    @test NeXLCore._mp_level4("Ca5(PO4)3OH") == Dict{Element, Int}(n"P" => 3, n"Ca" => 5, n"H" => 1, n"O" => 13)
    @test NeXLCore._mp_level4("(Ca)5(PO4)3OH") == Dict{Element, Int}(n"P" => 3, n"Ca" => 5, n"H" => 1, n"O" => 13)
    @test NeXLCore._mp_level4("(Ca5)(PO4)3(OH)")  == Dict{Element, Int}(n"P" => 3, n"Ca" => 5, n"H" => 1, n"O" => 13)
    @test NeXLCore._mp_level4("(Ca5)(PO4)3OH")  == Dict{Element, Int}(n"P" => 3, n"Ca" => 5, n"H" => 1, n"O" => 13)
    @test_throws ErrorException NeXLCore._mp_level4("(Ca5(PO4)3OH")
    @test_throws ErrorException NeXLCore._mp_level4("Ca5((PO4)3OH")
    @test_throws ErrorException NeXLCore._mp_level4("Ca5)(PO4)3OH")

    @test NeXLCore._mp_level3("Ca5(PO4)3⋅OH") == Dict{Element, Int}(n"P" => 3, n"Ca" => 5, n"H" => 1, n"O" => 13)
    @test NeXLCore._mp_level3("Ca5(PO4)3⋅1OH") == Dict{Element, Int}(n"P" => 3, n"Ca" => 5, n"H" => 1, n"O" => 13)

    function cda(d1, d2)
        all(isapprox(uv(get(d1,el, zero(valtype(d1)))), uv(get(d2,el,zero(valtype(d2)))), atol=1.0e-4) for el in intersect(keys(d1), keys(d2)))
    end

    # Check a couple of materials (Compare to DTSA-II)
    @test cda(NeXLCore.mmm_arsenopyrite.massfraction, Dict(n"S"=>0.1928, n"Fe"=>0.3416,n"Co"=>0.0002,n"As"=>0.4610))
    @test cda(NeXLCore.mmm_biotite.massfraction, Dict(n"H"=>0.0043,n"O"=>0.4372,n"Mg"=>0.1186,n"Al"=>0.0846,n"Si"=>0.1791,n"K"=>0.0835,n"Ti"=>0.0094,n"Fe"=>0.0780))
    @test cda(Dict(elm=>value(v) for (elm, v) in NeXLCore.srm470_k412.massfraction), Dict(n"O"=>0.42758,n"Mg"=>0.116567,n"Al"=>0.0490615,n"Si"=>0.211982,n"Ca"=>0.10899,n"Fe"=>0.077419))

    @test cda(NeXLCore._mp_level1("Ca5(PO4)3⋅1OH",lookup=s->nothing), Dict{Element, Float64}(n"P" => 0.184989, n"Ca" => 0.398942, n"H" => 0.00200674, n"O" => 0.414062))

    # Test the lookup function for named materials
    function luf(s)
        d = Dict(
            "K412"=>NeXLCore.srm470_k412.massfraction,
            "K411"=>NeXLCore.srm470_k411.massfraction,
            "Arsenopyrite"=>NeXLCore.mmm_arsenopyrite.massfraction
        )
        return get(d, s, nothing)
    end
    res1 = Dict(elm=>0.8*val for (elm,val) in srm470_k412.massfraction)
    res2 = Dict(elm=>0.2*val for (elm,val) in srm470_k411.massfraction)
    res = Dict(elm=>sum(uv(get(res1, elm, zero(valtype(res1)))), uv(get(res2, elm, zero(valtype(res2))))) for elm in union(keys(res1),keys(res2)))
    @test cda(parse(Material, "0.8*K412+0.2*K411", lookup=luf).massfraction, res)

    res1 = Dict(elm=>0.34*val for (elm,val) in NeXLCore.mmm_arsenopyrite.massfraction)
    res2 = Dict(elm=>0.62*val for (elm,val) in mat"NaAlSi3O8".massfraction)
    res = Dict(elm=>sum(uv(get(res1, elm, zero(valtype(res1)))), uv(get(res2, elm, zero(valtype(res2))))) for elm in union(keys(res1),keys(res2)))
    @test cda(parse(Material, "0.34*Arsenopyrite+0.62*NaAlSi3O8", lookup=luf).massfraction, res)

    @test isapprox(mat"RbTiO.PO4"[n"P"], 0.1268, atol=0.0001)
    @test isapprox(mat"Ca5(PO4)3⋅1OHNa"[n"H"], 0.0019, atol=0.0001)
    @test isapprox(mat"Ca5(PO4)3⋅1OHNa"[n"O"], 0.3959, atol=0.0001)

end


text:
@testset "Testing Allocators" begin
    verify_alignment(p) = @test iszero(mod(convert(UInt, p), 64))

    #####
    ##### Aligned Allocator
    #####

    allocator = CachedArrays.AlignedAllocator()

    for _ in 1:10000
        ptr, token = CachedArrays.allocate(allocator, rand(1:1_000_000))
        verify_alignment(ptr)
        CachedArrays.free(ptr, token)
    end
end


text:
function kernel!(F, v_m, v_a,
                                ybus_re_nzval, ybus_re_colptr, ybus_re_rowval,
                                ybus_im_nzval, ybus_im_colptr, ybus_im_rowval,
                                pinj, qinj, pv, pq, nbus, ::CPUBackend)
    npv = length(pv)
    npq = length(pq)
    n = npq + npv
    # i = @index(Global, Linear)
    for i in 1:n
        # REAL PV: 1:npv
        # REAL PQ: (npv+1:npv+npq)
        # IMAG PQ: (npv+npq+1:npv+2npq)
        fr = (i <= npv) ? pv[i] : pq[i - npv]
        F[i] -= pinj[fr]
        if i > npv
            F[i + npq] -= qinj[fr]
        end
        @inbounds for c in ybus_re_colptr[fr]:ybus_re_colptr[fr+1]-1
            to = ybus_re_rowval[c]
            aij = v_a[fr] - v_a[to]
            # f_re = a * cos + b * sin
            # f_im = a * sin - b * cos
            coef_cos = v_m[fr]*v_m[to]*ybus_re_nzval[c]
            coef_sin = v_m[fr]*v_m[to]*ybus_im_nzval[c]
            cos_val = cos(aij)
            sin_val = sin(aij)
            F[i] += coef_cos * cos_val + coef_sin * sin_val
            if i > npv
                F[npq + i] += coef_cos * sin_val - coef_sin * cos_val
            end
        end
    end
end

function kernel_roc!(F, v_m, v_a,
                                ybus_re_nzval, ybus_re_colptr, ybus_re_rowval,
                                ybus_im_nzval, ybus_im_colptr, ybus_im_rowval,
                                pinj, qinj, pv, pq, nbus)
    npv = length(pv)
    npq = length(pq)
    n = npq + npv
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    for i in index:stride:n
        # REAL PV: 1:npv
        # REAL PQ: (npv+1:npv+npq)
        # IMAG PQ: (npv+npq+1:npv+2npq)
        fr = (i <= npv) ? pv[i] : pq[i - npv]
        F[i] -= pinj[fr]
        if i > npv
            F[i + npq] -= qinj[fr]
        end
        for c in ybus_re_colptr[fr]:ybus_re_colptr[fr+1]-1
            to = ybus_re_rowval[c]
            aij = v_a[fr] - v_a[to]
            # f_re = a * cos + b * sin
            # f_im = a * sin - b * cos
            coef_cos = v_m[fr]*v_m[to]*ybus_re_nzval[c]
            coef_sin = v_m[fr]*v_m[to]*ybus_im_nzval[c]
            cos_val = AMDGPU.GPUArrays.cos(aij)
            sin_val = AMDGPU.GPUArrays.sin(aij)
            F[i] += coef_cos * cos_val + coef_sin * sin_val
            if i > npv
                F[npq + i] += coef_cos * sin_val - coef_sin * cos_val
            end
        end
    end
    return nothing
end

function kernel_oneapi!(F, v_m, v_a,
                                ybus_re_nzval, ybus_re_colptr, ybus_re_rowval,
                                ybus_im_nzval, ybus_im_colptr, ybus_im_rowval,
                                pinj, qinj, pv, pq, nbus)
    npv = length(pv)
    npq = length(pq)
    n = npq + npv
    i = get_global_id()
    # REAL PV: 1:npv
    # REAL PQ: (npv+1:npv+npq)
    # IMAG PQ: (npv+npq+1:npv+2npq)
    fr = (i <= npv) ? pv[i] : pq[i - npv]
    F[i] -= pinj[fr]
    if i > npv
        F[i + npq] -= qinj[fr]
    end
    @inbounds for c in ybus_re_colptr[fr]:ybus_re_colptr[fr+1]-1
        to = ybus_re_rowval[c]
        aij = v_a[fr] - v_a[to]
        # f_re = a * cos + b * sin
        # f_im = a * sin - b * cos
        coef_cos = v_m[fr]*v_m[to]*ybus_re_nzval[c]
        coef_sin = v_m[fr]*v_m[to]*ybus_im_nzval[c]
        # TODO: support cos and sin in oneAPI
        #cos_val = oneAPI.GPUArrays.cos(aij)
        #sin_val = oneAPI.GPUArrays.sin(aij)
        cos_val = aij
        sin_val = aij
        F[i] += coef_cos * cos_val + coef_sin * sin_val
        if i > npv
            F[npq + i] += coef_cos * sin_val - coef_sin * cos_val
        end
    end
    return nothing
end

function kernel_cuda!(F, v_m, v_a,
                                ybus_re_nzval, ybus_re_colptr, ybus_re_rowval,
                                ybus_im_nzval, ybus_im_colptr, ybus_im_rowval,
                                pinj, qinj, pv, pq, nbus)
    npv = length(pv)
    npq = length(pq)
    n = npq + npv
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    for i in index:stride:n
        # REAL PV: 1:npv
        # REAL PQ: (npv+1:npv+npq)
        # IMAG PQ: (npv+npq+1:npv+2npq)
        fr = (i <= npv) ? pv[i] : pq[i - npv]
        F[i] -= pinj[fr]
        if i > npv
            F[i + npq] -= qinj[fr]
        end
        @inbounds for c in ybus_re_colptr[fr]:ybus_re_colptr[fr+1]-1
            to = ybus_re_rowval[c]
            aij = v_a[fr] - v_a[to]
            # f_re = a * cos + b * sin
            # f_im = a * sin - b * cos
            coef_cos = v_m[fr]*v_m[to]*ybus_re_nzval[c]
            coef_sin = v_m[fr]*v_m[to]*ybus_im_nzval[c]
            cos_val = CUDA.cos(aij)
            sin_val = CUDA.sin(aij)
            F[i] += coef_cos * cos_val + coef_sin * sin_val
            if i > npv
                F[npq + i] += coef_cos * sin_val - coef_sin * cos_val
            end
        end
    end
    return nothing
end


text:
"""
    @exception(
        docstring::AbstractString,
        module_name::AbstractString,
        exception_name::AbstractString,
        error_message::AbstractString,
    )

Construct an exception.

# Usage

```julia
using Scats.Internal.Extras: @exception
@exception(
    "Exception thrown when the file is not found.",
    "Gen",
    "ScatsGenNotAFile",
    "The file is not found (\\\"\", e.file, \"\\\").",
)
```

This macro call is equivalent to this exception definition:

```jldoctest; output = false
\"\"\"
    ScatsGenNotAFile <: Exception

Exception thrown when the file is not found.

\"\"\"
mutable struct ScatsGenNotAFile <: Exception
    file::AbstractString
    ScatsGenNotAFile(file::AbstractString) = new(file)
end

Base.showerror(io::IO, e::ScatsGenNotAFile) =
print(
    io, string(
        "\\n\\n",
        "Scats.Internal.Gen.ScatsGenNotAFile:\\n",
        "The file is not found (\\\"\", e.file, \"\\\").\\n"
    )
)

# output

```
"""
macro exception(
    docstring::AbstractString,
    module_name::AbstractString,
    exception_name::AbstractString,
    error_message::AbstractString,
)
    return esc(
        Meta.parse("""
        begin
        \"\"\"
            $exception_name <: Exception

        $docstring

        \"\"\"
        mutable struct $exception_name <: Exception
            file::AbstractString
            $exception_name(file::AbstractString) = new(file)
        end

        Base.showerror(io::IO, e::$exception_name) =
        print(
            io, string(
                "\n\n",
                "Scats.Internal.$module_name.$exception_name:\n",
                "$error_message\n"
            )
        )

        end""")
    )
end

"""
    @unpack(
        struct_instance::AbstractString,
        fields::AbstractString...,
    )

Unpack an instance of a structure on the specified fields.

# Usage

```jldoctest; output = false
using Scats.Internal.Extras: @unpack
using Scats.Internal.Gen: GenStruct
Gen = GenStruct()

@unpack("Gen", "N", "Δt", "q", "α", "β", "r", "A", "ν", "ϕ", "γ") ===
begin
N = Gen.N
Δt = Gen.Δt
q = Gen.q
α = Gen.α
β = Gen.β
r = Gen.r
A = Gen.A
ν = Gen.ν
ϕ = Gen.ϕ
γ = Gen.γ
end

# output

true
```
"""
macro unpack(
    struct_instance::AbstractString,
    fields::AbstractString...,
)
    result_string = string(
        "begin\n",
        join(
            [ "$field = $struct_instance.$field\n" for field in fields ]
        ),
        "end"
    )
    return esc(Meta.parse(result_string))
end

"""
    @pack(
        struct_instance::AbstractString,
        fields::AbstractString...,
    )

Pack the specified variables in the structure instance fields of the same name.

# Usage

```jldoctest; output = false
using Scats.Internal.Extras: @pack
using Scats.Internal.Input: InputStruct
Input = InputStruct()

N = 230
Δt = 1
q = 0.01

@pack("Input", "N", "Δt", "q") ===
begin
Input.N = N
Input.Δt = Δt
Input.q = q
end

# output

true
```
"""
macro pack(
    struct_instance::AbstractString,
    fields::AbstractString...,
)
    result_string = string(
        "begin\n",
        join(
            [ "$struct_instance.$field = $field\n" for field in fields ]
        ),
        "end"
    )
    return esc(Meta.parse(result_string))
end


text:
"""
mutable struct InverseDistanceSurrogate{X,Y,P,L,U} <: AbstractSurrogate

The square polynomial model can be expressed by 𝐲 = 𝐗β + ϵ, with β = 𝐗ᵗ𝐗⁻¹𝐗ᵗ𝐲
"""
mutable struct SecondOrderPolynomialSurrogate{X,Y,B,L,U} <: AbstractSurrogate
    x::X
    y::Y
    β::B
    lb::L
    ub::U
end

function SecondOrderPolynomialSurrogate(x,y,lb::Number,ub::Number)
    n = length(x)
    d = 1
    X = ones(eltype(x[1]),n,3*d)
    X[:,2] = x
    X[:,3] = x.^2
    β = X\y
    return SecondOrderPolynomialSurrogate(x,y,β,lb,ub)
end

function (sec_ord::SecondOrderPolynomialSurrogate)(val::Number)
    return sec_ord.β[1] + sec_ord.β[2]*val + sec_ord.β[3]*val^2
end

function SecondOrderPolynomialSurrogate(x,y,lb,ub)
    n = length(x)
    d = length(lb)
    X = ones(eltype(x[1]),n,3*d)
    for j = 1:d
        X[:,j+1] =[x[i][j] for i=1:n]
    end
    for j = 1:d-1
        X[:,j+d+1] = [x[i][j]*x[i][j+1] for i = 1:n]
    end
    for j = 1:d
        X[:,j+2*d] = [x[i][j]^2 for i=1:n]
    end
    β = X\y
    return SecondOrderPolynomialSurrogate(x,y,β,lb,ub)
end

function (my_second_ord::SecondOrderPolynomialSurrogate)(val)
    #just create the val vector as X and multiply
    d = length(val)
    X = [[one(eltype(val[1]))]; [val[j] for j =1:d]; [val[j]*val[j+1] for j = 1:d-1]; [val[j]^2 for j = 1:d]]
    return my_second_ord.β'*X
end

function add_point!(my_second::SecondOrderPolynomialSurrogate,x_new,y_new)
    d = length(my_second.lb)
    if d == 1
        #1D
        my_second.x = vcat(my_second.x,x_new)
        my_second.y = vcat(my_second.y,y_new)
        n = length(my_second.x)
        X = ones(eltype(my_second.x[1]),n,3*d)
        X[:,2] = my_second.x
        X[:,3] = my_second.x.^2
        my_second.β = X\my_second.y
    else
        #ND
        my_second.x = vcat(my_second.x,x_new)
        my_second.y = vcat(my_second.y,y_new)
        n = length(my_second.x)
        X = ones(eltype(my_second.x[1]),n,3*d)
        for j = 1:d
            X[:,j+1] =[my_second.x[i][j] for i=1:n]
        end
        for j = 1:d-1
            X[:,j+d+1] = [my_second.x[i][j]*my_second.x[i][j+1] for i = 1:n]
        end
        for j = 1:d
            X[:,j+2*d] = [my_second.x[i][j]^2 for i=1:n]
        end
        my_second.β = X\my_second.y
    end
    nothing
end


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using GLMakie
using GeometryBasics, LinearAlgebra, StatsBase
using Makie: get_dim, surface_normals
GLMakie.activate!() # HIDE
let
    ### The bottle
    u = LinRange(0,π,100)
    v = LinRange(0,2π,100)
    x = [-2/15 * cos(u)*(3*cos(v) - 30*sin(u) + 90*cos(u)^4 *sin(u)
            - 60*cos(u)^6 * sin(u) + 5*cos(u)*cos(v)*sin(u)) for u in u, v in v]
    y = [-1/15 * sin(u)*(3*cos(v) - 3*cos(u)^2*cos(v) - 48*cos(u)^4 *cos(v) + 48*cos(u)^6 *cos(v)
            - 60*sin(u) + 5*cos(u)*cos(v)*sin(u) - 5*cos(u)^3*cos(v)*sin(u)
            -80*cos(u)^5 * cos(v)*sin(u) + 80*cos(u)^7*cos(v)*sin(u)) for u in u, v in v]
    z = [2/15 * (3 + 5*cos(u)*sin(u))*sin(v) for u in u, v in v]

    function getMesh(x,y,z)
        positions = vec(map(CartesianIndices(z)) do i
        GeometryBasics.Point{3, Float32}(
            get_dim(x, i, 1, size(z)),
            get_dim(y, i, 2, size(z)),
            z[i])
        end)
        faces = decompose(GLTriangleFace, Rect2D(0f0, 0f0, 1f0, 1f0), size(z))
        normals = surface_normals(x, y, z)
        vertices = GeometryBasics.meta(positions; normals=normals)
        meshObj = GeometryBasics.Mesh(vertices, faces)
        meshObj
    end
    meshKlein = getMesh(x,y,z)
    fig = Figure(resolution = (1200, 800), backgroundcolor = :white)
    fig.layout.alignmode[] = Outside(0)
    ax = LScene(fig, scenekw = (camera = cam3d!, show_axis = false))
    mesh!(ax,meshKlein, color = [norm(v) for v in coordinates(meshKlein)],
        colormap = (:Spectral_11,0.95))
    rotate_cam!(ax.scene, (-1.5,0,0))
    fig[1,1] = ax
    fig
    save(joinpath(@__DIR__, "output", "KleinBottle.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["GLMakie", "GeometryBasics", "LinearAlgebra", "StatsBase", "Makie"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using RootsNewton
using Test

#test to see if returns nothing for functions without real roots
f(x)=x^2+2
fp(x)=2*x
@test newtonroot(f,fp;iv=1.0, tolerance=1E-7, maxiter=1000).value==nothing

#Tests of functions with known solutions
f(x)=x-2
fp(x)=1
@test newtonroot(f,fp;iv=1.0, tolerance=1E-7, maxiter=1000).value ≈ 2

f(x)=x^3-1
fp(x)=3*x^2
@test newtonroot(f,fp;iv=1.0, tolerance=1E-7, maxiter=1000).value ≈ 1

f(x)=log(x)
fp(x)=1/x
@test newtonroot(f,fp;iv=1.0, tolerance=1E-7, maxiter=1000).value ≈ 1

f(x)=cos(x)
fp(x)=-sin(x)
@test newtonroot(f,fp;iv=1.0, tolerance=1E-7, maxiter=1000).value ≈ pi/2

#Tests using automatic differentiation

f(x)=x-2
@test newtonroot(f;iv=1.0, tolerance=1E-7, maxiter=1000).value ≈ 2

f(x)=x^3-1
@test newtonroot(f;iv=2.0, tolerance=1E-7, maxiter=1000).value ≈ 1

f(x)=log(x)
@test newtonroot(f;iv=2.0, tolerance=1E-7, maxiter=1000).value ≈ 1

f(x)=cos(x)
@test newtonroot(f;iv=1.0, tolerance=1E-7, maxiter=1000).value ≈ pi/2

#Tests with BigFloat
f(x)=x-2
@test newtonroot(f;iv=BigFloat(1.0), tolerance=1E-7, maxiter=1000).value ≈ 2

#Tests that maxiter is working
f(x)=x^3-1
@test newtonroot(f;iv=100, tolerance=1E-7, maxiter=3).iter <= 3

#Tests that tolerance is working
f(x)=x^3-1
@test 1-2 <= newtonroot(f;iv=10, tolerance=2, maxiter=1000).value <= 1+2 || newtonroot(f;iv=10, tolerance=2, maxiter=1000).iter==1000

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


text:
struct Operator{T <: AbstractMatrix{Scalar}} <: AbstractQuantumMatrix
    data::T
    function Operator{T}(A::AbstractMatrix)
        new{T}(A)
    end
end

const Hamiltonian = Operator{Hermitian}
const Potential = Operator{Diagonal}
const KineticEnergy = Operator{Tridiagonal}

Operator(A::T) where {T <: AbstractMatrix{Scalar}} = Operator{T}(A)
Operator(N::Integer) = Operator(zeros(Scalar, N, N))

function (A::Operator)(ψ::Ket)
    check_dimensions_match(A,ψ)
    Ket(data(A) * data(ψ))
end

function (A::Operator)(ϕ::Bra)
    check_dimensions_match(A, ϕ)
    Bra(adjoint(data(A)) * data(ϕ))
end

function *(A::Operator, B::Operator)
    check_dimensions_match(A, B)
    Operator(data(A) * data(B))
end

*(A::Operator, ψ::Ket) = (A)(ψ)
*(ϕ::Bra, A::Operator) = (A)(ϕ)

function mul!(C::AbstractQuantumMatrix, A::AbstractQuantumMatrix, B::AbstractQuantumMatrix, α::Scalar, β::Scalar)
    mul!(data(C), data(A), data(B), α, β)
end


text:
using GeneralizedGenerated
using JuliaVariables
using Test
using BenchmarkTools
using DataStructures

rmlines = NGG.rmlines

@testset "no kwargs" begin

@gg function f1(a)
    quote
        x -> a + x
    end |> rmlines
end

@test f1(1)(2) == 3

@gg function f2(a)
    quote
        x -> begin
            a += 2
            x + a
        end
    end |> rmlines
end

@test f2(1)(2) == 5

@gg function f3(a)
    quote
        k = 20
        x -> begin
            a += 2
            x + a + k
        end
    end
end

@test f3(1)(2) == 25

end

@testset "kwargs" begin

@gg function f4(a)
    quote
        function z(x, k=1)
            x + 20 + a + k
        end
    end
end

@test f4(10)(2) == 33
end

@testset "namedtuple" begin

@gg function f5(a)
    quote
        function z(x, k=1)
            (x=x, k=k, q=a)
        end
    end
end

@test f5(10)(2) == (x=2, k=1, q=10)
end

@testset "mk funcs" begin

f_ = mk_function(:((x, y) -> x + y))
@test f_(1, 2) == 3

f_ = mk_function(:(function (x, y) x + y end))
@test f_(1, 2) == 3

end

@testset "type encoding more datatypes" begin

@gg function f5(a)
    tp = (1, 2, 3)
    quote
        function z(x, k=$tp)
            (x=x, k=k, q=a)
        end
    end
end

@test f5(10)(2) == (x=2, k=(1, 2, 3), q=10)

@gg function f5(a)
    tp = (a1=1, a2=2, a3=3)
    quote
        function z(x, k=$tp)
            (x=x, k=k, q=a)
        end
    end
end

@test f5(10)(2) == (x=2, k=(a1=1, a2=2, a3=3), q=10)

@gg function f5(a)
    tp = "233"
    quote
        function z(x, k=$tp)
            (x=x, k=k, q=a)
        end
    end
end

@test f5(10)(2) == (x=2, k="233", q=10)

@gg function f5(a)
    tp = list(1, 2, 3)
    quote
        function z(x; k=$tp)
            (x=x, k=k, q=a)
        end
    end
end

@test f5(10)(2) == (x=2, k=list(1, 2, 3), q=10)
@test f5(10)(2; k=10) == (x=2, k=10, q=10)

end

@testset "runtime eval" begin

a = to_type(:(1 + 2))
@test :(1 + 2) == from_type(a)
@test string(from_type(a)) == string(:(1 + 2))

@test runtime_eval(1) == 1
@test mk_function(:(
    x -> x + 1
))(2) == 3

@test_throws Any mk_function(quote
    x -> x + 1
end)

@test runtime_eval(quote
    x -> x + 1
end)(1) == 2

end

@testset "self recursive" begin
    to_test = quote
        g(x, r=0) = x === 0 ? r : begin
            g = g # required for self recur
            g(x-1, r + x)
        end
        g(10)
    end |> runtime_eval

    g(x, r=0) = x === 0 ? r : g(x-1, r + x)
    expected = g(10)
    @test expected == to_test
end

@testset "self recursive" begin
    to_test = quote
        g(x, r=0) = x === 0 ? r : begin
            g = g # required for self recur
            g(x-1, r + x)
        end
        g(10)
    end |> runtime_eval

    g(x, r=0) = x === 0 ? r : g(x-1, r + x)
    expected = g(10)
    @test expected == to_test
end

@testset "support where clauses and return type annotations for @gg" begin
    @gg function foo(x::T) where T
        :(x, T)
    end
    @test foo(1) == (1, Int)
    @gg function bar(x::T) where T
        quote
            g = x + 20
            x = 10
            () -> begin
                x = g
                x
            end
        end
    end
    @test bar(2)() == 2 + 20

    @gg function foobar(x::T, y::A) where {T <: Number, A <: AbstractArray{T}}
        quote
            g = x + 20
            x = 10
            () -> begin
                x = g
                (A, x + y[1])
            end
        end
    end
    @test foobar(2, [3])() == (Vector{Int}, 2 + 20 + 3)
end

@testset "support default arguments" begin
    @gg function h(x, c)
        quote
            d = x + 10
            function g(x, y=c)
                x + y + d
            end
        end
    end
    @test h(1, 2)(3) == 16
end

module S
    run(y) = y + 1
end

struct K
    f1::Function
    f2::Function
end
@testset "specifying evaluation modules" begin
    @gg m function g(m::Module, y) :(run(y)) end
    @test g(S, 1) == 2

    @gg m function h(m, y)
        quote
        c = m.f1(y)
        () -> begin c = m.f2(c) end
        end
    end
    k = K(x -> x + 1, x -> x * 9)
    next = h(k, 1)
    @test next() == 18
    @test next() == 18 * 9
end

@testset "test free variables of let bindings" begin
    @gg function test_free_of_let()
        quote
            let x = 1
                f = () -> begin
                    x * 3
                end
                x = 2
                f
            end
        end
    end
    @test test_free_of_let()() == 6
end

@testset "show something" begin
    f1 = mk_function(:(x -> x + 1))
    f2 = mk_function(:((x :: Int = 2, ) -> x + 1))
    @test f1(1) == 2
    @test f2() == 3
    println(f1)
    println(f2)
end

@testset "omit func argname: #34" begin
   f1 = mk_function(:( (:: Int) -> 0 ))
   @test f1(1) == 0
   @test_throws MethodError f1("")
end

@testset "support for 1-d arrays' type encoding: #42" begin
    from_type(to_type([1, 2, 3])) == [1, 2, 3]
    from_type(to_type([1, 2, "3"])) == [1, 2, "3"]
end


text:
using PhaseRec
using XUnit
using LinearAlgebra

@testset "Test PhaseRec" begin include("phaserec.jl") end


text:
#ДАНО: Робот - в произвольной клетке ограниченного прямоугольного поля (без внутренних перегородок)
#РЕЗУЛЬТАТ: Робот - в исходном положении,
#в клетке с роботом стоит маркер, и все остальные клетки поля промаркированы в шахматном порядкe

function moves!(r,side) #двигает робота и возвращает количество шагов
    num = 0
    while isborder(r,side) == false
        move!(r,side)
        num+=1
    end
    return num
end

function reverse(side)
    for i=0:3
        if side == HorizonSide(i)
            return HorizonSide(mod(i+2, 4))
        end
    end
end

function move_back!(r,side,num) #передвигает робота в сторону на число шагов
    for i = 1:num
        move!(r,side)
    end
end

function mark_chek(r, x_coord, y_coord, size) #проверка координаты на наличие маркера(size-размер отступа)
    if ((mod(x_coord,2*size) in 0:size-1)&&(mod(y_coord,2*size) in 0:size-1))
        putmarker!(r)
    end
    if ((mod(x_coord,2*size) in size:2*size-1)&&(mod(y_coord,2*size) in size:2*size-1))
        putmarker!(r)
    end
end

function chestt(r)
    side = Ost
    y = moves!(r,West)
    x = moves!(r,Sud)
    num = 1
    while isborder(r,Nord) == false
        mark_chek(r,abs(x),abs(y),1)
        while isborder(r,side) == false
            move!(r,side)
            y -= num
            mark_chek(r,abs(x),abs(y),1)
        end
        move!(r,Nord)
        side = inverse(side)
        x -= 1
        num *= -1
    end
    mark_chek(r,abs(x),abs(y),1)
    while isborder(r,side)==false
        move!(r,side)
        y -= num
        mark_chek(r,abs(x),abs(y),1)
    end
    move_back(r,Nord,x)
    move_back(r,Ost,y)
end

max_stars_repo_name: gardend/mirea-progs

max_issues_repo_name: gardend/mirea-progs

max_forks_repo_name: gardend/mirea-progs


text:
# Autogenerated wrapper script for TetGen_jll for x86_64-w64-mingw32
export libtet

JLLWrappers.@generate_wrapper_header("TetGen")
JLLWrappers.@declare_library_product(libtet, "libtet.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libtet,
        "bin\\libtet.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Test, MLToolkit.DistributionsX

@testset "DistributionsX" begin
    tests = [
        "noise.jl",
        "normal.jl",
        "gumbel.jl",
        "bernoulli.jl",
        # "displaced_poisson.jl",
        # "ibp.jl",
        # "beta.jl",
        # "npd.jl",
    ]
    foreach(include, tests)

    @warn "`DistributionsX.jl` is not tested."
end


text:
using DrWatson
@quickactivate "QMC"

using QMC

# main.jl
#
# A projector QMC program for the TFIM

using Random
Random.seed!(1234)

using ProgressMeter

using Measurements
using Statistics
using FFTW

using DelimitedFiles
using JLD2
using Printf
using Lattices

using DataStructures
using ArgParse

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

function init_mc_cli(parsed_args)
    PBC = parsed_args["periodic"]
    h = parsed_args["field"]
    J = parsed_args["interaction"]

    Dim = length(parsed_args["dims"])
    nX = parsed_args["dims"]

    BC_name = PBC ? "PBC" : "OBC"

    if Dim == 1
        nX = nX[1]
        bond_spin, Ns, Nb = lattice_bond_spins(nX, PBC)
    elseif Dim == 2
        bond_spin, Ns, Nb = lattice_bond_spins(nX, PBC)
        nX = nX[1]
    else
        error("Unsupported number of dimensions")
    end

    # MC parameters
    M = parsed_args["M"] # length of the operator_list is 2M
    MCS = parsed_args["measurements"] # the number of samples to record
    EQ_MCS = div(MCS, 10)
    skip = parsed_args["skip"]  # number of MC steps to perform between each msmt

    # path = "$(Dim)D/$(nX)/$(BC_name)/J$(J)/h$(h)/skip$(skip)/"
    d = @ntuple Dim nX BC_name J h skip M
    mc_opts = @ntuple M MCS EQ_MCS skip

    # Ns, Nb = length(lattice), length(bond_spin)
    H = TFIM(bond_spin, Dim, Ns, Nb, h, J)
    if haskey(parsed_args, "beta")
        qmc_state = BinaryThermalState(H, 2M)
    else
        qmc_state = BinaryGroundState(H, M)
    end

    return H, qmc_state, savename(d; digits = 4), mc_opts
end

function make_info_file(info_file, samples_file, mc_opts, op_list_length, observables, corr_time)
    M, MCS, EQ_MCS, skip = mc_opts
    mag, abs_mag, mag_sqr, energy = observables

    open(info_file, "w") do file_io
        streams = [Base.stdout, file_io]

        for io in streams
            @printf(io, "⟨M⟩   = % .16f +/- %.16f\n", mag.val, mag.err)
            @printf(io, "⟨|M|⟩ = % .16f +/- %.16f\n", abs_mag.val, abs_mag.err)
            @printf(io, "⟨M^2⟩ = % .16f +/- %.16f\n", mag_sqr.val, mag_sqr.err)
            @printf(io, "⟨E⟩   = % .16f +/- %.16f\n\n", energy.val, energy.err)

            println(io, "Correlation time: $(corr_time)\n")

            println(io, "Initial Operator list length: $(2*M)")
            println(io, "Final Operator list length: $(op_list_length)")
            println(io, "Number of MC measurements: $(MCS)")
            println(io, "Number of equilibration steps: $(EQ_MCS)")
            println(io, "Number of skips between measurements: $(skip)\n")

            println(io, "Samples outputted to file: $(samples_file)")
        end
    end
end

function save_data(path, mc_opts, qmc_state, measurements, observables, corr_time)
    info_file = path * "_info.txt"
    samples_file = path * "_samples.txt"
    qmc_state_file = path * "_state.jld2"

    open(samples_file, "w") do io
        writedlm(samples_file, measurements, " ")
    end

    @time @save qmc_state_file qmc_state

    M = length(qmc_state.operator_list)

    make_info_file(info_file, samples_file, mc_opts, M, observables, corr_time)
end

function mixedstate(parsed_args)
    H, qmc_state, sname, mc_opts = init_mc_cli(parsed_args)
    beta = parsed_args["beta"]

    M, MCS, EQ_MCS, skip = mc_opts

    mkpath(datadir("sims", "mixedstate"))
    path = datadir("sims", "mixedstate", savename((@ntuple beta), sname; digits = 4))

    measurements = zeros(Int, MCS, nspins(H))
    mags = zeros(MCS)
    ns = zeros(MCS)

    max_ns = maximum(@showprogress "Warm up..." [mc_step_beta!(qmc_state, H, beta; eq=true) for i in 1:EQ_MCS])

    resize_op_list!(qmc_state, round(Int, (3//2)*max_ns, RoundUp))

    @showprogress "MCMC...   " for i in 1:MCS # Monte Carlo Steps
        ns[i] = mc_step_beta!(qmc_state, H, beta) do cluster_data, qmc_state, H
            spin_prop = qmc_state.left_config
            measurements[i, :] = spin_prop
            mags[i] = magnetization(spin_prop)
        end

        for _ in 1:skip
            mc_step_beta!(qmc_state, H, beta)
        end
    end

    mag = mean_and_stderr(mags)
    abs_mag = mean_and_stderr(abs, mags)
    mag_sqr = mean_and_stderr(abs2, mags)

    energy = mean_and_stderr(x -> -x/beta, ns) + H.J*nbonds(H) + H.h*nspins(H)
    energy /= nspins(H)

    observables = (mag, abs_mag, mag_sqr, energy)

    # measure correlation time from equilibriation samples
    @time corr_time = correlation_time(mags .^ 2)

    save_data(path, mc_opts, qmc_state, measurements, observables, corr_time)
end

function groundstate(parsed_args)
    H, qmc_state, sname, mc_opts = init_mc_cli(parsed_args)

    M, MCS, EQ_MCS, skip = mc_opts

    mkpath(datadir("sims", "groundstate"))
    path = datadir("sims", "groundstate", sname)

    measurements = zeros(Int, MCS, nspins(H))
    mags = zeros(MCS)
    ns = zeros(MCS)

    @showprogress "Warm up..." for i in 1:EQ_MCS
        mc_step!(qmc_state, H)
    end

    @showprogress "MCMC...   " for i in 1:MCS # Monte Carlo Production Steps
        mc_step!(qmc_state, H) do cluster_data, qmc_state, H
            spin_prop = sample(qmc_state)
            measurements[i, :] = spin_prop

            ns[i] = num_single_site_diag(qmc_state.operator_list)
            mags[i] = magnetization(spin_prop)
        end

        for _ in 1:skip
            mc_step!(qmc_state, H)
        end
    end

    mag = mean_and_stderr(mags)
    abs_mag = mean_and_stderr(abs, mags)
    mag_sqr = mean_and_stderr(abs2, mags)

    @time energy = jackknife(ns) do n
        if H.h != 0
            (-H.h * ((1.0 / n) - 1)) + H.J * (nbonds(H) / nspins(H))
        else
            H.J * (nbonds(H) / nspins(H))
        end
    end

    observables = (mag, abs_mag, mag_sqr, energy)

    # measure correlation time from equilibriation samples
    @time corr_time = correlation_time(mags .^ 2)

    save_data(path, mc_opts, qmc_state, measurements, observables, corr_time)
end

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

s = ArgParseSettings()

@add_arg_table! s begin
    "groundstate"
        help = "Use Projector SSE to simulate the ground state"
        action = :command
    "mixedstate"
        help = "Use vanilla SSE to simulate the system at non-zero temperature"
        action = :command
end

@add_arg_table! s["groundstate"] begin
    "dims"
        help = "The dimensions of the square lattice"
        required = true
        arg_type = Int
        nargs = '+'
    "--periodic", "-p"
        help = "Periodic BCs"
        action = :store_true
    "--field"
        help = "Strength of the transverse field"
        arg_type = Float64
        default = 1.0
    "--interaction", "-J"
        help = "Strength of the interaction"
        arg_type = Float64
        default = 1.0

    "-M"
        help = "Half-size of the operator list"
        arg_type = Int
        default = 1000

    "--measurements", "-n"
        help = "Number of samples to record"
        arg_type = Int
        default = 100_000
    "--skip", "-s"
        help = "Number of MC steps to perform between each measurement"
        arg_type = Int
        default = 0
end

import_settings!(s["mixedstate"], s["groundstate"])

@add_arg_table! s["mixedstate"] begin
    "--beta"
        help = "The inverse-temperature parameter for the simulation"
        arg_type = Float64
        default = 10.0
end

parsed_args = parse_args(ARGS, s)

if parsed_args["%COMMAND%"] == "groundstate"
    groundstate(parsed_args["groundstate"])
else
    mixedstate(parsed_args["mixedstate"])
end


text:
crmatrices=Dict(
(l=1, m=-1, mu=-1, i=1) => SVector{3,ComplexF64}(0.1666666666666667+0.0im, 0.0+0.1666666666666667im, 0.0+0.0im),
(l=1, m=-1, mu=0, i=1) => SVector{3,ComplexF64}(0.0+0.0im, 0.0+0.0im, 0.0+0.0im),
(l=1, m=-1, mu=1, i=1) => SVector{3,ComplexF64}(-0.1666666666666667+0.0im, 0.0+-0.1666666666666667im, 0.0+0.0im),
(l=1, m=0, mu=-1, i=1) => SVector{3,ComplexF64}(0.0+0.0im, 0.0+0.0im, -0.2357022603955158+0.0im),
(l=1, m=0, mu=0, i=1) => SVector{3,ComplexF64}(0.0+0.0im, 0.0+0.0im, 0.0+0.0im),
(l=1, m=0, mu=1, i=1) => SVector{3,ComplexF64}(0.0+0.0im, 0.0+0.0im, 0.2357022603955158+0.0im),
(l=1, m=1, mu=-1, i=1) => SVector{3,ComplexF64}(-0.1666666666666667+0.0im, 0.0+0.1666666666666667im, 0.0+0.0im),
(l=1, m=1, mu=0, i=1) => SVector{3,ComplexF64}(0.0+0.0im, 0.0+0.0im, 0.0+0.0im),
(l=1, m=1, mu=1, i=1) => SVector{3,ComplexF64}(0.1666666666666667+0.0im, 0.0+-0.1666666666666667im, 0.0+0.0im),
(l=1, m=-1, mu=-1, i=2) => SVector{3,ComplexF64}(0.0+-0.1666666666666667im, 0.1666666666666667+0.0im, 0.0+0.0im),
(l=1, m=-1, mu=0, i=2) => SVector{3,ComplexF64}(0.0+0.0im, 0.0+0.0im, 0.0+0.0im),
(l=1, m=-1, mu=1, i=2) => SVector{3,ComplexF64}(0.0+-0.1666666666666667im, 0.1666666666666667+0.0im, 0.0+0.0im),
(l=1, m=0, mu=-1, i=2) => SVector{3,ComplexF64}(0.0+0.0im, 0.0+0.0im, 0.0+0.2357022603955158im),
(l=1, m=0, mu=0, i=2) => SVector{3,ComplexF64}(0.0+0.0im, 0.0+0.0im, 0.0+0.0im),
(l=1, m=0, mu=1, i=2) => SVector{3,ComplexF64}(0.0+0.0im, 0.0+0.0im, 0.0+0.2357022603955158im),
(l=1, m=1, mu=-1, i=2) => SVector{3,ComplexF64}(0.0+0.1666666666666667im, 0.1666666666666667+0.0im, 0.0+0.0im),
(l=1, m=1, mu=0, i=2) => SVector{3,ComplexF64}(0.0+0.0im, 0.0+0.0im, 0.0+0.0im),
(l=1, m=1, mu=1, i=2) => SVector{3,ComplexF64}(0.0+0.1666666666666667im, 0.1666666666666667+0.0im, 0.0+0.0im),
(l=1, m=-1, mu=-1, i=3) => SVector{3,ComplexF64}(0.0+0.0im, 0.0+0.0im, 0.0+0.0im),
(l=1, m=-1, mu=0, i=3) => SVector{3,ComplexF64}(-0.2357022603955158+0.0im, 0.0+-0.2357022603955158im, 0.0+0.0im),
(l=1, m=-1, mu=1, i=3) => SVector{3,ComplexF64}(0.0+0.0im, 0.0+0.0im, 0.0+0.0im),
(l=1, m=0, mu=-1, i=3) => SVector{3,ComplexF64}(0.0+0.0im, 0.0+0.0im, 0.0+0.0im),
(l=1, m=0, mu=0, i=3) => SVector{3,ComplexF64}(0.0+0.0im, 0.0+0.0im, 0.3333333333333333+0.0im),
(l=1, m=0, mu=1, i=3) => SVector{3,ComplexF64}(0.0+0.0im, 0.0+0.0im, 0.0+0.0im),
(l=1, m=1, mu=-1, i=3) => SVector{3,ComplexF64}(0.0+0.0im, 0.0+0.0im, 0.0+0.0im),
(l=1, m=1, mu=0, i=3) => SVector{3,ComplexF64}(0.2357022603955158+0.0im, 0.0+-0.2357022603955158im, 0.0+0.0im),
(l=1, m=1, mu=1, i=3) => SVector{3,ComplexF64}(0.0+0.0im, 0.0+0.0im, 0.0+0.0im)
)


text:
using DataFrames
using CSV

include("btgDerivatives.jl")
include("integration.jl")
using Polynomials

#f = x -> sin(x^2)
#df = x -> 2*x*cos(x^2)
#df2 = x -> 2*cos(x^2)-4*x^2*sin(x^2)

#f = x -> x^14
#df = x-> 14 * x^13
#df2 = x -> 182 * x^12

f =  x -> x .^47 + x.^24
df = x -> 47 * x .^46 + 24*x.^23
df2 = x -> 2162*x .^45 + 552*x.^22

#f = x -> sin.(exp.(3 .-x))
#df = x -> -exp.(3 .-x) .* cos.(exp.(3 .-x))
#df2 = x -> exp.(3 .-2*x) .*(exp.(x) .* cos.(exp.(3 .-x)).-exp(3) .* sin.(exp.(3 .-x)))

#f = x-> sin.(x .+pi/2)
#df = x-> cos.(x .+pi/2)
#df2 = x -> -sin.(x .+pi/2)

(nodes, weights) = getTuranData()

(h, A) = checkDerivative(df, df2, .5)
println(polyfit(h, A, 1))

function eval_Gauss_Turan(f, df, df2, nodes, weights)
    fn = f.(nodes)
    dfn = df.(nodes)
    df2n = df2.(nodes)
    return fn'*weights[:, 1] + dfn'*weights[:, 2] + df2n'*weights[:, 3]
end

res1 = int1D(f, [-1, 1.1])
println("Gauss error: ", abs(res1 - 2.47375))

F = x -> [f(x); df(x); df2(x)]
res = Gauss_Turan(F, [-1, 1.1], nodes, weights)
println("Turan error: ", abs(res - 2.47375))

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
module Generate

import ..Git, ..Read

copyright_year() = readchomp(`date +%Y`)
copyright_name() = readchomp(`git config --global --get user.name`)
github_user() = readchomp(ignorestatus(`git config --global --get github.user`))

function git_contributors(dir::String, n::Int=typemax(Int))
    contrib = Dict()
    tty = @windows? "CON:" : "/dev/tty"
    for line in eachline(tty |> Git.cmd(`shortlog -nes`, dir=dir))
        m = match(r"\s*(\d+)\s+(.+?)\s+\<(.+?)\>\s*$", line)
        m == nothing && continue
        commits, name, email = m.captures
        if haskey(contrib,email)
            contrib[email][1] += int(commits)
        else
            contrib[email] = [int(commits),name]
        end
    end
    names = Dict()
    for (commits,name) in values(contrib)
        names[name] = get(names,name,0) + commits
    end
    names = sort!(collect(keys(names)),by=name->names[name],rev=true)
    length(names) <= n ? names : [names[1:n], "et al."]
end

function package(
    pkg::String,
    license::String;
    force::Bool = false,
    authors::String = "",
    years::Union(Int,String) = copyright_year(),
    user::String = github_user(),
)
    isnew = !ispath(pkg)
    try
        if isnew
            url = isempty(user) ? "" : "git://github.com/$user/$pkg.jl.git"
            Generate.init(pkg,url)
        else
            Git.dirty(dir=pkg) && error("$pkg is dirty – commit or stash your changes")
        end
        Git.transact(dir=pkg) do
            if isempty(authors)
                authors = isnew ? copyright_name() : git_contributors(pkg,5)
            end
            Generate.license(pkg,license,years,authors,force=force)
            Generate.readme(pkg,user,force=force)
            Generate.entrypoint(pkg,force=force)
            Generate.tests(pkg,force=force)
            Generate.travis(pkg,force=force)

            msg = """
            $pkg.jl $(isnew ? "generated" : "regenerated") files.

                license:  $license
                authors:  $(join([authors],", "))
                years:    $years
                user:     $user

            Julia Version $VERSION [$(Base.GIT_VERSION_INFO.commit_short)]
            """

            if isnew
                info("Committing $pkg generated files")
                Git.run(`commit -q -m $msg`, dir=pkg)
            elseif Git.dirty(dir=pkg)
                Git.run(`reset -q --`, dir=pkg)
                info("Regenerated files left unstaged, use `git add -p` to select")
                open(io->print(io,msg), joinpath(Git.dir(pkg),"MERGE_MSG"), "w")
            else
                info("Regenerated files are unchanged")
            end
        end
    catch
        isnew && rm(pkg, recursive=true)
        rethrow()
    end
end

function init(pkg::String, url::String="")
    if !ispath(pkg)
        info("Initializing $pkg repo: $(abspath(pkg))")
        Git.run(`init -q $pkg`)
        Git.run(`commit -q --allow-empty -m "initial empty commit"`, dir=pkg)
    end
    isempty(url) && return
    info("Origin: $url")
    Git.run(`remote add origin $url`,dir=pkg)
    Git.set_remote_url(url,dir=pkg)
end

function license(pkg::String, license::String,
                 years::Union(Int,String),
                 authors::Union(String,Array);
                 force::Bool=false)
    genfile(pkg,"LICENSE.md",force) do io
        if !haskey(LICENSES,license)
            licenses = join(sort!([keys(LICENSES)...], by=lowercase), ", ")
            error("$license is not a known license choice, choose one of: $licenses.")
        end
        print(io, LICENSES[license](pkg, string(years), authors))
    end || info("License file exists, leaving unmodified; use `force=true` to overwrite")
end

function readme(pkg::String, user::String=""; force::Bool=false)
    genfile(pkg,"README.md",force) do io
        println(io, "# $pkg")
        isempty(user) && return
        url = "https://travis-ci.org/$user/$pkg.jl"
        println(io, "\n[![Build Status]($url.svg?branch=master)]($url)")
    end
end

function tests(pkg::String; force::Bool=false)
    genfile(pkg,"test/runtests.jl",force) do io
        print(io, """
        using $pkg
        using Base.Test

        # write your own tests here
        @test 1 == 1
        """)
    end
end

function travis(pkg::String; force::Bool=false)
    genfile(pkg,".travis.yml",force) do io
        print(io, """
        language: cpp
        compiler:
          - clang
        notifications:
          email: false
        env:
          matrix:
            - JULIAVERSION="juliareleases"
            - JULIAVERSION="julianightlies"
        before_install:
          - sudo add-apt-repository ppa:staticfloat/julia-deps -y
          - sudo add-apt-repository ppa:staticfloat/\${JULIAVERSION} -y
          - sudo apt-get update -qq -y
          - sudo apt-get install libpcre3-dev julia -y
          - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
        script:
          - julia -e 'Pkg.init(); Pkg.clone(pwd()); Pkg.test("$pkg")'
        """)
    end
end

function entrypoint(pkg::String; force::Bool=false)
    genfile(pkg,"src/$pkg.jl",force) do io
        print(io, """
        module $pkg

        # package code goes here

        end # module
        """)
    end
end

function genfile(f::Function, pkg::String, file::String, force::Bool=false)
    path = joinpath(pkg,file)
    if force || !ispath(path)
        info("Generating $file")
        mkpath(dirname(path))
        open(f, path, "w")
        Git.run(`add $file`, dir=pkg)
        return true
    end
    return false
end

copyright(years::String, authors::String) = "> Copyright (c) $years: $authors."

function copyright(years::String, authors::Array)
    text = "> Copyright (c) $years:"
    for author in authors
        text *= "\n>  * $author"
    end
    return text
end

mit(pkg::String, years::String, authors::Union(String,Array)) =
"""
The $pkg.jl package is licensed under the MIT "Expat" License:

$(copyright(years,authors))
>
> 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.
"""

bsd(pkg::String, years::String, authors::Union(String,Array)) =
"""
The $pkg.jl package is licensed under the Simplified "2-clause" BSD License:

$(copyright(years,authors))
>
> Redistribution and use in source and binary forms, with or without
> modification, are permitted provided that the following conditions are
> met:
>
> 1. Redistributions of source code must retain the above copyright
>    notice, this list of conditions and the following disclaimer.
> 2. Redistributions in binary form must reproduce the above copyright
>    notice, this list of conditions and the following disclaimer in the
>    documentation and/or other materials provided with the distribution.
>
> THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
> "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
> LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
> A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
> OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
> SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
> LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
> DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
> THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
> (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
> OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
"""

const LICENSES = [ "MIT" => mit, "BSD" => bsd ]

end # module

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
""" `train!(model::AbstractGP;iterations::Integer=100,callback=0,conv_function=0)`

Function to train the given GP `model`.

**Keyword Arguments**

there are options to change the number of max iterations,
- `iterations::Int` : Number of iterations (not necessarily epochs!)for training
- `callback::Function` : Callback function called at every iteration. Should be of type `function(model,iter) ...  end`
- `conv_function::Function` : Convergence function to be called every iteration, should return a scalar and take the same arguments as `callback`
"""
function train!(model::AbstractGP;iterations::Integer=100,callback=0,Convergence=0)
    if model.verbose > 0
      println("Starting training $model with $(model.nSample) samples with $(size(model.X,2)) features and $(model.nLatent) latent GP"*(model.nLatent > 1 ? "s" : ""))# using the "*model.Name*" model")
    end

    @assert iterations > 0  "Number of iterations should be positive"
    # model.evol_conv = [] #Array to check on the evolution of convergence
    local_iter::Int64 = 1; conv = Inf;

    while true #loop until one condition is matched
        try #Allow for keyboard interruption without losing the model
            update_parameters!(model) #Update all the variational parameters
            model.Trained = true
            if callback != 0
                callback(model,model.inference.nIter) #Use a callback method if put by user
            end
            if model.Autotuning && (model.inference.nIter%model.atfrequency == 0) && model.inference.nIter >= 3
                update_hyperparameters!(model) #Update the hyperparameters
            end
            ### Print out informations about the convergence
            if model.verbose > 2 || (model.verbose > 1  && local_iter%10==0)
                print("Iteration : $local_iter ")
                 print("ELBO is : $(ELBO(model))")
                 print("\n")
             end
            local_iter += 1; model.inference.nIter += 1
            (local_iter <= iterations) || break; #Verify if the number of maximum iterations has been reached
            # (iter < model.nEpochs && conv > model.ϵ) || break; #Verify if any condition has been broken
        catch e
            if isa(e,InterruptException)
                println("Training interrupted by user at iteration $local_iter");
                break;
            else
                rethrow(e)
            end
        end
    end
    if model.verbose > 0
      println("Training ended after $local_iter iterations. Total number of iterations $(model.inference.nIter)")
    end
    computeMatrices!(model) #Compute final version of the matrices for prediction
    post_process!(model)
    model.Trained = true
end

function update_parameters!(model::GP)
    computeMatrices!(model); #Recompute the matrices if necessary (when hyperparameters have been updated)
end

function update_parameters!(model::VGP)
    computeMatrices!(model); #Recompute the matrices if necessary (always for the stochastic case, or when hyperparameters have been updated)
    variational_updates!(model);
end

function update_parameters!(model::SVGP)
    if model.inference.Stochastic
        model.inference.MBIndices = StatsBase.sample(1:model.inference.nSamples,model.inference.nSamplesUsed,replace=false) #Sample nSamplesUsed indices for the minibatches
    end
    computeMatrices!(model); #Recompute the matrices if necessary (always for the stochastic case, or when hyperparameters have been updated)
    variational_updates!(model);
end

function computeMatrices!(model::GP{<:Likelihood,<:Inference,T}) where {T<:Real}
    if model.inference.HyperParametersUpdated
        model.Knn .= Symmetric.(KernelModule.kernelmatrix.([model.X],model.kernel) .+ model.likelihood.ϵ.*[I])
        model.invKnn .= inv.(model.Knn)
    end
end

function computeMatrices!(model::VGP{<:Likelihood,<:Inference,T}) where {T<:Real}
    if model.inference.HyperParametersUpdated
        model.Knn .= Symmetric.(KernelModule.kernelmatrix.([model.X],model.kernel) .+ convert(T,Jittering()).*[I])
        model.invKnn .= inv.(model.Knn)
    end
end

function computeMatrices!(model::SVGP{<:Likelihood,<:Inference,T}) where {T<:Real}
    if model.inference.HyperParametersUpdated
        model.Kmm .= broadcast((Z,kernel)->Symmetric(KernelModule.kernelmatrix(Z,kernel)+getvariance(kernel)*convert(T,Jittering())*I),model.Z,model.kernel)
        model.invKmm .= inv.(model.Kmm)
    end
    #If change of hyperparameters or if stochatic
    if model.inference.HyperParametersUpdated || model.inference.Stochastic
        KernelModule.kernelmatrix!.(model.Knm,[model.X[model.inference.MBIndices,:]],model.Z,model.kernel)
        model.κ .= model.Knm.*model.invKmm
        model.K̃ .= kerneldiagmatrix.([model.X[model.inference.MBIndices,:]],model.kernel) .+ [convert(T,Jittering())*ones(T,model.inference.nSamplesUsed)] - opt_diag.(model.κ,model.Knm)
        @assert sum(count.(broadcast(x->x.<0,model.K̃)))==0 "K̃ has negative values"
    end
    model.inference.HyperParametersUpdated=false
end


text:
#######
## Common Regressor methods
########
struct LinearFitresult{T, F<:Real, M<:AbstractArray{F}} <: MMI.MLJType
    sol_matrix::M
    bias::Bool
    names::T
end

_convert(common_type, x::AbstractVector) = convert(AbstractVector{common_type}, x)
_convert(common_type, x::AbstractMatrix) = convert(AbstractMatrix{common_type}, MMI.matrix(x))
matrix_(X::AbstractVector) = X
matrix_(X) = MMI.matrix(X)
_names(y::AbstractVector) = nothing
_names(Y) = collect(MMI.schema(Y).names)

function _matrix(X, target)
    Xmatrix_ = MMI.matrix(X)
    Y_ = matrix_(target)
    common_type = promote_type(eltype(Xmatrix_), eltype(Y_))
    Xmatrix = _convert(common_type, Xmatrix_)
    Y = _convert(common_type, Y_)
    return Xmatrix, Y, _names(target)
end

function _regressor_fitted_params(fr::LinearFitresult{Nothing, <:Real, <:AbstractVector})
    return (
        coefficients=fr.sol_matrix[1:end-Int(fr.bias)],
        intercept=ifelse(fr.bias, fr.sol_matrix[end], nothing)
    )
end

function _regressor_fitted_params(fr::LinearFitresult{<:Vector, <:Real, <:AbstractMatrix})
    return (
        coefficients=fr.sol_matrix[1:end-Int(fr.bias), :],
        intercept=fr.bias ? fr.sol_matrix[end, :] : nothing
    )
end

function _predict_regressor(
    fr::LinearFitresult{Nothing, <:Real, <:AbstractVector},
    Xmat_new::AbstractMatrix,
    prototype
)
    if fr.bias
        return @views Xmat_new * fr.sol_matrix[1:end-1] .+ transpose(fr.sol_matrix[end])
    else
        return @views Xmat_new * fr.sol_matrix[1:end-1]
    end
end

function _predict_regressor(
    fr::LinearFitresult{<:Vector, <:Real, <:AbstractMatrix},
    Xmat_new::AbstractMatrix,
    prototype
)
    if fr.bias
        return MMI.table(
            Xmat_new * @view(fr.sol_matrix[1:end-1, :]) .+ transpose(
                @view(fr.sol_matrix[end, :])
            );
            names=fr.names,
            prototype=prototype
        )
    else
        return MMI.table(
            Xmat_new * @view(fr.sol_matrix[1:end-1, :]);
            names=fr.names,
            prototype=prototype
        )
    end
end

####
#### LinearRegressor & MultitargetLinearRegressor
####

"""
    LinearRegressor(; bias::Bool=true)

$LinearRegressor_DESCR

# Keyword Parameters

- `bias::Bool=true`: if true includes a bias term else fits without bias term.
"""
@mlj_model mutable struct LinearRegressor <: MMI.Deterministic
    bias::Bool = true
end

"""
    MultitargetLinearRegressor(; bias::Bool=true)

$MultitargetLinearRegressor_DESCR

# Keyword Parameters

- `bias::Bool=true`: if true includes a bias term else fits without bias term.
"""
@mlj_model mutable struct MultitargetLinearRegressor <: MMI.Deterministic
    bias::Bool = true
end

const LINREG = Union{LinearRegressor, MultitargetLinearRegressor}

function MMI.fit(model::LINREG, verbosity::Int, X, y)
    Xmatrix, y_, target_header= _matrix(X, y)
    θ = MS.llsq(Xmatrix, y_; bias=model.bias)
    fitresult = LinearFitresult(θ, model.bias, target_header)
    report = NamedTuple()
    cache = nothing
    return fitresult, cache, report
end

function MMI.fitted_params(::LINREG, fr)
    return _regressor_fitted_params(fr)
end

function MMI.predict(::LINREG, fr, Xnew)
    Xmat_new = MMI.matrix(Xnew)
    return _predict_regressor(fr, Xmat_new, Xnew)
end

####
#### RidgeRegressor & MultitargetRidgeRegressor
####

_check_typeof_lambda(x)= x isa AbstractVecOrMat || (x isa Real && x ≥ 0)

"""
    RidgeRegressor(; lambda::Union{Real, AbstractVecOrMat}=1.0, bias::Bool=true)

$RidgeRegressor_DESCR

# Keyword Parameters

- `lambda::Union{Real, AbstractVecOrMat}=1.0`: non-negative parameter for the
    regularization strength.
- `bias::Bool=true`: if true includes a bias term else fits without bias term.
"""
@mlj_model mutable struct RidgeRegressor <: MMI.Deterministic
    lambda::Union{Real, AbstractVecOrMat} = 1.0::(_check_typeof_lambda(_))
    bias::Bool = true
end

"""
    MultitargetRidgeRegressor(; lambda::Union{Real, AbstractVecOrMat}=1.0, bias::Bool=true)

$MultitargetRidgeRegressor_DESCR

# Keyword Parameters

- `lambda::Union{Real, AbstractVecOrMat}=1.0`: non-negative parameter for the
    regularization strength.
- `bias::Bool=true`: if true includes a bias term else fits without bias term.
"""
@mlj_model mutable struct MultitargetRidgeRegressor <: MMI.Deterministic
    lambda::Union{Real, AbstractVecOrMat} = 1.0::(_check_typeof_lambda(_))
    bias::Bool = true
end

const RIDGEREG = Union{RidgeRegressor, MultitargetRidgeRegressor}

function MMI.fit(model::RIDGEREG, verbosity::Int, X, y)
    Xmatrix, y_, target_header = _matrix(X, y)
    θ = MS.ridge(Xmatrix, y_, model.lambda; bias=model.bias)
    fitresult = LinearFitresult(θ, model.bias, target_header)
    report = NamedTuple()
    cache = nothing
    return fitresult, cache, report
end

function MMI.fitted_params(::RIDGEREG, fr)
    return _regressor_fitted_params(fr)
end

function MMI.predict(::RIDGEREG, fr, Xnew)
    Xmat_new = MMI.matrix(Xnew)
    return _predict_regressor(fr, Xmat_new, Xnew)
end

############
### Models Metadata
############
metadata_model(
    LinearRegressor,
    input=Table(Continuous),
    target=AbstractVector{Continuous},
    weights=false,
    descr=LinearRegressor_DESCR,
    path="$(PKG).LinearRegressor"
)

metadata_model(
    MultitargetLinearRegressor,
    input=Table(Continuous),
    target=Table(Continuous),
    weights=false,
    descr=MultitargetLinearRegressor_DESCR,
    path="$(PKG).MultitargetLinearRegressor"
)

metadata_model(
    RidgeRegressor,
    input=Table(Continuous),
    target=AbstractVector{Continuous},
    weights=false,
    descr=RidgeRegressor_DESCR ,
    path="$(PKG).RidgeRegressor"
)

metadata_model(
    MultitargetRidgeRegressor,
    input=Table(Continuous),
    target=Table(Continuous),
    weights=false,
    descr=MultitargetRidgeRegressor_DESCR,
    path="$(PKG).MultitargetRidgeRegressor"
)


text:
function chipher(str)
  result = ""
  for i in  1:length(str)
    if ismatch(r"[0-9A-Za-z]", str[i:i])
      result *= string(219 - int(str[i]))
    else
      result *= str[i:i]
    end
  end
  result
end

println(chipher("Today is a good day to die"))


text:
"""
    AbstractInverseRetractionMethod

Abstract type for methods for inverting a retraction (see [`inverse_retract`](@ref)).
"""
abstract type AbstractInverseRetractionMethod end

"""
    AbstractRetractionMethod

Abstract type for methods for [`retract`](@ref)ing a tangent vector to a manifold.
"""
abstract type AbstractRetractionMethod end

"""
    ApproximateInverseRetraction <: AbstractInverseRetractionMethod

An abstract type for representing approximate inverse retraction methods.
"""
abstract type ApproximateInverseRetraction <: AbstractInverseRetractionMethod end

"""
    ApproximateRetraction <: AbstractInverseRetractionMethod

An abstract type for representing approximate retraction methods.
"""
abstract type ApproximateRetraction <: AbstractRetractionMethod end

"""
    ExponentialRetraction

Retraction using the exponential map.
"""
struct ExponentialRetraction <: AbstractRetractionMethod end

"""
    PolarRetraction <: AbstractRetractionMethod

Retractions that are based on singular value decompositions of the matrix / matrices
for point and tangent vector on a [`AbstractManifold`](@ref)
"""
struct PolarRetraction <: AbstractRetractionMethod end

"""
    ProjectionRetraction <: AbstractRetractionMethod

Retractions that are based on projection and usually addition in the embedding.
"""
struct ProjectionRetraction <: AbstractRetractionMethod end

"""
    QRRetraction <: AbstractRetractionMethod

Retractions that are based on a QR decomposition of the
matrix / matrices for point and tangent vector on a [`AbstractManifold`](@ref)
"""
struct QRRetraction <: AbstractRetractionMethod end

"""
    LogarithmicInverseRetraction

Inverse retraction using the [`log`](@ref)arithmic map.
"""
struct LogarithmicInverseRetraction <: AbstractInverseRetractionMethod end

"""
    PolarInverseRetraction <: AbstractInverseRetractionMethod

Inverse retractions that are based on a singular value decomposition of the
matrix / matrices for point and tangent vector on a [`AbstractManifold`](@ref)
"""
struct PolarInverseRetraction <: AbstractInverseRetractionMethod end

"""
    ProjectionInverseRetraction <: AbstractInverseRetractionMethod

Inverse retractions that are based on a projection (or its inversion).
"""
struct ProjectionInverseRetraction <: AbstractInverseRetractionMethod end

"""
    QRInverseRetraction <: AbstractInverseRetractionMethod

Inverse retractions that are based on a QR decomposition of the
matrix / matrices for point and tangent vector on a [`AbstractManifold`](@ref)
"""
struct QRInverseRetraction <: AbstractInverseRetractionMethod end

"""
    NLsolveInverseRetraction{T<:AbstractRetractionMethod,TV,TK} <:
        ApproximateInverseRetraction

An inverse retraction method for approximating the inverse of a retraction using `NLsolve`.

# Constructor

    NLsolveInverseRetraction(
        method::AbstractRetractionMethod[, X0];
        project_tangent=false,
        project_point=false,
        nlsolve_kwargs...,
    )

Constructs an approximate inverse retraction for the retraction `method` with initial guess
`X0`, defaulting to the zero vector. If `project_tangent` is `true`, then the tangent
vector is projected before the retraction using `project`. If `project_point` is `true`,
then the resulting point is projected after the retraction. `nlsolve_kwargs` are keyword
arguments passed to `NLsolve.nlsolve`.
"""
struct NLsolveInverseRetraction{TR<:AbstractRetractionMethod,TV,TK} <:
       ApproximateInverseRetraction
    retraction::TR
    X0::TV
    project_tangent::Bool
    project_point::Bool
    nlsolve_kwargs::TK
    function NLsolveInverseRetraction(m, X0, project_point, project_tangent, nlsolve_kwargs)
        return new{typeof(m),typeof(X0),typeof(nlsolve_kwargs)}(
            m,
            X0,
            project_point,
            project_tangent,
            nlsolve_kwargs,
        )
    end
end
function NLsolveInverseRetraction(
    m,
    X0 = nothing;
    project_tangent::Bool = false,
    project_point::Bool = false,
    nlsolve_kwargs...,
)
    return NLsolveInverseRetraction(m, X0, project_point, project_tangent, nlsolve_kwargs)
end

"""
    default_inverse_retraction_method(M::AbstractManifold)

The [`AbstractInverseRetractionMethod`](@ref) that is used when calling
[`inverse_retract`](@ref) without specifying the inverse retraction method.
By default, this is the [`LogarithmicInverseRetraction`](@ref).
"""
default_inverse_retraction_method(::AbstractManifold) = LogarithmicInverseRetraction()

"""
    default_retraction_method(M::AbstractManifold)

The [`AbstractRetractionMethod`](@ref) that is used when calling [`retract`](@ref) without
specifying the retraction method. By default, this is the [`ExponentialRetraction`](@ref).
"""
default_retraction_method(::AbstractManifold) = ExponentialRetraction()

"""
    inverse_retract!(M::AbstractManifold, X, p, q[, method::AbstractInverseRetractionMethod])

Compute the inverse retraction, a cheaper, approximate version of the
[`log`](@ref)arithmic map), of points `p` and `q` on the [`AbstractManifold`](@ref) `M`.
Result is saved to `X`.

Inverse retraction method can be specified by the last argument, defaulting to
[`default_inverse_retraction_method`](@ref)`(M)`. See the documentation of respective manifolds for
available methods.

See also [`retract!`](@ref).
"""
function inverse_retract!(M::AbstractManifold, X, p, q)
    return inverse_retract!(M, X, p, q, default_inverse_retraction_method(M))
end
function inverse_retract!(
    M::AbstractManifold,
    X,
    p,
    q,
    method::LogarithmicInverseRetraction,
)
    return log!(M, X, p, q)
end
function inverse_retract!(
    M::AbstractManifold,
    X,
    p,
    q,
    method::AbstractInverseRetractionMethod,
)
    return error(
        manifold_function_not_implemented_message(M, inverse_retract!, X, p, q, method),
    )
end

"""
    inverse_retract(M::AbstractManifold, p, q)
    inverse_retract(M::AbstractManifold, p, q, method::AbstractInverseRetractionMethod

Compute the inverse retraction, a cheaper, approximate version of the
[`log`](@ref)arithmic map), of points `p` and `q` on the [`AbstractManifold`](@ref) `M`.

Inverse retraction method can be specified by the last argument, defaulting to
[`default_inverse_retraction_method`](@ref)`(M)`.
For available inverse retractions on certain manifolds see the documentation on the
corresponding manifold.

See also [`retract`](@ref).
"""
function inverse_retract(M::AbstractManifold, p, q)
    X = allocate_result(M, inverse_retract, p, q)
    inverse_retract!(M, X, p, q)
    return X
end
function inverse_retract(M::AbstractManifold, p, q, method::AbstractInverseRetractionMethod)
    X = allocate_result(M, inverse_retract, p, q)
    inverse_retract!(M, X, p, q, method)
    return X
end

@doc raw"""
    retract(M::AbstractManifold, p, X)
    retract(M::AbstractManifold, p, X, t::Real=1)
    retract(M::AbstractManifold, p, X, method::AbstractRetractionMethod)
    retract(M::AbstractManifold, p, X, t::Real=1, method::AbstractRetractionMethod)

Compute a retraction, a cheaper, approximate version of the [`exp`](@ref)onential map,
from `p` into direction `X`, scaled by `t`, on the [`AbstractManifold`](@ref) `M`.

A retraction ``\operatorname{retr}_p: T_p\mathcal M → \mathcal M`` is a smooth map that fulfills

1. ``\operatorname{retr}_p(0) = p``
2. ``D\operatorname{retr}_p(0): T_p\mathcal M \to T_p\mathcal M`` is the identity map, i.e. ``D\operatorname{retr}_p(0)[X]=X``,

where ``D\operatorname{retr}_p`` denotes the differential of the retraction

The retraction is called of second order if for all ``X`` the curves ``c(t) = R_p(tX)``
have a zero acceleration at ``t=0``, i.e. ``c''(0) = 0``.

Retraction method can be specified by the last argument, defaulting to
[`default_retraction_method`](@ref)`(M)`. For further available retractions see the documentation of respective manifolds.

Locally, the retraction is invertible. For the inverse operation, see [`inverse_retract`](@ref).
"""
function retract(M::AbstractManifold, p, X)
    q = allocate_result(M, retract, p, X)
    retract!(M, q, p, X)
    return q
end
retract(M::AbstractManifold, p, X, t::Real) = retract(M, p, t * X)
function retract(M::AbstractManifold, p, X, method::AbstractRetractionMethod)
    q = allocate_result(M, retract, p, X)
    retract!(M, q, p, X, method)
    return q
end
function retract(M::AbstractManifold, p, X, t::Real, method::AbstractRetractionMethod)
    return retract(M, p, t * X, method)
end

"""
    retract!(M::AbstractManifold, q, p, X)
    retract!(M::AbstractManifold, q, p, X, t::Real=1)
    retract!(M::AbstractManifold, q, p, X, method::AbstractRetractionMethod)
    retract!(M::AbstractManifold, q, p, X, t::Real=1, method::AbstractRetractionMethod)

Compute a retraction, a cheaper, approximate version of the [`exp`](@ref)onential map,
from `p` into direction `X`, scaled by `t`, on the [`AbstractManifold`](@ref) manifold `M`.
Result is saved to `q`.

Retraction method can be specified by the last argument, defaulting to
[`default_retraction_method`](@ref)`(M)`. See the documentation of respective manifolds for available
methods.

See [`retract`](@ref) for more details.
"""
retract!(M::AbstractManifold, q, p, X) = retract!(M, q, p, X, default_retraction_method(M))
retract!(M::AbstractManifold, q, p, X, t::Real) = retract!(M, q, p, t * X)
retract!(M::AbstractManifold, q, p, X, ::ExponentialRetraction) = exp!(M, q, p, X)
function retract!(M::AbstractManifold, q, p, X, t::Real, method::AbstractRetractionMethod)
    return retract!(M, q, p, t * X, method)
end
function retract!(M::AbstractManifold, q, p, X, method::AbstractRetractionMethod)
    return error(manifold_function_not_implemented_message(M, retract!, q, p, method))
end


text:
module SelectionProbit

using WGLMakie, GLMakie, Distributions, Colors, JSServe

include("plot.jl")

export showplot

end # module


text:
# Functions of reading phenotype, weights files and adding them to the model
#Read adjusted phenotypes
function read_yk(ykfile::String)
    yk = readdlm(ykfile, '\t', Float64, '\n'; header = false)
    return Phenotypes(map(x->string(Int64(x)), yk[:,1]), yk[:,2:end])
end
# LOSO adjusted
function setup_phenotypes(model, prefix, weightsfilename)
    model.phenotypesprefix = prefix
    model.phenotypesfilename =  string(model.phenotypesprefix, ".layer1.n1.tsv")
    model.phenotypes = read_yk(model.phenotypesfilename)

    model.weightsfilename = weightsfilename
    model.weights = read_weights(model.weightsfilename)

    model.phenoOrder, model.weightOrder = get_orders(model.phenotypes.IID, model.weights.IID)
    get_phenotype_intersect!(model.phenotypes, model.phenoOrder)
    get_weight_intersect!(model.weights, model.weightOrder)
end
#File with IID, Rinv, same order as phenotype file
function read_weights(weightsfile::String)
    w  = readdlm(weightsfile, ' ', Float64, '\n'; header = false)
    return Weights(map(x->string(Int64(x)), w[:,1]), w[:,2])
end
# Full MA adjusted (used full phenotype filename, not prefix)
function setup_phenotypes_full(model, phenotypesfilename, weightsfilename)
    model.phenotypesfilename =  phenotypesfilename
    model.phenotypes = read_yk(model.phenotypesfilename)

    model.weightsfilename = weightsfilename
    model.weights = read_weights(model.weightsfilename)

    model.phenoOrder, model.weightOrder = get_orders(model.phenotypes.IID, model.weights.IID)
    get_phenotype_intersect!(model.phenotypes, model.phenoOrder)
    get_weight_intersect!(model.weights, model.weightOrder)
end

max_stars_repo_name: egmreynolds/MANA

max_issues_repo_name: egmreynolds/MANA

max_forks_repo_name: egmreynolds/MANA


text:
using BinaryBuilder

# Collection of sources required to build Blosc
sources = [
    "https://github.com/Blosc/c-blosc.git" =>
    "88864b5fb759d12a0b32c47ad9f0436d00adea88", # v1.14.3
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/c-blosc

mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DBUILD_TESTS=Off -DBUILD_BENCHMARKS=Off -DCMAKE_SYSTEM_PROCESSOR=$(uname -p) ..
make && make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms() # build on all supported platforms

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

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "Blosc", v"1.14.3", sources, script, platforms, products, dependencies)

println("Contents of ", pwd(), " = ", readdir("."))
if isdir("products")
    println("PRODUCTS = ", readdir("products"))
else
    println("products/ directory not found")
end


text:
# TSP

using JuMP
using CPLEX

m = Model(solver=CplexSolver(CPX_PARAM_MIPDISPLAY=1, CPX_PARAM_MIPINTERVAL=1))

n = 6
cities = [ 50 200;
          100 100;
          100 300;
          500 100;
          500 300;
          550 200]
dist = zeros(n, n)
for i = 1:n
  for j = i:n
    d = norm(cities[i,1:2]-cities[j,1:2])
    dist[i,j] = d
    dist[j,i] = d
  end
end
#println(dist)

@variable(m, x[1:n,1:n], Bin)

@objective(m, Min, sum(dist[i,j]*x[i,j] for i=1:n,j=i:n))

# Make symmetric
for i = 1:n
  @constraint(m, x[i,i] == 0)
  for j = (i+1):n
    @constraint(m, x[i,j] == x[j,i])
  end
end

# Sum in to city = 1
for i = 1:n
  @constraint(m, sum(x[i,j] for j=1:n) == 2)
end

# Sum out of city = 1
for j = 1:n
  @constraint(m, sum(x[i,j] for i=1:n) == 2)
end

function subtour(cb)
  println("In subtour")
  cur_sol = getvalue(x)
  # println(cur_sol)

  # Find any subtour
  in_subtour = fill(false,n)
  cur_node = 1
  in_subtour[cur_node] = true
  subtour_length = 1
  while true
    # Find next unvisited node
    found_node = false
    for j = 1:n
      if !in_subtour[j]
        if cur_sol[cur_node,j] >= 0.9
          # Arc to unvisited node
          cur_node = j
          in_subtour[j] = true
          found_node = true
          subtour_length += 1
          break
        end
      end
    end
    if !found_node
      println("Done exploring")
      # Completely explored this subtour
      if subtour_length == n
        println("found tour!")
        # Done!
        break
      else
        # Add lazy constraint
        expr = AffExpr()
        for i = 1:n
          if !in_subtour[i]
            continue
          end
          # i is in S
          for j = 1:n
            if i == j
              continue
            end
            if in_subtour[j]
              # Both ends in subtour, bad
            else
              # j isn't in subtour
              expr += x[i,j]
            end
          end
        end
        println(expr)
        readline(STDIN)
        @lazyconstraint(cb, expr >= 2)
        break
      end
    end
  end

end

addlazycallback(m, subtour)
stat = solve(m)

println(getvalue(x))
println(stat)
println(getobjectivevalue(m))


text:
import Pkg
Pkg.add("Documenter")
using Documenter
using SPGBox
push!(LOAD_PATH, "../src/")
makedocs(
    modules = [SPGBox],
    sitename = "SPGBox.jl",
    pages = [
        "Installation" => "index.md",
        "User guide" => "usage.md",
        "Options" => "options.md",
        "Reference" => "reference.md",
    ],
)
deploydocs(
    repo = "github.com/m3g/SPGBox.jl.git",
    target = "build",
    branch = "gh-pages",
    versions = ["stable" => "v^", "v#.#"],
)


text:
# using Pkg
# Pkg.add("NLopt")
# Pkg.add("RecipesBase")

using LikelihoodProfiler, Test

include("./cases_func.jl")
methods_list = [:CICO_ONE_PASS, :LIN_EXTRAPOL, :QUADR_EXTRAPOL]

println("Starting tests for get_interval")
@testset "get_interval" begin include("test_get_interval.jl") end

println("Starting tests for get_endpoint")
@testset "get_endpoint CICO_ONE_PASS parameters" begin include("test_get_endpoint_parameters.jl") end
@testset "get_endpoint CICO_ONE_PASS function" begin include("test_get_endpoint_scan_func.jl") end
@testset "get_endpoint errors and warnings" begin include("test_get_endpoint_errors.jl") end

println("Starting tests for all all methods of get_right_endpoint")
@testset "get_right_endpoint methods" begin include("test_methods.jl") end
@testset "get_right_endpoint CICO_ONE_PASS" begin include("test_cico_one_pass.jl") end

println("Starting tests for profile")
@testset "profile" begin include("test_profile.jl") end

println("Starting tests for different fitting alg")
@testset "local_alg" begin include("test_algorithms.jl") end

println("Starting tests for Plot @recipe")
@testset "PLOT_INTERVAL" begin include("test_plots.jl") end

println("Starting tests for loss error")
@testset "LOSS_ERROR" begin include("test_loss_error.jl") end

println("Starting tests for bands")
@testset "Bands" begin include("test_bands.jl") end

println("Starting tests for get_optimal")
@testset "get_optimal(...)" begin include("test_get_optimal.jl") end

# experimental tests

@testset "testing derivative-free algorithms" begin include("test_deriv_free_algs.jl") end
@testset "gradient-based algorithms" begin include("test_grad_algs.jl") end


text:
using RiskAdjustedLinearizations, Test, JLD2
include(joinpath(dirname(@__FILE__), "../../examples/wachter_disaster_risk/wachter.jl"))

# Get stochastic steady state
sssout = JLD2.jldopen(joinpath(dirname(@__FILE__), "../reference/iterative_sss_output.jld2"), "r")
z = vec(sssout["z"])
y = vec(sssout["y"])
Ψ = sssout["Psi"]

# Create model and update it
m = WachterDisasterRisk()
ral = inplace_wachter_disaster_risk(m)
RiskAdjustedLinearizations.update!(ral, z, y, Ψ)
@info "The following message about Blanchard-Kahn conditions is expected."
@test RiskAdjustedLinearizations.blanchard_kahn(ral)
@test RiskAdjustedLinearizations.blanchard_kahn(ral; verbose = :none)

nothing


text:
using CommonEq
using Test

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


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

JLLWrappers.@generate_main_file_header("libigc")
JLLWrappers.@generate_main_file("libigc", UUID("94295238-5935-5bd7-bb0f-b00942e9bdd5"))
end  # module libigc_jll


text:
# Copyright (c) 2013: Clp.jl contributors
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module TestMOIWrapper

using Test
using MathOptInterface
import Clp

const MOI = MathOptInterface

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

function test_SolverName()
    @test MOI.get(Clp.Optimizer(), MOI.SolverName()) == "Clp"
    return
end

function test_supports_default_copy_to()
    @test !MOI.supports_incremental_interface(Clp.Optimizer())
    return
end

function test_runtests()
    # This is what JuMP would construct
    model = MOI.Utilities.CachingOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
        MOI.instantiate(Clp.Optimizer; with_bridge_type = Float64),
    )
    @test model.optimizer.model.model_cache isa
          MOI.Utilities.UniversalFallback{Clp.OptimizerCache}
    # `Variable.ZerosBridge` makes dual needed by some tests fail.
    MOI.Bridges.remove_bridge(
        model.optimizer,
        MathOptInterface.Bridges.Variable.ZerosBridge{Float64},
    )
    MOI.set(model, MOI.Silent(), true)
    MOI.Test.runtests(
        model,
        MOI.Test.Config(
            exclude = Any[MOI.DualObjectiveValue, MOI.ObjectiveBound],
        ),
        exclude = [
            # Unable to prove infeasibility
            "test_conic_NormInfinityCone_INFEASIBLE",
            "test_conic_NormOneCone_INFEASIBLE",
        ],
    )
    return
end

function test_Nonexistant_unbounded_ray()
    inner = Clp.Optimizer()
    model =
        MOI.Utilities.CachingOptimizer(MOI.default_cache(inner, Float64), inner)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variables(model, 5)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.optimize!(model)
    status = MOI.get(model, MOI.TerminationStatus())
    @test status == MOI.DUAL_INFEASIBLE
    return
end

function test_RawOptimizerAttribute()
    model = Clp.Optimizer()
    MOI.set(model, MOI.RawOptimizerAttribute("LogLevel"), 1)
    @test MOI.get(model, MOI.RawOptimizerAttribute("LogLevel")) == 1
    MOI.set(model, MOI.RawOptimizerAttribute("LogLevel"), 2)
    @test MOI.get(model, MOI.RawOptimizerAttribute("LogLevel")) == 2

    MOI.set(model, MOI.RawOptimizerAttribute("SolveType"), 1)
    @test MOI.get(model, MOI.RawOptimizerAttribute("SolveType")) == 1
    MOI.set(model, MOI.RawOptimizerAttribute("SolveType"), 4)
    @test MOI.get(model, MOI.RawOptimizerAttribute("SolveType")) == 4

    MOI.set(model, MOI.RawOptimizerAttribute("PresolveType"), 1)
    @test MOI.get(model, MOI.RawOptimizerAttribute("PresolveType")) == 1
    MOI.set(model, MOI.RawOptimizerAttribute("PresolveType"), 0)
    @test MOI.get(model, MOI.RawOptimizerAttribute("PresolveType")) == 0
    return
end

function test_All_parameters()
    model = Clp.Optimizer()
    param = MOI.RawOptimizerAttribute("NotAnOption")
    @test !MOI.supports(model, param)
    @test_throws MOI.UnsupportedAttribute(param) MOI.get(model, param)
    @test_throws MOI.UnsupportedAttribute(param) MOI.set(model, param, false)
    for key in Clp.SUPPORTED_PARAMETERS
        @test MOI.supports(model, MOI.RawOptimizerAttribute(key))
        value = MOI.get(model, MOI.RawOptimizerAttribute(key))
        MOI.set(model, MOI.RawOptimizerAttribute(key), value)
        @test MOI.get(model, MOI.RawOptimizerAttribute(key)) == value
    end
    return
end

function test_copy_to_bug()
    model = MOI.Utilities.Model{Float64}()
    x = MOI.add_variable(model)
    con = [
        MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
            MOI.EqualTo(1.0),
        ) for i in 1:2
    ]
    clp = Clp.Optimizer()
    index_map = MOI.copy_to(clp, model)
    @test index_map[con[1]] != index_map[con[2]]
    return
end

function test_options_after_empty!()
    model = Clp.Optimizer()
    @test MOI.get(model, MOI.Silent()) == false
    MOI.set(model, MOI.Silent(), true)
    @test MOI.get(model, MOI.Silent()) == true
    MOI.empty!(model)
    @test MOI.get(model, MOI.Silent()) == true
    return
end

end  # module TestMOIWrapper

TestMOIWrapper.runtests()


text:
push!(LOAD_PATH, "./koans")

using Test

using ArithmeticKoans
@testset "Arithmetic Operators" begin
    @test ArithmeticKoans.sum(1, 2) == 3
    @test ArithmeticKoans.difference(8, 1) == 7
    @test ArithmeticKoans.multiplication(10, 2) == 20
    @test ArithmeticKoans.division(35, 5) == 7.0
    @test ArithmeticKoans.integer_division(35, 5) == 7
    @test ArithmeticKoans.power(2, 2) == 4
    @test ArithmeticKoans.modulo(12, 10) == 2
end

using LogicalKoans
@testset "Logical Operators" begin
    @test LogicalKoans.not(true) == false
    @test LogicalKoans.not(false) == true
    @test LogicalKoans.equals(true, true) == true
    @test LogicalKoans.equals(false, false) == true
    @test LogicalKoans.not_equals(true, false) == true
    @test LogicalKoans.less_than(1, 10) == true
    @test LogicalKoans.greater_than(10, 1) == true
    @test LogicalKoans.less_than_or_equal_to(2, 2) == true
    @test LogicalKoans.less_than_or_equal_to(1, 10) == true
    @test LogicalKoans.less_than_or_equal_to(10, 1) == false
    @test LogicalKoans.greater_than_or_equal_to(2, 2) == true
    @test LogicalKoans.greater_than_or_equal_to(10, 1) == true
    @test LogicalKoans.greater_than_or_equal_to(1, 10) == false
end

using BitwiseKoans
@testset "Bitwise Operators" begin
    @test BitwiseKoans.not(2) == -3
    @test BitwiseKoans.and(3, 5) == 1
    @test BitwiseKoans.or(2, 4) == 6
    @test BitwiseKoans.xor_koan(2, 4) == 6
    @test BitwiseKoans.right_shift(2, 1) == 1
    @test BitwiseKoans.left_shift(2, 1) == 4
end

using ComplexNumberKoans
@testset "Complex Numbers" begin
end

using StringKoans
@testset "Strings" begin
end

using ConversionKoans
@testset "Type Conversions" begin
end

using ArrayKoans
@testset "Arrays" begin
end

using VectorOperationKoans
@testset "Vector Operators" begin
end

using MultiDimensionalArrayKoans
@testset "Multi-dimensional Arrays" begin
end

using IndexingKoans
@testset "Indexing" begin
end

using IterationKoans
@testset "Iteration" begin
end

using ModuleKoans
@testset "Modules" begin
end

using StructKoans
@testset "Structs" begin
end

using UnionKoans
@testset "Unions" begin
end

using FunctionKoans
@testset "Functions" begin
end

using HighOrderFunctionKoans
@testset "High Order Functions" begin
    f = (x, y) -> x + y
    @test HighOrderFunctionKoans.round_list([1.5, 2.3, 4.7]) == [2.0, 2.0, 5.0]
    @test HighOrderFunctionKoans.apply_fn(f, 7) == 12
    @test HighOrderFunctionKoans.add_n([1, 2, 4, 5],  3) == [4, 5, 7, 8]
    @test HighOrderFunctionKoans.remove_multiples([10, 12, 15, 18, 20], 3) == [10, 20]
end

using MethodKoans
@testset "Methods" begin
end

using DestructuringKoans
@testset "Destructuring" begin
    @test DestructuringKoans.add_separately(1+2im, 3-7im) == (4, -5)
    @test DestructuringKoans.add_coordinates((2, 3), (4, -1)) == (6, 2)
end

using MultipleDispatchKoans
@testset "Multiple Dispatch" begin
    dog = MultipleDispatchKoans.Dog("goodboi")
    cat = MultipleDispatchKoans.Cat("Sir Catcelot")
    @test MultipleDispatchKoans.hello() == "hello stranger"
    @test MultipleDispatchKoans.hello("julia") == "hello julia"
    @test MultipleDispatchKoans.pet_talk(dog) == "$(dog.name): bark"
    @test MultipleDispatchKoans.pet_talk(cat) == "$(cat.name): meow"
end

using ValueTypeKoans
@testset "Value Types" begin
end

using PythonInterop
@testset "Python Interop" begin
    @test PythonInterop.py_sinpi(0) == 0
    @test round(PythonInterop.py_sinpi(1/pi), Base.Rounding.RoundToZero, digits=2) == 0.84

    pyfib = PythonInterop.get_pyfib()
    @test PythonInterop.jlfib(1, PythonInterop.jlfib) == 1
    @test PythonInterop.jlfib(5, PythonInterop.jlfib) == 5
    @test PythonInterop.jlfib(5, pyfib) == 5
    @test PythonInterop.jlfib(7, pyfib) == 13
end


text:
function read_instructions(lines)
    instructions = Tuple{AbstractString, Int64}[]
    for line in lines
        op, value = split(line)
        push!(instructions, (op, parse(Int64, value)))
    end
    instructions
end

function find_possibly_corrupted_lines(instructions)
    indices = Int64[]
    for (ix, (op, _)) in enumerate(instructions)
        if op == "jmp" || op == "nop"
            push!(indices, ix)
        end
    end
    indices
end

function execute_correct_program(instructions)
    indices = find_possibly_corrupted_lines(instructions)
    for ix in indices
        modified_instructions = copy(instructions)
        op, k = modified_instructions[ix]
        if op == "jmp"
            modified_instructions[ix] = ("nop", k)
        elseif op == "nop"
            modified_instructions[ix] = ("jmp", k)
        end
        accumulator = run_program(modified_instructions)
        if accumulator !== nothing
            return accumulator
        end
    end
    return nothing
end

function run_program(instructions)
    n = length(instructions)
    visited = falses(n)
    ix = 1
    accumulator = 0
    while !visited[ix]
        (op, k) = instructions[ix]
        visited[ix] = true
        if op == "nop"
            ix += 1
        elseif op == "acc"
            accumulator += k
            ix += 1
        elseif op == "jmp"
            ix += k
        end

        if ix > n
            return accumulator
        end
    end
    nothing
end

function main()
    lines = readlines("../resources/input_08.txt")
    instructions = read_instructions(lines)
    result = execute_correct_program(instructions)
    println(result)
end

main()


text:
import Base: show, sum, conj, adjoint, isequal, (+)

export IntegerPartition, Ferrers

"""
An `IntegerPartition` is a decreasing list of positive integers.
Construct in the following ways:
* `IntegerParition([a,b,c,...])`
* `IntegerPartition(a,b,c,...)`
"""
struct IntegerPartition
    parts::Vector{Int}
    val::Int

    function IntegerPartition()
        parts = Vector{Int}()
        val   = 0
        new(parts,val)
    end

    function IntegerPartition(pts::Vector{Int})
        if any(pts .<= 0)
            error("All parts of an IntegerPartition must be positive")
        end
        parts = sort(pts, rev=true)
        value = sum(parts)
        new(parts, value)
    end

    function IntegerPartition(pts...)
        return IntegerPartition(collect(pts))
    end

end

parts(P::IntegerPartition) = copy(P.parts)
num_parts(P::IntegerPartition) = length(P.parts)
sum(P::IntegerPartition) = P.val

function show(io::IO, P::IntegerPartition)
    str = "("
    np = num_parts(P)
    for k=1:np
        str *= string(P.parts[k])
        if k < np
            str *= "+"
        end
    end
    str *= ")"
    print(io, str)
end

"""
`Ferrers(P::IntegerParition)` prints a graphical representation of the
partition `P` in the form of a Ferrer's diagram.
"""
function Ferrers(P::IntegerPartition, sym::Char='*')
    np = num_parts(P)
    for i=1:np
        println(sym^P.parts[i])
    end
    nothing
end

"""
`conj(P::IntegerPartition)` returns the Ferrer's conjugate of `P`.
Also available as `P'`.
"""
function conj(P::IntegerPartition)::IntegerPartition
    np = num_parts(P)
    if np == 0
        return IntegerPartition()
    end
    big = P.parts[1]  # largest part

    new_parts = Array{Int,1}(undef,big)   # Vector{Int}(big)
    for k=1:big
        new_parts[k] = count(P.parts .>= k)
    end
    return IntegerPartition(new_parts)
end

adjoint(P::IntegerPartition) = conj(P)

"""
The sum of `IntegerPartition`s is their concatenation (multiset union).
"""
function (+)(P::IntegerPartition,Q::IntegerPartition)
    pts = [parts(P);parts(Q)]
    return IntegerPartition(pts)
end

isequal(P::IntegerPartition,Q::IntegerPartition) = isequal(P.parts,Q.parts)
==(P::IntegerPartition, Q::IntegerPartition) = isequal(P,Q)


text:
const k8 = [ 4, 10,  9,  2, 13,  8,  0, 14,  6, 11,  1, 12,  7, 15,  5,  3]
const k7 = [14, 11,  4, 12,  6, 13, 15, 10,  2,  3,  8,  1,  0,  7,  5,  9]
const k6 = [ 5,  8,  1, 13, 10,  3,  4,  2, 14, 15, 12,  7,  6,  0,  9, 11]
const k5 = [ 7, 13, 10,  1,  0,  8,  9, 15, 14,  4,  6, 12, 11,  2,  5,  3]
const k4 = [ 6, 12,  7,  1,  5, 15, 13,  8,  4, 10,  9, 14,  0,  3, 11,  2]
const k3 = [ 4, 11, 10,  0,  7,  2,  1, 13,  3,  6,  8,  5,  9, 12, 15, 14]
const k2 = [13, 11,  4,  1,  3, 15,  5,  9,  0, 10, 14,  7,  6,  8,  2, 12]
const k1 = [ 1, 15, 13,  0,  5,  7, 10,  4,  9,  2,  3, 14,  6, 11,  8, 12]
const k87 = zeros(UInt32,256)
const k65 = zeros(UInt32,256)
const k43 = zeros(UInt32,256)
const k21 = zeros(UInt32,256)
for i in 1:256
    j = (i-1) >> 4 + 1
    k = (i-1) & 15 + 1
    k87[i] = (k1[j] << 4) | k2[k]
    k65[i] = (k3[j] << 4) | k4[k]
    k43[i] = (k5[j] << 4) | k6[k]
    k21[i] = (k7[j] << 4) | k8[k]
end

function f(x)
    y = (k87[(x>>24) & 0xff + 1] << 24) | (k65[(x>>16) & 0xff + 1] << 16) |
        (k43[(x>> 8) & 0xff + 1] <<  8) | k21[x & 0xff + 1]
    (y << 11) | (y >> (32-11))
end

bytes2int(arr) = arr[1] + (UInt32(arr[2]) << 8) + (UInt32(arr[3]) << 16) + (UInt32(arr[4])) << 24
int2bytes(x) = [UInt8(x&0xff), UInt8((x&0xff00)>>8), UInt8((x&0xff0000)>>16), UInt8(x>>24)]

function mainstep(inputbytes, keybytes)
    intkey = bytes2int(keybytes)
    lowint = bytes2int(inputbytes[1:4])
    topint = bytes2int(inputbytes[5:8])
    xorbytes = f(UInt32(intkey) + UInt32(lowint)) ⊻ topint
    vcat(int2bytes(xorbytes), inputbytes[1:4])
end

const input = [0x21, 0x04, 0x3B, 0x04, 0x30, 0x04, 0x32, 0x04]
const key = [0xF9, 0x04, 0xC1, 0xE2]
println("The encoded bytes are $(mainstep(input, key))")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module MaximumLikelihoodPower

"""
    module MaximumLikelihoodPower

Maximum likelihood calculations for power laws.

Functions: mle, KStatistic, scanmle, scanKS
"""
MaximumLikelihoodPower

export mle, KSstatistic, scanmle, scanKS, mleKS

include("mle.jl")
include("aux.jl")

end # module


text:
# ZJ: currenlty there is no support for distributed array and given the speed below it's not ready
using DistributedArrays
@time addprocs()
@everywhere using DistributedArrays
@everywhere using FastGroupBy

@time did6 = distribute(id6)
@time dv1 = distribute(v1)

function sumby_da(by, val)
  T = eltype(by)
  S = eltype(val)
  res = [@spawnat k sumby(localpart(by), localpart(val)) for k = 2:nprocs()]
  # algorithms to collate all dicts
  fnl_res = fetch(res[1])
  szero = zero(S)
  for i = 2:length(res)
    next_res = fetch(res[i])
    for k = keys(next_res)
      fnl_res[k] = get(fnl_res, k, szero) + next_res[k]
    end
  end
  fnl_res
end

@time sumby_da(did6, dv1)

@time pool1 = [@sprintf "id%010d" k for k in 1:(N/K)]
@time id3 = rand(pool1, N)

@time did3 = distribute(id3)

@time res = [@spawnat k sumby(hash.(localpart(did3)), localpart(dv1)) for k = 2:nprocs()]
@time res = [@spawnat k sumby_dict(localpart(did3), localpart(dv1)) for k = 2:nprocs()]
@time fetch.(res)

using CategoricalArrays
import CategoricalArrays.CategoricalArray
import CategoricalArrays.CategoricalPool
@time a = CategoricalArray(id3, ordered=false)

@time a = CategoricalArray(rand(UInt32(1):UInt32(N/K),N), CategoricalPool(pool1,false))


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct GAMSInput <: File
        name::String
        col::Array{Symbol,1}
    end

Read .map or .set file

# Arguments
- `name::String`: input file name
- `col::Array{Symbol,1}`: column names
"""
mutable struct GAMSInput <: File
    "input file name"
    name::String
    "column names"
    col::Array{Symbol,1}
end

function GAMSInput(; name, col, )
    GAMSInput(name, col, )
end

"""Get [`GAMSInput`](@ref) `name`."""
get_name(value::GAMSInput) = value.name
"""Get [`GAMSInput`](@ref) `col`."""
get_col(value::GAMSInput) = value.col

"""Set [`GAMSInput`](@ref) `name`."""
set_name!(value::GAMSInput, val) = value.name = val
"""Set [`GAMSInput`](@ref) `col`."""
set_col!(value::GAMSInput, val) = value.col = val

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
@testset "Ellipsoid" begin
    e = Ellipsoid([0,0], [1,2])
    @test e == deepcopy(e)
    @test hash(e) == hash(deepcopy(e))
    @test [0.3,2*sqrt(1 - 0.3^2)-0.01] ∈ e
    @test [0.3,2*sqrt(1 - 0.3^2)+0.01] ∉ e

    @test ((x,nout) = surfpt_nearby([0,0],e); (x≈[1,0] && nout≈[1,0]) || (x≈[0,2] && nout≈[0,1]))  # handle point at center properly
    @test all([surfpt_nearby([ρ*1sx,0],e)[1] ≈ [1sx,0] for ρ = (one⁻⁻,one⁺⁺), sx = (-1,1)])
    @test all([surfpt_nearby([0,ρ*2sy],e)[1] ≈ [0,2sy] for ρ = (one⁻⁻,one⁺⁺), sy = (-1,1)])
    @test all([(p = [1sx,0]; surfpt_nearby(p,e) ≈ (p,[sx,0])) for sx = (-1,1)])  # handle point on boundary properly
    @test all([(p = [0,2sy]; surfpt_nearby(p,e) ≈ (p,[0,sy])) for sy = (-1,1)])  # handle point on boundary properly

    @test normal([1.1,0],e) ≈ [1,0]
    @test normal([-1.1,0],e) ≈ [-1,0]
    @test normal([0,2.01],e) ≈ [0,1]
    @test bounds(e) == ([-1,-2],[1,2])
    @test checkbounds(e)
    @test checkbounds(Ellipsoid([0,0], [1,2], [1 1; 1 -1]))

    @test (∆ = rand(2); translate(e,∆) ≈ Ellipsoid([0,0]+∆, [1,2]))

    b = Box([0,0], [2,4])
    eb = Ellipsoid(b)
    @test e == eb
    @test hash(e) == hash(eb)
end  # @testset "Ellipsoid"

@testset "Ellipsoid, rotated" begin
    θ = π/3
    R = [cos(θ) sin(θ); sin(θ) -cos(θ)]
    er = Ellipsoid([0,0], [2,3], R)
    bp = GeometryPrimitives.boundpts(er)

    bp1, bp2 = bp[:,1], bp[:,2]

    # Test the two bounding points are on the ellipsoid perimeter.
    @test er == deepcopy(er)
    @test hash(er) == hash(deepcopy(er))
    @test (one⁻ * bp1 ∈ er) && (one⁻ * bp2 ∈ er)
    @test (one⁺ * bp1 ∉ er) && (one⁺ * bp2 ∉ er)

    @test ((x,nout) = surfpt_nearby([0,0],er); (x≈2*R[:,1] && nout≈R[:,1]) || (x≈3*R[:,2] && nout≈R[:,2]))  # handle point at center properly
    @test all([surfpt_nearby(R*[ρ*2sx,0],er)[1] ≈ R*[2sx,0] for ρ = (one⁻⁻,one⁺⁺), sx = (-1,1)])
    @test all([surfpt_nearby(R*[0,ρ*3sy],er)[1] ≈ R*[0,3sy] for ρ = (one⁻⁻,one⁺⁺), sy = (-1,1)])
    @test all([(p = R*[2sx,0]; surfpt_nearby(p,er) ≈ (p,R*[sx,0])) for sx = (-1,1)])  # handle point on boundary properly
    @test all([(p = R*[0,3sy]; surfpt_nearby(p,er) ≈ (p,R*[0,sy])) for sy = (-1,1)])  # handle point on boundary properly

    # Test the normal vector at the two bounding points are the x- and y-directions.
    @test normal(bp1, er) ≈ [1,0]
    @test normal(bp2, er) ≈ [0,1]

    xmax, ymax = bp1[1], bp2[2]
    @test bounds(er) == ([-xmax, -ymax], [xmax, ymax])
    @test checkbounds(er)

    @test (∆ = rand(2); translate(er,∆) ≈ Ellipsoid([0,0]+∆, [2,3], R))

    br = Box([0,0], 2*[2,3], R)
    ebr = Ellipsoid(br)
    @test er == ebr
    @test hash(er) == hash(ebr)
end  # @testset "Ellipsoid, rotated"