text:
module wrapper
    ################################################################
    # Original generator script found at:
    # https://github.com/JuliaInterop/Clang.jl#build-a-custom-c-bindings-generator
    ################################################################
    using Clang

    ################################################################
    # Wraper generator
    # Loop on all headerfiles files and generate 2 julia files:
    #   - API_FILE: functions wrapper
    #   - COMMON_FILE: data wrapper
    ################################################################
    function generate(headers, includes, libname, apifile, commonfile, ignore)
        # create a work context
        ctx = DefaultContext()

        # parse headers
        parse_headers!(ctx, headers, includes=includes)

        # settings
        ctx.libname = libname
        ctx.options["is_function_strictly_typed"] = false
        ctx.options["is_struct_mutable"] = false

        # write output
        api_dir = joinpath(@__DIR__, "..", "src", "bindings")
        if !isdir(api_dir) mkdir(api_dir) end
        api_file = joinpath(api_dir, apifile)
        api_stream = open(api_file, "w")

        let api_buffer = Array{Expr}(undef,0)

            for trans_unit in ctx.trans_units
                root_cursor = getcursor(trans_unit)
                push!(ctx.cursor_stack, root_cursor)
                header = spelling(root_cursor)
                @info "+ Wrapping header: $header ..."
                # loop over all of the child cursors and wrap them, if appropriate.
                ctx.children = children(root_cursor)
                for (i, child) in enumerate(ctx.children)
                    child_name = name(child)
                    child_header = filename(child)
                    ctx.children_index = i
                    # choose which cursor to wrap
                    startswith(child_name, "__") && continue  # skip compiler definitions
                    child_name in keys(ctx.common_buffer) && continue  # already wrapped
                    child_name in keys(ctx.api_buffer) && continue  # already wrapped
                    child_header != header && continue  # skip if cursor filename is not in the headers to be wrapped
                    isinlined(child) && continue # skip inlined functions (static inline procedures are not in library ABI)
                    child_name in ignore && continue # explicitly ignore structs/functions
                    @info "  | + $child_name"
                    wrap!(ctx, child)
                end
                # Remove duplicated Expr already present in a different header
                @info "- Removing duplicates from: $header ..."
                for (i,element) in enumerate(ctx.api_buffer)
                    if element in api_buffer
                        name = string(element.args[1])
                        @info "  | - $name"
                        deleteat!(ctx.api_buffer, i)
                    end
                end
                api_buffer = vcat(api_buffer,ctx.api_buffer)
                # Write down context api_buffer
                @info "Writing $(api_file)"
                println(api_stream, "# Julia wrapper for header: $(basename(header))")
                println(api_stream, "# Automatically generated using Clang.jl\n")
                print_buffer(api_stream, ctx.api_buffer)
                empty!(ctx.api_buffer)  # clean up api_buffer for the next header
            end
            empty!(api_buffer)
            close(api_stream)
        end

        # write "common" definitions: types, typealiases, etc.
        common_file = joinpath(api_dir, commonfile)
        @info "Writing $(common_file)"
        open(common_file, "w") do f
            println(f, "# Automatically generated using Clang.jl\n")
            print_buffer(f, dump_to_buffer(ctx.common_buffer))
        end

        ################################################################
        # Generate default main file at ../src
        ################################################################
        # uncomment the following code to generate dependency and template files
        # copydeps(dirname(api_file))
        # print_template(joinpath(dirname(api_file), "LibTemplate.jl"))

    end

end

################################################################
# Load sc/p4est header filenames
################################################################
headers_file = "headers.jl"
if !isfile(headers_file)
    error("Headers file (\"headers.jl\") was not found.")
end
include(headers_file)

################################################################
# User define p4est and MPI include paths
################################################################
P4EST_ROOT_DIR  = haskey(ENV,"P4EST_ROOT_DIR") ? ENV["P4EST_ROOT_DIR"] : "/usr"
MPI_INCLUDE_DIR = haskey(ENV,"MPI_INCLUDE_DIR") ? ENV["MPI_INCLUDE_DIR"] : "/usr/include/x86_64-linux-gnu/openmpi"

################################################################
# Build absolute paths
################################################################
const P4EST_INCLUDE     = joinpath(dirname(P4EST_ROOT_DIR), "include") |> normpath
const SC_HEADERS        = [joinpath(P4EST_INCLUDE, header) for header in SC_HEADER_FILES if endswith(header, ".h")]
const P4EST_HEADERS     = [joinpath(P4EST_INCLUDE, header) for header in P4EST_HEADER_FILES if endswith(header, ".h")]
const P6EST_HEADERS     = [joinpath(P4EST_INCLUDE, header) for header in P6EST_HEADER_FILES if endswith(header, ".h")]
const P8EST_HEADERS     = [joinpath(P4EST_INCLUDE, header) for header in P8EST_HEADER_FILES if endswith(header, ".h")]
const MPI_INCLUDES      = [MPI_INCLUDE_DIR]

################################################################
# Explicitly ignore some structs/functions during automatic generation
################################################################
const SC_IGNORE         = ["sc_array"]
const PXEST_IGNORE      = ["sc_extern_c_hack_3()", "sc_extern_c_hack_4()"]

################################################################
# Generate API and COMMON files
################################################################
includes = vcat(MPI_INCLUDES, P4EST_INCLUDE)
wrapper.generate(SC_HEADERS, includes, "p4est_lib", "sc_api.jl", "sc_common.jl", SC_IGNORE)
wrapper.generate(P4EST_HEADERS, includes, "p4est_lib", "p4est_api.jl", "p4est_common.jl", PXEST_IGNORE)
wrapper.generate(P6EST_HEADERS, includes, "p4est_lib", "p6est_api.jl", "p6est_common.jl", PXEST_IGNORE)
wrapper.generate(P8EST_HEADERS, includes, "p4est_lib", "p8est_api.jl", "p8est_common.jl", PXEST_IGNORE)


text:
function  mkv2t(t,ne)

sz = size(t);
if sz[1]==ne
    t = t';
end

ndof = maximum(t[:]);
re = zeros(Int,ndof,1); # store number of neighboring elements for each entity
for i = 1:ne
    k = t[:,i];
    k = k[k.>0];
    re[k] = re[k] .+ 1;
end
re = [0; cumsum(re,dims=1)];

ce = zeros(Int,re[end],1);  # store neighboring-element indices for each entity
in = ones(Int,ndof,1);
for i = 1:ne
    k = t[:,i];   # entities on element i
    k = k[k.>0];
    # re(k): pointer to the list of entities k
    ce[re[k]+in[k]] .= i;
    in[k] = in[k] .+ 1;
end

return re,ce

end

ext: jl

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
function mass_matrix_converter_entries!(
    mass_matrix,
    converter::C,
    global_index::Base.ImmutableDict{Symbol, Int64},
) where {C <: PSY.Converter}
    @debug "Using default mass matrix entries $C"
end

function mdl_converter_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    dynamic_device::DynamicWrapper{
        PSY.DynamicInverter{PSY.AverageConverter, O, IC, DC, P, F},
    },
) where {
    O <: PSY.OuterControl,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    P <: PSY.FrequencyEstimator,
    F <: PSY.Filter,
}

    #Obtain inner variables for component
    md = inner_vars[md_var]
    mq = inner_vars[mq_var]
    Vdc = inner_vars[Vdc_var]
    θ_oc = inner_vars[θ_oc_var]

    #Transform reference frame to grid reference frame
    m_ri = dq_ri(θ_oc + pi / 2) * [md; mq]

    #Update inner_vars
    inner_vars[Vr_cnv_var] = m_ri[R] * Vdc
    inner_vars[Vi_cnv_var] = m_ri[I] * Vdc
end

function mdl_converter_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    dynamic_device::DynamicWrapper{
        PSY.DynamicInverter{PSY.RenewableEnergyConverterTypeA, O, IC, DC, P, PSY.RLFilter},
    },
) where {
    O <: PSY.OuterControl,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    P <: PSY.FrequencyEstimator,
}
    #Obtain inner variables for component
    V_R = inner_vars[Vr_inv_var]
    V_I = inner_vars[Vi_inv_var]
    V_t = sqrt(V_R^2 + V_I^2)
    θ = atan(V_I / V_R)
    Ip_cmd = inner_vars[Id_ic_var]
    Iq_cmd = inner_vars[Iq_ic_var]

    #Get Converter parameters
    converter = PSY.get_converter(dynamic_device)
    T_g = PSY.get_T_g(converter)
    Rrpwr = PSY.get_Rrpwr(converter)
    Brkpt = PSY.get_Brkpt(converter)
    Zerox = PSY.get_Zerox(converter)
    Lvpl1 = PSY.get_Lvpl1(converter)
    Vo_lim = PSY.get_Vo_lim(converter)
    Lv_pnt0, Lv_pnt1 = PSY.get_Lv_pnts(converter)
    Io_lim = PSY.get_Io_lim(converter)
    T_fltr = PSY.get_T_fltr(converter)
    K_hv = PSY.get_K_hv(converter)
    Iqr_min, Iqr_max = PSY.get_Iqr_lims(converter)
    Lvpl_sw = PSY.get_Lvpl_sw(converter)
    Q_ref = PSY.get_Q_ref(converter)
    R_source = PSY.get_R_source(converter)
    X_source = PSY.get_X_source(converter)

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(dynamic_device, PSY.RenewableEnergyConverterTypeA)
    #Define internal states for Converter
    internal_states = @view device_states[local_ix]
    Ip = internal_states[1]
    Iq = internal_states[2]
    Vmeas = internal_states[3]

    # Compute additional variables
    # Active Power Part
    #Update Ip ramp limits
    Rp_up = Ip >= 0 ? Rrpwr : Inf
    Rp_dn = Ip >= 0 ? -Inf : -Rrpwr
    #Saturate Ip if LVPL is active
    Ip_sat = Ip
    Ip_binary = 1.0
    if Lvpl_sw == 1
        LVPL = get_LVPL_gain(Vmeas, Zerox, Brkpt, Lvpl1)
        Ip_sat = Ip <= LVPL ? Ip : LVPL
        Ip_binary = Ip <= LVPL ? 1.0 : 0.0
    end
    Ip_in = clamp(Ip_cmd - Ip_sat, Rp_dn, Rp_up)
    #Get Low Voltage Active Current Management Gain
    G_lv = get_LV_current_gain(V_t, Lv_pnt0, Lv_pnt1)

    # Reactive Power Part
    #Update Iq ramp limits
    Rq_up = Q_ref >= 0 ? Iqr_max : Inf
    Rq_dn = Q_ref >= 0 ? -Inf : Iqr_min
    Iq_in = clamp(Iq_cmd - Iq, Rq_dn, Rq_up)
    Iq_extra = max(K_hv * (V_t - Vo_lim), 0.0)
    Id_cnv = G_lv * Ip_sat
    Iq_cnv = -Iq - Iq_extra
    #Reference Transformation
    Ir_cnv = Id_cnv * cos(θ) - Iq_cnv * sin(θ)
    Ii_cnv = Id_cnv * sin(θ) + Iq_cnv * cos(θ)

    ###
    #Obtain parameters
    filt = PSY.get_filter(dynamic_device)
    rf = PSY.get_rf(filt)
    lf = PSY.get_lf(filt)

    function V_cnv_calc(Ir_cnv, Ii_cnv, Vr_inv, Vi_inv)
        if lf != 0.0 || rf != 0.0
            Z_source_mag_sq = R_source^2 + X_source^2
            Zf_mag_sq = rf^2 + lf^2
            r_total_ratio = rf / Zf_mag_sq + R_source / Z_source_mag_sq
            l_total_ratio = lf / Zf_mag_sq + X_source / Z_source_mag_sq
            denom = 1.0 / (r_total_ratio^2 + l_total_ratio^2)
            rf_ratio = rf / Zf_mag_sq
            lf_ratio = lf / Zf_mag_sq
            Vr_cnv =
                denom * (
                    (Ir_cnv + Vi_inv * lf_ratio + Vr_inv * rf_ratio) * r_total_ratio -
                    (Ii_cnv + Vi_inv * rf_ratio - Vr_inv * lf_ratio) * l_total_ratio
                )
            Vi_cnv =
                denom * (
                    (Ir_cnv + Vi_inv * lf_ratio + Vr_inv * rf_ratio) * l_total_ratio +
                    (Ii_cnv + Vi_inv * rf_ratio - Vr_inv * lf_ratio) * r_total_ratio
                )
        else
            Vr_cnv = Vr_inv
            Vi_cnv = Vi_inv
        end
        return Vr_cnv, Vi_cnv
    end

    #Compute converter voltage
    Vr_cnv, Vi_cnv = V_cnv_calc(Ir_cnv, Ii_cnv, V_R, V_I)

    #Update ODEs
    output_ode[local_ix[1]] = Ip_binary * (1.0 / T_g) * Ip_in #(Ip_cmd - Ip)
    output_ode[local_ix[2]] = (1.0 / T_g) * Iq_in # (Iq_cmd - Iq)
    output_ode[local_ix[3]] = (1.0 / T_fltr) * (V_t - Vmeas)

    #Update inner_vars
    inner_vars[Ir_cnv_var] = Ir_cnv
    inner_vars[Ii_cnv_var] = Ii_cnv
    inner_vars[Vr_cnv_var] = Vr_cnv
    inner_vars[Vi_cnv_var] = Vi_cnv
