text:
#
# Initial set up. This stuff should mean that this script will run on a fresh install.
#

using Pkg
Pkg.activate(".");
Pkg.instantiate();

#
# Load the various packages on which this script depends.
#

using Revise
using DrWatson, Stheno, BenchmarkTools, PGFPlotsX, ProgressMeter, TemporalGPs, Random,
    DataFrames, Zygote

using DrWatson: @dict, @tagsave

using LinearAlgebra
LinearAlgebra.BLAS.set_num_threads(1)

const exp_dir_name = "single_output_gps"
const data_dir = joinpath(datadir(), exp_dir_name)

build_gp(k_base, σ², l) = GP(σ² * transform(k_base, 1 / l), GPC())

# Naive implementation.
function build(::Val{:naive}, k_base, σ², l, x, σ²_n)
    f = build_gp(k_base, σ², l)
    return f(x, σ²_n)
end

# Basic lgssm implementation.
function build(::Val{:heap}, k_base, σ², l, x, σ²_n)
    f = build_gp(k_base, σ², l)
    return to_sde(f)(x, σ²_n)
end

# Stack-allocated arrays lgssm implementation.
function build(::Val{:stack}, k_base, σ², l, x, σ²_n)
    f = build_gp(k_base, σ², l)
    return to_sde(f, SArrayStorage(Float64))(x, σ²_n)
end

# Basic lgssm implementation.
function build(::Val{:heapF32}, k_base, σ², l, x, σ²_n)
    f = build_gp(k_base, σ², l)
    return to_sde(f, ArrayStorage(Float32))(x, Float32(σ²_n))
end

# Stack-allocated arrays lgssm implementation.
function build(::Val{:stackF32}, k_base, σ², l, x, σ²_n)
    f = build_gp(k_base, σ², l)
    return to_sde(f, SArrayStorage(Float32))(x, Float32(σ²_n))
end

# Generic logpdf computation.
function build_and_logpdf(val::Val, k_base, σ², l, x, σ²_n, y)
    return logpdf(build(val, k_base, σ², l, x, σ²_n), y)
end

# Specify which experiments to run.
tagsave(
    joinpath(datadir(), exp_dir_name, "meta", "settings.bson"),
    Dict(
        :Ns => [
            2, 5, 10,
            20, 50, 100,
            200, 500, 1_000,
            2_000, 5_000, 10_000,
            20_000, 50_000, 100_000,
            200_000, 500_000, 1_000_000,
            2_000_000, 5_000_000, 10_000_000,
        ],
        :kernels => [
            # (k=Matern12Kernel(), sym=:Matern12Kernel, name="Matern12Kernel"),
            # (k=Matern32Kernel(), sym=:Matern32Kernel, name="Matern32Kernel"),
            (k=Matern52Kernel(), sym=:Matern52Kernel, name="Matern52Kernel"),
        ],
        :implementations => [
            (
                name="naive",
                val=Val(:naive),
                colour="blue",
                marker="square*",
            ),
            (
                name="heap",
                val=Val(:heap),
                colour="black",
                marker="*",
            ),
            (
                name="stack",
                val=Val(:stack),
                colour="red",
                marker="triangle*",
            ),
            (
                name="heap-F32",
                val=Val(:heapF32),
                colour="green",
                marker="*",
            ),
            (
                name="stack-F32",
                val=Val(:stackF32),
                colour="purple",
                marker="triangle*",
            ),
        ]
    )
)

# Run to execute all benchmarking experiments specified in settings.bson.
let
    settings = load(joinpath(datadir(), exp_dir_name, "meta", "settings.bson"))
    settings_list = dict_list(settings)
    progress = Progress(length(settings_list))
    for (p, setting) in enumerate(settings_list)

        N = setting[:Ns]
        kernel = setting[:kernels]
        impl = setting[:implementations]

        if N > 10_000 && impl.name == "naive"
            println("skipping naive implementation with N=$N")
            continue
        end

        # Construct data for problem.
        x = range(-5.0; length=N, step=1e-2)
        σ², l, σ²_n = 1.0, 2.3, 0.5
        k = kernel.k
        rng = MersenneTwister(123456)
        # y = rand(rng, build(Val(:stack), k, σ², l, x, σ²_n))
        y = rand(rng, build(impl.val, k, σ², l, x, σ²_n))

        # Generate results including construction of GP.
        results = @benchmark(build_and_logpdf($(impl.val), $k, $σ², $l, $x, $σ²_n, $y))
        grad_results = @benchmark(Zygote.gradient(
            $build_and_logpdf, $(impl.val), $k, $σ², $l, $x, $σ²_n, $y,
        ))

        # Generate results excluding construction of GP.
        fx = build(impl.val, k, σ², l, x, σ²_n)
        no_build_results = @benchmark(logpdf($fx, $y))
        no_build_grad_results = @benchmark(Zygote.gradient(logpdf, $fx, $y))

        # Save results in predictable location.
        @tagsave(
            joinpath(
                datadir(),
                exp_dir_name,
                savename(Dict(:N=>N, :k=>kernel.sym, :impl=>impl.name), "bson"),
            ),
            Dict(
                "setting" => setting,
                "results" => results,
                "grad_results" => grad_results,
                "no_build_results" => no_build_results,
                "no_build_grad_results" => no_build_grad_results,
                "N" => N,
            ),
        )

        next!(progress)
    end
end

# Run to generate figures specified in settings.bson.
let

    # Load + extract experimental settings.
    settings = load(joinpath(data_dir, "meta", "settings.bson"))
    colour_map = Dict(impl.name => impl.colour for impl in settings[:implementations])

    # Load and post-process results.
    results = DataFrame(
        map(eachrow(collect_results(data_dir))) do row
            row = merge(
                copy(row),
                (
                    kernel = row.setting[:kernels].name,
                    impl = row.setting[:implementations].name,
                    total_time = time(row.results) / 1e9,
                    grad_total_time = time(row.grad_results) / 1e9,
                    run_time = time(row.no_build_results) / 1e9,
                    grad_run_time = time(row.no_build_grad_results) / 1e9,
                ),
            )
            row = merge(
                row,
                (
                    build_time = row.total_time - row.run_time,
                    grad_build_time = row.grad_total_time - row.grad_run_time,
                    total_rate = row.total_time / row.N,
                    grad_total_rate = row.grad_total_time / row.N,
                    run_rate = row.run_time / row.N,
                    grad_run_rate = row.grad_run_time / row.N,
                )
            )
            row = merge(
                row,
                (
                    grad_ratio_total = row.grad_total_time / row.total_time,
                    grad_ratio_build = row.grad_build_time / row.build_time,
                    grad_ratio_run = row.grad_run_time / row.run_time,
                    build_fraction = row.build_time / row.total_time,
                    grad_build_fraction = row.grad_build_time / row.grad_total_time,
                ),
            )
            return row
        end
    )

    # List the things to generate plots for.
    plot_keys = [:total_time, :grad_total_time, :run_time, :grad_run_time, :build_time,
        :grad_build_time, :build_fraction, :grad_build_fraction, :grad_ratio_total,
        :grad_ratio_run, :grad_ratio_build, :total_rate, :grad_total_rate, :run_rate,
        :grad_run_rate,
    ]

    plots = by(results, [:kernel, :impl]) do result
        result = sort(result, :N)
        impl, kernel = first(result).impl, first(result).kernel
        return (
            label = "\\label {$(impl)-$(kernel)}",
            legend_entry = LegendEntry(kernel * " -- " * impl),
            Pair.(
                plot_keys,
                map(plot_keys) do key
                    @pgf PlotInc(
                        {
                            color = colour_map[impl],
                            mark_options={fill=colour_map[impl]}
                        },
                        Coordinates(result[!, :N], result[!, key]),
                    )
                end,
            )...,
        )
    end

    #
    # Plot stuff using the resources generated above.
    #

    # Place to store things to plot.
    plots_to_save = Vector{NamedTuple{(:name, :data)}}(undef, 0)

    x_ticks = 10.0 .^ (0:2:8)
    horz_sep = "20pt"

    # Generate total time / grad total time 1 x 2 group plot.
    for (field, grad_field, yticks) in (
        (:total_time, :grad_total_time, 10.0 .^ (-8:2:2)),
        (:run_time, :grad_run_time, 10.0 .^ (-8:2:2)),
        (:build_time, :grad_build_time, 10.0 .^ (-8:2:2)),
        (:total_rate, :grad_total_rate, 10.0 .^ (-8:2:-1)),
        (:run_rate, :grad_run_rate, 10.0 .^ (-8:2:-1)),
    )
        for (xmode, ymode, name_suffix) in (
            ("normal", "log", ""),
            ("log", "log", "_log"),
        )
            push!(
                plots_to_save,
                (
                    name=string(field) * name_suffix,
                    data=@pgf GroupPlot( # Generate group plot
                        {
                            group_style =
                            {
                                group_size="2 by 1",
                                horizontal_sep=horz_sep,
                                yticklabels_at="edge left",
                            },
                        },
                        { # Style for lhs plot
                            xmode=xmode,
                            ymode=ymode,
                            xlabel="N",
                            ylabel="time (s)",
                            # legend_to_name="legend:$(string(field))",
                            xmajorgrids,
                            ymajorgrids,
                            width="0.48\\textwidth",
                            ytick=yticks,
                            ymin=minimum(yticks),
                            ymax=maximum(yticks),
                            xtick=x_ticks,
                            xmin=0,
                            xmax=1e8,
                        },
                        plots[!, field], # lhs plot
                        { # style for rhs plot
                            xmode=xmode,
                            ymode=ymode,
                            xlabel="N",
                            xmajorgrids,
                            ymajorgrids,
                            width="0.48\\textwidth",
                            ytick=yticks,
                            ymin=minimum(yticks),
                            ymax=maximum(yticks),
                            xtick=x_ticks,
                            xmin=0,
                            xmax=1e8,
                            legend_style =
                            {
                                font = "\\tiny",
                                at = Coordinate(1 + 0.1, 1),
                                anchor = "north west",
                                legend_columns = 1,
                            },
                            transpose_legend,
                        },
                        plots[!, grad_field], # rhs plot
                        plots[!, :legend_entry],
                    )
                ),
            )
        end
    end

    # Plot proportion of time spent setting up in both forwards and reverse passes.
    let
        yticks = [0.0, 0.5, 1.0]
        push!(
            plots_to_save,
            (
                name = "build_fraction",
                data = @pgf GroupPlot( # Generate group plot
                    {
                        group_style =
                        {
                            group_size="2 by 1",
                            horizontal_sep=horz_sep,
                            xticklabels_at="edge bottom",
                            yticklabels_at="edge left",
                        },
                    },
                    { # Style for lhs plot
                        xmode="log",
                        ymode="normal",
                        xlabel="N",
                        legend_to_name="legend:build-fraction",
                        legend_columns=3,
                        transpose_legend,
                        xmajorgrids,
                        ymajorgrids,
                        width="0.48\\textwidth",
                        xtick=x_ticks,
                        ytick=yticks,
                        ymin=minimum(yticks),
                        ymax=maximum(yticks),
                    },
                    plots[!, :build_fraction], # lhs plot
                    plots[!, :legend_entry],
                    { # style for rhs plot
                        xmode="log",
                        ymode="normal",
                        xlabel="N",
                        xmajorgrids,
                        ymajorgrids,
                        width="0.48\\textwidth",
                        xtick=x_ticks,
                        ytick=yticks,
                        ymin=minimum(yticks),
                        ymax=maximum(yticks),
                    },
                    plots[!, :grad_build_fraction], # rhs plot
                )
            ),
        )
    end

    # Plot various ratios of logpdf to logpdf + gradient of logpdf evaluations.
    for (field, yticks) in [
        (:grad_ratio_total, 0:5:50),
        (:grad_ratio_run, vcat(0:5:50, 2)),
        (:grad_ratio_build, 0:5:50),
    ]
        push!(
            plots_to_save,
            (
                name=string(field),
                data = @pgf Axis(
                    {
                        xmode="log",
                        ymode="normal",
                        xlabel="N",
                        width="0.45\\textwidth",
                        xmajorgrids,
                        ymajorgrids,
                        xtick=x_ticks,
                        ytick=yticks,
                    },
                    plots[!, field],
                    # legend_entries,
                )
            )
        )
    end

    # Write results plots to results directory.
    plots_dir = joinpath(plotsdir(), exp_dir_name)
    mkpath(plots_dir)
    for plot in plots_to_save
        pgfsave(
            joinpath(plots_dir, plot.name * ".tex"),
            plot.data;
            include_preamble=true,
        )
    end
