text:
#
# Sequencing views
# ================

# Abstract type and high level behaviour
# --------------------------------------

abstract type AbstractSequencingView end

# Specific concrete view types
# ----------------------------

# A basic view type - used most often.

struct BasicSequencingView <: AbstractSequencingView
    sequence::UInt64
    start::UInt32
    stop::UInt32
end

function BasicSequencingView(id::Integer, start::Integer, stop::Integer)
    return BasicSequencingView(convert(UInt64, id), convert(UInt32, start), convert(UInt32, stop))
end

@inline function BasicSequencingView(bsv::BasicSequencingView, start::Integer, stop::Integer)
    return BasicSequencingView(seqid(bsv), start, stop)
end

@inline seqid(x::BasicSequencingView) = x.sequence
@inline Base.first(x::BasicSequencingView) = x.start
@inline Base.last(x::BasicSequencingView) = x.stop

# A view type that includes meta info in the form of an additional UInt64 tag.

struct TaggedSequencingView <: AbstractSequencingView
    sequence::UInt64
    start::UInt32
    stop::UInt32
    tag::UInt64
end

function TaggedSequencingView(id::Integer, start::Integer, stop::Integer, tag::UInt64 = zero(UInt64))
    return TaggedSequencingView(id, convert(UInt32, start), convert(UInt32, stop), tag)
end

@inline function TaggedSequencingView(tsv::TaggedSequencingView, start::Integer, stop::Integer)
    return TaggedSequencingView(seqid(tsv), start, stop, tag(tsv))
end

@inline seqid(x::TaggedSequencingView) = x.sequence
@inline Base.first(x::TaggedSequencingView) = x.start
@inline Base.last(x::TaggedSequencingView) = x.stop
@inline tag(x::TaggedSequencingView) = x.tag

const BasicViews = Vector{BasicSequencingView}
const TaggedViews = Vector{TaggedSequencingView}

# Generic sequencing view methods
# -------------------------------

@inline isfwd(x::AbstractSequencingView) = first(x) < last(x)
@inline isbwd(x::AbstractSequencingView) = first(x) > last(x)

@inline function Base.length(x::AbstractSequencingView)
    return ifelse(isfwd(x), last(x) - first(x) + 1, first(x) - last(x) + 1)
end

@inline Base.maximum(x::AbstractSequencingView) = max(first(x), last(x))
@inline Base.minimum(x::AbstractSequencingView) = min(first(x), last(x))

@inline Base.UnitRange(x::AbstractSequencingView) = minimum(x):maximum(x)
@inline Base.convert(::Type{Base.UnitRange}, x::AbstractSequencingView) = Base.UnitRange(x)

@inline Base.firstindex(x::AbstractSequencingView) = 1
@inline Base.lastindex(x::AbstractSequencingView) = length(x)
@inline window(x::AbstractSequencingView) = first(x):last(x)
@inline Base.eachindex(x::AbstractSequencingView) = Base.OneTo(lastindex(x))

function Base.show(io::IO, x::AbstractSequencingView)
    print(io, "A ", length(x), "bp view over a sequence ")
    if isbwd(x)
        print(io, '(', last(x), " <=== ", first(x), ')')
    else
        print(io, '(', first(x), " ===> ", last(x), ')')
    end
end

@inline function Base.checkbounds(x::AbstractSequencingView, i::Integer)
    if 1 ≤ i ≤ length(x)
        return true
    end
    throw(BoundsError(x, i))
end

@inline function Base.getindex(x::AbstractSequencingView, i::Integer)
    @boundscheck checkbounds(x, i)
    return ifelse(isfwd(x), first(x) + i - 1, first(x) - i + 1)
end

@inline function Base.reverse(x::T) where {T<:AbstractSequencingView}
    return T(x, last(x), first(x))
end

@inline function subview(x::T, start::Integer, stop::Integer) where {T<:AbstractSequencingView}
    return T(x, x[start], x[stop])
end

function take_paired_ends(sv::AbstractSequencingView, flen::Int, rlen::Int)
    return subview(sv, firstindex(sv), flen), subview(sv, lastindex(sv), lastindex(sv) - rlen + 1)
end

function take_single_end(sv::AbstractSequencingView, len::Int = length(sv))
    if (len > length(sv)) | (len < 0)
        throw(ArgumentError("len must be 0 < len ≤ length(sv)"))
    end
    forwards = rand(Bool)
    start = ifelse(forwards, firstindex(sv), lastindex(sv))
    stop = ifelse(forwards, len, lastindex(sv) - len + 1)
    return subview(sv, start, stop)
end

function extract_sequence(ref::BioSequence, view::AbstractSequencingView)
    if isbwd(view)
        subseq = ref[last(view):first(view)]
        reverse_complement!(subseq)
    else
        subseq = ref[first(view):last(view)]
    end
    return subseq
end

function extract_sequence(genome::Vector{<:BioSequence}, view::AbstractSequencingView)
    return extract_sequence(genome[seqid(view)], view)
end

# Utils for dealing with multiple views
# -------------------------------------

# Transformations

## The amplify transformation

function amplify(views::Vector{T}, n::Integer) where {T<:AbstractSequencingView}
    newviews = Vector{T}(undef, length(views) * n)
    i = 0
    for view in views
        @inbounds for _ in 1:n
            i += 1
            newviews[i] = view
        end
    end
    return newviews
end

# Utility function, randomly select values from a range, without replacement.
function sample_values(range::UnitRange, n::Int)
    @assert n <= length(range)
    s = Set{eltype(range)}()
    while length(s) < n
        push!(s, rand(range))
    end
    return collect(s)
end

## The subsample transformation

function subsample(vs::Vector{T}, n::Int) where {T<:AbstractSequencingView}
    if length(vs) < n
        throw(ArgumentError("number of views to sample is greater than total number of views."))
    end
    p = sample_values(1:length(vs), n)
    return vs[p]
end

## The selection transformation

function select(fn::Function, vs::Vector{T}) where {T<:AbstractSequencingView}
    nselected = 0
    newviews = Vector{T}(undef, length(vs))
    @inbounds for oldview in vs
        prob = fn(oldview)
        if rand() < prob
            nselected += 1
            newviews[nselected] = oldview
        end
    end
    resize!(newviews, nselected)
    return newviews
end

## The tag transformation

function tag(vs::BasicViews, npools::Int)
    newviews = TaggedViews(undef, length(vs))
    tags = rand(one(UInt64):UInt64(npools), length(vs))
    for i in eachindex(vs)
        view = vs[i]
        newviews[i] = TaggedSequencingView(seqid(view), first(view), last(view), tags[i])
    end
    return newviews
end

## The fragment transformation

struct BreakpointFragmenter{T<:AbstractSequencingView}
    seq::T
    nbreak::Int
end

Base.IteratorSize(::Type{BreakpointFragmenter{T}}) where {T} = Base.HasLength()
Base.IteratorEltype(::Type{BreakpointFragmenter{T}}) where {T} = Base.HasEltype()
Base.eltype(::Type{BreakpointFragmenter{T}}) where T = T
Base.length(it::BreakpointFragmenter) = it.nbreak + 1

function Base.iterate(it::BreakpointFragmenter{T}) where {T}
    breakpoints = sort!(sample_values(1:length(it.seq), it.nbreak))
    return iterate(it, (1, 0, breakpoints))
end

function Base.iterate(it::BreakpointFragmenter{T}, state::Tuple{Int,Int,Vector{Int}}) where {T}
    i, lastbp, breakpoints = state
    lastbp == lastindex(it.seq) && return nothing
    newstart = lastbp + 1
    newstop = i > it.nbreak ? lastindex(it.seq) : breakpoints[i]
    return subview(it.seq, newstart, newstop), (i + 1, newstop, breakpoints)
end

fragments(seq::AbstractSequencingView, nbreak::Int) = BreakpointFragmenter{typeof(seq)}(seq, nbreak)

function fragment(vs::Vector{T}, meansize::Int) where {T<:AbstractSequencingView}
    newviews = Vector{T}()
    for view in vs
        nbreak = div(length(view), meansize) - 1
        for f in fragments(view, nbreak)
            push!(newviews, f)
        end
    end
    return newviews
end

## End-taking transformations

function take_paired_ends(vs::Vector{T}, flen::Int, rlen::Int) where {T<:AbstractSequencingView}
    i = 0
    newviews = Vector{T}(undef, length(vs) * 2)
    for view in vs
        vl = length(view)
        if vl > flen && vl > rlen # Guard against the sloppy fragmenter.
            fwdview, revview = take_paired_ends(view, flen, rlen)
            newviews[i += 1] = fwdview
            newviews[i += 1] = revview
        end
    end
    resize!(newviews, i)
    return newviews
end

function take_single_ends(vs::Vector{T}, len::Int) where {T<:AbstractSequencingView}
    newviews = Vector{T}(undef, length(vs))
    i = 0
    for view in vs
        vl = length(view)
        if vl > len
            newviews[i += 1] = take_single_end(view, len)
        end
    end
    resize!(newviews, i)
    return newviews
end

function take_single_ends(vs::Vector{T}, len::Nothing) where {T<:AbstractSequencingView}
    newviews = Vector{T}(undef, length(vs))
    for (i, view) in enumerate(vs)
        newviews[i] = take_single_end(view)
    end
    return newviews
end

## Flipping transformation
#=
function flip_views!(vs::Vector{<:AbstractSequencingView})
    willflip = rand(Bool, length(vs))
    @inbounds for (i, view) in enumerate(vs)
        if willflip[i]
            vs[i] = reverse(vs[i])
        end
    end
    return vs
end

flip_views(vs::Vector{<:AbstractSequencingView}) = flip_views!(deepcopy(vs))
=#
# Summarizing functions
# ---------------------

function summarize_lengths(vs::Vector{<:AbstractSequencingView})::Tuple{Int, Int, Int}
    lens = length.(vs)
    return maximum(lens), div(sum(lens), length(vs)), minimum(lens)
end

function summarize_tags(vs::Vector{TaggedSequencingView})
    tagdict = Dict{UInt64, Int}()
    for v in vs
        t = tag(v)
        if t > 0
            tagdict[t] = get(tagdict, t, 0) + 1
        end
    end
    return tagdict
end


text:
module ServiceMesh

using ReinforcementLearningBase
using Random
using IntervalSets
using Distributions
using StableRNGs

# From RLBase
# export state, reward, reset!, action_space, state_space, is_terminated

# Others
export ServiceMeshEnv, JobParams, SinusArrival, SingleArrival, ConstantArrival, PoissonWrapper, StreamArrival, FlippingArrival

include("job.jl")
include("microservice.jl")
include("servicemeshenv.jl")

end # module

ext: jl

max_stars_repo_name: albheim/ServiceMeshControl

max_issues_repo_name: albheim/ServiceMeshControl

max_forks_repo_name: albheim/ServiceMeshControl


text:
using Symbolics, Test

@variables a b::Real z::Complex (Z::Complex)[1:10]

@testset "types" begin
    @test a isa Num
    @test b isa Num
    @test eltype(Z) <: Complex{Real}

    for x in [z, Z[1], z+a, z*a, z^2, z/z] # z/z is sus
        @test x isa Complex{Num}
        @test real(x) isa Num
        @test imag(x) isa Num
        @test conj(x) isa Complex{Num}
    end
end

@testset "repr" begin
    @test repr(z) == "z"
    @test repr(a + b*im) == "a + b*im"
end


text:
var a = "global";
{
  fun showA() {
    print a;
  }

  showA();
  var a = "block";
  showA();
}


text:
using DataFrames
using HDF5

function main()
    table_path = ARGS[1]
    id = h5read(table_path,"id")
    x = h5read(table_path,"x")
    y = h5read(table_path,"y")

    id1 = h5read(table_path,"id")
    x1 = h5read(table_path,"x")
    y1 = h5read(table_path,"y")

    df1 = DataFrame()
    df1[:id] = id
    df1[:x] = x
    df1[:y] = y

    df2 = DataFrame()
    df2[:id] = id1
    df2[:x1] = x1
    df2[:y1] = y1

    tic()
    join(df1,df2,on=:id)
    t = toq()

    println("DataFrames.jl join took: ", t)
end

main()

max_stars_repo_name: ehsantn/data-analytics-benchmarks

max_issues_repo_name: ehsantn/data-analytics-benchmarks

max_forks_repo_name: ehsantn/data-analytics-benchmarks


text:
function add1!(x,y)
    x[y] .+= 1
end

function devec_add1!(x,y)
	for i=1:length(y)
        x[y[i]] += 1
    end
end

function devec_add1_logical!(x,y)
	for i=1:length(y)
        if y[i]
            x[i] += 1
        end
    end
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
"""
複数の変数に対して，任意の q 値に対するパーセンタイル値を求めるには，以下のような関数を定義すればよい。

missing を除去して計算される。
"""