end


text:
function FieldGammaMnsFunction(arg0::jint, arg1::RealFieldElement, arg2::jint, arg3::Field)
    return FieldGammaMnsFunction((jint, RealFieldElement, jint, Field), arg0, arg1, arg2, arg3)
end

function get_derivative(obj::FieldGammaMnsFunction, arg0::jint, arg1::jint, arg2::jint)
    return jcall(obj, "getDerivative", RealFieldElement, (jint, jint, jint), arg0, arg1, arg2)
end

function get_value(obj::FieldGammaMnsFunction, arg0::jint, arg1::jint, arg2::jint)
    return jcall(obj, "getValue", RealFieldElement, (jint, jint, jint), arg0, arg1, arg2)
end


text:
using BlockArrays
using BenchmarkTools
using FileIO

include("generate_report.jl")

const SUITE = BenchmarkGroup()

g = addgroup!(SUITE, "indexing")
# g_block = addgroup!(SUITE, "blockindexing")
g_size = addgroup!(SUITE, "size")

for n = (5,)
    for BT in (BlockArray, PseudoBlockArray)
        block_vec = BT(rand(n),       [1,3,1])
        block_mat = BT(rand(n,n),     [1,3,1], [4,1])
        block_arr = BT(rand(n,n,n),   [1,3,1], [4,1], [3, 2])
        g["getindex", BT.name.name, "vector", n] = @benchmarkable getindex($block_vec, 3)
        g["getindex", BT.name.name, "matrix", n] = @benchmarkable getindex($block_mat, 3, 2)
        g["getindex", BT.name.name, "rank3", n]  = @benchmarkable getindex($block_arr, 3, 2 ,3)

        g["setindex!", BT.name.name, "vector", n] = @benchmarkable setindex!($block_vec, 3)
        g["setindex!", BT.name.name, "matrix", n] = @benchmarkable setindex!($block_mat, 3, 2)
        g["setindex!", BT.name.name, "rank3", n]  = @benchmarkable setindex!($block_arr, 3, 2 ,3)

        g_size[BT.name.name, "vector", n] = @benchmarkable size($block_vec)
        g_size[BT.name.name, "matrix", n] = @benchmarkable size($block_mat)
        g_size[BT.name.name, "rank3", n]  = @benchmarkable size($block_arr)
    end
end

function run_benchmarks(name, tagfilter = @tagged ALL)
    const paramspath = joinpath(dirname(@__FILE__), "params.jld")
    if !isfile(paramspath)
        println("Tuning benchmarks...")
        tune!(SUITE, verbose=true)
        JLD.save(paramspath, "SUITE", params(SUITE))
    end
    loadparams!(SUITE, JLD.load(paramspath, "SUITE"), :evals, :samples)
    results = run(SUITE[tagfilter], verbose = true, seconds = 2)
    JLD.save(joinpath(dirname(@__FILE__), name * ".jld"), "results", results)
end

function generate_report(v1, v2)
    v1_res = load(joinpath(dirname(@__FILE__), v1 * ".jld"), "results")
    v2_res = load(joinpath(dirname(@__FILE__), v2 * ".jld"), "results")
    open(joinpath(dirname(@__FILE__), "results_$(v1)_$(v2).md"), "w") do f
        printreport(f, judge(minimum(v1_res), minimum(v2_res)); iscomparisonjob = true)
    end
end

function generate_report(v1)
    v1_res = load(joinpath(dirname(@__FILE__), v1 * ".jld"), "results")
    open(joinpath(dirname(@__FILE__), "results_$(v1).md"), "w") do f
        printreport(f, minimum(v1_res); iscomparisonjob = false)
    end
end


text:
# Function to load data from CSV file. Empty lines indicate end of a sequence.
function load(path::ASCIIString)
    stream = open(path, "r")
    X = Array{XT,1}[ XT[ ] ]
    Y = Array{YT,1}[ YT[ ] ]
    for line in map(strip, eachline(stream))
        if isempty(line)
            push!(X, XT[ ])
            push!(Y, YT[ ])
        else
            x1, x2, y = split(line, ',')
            push!(X[end], [ float(x1), float(x2) ])
            push!(Y[end], y)
        end
    end
    close(stream)
    return X, Y
end

function colorize(col::Symbol, str::ASCIIString)
    get(Base.text_colors, col, :normal) * str * get(Base.text_colors, :normal, :normal)
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# This method is called when a SemidiscretizationHyperbolic is constructed.
# It constructs the basic `cache` used throughout the simulation to compute
# the RHS etc.
function create_cache(mesh::StructuredMesh, equations::AbstractEquations, dg::DG, ::Any, ::Type{uEltype}) where {uEltype<:Real}
  elements = init_elements(mesh, equations, dg.basis, uEltype)

  cache = (; elements)

  # Add specialized parts of the cache required to compute the volume integral etc.
  cache = (;cache..., create_cache(mesh, equations, dg.volume_integral, dg, uEltype)...)

  return cache
end

# Extract contravariant vector Ja^i (i = index) as SVector
@inline function get_contravariant_vector(index, contravariant_vectors, indices...)
  SVector(ntuple(@inline(dim -> contravariant_vectors[dim, index, indices...]), Val(ndims(contravariant_vectors) - 3)))
end

@inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation,
                                                  boundary_condition::BoundaryConditionPeriodic,
                                                  mesh::StructuredMesh, equations,
                                                  surface_integral, dg::DG, cache,
                                                  direction, node_indices, surface_node_indices, element)
  @assert isperiodic(mesh, orientation)
end

@inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation,
                                                  boundary_condition,
                                                  mesh::StructuredMesh, equations,
                                                  surface_integral, dg::DG, cache,
                                                  direction, node_indices, surface_node_indices, element)
  @unpack node_coordinates, contravariant_vectors, inverse_jacobian = cache.elements
  @unpack surface_flux = surface_integral

  u_inner = get_node_vars(u, equations, dg, node_indices..., element)
  x = get_node_coords(node_coordinates, equations, dg, node_indices..., element)

  # If the mapping is orientation-reversing, the contravariant vectors' orientation
  # is reversed as well. The normal vector must be oriented in the direction
  # from `left_element` to `right_element`, or the numerical flux will be computed
  # incorrectly (downwind direction).
  sign_jacobian = sign(inverse_jacobian[node_indices..., element])

  # Contravariant vector Ja^i is the normal vector
  normal = sign_jacobian * get_contravariant_vector(orientation, contravariant_vectors,
                                                    node_indices..., element)

  # If the mapping is orientation-reversing, the normal vector will be reversed (see above).
  # However, the flux now has the wrong sign, since we need the physical flux in normal direction.
  flux = sign_jacobian * boundary_condition(u_inner, normal, direction, x, t, surface_flux, equations)

  for v in eachvariable(equations)
    surface_flux_values[v, surface_node_indices..., direction, element] = flux[v]
  end
end

include("containers.jl")
include("dg_1d.jl")
include("dg_2d.jl")
include("dg_3d.jl")

include("indicators_1d.jl")
include("indicators_2d.jl")
include("indicators_3d.jl")

# Specialized implementations used to improve performance
include("dg_2d_compressible_euler.jl")
include("dg_3d_compressible_euler.jl")

end # @muladd


text:
function import_panel(name::AbstractString, bedfile::AbstractString)
    open(rel".." * "/data/panel/$name.panel", "w") do f
        chr, buf = "", IOBuffer()

        for line in eachline(bedfile) @when length(line) > 1
            c, a, b = split(line)
            if c != chr
                if chr != ""
                    data = takebuf_array(buf)
                    write(f, chr, '\0', length(data)>>3, data)
                end
                chr = c
            end

            a, b = parse(i32, a)+i32(1), parse(i32, b)
            write(buf, a, b)
        end

        data = takebuf_array(buf)
        write(f, chr, '\0', length(data)>>3, data)
    end
end

function load_panel(name::AbstractString)
    panel = Dict{String, IntRangeSet{i32}}()
    open(rel".." * "/data/panel/$name.panel") do f
        while !eof(f)
            chr = readuntil(f, '\0')[1:end-1]
            len = f >> i64
            buf = read(f, i32, 2len)

            p = panel[chr] = IntRangeSet{i32}()

            for i in 1:len
                push!(p, buf[2i-1]:buf[2i])
            end
        end
    end
    panel
end


text:
##### Problem Set 1 #####

using JLD2
using Random
using LinearAlgebra
using Statistics
using CSV
using DataFrames
using FreqTables
using Distributions

# NOTE: Remember to set path
cd("/Users/matthew/Documents/Summer University of Oklahoma Course/fall-2020/ProblemSetSolutions")

##### Problem 1
### Set the seed
Random.seed!(1234)

### Part a ###
## i.) Create 10x7 matrix of numbers distr. U~[-5, 10]
A = rand(Uniform(-5, 10), 10, 7)

## ii.) Create 10x7 matrix of numbers distr. N~[-2, 15]
B = rand(Normal(-2, 15), 10, 7)