end


text:
import Base.show

#struct LC{T}
#    map::Dict{T,Float64}
#end

LC{T} = Dict{T,Float64}

function LC(dt::DataType)
     #return LC(Dict{dt,Float64}())
     return LC{dt}()
end

function parseLC(lcs::String)
    debug=0
    #lc = LC{String}(Dict{String,Float64}())
    lc = LC{String}()
    ss = split(lcs,"+")
    for si in ss
        sis= strip(si)
        #sx = split( sis , r"\s+" )
        sx = split( sis , r"[\s\*]+" )
        vx::String = ""
        nx::Float64 = 1.
        if(size(sx,1)==2)
            if(debug>0); println("$(sis) --> $(sx[1]) $(sx[2])") ; end
            nx = float(sx[1])
            vx = strip(sx[2])
        else
            if(debug>0); println("$(sis) --> $(sx[1])") ; end
            vx = strip(sx[1])
        end
        #lci = LC{String}( Dict(vx => nx) )
        lci = LC{String}( vx => nx )
        #lci.map[vx] = nx
        lc = add(lc,lci)
    end
    return lc
end

function add{T}(a::LC{T},b::LC{T})
    #keys_a::Array{T,1} = collect(keys(a.map))
    #keys_b::Array{T,1} = collect(keys(b.map))
    keys_a::Array{T,1} = collect(keys(a))
    keys_b::Array{T,1} = collect(keys(b))

    keys_ab::Array{T,1} = unique([keys_a;keys_b])
    #map_ab  = Dict{T,Float64}()
    map_ab  = LC{T}()
    for x in keys_ab
        #vi::Float64 = ( haskey(a.map,x) ) ? a.map[x] : 0
        #vi         += ( haskey(b.map,x) ) ? b.map[x] : 0
        vi::Float64 = ( haskey(a,x) ) ? a[x] : 0
        vi         += ( haskey(b,x) ) ? b[x] : 0
        map_ab[x] = vi
    end
    #return LC{T}(map_ab)
    return map_ab
end

function add{T}(a::LC{T},b::T)
    return add(a,b,1.0)
end

function add{T}(a::LC{T},b::T,count::Float64)
    lcb = LC{T}( b => count )
    return add(a,lcb)
end

function mult{T}(a::LC{T},b::Float64)
    lc_new = LC{T}()
    for ki in keys(a)
        lc_new[ki] = b * a[ki]
    end
    return lc_new
end

"""
LinConstr has fields constr::Char, lhs::LC{String} and rhs::LC{String}
and represents a linear constraint
.constr: can be "=", "<" , ">"
convention: "1" in lhs/rhs is the key for the non-variable-multiplied value
"""
struct LinConstr
    constr::Char
    lhs::LC{String}
    rhs::LC{String}
end

function Base.show(io::IO,lc::LC{String})
    cnt=0
    for ki in keys(lc)
        cnt+=1
        #if(lc[ki]==0)
            # skip entry
        #else
            #print( io , @sprintf("%5.3d",lc[ki]) )
            print( io , lc[ki])
            print( io , string( " ",ki ) )
            if(cnt<length(lc)); print(io," + "); end
        #end
    end
end
function Base.show(io::IO,lc::LinConstr)
    show(io,lc.lhs)
    print(io,string(" ",lc.constr," "))
    show(io,lc.rhs)
end

function resolve_lc_s(n::Int64,d_s::Dict{String,Int64},lc::LC{String})
    lcrow::Array{Float64,1} = zeros(n)
    for ei in keys(lc)
        if(!haskey(d_s,ei)); error("LC::resolve_lc : error : no element $(ei)"); end
        lcrow[d_s[ei]] = lc[ei]
    end
    return lcrow
end

function resolve_lc{T}(n::Int64,d_s::Dict{T,Int64},lc::LC{T})
    lcrow::Array{Float64,1} = zeros(n)
    for ei in keys(lc)
        if(!haskey(d_s,ei)); error("LC::resolve_lc : error : no element $(ei)"); end
        lcrow[d_s[ei]] = lc[ei]
    end
    return lcrow
end

"""
parses linear constraints
comparison operators are: '=' , '<' , '>'
on left and right hand side there must be strings which can be parse
by parseLC
"""
function parse_linconstr(str)
    iseq = contains(str,"=")
    isle = contains(str,"<")
    isge = contains(str,">")

    n_elg = (iseq?1:0)+(isle?1:0)+(isge?1:0)
    if(n_elg!=1); error(string("Unparsable constraint: ",str)); end

    op::Char = '?'

    if(iseq) ; op ='=' ; end
    if(isle) ; op ='<' ; end
    if(isge) ; op ='>' ; end

    str_ab = split(str,string(op))
    a = String( str_ab[1] )
    b = String( str_ab[2] )

    #lhs = parseLC(a)
    #rhs = parseLC(b)
    lhs = parse_lc_expr(strip(a))
    rhs = parse_lc_expr(strip(b))

    return LinConstr(op,lhs,rhs)
end

"""
parse_lc_expr(str::String)
returns LC{String}

parses formulas which are sums like:
" 5 va + 4 bxy - 2 + 1 va "
" 3 "
" 4 v "
" 3 * x + 2 * abc + 3"
"""
function parse_lc_expr(str::String)
    # enry "1" maps to the "not-multiplied-with-var" part
    val = LC{String}("1" => 0)

    str = strip(str);
    if(isempty(str)); return val; end

    if(contains(str,"+"))
        split_pos = search(str,"+")[1]
        va = parse_lc_expr( String( str[1:split_pos-1] ) )
        vb = parse_lc_expr( String( str[split_pos+1:length(str)] ) )
        return add(va,vb)
    else
        if(contains(str,"-"))
            split_pos = search(str,"-")[1]
            va = parse_lc_expr( String( str[1:split_pos-1] ) )
            vb = parse_lc_expr( String( str[split_pos+1:length(str)] ) )
            return add(va, mult( vb , -1.0 ) )
        end
    end

    ## now down to parsing literals..
    #str_s = split(str,"\s+")
    str_s = split(str,r"[\s\*]+")
    if(length(str_s)==1)
        str_s = strip(string(str_s[1]))
        #str_s = String(str_s)
        #rmatch = ismatch(r"\d*\.?\d+",str_s)
        parsed_float = false
        try
            pf = parse(Float64,str_s)
            val["1"] = pf
            parsed_float = true
        end
        if(!parsed_float)
            val[str_s] = 1.0
        end
        return val
    else
        if(length(str_s)>2); error("Error parsing constraint: "+str); end
        v_num = parse(Float64,str_s[1])
        v_var = strip(String(str_s[2]))
        val[v_var] = v_num
        return val
    end
end

ext: jl


text:
using Sherlogs, PyPlot, StatsBase, Printf

A = Sherlog64.(rand(1000,1000))
b = Sherlog64.(rand(1000))

# trigger logging of A,b
A .+ 0
b .+ 0
lb0 = get_logbook()

# LU decomposition
x = A\b
lb1 = get_logbook()

H0 = entropy(lb0.logbook/sum(lb0),2)
H1 = entropy(lb1.logbook/sum(lb1),2)

## PLOT
H0s = @sprintf "%.2f" H0
H1s = @sprintf "%.2f" H1

f = [1/512,1/4,1,4,512,Inf,-1/512,-1/4,-1,-4,-512,-Inf]
fs = ["⅟512","¼","1","4","512","∞","-⅟512","-¼","-1","-4","-512","-∞"]
fi = Int.(reinterpret.(UInt16,Float16.(f)))

fig,ax = subplots(1,1,figsize=(8,3))

ax.plot(lb0.logbook/sum(lb0)*100,label="U(0,1), H=$H0s bit")
ax.plot(lb1.logbook/sum(lb1)*100,label="A\\b,      H=$H1s bit")
legend()
ax.set_title("Float16 bitpattern histogram",loc="left")
ax.set_ylabel("[%]")
ax.set_xticks(fi)
ax.set_xticklabels(fs)
ax.set_xlim(0,2^16)

tight_layout()
savefig("matrixsolve.png")
close(fig)


text:
#
# Function to perform the most important selections that are used in solute-solvent
# analysis
#

using Parameters

# Main function: receives the atoms vector and a julia function to select

function select(set::AbstractVector{T}; by = all) where {T}
    selected = similar(set, 0)
    for el in set
        if by(el)
            push!(selected, el)
        end
    end
    return selected
end

# Given a selection string

function select(set::AbstractVector{T}, selection::String) where {T}
    query = parse_query(selection)
    return select(set, by = el -> apply_query(query, el))
end

#
# Return indexes only
#
function selindex(set::AbstractVector{T}; by = all) where {T}
    indexes = Vector{Int}(undef, 0)
    for i in eachindex(set)
        if by(set[i])
            push!(indexes, i)
        end
    end
    return indexes
end

function selindex(set::AbstractVector{T}, selection::String) where {T}
    query = parse_query(selection)
    return selindex(set, by = el -> apply_query(query, el))
end

# Comparison operators

const operators = (
    "=" => (x, y) -> isequal(x, y),
    "<" => (x, y) -> isless(x, y),
    ">" => (x, y) -> isless(y, x),
    "<=" => (x, y) -> (!isless(y, x)),
    ">=" => (x, y) -> (!isless(x, y)),
)

#
# Keywords
#

struct Keyword
    ValueType::Type
    name::String
    field::Symbol
    operators::Tuple
end

function (key::Keyword)(s::AbstractVector{<:AbstractString})
    @unpack name, field, operators = key
    for op in operators
        if (i = has_key(op.first, s)) > 0
            return el -> op.second(getfield(el, field), parse_to_type(key, s[i+1]))
        end
    end
    # If no operator was found, assume that `=` was intended
    i = has_key(name, s)
    return el -> isequal(getfield(el, field), parse_to_type(key, s[i+1]))
end

"""
```
FunctionalKeyword{T}
```

This is a structure that will store a keyword that depends on an external function
requiring an operator and an argument.

## Example:

```
element_keyword = FunctionalKeyword(String,
                                    "element",
                                    element,
                                    ("=",isequal))

```
will define a keyword "element" to be used as `element C`, which will return
`true` if there is an `element` function such that `element(atom) == C`.

"""
struct FunctionalKeyword{FunctionType}
    ValueType::Type
    name::String
    by::FunctionType
    operators::Tuple
end

function (key::FunctionalKeyword)(s::AbstractVector{<:AbstractString})
    @unpack name, by, operators = key
    for op in operators
        if (i = has_key(op.first, s)) > 0
            return el -> op.second(by(el), parse_to_type(key, s[i+1]))
        end
    end
    # If no operator was found, assume that `=` was intended
    i = has_key(name, s)
    return el -> isequal(by(el), parse_to_type(key, s[i+1]))
end

#
# Macro keywords (functions without parameters)
#

struct MacroKeyword{FunctionType}
    name::String
    by::FunctionType
end

function (key::MacroKeyword)(s::AbstractVector{<:AbstractString})
    return key.by
end

"""

`parse_to_type(key::Keyword, val::String)`

Tries to parse `val` into the type of value expected by `key.ValueType`.

"""
function parse_to_type(key::Union{Keyword,FunctionalKeyword}, val)
    if key.ValueType == String
        return val
    end
    try
        val = parse(key.ValueType, val)
        return val
    catch
        parse_error(
            "Could not parse $val for keyword $(key.name), expected $(key.ValueType)",
        )
    end
end

#
# Keywords for PDBTools
#