function quantile2(df, qs=[0, 0.25, 0.5, 0.75, 1])
    Names =["q$(Int(100q))" for q in qs]
    Names = replace(Names, "q0" => "min", "q50" => "median", "q100" => "max")
    qtile = map(x -> quantile(skipmissing(x), qs), eachcol(df))
    df2 = DataFrame(Matrix(reshape(vcat(qtile...), length(qs), :)'), Names)
    n = map(x -> length(filter(!ismissing, x)), eachcol(df))
    insertcols!(df2, 1, :variable => names(df), :n => n)
    df2
end

"""
使用例

julia> using RDatasets, Statistics

julia> iris = dataset("datasets", "iris");

julia> quantile2(iris[:, 1:4])
4×7 DataFrame
 Row │ variable     n      min      q25      median   q75      max
     │ String       Int64  Float64  Float64  Float64  Float64  Float64
─────┼─────────────────────────────────────────────────────────────────
   1 │ SepalLength    150      4.3      5.1     5.8       6.4      7.9
   2 │ SepalWidth     150      2.0      2.8     3.0       3.3      4.4
   3 │ PetalLength    150      1.0      1.6     4.35      5.1      6.9
   4 │ PetalWidth     150      0.1      0.3     1.3       1.8      2.5
"""

max_stars_repo_name: r-de-r/stats

max_issues_repo_name: r-de-r/stats

max_forks_repo_name: r-de-r/stats


text:
mutable struct SimpleExplainer
    model
    vp
    nb
    mask

    function SimpleExplainer(g, mg, vertex2array, frozen_model)
        vp = reduce(hcat, [vertex2array(mg, i) for i in vertices(g)])
        nvp = ArrayNode(Flux.param(vp))
        new(frozen_model, vp, graph2bag(TreeNode((nvp, BagNode(nvp,  neighbors_bags(g))))), Flux.param(ones(Float32, 1, nobs(vp)) .* 3))
    end
end

Flux.@treelike(SimpleExplainer)

function (se::SimpleExplainer)()
    nvp = ArrayNode(se.vp .* σ.(se.mask))
    g = graph2bag(TreeNode((nvp, BagNode(nvp,  se.nb))))
    se.model(g)
end


text:
module Metrics
    using ..Representation
    function calc_population_fitness!(population::Representation.Population,
                                    true_values::Array{Float64, 1},
                                    error_function, extra_params=nothing)
        for memeber in population.members
            pred_values = [kernels_to_values(x, memeber.gauss_kernels) for x in true_values]
            if extra_params != nothing && error_function == weighted_mean_abs_error
                (error, hit_ratio) = error_function(true_values,
                                        pred_values,
                                        extra_params["threshold"],
                                        extra_params["lower_weight"],
                                        extra_params["upper_weight"])
            else
                (error, hit_ratio) = error_function(true_values, pred_values)
            end
            memeber.fitness = error
            memeber.hit_ratio = hit_ratio
        end
    end

    function mean_sqr_error(true_values::Array{Float64, 1}, pred_values::Array{Float64, 1})
        error = sum((true_values - pred_values).^2)/length(true_values)
        hit_ratio = sum(true_values .== pred_values)
        return error, hit_ratio
    end

    function mean_absolute_error(true_values::Array{Float64, 1}, pred_values::Array{Float64, 1})
        error = sum(abs.(true_values - pred_values))/length(true_values)
        hit_ratio = sum(true_values .== pred_values)
        return error, hit_ratio
    end

    function weighted_mean_abs_error(true_values::Array{Float64, 1},
                                    pred_values::Array{Float64, 1},
                                    threshold::Float64,
                                    lower_weight::Float64,
                                    upper_weight::Float64)

        abs_error = abs.(true_values - pred_values)
        mask = [apply_threshold(value, threshold,
                        lower_weight, upper_weight) for value in abs_error]
        hit_ratio = sum(mask .== 1.0)
        return sum(abs_error.*mask)/length(true_values), hit_ratio

    end

    function apply_threshold(value::Float64, threshold::Float64,
                                lower_weight::Float64, upper_weight::Float64)
        return value <= threshold ? lower_weight : upper_weight
    end

    function kernels_to_values(x::Float64,
                                kernels::Array{Tuple{Float64, Float64, Float64}, 1})
        sum = Float32(0.0)
        for kernel in kernels
            (weight, c, γ) = kernel
            sum += weight * exp(-γ*(c - x)^2)
        end
        return sum
    end

end


text:
"""
    vcf2dataframe(vcf_data_array)

Convert a vcf data array from the readbigvcf.jl function
and return a standard data frame.

Requires: DataFrames
"""
function vcf2dataframe(vcf_data_array)

vcf = vcf_data_array
hapmat = []
header = []
dataline = []
lup = Dict()

for i in eachindex(vcf)

    if occursin(r"^##", vcf[i])

        continue

    elseif occursin(r"^#CHROM", vcf[i])

        ids = split(vcf[i], r"\t")
        ids[1] = replace(ids[1], "#" => "")

        push!(header, ids[1], ids[2])

        for a = 10:length(ids)
            push!(header, ids[a] * "p", ids[a] * "q")
        end

        hapmat = [hapmat; header]

    else

        g = split(vcf[i], '\t')

        push!(dataline, g[1], g[2])

        for i = 10:length(g)
            a = split(g[i], r"[\||\:|\/]")
            push!(dataline, a[1], a[2])
        end

        dataline = replace.(dataline, "." => "-1") #missing to -1

        addbases = false #routine to add bases calls

        if addbases == true
            alts = enumerate(split(g[5], ","))

            if (length(g[4]) +  length(g[5])) > 2 #indels remain numeric
                continue
            end

            dataline = replace.(dataline, r"0" => g[4])

            for (i,v) = alts
                dataline = replace.(dataline, Regex("$i") => v)
            end
        end

        hapmat = [hapmat dataline]

        dataline = []
    end
end

p = permutedims(hapmat)
pm = p[2:end, :]
h = p[1, :]
t = DataFrame(pm, :auto)
tt = rename!(t, Symbol.(h))

return tt

end

max_stars_repo_name: ScottWatkins/STRdiff

max_issues_repo_name: ScottWatkins/STRdiff

max_forks_repo_name: ScottWatkins/STRdiff


text:
module test_octo_schema

using Test # @test_throws
using Octo.Adapters.SQL # Schema.changeset Repo
using Octo.Schema # InvalidChangesetError validate_length

struct User
end
Schema.model(User, table_name="users")

Schema.changeset(User) do model
    validate_length(model, :username, min=6)
end

Repo.connect(adapter=Octo.Adapters.SQL)

Repo.insert!(User, (username="abcdef",))
@test_throws InvalidChangesetError Repo.insert!(User, (username="a",))

Schema.model(User, table_name="users", primary_key=[:pk1, :pk2])
Repo.get(User, (pk1=1, pk2=2))

Schema.model(User, table_name="users")
Schema.model(User, table_name="users", primary_key=nothing)
Schema.model(User, table_name="users", primary_key=:pk1)
Schema.model(User, table_name="users", primary_key=[:pk1, :pk2])
Schema.model(User, table_name="users", primary_key="pk1")
Schema.model(User, table_name="users", primary_key=["pk1", "pk2"])

end # module test_octo_schema


text:
struct MultiDiscreteSpace{N} <:AbstractSpace
    counts::Array{Int, N}
end

size(s::MultiDiscreteSpace) = size(s.counts)

"To compat with Python, here we start with 0"
sample(s::MultiDiscreteSpace) = map(x -> rand(0:x-1), s.counts)
occursin(x::Int, s::MultiDiscreteSpace) = all(e -> 0 ≤ x < e, s.counts)
occursin(xs::Array{Int}, s::MultiDiscreteSpace) = size(s) == size(xs) &&
    all(map((e, x) -> 0 ≤ x < e , s.counts, xs))
==(x::MultiDiscreteSpace, y::MultiDiscreteSpace) = x.counts == y.counts


text:
# AUTO GENERATED FILE - DO NOT EDIT

export datepicker

"""
    datepicker(;kwargs...)

A DatePicker component.
DatePicker
Keyword arguments:
- `allowInput` (Bool; optional): flatpickr prop passthrough. Allows the user to enter a date directly
into the input field
- `className` (String; optional): The CSS class names.
- `id` (String; required): The ID of this component, used to identify dash components
in callbacks. The ID needs to be unique across all of the
components in an app.
- `dateFormat` (String; optional): The date format.
- `light` (Bool; optional): True to use the light version.
- `locale` (String; optional): The language locale used to format the days of the week, months, and numbers.
The full list of supported locales can be found here https://github.com/flatpickr/flatpickr/tree/master/src/l10n
- `maxDate` (String; optional): The maximum date that a user can pick to.
- `minDate` (String; optional): The minimum date that a user can start picking from.
- `short` (Bool; optional): True to use the short version.
- `value` (String | Array of Strings; optional): The value of the date value provided to flatpickr
- `placeholder` (String; required): Placeholder
- `labelText` (String; required): Label text
"""
function datepicker(; kwargs...)
        available_props = Symbol[:allowInput, :className, :id, :dateFormat, :light, :locale, :maxDate, :minDate, :short, :value, :placeholder, :labelText]
        wild_props = Symbol[]
        return Component("datepicker", "DatePicker", "dash_carbon_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: Matheus-Rangel/dash-carbon-components

max_issues_repo_name: Matheus-Rangel/dash-carbon-components

max_forks_repo_name: Matheus-Rangel/dash-carbon-components


text:
using Base.Cartesian

"""
deleteat{N,T,i}(x::NTuple{N,T}, ::Type{Val{i}})

Return copy of x with ith entry ommited.
"""
@generated function deleteat{N,T,i}(x::NTuple{N,T}, ::Type{Val{i}})
    (1 <= i <= N) || throw(MethodError(drop_index, (x,Val{i})))
    args = [:(x[$j]) for j in deleteat!([1:N...], i)]
    Expr(:tuple, args...)
end

"""
deleteat{N,T}(x::NTuple{N,T}, i::Int)
"""
@generated deleteat{N,T}(x::NTuple{N,T}, i::Int) = quote
    # would be nice to eliminate boundscheck
    (1 <= i <= N) || throw(BoundsError(x,i))
    @nif $(N) d->(i==d) d-> deleteat(x, Val{d})
end

"""
best_arg, best_val = argmax(f, iter)

Computes the first element of iterator, on which f takes its maximum.
"""
function argmax(f, iter)
    state = start(iter)
    best_arg, state = next(iter, state)
    best_val = f(best_arg)
    while !done(iter,state)
        arg, state = next(iter,state)
        val = f(arg)
        if val > best_val
            best_val = val
            best_arg = arg
        end
    end
    best_arg, best_val
end


text:
module nng

const LIB ="libnng.so.1"

mutable struct nng_socket
    id::UInt32
end

mutable struct nng_listener
    id::UInt32
end

mutable struct nng_dialer
    id::UInt32
end

mutable struct nng_aio end

const NNG_OPT_SUB_SUBSCRIBE = "sub:subscribe"
const NNG_OPT_SUB_UNSUBSCRIBE = "sub:unsubscribe"

ERROR_CODES = Dict(
1  => "NNG_EINTR",
2  => "NNG_ENOMEM",
3  => "NNG_EINVAL",
4  => "NNG_EBUSY",
5  => "NNG_ETIMEDOUT",
6  => "NNG_ECONNREFUSED",
7  => "NNG_ECLOSED",
8  => "NNG_EAGAIN",
9  => "NNG_ENOTSUP",
10 => "NNG_EADDRINUSE",
11 => "NNG_ESTATE",
12 => "NNG_ENOENT",
13 => "NNG_EPROTO",
14 => "NNG_EUNREACHABLE",
15 => "NNG_EADDRINVAL",
16 => "NNG_EPERM",
17 => "NNG_EMSGSIZE",
18 => "NNG_ECONNABORTED",
19 => "NNG_ECONNRESET",
20 => "NNG_ECANCELED",
21 => "NNG_ENOFILES",
22 => "NNG_ENOSPC",
23 => "NNG_EEXIST",
24 => "NNG_EREADONLY",
25 => "NNG_EWRITEONLY",
26 => "NNG_ECRYPTO",
27 => "NNG_EPEERAUTH",
28 => "NNG_ENOARG",
29 => "NNG_EAMBIGUOUS",
30 => "NNG_EBADTYPE",
31 => "NNG_ECONNSHUT",
1000 => "NNG_EINTERNAL"
)

@enum SOCKET_TYPES begin
    PUSH0
    PULL0
    PUB0
    SUB0
    PAIR0
    REQ0
    REP0
    SURVEYOR0
    RESPONDENT0
    BUS0
end

function _handle_err(err:: Int32)::Int32
    if err != 0
        throw(error("NNG error: $(ERROR_CODES[err])"))
    end
    return err
end

"""
Low level function to get the socket ID.

In many occasions, this returns -1 (Error), but the socket works anyway
"""
_nng_socket_id(socket::nng_socket) =
    ccall((:nng_socket_id, LIB), Cint, (Ref{nng_socket},), Ref(socket))

"""
Low level function to set a string as an option
"""
_nng_setopt(socket::nng_socket, option::String, value::String) =
    ccall((:nng_setopt, LIB), Cint, (nng_socket, Cstring, Cstring, Csize_t), socket, option, value, length(value))

"""
Low level listen
"""
_nng_listen(socket::nng_socket, url::String) = ccall(
    (:nng_listen, LIB),
    Cint,
    (nng_socket, Cstring, Ref{nng_listener}, Cint),
    socket,
    url,
    Ref(nng_listener(0)),
    0,
)
"""
Low level dial
"""
_nng_dial(socket::nng_socket, url::String) = ccall(
    (:nng_dial, LIB),
    Cint,
    (nng_socket, Cstring, Ref{nng_dialer}, Cint),
    socket,
    url,
    Ref(nng_dialer(0)),
    0,
)

"""
Low level PULL
"""
_nng_pull0_open(socket::nng_socket) =
    ccall((:nng_pull0_open, LIB), Cint, (Ref{nng_socket},), Ref(socket))

"""
Low level PUSH
"""
_nng_push0_open(socket::nng_socket) =
    ccall((:nng_push0_open, LIB), Cint, (Ref{nng_socket},), Ref(socket))

"""
Low level PUB
"""
_nng_pub0_open(socket::nng_socket) =
    ccall((:nng_pub0_open, LIB), Cint, (Ref{nng_socket},), Ref(socket))

"""
Low level SUB
"""
_nng_sub0_open(socket::nng_socket) =
    ccall((:nng_sub0_open, LIB), Cint, (Ref{nng_socket},), Ref(socket))

"""
Low level PAIR
"""
_nng_pair0_open(socket::nng_socket) =
    ccall((:nng_pair0_open, LIB), Cint, (Ref{nng_socket},), Ref(socket))

"""
Low level REQ
"""
_nng_req0_open(socket::nng_socket) =
ccall((:nng_req0_open, LIB), Cint, (Ref{nng_socket},), Ref(socket))

"""
Low level REP
"""
_nng_rep0_open(socket::nng_socket) =
ccall((:nng_rep0_open, LIB), Cint, (Ref{nng_socket},), Ref(socket))

"""
Low level SURVEYOR
"""
_nng_surveyor0_open(socket::nng_socket) =
ccall((:nng_surveyor0_open, LIB), Cint, (Ref{nng_socket},), Ref(socket))

"""
Low level RESPONDENT
"""
_nng_respondent0_open(socket::nng_socket) =
ccall((:nng_respondent0_open, LIB), Cint, (Ref{nng_socket},), Ref(socket))

"""
Low level BUS
"""
_nng_bus0_open(socket::nng_socket) =
ccall((:nng_bus0_open, LIB), Cint, (Ref{nng_socket},), Ref(socket))

"""
Low level CLOSE
"""
_nng_close(socket::nng_socket) = ccall((:nng_close, LIB), Cint, (nng_socket,), socket)

"""
Low level SEND.
"""
function _nng_send(socket::nng_socket, message::AbstractString)
    return ccall(
        (:nng_send, LIB),
        Cint,
        (nng_socket, Ptr{UInt8}, Csize_t, Cint),
        socket,
        pointer(message),
        sizeof(message) + 1,
        0,
    )
end

"""
Low level RECV.
"""
function _nng_recv(socket::nng_socket)
    buf = Vector{Ptr{UInt8}}(undef, 1)
    size = Csize_t(0)
    rval = ccall(
        (:nng_recv, LIB),
        Cint,
        (nng_socket, Ptr{UInt8}, Ref{Csize_t}, Cint),
        socket,
        pointer(buf),
        size,
        1,
    )

    msg = unsafe_string(buf[1])
    ccall((:nng_free, LIB), Cvoid, (Ptr{Cvoid}, Csize_t), buf[1], sizeof(msg))

    return msg
end

function_mapping = Dict{SOCKET_TYPES, Function}(
PULL0 => _nng_pull0_open,
PUSH0 => _nng_push0_open,
PUB0 => _nng_pub0_open,
SUB0 => _nng_sub0_open,
PAIR0 => _nng_pair0_open,
REQ0 => _nng_req0_open,
REP0 => _nng_rep0_open,
SURVEYOR0 => _nng_surveyor0_open,
RESPONDENT0 => _nng_respondent0_open,
BUS0 => _nng_bus0_open
)

"""
NNG listen to a socket for other connections
"""
function listen(addr::String, proto::SOCKET_TYPES)::nng_socket
    if haskey(function_mapping, proto) == false
        throw(error("Not one of accepted socket types"))
    end
    socket = nng_socket(0)
    err_val = _handle_err(function_mapping[proto](socket))
    err_val = _handle_err(_nng_listen(socket, addr))
    return socket
end

"""
NNG dial to a socket
"""
function dial(addr::String, proto::SOCKET_TYPES)::nng_socket
    if haskey(function_mapping, proto) == false
        throw(error("Not one of accepted socket types"))
    end
    socket = nng_socket(0)
    err_val = _handle_err(function_mapping[proto](socket))
    err_val = _handle_err(_nng_dial(socket, addr))
    return socket
end

"""
NNG dial to a socket and subscribe to a topic
"""
function dial(addr::String, proto::SOCKET_TYPES, topic::String)::nng_socket
    if haskey(function_mapping, proto) == false
        throw(error("Not one of accepted socket types"))
    end
    socket = nng_socket(0)
    err_val = _handle_err(function_mapping[proto](socket))
    err_val = _handle_err(_nng_setopt(socket, NNG_OPT_SUB_SUBSCRIBE, topic))
    err_val = _handle_err(_nng_dial(socket, addr))
    return socket
end

"""
NNG dial with a socket
"""
function dial(socket::nng_socket, addr::String)::nng_socket
    err_val = _handle_err(_nng_dial(socket, addr))
    return socket
end

"""
NNG close a socket
"""
function close(socket::nng_socket)::Int32
    return _handle_err(_nng_close(socket))
end

"""
NNG subscribe to a topic
"""
function subscribe(socket::nng_socket, topic::String)::Int32
    return _handle_err(_nng_setopt(socket, NNG_OPT_SUB_SUBSCRIBE, topic))
end

"""
NNG unsubscribe to a topic
"""
function unsubscribe(socket::nng_socket, topic::String)::Int32
    return _handle_err(_nng_setopt(socket, NNG_OPT_SUB_UNSUBSCRIBE, topic))
end

"""
NNG send an abstract string from a socket
"""
function send(socket::nng_socket, msg::AbstractString)::Int32
    return _handle_err(_nng_send(socket, msg))
end

"""
NNG receive an abstract string to a socket
"""
function recv(socket::nng_socket)::AbstractString
    return _nng_recv(socket)
end

end


text:
"""
This is the module where the actual interaction with typescript happens.

It's mostly wrapped up in the Semagrams struct; see the documentation there.
"""
module UI

export Semagram, get_acset, serve_semagram, save, load, close!, export_html, export_pdf

using JSExpr
using WebIO
using Mux
import JSON
import DefaultApplication
using ..Schema
using ..Muesli
using ..Data
using Catlab.CSetDataStructures

"""
Basically the "control wires" for a semagram.

Contains:
- `divid`, the DOM ID of the div where the semagram should be anchored.
- `scope`, the WebIO scope in which the semagram lives
- `handle`, the WebIO handle in which data is passed *back* from the semagram.
   Note: there is currently no way of passing data *to* the semagram, but when
   we add this, it should be a parameter here.
- `ws`, the schema of the semagram. This gets passed to typescript, and also
   tells us how to interpret the data that we are getting back from typescript.
"""
struct Semagram{T <: ACSet}
  scope::Scope
  receiving::Observable{Dict{String,Any}}
  sending::Observable{Dict{String,Any}}
  exported::Observable{String}
  function Semagram{T}(ls::LocatedSemagramData) where {T <: ACSet}
    try
      T()
    catch e
      error("""
        error: all type parameters to the acset type must be fully instantiated.

        Right: Semagram{WeightedGraph{Int}}(WeightedGraphSema)
        Wrong: Semagram{WeightedGraph}(WeightedGraphSema)
      """)
    end
    deps = [
      "semagrams" => joinpath(@__DIR__, "..", "deps", "bundles", "app.bundle.js")
    ]
    scope = Scope(imports=deps, dom=dom"div"())
    ls_json = to_json(ls)
    receiving = Observable(scope, "receiving", ls_json)
    sending = Observable(scope, "sending", ls_json)
    exported = Observable(scope, "exported", "")
    onjs(sending, @js function (newls)
           this.state.resetWith(newls)
         end)
    on((newls) -> receiving[] = newls, sending)
    mountfn = @js function ()
      @var semagrams = System.registry.get(System.resolveSync("semagrams"))
      @var scopeobj = this
      setTimeout(() ->
        semagrams.main(
          $receiving[],
          scopeobj,
          (x) -> $receiving[] = x,
          (x) -> $exported[] = x), 20)
    end
    onmount(scope, mountfn)
    new{T}(scope, receiving, sending, exported)
  end
end

function Semagram{T}(s::SemagramSchema) where {T <: ACSet}
  ls = LocatedSemagramData(s)
  Semagram{T}(ls)
end

@WebIO.register_renderable(Semagram) do sg
    return WebIO.render(sg.scope)
end

function get_acset(sema::Semagram{T}) where {T}
  to_acset(save(sema), T)
end

function save(sema::Semagram)
  from_json(sema.receiving[], LocatedSemagramData)
end

function save(sema::Semagram, fn::String)
  open(fn, "w") do f
    write(f, JSON.json(sema.receiving[]))
  end
end

function load(sema::Semagram, fn::String)
  load(sema, load(fn))
end

function load(sema::Semagram, ls::LocatedSemagramData)
  sema.sending[] = to_json(ls)
end

function load(fn::String)
  v = Dict()
  open(fn, "r") do f
    v = JSON.parse(read(f, String))
  end
  from_json(v, LocatedSemagramData)
end

# NOTE: Not sure how to catch all the errors associated with WebSockets

function close!(sema_task::Task)
  try
    schedule(sema_task,ErrorException("stop"),error=true)
  catch e
    @warn "Semagram server closed"
  end
end

HTML_WRAPPING = """
<!doctype html>

<html>
<head>
  <meta charset="UTF-8">
</head>
<body>
  @SVG@
</body>
</html>
"""

function export_html(sema::Semagram;
                     standalone::Union{Bool,Nothing}=nothing,
                     file::Union{String,Nothing}=nothing)
  # standalone defaults to true if we are writing to a file,
  # false otherwise
  standalone = standalone == nothing ? file != nothing : standalone
  body = sema.exported[]
  html = if standalone
    replace(HTML_WRAPPING, "@SVG@" => body)
  else
    body
  end
  if file == nothing
    HTML(html)
  else
    open(file, "w") do f
      write(f, html)
    end;
  end
end

function export_pdf(sema::Semagram, file::String)
  htmlfile = tempname() * ".html"
  export_html(sema; file=htmlfile)
  run(`wkhtmltopdf -s A6 -O Landscape $htmlfile $file`)
end

"""
This serves a semagram in a standalone page (as opposed to having it displayed
inline in a Jupyter notebook).

TODO: I currently don't know how to kill a task spawned by this...
"""
function serve_semagram(s::Semagram; port=8000)
  task = webio_serve(page("/", req -> node(:div, s)), port)
  sleep(1)
  DefaultApplication.open("http://localhost:$port")
  task
end

end # module UI


text:
module XVGReader

    include("read.jl")
    export read_xvg,f

end


text:
differentiate(f::Fun{OscLaurent{DD,RR}}) where {DD,RR} = (f.space.exp ≈ 0.0) ?  Fun(f.space,1/(domain(f.space).L)*mobiusdiff(f.coefficients)) : Fun(f.space,vcat(1.0im*f.space.exp*f.coefficients,[0.0,0.0]) + 1/(domain(f.space).L)*mobiusdiff(f.coefficients))

function Lag(n::Int64,x::Float64) # evaluate Laguerre polynomials at x, α = 1
  c = zeros(Float64,n+2)
  c[2] = exp.(-x/2)
  k = 0
  for i = 3:n+2 #k = 1 is the zeroth order, gives (k-1)th order
    c[i] = (2*k+2-x)*c[i-1]/(k+1) - c[i-2]
    k = k + 1
  end
  return c
end

Base.sum(f::Fun{T}) where T <: ArraySpace = map(sum,f)

# Old version, still here for testing
function sum2(f::Fun{OscLaurent{DD,RR}}) where {DD,RR}
    if space(f).exp ≈ 0.
        data = f.coefficients[2:2:end]
        data = data.*(j*(-1)^j for j in 1:length(data))
        s = sum(data)
        data = f.coefficients[3:2:end]
        data = data.*(j*(-1)^j for j in 1:length(data))
        return -2*pi*(s + sum(data))*domain(f).L
    else
        m = convert(Int64,floor(length(coefficients(f))/2))
        ex = space(f).exp
        β = domain(f).L
        z = 2*β*abs(ex)
        c = Lag(m-1,z)[2:end]
        #cc = v.coefs
        #if abs(v.exp) == 0
        #  c = c.*(N(2*m)[2:length(c)+1])
        #end # if v.exp > 0, only negative terms contribute
        if ex > 0.
          data = f.coefficients[2:2:end]
        elseif ex < 0.
            data = f.coefficients[3:2:end]
        end
        data = data.*((-1)^j for j in 1:length(data))
        return -4*pi*β*sum(data.*c[1:length(data)])
    end
end

function Base.sum(f::Fun{OscLaurent{DD,RR}}) where {DD,RR <: BigFloat}
    if space(f).exp ≈ 0
        s = BigFloat("0")*1im
        o = 1
        j = 0
        for i = 2:2:length(f.coefficients)
            o *= -1
            j += 1
            @inbounds s += f.coefficients[i]*j*o
        end
        o = 1
        j = 0
        for i = 3:2:length(f.coefficients)
            o *= -1
            j += 1
            @inbounds s += f.coefficients[i]*j*o
        end
        return -2*big(pi)*s*BigFloat(domain(f).L)
    else
        #m = convert(Int64,floor(length(coefficients(f))/2))
        ex = space(f).exp
        β = BigFloat(domain(f).L)
        z = 2*β*abs(ex)
        c0 = 0
        c1 = exp(-z/2)
        #c = Lag(m-1,z)[2:end]
        #cc = v.coefs
        #if abs(v.exp) == 0
        #  c = c.*(N(2*m)[2:length(c)+1])
        #end # if v.exp > 0, only negative terms contribute
        if real(ex) > 0
          #data = f.coefficients[2:2:end]
          j = 2
        else
            #data = f.coefficients[3:2:end]
            j = 3
        end
        s = BigFloat("0")*1im
        o = 1
        k = 0
        for i = j:2:length(f.coefficients)
            o *= -1
            c2 = (2*k +2 -z)*c1/(k+1) - c0
            @inbounds s += o*c1*f.coefficients[i]
            c0 = c1
            c1 = c2
            k += 1
        end

        #data = data.*((-1)^j for j in 1:length(data))
        return -4*big(pi)*β*s
    end
end

function Base.sum(f::Fun{OscLaurent{DD,RR}}) where {DD,RR}
    if space(f).exp ≈ 0.
        s = 0.0im
        o = 1
        j = 0
        for i = 2:2:length(f.coefficients)
            o *= -1
            j += 1
            @inbounds s += f.coefficients[i]*j*o
        end
        o = 1
        j = 0
        for i = 3:2:length(f.coefficients)
            o *= -1
            j += 1
            @inbounds s += f.coefficients[i]*j*o
        end
        return -2*pi*s*domain(f).L
    else
        #m = convert(Int64,floor(length(coefficients(f))/2))
        ex = space(f).exp
        β = domain(f).L
        z = 2*β*abs(ex)
        c0 = 0.
        c1 = exp(-z/2)
        #c = Lag(m-1,z)[2:end]
        #cc = v.coefs
        #if abs(v.exp) == 0
        #  c = c.*(N(2*m)[2:length(c)+1])
        #end # if v.exp > 0, only negative terms contribute
        if real(ex) > 0.
          #data = f.coefficients[2:2:end]
          j = 2
        else
            #data = f.coefficients[3:2:end]
            j = 3
        end
        s = 0.0im
        o = 1
        k = 0
        for i = j:2:length(f.coefficients)
            o *= -1
            c2 = (2*k +2 -z)*c1/(k+1) - c0
            @inbounds s += o*c1*f.coefficients[i]
            c0 = c1
            c1 = c2
            k += 1
        end

        #data = data.*((-1)^j for j in 1:length(data))
        return -4*pi*β*s
    end
end


text:
#=
-------------------------------------------------------------------------------------------------------------
This file runs the 1D Wave Equation using a unified CG/DG method using the AGGP storage
introduced in F.X. Giraldo's Introduction to Element-based Galerkin Methods using
Tensor-Product Bases: Analysis, Algorithms, and Applications.
This code constructs global matrices, which although not the most efficient implementation,
allows for the analysis of the resulting global matrices.

Written by F.X. Giraldo on July 5, 2021.
Department of Applied Mathematics
Naval Postgraduate School
Monterey, CA 93943

For time-integration, we use SSP methods of order 2 and 3; and LSRK45 (4th order, 5-stage).

The interpolation points used are the following:
ipoints=1: Lobatto
ipoints=2: Legendre
ipoints=3: Chebyshev
ipoints=4: Equi-spaced

The integration points used are the following:
qpoints=1: Lobatto
qpoints=2: Legendre
-------------------------------------------------------------------------------------------------------------
=#

using Plots

include("QuadraturePoints.jl")
include("element_matrices.jl")
include("create_grid.jl")
include("global_matrices.jl")
include("exact_solution.jl")
include("compute_norms.jl")
include("time_integration.jl")

#Some Constants
DFloat = Float64
machine_zero=eps(DFloat)

#{{{ Main
function main()

    #Fix Interpolation and Integration Order
    N=4
    Q=N+1
    Np=N+1
    Nq=Q+1
    ipoints=1
    qpoints=1
    space_method="DG"
    #space_method="DG"
    #ti_type="SSP"
    ti_type="LSRK"
    case=1 #1=exponential, 2=square wave
    stages=3 #1,2, or 3 only
    if (ti_type == "LSRK")
        stages=5
    end
    Ne=20
    Npoin_dg=Ne*Np
    Npoin_cg=Ne*N + 1
    dt=DFloat(0.0001)
    time_final=DFloat(2.0)
    plot_movie=true
    Courant_max=0.5
    ntime = ceil(Int64, time_final / dt)
    @show(N,Q,Ne,Npoin_cg,Npoin_dg,dt,ntime,case,ti_type,stages)

    #Select Interpolation Points
    if ipoints == 1
        (ξ,ω) = QuadraturePoints.lobatto_gauss(Np,DFloat) #Lobatto
    elseif ipoints == 2
        (ξ,ω) = QuadraturePoints.legendre_gauss(Np,DFloat) #Legendre
    elseif ipoints == 3
        (ξ,ω) = QuadraturePoints.chebyshev_gauss(Np,DFloat) #Chebyshev
    elseif ipoints == 4
        (ξ,ω) = QuadraturePoints.equispaced_points(Np,DFloat) #Equi-spaced
    end

    #Select Quadrature Points
    if qpoints == 1
        (ξq,ωq) = QuadraturePoints.lobatto_gauss(Nq,DFloat) #Lobatto
    elseif qpoints == 2
        (ξq,ωq) = QuadraturePoints.legendre_gauss(Nq,DFloat) #Legendre
    end

    #Construct Lagrange Polynomials
    (ψ,dψ) = QuadraturePoints.lagrange_basis(Np,Nq,ξ,ξq,DFloat)

    #Construct Local Mass and Differentiation Matrices
    (Me,De,Dwe,Fe) = element_matrices(ψ,dψ,Np,Nq,ωq,DFloat)

    #Construct Grid and (i,e)->I map
    (coord_cg,coord_dg,intma_cg,intma_dg,periodicity_cg,periodicity_dg) = create_grid(Np,Npoin_cg,Npoin_dg,Ne,ξ,DFloat)

    #=
    @show(coord_cg, intma_cg)
    @show( size(intma_cg) )
    @show(coord_dg, intma_dg)
    @show( size(intma_dg) )
    =#

    #CG or DG?
    if (space_method == "CG")
        Npoin=Npoin_cg
        coord=copy(coord_cg)
        intma=copy(intma_cg)
        periodicity=copy(periodicity_cg)
    elseif (space_method == "DG")
        Npoin=Npoin_dg
        coord=copy(coord_dg)
        intma=copy(intma_dg)
        periodicity=copy(periodicity_dg)
    end
    @show(space_method,Npoin,Ne,N,Q)

    #Construct Global Mass and Differentiation Matrices
    (M,D,Dw,F) = global_matrices(Me,De,Dwe,Fe,intma,coord,Npoin,Ne,Np,periodicity,DFloat)
    #Dhat=-M\D
    Dhat=M\Dw
    Fhat=M\F

    #Compute Initial Condition
    time=DFloat(0.0)
    (qe,u) = exact_solution(coord,Npoin,time,case,DFloat)
    q0=copy(qe)
    dx_min=coord[2]-coord[1]

    #Recompute DT, NTIME
    dt=Courant_max*dx_min/u
    courant=u*dt/dx_min
    ntime = ceil(Int64, time_final / dt)
    println("Courant = ",courant," dt = ",dt," ntime = ",ntime)

    #Time Integration
    q0 = time_integration!(q0,u,coord,Dhat,Fhat,intma,periodicity,time,ntime,dt,stages,DFloat,ti_type,space_method,plot_movie)
    println(" itime = ",ntime," time = ",time," qmax = ",maximum(q0)," qmin = ",minimum(q0))

    #Compute L2 Norm
    (qe,u) = exact_solution(coord,Npoin,time,case,DFloat)
    norms = compute_norms(q0,qe,Npoin,DFloat)
    println(" norm.L1 = ",norms[1]," norm.L2 = ",norms[2]," norm.L∞ = ",norms[3])

    #Plot Interpolation
    qarray=zeros(DFloat,Npoin,2)
    qarray[:,1] .= q0
    qarray[:,2] .= qe
    plot_handle=plot(coord,qarray,xlabel="x",ylabel="q(x,t)",legend=true,lw=3,label=[space_method "Exact"],title=space_method)
    display(plot_handle)
#    savefig(plot_handle,"Project2.png")

    println("Done") #output

end
#}}} Main

#----------------------------------#
# Run the main function
#----------------------------------#
main()


text:
#Module for the kernel type
#

module KernelFunctions

export Kernel
export CreateKernelMatrix, CreateDiagonalKernelMatrix
export delta_kroenecker

type Kernel
    kernel_function::Function # Kernel function
    coeff::Float64 #Weight for the kernel
    derivative_kernel::Function # Derivative of the kernel function (used for hyperparameter optimization)
    param::Float64 #Hyperparameter for the kernel function (depends of the function)
    Nparams::Int64 #Number of hyperparameters
    compute::Function #General computational function
    compute_deriv::Function #General derivative function
    #Constructor
    function Kernel(kernel, coeff::Float64; params=0)
      this = new()
      this.coeff = coeff
      this.Nparams = 1
      if kernel=="rbf"
        this.kernel_function = rbf
        this.derivative_kernel = deriv_rbf
        this.param = params
      elseif kernel=="quadra"
        this.kernel_function = quadratic
        this.derivative_kernel = deriv_quadratic
        this.param = params
      elseif kernel=="linear"
        this.kernel_function = linear
        this.Nparams = 0
      elseif kernel=="laplace"
        this.kernel_function = laplace
        this.derivative_kernel = deriv_laplace
        this.param = params
      elseif kernel=="abel"
        this.kernel_function = abel
        this.derivative_kernel = deriv_abel
        this.param = params
      elseif kernel=="imk"
        this.kernel_function = imk
        this.derivative_kernel = deriv_imk
        this.param = params
      else
        error("Kernel function $(kernel_list[i]) not available, options are : rbf, linear, laplace, abel, imk")
      end
      if this.Nparams > 0
        this.compute = function(X1,X2)
            this.kernel_function(X1,X2,this.param)
          end
        this.compute_deriv = function(X1,X2)
            this.derivative_kernel(X1,X2,this.param)
          end
      else
        this.compute = function(X1,X2)
            this.kernel_function(X1,X2)
          end
      end
      return this
    end
end

function CreateKernelMatrix(X1,kernel_function;X2=0) #Create the kernel matrix from the training data or the correlation matrix between two set of data
  if X2 == 0
    ksize = size(X1,1)
    K = zeros(ksize,ksize)
    for i in 1:ksize
      for j in 1:i
        K[i,j] = kernel_function(X1[i,:],X1[j,:])
        if i != j
          K[j,i] = K[i,j]
        end
      end
    end
    return K
  else
    ksize1 = size(X1,1)
    ksize2 = size(X2,1)
    K = zeros(ksize1,ksize2)
    for i in 1:ksize1
      for j in 1:ksize2
        K[i,j] = kernel_function(X1[i,:],X2[j,:])
      end
    end
    return K
  end
end

function CreateDiagonalKernelMatrix(X,kernel_function;MatrixFormat=false)
  n = size(X,1)
  kermatrix = zeros(n)
  for i in 1:n
    kermatrix[i] = kernel_function(X[i,:],X[i,:])
  end
  if MatrixFormat
    return diagm(kermatrix)
  else
    return kermatrix
  end
end

function delta_kroenecker(X1::Array{Float64,1},X2::Array{Float64,1})
  return X1==X2 ? 1 : 0
end

#Gaussian (RBF) Kernel
function rbf(X1::Array{Float64,1},X2::Array{Float64,1},Θ)
  if X1 == X2
    return 1
  end
  exp(-(norm(X1-X2))^2/(Θ^2))
end

function deriv_rbf(X1::Array{Float64,1},X2::Array{Float64,1},Θ)
  a = norm(X1-X2)
  if a != 0
    return 2*a^2/(Θ^3)*exp(-a^2/(Θ^2))
  else
    return 0
  end
end

#Laplace Kernel
function laplace(X1::Array{Float64,1},X2::Array{Float64,1},Θ)
  if X1 == X2
    return 1
  end
  exp(-Θ*norm(X1-X2,2))
end

function deriv_laplace(X1::Array{Float64,1},X2::Array{Float64,1},Θ)
  a = norm(X1-X2,2)
  -a*exp(-Θ*a)
end

#Abel Kernel
function abel(X1::Array{Float64,1},X2::Array{Float64,1},Θ)
  if X1==X2
    return 1
  end
  exp(-Θ*norm(X1-X2,1))
end

function deriv_abel(X1::Array{Float64,1},X2::Array{Float64,1},Θ)
  a = norm(X1-X2,1)
  -a*exp(-Θ*a)
end

function imk(X1::Array{Float64,1},X2::Array{Float64,1},Θ)
  1/(sqrt(norm(X1-X2)^2+Θ))
end

function deriv_imk(X1::Array{Float64,1},X2::Array{Float64,1},Θ)
  -1/(2*(sqrt(norm(X1-X2)+Θ))^3)
end

#Linear Kernel
function linear(X1::Array{Float64,1},X2::Array{Float64,1})
  dot(X1,X2)
end

#Quadratic Kernel (special case of polynomial kernel)
function quadratic(X1::Array{Float64,1},X2::Array{Float64,1},Θ)
  (dot(X1,X2)+Θ)^2
end

function deriv_quadratic(X1::Array{Float64,1},X2::Array{Float64,1},Θ)
  2*(dot(X1,X2)+Θ)
end

end #end of module


text:
### The implementation of AM in Klara uses the mixture proposal of the following article:
### Gareth O. Roberts and Jeffrey S. Rosenthal
### Examples of Adaptive MCMC
### Journal of Computational and Graphical Statistics, 2012, 18 (2), pp 349-367

### The original adaptive Metropolis algorithm was introduced by the following article:
### Heikki Haario, Eero Saksman and Johanna Tamminen
### An Adaptive Metropolis Algorithm
### Bernoulli, 2001, 7 (2), pp 223-242

### Abstract AM state

abstract type AMState{F<:VariateForm} <: MHSamplerState{F} end

### AM state subtypes

## UnvAMState holds the internal state ("local variables") of the AM sampler for univariate parameters

mutable struct UnvAMState <: AMState{Univariate}
  proposal::Union{UnivariateGMM, RealNormal}
  pstate::ParameterState{Continuous, Univariate} # Parameter state used internally by AM
  tune::MCTunerState
  ratio::Real # Acceptance ratio
  lastmean::Real
  secondlastmean::Real
  C::Real
  sqrtminorscale::Real
  w::RealVector
  count::Integer
  diagnosticindices::Dict{Symbol, Integer}

  function UnvAMState(
    proposal::Union{UnivariateGMM, RealNormal},
    pstate::ParameterState{Continuous, Univariate},
    tune::MCTunerState,
    ratio::Real,
    lastmean::Real,
    secondlastmean::Real,
    C::Real,
    sqrtminorscale::Real,
    w::RealVector,
    count::Integer,
    diagnosticindices::Dict{Symbol, Integer}
  )
    if !isnan(ratio)
      @assert ratio > 0 "Acceptance ratio should be positive"
    end
    if !isnan(sqrtminorscale)
      @assert sqrtminorscale >= 0 "Scaling of stabilizing covariance should be non-negative"
    end
    if !isnan(w[1])
      @assert w[1] > 0 "Weight of core mixture component must be positive"
    end
    if !isnan(w[2])
      @assert w[2] >= 0 "Weight of minor mixture component must be non-negative"
    end
    @assert count >= 0 "Number of iterations (count) should be non-negative"
    new(proposal, pstate, tune, ratio, lastmean, secondlastmean, C, sqrtminorscale, w, count, diagnosticindices)
  end
end

UnvAMState(
  proposal::Union{UnivariateGMM, RealNormal},
  pstate::ParameterState{Continuous, Univariate},
  C::Real,
  sqrtminorscale::Real,
  w::RealVector,
  tune::MCTunerState=BasicMCTune(),
  lastmean::Real=NaN
) =
  UnvAMState(proposal, pstate, tune, NaN, lastmean, lastmean, C, sqrtminorscale, w, 0, Dict{Symbol, Integer}())

## MuvAMState holds the internal state ("local variables") of the AM sampler for multivariate parameters

mutable struct MuvAMState <: AMState{Multivariate}
  proposal::Union{MultivariateGMM, AbstractMvNormal}
  pstate::ParameterState{Continuous, Multivariate} # Parameter state used internally by AM
  tune::MCTunerState
  ratio::Real # Acceptance ratio
  lastmean::RealVector
  secondlastmean::RealVector
  C::RealMatrix
  sqrtminorscale::Real
  w::RealVector
  count::Integer
  diagnosticindices::Dict{Symbol, Integer}

  function MuvAMState(
    proposal::Union{MultivariateGMM, AbstractMvNormal},
    pstate::ParameterState{Continuous, Multivariate},
    tune::MCTunerState,
    ratio::Real,
    lastmean::RealVector,
    secondlastmean::RealVector,
    C::RealMatrix,
    sqrtminorscale::Real,
    w::RealVector,
    count::Integer,
    diagnosticindices::Dict{Symbol, Integer}
  )
    if !isnan(ratio)
      @assert ratio > 0 "Acceptance ratio should be positive"
    end
    if !isnan(sqrtminorscale)
      @assert sqrtminorscale >= 0 "Scaling of stabilizing covariance should be non-negative"
    end
    if !isnan(w[1])
      @assert w[1] > 0 "Weight of core mixture component must be positive"
    end
    if !isnan(w[2])
      @assert w[2] >= 0 "Weight of minor mixture component must be non-negative"
    end
    @assert count >= 0 "Number of iterations (count) should be non-negative"
    new(proposal, pstate, tune, ratio, lastmean, secondlastmean, C, sqrtminorscale, w, count, diagnosticindices)
  end
end

MuvAMState(
  proposal::Union{MultivariateGMM, AbstractMvNormal},
  pstate::ParameterState{Continuous, Multivariate},
  C::RealMatrix,
  sqrtminorscale::Real,
  w::RealVector,
  tune::MCTunerState=BasicMCTune(),
  lastmean::RealVector=Array{eltype(pstate)}(pstate.size)
) =
  MuvAMState(proposal, pstate, tune, NaN, lastmean, lastmean, C, sqrtminorscale, w, 0, Dict{Symbol, Integer}())

### Adaptive Metropolis (AM) sampler

struct AM <: MHSampler
  C0::RealMatrix # Initial covariance of proposal according to best prior knowledge
  corescale::Real # Scaling factor of covariance matrix of the core mixture component
  minorscale::Real # Scaling factor of covariance matrix of the stabilizing mixture component
  c::Real # Non-negative constant with relative small value that determines the mixture weight of the stabilizing component
  t0::Integer

  function AM(C0::RealMatrix, corescale::Real, minorscale::Real, c::Real, t0::Integer)
    @assert 0 < corescale "Constant corescale must be positive, got $corescale"
    @assert 0 < minorscale "Constant minorscale must be positive, got $minorscale"
    @assert 0 <= c "Constant c must be non-negative"
    @assert t0 > 0 "t0 is not positive"
    new(C0, corescale, minorscale, c, t0)
  end
end

AM(C0::RealMatrix; corescale::Real=1., minorscale::Real=1., c::Real=0.05, t0::Integer=10) =
  AM(C0, corescale, minorscale, c, t0)

AM(C0::RealVector; corescale::Real=1., minorscale::Real=1., c::Real=0.05, t0::Integer=10) =
  AM(diagm(C0), corescale=corescale, minorscale=minorscale, c=c, t0=t0)

AM(C0::Real, d::Integer=1; corescale::Real=1., minorscale::Real=1., c::Real=0.05, t0::Integer=10) =
  AM(fill(C0, d), corescale=corescale, minorscale=minorscale, c=c, t0=t0)

### Initialize AM sampler

## Initialize parameter state

function initialize!(
  pstate::ParameterState{Continuous, F},
  parameter::Parameter{Continuous, F},
  sampler::AM,
  outopts::Dict
) where F<:VariateForm
  parameter.logtarget!(pstate)
  @assert isfinite(pstate.logtarget) "Log-target not finite: initial value out of support"

  if !isempty(outopts[:diagnostics])
    pstate.diagnostickeys = copy(outopts[:diagnostics])
    pstate.diagnosticvalues = Array{Any}(length(pstate.diagnostickeys))
  end
end

## Initialize AM state

set_gmm(sampler::AM, pstate::ParameterState{Continuous, Univariate}, C::Real, sqrtminorscale::Real, w::RealVector) =
  UnivariateGMM([pstate.value, pstate.value], Float64[sqrt(sampler.corescale*C), sqrtminorscale], Categorical(w))

set_gmm!(sstate::UnvAMState, sampler::AM, pstate::ParameterState{Continuous, Univariate}) =
  sstate.proposal = set_gmm(sampler, pstate, sstate.C, sstate.sqrtminorscale, sstate.w)

set_normal(sampler::AM, pstate::ParameterState{Continuous, Univariate}, sqrtminorscale::Real) =
  Normal(pstate.value, sqrtminorscale)

set_normal!(sstate::UnvAMState, sampler::AM, pstate::ParameterState{Continuous, Univariate}) =
  sstate.proposal = set_normal(sampler, pstate, sstate.sqrtminorscale)

function sampler_state(
  parameter::Parameter{Continuous, Univariate},
  sampler::AM,
  tuner::MCTuner,
  pstate::ParameterState{Continuous, Univariate},
  vstate::VariableStateVector,
  diagnostickeys::Vector{Symbol}
)
  sqrtminorscale = sqrt(sampler.minorscale)
  w = [1-sampler.c, sampler.c]

  sstate = UnvAMState(
    set_normal(sampler, pstate, sqrtminorscale),
    generate_empty(pstate),
    sampler.C0[1, 1],
    sqrtminorscale,
    w,
    tuner_state(parameter, sampler, tuner),
    pstate.value
  )

  set_diagnosticindices!(sstate, [:accept], diagnostickeys)

  sstate
end

set_gmm(sampler::AM, pstate::ParameterState{Continuous, Multivariate}, C::RealMatrix, sqrtminorscale::Real, w::RealVector) =
  MixtureModel([MvNormal(pstate.value, sampler.corescale*C), MvNormal(pstate.value, sqrtminorscale)], w)

set_gmm!(sstate::MuvAMState, sampler::AM, pstate::ParameterState{Continuous, Multivariate}) =
  sstate.proposal = set_gmm(sampler, pstate, sstate.C, sstate.sqrtminorscale, sstate.w)

set_normal(sampler::AM, pstate::ParameterState{Continuous, Multivariate}, sqrtminorscale::Real) =
  MvNormal(pstate.value, sqrtminorscale)

set_normal!(sstate::MuvAMState, sampler::AM, pstate::ParameterState{Continuous, Multivariate}) =
  sstate.proposal = set_normal(sampler, pstate, sstate.sqrtminorscale)

function sampler_state(
  parameter::Parameter{Continuous, Multivariate},
  sampler::AM,
  tuner::MCTuner,
  pstate::ParameterState{Continuous, Multivariate},
  vstate::VariableStateVector,
  diagnostickeys::Vector{Symbol}
)
  sqrtminorscale = sqrt(sampler.minorscale)
  w = [1-sampler.c, sampler.c]

  sstate = MuvAMState(
    set_normal(sampler, pstate, sqrtminorscale),
    generate_empty(pstate),
    copy(sampler.C0),
    sqrtminorscale,
    w,
    tuner_state(parameter, sampler, tuner),
    copy(pstate.value)
  )

  set_diagnosticindices!(sstate, [:accept], diagnostickeys)

  sstate
end

## Reset parameter state

function reset!(
  pstate::ParameterState{Continuous, Univariate},
  x::Real,
  parameter::Parameter{Continuous, Univariate},
  sampler::AM
)
  pstate.value = x
  parameter.logtarget!(pstate)
end

function reset!(
  pstate::ParameterState{Continuous, Multivariate},
  x::RealVector,
  parameter::Parameter{Continuous, Multivariate},
  sampler::AM
)
  pstate.value = copy(x)
  parameter.logtarget!(pstate)
end

## Reset sampler state

function reset!(
  sstate::UnvAMState,
  pstate::ParameterState{Continuous, Univariate},
  parameter::Parameter{Continuous, Univariate},
  sampler::AM,
  tuner::MCTuner
)
  set_normal!(sstate, sampler, pstate)
  reset!(sstate.tune, sampler, tuner)
  sstate.lastmean = pstate.value
  sstate.C0 = sampler.C0[1, 1]
  sstate.C = sstate.C0
  sstate.count = 0
end

function reset!(
  sstate::MuvAMState,
  pstate::ParameterState{Continuous, Multivariate},
  parameter::Parameter{Continuous, Multivariate},
  sampler::AM,
  tuner::MCTuner
)
  set_normal!(sstate, sampler, pstate)
  reset!(sstate.tune, sampler, tuner)
  sstate.lastmean = copy(pstate.value)
  sstate.C = copy(sampler.C0)
  sstate.count = 0
end

show(io::IO, sampler::AM) = print(
  io,
  "AM sampler: scaling of core covariance = ",
  sampler.corescale,
  ", scaling of stabilizing covariance = ",
  sampler.minorscale,
  ", weight of stabilizing mixture component = ",
  sampler.c
)


text:
push!(LOAD_PATH, ".")

module GACrossover

using Debug
using Individual
using RealIndividual
using Parallel

abstract type PointCrossover end
abstract type AritmeticCrossover end

function crossover(p1::IndType, p2::IndType, ::Type{PointCrossover}, cr::Float64) where {IndType <: Individual.AbstractIndividual}
  gene_num = Individual.getNumGenes(p1)
  c1, c2 = deepcopy(p1), deepcopy(p2)

  Debug.ga_debug && println("----- Point Crossover -----\n")

  pr = Parallel.threadRand()

  Debug.ga_debug && println("Test: ", pr, " < ", cr, "\n")

  if pr < cr

    if Debug.ga_debug
      println("Parent 1: ", c1[:])
      println("Parent 2: ", c2[:], "\n")
    end

    point = Parallel.threadRand(1:(gene_num-1))

    Debug.ga_debug && println("Cut point: ", point, "\n")

    cuts = [1:point, (point + 1):gene_num]
    chosen_cut = Parallel.threadRand(cuts)
    c1[chosen_cut], c2[chosen_cut] = c2[chosen_cut], c1[chosen_cut]

    if Debug.ga_debug
      println("Child 1: ", c1[:])
      println("Child 2: ", c2[:], "\n")
    end

  end

  Debug.ga_debug && println("----- Point Crossover End -----\n")

  return (c1, c2)
end

function crossover(p1::IndType, p2::IndType, ::Type{AritmeticCrossover}, cr::Float64) where {IndType <: RealIndividual.AbstractRealIndividual}
  c1, c2 = deepcopy(p1), deepcopy(p2)

  pr = Parallel.threadRand()
  if pr < cr
    α = Parallel.threadRand()
    c1[:] = α * c1 + (1 - α) * c2
    c2[:] = α * c2 + (1 - α) * c1
  end
  return (c1, c2)
end

end

max_stars_repo_name: Markhyz/Julia-Genetic-Algorithm

max_issues_repo_name: Markhyz/Julia-Genetic-Algorithm

max_forks_repo_name: Markhyz/Julia-Genetic-Algorithm


text:
using Test

include("_includes/julia/preamble_using.jl")

macro test_example(name, test_code)
    file = "_includes/julia/$name.jl"
    quote
        @testset "$($name)" begin
            include("_includes/julia/preamble.jl")
            include($file)
            $test_code
        end
    end
end

@test_example "example1" begin
    @test termination_status(model) == MOI.OPTIMAL
end
@test_example "example1_motzkin" begin
    @test termination_status(model) == MOI.INFEASIBLE
end
@test_example "example2" begin
    @test termination_status(model) == MOI.OPTIMAL
    @test value(s) ≈  4.2956  rtol=1e-4
    @test value(t) ≈  4.2956  rtol=1e-4
end
@test_example "example2_objective" begin
    @test value(s) ≈  1.2032  rtol=1e-4
    @test value(t) ≈  1.2032  rtol=1e-4
end
@test_example "example2_both" begin
    @test value(t) ≈  1/4     rtol=1e-6
end
@test_example "example3" begin
    @test value(t) ≈ -0.1778  rtol=1e-3
end
@test_example "example4" begin
    @test value(t) ≈ -1       rtol=1e-6
end
@test_example "example4_inequalities" begin
    @test value(t) ≈  1.3911  rtol=1e-4
end

max_stars_repo_name: sums-of-squares/tools

max_issues_repo_name: sums-of-squares/tools

max_forks_repo_name: sums-of-squares/tools


text:
using BioMedQuery.Processes

const email= ""
#Enviroment variable that need to be setup
const umls_user = get(ENV, "UMLS_USER", "")
const umls_pswd = get(ENV, "UMLS_PSSWD", "")
const search_term="(obesity[MeSH Major Topic]) AND (\"2010\"[Date - Publication])"
const max_articles = 2
const verbose = false

# global dfs = DataFrame()

@testset "Search and Parse" begin
    println("-----------------------------------------")
    println("       Testing Search and Parse")

    global dfs = pubmed_search_and_parse(email, search_term, max_articles, verbose)

    @test size(dfs["basic"])[1] == max_articles

end

credentials_set = get(ENV, "TRAVIS_SECURE_ENV_VARS", "true")=="true" && umls_user != ""

@testset "UMLS DataFrames" begin
    if credentials_set
        println("-----------------------------------------")
        println("       Testing MESH2UMLS")
        m2u = map_mesh_to_umls_async(dfs["mesh_desc"], umls_user, umls_pswd)

        @test size(m2u)[1] > length(dfs["mesh_desc"])

        println("-----------------------------------------")
        println("       Testing Occurences")

        umls_concept = "Disease or Syndrome"

        labels2ind, occur = umls_semantic_occurrences(dfs, m2u, umls_concept)

        @test length(keys(labels2ind)) > 0
        @test length(findall(x->x=="Obesity", collect(keys(labels2ind)))) ==1
    else
        @warn "Skipping UMLS tests as no credentials provided"
    end
end


text:
include("../server/simulation/simulation.jl")

using .Simulation

function test_tiny_delay()
    println("Testing tiny delay")

    stream = Simulation.PoissonStream(UInt64(13163), 0.2)

    delay = Simulation.TinyDelay(stream)

    @assert Simulation.output(delay) == 0 "Delay output not zero at start"

    expected_output = [0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,1,0,0,0,0]
    for i in 1:20
        Simulation.pre!(delay)
        print(Simulation.output(delay), ",")
        # @assert Simulation.output(delay) == expected_output[i] "Expected output did not match: $i"
        Simulation.post!(delay)
    end
    println("")

    println("Passed")
end

function test_tiny_delay3()
    println("Testing tiny delay of 3")

    stream = Simulation.PoissonStream(UInt64(13163), 0.2)

    delay = Simulation.TinyDelay(stream, 3)

    @assert Simulation.output(delay) == 0 "Delay output not zero at start"

    # Not shifted     [0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,1,0,0,0,0]
    expected_output = [0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,1,0] # delayed
    for i in 1:20
        Simulation.pre!(delay)
        # print(Simulation.output(delay), ",")
        @assert Simulation.output(delay) == expected_output[i] "Expected output did not match: $i"
        Simulation.post!(delay)
    end
    # println("")

    println("Passed")
end

function test_small_delay()
    println("Testing small delay")

    stream = Simulation.PoissonStream(UInt64(13163), 0.2)

    delay = Simulation.SmallDelay(stream, 10)

    @assert Simulation.output(delay) == 0 "Delay output not zero at start"

    expected_output = [0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,1,0,0,0,1,0,1,0,0,0,1,0]
    for i in 1:100
        Simulation.pre!(delay)
        print(Simulation.output(delay), ",")
        # @assert Simulation.output(delay) == expected_output[i] "Expected output did not match: $i"
        Simulation.post!(delay)
    end
    println("")
    println("Passed")
end

function test_medium_delay()
    println("Testing medium delay")

    stream = Simulation.PoissonStream(UInt64(13163), 0.2)

    delay = Simulation.MediumDelay(stream, 20)

    @assert Simulation.output(delay) == 0 "Delay output not zero at start"

    expected_output = [0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,1,0,0,0,1,0,1,0,0,0,1,0]
    for i in 1:100
        Simulation.pre!(delay)
        print(Simulation.output(delay), ",")
        # @assert Simulation.output(delay) == expected_output[i] "Expected output did not match: $i"
        Simulation.post!(delay)
    end
    println("")
    println("Passed")
end

test_medium_delay()


text:
abstract type ModelAtmosphere end

struct PlanarAtmosphereLayer{F}
    colmass::F                 #g
    temp::F                    #K
    electron_number_density::F #cm^-3
    number_density::F          #cm^-3
    density::F                 #g cm^-3
end

struct PlanarAtmosphere{F} <: ModelAtmosphere
    layers::Vector{PlanarAtmosphereLayer{F}}
end

struct ShellAtmosphereLayer{F}
    colmass::F                 #g
    temp::F                    #K
    electron_number_density::F #cm^-3
    number_density::F          #cm^-3
    density::F                 #g cm^-3
    r::F                       #cm
end

struct ShellAtmosphere{F} <: ModelAtmosphere
    layers::Vector{ShellAtmosphereLayer{F}}
end

"""
    PlanarAtmosphere(atm::ShellAtmosphere)

Construct a planar atmosphere with the data from a shell atmosphere.  Mostly useful for testing.
"""
function PlanarAtmosphere(atm::ShellAtmosphere)
    PlanarAtmosphere([PlanarAtmosphereLayer(l.colmass, l.temp, l.electron_number_density,
                                            l.number_density, l.density) for l in atm.layers])
end

"""
    ShellAtmosphere(atm::PlanarAtmosphere, R)

Construct a shell atmosphere with the data from a planar atmosphere and an outer radius.  Mostly
useful for testing.
"""
function ShellAtmosphere(atm::PlanarAtmosphere, R)
    Δcolmass = diff((l->l.colmass).(atm.layers))
    Δs = 0.5([0 ; Δcolmass] + [Δcolmass; Δcolmass[end]]) ./ (l->l.density).(atm.layers)
    rs = R .- cumsum(Δs)

    ShellAtmosphere([ShellAtmosphereLayer(l.colmass, l.temp, l.electron_number_density,
                                          l.number_density, l.density, r)
                     for (l, r) in zip(atm.layers, rs)])
end

#pretty-printing
function Base.show(io::IO, m::MIME"text/plain", atm::A) where A <: ModelAtmosphere
    print(io, "$(A) with $(length(atm.layers)) layers")
end

"""
    read_model_atmosphere(filename; truncate_at_10000K=true)

Parse the provided model atmosphere file in
["Kurucz" format](https://marcs.astro.uu.se/krz_format.html). Returns either a `PlanarAtmosphere` or
a `ShellAtmosphere`.

When `truncate_at_10000K` is true, layers with temperatures greater than 10,000 Kelvin will be
elided.  These are typically at the deepest level, where optical depth is very large, and have
minimal impact on the surface spectrum.  They are ignored by default because Korg's default
partition and molecular equillibrium functions are tabulated only up to that tempurature.
"""
function read_model_atmosphere(fname::AbstractString; truncate_at_10000K=true)
    open(fname) do f
        #these files are small, so it's not a big deal to load them entirely into memory
        lines = collect(eachline(f))

        toks = strip.(split(lines[2], ' '))
        spherical = "RADIUS=" in toks
        R_header = spherical ? parse(Float64, toks[end-1]) : NaN

        layers = map(lines[14:end]) do line
            toks = strip.(split(line, ','))
            if toks[end] == ""
                toks = toks[1:end-1]
            end
            numbers = parse.(Float64, toks)
            if spherical
                #spherical shell atmospheres have a sixth column with distance from reference height
                numbers[6] += R_header
                ShellAtmosphereLayer(numbers...)
            else
                PlanarAtmosphereLayer(numbers...)
            end
         end

        if truncate_at_10000K
            filter!(layers) do layer
                layer.temp < 10_000
            end
        end

         if spherical
             ShellAtmosphere(layers)
         else
             PlanarAtmosphere(layers)
         end
    end
end


text:
module CStructures

include("layout.jl")
include("cstruct.jl")

end


text:
# Internal prediction function. Used at the end of shap().
using Statistics

function _predict(;reference::DataFrame, data_predict::DataFrame, model, predict_function::Function, n_features::Integer,
                  n_target_features::Integer, n_instances_explain::Integer, sample_size::Integer, precision::Union{Integer, Nothing},
                  chunk::Bool, reconcile_instance::Bool, explain::DataFrame)

  if !chunk
    # For robustness across prediction functions from other languages (R and Python),
    # the modeling dataset needs to be created before it's passed in predict_function().
    # This is data-to-be-predicted but the name is kept the same throughout to reduce memory.
    data_predicted = data_predict[:, 1:n_features]

    # User-defined predict() function.
    data_predicted = predict_function(model, data_predicted)

    # A dataset with meta-data and the predictions in the last column.
    data_predicted = hcat(data_predict[:, (n_features + 1):size(data_predict, 2)], data_predicted, copycols = false)

  else

    data_predicted = copy(data_predict)
  end
  #----------------------------------------------------------------------------
  # Cast the DataFrame to, for each random sample, take the difference between the Frankenstein
  # instances.
  user_fun_y_pred_name = names(data_predicted)[end]

  # Returns a length 1 numeric vector of the average prediction--i.e., intercept--from the reference group.
  intercept = Statistics.mean(data_predicted[:, user_fun_y_pred_name])

  data_predicted = DataFrames.unstack(data_predicted, [:index, :sample, :feature_name], :feature_group, user_fun_y_pred_name)

  # Shapley value for each Monte Carlo sample for each instance.
  data_predicted.shap_effect = data_predicted.real_target - data_predicted.fake_target
  #----------------------------------------------------------------------------
  DataFrames.select!(data_predicted, [:index, :sample, :feature_name, :shap_effect])
  #----------------------------------------------------------------------------
  # Final Shapley value calculation collapsed across Monte Carlo samples using a
  # custom _aggregate() function for speed and reduced memory. The dimensions of
  # the output DataFrame, 'data_predicted', are ('n instances' * 'n features') by 5.
  if reconcile_instance
    data_predicted = _aggregate(data_predicted, sample_size, n_instances_explain, n_target_features, reconcile_instance)
  else
    data_predicted = _aggregate(data_predicted, sample_size, n_instances_explain, n_target_features, reconcile_instance)
  end
  #----------------------------------------------------------------------------
  # Adjust the instance-level Shapley values so that the sum across features equals the model prediction.
  if reconcile_instance
    # User-defined predict() function.
    data_model_pred = predict_function(model, explain)

    rename!(data_model_pred, Dict(user_fun_y_pred_name => "model_pred"))

    data_model_pred[:, :index] .= 1:n_instances_explain

    # data_predicted = join(data_predicted, data_model_pred, on = [:index], kind = :left)
    data_predicted = leftjoin(data_predicted, data_model_pred, on = [:index])
    data_shap_pred = DataFrames.combine(DataFrames.groupby(data_predicted, [:index]), :shap_effect => (x -> sum(x)) => :shap_pred)
    data_shap_pred.shap_pred .= data_shap_pred.shap_pred .+ intercept

    # data_predicted = join(data_predicted, data_shap_pred, on = [:index], kind = :left)
    data_predicted = leftjoin(data_predicted, data_shap_pred, on = [:index])

    data_predicted[:, :residual] = data_predicted.model_pred - data_predicted.shap_pred

    data_scale = DataFrames.combine(DataFrames.groupby(data_predicted, [:index]), :shap_effect_var => (x -> x / maximum(x)) => :variance_scale)

    data_predicted[:, :variance_scale] .= data_scale.variance_scale

    data_scale = DataFrames.combine(DataFrames.groupby(data_predicted, [:index]), :variance_scale => (x -> sum(x)) => :variance_scale_total)

    # data_predicted = join(data_predicted, data_scale, on = [:index], kind = :left)
    data_predicted = leftjoin(data_predicted, data_scale, on = [:index])

    data_predicted[:, :shap_effect] = data_predicted.shap_effect .+ (data_predicted.residual .* (data_predicted.variance_scale .- (data_predicted.variance_scale .* data_predicted.variance_scale_total) ./
                                                                 (1 .+ data_predicted.variance_scale_total)))

    DataFrames.select!(data_predicted, [:index, :feature_name, :shap_effect, :shap_effect_sd])
  end
  #----------------------------------------------------------------------------
  data_predicted.intercept = repeat([intercept], size(data_predicted, 1))

  # Optional rounding of Shapley results to reduce dataset size.
  if precision !== nothing
    data_predicted.intercept .= round.(data_predicted.intercept, digits = precision)
    data_predicted.shap_effect .= round.(data_predicted.shap_effect, digits = precision)
    data_predicted.shap_effect_sd .= round.(data_predicted.shap_effect_sd, digits = precision)
  end

  return data_predicted
end  # End _predict().

max_stars_repo_name: nredell/shapML

max_issues_repo_name: nredell/shapML

max_forks_repo_name: nredell/shapML


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

function inflate_ir(ci::CodeInfo, linfo::MethodInstance)
    sptypes = sptypes_from_meth_instance(linfo)
    if ci.inferred
        argtypes, _ = matching_cache_argtypes(linfo, nothing, false)
    else
        argtypes = Any[ Any for i = 1:length(ci.slotflags) ]
    end
    return inflate_ir(ci, sptypes, argtypes)
end

function inflate_ir(ci::CodeInfo, sptypes::Vector{Any}, argtypes::Vector{Any})
    code = copy_exprargs(ci.code) # TODO: this is a huge hot-spot
    cfg = compute_basic_blocks(code)
    for i = 1:length(code)
        stmt = code[i]
        # Translate statement edges to bb_edges
        if isa(stmt, GotoNode)
            code[i] = GotoNode(block_for_inst(cfg, stmt.label))
        elseif isa(stmt, GotoIfNot)
            code[i] = GotoIfNot(stmt.cond, block_for_inst(cfg, stmt.dest))
        elseif isa(stmt, PhiNode)
            code[i] = PhiNode(Int32[block_for_inst(cfg, Int(edge)) for edge in stmt.edges], stmt.values)
        elseif isa(stmt, Expr) && stmt.head === :enter
            stmt.args[1] = block_for_inst(cfg, stmt.args[1])
            code[i] = stmt
        else
            code[i] = stmt
        end
    end
    nstmts = length(code)
    ssavaluetypes = let ssavaluetypes = ci.ssavaluetypes
        ssavaluetypes isa Vector{Any} ? copy(ssavaluetypes) : Any[ Any for i = 1:(ssavaluetypes::Int) ]
    end
    stmts = InstructionStream(code, ssavaluetypes, Any[nothing for i = 1:nstmts], copy(ci.codelocs), copy(ci.ssaflags))
    ir = IRCode(stmts, cfg, collect(LineInfoNode, ci.linetable), argtypes, Any[], sptypes)
    return ir
end

function replace_code_newstyle!(ci::CodeInfo, ir::IRCode, nargs::Int)
    @assert isempty(ir.new_nodes)
    # All but the first `nargs` slots will now be unused
    resize!(ci.slotflags, nargs)
    stmts = ir.stmts
    ci.code, ci.ssavaluetypes, ci.codelocs, ci.ssaflags, ci.linetable =
        stmts.inst, stmts.type, stmts.line, stmts.flag, ir.linetable
    for metanode in ir.meta
        push!(ci.code, metanode)
        push!(ci.codelocs, 1)
        push!(ci.ssavaluetypes::Vector{Any}, Any)
        push!(ci.ssaflags, IR_FLAG_NULL)
    end
    # Translate BB Edges to statement edges
    # (and undo normalization for now)
    for i = 1:length(ci.code)
        stmt = ci.code[i]
        if isa(stmt, GotoNode)
            stmt = GotoNode(first(ir.cfg.blocks[stmt.label].stmts))
        elseif isa(stmt, GotoIfNot)
            stmt = GotoIfNot(stmt.cond, first(ir.cfg.blocks[stmt.dest].stmts))
        elseif isa(stmt, PhiNode)
            stmt = PhiNode(Int32[last(ir.cfg.blocks[edge].stmts) for edge in stmt.edges], stmt.values)
        elseif isa(stmt, Expr) && stmt.head === :enter
            stmt.args[1] = first(ir.cfg.blocks[stmt.args[1]::Int].stmts)
        end
        ci.code[i] = stmt
    end
end

# used by some tests
inflate_ir(ci::CodeInfo) = inflate_ir(ci, Any[], Any[ Any for i = 1:length(ci.slotflags) ])

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
using StaticArrays
const Vec = SVector

function writefbodyline(f, indent, s)
    for i in 1:indent
        write(f, "\t")
    end
    return write(f, "$s\n")
end

function fort2julia(f_type)
    f_type = lowercase(f_type)
    if f_type == "real"
        return Float32
    elseif f_type == "real(kind=dp)"
        return Float64
    elseif f_type == "complex(kind=dp)"
        return Complex{Float64}
    elseif occursin("character", f_type)
        return String
    elseif f_type == "string"
        return String
    elseif f_type == "integer"
        return Int
    elseif f_type == "logical"
        return Bool
    elseif occursin(".D", f_type)
        return replace(f_type, "D" => "e")
    else
        return Nothing
    end
end
indentabs(indent) = prod(["\t" for i in 1:indent])

open(joinpath(@__DIR__, "mpirunflags.jl"), "w") do wf
    write(wf, "_MPI_FLAGS() = ExecFlag[\n")
    # writefbodyline(1, "flags = ")
    open(joinpath(@__DIR__, "..", "assets", "mpirun_man.txt"), "r") do f
        line = readline(f)
        while line != "OPTIONS"
            line = readline(f)
        end
        while line != "Environment Variables"
            line = strip(readline(f))
            if !isempty(line) && line[1] == '-'
                name        = ""     #--
                symbols     = Symbol[] #-
                type        = Nothing
                description = ""
                sline       = strip.(split(line), ',')
                for s in sline
                    if s[2] == '-' #--npernode
                        name = strip(s, '-')
                    elseif occursin('<', s) # <#persocket>
                        type = if occursin("#", s)
                            occursin(',', s) ? Vector{Int} : Int
                        elseif occursin("ppr", s) #ppr:N:<object>
                            Pair{Int,String}
                        else
                            occursin(',', s) ? Vector{String} : String
                        end
                        break
                    else  #-np
                        push!(symbols, Symbol(strip(s, '-')))
                    end
                end
                line = strip(readline(f))
                while !isempty(line)
                    description *= join(split(line), " ") * "\n" * indentabs(2)
                    line = strip(readline(f))
                end
                description = replace(strip(description), "\"" => "'")
                if name != "" && isempty(symbols)
                    symbols = [Symbol(name)]
                end
                if type <: Number
                    val = zero(type)
                elseif type <: Pair
                    val = 0 => ""
                elseif type === Nothing
                    val = nothing
                else
                    val = "\"\""
                end
                for symbol in symbols
                    writefbodyline(wf, 1,
                                   """ExecFlag(Symbol("$symbol"), "$name", "$description", $val, 1),""")
                end
            end
        end
    end
    return writefbodyline(wf, 0, "]")
end

open(joinpath(@__DIR__, "wannier90flags.jl"), "w") do wf
    write(wf, "_WAN_FLAGS() = Dict{Symbol, Type}(\n")
    open(joinpath(@__DIR__, "..", "assets", "calculations", "wannier",
                  "calculation_flags.txt"), "r") do f
        while !eof(f)
            line = readline(f)
            if isempty(line) || line[1] == '!'
                continue
            else
                s_line = split(line, "::")
                flagpart = s_line[end]
                flag = Symbol(strip(strip(split(split(split(s_line[end], "=")[1], "(")[1],
                                                "!")[1], ':')))
                fl_type = fort2julia(strip(split(s_line[1])[1], ','))
                flagstring = occursin("(", flagpart) && occursin(")", flagpart) ?
                             """Symbol("$flag") => Vector{$fl_type},""" :
                             """Symbol("$flag") => $fl_type,"""
                writefbodyline(wf, 1, flagstring)
            end
        end
    end
    writefbodyline(wf, 1, ":preprocess => Bool") #temporary hack
    return write(wf, ")")
end

strip_split(line, args...) = strip.(split(line, args...))
function read_block(f, startstr::String, endstr::String)
    block = [startstr]
    line = readline(f)
    while !eof(f)
        while !occursin(endstr, line)
            line = readline(f)
            push!(block, line)
        end
        return block
    end
    return error("Block not found: start = $startstr, end = $endstr.")
end

function vardim(line)
    if all(occursin.(['(', ')'], (line,)))
        dim = length(split(split(split(line, '(')[2], ')')[1], ','))
    else
        dim = 0
    end
    return dim
end

function qe_write_variable(wf, indent, lines, i)
    name = gensym()
    var_i = i
    i += 2
    typ = fort2julia(strip_split(lines[i])[2])
    description = ""
    # default = nothing
    i += 1
    line = lines[i]
    while !occursin("+------", line)
        # if occursin("Default", line)
        #     _t = strip_split(line)[2]
        #     _t = strip(strip(_t,'('),')')
        #     if occursin("D", _t)
        #         default = Meta.parse(typ, replace(_t,"D" => "e"))
        #     else
        #         _t = occursin("=",_t) ?split(_t,"=")[end] : _t
        #         default = typ ==String ? _t : Meta.parse(_t)
        #         println(_t)
        #         if typeof(default) != Symbol
        #             default = convert(typ, default)
        #         end
        #     end
        if occursin("Description", line)
            description *= strip_split(line, ":")[2] * "\n"
            i += 1
            line = lines[i]
            while !occursin("+------", line)
                description *= strip(lines[i]) * "\n"
                i += 1
                line = lines[i]
            end
            @goto break_label
        end
        i += 1
        line = lines[i]
    end
    @label break_label
    line = lines[var_i]
    dim = vardim(line)
    typ = dim == 0 ? typ : Array{typ,dim}
    description = replace(replace(replace(replace(replace(description, "\"" => "'"),
                                                  "\\" => "\\\\"), "\$" => "\\\$"),
                                  "\\\\t" => "\\t"), "\\\\n" => "\\n")
    if occursin("Variables", line)
        spl = [split(x, "(")[1]
               for x in strip.(filter(x -> !occursin("=", x), split(line)[2:end]), ',')]
        names = Symbol.(spl)
        for name in names
            writefbodyline(wf, indent,
                           """QEFlagInfo{$typ}(Symbol("$name"), "$description"),""")
        end
        return i
    else
        if occursin("(", line) && occursin(")", line)
            name = Symbol(split(strip_split(line, ":")[end], "(")[1])
        else
            name = Symbol(strip_split(line, ":")[end])
        end
        writefbodyline(wf, indent, """QEFlagInfo{$typ}(Symbol("$name"), "$description"),""")
        return i
    end
end

function write_QEControlBlockInfo(wf, indent, lines)
    name = Symbol(lowercase(strip_split(lines[1], "&")[2]))
    writefbodyline(wf, indent, """QEControlBlockInfo(Symbol("$name"), [""")
    for i in 1:length(lines)
        line = lines[i]
        if occursin("Variable", line)
            i += qe_write_variable(wf, indent + 1, lines, i)
        end
    end
    return writefbodyline(wf, indent + 1, "]),")
    # writefbodyline(wf, indent, "),")
end

function write_QEDataBlockInfo(wf, indent, lines)
    spl                 = split(strip(lines[1]))
    name                = length(spl) > 1 ? lowercase(spl[2]) : "noname"
    options             = Symbol.(spl[4:2:end])
    description         = ""
    options_description = ""
    istop               = findfirst(x -> occursin("DESCRIPTION", x), lines)
    if istop === nothing
        return
    end
    for i in 1:istop-1
        line = strip(replace(replace(lines[i], "_" => ""), "-" => ""))
        if isempty(line)
            continue
        end
        description *= strip(line) * "\n" * indentabs(indent + 2)
    end
    description = replace(description, "\"" => "'")

    istart = findnext(x -> occursin("Card's flags", x), lines, istop)
    if istart != nothing
        istop = findnext(x -> occursin("+------", x), lines, istart)
        for i in istart+1:istop-1
            line = lines[i]
            if occursin("Default", line) || occursin("Description", line)
                continue
            end
            options_description *= strip(line) * "\n" * indentabs(indent + 2)
        end
    end

    options_description = replace(options_description, "\"" => "'")
    writefbodyline(wf, indent + 1,
                   """QEDataBlockInfo(Symbol("$name"),"$description", $options, "$options_description", [""")
    i = istop
    while i <= length(lines) - 1
        line = strip(lines[i])
        if occursin("Variable", line)
            i = qe_write_variable(wf, indent + 2, lines, i)
        end
        i += 1
    end
    return writefbodyline(wf, indent + 1, "]),")
end
function write_QECalculationInfo(wf, filename, indent, exec)
    writefbodyline(wf, indent, """QECalculationInfo("$exec", [""")
    allcontrolwritten = false
    anydatawritten = false
    open(filename, "r") do f
        while !eof(f)
            line = readline(f)
            if occursin("NAMELIST", line)
                write_QEControlBlockInfo(wf, indent + 1,
                                         read_block(f, line, "END OF NAMELIST"))
            elseif occursin("CARD:", line)
                anydatawritten = true
                if !allcontrolwritten
                    writefbodyline(wf, indent + 1, "],")
                    writefbodyline(wf, indent + 1, "[")
                    allcontrolwritten = true
                end
                write_QEDataBlockInfo(wf, indent + 1, read_block(f, line, "END OF CARD"))
            end
        end
    end
    if !anydatawritten
        if !allcontrolwritten
            writefbodyline(wf, indent + 1, "],")
        end
        writefbodyline(wf, indent + 1, "QEDataBlockInfo[]")
    else
        writefbodyline(wf, indent + 1, "]")
    end
    return writefbodyline(wf, indent, "),")
end

searchdir(path::String, key) = filter(x -> occursin(key, x), readdir(path))
open(joinpath(@__DIR__, "qeflags.jl"), "w") do wf
    write(wf, "_QEINPUTINFOS() = QECalculationInfo[\n")
    calculation_files = searchdir(joinpath(@__DIR__, "..", "assets", "calculations", "qe"),
                                  "INPUT")
    filepaths = joinpath.(Ref(joinpath(@__DIR__, "..", "assets", "calculations", "qe")),
                          calculation_files)
    for _f in filepaths
        exec_name = join([lowercase(splitext(split(_f, "_")[end])[1]), ".x"], "")
        write_QECalculationInfo(wf, _f, 1, exec_name)
    end
    return write(wf, "]")
end

open(joinpath(@__DIR__, "abinitflags.jl"), "w") do wf
    flaglines = split.(readlines(joinpath(@__DIR__, "..", "assets", "calculations",
                                          "abinit", "calculation_variables.txt")))
    write(wf, "_ABINITFLAGS() = Dict{Symbol, Type}(\n")
    for (fl, typ) in flaglines
        write(wf, "\t:$fl => $typ,\n")
    end
    return write(wf, ")")
end

#all this so I can be lazy and just use repr
struct ElkFlagInfo{T}
    name::Symbol
    default::Union{T,Nothing}  #check again v0.7 Some
    description::String
end
ElkFlagInfo() = ElkFlagInfo{Nothing}(:error, "")
Base.eltype(x::ElkFlagInfo{T}) where {T} = T

struct ElkControlBlockInfo
    name::Symbol
    flags::Vector{<:ElkFlagInfo}
    description::String
end

function closing(c::Char)
    if c == '{'
        return '}'
    elseif c == '('
        return ')'
    elseif c == '['
        return ']'
    end
end

function readuntil_closing(io::IO, opening_char::Char)
    out = ""
    counter = 1
    closing_char = closing(opening_char)
    while counter > 0
        c = read(io, Char)
        out *= c
        if c == opening_char
            counter += 1
        elseif c == closing_char
            counter -= 1
        end
    end
    return out
end

function elk2julia_type(s::AbstractString)
    if s == "integer"
        return Int
    elseif s == "real"
        return Float64
    elseif s == "logical"
        return Bool
    elseif s == "string"
        return String
    elseif s == "complex"
        return ComplexF64
    end
end

function replace_multiple(str, replacements::Pair{String,String}...)
    tstr = deepcopy(str)
    for r in replacements
        tstr = replace(tstr, r)
    end
    return tstr
end

function parse_elk_default(::Type{T}, s) where {T}
    s_cleaned = strip(strip(strip(replace_multiple(s, "{" => "", "}" => "", "." => "",
                                                   " " => ""), '$'), ')'), '(')
    if occursin("-", s_cleaned)
        return nothing
    elseif occursin(",", s_cleaned)
        s_t = split(s_cleaned, ',')
        return parse.(eltype(T), s_t)
    else
        return parse(T, s_cleaned)
    end
end

function blockflags(s::String)
    flag_regex = r"\{([^.}]*)\}"
    flags = ElkFlagInfo[]
    lines = split(s, "\\\\")
    for l in lines
        sline = strip_split(l, '&')
        flag = match(flag_regex, sline[1])
        descr = sline[2]
        flag_i = split(replace_multiple(sline[3], "(" => " ", ")" => " "))
        typ = length(flag_i) > 1 && flag_i[1] != "string" ?
              Vec{parse(Int, flag_i[2]),elk2julia_type(flag_i[1])} :
              elk2julia_type(flag_i[1])
        #TODO default
        default = typ == String ? sline[4] : parse_elk_default(typ, sline[4])
        flagname = Symbol(strip_split(flag.captures[1], '(')[1])
        if flagname == :wann_bands
            push!(flags, ElkFlagInfo{UnitRange{Int}}(:wann_bands, nothing, string(descr)))
        elseif flagname == :wann_projections
            push!(flags,
                  ElkFlagInfo{Vector{String}}(:wann_projections, nothing, string(descr)))
        elseif flagname == :wann_seedname
            push!(flags, ElkFlagInfo{Symbol}(:wann_seedname, nothing, string(descr)))
        else
            push!(flags, ElkFlagInfo{typ}(flagname, default, string(descr)))
        end
    end
    return flags
end

function replace_latex_symbols(s::String)
    return replace_multiple(s, "\\_" => "_", "\\hline" => "", "\\tt " => "", "^" => "e")
end

function read_elk_doc(fn::String)
    blocks = ElkControlBlockInfo[]
    blockname_regex = r"\{(.*)\}"
    open(fn, "r") do f
        while !eof(f)
            line = readline(f)
            if occursin("\\block", line)
                m = match(blockname_regex, line)
                blockname = Symbol(m.captures[1])
                block_text = replace_latex_symbols(readuntil_closing(f, '{'))
                flags = blockflags(block_text)
                push!(blocks,
                      ElkControlBlockInfo(blockname, flags, readuntil(f, "\\block")))
                seek(f, Base.position(f) - 6) #\\block offset
            end
        end
    end
    return blocks
end

ELK_CONTROLBLOCKS = read_elk_doc(joinpath(@__DIR__, "..", "assets", "calculations", "elk",
                                          "elk.txt"))

open(joinpath(@__DIR__, "elkflags.jl"), "w") do f
    return write(f, "_ELK_CONTROLBLOCKS() = $(repr(ELK_CONTROLBLOCKS))")
end


text: f(l,s=length(l))=replace("x z typing...","x"=>l[1],"z"=>"$(s>1 ? "and $(s>2 ? "$s others" : l[2]) are" : "is")")

max_stars_repo_name: pgrepds/Codegolf

max_issues_repo_name: pgrepds/Codegolf

max_forks_repo_name: pgrepds/Codegolf


text:
module Test
export generate_dot, generate_sum

using LinearAlgebra, Random
using ..EFT: two_prod

"""
    (x, y, d, c) = generate_dot(n, c)

Generate two Float64 vectors whose dot product is ill-conditioned.

Inputs:
  n   -- vectors size
  c   -- target condition number (must be >1)
  rng -- (pseudo-)random number generator

Results:
  x, y -- vectors of size n
  d    -- accurate dot product, rounded to nearest
  c    -- actual condition number of the dot product
"""
function generate_dot(n, c::T; rng=Random.GLOBAL_RNG) where {T}
    generate(c, 100, "dot product") do c1
        generate_dot_(n, T(c1), rng)
    end
end

"""
    (x, s, c) = generate_sum(n, c)

Generate a Float64 vectors whose sum is ill-conditioned.

Inputs:
  n   -- vectors size
  c   -- target condition number (must be >1)
  rng -- (pseudo-)random number generator

Results:
  x -- vector of size n
  s -- accurate sum, rounded to nearest
  c -- actual condition number of the sum
"""
function generate_sum(n, c::T; rng=Random.GLOBAL_RNG) where {T}
    generate(c, 10, "sum") do c1
        generate_sum_(n, T(c1), rng)
    end
end

function generate(f, c, cmin, title)
    c1 = c
    c = max(c, cmin)
    for i in 1:100
        res = f(c1)
        c_ = last(res)

        # println("$i -> $c1 \t $c_")

        if     c_ > 3*c
            c1 = max(1.01, 0.8*c1)
        elseif c_ < c/3
            c1 *= 1.1
        else
            return res
        end
    end
    @error "Could not generate $title with requested condition number"
end

function generate_dot_(n, c::T, rng) where {T}
    R = Rational{BigInt}

    # Initialization
    x = zeros(T, n)
    y = zeros(T, n)

    # First half of the vectors:
    #   random numbers within a large exponent range
    n2 = div(n, 2)
    b = log2(c)

    e = rand(rng, n2) .* b/2
    e[1]  = b/2 + 1           # Make sure exponents b/2
    e[n2] = 0                 # and 0 actually occur
    for i in 1:n2
        x[i] = (2*rand(rng, T)-1) * 2^(e[i])
        y[i] = (2*rand(rng, T)-1) * 2^(e[i])
    end

    # Second half of the vectors such that
    #   (*) log2( dot (x[1:i], y[1:i]) ) decreases from b/2 to 0
    δe = -b/(2*(n-n2-1))
    e = b/2:δe:0
    for i in eachindex(e)
        # Random x[i]
        cx = (2*rand(rng)-1) * 2^(e[i])
        x[i+n2] = cx

        # y[i] chosen according to (*)
        cy = (2*rand(rng)-1) * 2^(e[i])
        y[i+n2] = (cy - T(dot(R.(x), R.(y)))) / cx
    end

    # Random permutation of x and y
    perm = randperm(rng, n)
    X = x[perm]
    Y = y[perm]

    # Dot product, rounded to nearest
    d = T(dot(R.(X), R.(Y)))

    # Actual condition number
    c_ = 2 * dot(abs.(X), abs.(Y)) / abs(d)

    (X,Y,d,c_)
end

function generate_sum_(n, c::T, rng) where {T}
    R = Rational{BigInt}

    (x, y, _, _) = generate_dot_(n÷2, c, rng)

    z = (two_prod.(x, y)
         |> Iterators.flatten
         |> collect)

    # Complete if necessary
    if length(z) < n
        push!(z, rand(rng))
    end

    z = shuffle(z)

    # Sum, rounded to nearest
    s = T(sum(R.(z)))

    # Actual condition number
    c_ = sum(abs.(z)) / abs(s)

    (z, s, c_)
end

end


text:
import Base: *, \

# Identity operator

struct Identity end

(*)(::Identity,s::GridData) = s

# Lots of other operators

include("innerproducts.jl")
include("convolution.jl")
include("lgf.jl")
include("lgf-helmholtz.jl")
include("ddf.jl")
include("laplacian.jl")
include("helmholtz.jl")
include("intfact.jl")
include("diffcalculus.jl")
include("shift.jl")
include("nlcalculus.jl")
include("regularization.jl")


text:
using test2
using Test

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


text:
include("Graph.jl")
include("Selector.jl")

function testRoleSim(G, k)
    T = time()
    f = RoleSimSelector(G)
    foreach(i -> x = f(i, k), 1 : G.n)
    return time()-T
end

function testStructSim(G, k)
    T = time()
    f = StructSimSelector(G)
    foreach(i -> x = f(i, k), 1 : G.n)
    return time()-T
end

function testForestSim(G, k)
    T = time()
    f = ForestSimSelector(G)
    foreach(i -> x = f(i, k), 1 : G.n)
    return time()-T
end

function testForestSimApprox(G, k)
    T = time()
    f = ForestSimApproxSelector(G)
    foreach(i -> x = f(i, k), 1 : G.n)
    return time()-T
end


text:
using Test

include("hello-world.jl")

@testset "Say Hi!" begin
    @test hello() == "Hello, World!"
end


text:
function mark_stairs!(r)
    side1 = Sud
    side2 = Ost
    h, w = go_to_corner!(r, side1, side2)
    s = West
    steps = 0
    while true
        if s == West
            steps = mark_line!(r, s)
        else
            move_mark_some_steps!(r, s, steps - 1)
            if isborder(r, West)
                break
            end
            move!(r, West)
        end
        if isborder(r, Nord)
            break
        end
        move!(r, Nord)
        s = inverse_side(s)
    end
    go_to_corner!(r, side1, side2)
    move_back(r, inverse_side(side1), inverse_side(side2), h, w)
end

max_stars_repo_name: dmitrii-a-piskunov/mirea-progs

max_issues_repo_name: dmitrii-a-piskunov/mirea-progs

max_forks_repo_name: dmitrii-a-piskunov/mirea-progs


text:
function visualize!(vis, model, q;
	Δt = 0.1, r = 0.1)

	default_background!(vis)
    setobject!(vis["particle"],
		Rect(Vec(0, 0, 0),Vec(2r, 2r, 2r)),
		MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0, 1.0)))

    anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

    for t = 1:length(q)
        MeshCat.atframe(anim, t) do
            settransform!(vis["particle"], Translation(q[t][1:3]...))
        end
    end

    MeshCat.setanimation!(vis, anim)
end