## iii.) Create 5x7 matrix: first 5r and 5c from A and last two columns of B
C = hcat(A[1:5, 1:5], B[1:5, 6:7])

## iv.) D_i,j = A_i,j if A_i,j <= 0, or 0 otherwise
D = min.(A, 0)

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

### Part b ###
## Use built-in Julia function to list the number of elements of A
length(A)
##################################

### Part c ###
## Get number of unique elements of D
length(unique(D))
##################################

### Part d ###
## Use reshape() to create matrix E which is vec operator applied to B
## Can you find an easier way to accomplish this?
E = reshape(B, (length(B), 1))
## Alternative way
vec(B)
##################################

### Part e ###
## Create F, 3-dim that contains A in col 1 of 3rd dim and B in col 2 of 3rd dim
# Can use cat
F = cat(A, B; dims = 3)
##################################

### Part f ###
## Use permutedims() to twist F into a 2x10x7
F = permutedims(F, [3, 2, 1])
##################################

### Part g ###
## Create G, the kronecker product of B and C
## What happens when you try kronecker product of C and F?
G = kron(B, C)
kron(C, F) # This results in an error since F is 3 dimensional
##################################

### Part h ###
## Save matrices A, B, C, D, E, F, and G as .jld file called firstmatrix
# @save from JLD2 package
@save "matrixpractice.jld" A, B, C, D, E, F, G
##################################


text:
function make_regression2(n=100, p=3, c=5;intercept=true, noise=0.1, rng=nothing)
    if rng === nothing
        rng = Random.GLOBAL_RNG
    elseif rng isa Integer
        rng = Random.MersenneTwister(rng)
    end

    X = randn(rng, n, p)
    X = MLJBase.augment_X(randn(rng, n, p), intercept)
    A = randn(rng, p + Int(intercept), c)
    Y = (X * A) .+ (noise .* randn(rng, n, c))
    return MLJBase.table(X), MLJBase.table(Y, names=[Symbol("y$(i)") for i = 1:c])
end

function test_regression(model, X, y)
    fitresult, report, cache = fit(model, 0, X, y)
    yhat = predict(model, fitresult, X)
    fr = fitted_params(model, fitresult)
    return yhat, fr
end

function test_composition_model(ms_model, mlj_model, X, X_array)
    mlj_model_type = typeof(mlj_model)
    Xtr_ms = permutedims(
        MultivariateStats.transform(ms_model, permutedims(X_array))
    )
    fitresult, _, _ = fit(mlj_model, 1, X)
    Xtr_mlj_table = transform(mlj_model, fitresult, X)
    Xtr_mlj = matrix(Xtr_mlj_table)
    # Compare MLJ and MultivariateStats transformed matrices
    @test Xtr_mlj ≈ Xtr_ms
    # test metadata
    d = info_dict(mlj_model_type)
    @test d[:input_scitype] == Table(Continuous)
    @test d[:output_scitype] == Table(Continuous)
    @test d[:name] == string(mlj_model_type)
end


text:
export
  multiplyDistributions,
  status

function status(cfg, fg, requestDict)
  @info "got status request"
  Dict{String, Any}("status" => "OK")
end

function multiplyDistributions(cfg, fg, requestDict)::Dict{String, Any}
  requestVal = Unmarshal.unmarshal(MultiplyDistributionsRequest, requestDict["payload"])

  funcs = BallTreeDensity[]

  for pbd in requestVal.weights
    pts = pbd.points
    @warn "recalculating new bandwidth for multiplyDistributions"
    # TODO assuming 1D Euclidean
    p = AMP.manikde!( reshape(pts, 1, :), (:Euclid,))
    push!(funcs, p)
  end

  pqr = manifoldProduct(funcs, (:Euclid,))

  bw = vec(getBW(pqr))
  n = Npts(pqr)
  rpts = vec(getPoints(pqr))
  resObj = Packed_BallTreeDensity(1, "Gaussian", bw, 1/n*ones(n), rpts)

  Dict{String, Any}("status" => "OK", "payload" => resObj)
end

# t2 = "{\"payload\":{\"weights\":[{\"distType\":\"BallTreeDensity\",\"kernelType\":\"Gaussian\",\"weights\":[1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0],\"bandwidth\":[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],\"points\":[-0.8272759752535751,1.90340735325013,0.5617238720356695,-1.2610130840544629,-0.5196580171571705,0.47808234574398156,-0.5034354386657888,-0.1289061606999511,-0.6283440136328498,-0.04520468574277743],\"dim\":1},{\"distType\":\"BallTreeDensity\",\"kernelType\":\"Gaussian\",\"weights\":[1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0],\"bandwidth\":[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],\"points\":[0.6986754683541115,0.4206784705438769,-0.01798641761511103,0.2862258908507545,0.5021835960347083,-1.4424531814460273,0.7667313775874391,0.13048774810077268,-0.004778323705883176,-0.6408606329686699],\"dim\":1}]},\"request\":\"multiplyDistributions\"}"


text:
# Stack ops are things like Dup, Swap, Flush etc
abstract type PushStackOp{T} <: PushInstructionA1tA1{T} end

struct Dup{T} <: PushStackOp{T} end
@inline function execute!(i::Dup{T}, interp::AbstractPushInterpreter) where {T}
    stackhas(interp, T, 1) || return(nothing)
    s = stack(interp, T)
    push!(s, last(s))
end

struct Swap{T} <: PushStackOp{T} end
@inline function execute!(i::Swap{T}, interp::AbstractPushInterpreter) where {T}
    stackhas(interp, T, 2) || return(nothing)
    a, b = pop2!(interp, T) # a was last and b next to last
    push!(interp, a)
    push!(interp, b) # b is now last so a is next to last
end

struct Flush{T} <: PushStackOp{T} end
@inline function execute!(i::Flush{T}, interp::AbstractPushInterpreter) where {T}
    empty!(stack(interp, T))
end

struct StackDepth{T} <: PushStackOp{T} end
@inline function execute!(i::StackDepth{T}, interp::AbstractPushInterpreter) where {T}
    push!(stack(interp, Int), length(stack(interp, T)))
end


text:
"""
    StdDPP(kernel::Kernel)

Standard DPP (Determinantal Point Process) sampling given `kernel`.
The size of the returned `Z` is not fixed (but cannot be empty unlike in a classical DPP)
"""
struct StdDPP{K<:Kernel} <: OffIPSA
    kernel::K
    function StdDPP(kernel::K) where {K<:Kernel}
        return new{K}(kernel)
    end
end

function inducingpoints(rng::AbstractRNG, alg::StdDPP, X::AbstractVector; kwargs...)
    dpp = DPP(kernel, X)
    samp = rand(rng, dpp)
    while isempty(samp) # Sample from the DPP until there is a non-empty set
        samp = rand(rng, dpp)
    end
    return X[samp]
end


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2020
#
# This is a Julia package for discrete event simulation
#

"""
```
event!([clk], ex, t; <keyword arguments>)
event!([clk], ex, T, t; <keyword arguments>)
```
Schedule an event for a given simulation time.

# Arguments
- `clk<:AbstractClock`: clock, it not supplied, the event is scheduled to 𝐶,
- `ex<:Action`: an expression or function or a tuple of them,
- `T::Timing`: a timing, one of `at`, `after` or `every`,
- `t<:Number`: simulation time, if t < clk.time set t = clk.time,

# Keyword arguments
- `cycle<:Number=0.0`: repeat cycle time for an event,
- `cid::Int=clk.id`: if cid ≠ clk.id, assign the event to the parallel clock
    with id == cid. This overrides `spawn`,
- `spawn::Bool=false`: if true, spawn the event at other available threads,
- `sync::Bool=false`: if true, force a synchronization of all parallel clocks
    before executing the event.

# Examples
```jldoctest
julia> using DiscreteEvents, Unitful

julia> import Unitful: s, minute, hr

julia> myfunc(a, b) = a+b
myfunc (generic function with 1 method)

julia> event!(𝐶, fun(myfunc, 1, 2), 1) # a 1st event to 1

julia> event!(𝐶, fun(myfunc, 2, 3), 1) #  a 2nd event to the same time

julia> event!(𝐶, fun(myfunc, 3, 4), 1s)
Warning: clock has no time unit, ignoring units

julia> setUnit!(𝐶, s)
0.0 s

julia> event!(𝐶, fun(myfunc, 4, 5), 1minute)

julia> event!(fun(myfunc, 5, 6), after, 1hr)

```
"""
function event!(clk::CL, ex::A, t::U; cycle::V=0.0,
                cid::Int=clk.id, spawn::Bool=false, sync::Bool=false) where {CL<:AbstractClock,A<:Action,U<:Number,V<:Number}
    t = _tadjust(clk, t)
    cycle = _tadjust(clk, cycle)
    t = max(t, _tadjust(clk, tau(clk)))

    if cid == clk.id && spawn  # evaluate spawn only if cid == clk.id
        cid = _spawnid(clk)
    end
    _assign(clk, DiscreteEvent(ex, t, cycle), cid)
end
function event!(clk::CL, ex::A, T::Timing, t::U;
                cid::Int=clk.id, spawn::Bool=false, sync::Bool=false) where {CL<:AbstractClock,A<:Action,U<:Number}
    t = _tadjust(clk, t)
    if T == after
        event!(clk, ex, t+clk.time, cid=cid, spawn=spawn, sync=sync)
    elseif T == every
        event!(clk, ex, clk.time, cycle=t, cid=cid, spawn=spawn, sync=sync)
    else
        event!(clk, ex, t, cid=cid, spawn=spawn, sync=sync)
    end
end
event!(ex::A, t::N; kw...) where {A<:Action,N<:Number} = event!(𝐶, ex, t; kw...)
event!(ex::A, T::Timing, t::N; kw...) where {A<:Action,N<:Number} = event!(𝐶, ex, T, t; kw...)

"""
    event!([clk], ex, cond; <keyword arguments>)

Schedule ex as a conditional event, conditions cond get evaluated at each clock tick.

# Arguments
- `clk<:AbstractClock`: if no clock is supplied, the event is scheduled to 𝐶,
- `ex<:Action`: an expression or function or a tuple of them,
- `cond<:Action`: a condition is true if all functions or expressions therein return true,
- `cid::Int=clk.id`: assign the event to the parallel clock cid. This overrides `spawn`,
- `spawn::Bool=false`: if true, spawn the event at other available threads.

# Examples
```jldoctest
julia> using DiscreteEvents

julia> c = Clock()   # create a new clock
Clock 0, thrd 1 (+ 0 ac): state=DiscreteEvents.Undefined(), t=0.0 , Δt=0.0 , prc:0
  scheduled ev:0, cev:0, sampl:0

julia> event!(c, fun((x)->println(tau(x), ": now I'm triggered"), c), fun(>=, fun(tau, c), 5))

julia> c                       # a conditional event turns sampling on  ⬇
Clock 0, thrd 1 (+ 0 ac): state=DiscreteEvents.Undefined(), t=0.0 , Δt=0.01 , prc:0
  scheduled ev:0, cev:1, sampl:0

julia> run!(c, 10)   # sampling is not exact, so it takes 501 sample steps to fire the event
5.009999999999938: now I'm triggered
"run! finished with 0 clock events, 501 sample steps, simulation time: 10.0"

julia> c           # after the event sampling is again switched off ⬇
Clock 0, thrd 1 (+ 0 ac): state=DiscreteEvents.Idle(), t=10.0 , Δt=0.0 , prc:0
  scheduled ev:0, cev:0, sampl:0
```
"""
function event!(clk::T, ex::A, cond::C;
                cid::Int=clk.id, spawn=false) where {T<:AbstractClock,A<:Action,C<:Action}
    if _busy(clk) && all(_evaluate(cond))   # all conditions met
        _evaluate(ex)                      # execute immediately
    else
        if cid == clk.id && spawn  # evaluate spawn only if cid == clk.id
            cid = _spawnid(clk)
        end
        _assign(clk, DiscreteCond(cond, ex), cid)
    end
    return nothing