keywords = [
    Keyword(Int, "index", :index, operators),
    Keyword(Int, "index_pdb", :index_pdb, operators),
    Keyword(Int, "resnum", :resnum, operators),
    Keyword(Int, "residue", :residue, operators),
    Keyword(Float64, "beta", :beta, operators),
    Keyword(Float64, "occup", :occup, operators),
    Keyword(Int, "model", :model, operators),
    Keyword(String, "name", :name, operators),
    Keyword(String, "segname", :segname, operators),
    Keyword(String, "resname", :resname, operators),
    Keyword(String, "chain", :chain, operators),
]

macro_keywords = [
    MacroKeyword("water", iswater),
    MacroKeyword("protein", isprotein),
    MacroKeyword("polar", ispolar),
    MacroKeyword("nonpolar", isnonpolar),
    MacroKeyword("basic", isbasic),
    MacroKeyword("acidic", isacidic),
    MacroKeyword("charged", ischarged),
    MacroKeyword("aliphatic", isaliphatic),
    MacroKeyword("aromatic", isaromatic),
    MacroKeyword("hydrophobic", ishydrophobic),
    MacroKeyword("neutral", isneutral),
    MacroKeyword("backbone", isbackbone),
    MacroKeyword("sidechain", issidechain),
    MacroKeyword("all", all),
]

functional_keywords = [FunctionalKeyword(String, "element", element, operators)]

#
# parse_query and apply_query are a very gentle contribution given by
# CameronBieganek in https://discourse.julialang.org/t/parsing-selection-syntax/43632/9
# while explaining to me how to creat a syntex interpreter
#
"""
```
has_key(key::String, s::AbstractVector{<:AbstractString})
```

Returns the first index of the vector `s` in which where `key` is found, or 0.

## Example:

```julia-repl

julia> PDBTools.has_key("or",["name","CA","or","index","1"])
3

julia> PDBTools.has_key("and",["name","CA","or","index","1"])
0

```

"""
function has_key(key::String, s::AbstractVector{<:AbstractString})
    i = findfirst(isequal(key), s)
    if i == nothing
        0
    else
        i
    end
end

"""
```
parse_query(selection:String)
```

Calls `parse_query_vector` after splitting the selection string.

"""
parse_query(selection::String) = parse_query_vector(split(selection))

"""

```
parse_query_vector(s::AbstractVector{<:AbstractString})
```

"""
function parse_query_vector(s)
    if (i = has_key("or", s)) > 0
        deleteat!(s, i)
        (|, parse_query_vector.((s[1:i-1], s[i:end]))...)
    elseif (i = has_key("and", s)) > 0
        deleteat!(s, i)
        (&, parse_query_vector.((s[1:i-1], s[i:end]))...)
    elseif (i = has_key("not", s)) > 0
        deleteat!(s, i)
        (!, parse_query_vector(s[i:end]))

        # keywords
    else
        for key in keywords
            if (i = has_key(key.name, s)) > 0
                deleteat!(s, i)
                return key(s)
            end
        end
        for key in macro_keywords
            if (i = has_key(key.name, s)) > 0
                deleteat!(s, i)
                return key(s)
            end
        end
        for key in functional_keywords
            if (i = has_key(key.name, s)) > 0
                deleteat!(s, i)
                return key(s)
            end
        end
        parse_error("Unable to parse selection string.")
    end
end

function apply_query(q, a)
    if !(q isa Tuple)
        q(a)
    else
        f, args = Iterators.peel(q)
        f(apply_query.(args, Ref(a))...)
    end
end

#
# Simple error message (https://discourse.julialang.org/t/a-julia-equivalent-to-rs-stop/36568/13)
#

struct NoBackTraceException
    exc::Exception
end

function Base.showerror(io::IO, ex::NoBackTraceException, bt; backtrace = true)
    Base.with_output_color(get(io, :color, false) ? Base.error_color() : :nothing, io) do io
        showerror(io, ex.exc)
    end
end

parse_error(str) = throw(NoBackTraceException(ErrorException(str)))

max_stars_repo_name: mcubeg/PDBTools

max_issues_repo_name: mcubeg/PDBTools

max_forks_repo_name: mcubeg/PDBTools


text:
__precompile__()
module Day3_module

export get_positions

dirs = Dict(
    'U' => [0, +1],
    'D' => [0, -1],
    'R' => [+1, 0],
    'L' => [-1, 0],
)

function get_positions(steps)
print(typeof(steps))
    current = [0, 0]
    pos = Array{Tuple{Int64, Int64}}(undef, 0)
    for next in steps
        direction = dirs[next[1]]
        distance = parse(Int64, next[2:end])
        for i=1:distance
            current += direction
            push!(pos, tuple(current...))
        end
    end
    pos
end

precompile(get_positions, (Array{SubString{String},1},))

end

max_stars_repo_name: st--/advent-of-code

max_issues_repo_name: st--/advent-of-code

max_forks_repo_name: st--/advent-of-code


text:
using BinaryBuilder, Dates, Pkg, Base.BinaryPlatforms
include("../common.jl")

# Don't mount any shards that you don't need to
Core.eval(BinaryBuilder, :(bootstrap_list = [:rootfs]))

compiler_target = try
    parse(Platform, ARGS[end])
catch
    error("This is not a typical build_tarballs.jl!  Must provide exactly one platform as the last argument!")
end
deleteat!(ARGS, length(ARGS))
name = "PlatformSupport"
version = VersionNumber("$(year(today())).$(month(today())).$(day(today()))")

# Add check that FreeBSD and MacOS shards have version numbers embedded
if (Sys.isapple(compiler_target) || Sys.isfreebsd(compiler_target)) && os_version(compiler_target) === nothing
    error("Compiler target $(triplet(compiler_target)) does not have an `os_version` tag!")
end

sources = [
    ArchiveSource("https://mirrors.edge.kernel.org/pub/linux/kernel/v4.x/linux-4.20.9.tar.xz",
                  "b5de28fd594a01edacd06e53491ad0890293e5fbf98329346426cf6030ef1ea6"),
    ArchiveSource("https://sourceforge.net/projects/mingw-w64/files/mingw-w64/mingw-w64-release/mingw-w64-v7.0.0.tar.bz2",
                  "aa20dfff3596f08a7f427aab74315a6cb80c2b086b4a107ed35af02f9496b628"),
    ArchiveSource("https://download.freebsd.org/ftp/releases/amd64/11.4-RELEASE/base.txz",
                  "3bac8257bdd5e5b071f7b80cc591ebecd01b9314ca7839a2903096cbf82169f9"),
    ArchiveSource("https://github.com/llvm/llvm-project/releases/download/llvmorg-8.0.1/libcxx-8.0.1.src.tar.xz",
                  "7f0652c86a0307a250b5741ab6e82bb10766fb6f2b5a5602a63f30337e629b78"),
    DirectorySource("./bundled"),
]

macos_sdk = if Sys.isapple(compiler_target) && arch(compiler_target) == "aarch64"
    # MacOSX11.1 is not yet available on phracker/MacOSX-SDKs
    # https://github.com/phracker/MacOSX-SDKs/pull/32#issuecomment-749230532
    ArchiveSource("https://github.com/larskanis/MacOSX-SDKs/releases/download/11.1/MacOSX11.1.sdk.tar.xz",
                  "97f44b22949cea4522408ccca9a8d87f2d09779b2878423d2d7a2cb805c3d42d")
else
    ArchiveSource("https://github.com/phracker/MacOSX-SDKs/releases/download/10.15/MacOSX10.12.sdk.tar.xz",
                  "6852728af94399193599a55d00ae9c4a900925b6431534a3816496b354926774")
end

push!(sources, macos_sdk)

script = "COMPILER_TARGET=$(BinaryBuilder.aatriplet(compiler_target))\n" * raw"""
## Function to take in a target such as `aarch64-linux-gnu`` and spit out a
## linux kernel arch like "arm64".
target_to_linux_arch()
{
    case "$1" in
        arm*)
            echo "arm"
            ;;
        aarch64*)
            echo "arm64"
            ;;
        powerpc*)
            echo "powerpc"
            ;;
        i686*)
            echo "x86"
            ;;
        x86*)
            echo "x86"
            ;;
    esac
}

## sysroot is where most of this stuff gets plopped
sysroot=${prefix}/${COMPILER_TARGET}/sys-root

# Install kernel headers
case "${COMPILER_TARGET}" in
    *-linux-*)
        cd $WORKSPACE/srcdir/linux-*/

        # Grumble, grumble, need gcc just to install some headers...
        apk add gcc musl-dev

        # The kernel make system can't deal with spaces (for things like ccache) very well
        KERNEL_FLAGS="ARCH=$(target_to_linux_arch ${COMPILER_TARGET}) -j${nproc}"
        make ${KERNEL_FLAGS} mrproper V=1
        make ${KERNEL_FLAGS} headers_check V=1
        make ${KERNEL_FLAGS} INSTALL_HDR_PATH=${sysroot}/usr V=1 headers_install

        # Move case-sensitivity issues, breaking netfilter without a patch
        NF="${prefix}/${COMPILER_TARGET}/sys-root/usr/include/linux/netfilter"
        for NAME in CONNMARK DSCP MARK RATEEST TCPMSS; do
            mv "${NF}/xt_${NAME}.h" "${NF}/xt_${NAME}_.h"
        done

        for NAME in ECN TTL; do
            mv "${NF}_ipv4/ipt_${NAME}.h" "${NF}_ipv4/ipt_${NAME}_.h"
        done
        mv "${NF}_ipv6/ip6t_HL.h" "${NF}_ipv6/ip6t_HL_.h"
        ;;

    *-mingw*)
        cd $WORKSPACE/srcdir/mingw-*/mingw-w64-headers
        ./configure --prefix=/ \
            --enable-sdk=all \
            --enable-secure-api \
            --host=${COMPILER_TARGET}

        make install DESTDIR=${sysroot}
        ;;

    *-freebsd*)
        mkdir -p "${sysroot}/usr"
        mv usr/include "${sysroot}/"
        ln -sf "../include" "${sysroot}/usr/include"
        ;;

    *-apple-*)
        cd ${WORKSPACE}/srcdir/MacOSX*.sdk
        mkdir -p "${sysroot}/usr"
        mv usr/include "${sysroot}/usr"
        mv System "${sysroot}/"

        # Grumble, grumble, need gcc just to install some headers...
        apk add gcc g++ musl-dev

        # Also deploy libcxx headers
        cd ${WORKSPACE}/srcdir/libcxx*
        mkdir build && cd build
        PREFIX="${sysroot}" cmake .. -DLLVM_ENABLE_PROJECTS='libcxx' \
                 -DCMAKE_INSTALL_PREFIX="${sysroot}/usr" \
                 -DCMAKE_CROSSCOMPILING=True \
                 -DLLVM_HOST_TRIPLE=${COMPILER_TARGET} \
                 -DDARWIN_macosx_CACHED_SYSROOT:STRING="${sysroot}"
        make install-cxx-headers
        ;;
    *)
        echo "ERROR: Unmatched platform!" >&2
        exit 1
        ;;
esac

# We create a link from ${COMPILER_TARGET}/sys-root/usr/local to ${prefix}.
# This is the most reliable way for our sysroot'ed compilers to find destination
# libraries so far, hopefully this changes in the future.
mkdir -p "${sysroot}/usr"
ln -s "${prefix}" "${sysroot}/usr/local"
"""

# Build the artifacts
ndARGS, deploy, deploy_target = find_deploy_arg(ARGS)
build_info = build_tarballs(ndARGS, "$(name)-$(triplet(compiler_target))", version, sources, script, [host_platform], Product[], []; skip_audit=true)
if deploy
    upload_and_insert_shards(deploy_target, name, version, build_info; target=compiler_target)
end

max_stars_repo_name: allixender/Yggdrasil

max_issues_repo_name: allixender/Yggdrasil

max_forks_repo_name: allixender/Yggdrasil


text:
using BlackBoxOptim