end
event!( ex::A, cond::C; kw...) where {A<:Action,C<:Action} = event!(𝐶, ex, cond; kw...)

"""
    periodic!([clk], ex, Δt; spawn)

Register a function or expression for periodic execution at the clock`s sample rate.

# Arguments
- `clk<:AbstractClock`: if not supplied, it registers on 𝐶,
- `ex<:Action`: an expression or function or a tuple of them,
- `Δt<:Number=clk.Δt`: set the clock's sampling rate, if no Δt is given, it takes
    the current sampling rate, if that is 0, it calculates one,
- `spawn::Bool=false`: if true, spawn the periodic event to other available threads.
"""
function periodic!(clk::Clock, ex::T, Δt::U=clk.Δt;
                   spawn=false) where {T<:Action,U<:Number}
    clk.Δt = Δt == 0 ? _scale(clk.end_time - clk.time)/100 : Δt
    _assign(clk, Sample(ex), spawn ? _spawnid(clk) : 0)
end
periodic!(ex::T, Δt::U=𝐶.Δt; kw...) where {T<:Action,U<:Number} = periodic!(𝐶, ex, Δt; kw...)
periodic!(ac::ActiveClock, ex::T, Δt::U=ac.clock.Δt; kw...) where {T<:Action,U<:Number} = periodic!(ac.clock, ex, Δt; kw...)

# Return a random number out of the thread ids of all available parallel clocks.
# This is used for `spawn`ing tasks or events to them.
function _spawnid(clk::Clock) :: Int
    if isempty(clk.ac)
        return 0
    else
        return rand(rng, (0, (i.id for i in clk.ac)...))
    end
end

# calculate the scale from a given number
function _scale(n::T)::Float64 where {T<:Number}
    if n > 0
        i = 1.0
        while !(10^i ≤ n < 10^(i+1))
            n < 10^i ? i -= 1 : i += 1
        end
        return 10^i
    else
        return 1.0
    end
end

# ---------------------------------------------------------------
# assign and register events and samples to a clock
# ---------------------------------------------------------------
# There are several ways to do it:
# 1. assign it directly to a clock or an active clock or
# 2. assign it via a channel to another one, if given a different id.
#    In this case the event is sent over the channel to the target clock.
#    - The master clock (id=0) can directly send to an active clock.
#    - An active clock can send directly to master.
#    - An active clock can only send via master to another active clock.
function _assign(c::S, ev::T, id::Int=c.id) where {S<:AbstractClock, T<:AbstractEvent}
    if id == c.id
        _register!(c, ev)       # 1: register it to yourself
    else
        _register(c, ev, id)    # 2: register it to another clock
    end
end

# function barrier: register a concrete event directly to a clock.
function _register!(c::Clock, ev::DiscreteEvent)
    t = ev.t
    while any(i->i==t, values(c.sc.events)) # in case an event at that time exists
        t = nextfloat(float(t))                  # increment scheduled time
    end
    c.sc.events[ev] = t
    return
end
function _register!(c::Clock, cond::DiscreteCond)
    (c.Δt == 0) && (c.Δt = _scale(c.end_time - c.time)/100)
    push!(c.sc.cevents, cond)
    return
end
_register!(c::Clock, sp::Sample) = ( push!(c.sc.samples, sp); nothing )
_register!(ac::ActiveClock, ev::T) where {T<:AbstractEvent} = _register!(ac.clock, ev)

# Register an event to another clock via a channel.
# - c::Clock: a master clock can forward events to active clocks,
# - ac::ActiveClock: active clocks can forward events only through master,
#                    he then does the distribution for them,
# - ev<:AbstractEvent: the event to register,
# - id::Int: the id of the clock it should get registered to.
function _register(c::Clock, ev::T, id::Int) where {T<:AbstractEvent}
    if id > 0 && threadid() == 1       # only master can forward events
        if id ≤ length(c.ac)
            @inbounds put!(c.ac[id].forth, Register(ev))
            return
        end
    end
    _register!(c, ev)                   # otherwise take it yourself
end
function _register(ac::ActiveClock, ev::T, id::Int) where {T<:AbstractEvent}
    if ac.id == id                     # if id is your id
        _register!(ac.clock, ev)       # take it yourself
    elseif ac.clock.state == Busy()    # if master handles its channels
        put!(ac.back, Forward(ev, id)) # put it over the channel
        return
    else
        _register(ac.master[], ev, id) # otherwise call him directly
    end
end


text:
abstract type RestartScheme end
struct NoRestarts <: RestartScheme

end

struct FixedFrequencyRestarts <: RestartScheme
  restart_length::Int64
  always_reset_to_average::Bool # do we alway reset to average?
end

struct AdaptiveRestarts <: RestartScheme
  beta::Float64
  always_reset_to_average::Bool # do we alway reset to average?
end

mutable struct RestartInfo
  last_restart_primal_solution::Vector{Float64}
  last_restart_dual_solution::Vector{Float64}
  normalized_duality_gap::Float64
end

function initialize_restart_info(primal_size::Int64, dual_size::Int64)
  last_restart_primal_solution = zeros(primal_size)
  last_restart_dual_solution = zeros(dual_size)
  normalized_duality_gap = Inf
  return RestartInfo(
    last_restart_primal_solution,
    last_restart_dual_solution,
    normalized_duality_gap,
  )
end

function compute_normalized_duality_gap(
  problem::FirstOrderLp.QuadraticProgrammingProblem,
  primal_solution::Vector{Float64},
  dual_solution::Vector{Float64},
  primal_weight::Float64,
  distance_travelled::Float64,
)
  local_duality_gap = FirstOrderLp.bound_optimal_objective(
    problem,
    primal_solution,
    dual_solution,
    ones(length(primal_solution)) * primal_weight, # primal_norm_params
    ones(length(dual_solution)) / primal_weight, # dual_norm_params
    distance_travelled,
    FirstOrderLp.EUCLIDEAN_NORM, # norm
    solve_approximately = false,
  )
  return FirstOrderLp.get_gap(local_duality_gap) / distance_travelled
end

mutable struct NormalizedDualityGaps
  average_distance_travelled::Float64
  average_gap::Float64
  current_distance_travelled::Float64
  current_gap::Float64
end

function compute_normalized_duality_gaps(
  problem::FirstOrderLp.QuadraticProgrammingProblem,
  restart_info::RestartInfo,
  current_primal_solution::Vector{Float64},
  current_dual_solution::Vector{Float64},
  avg_primal_solution::Vector{Float64},
  avg_dual_solution::Vector{Float64},
  primal_weight::Float64,
)
  distance_travelled_by_avg = norm([
    avg_primal_solution - restart_info.last_restart_primal_solution
    avg_dual_solution - restart_info.last_restart_dual_solution
  ])
  avg_normalized_duality_gap = compute_normalized_duality_gap(
    problem,
    avg_primal_solution,
    avg_dual_solution,
    primal_weight,
    distance_travelled_by_avg,
  )
  distance_travelled_by_current = norm([
    current_primal_solution - restart_info.last_restart_primal_solution
    current_dual_solution - restart_info.last_restart_dual_solution
  ])
  current_normalized_duality_gap = compute_normalized_duality_gap(
    problem,
    current_primal_solution,
    current_dual_solution,
    primal_weight,
    distance_travelled_by_current,
  )
  return NormalizedDualityGaps(
    distance_travelled_by_avg,
    avg_normalized_duality_gap,
    distance_travelled_by_current,
    current_normalized_duality_gap,
  )
end

function select_best_iterate(
  gaps::NormalizedDualityGaps,
  current_primal_solution::Vector{Float64},
  current_dual_solution::Vector{Float64},
  avg_primal_solution::Vector{Float64},
  avg_dual_solution::Vector{Float64},
  always_reset_to_average::Bool,
)
  restart_to_current = false
  if !always_reset_to_average
    if gaps.current_gap < gaps.average_gap
      restart_to_current = true
    end
  end
  if restart_to_current
    candidate_normalized_duality_gap = gaps.current_gap
    candidate_primal_solution = current_primal_solution
    candidate_dual_solution = current_dual_solution
  else
    candidate_normalized_duality_gap = gaps.average_gap
    candidate_primal_solution = avg_primal_solution
    candidate_dual_solution = avg_dual_solution
  end

  return candidate_normalized_duality_gap,
  candidate_primal_solution,
  candidate_dual_solution
end

function run_restart_scheme(
  restart_scheme::RestartScheme,
  problem::FirstOrderLp.QuadraticProgrammingProblem,
  current_primal_solution::Vector{Float64},
  current_dual_solution::Vector{Float64},
  primal_weight::Float64,
  solution_weighted_avg::FirstOrderLp.SolutionWeightedAverage,
  restart_info::RestartInfo,
  verbose::Bool,
)

  do_fixed_frequency_restart =
    isa(restart_scheme, FixedFrequencyRestarts) &&
    solution_weighted_avg.sum_primal_solutions_count >=
    restart_scheme.restart_length
  consider_adaptive_restart = isa(restart_scheme, AdaptiveRestarts)

  primal_avg, dual_avg = FirstOrderLp.compute_average(solution_weighted_avg)
  gaps = compute_normalized_duality_gaps(
    problem,
    restart_info,
    current_primal_solution,
    current_dual_solution,
    primal_avg,
    dual_avg,
    primal_weight,
  )

  do_restart = false
  if do_fixed_frequency_restart || consider_adaptive_restart
    candidate_normalized_duality_gap,
    candidate_primal_solution,
    candidate_dual_solution = select_best_iterate(
      gaps,
      current_primal_solution,
      current_dual_solution,
      primal_avg,
      dual_avg,
      restart_scheme.always_reset_to_average,
    )
    if do_fixed_frequency_restart
      do_restart = true
    elseif consider_adaptive_restart
      if candidate_normalized_duality_gap <=
         restart_scheme.beta * restart_info.normalized_duality_gap
        do_restart = true
      end
    end
  end
  if do_restart
    if verbose
      println("restart algorithm after $(solution_weighted_avg.sum_primal_solutions_count) iterations")
    end
    restart_info.normalized_duality_gap = candidate_normalized_duality_gap
    current_primal_solution .= candidate_primal_solution
    current_dual_solution .= candidate_dual_solution
    restart_info.last_restart_primal_solution .= current_primal_solution
    restart_info.last_restart_dual_solution .= current_dual_solution
    FirstOrderLp.reset_solution_weighted_average(solution_weighted_avg)
  end

  return do_restart, gaps