function get_metric(this::DSB2018.Data.Image, res::Matrix{Float64})

    adaptive_threshold = DSB2018.Model.skimage_filters.threshold_local(res, 11, offset=0.0)
    global_threshold   = DSB2018.Model.skimage_filters.threshold_otsu(res)

    w_adaptve = 0.1
    threshold = adaptive_threshold * w_adaptve .+ global_threshold * (1.0 - w_adaptve)

    img_t = res .> threshold
    img_t = DSB2018.Model.skimage_morphology.remove_small_objects(img_t, min_size=10)
    img_t = DSB2018.Model.ndi.binary_fill_holes(img_t)

    labeled_prediction = Images.label_components(img_t)
    prediction_vector  = DSB2018.Model.split_labeled_matrix(labeled_prediction)
    truth_images = getTrainImageMasks(this)

    return metric(prediction_vector, truth_images)
end

ids = DSB2018.Data.getTrainImageIds()
img = DSB2018.Data.loadTrainImage(ids[1])

x̂ = Float64.(Gray.(img.image))
ŷ = Float64.(img.mask .> 0)

function mask(c::T, x::Float64, y::Float64) where {T}
    μx = c[1]
    μy = c[2]
    σx = c[3]
    σy = c[4]
    ρ  = c[5]
    w  = c[6]

    h = (x - μx)^2.0/σx^2. + (y - μy)^2.0/σy^2. - 2.*ρ*(x-μx)*(y-μy)/σx/σy
    w * exp( -h/2./(1.0 - ρ^2.) )
end

function bbpredict(c::Vector{Float64}, sz::Tuple{Int64, Int64})
    τ = Matrix{Float64}(sz)
    r = CartesianRange(CartesianIndex(1, 1), CartesianIndex(size(τ)))
    for iz in r
        x, y = iz[1]/sz[1], iz[2]/sz[2]
        τ[iz] = mask(c, x, y)
    end
    τ
end

function bbpredict!(c::Vector{Float64}, τ::Matrix{Float64})
    sz = size(τ)
    r = CartesianRange(CartesianIndex(1, 1), CartesianIndex(size(τ)))
    for iz in r
        x, y = iz[1]/sz[1], iz[2]/sz[2]
        τ[iz] = mask(c, x, y)
    end
    τ
end

c = [0.5, 0.1, 0.3, 0.6, 0.9, 1.]

function bbloss(c::Vector{Float64}, img::Matrix{Float64})
    out = 0.0
    for iz in CartesianRange(CartesianIndex(1, 1), CartesianIndex(size(img)))
        x, y = iz[1]/size(img)[1], iz[2]/size(img)[2]
        out += (img[iz] - mask(c, x, y)).^2
    end
    out / mean(img).^2
end

sr = Vector{Tuple{Float64, Float64}}(6)
sr[1] = (0, 1)
sr[2] = (0, 1)
sr[3] = (0.001, 0.05)
sr[4] = (0.001, 0.05)
sr[5] = (-0.5, 0.5)
sr[6] = (0.001, 2.0)

res = bboptimize(x -> bbloss([c0[1:2]; x], x̂),
                NumDimensions = 4,
                SearchRange = sr,
                MaxSteps = 100000,
                PopulationSize = 10,
                Method = :random_search)

c0 = best_candidate(res)
imshow(bbpredict(c0, size(ŷ)))
imshow(x̂)


text:
using Colorsys
using Test

@testset "Colorsys.jl" begin
    @testset "_v" begin
        @test Colorsys._v(1,2,0.1) ≈ 1.6
        @test Colorsys._v(1,2,0.4) == 2
        @test Colorsys._v(1,2,0.6) ≈ 1.4
        @test Colorsys._v(1,2,0.8) == 1
    end
    y,i,q = rgb2yiq(120, 130, 140)
    @test y≈0.5023529411 && i≈-0.0360784313 && q≈0.0039215686
    r,g,b = yiq2rgb(0.8, 0.6, 0.4)
    @test (r,g,b) == (255,96,211)
    @testset "rgb2hls" begin
        h,l,s = rgb2hls(120, 130, 140)
        @test h≈210.0000000000 && l≈0.5098039215 && s≈0.0800000000
        h,l,s = rgb2hls(120, 120, 120)
        @test h≈0.0000000000 && l≈0.4705882352 && s≈0.0000000000
    end
    @testset "hls2rgb" begin
        r,g,b = hls2rgb(120, 0.6, 0.4)
        @test (r,g,b) == (112,193,112)
        r,g,b = hls2rgb(120, 0.6, 0.0)
        @test (r,g,b) == (0,0,0)
    end
    @testset "rgb2hsv" begin
        h,s,v = rgb2hsv(120, 130, 140)
        @test h≈210.0000000000 && s≈0.1428571428 && v≈0.5490196078
        h,s,v = rgb2hsv(120, 120, 120)
        @test h≈0.0 && s≈0.0 && v≈0.4705882352
        h,s,v = rgb2hsv(220, 130, 240)
        @test h≈289.0909090909 && s≈0.4583333333 && v≈0.9411764705
    end
    @testset "hsv2rgb" begin
        r,g,b = hsv2rgb(40, 0.6, 0.4)
        @test (r,g,b) == (102,81,40)
        r,g,b = hsv2rgb(80, 0.6, 0.4)
        @test (r,g,b) == (81,102,40)
        r,g,b = hsv2rgb(120, 0.6, 0.4)
        @test (r,g,b) == (40,102,40)
        r,g,b = hsv2rgb(200, 0.6, 0.4)
        @test (r,g,b) == (40,81,102)
        r,g,b = hsv2rgb(260, 0.6, 0.4)
        @test (r,g,b) == (61,40,102)
        r,g,b = hsv2rgb(320, 0.6, 0.4)
        @test (r,g,b) == (102,40,81)
        r,g,b = hsv2rgb(120, 0, 0.2)
        @test (r,g,b) == (51,51,51)
    end
    h,l,s = yiq2hls(0.8, 0.6, 0.4)
    @test h≈316.6037735849 && l≈0.6882352941 && s≈1.0000000000
    h,s,v = yiq2hsv(0.8, 0.6, 0.4)
    @test h≈316.6037735849 && s≈0.6235294117 && v≈1.0000000000
    y,i,q = hls2yiq(120, 0.6, 0.4)
    @test y≈0.6266274509 && i≈-0.0889411764 && q≈-0.1651764705
    h,s,v = hls2hsv(120, 0.6, 0.4)
    @test h≈120.0000000000 && s≈0.4196891191 && v≈0.7568627450
    y,i,q = hsv2yiq(120, 0.6, 0.4)
    @test y≈0.3003137254 && i≈-0.0680784313 && q≈-0.1264313725
    h,l,s = hsv2hls(120, 0.6, 0.4)
    @test h≈120.0000000000 && l≈0.2784313725 && s≈0.4366197183
end


text:
include("colorscheme.jl")
using PyPlot, DrWatson
using3D()
export rdspl, axis_zoomin!, add_identifiers!, nice_arrow!, figx, figy

DrWatson._wsave(s, fig::Figure) = fig.savefig(s, dpi = 600, transparent = false)

"""
    rdspl(x, n = 3)
Round `x` with `n` sigdigits for display purposes.
"""
rdspl(x::Real, n = 3) = round(x, sigdigits=n)
rdspl(x::AbstractVector, n = 3) = Tuple((round.(Float64.(x); sigdigits=n)))

PyPlot.rc("lines", lw = 3)
PyPlot.rc("errorbar", capsize = 6)
PyPlot.rc("axes", grid = true)
PyPlot.rc("grid", color = "0.75", alpha = 0.75)
PyPlot.rc("axes", axisbelow=true) # makes grid behind other plotting elements

MATHFONTSIZE = 34
PyPlot.rc("font", size = 28) # set default fontsize
PyPlot.rc("xtick", labelsize = 24)
PyPlot.rc("ytick", labelsize = 24)
PyPlot.rc("axes", labelsize = MATHFONTSIZE)
PyPlot.rc("legend", fontsize = 30, handlelength = 1)
# PyPlot.rc("font", family = "Times New Roman") # Serif main font
PyPlot.rc("font", family = "DejaVu Sans") # sans main font
# PyPlot.rc("mathtext", rm = "sanserif", fontset="dejavusans") # sans math font
PyPlot.rc("mathtext", rm = "serif", fontset="cm") # serif math font

for z in ("x", "y")
    PyPlot.rc("$(z)tick.major", size = 7, width = 1.2)
    PyPlot.rc("$(z)tick.minor", size = 3, visible = false)
end

PyPlot.rc("figure", figsize = (figx, figy))
PyPlot.rc("savefig", dpi = 600, transparent = true, format = "png")

# set default color cycle
PyPlot.rc("axes", prop_cycle = matplotlib.cycler(color=COLORS))

if false # test font
    figure(figsize=(10,8)); plot(rand(10), label = "\$a=\\int_0^\\infty xdx\$")
    xlabel("x label"); ylabel("\$H_2\$"); legend(); tight_layout()
end

if false # test color scheme
    fig = figure(figsize = (20, 15)) # show colors
    ax1 = subplot(231)
    ax2 = subplot(232)
    ax3 = subplot(233)
    ax4 = subplot(223)
    ax5 = subplot(224)
    lw = 60
    L = length(COLORSCHEME)
    for (i, c) in enumerate(COLORS)
        chsv = matplotlib.colors.rgb_to_hsv(matplotlib.colors.to_rgb(c))
        ax1.plot([0, 1], [0, 0] .+ i, color = c, lw = lw)
        ax1.set_title("color")
        ax2.plot([0, 1], [0, 0] .+ i, color = string(chsv[3]), lw = lw)
        ax2.set_title("brightness")
        ax3.plot([0, 1], [0, 0] .+ i, color = string(chsv[2]), lw = lw)
        ax3.set_title("saturation")
        x = 0:0.05:5π
        ax4.plot(x, rand(1:3) .* cos.(x .+ i/2) .+ rand(length(x))/5; color=c, lw = 2)
        ax5.bar(collect(1:4) .+ (i-1)/L, 0.5rand(4) .+ 0.5, 1/L; color=c)
    end
    fig = tight_layout()
end

bbox = Dict(:boxstyle => "round,pad=0.3", :facecolor=>"white", :alpha => 1.0)

"`add_identifiers!(fig = gcf(), axs = fig.get_axes(); xloc, yloc)`"
function add_identifiers!(fig = gcf(), axs = fig.get_axes();
        xloc = 0.975, yloc = 1
    )
    for (i, ax) in enumerate(axs)
        l = collect('a':'z')[i]
        if ax.name ≠ "3d"
            ax.text(xloc, yloc, "$(l)"; transform = ax.transAxes,
            bbox = bbox, zorder = 99, va = "top")
        end
    end
end

"""
    nice_arrow!(ax, xc, yc, xspan, yspan;
        style = "<->", tex = "", xo = 0.2, yo = -0.2, color = "C0"
    )
"""
function nice_arrow!(ax, xc, yc, xspan, yspan;
    style = "<->", tex = "", xo = 0.2, yo = -0.2, color = "C0")

    ax.annotate("";
        xy=(xc-xspan/2, yc - yspan/2), xycoords="data",
        xytext=(xc+xspan/2, yc + yspan/2), textcoords="data",
        arrowprops = Dict(
            :arrowstyle=>style,
            :connectionstyle=>"arc3",
            :lw=>1.5, :color => color),
        zorder = 99, color
    )
    if tex != ""
        ax.text(xc + xo, yc + yo, tex; va = :center, color)
    end
end

"""
    axis_zoomin!(zoomin, origin, zbox, rbox, co = "C0"; kwargs...)
Create a zoomin box connecting two axes, the `zoomin` and `origin`.
The zoom-in box is in the `origin` axis, while the `zoomin` axis is the
box. `rbox` is the enclosing box of the `zoomin` axes, while `zbox`
is the small "zoom-in" box of the `origin` axis. They must be in the form
((x1, y1), (x2, y2)). `co` is color `α` the alpha setting.
"""
function axis_zoomin!(zoomin, origin, zbox, rbox, co = "C0";
    connect_lines = true, lw = 2.0, α = 1.0, dir = :right, set_lims::Bool=true)
    # plot box in zoomin axis
    line, = zoomin.plot(
        [rbox[1][1], rbox[2][1], rbox[2][1], rbox[1][1], rbox[1][1]],
        [rbox[1][2], rbox[1][2], rbox[2][2], rbox[2][2], rbox[1][2]],
        color=co, lw = lw, alpha = α
    )
    line.set_clip_on(false)
    # plot box in origin axis
    line, = origin.plot(
        [zbox[1][1], zbox[2][1], zbox[2][1], zbox[1][1], zbox[1][1]],
        [zbox[1][2], zbox[1][2], zbox[2][2], zbox[2][2], zbox[1][2]],
        color=co, lw = lw, alpha = α
    )
    line.set_clip_on(false)

    # xyA is zoomin axis, xyB is origin axis
    if connect_lines
    if dir == :right
        for e in 1:2
            con = matplotlib.patches.ConnectionPatch(
            xyA = (rbox[1][1], rbox[e][2]), xyB=(zbox[2][1], zbox[e][2]),
            coordsA="data", coordsB="data",
            axesA = zoomin, axesB=origin, color=co, lw = lw, alpha = α)
            con.set_clip_on(false)
            zoomin.add_artist(con)
        end
    elseif dir == :top
        for e in 1:2
            con = matplotlib.patches.ConnectionPatch(
            xyA = (rbox[e][1], rbox[1][2]), xyB=(zbox[e][1], zbox[2][2]),
            coordsA="data", coordsB="data",
            axesA = zoomin, axesB=origin, color=co, lw = lw, alpha = α)
            con.set_clip_on(false)
            zoomin.add_artist(con)
        end
    end
    end
    # Set limits for zoom axis
    if set_lims
        zoomin.set_xlim(zbox[1][1], zbox[2][1])
        zoomin.set_ylim(zbox[1][2], zbox[2][2])
    end
    # remove axis for zoomin
    zoomin.axis("off")
end

function darken_color(color, amount=0.5)
    c = matplotlib.colors.to_rgb(color)
    c = matplotlib.colors.rgb_to_hsv(c)
    c = (c[1], c[2], clamp(amount*c[3], 0, 1))
    return matplotlib.colors.hsv_to_rgb(c)
end

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
# VTK cell definitions

# Definitions copied from the vtkCellType.h file of the VTK source code:
# https://raw.githubusercontent.com/Kitware/VTK/master/Common/DataModel/vtkCellType.h

__precompile__()

module VTKCellTypes

import Compat.UTF8String

immutable VTKCellType
    vtk_name::UTF8String
    vtk_id::UInt8
    nodes::Int
end

# Linear cells
const VTK_EMPTY_CELL       = VTKCellType("VTK_EMPTY_CELL",      UInt8(0),  0)
const VTK_VERTEX           = VTKCellType("VTK_VERTEX",          UInt8(1),  1)
const VTK_POLY_VERTEX      = VTKCellType("VTK_POLY_VERTEX",     UInt8(2), -1)
const VTK_LINE             = VTKCellType("VTK_LINE",            UInt8(3),  2)
const VTK_POLY_LINE        = VTKCellType("VTK_POLY_LINE",       UInt8(4), -1)
const VTK_TRIANGLE         = VTKCellType("VTK_TRIANGLE",        UInt8(5),  3)
const VTK_TRIANGLE_STRIP   = VTKCellType("VTK_TRIANGLE_STRIP",  UInt8(6), -1)
const VTK_POLYGON          = VTKCellType("VTK_POLYGON",         UInt8(7), -1)
const VTK_PIXEL            = VTKCellType("VTK_PIXEL",           UInt8(8),  4)
const VTK_QUAD             = VTKCellType("VTK_QUAD",            UInt8(9),  4)
const VTK_TETRA            = VTKCellType("VTK_TETRA",           UInt8(10), 4)
const VTK_VOXEL            = VTKCellType("VTK_VOXEL",           UInt8(11), 8)
const VTK_HEXAHEDRON       = VTKCellType("VTK_HEXAHEDRON",      UInt8(12), 8)
const VTK_WEDGE            = VTKCellType("VTK_WEDGE",           UInt8(13), 6)
const VTK_PYRAMID          = VTKCellType("VTK_PYRAMID",         UInt8(14), 5)
const VTK_PENTAGONAL_PRISM = VTKCellType("VTK_PENTAGONAL_PRISM",UInt8(15),10)
const VTK_HEXAGONAL_PRISM  = VTKCellType("VTK_HEXAGONAL_PRISM", UInt8(16),12)

# Quadratic, isoparametric cells
const VTK_QUADRATIC_EDGE                   = VTKCellType("VTK_QUADRATIC_EDGE",                  UInt8(21), 3)
const VTK_QUADRATIC_TRIANGLE               = VTKCellType("VTK_QUADRATIC_TRIANGLE",              UInt8(22), 6)
const VTK_QUADRATIC_QUAD                   = VTKCellType("VTK_QUADRATIC_QUAD",                  UInt8(23), 8)
const VTK_QUADRATIC_POLYGON                = VTKCellType("VTK_QUADRATIC_POLYGON",               UInt8(36),-1)
const VTK_QUADRATIC_TETRA                  = VTKCellType("VTK_QUADRATIC_TETRA",                 UInt8(24),10)
const VTK_QUADRATIC_HEXAHEDRON             = VTKCellType("VTK_QUADRATIC_HEXAHEDRON",            UInt8(25),20)
const VTK_QUADRATIC_WEDGE                  = VTKCellType("VTK_QUADRATIC_WEDGE",                 UInt8(26),15)
const VTK_QUADRATIC_PYRAMID                = VTKCellType("VTK_QUADRATIC_PYRAMID",               UInt8(27),13)
const VTK_BIQUADRATIC_QUAD                 = VTKCellType("VTK_BIQUADRATIC_QUAD",                UInt8(28), 9)
const VTK_TRIQUADRATIC_HEXAHEDRON          = VTKCellType("VTK_TRIQUADRATIC_HEXAHEDRON",         UInt8(29),27)
const VTK_QUADRATIC_LINEAR_QUAD            = VTKCellType("VTK_QUADRATIC_LINEAR_QUAD",           UInt8(30), 6)
const VTK_QUADRATIC_LINEAR_WEDGE           = VTKCellType("VTK_QUADRATIC_LINEAR_WEDGE",          UInt8(31),12)
const VTK_BIQUADRATIC_QUADRATIC_WEDGE      = VTKCellType("VTK_BIQUADRATIC_QUADRATIC_WEDGE",     UInt8(32),18)
const VTK_BIQUADRATIC_QUADRATIC_HEXAHEDRON = VTKCellType("VTK_BIQUADRATIC_QUADRATIC_HEXAHEDRON",UInt8(33),24)
const VTK_BIQUADRATIC_TRIANGLE             = VTKCellType("VTK_BIQUADRATIC_TRIANGLE",            UInt8(34), 7)

# Cubic, isoparametric cell
const VTK_CUBIC_LINE                       = VTKCellType("VTK_CUBIC_LINE",UInt8(35),4)

# Special class of cells formed by convex group of points
const VTK_CONVEX_POINT_SET = VTKCellType("VTK_CONVEX_POINT_SET",UInt8(41),-1)

# Polyhedron cell (consisting of polygonal faces)
const VTK_POLYHEDRON = VTKCellType("VTK_POLYHEDRON",UInt8(42),-1)

# Higher order cells in parametric form
const VTK_PARAMETRIC_CURVE        = VTKCellType("VTK_PARAMETRIC_CURVE",       UInt8(51),-1)
const VTK_PARAMETRIC_SURFACE      = VTKCellType("VTK_PARAMETRIC_SURFACE",     UInt8(52),-1)
const VTK_PARAMETRIC_TRI_SURFACE  = VTKCellType("VTK_PARAMETRIC_TRI_SURFACE", UInt8(53),-1)
const VTK_PARAMETRIC_QUAD_SURFACE = VTKCellType("VTK_PARAMETRIC_QUAD_SURFACE",UInt8(54),-1)
const VTK_PARAMETRIC_TETRA_REGION = VTKCellType("VTK_PARAMETRIC_TETRA_REGION",UInt8(55),-1)
const VTK_PARAMETRIC_HEX_REGION   = VTKCellType("VTK_PARAMETRIC_HEX_REGION",  UInt8(56),-1)

# Higher order cells
const VTK_HIGHER_ORDER_EDGE        = VTKCellType("VTK_HIGHER_ORDER_EDGE",       UInt8(60),-1)
const VTK_HIGHER_ORDER_TRIANGLE    = VTKCellType("VTK_HIGHER_ORDER_TRIANGLE",   UInt8(61),-1)
const VTK_HIGHER_ORDER_QUAD        = VTKCellType("VTK_HIGHER_ORDER_QUAD",       UInt8(62),-1)
const VTK_HIGHER_ORDER_POLYGON     = VTKCellType("VTK_HIGHER_ORDER_POLYGON",    UInt8(63),-1)
const VTK_HIGHER_ORDER_TETRAHEDRON = VTKCellType("VTK_HIGHER_ORDER_TETRAHEDRON",UInt8(64),-1)
const VTK_HIGHER_ORDER_WEDGE       = VTKCellType("VTK_HIGHER_ORDER_WEDGE",      UInt8(65),-1)
const VTK_HIGHER_ORDER_PYRAMID     = VTKCellType("VTK_HIGHER_ORDER_PYRAMID",    UInt8(66),-1)
const VTK_HIGHER_ORDER_HEXAHEDRON  = VTKCellType("VTK_HIGHER_ORDER_HEXAHEDRON", UInt8(67),-1)

end # module


text:
using Random
using RandomNumbers.Xorshifts
using LinearAlgebra

rng = Xorshift1024Plus(21457)

mutable struct BoltzmannMachine
    W::AbstractArray
    rng::AbstractRNG
end

function init!(v, h; rng = nothing)
    w = Matrix{Float64}(undef, v, h)
    # @show w
    # Glorot initialization for weights
    low = -0.1 * sqrt(6.0 / (h + v))
    rand!(rng, w, [low, -low])

    # Insert weights for the bias units into the first row and first column.
    w = vcat(zeros(1, h), w)
    w = hcat(zeros(v+1), w)

    return BoltzmannMachine(w, rng)
end

sigmoid(x) = 1 / (1 + exp(-x))

to_one(x) = x == true ? 1 : 0

function kstep(W; k = 50, rng = nothing)
    (N, M) = size(W)
    samples = rand(rng, 1, N)
    neg_associations = Matrix{Float64}(undef, N, M)

    for i = 1:k-1
        hidden_activations = samples * W
        # @show size(hidden_activations)
        hidden_probs = sigmoid.(hidden_activations)
        # @show size(hidden_probs)
        hidden_states = hidden_probs .> rand(rng, (1, M))
        hidden_states = to_one.(hidden_states)
        hidden_states[:, 1] .= 1

        visible_activations = hidden_states * W'
        visible_probs = sigmoid.(visible_activations)
        # @show size(visible_probs)
        # @show size(visible_probs')
        # @show size(hidden_probs)
        neg_associations = visible_probs' * hidden_probs
        # @show size(neg_associations)
    end

    return neg_associations
end

function train(data, rbm; epochs = 3000, η = 0.001)
    (N, M) = size(data)

    # Add bias to data
    data_bias = hcat(ones(N), data)

    neg_associations = Matrix{Float64}(undef, N, M)

    for e = 1:epochs
        pos_hidden_activations = data_bias * rbm.W
        pos_hidden_probs = sigmoid.(pos_hidden_activations)
        # Fix the bias unit
        pos_hidden_probs[:, 1] .= 1
        pos_hidden_states = pos_hidden_probs .> rand(rbm.rng, size(pos_hidden_probs)...)
        pos_hidden_states = to_one.(pos_hidden_states)
        pos_associations = data_bias' * pos_hidden_probs
        # @show size(pos_associations)

        neg_associations = kstep(rbm.W; k = 5000, rng = rbm.rng)

        # Update weights.
        @. rbm.W += η * ((pos_associations - neg_associations) / N)
    end