end

max_stars_repo_name: pedersor/google-research

max_issues_repo_name: pedersor/google-research

max_forks_repo_name: pedersor/google-research


text:
function  getProblemSpecificParameters(nx,ny,nquadpoints,ThisICType,pointsxyz)

	if ThisICType=="linesource"
		SigmaS = 1*ones(ny+4,nx+4);
		SigmaA = 0*ones(ny+4,nx+4);
		SigmaT = SigmaS+SigmaA;
		x0 = -1.5; x1 =  1.5;  y0 = -1.5; y1 =  1.5
		dx = (x1-x0)/(nx)
		dy = (y1-y0)/(ny)
		tEnd = 1.0
		phi0 = setInitialConditionsGaussian(nx,ny,nquadpoints,dx,dy,x0,y0)
		Source = zeros(Float64,nquadpoints,ny+4,nx+4)
	elseif ThisICType=="beam"
		SigmaS = 0.0*ones(ny+4,nx+4);
		SigmaA = 0.0*ones(ny+4,nx+4);
		SigmaT = SigmaS+SigmaA;
		x0 = -1.5; x1 =  1.5;  y0 = -1.5; y1 =  1.5
		dx = (x1-x0)/(nx)
		dy = (y1-y0)/(ny)
		tEnd = 10.0

		beamdirection = [1/sqrt(2),1/sqrt(2),0]
		beamdirection = [0,1,0]
		Source = setInitialConditionsBeam(nx,ny,nquadpoints,dx,dy,x0,y0,pointsxyz,beamdirection)
		phi0= zeros(Float64,nquadpoints,ny+4,nx+4)
		v = 10.0
		for j=3:nx+2
			for i=3:ny+2
				y = y0+dy/2 + (i-3)*dy
				x = x0+dx/2 + (j-3)*dx
				if x>0.0 #&& x <0.1
					SigmaS[i,j] = v
				end

#				if x>0.4 && x<0.5
#					SigmaT[i,j] = v
#				end
#
#
#				if x>0.7 && x<0.8
#					SigmaT[i,j] = v
#				end
			end
		end
		SigmaT = SigmaS+SigmaA;

	elseif ThisICType=="checkerboard"
		x0 = 0.0; x1 = 7.0; y0 = 0.0; y1 = 7.0;
		dx = (x1-x0)/(nx)
		dy = (y1-y0)/(ny)
		SigmaS = zeros(ny+4,nx+4);
		SigmaT = zeros(ny+4,nx+4);
		phi0 = zeros(Float64,nquadpoints,ny+4,nx+4)
		tEnd = 3.2
		Source  = zeros(Float64,nquadpoints,ny+4,nx+4)

		for j=3:nx+2
			for i=3:ny+2
				y = y0+dy/2 + (i-3)*dy
				x = x0+dx/2 + (j-3)*dx
				cx = ceil(x);
				cy = ceil(y);
				g = float(ceil((cx+cy)/2.0)*2.0==(cx+cy))*((1.0<cx)&&(cx<7.)&&(1.0<cy)&&((cy-2.0*abs(cx-4.0))<4.0));
				SigmaS[i,j] = (1-g)*1+g*0;
                SigmaT[i,j] = (1-g)*0+g*10 + SigmaS[i,j]
				if ((3<x)&&(x<4)&&(3<y)&&(y<4))
					for q=1:nquadpoints
						Source[q,i,j] = 1.0;
					end
					SigmaS[i,j] = 0.0;
                    SigmaT[i,j] = 10.0;
				end
			end
		end
	elseif ThisICType=="dim1"
		SigmaS = 1*ones(ny+4,nx+4);
		SigmaA = 0*ones(ny+4,nx+4);
		SigmaT = SigmaS+SigmaA;
		x0 = -1.5; x1 =  1.5;  y0 = -1.5; y1 =  1.5
		dx = (x1-x0)/(nx)
		dy = (y1-y0)/(ny)
		tEnd = 10.0
		phi0 = setInitialConditionsGaussian1d(nx,ny,nquadpoints,dx,dy,x0,y0)
		Source = zeros(Float64,nquadpoints,ny+4,nx+4)
	else
		#warning('Not yet implemented')
	end
    return x0,x1,y0,y1,tEnd,SigmaS,SigmaT,Source,phi0
end

function setInitialConditionsGaussian(nx,ny,nq,dx,dy,domainx0,domainy0)
    # init to zero
    phi = zeros(Float64,nq,ny+4,nx+4)
    x0 = 0.0
    y0 = 0.0
    s2 = 0.03^2
    floor = 1e-4
    f(x,y) = max(floor,1.0/(4.0*pi*s2) *exp(-((x-x0)*(x-x0)+(y-y0)*(y-y0))/4.0/s2));
    for j=3:nx+2
        for i=3:ny+2
            y = domainy0+dy/2 + (i-3)*dy
            x = domainx0+dx/2 + (j-3)*dx
            for q=1:nq
                phi[q,i,j] = f(x,y)/4.0/pi;
            end
        end
    end
    return phi
end

function setInitialConditionsGaussian1d(nx,ny,nq,dx,dy,domainx0,domainy0)
    # init to zero
    phi = zeros(Float64,nq,ny+4,nx+4)
    x0 = 0.0
    y0 = 0.0
    s2 = 0.03^2
    floor = 1e-4
    f(x,y) = max(floor,1.0/(4.0*pi*s2) *exp(-((x-x0)*(x-x0)+(y-y0)*(y-y0))/4.0/s2));
    for j=3:nx+2
        for i=3:ny+2
            y = domainy0+dy/2 + (i-3)*dy
            x = domainx0+dx/2 + (j-3)*dx
            for q=1:nq
                phi[q,i,j] = f(x,0)/4.0/pi;
            end
        end
    end
    return phi
end

function setInitialConditionsBeam(nx,ny,nq,dx,dy,domainx0,domainy0,quadpoints,beamdirection)
    # init to zero
    phi = 1e-14 .*ones(Float64,nq,ny+4,nx+4)
	# beam postion:
    x0 = -0.50
    y0 =  0.00
	# beam width, height
	height = 0.05
	width =  0.05

    for j=3:nx+2
        for i=3:ny+2
            y = domainy0+dy/2 + (i-3)*dy
            x = domainx0+dx/2 + (j-3)*dx
			if x>x0-width && x<x0+width && y>y0-height && y<y0+height # beam should be here!
				# beam is a gaussian with pole pointing into beamdirection
				for q=1:nq
					Omega = quadpoints[q,:]
					eps = 1^2
					x =  dot(Omega,beamdirection)-1
					val = 1.0/sqrt(2.0*pi*eps) *
						exp(-x^2/(2*eps))
					phi[q,i,j] = val
				end
			end
        end
    end
    return phi
end

max_stars_repo_name: camminady/SN

max_issues_repo_name: camminady/SN

max_forks_repo_name: camminady/SN


text:
module Utils
export get_starttimelist, initlogo
using Dates
"""
get_starttimelist(st::DateTime, et::DateTime, unittime::Float64)
calculate start time list for parallel downloading

    st: start time
    et: end time
    unittime: unit time in Second

    this function returns
    stlist: list of start time

    e.g.
    st = DateTime(2019,1,1,0,0,0)
    et = DateTime(2019,1,1,12,0,0)
    unittime = 3600

    stlist = get_starttimelist(st, et, unittime)

"""
function get_starttimelist(st::DateTime, et::DateTime, unittime::Real)

    reftime = st
    stlist = []

    while reftime <= et
        push!(stlist, string(reftime))
        reftime += Dates.Second(float(unittime))
    end

    return stlist
end