end

function sample(rbm; k = 50, rng = nothing)
    (N, M) = size(rbm.W)
    samples = Array{Float64}(undef, k, N, N)
    # samples = Matrix{Float64}(undef, k, N)
    samples[1, :, :] = rand(rng, 1, N, N)

    for i = 1:k-1
        hidden_activations = samples[i, :, :] * rbm.W
        hidden_probs = sigmoid.(hidden_activations)
        hidden_probs[:, 1] .= 1
        hidden_states = hidden_probs .> rand(rng, 1, M)
        hidden_states = to_one.(hidden_states)

        visible_activations = hidden_states * rbm.W'
        # @show size(visible_activations)
        visible_probs = sigmoid.(visible_activations)
        samples[i+1, :, :] = visible_probs
    end

    return samples
end

data = [
    [1, 1, 1, 0, 0, 0];
    [1, 0, 1, 0, 0, 0];
    [1, 1, 1, 0, 0, 0];
    [0, 0, 1, 1, 1, 0];
    [0, 0, 1, 1, 0, 0];
    [0, 0, 1, 1, 1, 0]
]
data = reshape(data, (6, 6))

rbm = init!(6, 2; rng = rng)

train(data, rbm)

new_data = sample(rbm; k = 2000, rng = rng)
display(new_data[end, :, :])

max_stars_repo_name: edwinb-ai/IsingRBM

max_issues_repo_name: edwinb-ai/IsingRBM

max_forks_repo_name: edwinb-ai/IsingRBM


text:
"""
# Usage

    solve(prob::SingleTermFODEProblem, h, PIEX())

### References

```tex
@inproceedings{Garrappa2018NumericalSO,
  title={Numerical Solution of Fractional Differential Equations: A Survey and a Software Tutorial},
  author={Roberto Garrappa},
  year={2018}
}
```
"""
struct PIEX <: FractionalDiffEqAlgorithm end

function solve(FODE::SingleTermFODEProblem, h, ::PIEX)
    @unpack f, α, u0, tspan = FODE
    t0 = tspan[1]; T = tspan[2]
    N::Int64 = round(Int, (T-t0)/h)
    y = zeros(N+1)

    y[1]=u0
    for j in range(2, N+1, step=1)
        middle=0
        @turbo for i=0:j-1
            middle += bcoefficients(j-i, α)*f(t0+i*h, y[i+1])
        end
        middle = middle/gamma(α+1)
        y[j] = u0 + middle*h^α
    end
    t = collect(t0:h:T)
    return FODESolution(t, y)
end

function acoefficients(n, α)
    if n == 0
        return 1
    else
        return ((n-1)^(α+1)-2*n^(α+1)+(n+1)^(α+1))
    end
end

bcoefficients(n, α) = ((n+1)^α-n^α)


text:
@testset "PrimeIdealsSet" begin
  Qx, x = FlintQQ["x"]
  K, a = NumberField(x - 1, "a")
  O = maximal_order(K)

  S = @inferred PrimeIdealsSet(O, 2, 100)
  @test @inferred length(collect(S)) == 25

  S = @inferred PrimeIdealsSet(O, 2, 100, indexdivisors = false, ramified = false, degreebound = 100, coprimeto = 3*5*7)
  @test @inferred length(collect(S)) == 22

  K, a = NumberField(x^2 - 5, "a")
  O = maximal_order(K)

  S = @inferred PrimeIdealsSet(O, fmpz(2), fmpz(100))
  @test @inferred length(collect(S)) == 35
  S = @inferred PrimeIdealsSet(O, fmpz(2), fmpz(100), indexdivisors = false)
  @test @inferred length(collect(S)) == 34
  S = @inferred PrimeIdealsSet(O, fmpz(2), fmpz(100), indexdivisors = false, ramified = false)
  @test @inferred length(collect(S)) == 33
  S = @inferred PrimeIdealsSet(O, fmpz(2), fmpz(100), indexdivisors = false, ramified = false, degreebound = 1)
  @test @inferred length(collect(S)) == 20

  K, a = NumberField(x^5 - x + 1, "a")
  O = maximal_order(K)

  S = @inferred PrimeIdealsSet(O, fmpz(2), fmpz(100))
  @test @inferred length(collect(S)) == 51

  S = @inferred PrimeIdealsSet(O, fmpz(2), fmpz(100), degreebound = 1)
  @test @inferred length(collect(S)) == 18
  SS = collect(S)

  S = @inferred PrimeIdealsSet(O, fmpz(2), -1, degreebound = 1)
  z = 1
  T = NfOrdIdl[]
  for P in S
    push!(T, P)
    if z == 18
      break
    end
    z += 1
  end
  @test T == SS

  P = prime_decomposition(O, 2)[1][1]
  S = @inferred PrimeIdealsSet(O, fmpz(2), fmpz(100), coprimeto = P)
  @test @inferred length(collect(S)) == 50

  el = Hecke.find_elem_of_valuation_1(P, P^2)
  @test valuation(el, P) == 1

  S = @inferred PrimeIdealsSet(O, fmpz(2), fmpz(100), coprimeto = 2)
  @test @inferred length(collect(S)) == 49

  S = @inferred PrimeIdealsSet(O, fmpz(2), fmpz(100), coprimeto = fmpz(6))
  @test @inferred length(collect(S)) == 48

  S = @inferred PrimeIdealsSet(O, fmpz(2), fmpz(100), coprimeto = O(30))
  @test @inferred length(collect(S)) == 47

  @test_throws ErrorException PrimeIdealsSet(O, fmpz(-1), fmpz(1))
  @test_throws ErrorException PrimeIdealsSet(O, fmpz(1), -2)
  @test_throws ErrorException PrimeIdealsSet(O, fmpz(1), 2, coprimeto = "bla")
end

@testset "PrimeIdealsSet" begin
  l = Hecke.primes_up_to(20)
  for i = 2:length(l)
    d = l[i]
    K = quadratic_field(d)[1]
    OK = maximal_order(K)
    p = next_prime(fmpz(10)^70)
    while length(prime_decomposition_type(OK, p)) != 2
      p = next_prime(p)
    end
    P = prime_decomposition(OK, p)[1][1]
    @test valuation(P.gen_two, P) == 1
  end

  _, x = PolynomialRing(QQ, cached = false)
  K = number_field(x^4 + 2*x^3 - 35*x^2 - 36*x + 5, "a", cached = false)[1]
  OK = maximal_order(K)
  @assert length(prime_decomposition_type(OK, 3)) == 2
  @assert length(prime_decomposition_type(OK, 5)) == 4
end


text:
#-------------------------------------------------------------------------------
# Type: Characters
#-------------------------------------------------------------------------------
struct Characters{T <: Number, Tg <: AbstractFiniteGroup}
    χ::AbstractVector{T}
    R::Int
    g::Tg
end

function Characters(g::AbstractFiniteGroup, χ::AbstractVector{<:Number})
    R = check_real_rep(g, χ)
    Characters(isone(R) ? real.(χ) : χ, R, g)
end
#-------------------------------------------------------------------------------
function Base.display(c::Characters)
    reality = if c.R == 1
        "Real"
    elseif c.R == 0 Complex
        "Complex"
    else
        "Pseudo-real"
    end
    println("Characters of $reality representation of $(name(c.g)):")
    println(c.χ)
end
#-------------------------------------------------------------------------------
Base.getindex(c::Characters, i) = c.χ[i]
Base.length(c::Characters) = length(c.χ)
Base.eltype(::Characters{T,Tg}) where {T,Tg} = T
Base.conj(c::Characters) = Characters(conj(c.χ), c.R, c.g)
function Base.isless(c1::Characters, c2::Characters)
    @assert isequal(c1.g, c2.g) "Compared characters should belong to tha same group."
    order_chi(c1.χ, c2.χ)
end
function ==(c1::Characters, c2::Characters)
    isequal(c1.g, c2.g) && norm(c1.χ - c2.χ) < 1e-10
end
function dot(c1::Characters, c2::Characters)
    @assert isequal(c1.g, c2.g) "Dotted characters should belong to tha same group."
    g = c1.g
    p = sum(conj(c1.χ[k]) * c2.χ[k] * mult(g, k) for k = 1:length(class(g))) / order(g)
    pint = round(Int64, p)
    @assert abs(p-pint) < 1e-10 "Inner product = $p, not integer."
    pint
end
#-------------------------------------------------------------------------------
# Type: CharacterTable
#-------------------------------------------------------------------------------
struct CharacterTable
    tab::Vector{Characters}
    names::Vector{String}
end

function CharacterTable(tab::Vector{Characters})
    names = Vector{String}(undef, length(tab))
    n = 0
    for i = 1:length(tab)
        isassigned(names, i) && continue
        n += 1
        if iszero(tab[i].R)
            χc = conj(tab[i])
            for j = i+1:length(tab)
                if tab[j] == χc
                    names[i] = "$(i)+"
                    names[j] = "$(j)-"
                    break
                end
            end
        else
            names[i] = "$(i)"
        end
    end
    CharacterTable(tab, names)
end
#-------------------------------------------------------------------------------
function Base.display(ct::CharacterTable)
    g = ct.tab[1].g
    mats = Matrix(ct)
    n = size(mats, 1)+1
    out = Matrix{Any}(undef, n, n)
    out[1, 1] = ""
    out[1, 2:end] .= [name(g, cls[1]) for cls in class(g)]
    out[2:end, 1] .= ct.names
    out[2:end, 2:end] .= beautify.(mats)
    display(out)
end
Base.Matrix(ct::CharacterTable) = vcat((transpose(c.χ) for c in ct.tab)...)
Base.length(ct::CharacterTable) = length(ct.tab)
Base.size(ct::CharacterTable) = (length(ct.tab), length(ct.tab[1]))
Base.size(ct::CharacterTable, i::Integer) = isone(i) ? length(ct.tab) : (i == 2) ? length(ct.tab[1]) : 1
Base.getindex(ct::CharacterTable, i::Integer, j) = ct.tab[i][j]
Base.getindex(ct::CharacterTable, inds::AbstractVector{<:Integer}, j) = vcat((transpose(ct.tab[i][j]) for i in inds)...)
Base.getindex(ct::CharacterTable, c::Colon, j) = vcat((transpose(ct.tab[i][j]) for i in 1:length(ct))...)
Base.getindex(ct::CharacterTable, i) = ct.tab[i]
#-------------------------------------------------------------------------------
# Main output function
#-------------------------------------------------------------------------------
export charactertable
"""
    charactertable(g::AbstractFiniteGroup; method)

Return the character table of group `g`.
Use Burnside's method.
"""
function charactertable(g::AbstractFiniteGroup; tol::Real=1e-7)
    h = class_multab(g)
    burnside(g, h, tol=tol)
end

#-------------------------------------------------------------------------------
# Burnside's Mothod
#-------------------------------------------------------------------------------
"""
    burnside(g::FiniteGroup, h::AbstractArray{<:Integer, 3}; tol::Real=1e-7)

Get character table using Burnside's method.
"""
function burnside(
    g::AbstractFiniteGroup,
    h::AbstractArray{<:Integer, 3};
    tol::Real=1e-7
)
    NC = size(h, 1)
    v = [I(NC)]
    while true
        rh = sum(rand() * h[i,:,:] for i = 1:NC)
        v = vcat([vec_split(rh, vi, tol=tol) for vi in v]...)
        all(isone(size(vi, 2)) for vi in v) && break
    end
    characters = Vector{Characters}(undef, NC)
    @threads for i = 1:NC
        χ = beautify.(reshape(normalize_chi!(g, v[i]), :), tol=tol)
        characters[i] = Characters(g, χ)
    end
    tab = sort(characters)
    CharacterTable(tab)