"""
initlogo()
print initial logo
"""
function initlogo()

    print("

      _____        _       _____                          _                    _
     / ____|      (_)     |  __ \\                        | |                  | |
    | (___    ___  _  ___ | |  | |  ___ __      __ _ __  | |  ___    __ _   __| |
     \\___ \\  / _ \\| |/ __|| |  | | / _ \\\\ \\ /\\ / /| '_ \\ | | / _ \\  / _` | / _` |
     ____) ||  __/| |\\__ \\| |__| || (_) |\\ V  V / | | | || || (_) || (_| || (_| |
    |_____/  \\___||_||___/|_____/  \\___/  \\_/\\_/  |_| |_||_| \\___/  \\__,_| \\__,_|
                      _         _  _
                     | |       | |(_)           |
    __      __       | | _   _ | | _   __ _     | v1.0.0 (Last update 05/21/2019)
    \\ \\ /\\ / /   _   | || | | || || | / _` |    | © Kurama Okubo
     \\ V  V /_  | |__| || |_| || || || (_| |    |
      \\_/\\_/(_)  \\____/  \\__,_||_||_| \\__,_|    |

")

    println("Job start running at "*string(now())*"\n")

end

end

max_stars_repo_name: kura-okubo/SeisDownloadwJulia

max_issues_repo_name: kura-okubo/SeisDownloadwJulia

max_forks_repo_name: kura-okubo/SeisDownloadwJulia


text:
struct Ctypespec{T, CO<:Copaques, S<:Calignment, TS<:Tuple}
end

Ctypespec(::Type{T}) where {T} = Tuple{T}
Ctypespec(::Type{T}, ::Val{bits}) where {T, bits} = Tuple{T, bits}
Ctypespec(::Type{CO}) where {CO<:Copaques} = Ctypespec{concrete(CO), concrete(CO), strategy(CO), specification(CO)}

Ctypespec(::Type{CO}, ::Type{S}, ::Type{TS}) where {CO<:Copaques, S, TS} = Ctypespec{CO, CO, S, TS}

type(::Type{Ctypespec{T, CO, S, TS}}) where {T, CO, S, TS} = T
opaque(::Type{Ctypespec{T, CO, S, TS}}) where {T, CO, S, TS} = CO
strategy(::Type{Ctypespec{T, CO, S, TS}}) where {T, CO, S, TS} = S
specification(::Type{Ctypespec{T, CO, S, TS}}) where {T, CO, S, TS} = TS

concrete(::Type{T}) where {T} = T
concrete(::Type{CO}) where {CO<:Copaques} = error("Attempted to get type specification details of an opaque type `$(CO)`")
strategy(::Type{CO}) where {CO<:Copaques} = strategy(concrete(CO))
specification(::Type{CO}) where {CO<:Copaques} = specification(concrete(CO))

kind(::Type{TS}) where {TS<:Ctypespec} = kind(opaque(TS))
kind(::Type{CS}) where {CS<:Cstruct} = Cstruct
kind(::Type{CU}) where {CU<:Cunion} = Cunion
kind(::Type{CE}) where {CE<:Cenum} = Cenum


text:
# """
#     abstract type AbstractConfig{T <: SSDFloat} end
#
# Supertype of all detector/world/object configs.
#
# User defined geometries must be subtype of `AbstractConfig{T}`.
#
# There are a few functions which must be defined for a user config, e.g. `struct UserConfig{T} <: AbstractConfig{T}`:
#
# For cylindrical grids:
#
# * in(pt::CylindricalPoint{T}, config::UserConfig{T})::Bool where {T <: SSDFloat}
# * Grid(config::UserConfig{T})::CylindricalGrid{T} where {T <: SSDFloat}
# * get\\_ρ\\_and\\_ϵ(pt::CylindricalPoint{T}, config::UserConfig{T})::Tuple{T, T} where {T <: SSDFloat}
# * set\\_point\\_types\\_and\\_fixed\\_potentials!(point_types::Array{PointType, 3}, potential::Array{T, 3},
#         grid::CylindricalGrid{T}, config::UserConfig{T}; weighting\\_potential\\_contact\\_id::Union{Missing, Int} = missing)::Nothing where {T <: SSDFloat}
#
# For cartesian grids:
#
# * in(pt::CartesianPoint{3, T}, config::UserConfig{T})::Bool
# * Grid(config::UserConfig{T})::CartesianGrid3D{T} where {T <: SSDFloat}
# * get\\_ρ\\_and\\_ϵ(pt::CartesianPoint{3, T}, config::UserConfig{T})::Tuple{T, T} where {T <: SSDFloat}
# * set\\_point\\_types\\_and\\_fixed\\_potentials!(point_types::Array{PointType, 3}, potential::Array{T, 3},
#         grid::CartesianGrid3D{T}, config::UserConfig{T}; weighting\\_potential\\_contact\\_id::Union{Missing, Int} = missing)::Nothing where {T <: SSDFloat}
#
# """
abstract type AbstractConfig{T <: SSDFloat} end


text:
function DopamineDQNBaseline(env; minimal_action_set = false)
    # Replicate Dopamine's Results!

    γ=0.99
    batch_size=32
    buffer_size = 1000000
    tn_update_freq= 8000
    hist_length = 4
    update_freq = 4
    min_mem_size = 20000

    learning_rate = 0.00025
    momentum_term = 0.00
    squared_grad_term = 0.95
    min_grad_term = 1e-5

    example_state = MinimalRLCore.get_state(env)
    action_set = if minimal_action_set
        get_minimal_actions(env)
    else
        get_actions(env)
    end

    init_f = Flux.glorot_uniform

    model = Chain(
        Conv((8,8), hist_length=>32, relu, stride=4, init=init_f),
        Conv((4,4), 32=>64, relu, stride=2, init=init_f),
        Conv((3,3), 64=>64, relu, stride=1, init=init_f),
        (x)->reshape(x, :, size(x, 4)),
        Dense(3136, 512, relu, initW=init_f),
        Dense(512, length(action_set), identity, initW=init_f)) |> gpu

    target_network  = deepcopy(model)

    return DQNAgent(
        model,
        target_network,
        QLearningHuberLoss(γ),
        DeepRL.RMSPropTFCentered(learning_rate,
                                 squared_grad_term,
                                 momentum_term,
                                 min_grad_term),
        DeepRL.ϵGreedyDecay((1.0, 0.01), 250000, min_mem_size, action_set),
        buffer_size,
        hist_length,
        example_state,
        batch_size,
        tn_update_freq,
        update_freq,
        min_mem_size;
        hist_squeeze = Val{false}(),
        state_preproc = DeepRL.image_manip_atari,
        state_postproc = DeepRL.image_norm,
        rew_transform = (r)->clamp(Float32(r), -1.0f0, 1.0f0),
        device = Flux.use_cuda[] ? Val{:gpu}() : Val{:cpu}())
end

function RevisitingALEDQNBaseline(env)
    # Replicate Dopamine's Results!

    γ=0.99
    batch_size=32
    buffer_size = 1000000
    tn_update_freq= 10000
    hist_length = 4
    update_freq = 4
    min_mem_size = 50000

    learning_rate = 0.00025
    momentum_term = 0.00
    squared_grad_term = 0.95
    min_grad_term = 1e-5

    example_state = MinimalRLCore.get_state(env)

    init_f = Flux.glorot_uniform

    model = Chain(
        Conv((8,8), hist_length=>32, relu, stride=4, init=init_f),
        Conv((4,4), 32=>64, relu, stride=2, init=init_f),
        Conv((3,3), 64=>64, relu, stride=1, init=init_f),
        (x)->reshape(x, :, size(x, 4)),
        Dense(3136, 512, relu, initW=init_f),
        Dense(512, length(get_actions(env)), identity, initW=init_f)) |> gpu

    target_network  = deepcopy(model)

    return DQNAgent(
        model,
        target_network,
        QLearningHuberLoss(γ),
        DeepRL.RMSPropTFCentered(learning_rate,
                                 squared_grad_term,
                                 momentum_term,
                                 min_grad_term),
        DeepRL.ϵGreedyDecay((1.0, 0.01), 1000000, min_mem_size, get_actions(env)),
        buffer_size,
        hist_length,
        example_state,
        batch_size,
        tn_update_freq,
        update_freq,
        min_mem_size;
        hist_squeeze = Val{false}(),
        state_preproc = DeepRL.image_manip_atari,
        state_postproc = DeepRL.image_norm,
        rew_transform = (r)->clamp(Float32(r), -1.0f0, 1.0f0),
        device = Flux.use_cuda[] ? Val{:gpu}() : Val{:cpu}())
end


text:
#=
Created on Monday 12 August 2019
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Efficient ways of storing Kronecker powers.
=#

"""
Efficient way of storing Kronecker powers, e.g.

K = A ⊗ A ⊗ ... ⊗ A.
"""
struct KroneckerPower{T<:Any,TA<:AbstractMatrix{T}, N} <: AbstractKroneckerProduct{T}
   A::TA
   pow::Integer
   function KroneckerPower(A::AbstractMatrix{T}, pow::Integer) where {T}
      @assert pow ≥ 2 "KroneckerPower only makes sense for powers greater than 1"
      return new{eltype(A), typeof(A), pow}(A, pow)
    end
end

"""
    kronecker(A::AbstractMatrix, pow::Int)

Kronecker power, computes `A ⊗ A ⊗ ... ⊗ A`. Returns a lazy `KroneckerPower`
type.
"""
kronecker(A::AbstractMatrix, pow::Integer) = KroneckerPower(A, pow)

"""
    ⊗(A::AbstractMatrix, pow::Int)

Kronecker power, computes `A ⊗ A ⊗ ... ⊗ A`. Returns a lazy `KroneckerPower`
type.
"""
⊗(A::AbstractMatrix, pow::Integer) = kronecker(A, pow)

getmatrices(K::KroneckerPower{T,TA,N}) where {T,TA,N} = (K.A, KroneckerPower(K.A, K.pow-1))
getmatrices(K::KroneckerPower{T,TA,2}) where {T,TA} = (K.A, K.A)
getmatrices(K::KroneckerPower{T,TA,1}) where {T,TA} = (K.A, )

order(K::KroneckerPower) = K.pow
Base.size(K::KroneckerPower) = size(K.A).^K.pow
Base.eltype(K::KroneckerPower{T,TA,N}) where {T,TA,N} = T
issquare(K::KroneckerPower) = issquare(K.A)

# SCALAR EQUIVALENTS FOR AbstractKroneckerProduct

"""
    det(K::KroneckerPower)

Compute the determinant of a Kronecker power.
"""
function LinearAlgebra.det(K::KroneckerPower)
    squarecheck(K)
    A, pow = K.A, K.pow
    n = size(A, 1)
    return det(K.A)^(n * pow)
end

"""
    logdet(K::KroneckerPower)

Compute the logarithm of the determinant of a Kronecker power.
"""
function LinearAlgebra.logdet(K::KroneckerPower)
    squarecheck(K)
    A, pow = K.A, K.pow
    n = size(A, 1)
    return n * pow * logdet(K.A)
end

"""
    tr(K::KroneckerPower)

Compute the trace of a Kronecker power.
"""
function LinearAlgebra.tr(K::KroneckerPower)
    squarecheck(K)
    return tr(K.A)^K.pow
end

# Matrix operations

"""
    inv(K::KroneckerPower)

Compute the inverse of a Kronecker power.
"""
function Base.inv(K::KroneckerPower)
    squarecheck(K)
    return KroneckerPower(inv(K.A), K.pow)
end

"""
    adjoint(K::KroneckerPower)

Compute the adjoint of a Kronecker power.
"""
function Base.adjoint(K::KroneckerPower)
    return KroneckerPower(K.A', K.pow)
end

"""
    transpose(K::KroneckerPower)

Compute the transpose of a Kronecker power.
"""
function Base.transpose(K::KroneckerPower)
    return KroneckerPower(transpose(K.A), K.pow)
end

"""
    conj(K::KroneckerPower)

Compute the conjugate of a Kronecker power.
"""
function Base.conj(K::KroneckerPower)
    return KroneckerPower(conj(K.A), K.pow)
end

# mixed-product property
function Base.:*(K1::KroneckerPower{T,TA,N},
                        K2::KroneckerPower{S,TB,N}) where {T,TA,S,TB,N}
    if size(K1, 2) != size(K2, 1)
        throw(DimensionMismatch("Mismatch between K1 and K2"))
    end
    return KroneckerPower(K1.A * K2.A, N)
end


text:
function seriesTable(tabda, idfield; fmod = x->x)

  fieldvec = [ fmod(x[idfield]) for x in tabda];

  xfields = unique(fieldvec);

  if(size(xfields,1) == size(fieldvec,1))

    println("Series in $idfield not found");

    return [];

    end;

  n = size(tabda, 1);

  m = size(xfields,1);

  lookfield = Dict{Any, Int}();

  sizehint!(lookfield, m);

  [lookfield[x] = y for (y,x) in enumerate(xfields)];

  xcounts = zeros(Int, m);

  for i = 1:n

    xcounts[lookfield[fieldvec[i]]] += 1;

    end;

  q = maximum(xcounts);

  series = similar(tabda, m, q); ##Array{NamedTuple}(undef, m, q);

  ##series[:] .= (xundef = true, );

  xcounts[:] .= 0;

  for iter = 1:n

    id = lookfield[fieldvec[iter]];

    series[id, xcounts[id]+1] = tabda[iter];

    xcounts[id] += 1;

    end;

  return (series, xcounts);

  end;

function splitTable(tabda, xcounts)

  ret = Array{typeof(tabda[:,1])}(undef, 0); #[];

  n = size(tabda,1);

  for iter = 1:n

    push!(ret, tabda[iter, 1:xcounts[iter]]);

    end;

  return ret;

  end;


text:
"""
    mutable struct Blending
        direction::Tuple{Luxor.Point,Luxor.Point}
        hue::Tuple{Luxor.RGB,Luxor.RGB}
    end

Color gradient from the starting to stopping position of each waterfall in the cascade.
The gradient begins with the "natural" color (the color that would be used if a gradient
was not applied) and ends in the lightened color.

Blending is applied to `Plot{Vertical}`, for which sample information does not overalp.
If a waterfall contains both positive and negative values, gradients will be defined and
applied separately, so that **all** positive values have the same starting/stopping gradient
position and/or color and **all** negative values have the same starting/stopping gradient
position and/or color.

Arguments:
- `direction::Tuple{Luxor.Point,Luxor.Point}`, the gradient's starting and stopping points.
- `hue::Tuple{Luxor.RGB,Luxor.RGB}`, the gradient's starting and stopping colors.
"""
mutable struct Blending
    direction::Tuple{Luxor.Point,Luxor.Point}
    hue::Tuple{Luxor.RGB,Luxor.RGB}
end


text:
"""
    ShowText([io; kwargs...]) <: ShowBackend

A backend for showing block data using text for REPL use.
Text is displayed to `io` and `kwargs` are keyword arguments
for `PrettyTables.pretty_table`, which is used to display
collections of blocks.
"""
struct ShowText <: ShowBackend
    io
    kwargs
end

ShowText(io=stdout; hlines=:all, alignment=:l, kwargs...) = ShowText(
    io,
    (; hlines=hlines, alignment=alignment, kwargs...))

createhandle(backend::ShowText) = backend.io

function showblock!(io, backend::ShowText, (title, block)::Pair, obs)
    printstyled(io, title, bold=true)
    println(io)
    showblock!(io, backend, block, obs)
end

function showblock!(io, backend::ShowText, blocks::Tuple, obss::Tuple)
    header = [block isa Pair ? first(block) : "" for block in blocks]
    blocks = [block isa Pair ? last(block) : block for block in blocks]
    data = reshape([PrettyTables.AnsiTextCell(io -> showblock!(IOContext(io, :color => true), backend, block, obs))
                for (block, obs) in zip(blocks, obss)], 1, :)
    pretty_table(io, data; header=header, noheader=all(isempty, header), backend.kwargs...)
end

function showblocks!(io, backend::ShowText, blocks::Tuple, obsss::AbstractVector)
    header = [block isa Pair ? first(block) : "" for block in blocks]
    blocks = [block isa Pair ? last(block) : block for block in blocks]
    rows = []
    for obss in obsss
        row = reshape([AnsiTextCell(
                    io -> showblock!(IOContext(io, :color => true), backend, block, obs))
                for (block, obs) in zip(blocks, obss)], 1, :)
        push!(rows, row)
    end
    tabledata = reduce(vcat, rows)
    pretty_table(io, tabledata; header=header, noheader=all(isempty, header), backend.kwargs...)
end

showblocks!(io, backend::ShowText, block, obss::AbstractVector) =
    showblocks!(io, backend, (block,), map(obs -> (obs,), obss))

# Block implementations

function showblock!(io, ::ShowText, block::Label, obs)
    print(io, obs)
end

function showblock!(io, ::ShowText, block::Continuous, obs)
    print(io, obs)
end

function showblock!(io, ::ShowText, block::LabelMulti, obs)
    print(io, obs)
end

function showblock!(io, ::ShowText, block::OneHotLabel, obs)
    if !(sum(obs) ≈ 1)
        obs = softmax(obs)
    end
    obs = round.(obs; sigdigits=3)
    plot = UnicodePlots.barplot(block.classes, obs, width=20, compact=true)
    print(IOContext(io, :color => true), plot)
end


text:
############################################################################################
"""
$(TYPEDEF)

Default arguments for PMCMC constructor.

# Fields
$(TYPEDFIELDS)
"""
struct PMCMCDefault
    function PMCMCDefault()
        return new()
    end
end

############################################################################################
"""
$(TYPEDEF)

Particle MCMC container, container kernel and tuning.

# Fields
$(TYPEDFIELDS)
"""
struct PMCMC{M<:PMCMCKernel,N<:PMCMCTune} <: AbstractAlgorithm
    "PMCMC sampler"
    kernel::M
    "Tuning configuration for kernel."
    tune::N
    function PMCMC(kernel::M, tune::N) where {M<:PMCMCKernel,N<:PMCMCTune}
        return new{M,N}(kernel, tune)
    end
end
function PMCMC(
    _rng::Random.AbstractRNG,
    kernel::Type{M},
    pf::ParticleFilter,
    mcmc::MCMC,
    default::PMCMCDefault=PMCMCDefault(),
    info::BaytesCore.SampleDefault = BaytesCore.SampleDefault()
) where {M<:PMCMCKernel}
    ## Assign PMCMC Kernel
    pmcmc = kernel(pf, mcmc)
    ## Assign tuning container ~ Placeholder for now
    tune = PMCMCTune(mcmc.tune.tagged)
    ## Return struct
    return PMCMC(pmcmc, tune)
end

############################################################################################
"""
$(SIGNATURES)
Propose new parameter with pmcmc sampler. If update=true, objective function will be updated with input model and data.

# Examples
```julia
```

"""
function propose!(
    _rng::Random.AbstractRNG,
    pmcmc::PMCMC,
    model::ModelWrapper,
    data::D,
    temperature::F = model.info.reconstruct.default.output(1.0),
    update::U=BaytesCore.UpdateTrue(),
) where {D,F<:AbstractFloat, U<:BaytesCore.UpdateBool}
    ## Make PMCMC Proposal step
    val, diagnostics = propose!(_rng, pmcmc.kernel, model, data, temperature, update)
    ## Pack and return output
    return val, diagnostics
end

############################################################################################
#export
export PMCMC, PMCMCDefault, propose!


text:
# Generic evaluable

abstract type Evaluable{T} end

arguments(self::Evaluable) = ()
isconstant(::Any) = true
isconstant(self::Evaluable) = all(isconstant(arg) for arg in arguments(self))
iselconstant(::Any) = true
iselconstant(self::Evaluable) = all(iselconstant(arg) for arg in arguments(self))
restype(self::Evaluable{T}) where T = T
typetree(self::Evaluable) = [typeof(self).name.name, (typetree(arg) for arg in arguments(self))...]

# For internal use by the compiler
mutable struct FunctionData
    func :: Evaluable
    arginds :: Vector{Int}
    tgtsym :: Union{Symbol,Nothing}
    allocexprs :: Vector{Any}
    allocsyms :: Vector{Symbol}
    master_alloc :: Union{Symbol,Nothing}

    # Note: Many fields are left undefined
    function FunctionData(func::Evaluable)
        retval = new(func)
        retval.master_alloc = nothing
        retval
    end
end

function linearize(self::Evaluable)
    indices = OrderedDict{Evaluable,Int}()
    linearize!(indices, self)
    sequence = [FunctionData(func) for (func, index) in indices]
    for data in sequence
        data.arginds = Int[indices[arg] for arg in arguments(data.func)]
    end
    sequence
end

function linearize!(indices::OrderedDict{Evaluable,Int}, self::Evaluable)
    haskey(indices, self) && return
    for func in arguments(self)
        linearize!(indices, func)
    end
    indices[self] = length(indices) + 1
end

function ast(self::Evaluable)
    sequence = linearize(self)
    ast!(length(sequence), Set{Int}(), linearize(self))
end

function ast!(index::Int, seen::Set{Int}, sequence::Vector{FunctionData})
    index in seen && return string("%", string(index))

    s = ""
    argids = Stateful(sequence[index].arginds)
    for argid in argids
        bridge = isdone(argids) ? flatten((("└─",), repeated("  "))) : flatten((("├─",), repeated("│ ")))
        sublines = split(ast!(argid, seen, sequence), "\n")
        subtree = join((string(b, s) for (b, s) in zip(bridge, sublines)), "\n")
        s = string(s, "\n", subtree)
    end

    push!(seen, index)
    string("%", string(index), " = ", repr(sequence[index].func), s)
end

# Array functions

abstract type ArrayEvaluable{T,N} <: Evaluable{Array{T,N}} end

Base.eltype(::ArrayEvaluable{T}) where T = T

asarray(v::ArrayEvaluable) = v
asarray(v::Real) = Constant(v)
asarray(v::AbstractArray) = Constant(v)

Base.length(self::ArrayEvaluable) = prod(size(self))
Base.ndims(self::ArrayEvaluable{T,N}) where {T,N} = N :: Int
Base.size(self::ArrayEvaluable, dim::Int) = size(self)[dim]
Base.show(io::IO, self::ArrayEvaluable) = print(io, string(typeof(self).name.name), size(self))

isnormal(::ArrayEvaluable) = true
separate(self::ArrayEvaluable) = [(Tupl(Constant(collect(1:n)) for n in size(self)), self)]

# Compiled functions

struct Compiled{T}
    kernel :: Function
end

callable(self::Compiled) = Base.invokelatest(self.kernel)

struct CompiledArray{T,N}
    kernel :: Function
    shape :: Shape
end

Base.size(self::CompiledArray) = self.shape
callable(self::CompiledArray) = Base.invokelatest(self.kernel)

struct CompiledSparseArray{T,N}
    ikernel :: Compiled
    dkernel :: CompiledArray{T,N}
    shape :: Shape
    blockshape :: Shape
end

Base.size(self::CompiledSparseArray) = self.shape
callable(self::CompiledSparseArray) = (callable(self.ikernel), callable(self.dkernel))

# Compilation

function _compile(infunc::Evaluable, show::Bool)
    sequence = linearize(infunc)
    for data in sequence
        data.tgtsym = gensym()
        data.allocexprs = prealloc(data.func, [sequence[argid].master_alloc for argid in data.arginds])
        data.allocsyms = Symbol[gensym() for _ in 1:length(data.allocexprs)]
        if !isempty(data.allocsyms)
            data.master_alloc = data.allocsyms[end]
        end
    end

    # Code for allocation
    alloccode = Vector{Expr}()
    for data in sequence, (sym, expr) in zip(data.allocsyms, data.allocexprs)
        push!(alloccode, :($sym = $expr))
    end

    # Code for evaluation
    evalcode = Vector{Expr}()
    for data in sequence
        argsyms = Symbol[sequence[argid].tgtsym for argid in data.arginds]
        tgtsym = data.tgtsym
        push!(evalcode, :($(data.tgtsym) = $(codegen(data.func, argsyms, data.allocsyms))))
    end

    typeinfo = [(:point, Vector{Float64}), (:element, Element)]
    paramlist = (:($sym::$tp) for (sym, tp) in typeinfo)
    push!(alloccode, Expr(:function, Expr(:call, :evaluate, paramlist...), Expr(:block, evalcode...)))
    push!(alloccode, :(return evaluate))

    definition = Expr(:function, Expr(:call, :mkevaluate), Expr(:block, alloccode...))

    show && @show definition

    mod = Module()
    Core.eval(mod, :(import Jutils.Elements: Element))
    Core.eval(mod, :(import Jutils.Transforms: applytrans, applytrans_grad))
    Core.eval(mod, :(using LinearAlgebra))
    Core.eval(mod, :(import TensorOperations))
    Core.eval(mod, definition)

    mod.mkevaluate
end

compile(func::Evaluable{T}; show::Bool=false) where T = Compiled{T}(Base.invokelatest(_compile, func, show))

function compile(func::ArrayEvaluable{T,N}; show::Bool=false, dense::Bool=true) where {T,N}
    if dense
        CompiledArray{T,N}(Base.invokelatest(_compile, Normalize(func), show), size(func))
    else
        blocks = separate(func)
        @assert length(blocks) == 1
        inds, data = blocks[1]
        inds = Tupl(Normalize(k) for k in inds)
        ikernel = compile(inds; show=show)
        dkernel = compile(Normalize(data); show=show)
        CompiledSparseArray{T,N}(ikernel, dkernel, size(func), size(data))
    end
end

max_stars_repo_name: TheBB/Jutils

max_issues_repo_name: TheBB/Jutils

max_forks_repo_name: TheBB/Jutils


text:
"""
    autoimports(pkgdir::AbstractString)
    autoimports(paths::AbstractString...)
    autoimports(paths::AbstractVector{<:AbstractString})

# Keyword Arguments
- `project::AbstractString`
- `dry_run::Bool = false`: Do not write to files (imples `verbose=true`).
- `verbose::Bool = false`: Print operations to be performed.
- `sort::Bool = true`: Sort `using`s and `import`s.
"""
autoimports(pkgdir::AbstractString = find_pkgdir(); kwargs...) = autoimports(
    isdir(pkgdir) ? sort!(jlfiles(joinpath(pkgdir, "src"))) : [pkgdir];
    project = isdir(pkgdir) ? joinpath(pkgdir, "Project.toml") : Base.active_project(),
    kwargs...,
)

module __EmptyModule end

unbound_vars_for(vars) = setdiff(
    vars.right,
    vars.left,
    names(Base),
    names(Core),
    names(__EmptyModule; all = true),
)

function _sortnames!(formatted, paths)

#! format: off
    files = sorted_files!(
        let code = get(formatted, p, nothing)
            code === nothing ? SourceCode(p) : SourceCode(p, code)
        end
        for p in paths
    )
#! format: on

    for (p, code) in format_files_to_dict(files)
        formatted[p] = code
    end
    # Note: not using the result of `format_files_to_dict(files)`
    # directly to avoid skipping edits existing in `formatted`.

    return formatted
end

function _write_formatted(formatted; dry_run, verbose)
    if verbose
        @info "Updating files..." files = sort!(collect(keys(formatted)))
    end

    dry_run && return

    # Write results to files
    for (p, code) in formatted
        write(p, code)
    end
end

function autoimports(
    paths::AbstractArray{<:AbstractString};
    project::AbstractString = Base.active_project(),
    dry_run::Bool = false,
    verbose::Bool = false,
    sort::Bool = true,
)
    verbose = verbose | dry_run
    project = as_project_toml_path(project)

    if verbose
        @info "Analyzing and modifying $(length(paths)) files." project sort dry_run
    end
    if !occursin("AutoImports", ENV["JULIA_DEBUG"])
        ENV["JULIA_DEBUG"] = ENV["JULIA_DEBUG"] * ",AutoImports"
    end
    @debug "Debug logging enabled" ENV["JULIA_DEBUG"]

    # Find unbound variables
    vars = module_variables(paths)
    unbound_vars = unbound_vars_for(vars)
    @debug "Unbound variables found" unbound_vars = Text(join(unbound_vars, ", "))

    # List variables defined in deps
    depsvars = depsvars_from_tomlpath(project, imported_modules(paths))

    # Find packages that define `unbound_vars`
    tobeimported = Dict{Base.PkgId,Vector{Symbol}}()
    for v in unbound_vars
        for (pkg, dvars) in depsvars
            if v in dvars.left
                push!(get!(tobeimported, pkg, []), v)
                @goto ok
            end
        end
        @warn "Variable `$v` not found anywhere."
        @label ok
    end
    isempty(tobeimported) && return sortnames(paths; dry_run = dry_run, verbose = verbose)

    if verbose
        import_texts = sprint() do io
            for (pkg, vars) in sort!(
                collect(tobeimported),
                by = ((pkg, _),) -> pkg.name,
            )
                print(io, "using $(pkg.name): ")
                join(io, vars, ", ")
                println(io)
            end
        end
        @info """
        Required explicit imports:
        $import_texts"""
    end

    @debug "Adding explicit imports..." tobeimported
    files = add_using!(paths, tobeimported)

    formatted = format_files_to_dict(files)

    if sort
        _sortnames!(formatted, paths)
    end

    _write_formatted(formatted; dry_run = dry_run, verbose = verbose)
end

"""
    sortnames(path::AbstractString = ".")
    sortnames(paths::AbstractString...)
    sortnames(paths::AbstractVector{<:AbstractString})

# Keyword Arguments
- `dry_run::Bool = false`: Do not write to files (imples `verbose=true`).
- `verbose::Bool = false`: Print operations to be performed.
"""
sortnames(path::AbstractString = "."; kwargs...) =
    sortnames(isdir(path) ? sort!(jlfiles(path)) : [path]; kwargs...)

function sortnames(
    paths::AbstractArray{<:AbstractString};
    dry_run::Bool = false,
    verbose::Bool = false,
)
    verbose = verbose | dry_run
    formatted = _sortnames!(Dict{String,String}(), paths)
    _write_formatted(formatted; dry_run = dry_run, verbose = verbose)
end

@nospecialize
autoimports(paths::AbstractString...; kwargs...) = autoimports(collect(paths); kwargs...)
sortnames(paths::AbstractString...; kwargs...) = sortnames(collect(paths); kwargs...)
@specialize

sortnames_text(code::AbstractString) = sprint(sortnames_text, code)
function sortnames_text(io::IO, code::AbstractString)
    f = sort_names!(SourceCode("<string>", code))
    format_code(io, f)
end


text:
using LinearAlgebra.BLAS
using OpenBLASHighCoreCount
using OpenBLASHighCoreCount_jll
using Test

@testset "OpenBLASHighCoreCount test suite" begin
    libs = BLAS.get_config().loaded_libs
    idx = findfirst(l -> l.libname == OpenBLASHighCoreCount_jll.libopenblas_path, libs)
    @test idx !== nothing
    @test libs[idx].interface === :ilp64
    BLAS.set_num_threads(128)
    @test BLAS.get_num_threads() == 128
end


text:
@recipe function f(l::Line;numerate=false,m=false)
    # defaults
    grid --> true
    xticks := true
    legend := false
    @series begin
        linetype --> :path
        linecolor --> :black
        linewidth --> 1
        if m
            markershape --> :circle
            markerstrokecolor --> :black
            markercolor --> :white
            markersize --> 2
            if numerate
                series_annotations := ["$i" for i in 1:length(l.x)]
            end
        end
        (l.x,l.y)
    end
end

@recipe function f(x::Envelope; removed=false,numerate=false,m=false)

    # defaults
    grid --> true
    xticks := true
    legend := false

    # if line array exists, plot
    if length(x.L) > 0
        for l in x.L
            @series begin
                # subplot := 1
                linetype := :line
                linecolor := :black
                linewidth := 1
                if m
                    markershape := :circle
                    markerstrokecolor := :black
                    markercolor := :white
                    markersize := 2
                    if numerate
                        series_annotations := ["$i" for i in 1:length(l.x)]
                    end
                end
                (l.x,l.y)
            end
        end
    end
    # plot envelope, if exists
    if x.env_set
        @series begin
            # subplot := 1
            linetype := :line
            linecolor --> :red
            linewidth --> 4
            if m
                markershape := :circle
                markercolor := :white
                # markeralpha := 0.5
                markerstrokecolor := :black
                markersize := 3
                if numerate
                    series_annotations := ["$i" for i in 1:length(getx(x))]
                end
            end
            (getx(x),gety(x))
        end
        if removed
            for ir in x.removed
                if length(ir) > 0
                    @series begin
                        seriestype = :scatter
                        markershape := :rect
                        markersize := 3
                        markerstrokecolor := :black
                        markercolor := :white
                        markeralpha := 0.5
                        [ir[i].x for i in 1:length(ir)],[ir[i].y for i in 1:length(ir)]
                    end
                end
            end
        end
    end
end

function tplot1()
    n = 15
    x1 = collect(linspace(0,10,n))
    x2 = collect(linspace(-1,9,n))
    L1 = Line(x1,x1)
    L2 = Line(x2,ones(n)*5)
    e = Envelope([L1,L2])
    plot(e)
end

function tplot2()
    n = 15
    x1 = collect(linspace(0,10,n))
    x2 = collect(linspace(-1,9,n))
    L1 = Line(x1,x1)
    L2 = Line(x2,ones(n)*5)
    en = Envelope([L1,L2])
    p1 = plot(en)

    upper_env!(en)
    p2 = plot(en,removed=true)
    plot(p1,p2)
end

function f3a()

    fs = Function[x->ones(length(x));x->0.5x;x->x-2;x->2x-8]
    xs = [
         collect(linspace(-1,0.9,6)),
         collect(linspace(1,7,19)),
         collect(linspace(2,7,15)),
         collect(linspace(4,8,25))]
    ls = [Line(i[1],i[2](i[1])) for i in zip(xs,fs)]
    # create an envelope
    e = Envelope(ls)
    return e
end
function f3b()

    fs = Function[x->ones(length(x));x->0.5x;x->x-2;x->2x-8]
    xs = [
         collect(linspace(-1,3.1,6)),
         collect(linspace(1,7,19)),
         collect(linspace(2,7,15)),
         collect(linspace(4,8,25))]
    ls = [Line(i[1],i[2](i[1])) for i in zip(xs,fs)]
    # create an envelope
    e = Envelope(ls)
    return e
end
function f3c()

    fs = Function[x->ones(length(x));x->0.5x;x->x-2;x->2x-8.5]
    xs = [
         collect(linspace(-1,3.1,6)),
         collect(linspace(1,7,19)),
         collect(linspace(2,7,15)),
         collect(linspace(4,8,25))]
    ls = [Line(i[1],i[2](i[1])) for i in zip(xs,fs)]
    # create an envelope
    e = Envelope(ls)
    return e
end

function tplot3a()

    en = f3a()

    p1 = plot(en)

    upper_env!(en)
    p2 = plot(en,removed=true)
    plot(p1,p2)

end

function tplot3b()

    en = f3b()

    p1 = plot(en)

    upper_env!(en)
    p2 = plot(en,removed=true)

    plot(p1,p2)

end
function tplot3c()

    en = f3c()

    p1 = plot(en)

    upper_env!(en)
    p2 = plot(en,removed=true)

    plot(p1,p2)

end

function tplot4()

    fs = Function[x->ones(length(x));x->0.5x;x->x-2;x->2x-8]
    xs = [
         collect(linspace(0.1,1.5,5)),
         collect(linspace(1,7,19)),
         collect(linspace(2,7,15)),
         collect(linspace(4,8,25))]
    ls = [Line(i[1],i[2](i[1])) for i in zip(xs,fs)]
    e = Envelope(ls)

    p1 = plot(e)

    upper_env!(e)

    p2 = plot(e,removed=true)

    plot(p1,p2)
end

function splitf()
    x = [1,2,3,1.5,2.1,2.9]
    y = [1,1.5,1.7,1.2,1.8,2.1]
    L = Line(x,y)
    p1 = plot(L,title="original",numerate=true)
    e = splitLine(L)
    p2 = plot(e,title="split Line",numerate=true)
    plot(p1,p2)
end

max_stars_repo_name: floswald/MultiLine

max_issues_repo_name: floswald/MultiLine

max_forks_repo_name: floswald/MultiLine