end
#-------------------------------------------------------------------------------
function normalize_chi!(g::AbstractFiniteGroup, χ::AbstractArray)
    NC = length(χ)
    χ ./= χ[1]
    d2 = order(g) / sum( abs2(χ[j])/mult(g,j) for j=1:NC )
    for i = 1:NC
        χ[i] *= sqrt(d2)/mult(g,i)
    end
    χ
end
#-------------------------------------------------------------------------------
function order_chi(v1::AbstractVector, v2::AbstractVector; tol::Real=1e-7)
    for i = 1:length(v1)
        a, b = log(Complex(v1[i])), log(Complex(v2[i]))
        (real(a) < -10) && (real(b) < -10) && continue
        dr = real(a) - real(b)
        if dr > tol
            return false
        elseif dr < -tol
            return true
        else
            ai = imag(a) .< -tol ? 2π + imag(a) : imag(a)
            bi = imag(b) .< -tol ? 2π + imag(b) : imag(b)
            di = ai - bi
            if di > tol
                return false
            elseif di < -tol
                return true
            end
        end
    end
    println("Warning: two vector the same.")
    false
end
#-------------------------------------------------------------------------------
"""
Split the vector space according to the split of the eigen values.
"""
function vec_split(h::AbstractMatrix, vs::AbstractVecOrMat; tol::Real=1e-7)
    isone(size(vs, 2)) && return [vs]
    php = vs' * h * vs
    e, v = eigen(php)
    e_split = spectrum_split(e, tol=tol)
    # zero-energy degeneracy is numerically unstable!!!
    isone(length(e_split)) ? [vs] : [vs * v[:, ei] for ei in e_split]
end
#-------------------------------------------------------------------------------
"""
Group the spectrum into degenerate sets.
"""
function spectrum_split(vals::AbstractVector{<:Number}; tol::Real=1e-7)
    n = length(vals)
    NG = 0
    gvec = zeros(Int, n)
    for i = 1:n
        if iszero(gvec[i])
            NG += 1
            gvec[i] = NG
        else
            continue
        end
        target = vals[i]
        # Pay attention to the sorting of complex number:
        # a - b im, a + b im, ã - b̃ im, ã + b̃ im
        #tally = 0
        for j = i+1:n
            if abs(vals[j] - target) < tol
                gvec[j] = NG
            #    tally = 0
            #else
            #    tally += 1
            end
            #(tally > 2) && break
        end
    end
    collect_group(gvec, NG)
end

#-------------------------------------------------------------------------------
# Helper Functions
#-------------------------------------------------------------------------------
export class_multab
"""
    class_multab(multable, classes, class_vec)

Compute class multiplicity:

Cᵢ × Cⱼ = hᵢⱼₖ Cₖ

Outputs:
--------
h::Array{Int64, 3} : Class multiplicity.
"""
function class_multab(g::AbstractFiniteGroup)
    NC = class(g) |> length
    h = zeros(Int64, NC, NC, NC)
    for i = 1:NC
        Ci = class(g, i)
        for a in Ci, b in Ci
            k = inclass(g, g[a, b])
            h[i, i, k] += 1
        end
        for j = i+1:NC
            Cj = class(g, j)
            for a in Ci, b in Cj
                k = inclass(g, g[a, b])
                h[i, j, k] += 1
                h[j, i, k] += 1
            end
        end
    end
    for k=1:NC
        h[:, :, k] .÷= mult(g, k)
    end
    h
end
#-------------------------------------------------------------------------------

function beautify(a::Real; tol::Real=1e-7)
    a_int = round(a)
    abs(a-a_int) < tol ? a_int : a
end
function beautify(a::Number; tol::Real=1e-7)
    if abs(imag(a)) < tol
        convert(typeof(a), beautify(real(a), tol=tol))
    elseif abs(real(a)) < tol
        convert(typeof(a), 1im * beautify(imag(a), tol=tol))
    else
        a
    end
end


text:
export ⋆ₛ#, ⋆, UniformScalingExchange, J

"""
    ⋆ₛ(A::AbstractBlockMatrix, B::AbstractBlockMatrix)

Return the star product of two S-matrices (2x2 BlockMatrices)
See:
https://en.wikipedia.org/wiki/Redheffer_star_product#Connection_to_scattering_matrices
"""
function ⋆ₛ(A::AbstractMatrix, B::AbstractMatrix)
    d = Int(size(A, 1) // 2)
    @assert size(A) == size(B) == reverse(size(A)) "maps must be endomorphisms"
    invI_A₂₂B₁₁ = inv(I - A[(d+1):2d, (d+1):2d] * B[1:d, 1:d])
    invI_B₁₁A₂₂ = inv(I - B[1:d, 1:d] * A[(d+1):2d, (d+1):2d])
    [
        A[1:d, 1:d] + A[1:d, (d+1):2d] * invI_B₁₁A₂₂ * B[1:d, 1:d] * A[(d+1):2d, 1:d]        A[1:d, (d+1):2d] * invI_B₁₁A₂₂ * B[1:d, (d+1):2d];
        B[(d+1):2d, 1:d] * invI_A₂₂B₁₁ * A[(d+1):2d, 1:d]        B[(d+1):2d, (d+1):2d] + B[(d+1):2d, 1:d] * invI_A₂₂B₁₁ * A[(d+1):2d, (d+1):2d] * B[1:d, (d+1):2d];
    ]
end

"""
    ⋆ₛ(A::AbstractBlockMatrix, B::AbstractBlockMatrix)

Return the star product of two S-matrices (2x2 BlockMatrices)
See:
https://en.wikipedia.org/wiki/Redheffer_star_product#Connection_to_scattering_matrices
"""
function ⋆ₛ(A::AbstractBlockMatrix, B::AbstractBlockMatrix)
    # The inputs at the very least need to have conformable blocks
    # However, the Redheffer star is defined for blocks of all same shape
    invI_A₂₂B₁₁ = inv(I - A[Block(2, 2)] * B[Block(1, 1)])
    invI_B₁₁A₂₂ = inv(I - B[Block(1, 1)] * A[Block(2, 2)])
    mortar(
        (A[Block(1, 1)] + A[Block(1, 2)] * invI_B₁₁A₂₂ * B[Block(1, 1)] * A[Block(2, 1)],        A[Block(1, 2)] * invI_B₁₁A₂₂ * B[Block(1, 2)]),
        (B[Block(2, 1)] * invI_A₂₂B₁₁ * A[Block(2, 1)],        B[Block(2, 2)] + B[Block(2, 1)] * invI_A₂₂B₁₁ * A[Block(2, 2)] * B[Block(1, 2)]),
    )
end

"""
    ⋆ₛ(A::BlockMap, B::BlockMap)::BlockMap

Perform the star product on BlockMaps
"""
function ⋆ₛ(A::BlockMap, B::BlockMap)::BlockMap
    @assert (2, 2) == A.rows == B.rows "Both maps must be 2x2 block maps"
    @assert all(_has_regular_blocks.([A, B])) "Both maps must have regular block alignment"
    A₁₁, A₁₂, A₂₁, A₂₂ = A.maps
    B₁₁, B₁₂, B₂₁, B₂₂ = B.maps
    invI_A₂₂B₁₁ = LinearMap(x -> gmres(I - A₂₂ * B₁₁, x), size(A₂₂, 1), size(B₁₁, 2))
    invI_B₁₁A₂₂ = LinearMap(x -> gmres(I - B₁₁ * A₂₂, x), size(B₁₁, 1), size(A₂₂, 2))
    [
        A₁₁ + A₁₂ * invI_B₁₁A₂₂ * B₁₁ * A₂₁     A₁₂ * invI_B₁₁A₂₂ * B₁₂;
        B₂₁ * invI_A₂₂B₁₁ * A₂₁                 B₂₂ + B₂₁ * invI_A₂₂B₁₁ * A₂₂ * B₁₂
    ]
end

"""
    has_regular_blocks(M::BlockMap)::Bool

Returns true if the blocks in the BlockMap are partitioned by a row partition and a
column partition (i.e. all blocks in each column/row have the same width/height).
"""
function _has_regular_blocks(M::BlockMap)::Bool
    Ncol = M.rows[1]
    Nrow = length(M.rows)
    all(Ncol == M.rows[i] for i in 1:Nrow) &&
    all(M.colranges[i] == M.colranges[i+j*Ncol] for i in 1:Ncol, j in 0:(Nrow-1))
end

"""
    ⋆ₛ(A::LinearMap, B::LinearMap)::LinearMap

Take the scattering star product of two matrix-free scattering matrices.
Requires the inputs to be square matrices of the same size and assumes that the
matrix is divided into equal-size quadrants to obtain the blocks
"""
function ⋆ₛ(A::LinearMap, B::LinearMap)::LinearMap
    d = Int(size(A, 1) // 2)
    @assert size(A) == size(B) == reverse(size(A)) "maps must be endomorphisms"
    A₁₁ = LinearMap(x -> (A * vcat(x, zero(x)))[1:d], d)
    A₁₂ = LinearMap(x -> (A * vcat(zero(x), x))[1:d], d)
    A₂₂ = LinearMap(x -> (A * vcat(zero(x), x))[(1+d):2d], d)
    B₁₁ = LinearMap(x -> (B * vcat(x, zero(x)))[1:d], d)
    B₂₁ = LinearMap(x -> (B * vcat(x, zero(x)))[(1+d):2d], d)
    B₂₂ = LinearMap(x -> (B * vcat(zero(x), x))[(1+d):2d], d)
    LinearMap(2d) do x
        x₁ = @view x[1:d]
        x₂ = @view x[(d+1):2d]
        # This is slightly optimized to avoid repeating two linear mappings
        y = A * vcat(x₁, fill(zero(eltype(x)), d))
        y₁ = @view y[1:d]
        y₂ = gmres(I - A₂₂ * B₁₁, y[(d+1):2d])
        z = B * vcat(fill(zero(eltype(x)), d), x₂)
        z₁ = gmres(I - B₁₁ * A₂₂, z[1:d])
        z₂ = @view z[(d+1):2d]
        vcat(
            y₁ + A₁₂ * B₁₁ * y₂ + A₁₂ * z₁,
            B₂₁ * y₂ + z₂ + B₂₁ * A₂₂ * z₁
        )
    end
end

# """
#     exchangehalves(A::AbstractArray{<:Any,N}, dim::Integer) where N

# Exchange the first half of A along first dimension with the second half.
# Will fail if size(A, dim) is odd.
# """
# function exchangehalves(A::AbstractArray{<:Any, N}, dim::Integer) where N
#     if N == 0
#         return A
#     else
#         @assert 1 <= dim <= N
#         return _exchangehalves(A, dim)
#     end
# end
# function _exchangehalves(A::AbstractArray, dim)
#     d = Int(size(A, dim) // 2)
#     selectdim(A, dim, vcat((d+1):2d, 1:d))
# end
# # function _exchangehalves(A::BlockArray, dim)
# #     @assert length(blockaxes(A, dim)) == 2
# #     selectdim(A, dim, blockslice)
# # end

# """
#     ExchangeOperator

# Type equivalent to the permutation matrix

#     [
#         0 1;
#         1 0
#     ]

# with flexible size set so that each block is the same square size.
# """
# struct ExchangeOperator end

# """
#     J

# Object of type ExchangeOperator
# """
# const J = ExchangeOperator()

# """
#     ExchangeOperator()(d::Int)

# Returns a LinearMap version of the exchangehalves operator. `d` is the axis length
# and must be even.
# """
# (::ExchangeOperator, d::Int) = LinearMap{Bool}(d; issymmetric=true, ishermitian=true) do A
#     @assert size(A, 1) == d
#     exchangehalves(A, 1)
# end

# Base.:*(::ExchangeOperator, A::AbstractArray) = exchangehalves(A, 1)
# Base.:*(A::AbstractArray{<:Any,N}, ::ExchangeOperator) where N = exchangehalves(A, N)

# Base.:*(::ExchangeOperator, A::LinearMap) = J(size(A, 1)) * A
# Base.:*(A::LinearMap, ::ExchangeOperator) = A * J(size(A, 2))

# """
#     ⋆(A, B)

# Return the star product of two matrices
# See:
# https://en.wikipedia.org/wiki/Redheffer_star_product
# """
# ⋆(A, B) = J * ((J * A) ⋆ₛ (J * B))


text:
# workaround to return all variables, so that nothing get's optimized away,
# but don't use an array, since it will sprinkle the code with any's and grow_array

const syms = []

macro save(x)
    s = esc(gensym())
    push!(syms, s)
    :($s = $x)
end
macro ret()
    :(tuple($(syms...)))
end
immutable TestFunctor end
@compat (::TestFunctor)(a, b) = a+b

function use_operations()

    a = Vec{3, Float32}(0)
    b = Vec{3, Float32}(0.0)
    c = Vec{3, Float32}(0f0)

    a2 = Point{3}(0)
    b2 = Point{3}(0.0)
    c2 = Point{3}(0f0)

    a3 = Point(0f0, 1, 1.)
    b3 = Point{3, Float32}(0f0, 1, 1.)
    c3 = Point{3, Float32}((0f0, 1, 1.))

    m1 = @fsa([1 2;3 4;5 6])
    m2 = eye(Mat{3,3,Float32})

    @save dot(a, a)
    @save dot(a2, a2)

    @save a./b
    @save a.*b
    @save a == b
    @save m2*c
    @save m1.+m1
    @save m1./m1
    @save m2.*m2

    @save reduce(TestFunctor(), a3)
    @save reduce(TestFunctor(), m2)

    @save map(TestFunctor(), a3, a3)
    @save map(TestFunctor(), m2, m2)

    @save map(Float32, a3)
    @save map(Float64, m2)

    @save map(Float64, a3)
    @save map(Float32, m2)

    @ret
end

# This seems to be the easiest way to test,
# that all variables even if they come from inlining,
# are inferred as a concrete type.
# seems like it's hard to turn of inlining when doing coverage, though...
# so this must be executed locally to test for it.
context("type inference") do
    io = IOBuffer()
    use_operations()
    code_warntype(io, use_operations, ())
    str = takebuf_string(io)
    x = matchall(r"Any", str)
    @fact length(x) --> 0
end


text:
module SuperType

using BinaryTraits, Test

struct Duck end
abstract type Mobility end

@trait Fly as Mobility

function test()
    @test supertype(Fly) === Mobility
end

end # module

using Test
@testset "Super Type" begin
    import .SuperType
    SuperType.test()
end


text:
using DecisionTree

## Function that list all the leafs in a tree and the number of partition in each leaf
## It returns a dictionary.

function list_leafs!(leafs::Dict{Int64, Int64}, leaf::Leaf{T}, leafid::Int64, index::Int64) where {T}
    leafs[leafid] = index
    return index + 1
end

function list_leafs!(leafs::Dict{Int64, Int64}, node::Node{S, T}, leafid::Int64, index::Int64) where {S, T}
    index = list_leafs!(leafs, node.left, leafid << 1, index)
    index = list_leafs!(leafs, node.right, (leafid << 1) + 1, index)
    return index
end

function list_leafs(tree::DecisionTree.LeafOrNode{S, T}) where {S, T}
    leafs = Dict{Int64, Int64}()
    list_leafs!(leafs, tree, 1, 1)
    return leafs
end

## Function that returns a dictionnary
## Keys are the leafs names
## Values are the vectors of id of the features in the leaf

get_leaf_of_feature(leaf::Leaf{T}, feature::AbstractVector{S}, leafid::Int64) where {S, T} = leafid

function get_leaf_of_feature(tree::Node{S, T}, feature::AbstractVector{S}, leafid::Int64) where {S, T}
    if feature[tree.featid] < tree.featval
        return get_leaf_of_feature(tree.left, feature, leafid << 1)
    else
        return get_leaf_of_feature(tree.right, feature, (leafid << 1) + 1)
    end
end

function get_leaf_of_feature(tree::DecisionTree.LeafOrNode{S, T}, feature::AbstractVector{S}) where {S, T}
    leafnames = list_leafs(tree)
    return leafnames[get_leaf_of_feature(tree, feature, 1)]
end

function get_leaf_of_feature(tree::DecisionTree.LeafOrNode{S, T}, features::AbstractMatrix{S}) where {S, T}
    leafnames = list_leafs(tree)

    N = size(features, 1)
    leafs = Array{Int64, 1}(undef, N)
    for i in 1:N
        leafs[i] = leafnames[get_leaf_of_feature(tree, view(features, i, :), 1)]
    end

    return leafs
end

function get_repartition(tree::DecisionTree.LeafOrNode{S, T}, features::AbstractMatrix{S}, bootstrap::AbstractVector{Int64}) where {S, T}
    nleafs = length(tree)

    leafs_of_features = get_leaf_of_feature(tree, features[bootstrap, :])

    count_partition = zeros(Int64, nleafs)
    for idleaf in leafs_of_features
        count_partition[idleaf] += 1
    end

    partition = Array{Array{Int64, 1}, 1}(undef, nleafs)
    for ileaf in 1:nleafs
        partition[ileaf] = Array{Int64, 1}(undef, count_partition[ileaf])
        count_partition[ileaf] = 0
    end

    for i in 1:length(leafs_of_features)
        ileaf = leafs_of_features[i]
        count_partition[ileaf] += 1
        partition[ileaf][count_partition[ileaf]] = bootstrap[i]
    end

    return partition
end


text:
using Wflow
using Dates
using TOML

# edit existing TOML to get a short run in a temporary directory
tomlpath = joinpath(@__DIR__, "sbm_config.toml")
confdict = TOML.parsefile(tomlpath)
confdict["endtime"] = DateTime("2000-01-03T00:00:00")
# convert relative paths to absolute paths, to avoid having to copy input files
confdict["input"]["path_forcing"] = joinpath(@__DIR__, confdict["input"]["path_forcing"])
confdict["input"]["path_static"] = joinpath(@__DIR__, confdict["input"]["path_static"])
tmp_dir = mktempdir()
tomlpath = joinpath(tmp_dir, "config.toml")
open(tomlpath, "w") do io
    TOML.print(io, confdict)
end

Wflow.run_simulation(tomlpath)

#=
# test whether restarted runs get the same results as continuus ones, i.e. state is captured
tomlpath = joinpath(@__DIR__, "sbm_config.toml")
config = Wflow.Config(tomlpath)

# January at once, cold start
config.starttime = DateTime("2000-01-01T00:00:00")
config.endtime   = DateTime("2000-02-01T00:00:00")
config.model.reinit = true  # cold start
# note that this needs to be relative to the tomlpath
config.state.path_output = joinpath(dirname(tomlpath), "data/state-test/outstates-moselle-january.nc")
config.output.path = joinpath(dirname(tomlpath), "data/state-test/output-moselle-january.nc")
model = Wflow.initialize_sbm_model(config)
Wflow.run_simulation(model)

# first half of January, cold start
config.starttime = DateTime("2000-01-01T00:00:00")
config.endtime   = DateTime("2000-01-15T00:00:00")
config.model.reinit = true  # cold start
config.state.path_output = joinpath(dirname(tomlpath), "data/state-test/outstates-moselle-january-1of2.nc")
config.output.path = joinpath(dirname(tomlpath), "data/state-test/output-moselle-january-1of2.nc")
model = Wflow.initialize_sbm_model(config)
Wflow.run_simulation(model)

# second half of January, warm start
config.starttime = DateTime("2000-01-15T00:00:00")
config.endtime   = DateTime("2000-02-01T00:00:00")
config.model.reinit = false  # cold start
config.state.path_input = joinpath(dirname(tomlpath), "data/state-test/outstates-moselle-january-1of2.nc")
config.state.path_output = joinpath(dirname(tomlpath), "data/state-test/outstates-moselle-january-2of2.nc")
config.output.path = joinpath(dirname(tomlpath), "data/state-test/output-moselle-january-2of2.nc")
model = Wflow.initialize_sbm_model(config)
Wflow.run_simulation(model)
=#


text:
"""
    ComputGraphStyle <: Broadcast.BroadcastStyle

This style of broadcast will forward the broadcast expression
to be registered in a computation graph, rather than directly
calculate it.
"""
struct ComputGraphStyle <: Broadcast.BroadcastStyle end
Base.BroadcastStyle(::Type{<:Value}) = ComputGraphStyle()
Broadcast.BroadcastStyle(s::ComputGraphStyle, x::Broadcast.BroadcastStyle) = s

# # this enables method traits broadcast as a constant
Broadcast.broadcastable(x::Value) = x

function Broadcast.broadcasted(::ComputGraphStyle, f, args...)
    mt = Trait.Broadcasted(f)
    register(mt, args...)
end

Broadcast.materialize(x::Value) = register(Broadcast.materialize, x)

# NOTE: we don't need it anymore since we implemented backward for materialize
# gradient(::typeof(Broadcast.materialize), grad, output, x) = (grad, )

# directly forward to broadcasted
function backward(node::CachedNode, ::typeof(Broadcast.materialize), grad)
    backward_type_assert(node, grad)
    # TODO: replace with @assert when there is a compiler option for it
    @boundscheck backward_size_assert(node, grad)
    backward(arg(node, 1), grad)
end

#
function backward(node::CachedNode, ::Trait.Broadcasted, grad)
    grad_inputs = gradient(node, grad)
    for (each, each_grad) in zip(args(node), grad_inputs)
        backward(each, each_grad)
    end
    nothing
end

# arraymath.jl
for sym in (:(/), :(\), :*, :+, :-)
    f = Expr(:., :Base, QuoteNode(sym))

    if f != :/
        @eval ($f)(A::Number, B::Variable{<:AbstractArray}) = broadcast($f, A, B)
        @eval ($f)(A::Number, B::CachedNode{<:Node, <:AbstractArray}) = broadcast($f, A, B)
    end
    if f != :\
        @eval ($f)(A::Variable{<:AbstractArray}, B::Number) = broadcast($f, A, B)
        @eval ($f)(A::CachedNode{<:Node, <:AbstractArray}, B::Number) = broadcast($f, A, B)
    end
end

for sym in (:-, :conj, :real, :imag)
    f = Expr(:., :Base, QuoteNode(sym))
    @eval ($f)(A::Variable{<:AbstractArray}) = broadcast($f, A)
    @eval ($f)(A::CachedNode{<:Node, <:AbstractArray}) = broadcast($f, A)
end


text:
const ev_to_K = 11604.52500617u"K"

@testset "Frequencies" begin
    n_e = 1.0u"1/cm^3"
    ω_pe = electron_plasma_frequency(n_e)
    ω_pe_approx = 5.64e4 * sqrt(n_e |> ustrip) * u"s^-1"
    @test isapprox(ω_pe, ω_pe_approx, atol=1e2u"s^-1")

    n_i = 1.0u"1/cm^3"
    Z = 1.0
    m_i = 1.0u"g"
    μ = (m_i / ProtonMass) |> NoUnits
    ω_pi = ion_plasma_frequency(n_i, m_i, Z)
    ω_pi_approx = 1.32e3 * Z * μ^(-1/2) * sqrt(n_i |> ustrip) * u"s^-1"
    @test isapprox(ω_pi, ω_pi_approx, atol=1u"s^-1")
end

@testset "Velocities" begin
    T_eV = 1u"eV"
    TK = ustrip(u"eV", T_eV) * ev_to_K
    @test isapprox(electron_thermal_velocity(TK), 4.19e7u"cm/s", atol=1e5u"cm/s")
    @test isapprox(electron_thermal_velocity(T_eV), 4.19e7u"cm/s", atol=1e5u"cm/s")
end


text:
validate(square) = 1 <= square <= 64 || throw(DomainError("Invalid value: $square"))

"""Calculate the number of grains on square `square`."""
function on_square(square)
    validate(square)
    BigInt(1) * BigInt(2)^(square - 1)
    # exp(log(BigInt(1)) + log(2) * (square - 1))
end

"""Calculate the total number of grains after square `square`."""
function total_after(square)
    validate(square)
    sum(on_square(i) for i in 1:square)
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning