text:
using ProgressMeter
using DataFrames
using Statistics
using StatsBase
using Random

## Full EM. Probably will be needed when scRNA-seq is used as prior

function maximize_molecule_clusters!(cell_type_exprs::Matrix{Float64}, genes::Vector{Int}, confidence::Vector{Float64}, assignment_probs::Matrix{Float64};
        prior_exprs::Union{Matrix{Float64}, Nothing}=nothing, prior_stds::Union{Matrix{Float64}, Nothing}=nothing, add_pseudocount::Bool=false)
    cell_type_exprs .= 0.0;
    for i in 1:length(genes)
        t_gene = genes[i];
        t_conf = confidence[i];

        for j in 1:size(cell_type_exprs, 1)
            cell_type_exprs[j, t_gene] += t_conf * assignment_probs[j, i];
        end
    end

    if prior_exprs !== nothing
        mult = sum(cell_type_exprs, dims=2)
        cell_type_exprs .= adj_value_norm.(cell_type_exprs, prior_exprs .* mult, prior_stds .* mult)
    end

    if add_pseudocount
        cell_type_exprs .= (cell_type_exprs .+ 1) ./ (sum(cell_type_exprs, dims=2) .+ 1);
    else
        cell_type_exprs ./= sum(cell_type_exprs, dims=2);
    end
end

"""
Params:
- adjacent_weights: must be multiplied by confidence of the corresponding adjacent_point
"""
function expect_molecule_clusters!(assignment_probs::Matrix{Float64}, cell_type_exprs::Matrix{Float64}, genes::Vector{Int},
        adjacent_points::Vector{Vector{Int}}, adjacent_weights::Vector{Vector{Float64}})
    total_ll = 0.0
    for i in 1:length(genes)
        gene = genes[i]
        cur_weights = adjacent_weights[i]
        cur_points = adjacent_points[i]

        dense_sum = 0.0
        for ri in 1:size(assignment_probs, 1) # TODO: this function can probably be optimized by iterating by j first and then by ri
            c_d = 0.0
            for j in 1:length(cur_points) # TODO: can try to use sparsity to optimize it. Can store BitMatrix with info about a_p > 1e-10
                a_p = assignment_probs[ri, cur_points[j]]
                (a_p > 1e-5) || continue

                c_d += cur_weights[j] * a_p
            end

            assignment_probs[ri, i] = cell_type_exprs[ri, gene] * exp(c_d)
            dense_sum += assignment_probs[ri, i]
        end

        total_ll += log10(dense_sum)
        if dense_sum .> 1e-20
            assignment_probs[:, i] ./= dense_sum
        else
            assignment_probs[:, i] .= 1 / size(assignment_probs, 1)
        end
    end

    return total_ll
end

function cluster_molecules_on_mrf(df_spatial::DataFrame, adjacent_points::Vector{Vector{Int}}, adjacent_weights::Vector{Vector{Float64}};
        n_clusters::Int, confidence_threshold::Float64=0.95, kwargs...)

    cor_mat = pairwise_gene_spatial_cor(df_spatial.gene, df_spatial.confidence, adjacent_points, adjacent_weights; confidence_threshold=confidence_threshold);
    ct_exprs_init = nothing
    try
        ica_fit = fit(MultivariateStats.ICA, cor_mat, n_clusters, maxiter=10000);
        ct_exprs_init = copy((abs.(ica_fit.W) ./ sum(abs.(ica_fit.W), dims=1))')
    catch
        @warn "ICA did not converge, fall back to random initialization"
    end

    return cluster_molecules_on_mrf(df_spatial.gene, adjacent_points, adjacent_weights, df_spatial.confidence;
        cell_type_exprs=ct_exprs_init, n_clusters=n_clusters, kwargs...)
end

function cluster_molecules_on_mrf(genes::Vector{Int}, adjacent_points::Vector{Vector{Int}}, adjacent_weights::Vector{Vector{Float64}}, confidence::Vector{Float64};
        n_clusters::Int=1, tol::Float64=0.01, do_maximize::Bool=true, max_iters::Int=max(10000, div(length(genes), 200)), n_iters_without_update::Int=20,
        cell_type_exprs::Union{<:AbstractMatrix{Float64}, Nothing}=nothing, assignment::Union{Vector{Int}, Nothing}=nothing, assignment_probs::Union{Matrix{Float64}, Nothing}=nothing,
        verbose::Bool=true, progress::Union{Progress, Nothing}=nothing, weights_pre_adjusted::Bool=false, weight_mult::Float64=1.0, init_mod::Int=10000, kwargs...)
    # Initialization
    if cell_type_exprs === nothing
        if n_clusters <= 1
            if assignment !== nothing
                n_clusters = maximum(assignment)
            else
                error("Either n_clusters, assignment or cell_type_exprs must be specified")
            end
        end

        if init_mod < 0
            cell_type_exprs = copy(hcat(prob_array.(split(genes, rand(1:n_clusters, length(genes))), max_value=maximum(genes))...)')
        else
            gene_probs = prob_array(genes)
            cell_type_exprs = gene_probs' .* (0.95 .+ (hcat([[hash(x1 * x2^2) for x2 in 1:n_clusters] for x1 in 1:length(gene_probs)]...) .% init_mod) ./ 100000) # determenistic way of adding pseudo-random noise
        end
        # cell_type_exprs ./= sum(cell_type_exprs, dims=2)
        cell_type_exprs = (cell_type_exprs .+ 1) ./ (sum(cell_type_exprs, dims=2) .+ 1)
    else
        cell_type_exprs = deepcopy(Matrix(cell_type_exprs))
    end
    ## On the first iteration we don't have a pseudocount, but it shouldn't be a problem

    if !weights_pre_adjusted
        adjacent_weights = [weight_mult .* adjacent_weights[i] .* confidence[adjacent_points[i]] for i in 1:length(adjacent_weights)] # instead of multiplying each time in expect
    end

    if assignment_probs === nothing
        assignment_probs = zeros(size(cell_type_exprs, 1), length(genes));

        if assignment === nothing
            assignment_probs .= cell_type_exprs[:, genes]
            assignment_probs[:, vec(sum(assignment_probs, dims=1)) .< 1e-10] .= 1 / size(assignment_probs, 1)
            assignment_probs ./= sum(assignment_probs, dims=1)
        else
            assignment_probs[CartesianIndex.(assignment, 1:length(assignment))] .= 1.0;
        end
    else
        assignment_probs = deepcopy(assignment_probs)
    end

    assignment_probs_prev = deepcopy(assignment_probs)
    max_diffs = Float64[]
    change_fracs = Float64[]

    if verbose && progress === nothing
        progress = Progress(max_iters, 0.3)
    end

    # EM iterations
    n_iters = 0
    for i in 1:max_iters
        n_iters = i
        assignment_probs_prev .= assignment_probs

        expect_molecule_clusters!(assignment_probs, cell_type_exprs, genes, adjacent_points, adjacent_weights)

        if do_maximize
            maximize_molecule_clusters!(cell_type_exprs, genes, confidence, assignment_probs; add_pseudocount=true, kwargs...)
        end

        md, cf = estimate_difference_l0(assignment_probs, assignment_probs_prev, col_weights=confidence)
        push!(max_diffs, md)
        push!(change_fracs, cf)

        prog_vals = [("Iteration", i), ("Max. difference", md), ("Fraction of probs changed", cf)]
        if verbose
            next!(progress, showvalues=prog_vals)
        end

        if (i > n_iters_without_update) && (maximum(max_diffs[(end - n_iters_without_update):end]) < tol)
            if verbose
                finish!(progress, showvalues=prog_vals)
            end
            break
        end
    end

    if verbose
        @info "Algorithm stopped after $n_iters iterations. Error: $(round(max_diffs[end], sigdigits=3)). Converged: $(max_diffs[end] <= tol)."
    end

    if do_maximize
        maximize_molecule_clusters!(cell_type_exprs, genes, confidence, assignment_probs, add_pseudocount=false)
    end

    assignment = vec(mapslices(x -> findmax(x)[2], assignment_probs, dims=1));

    return (exprs=cell_type_exprs, assignment=assignment, diffs=max_diffs, assignment_probs=assignment_probs, change_fracs=change_fracs)
end

function filter_small_molecule_clusters(genes::Vector{Int}, confidence::Vector{Float64}, adjacent_points::Vector{Vector{Int}},
        assignment_probs::Matrix{Float64}, cell_type_exprs::Matrix{Float64}; min_mols_per_cell::Int, confidence_threshold::Float64=0.95)

    assignment = vec(mapslices(x -> findmax(x)[2], assignment_probs, dims=1));
    conn_comps_per_clust = get_connected_components_per_label(assignment, adjacent_points, 1;
        confidence=confidence, confidence_threshold=confidence_threshold)[1];
    n_mols_per_comp_per_clust = [length.(c) for c in conn_comps_per_clust];

    real_clust_ids = findall(maximum.(n_mols_per_comp_per_clust) .>= min_mols_per_cell)

    if length(real_clust_ids) == size(assignment_probs, 1)
        return assignment_probs, n_mols_per_comp_per_clust, real_clust_ids
    end

    assignment_probs = assignment_probs[real_clust_ids,:];
    cell_type_exprs = cell_type_exprs[real_clust_ids,:];

    for i in findall(vec(sum(assignment_probs, dims=1) .< 1e-10))
        assignment_probs[:, i] .= cell_type_exprs[:, genes[i]]
    end

    assignment_probs ./= sum(assignment_probs, dims=1);

    return assignment_probs, n_mols_per_comp_per_clust[real_clust_ids], real_clust_ids
end

## Utils

"""
    Adjust value based on prior. Doesn't penalize values < σ, penalize linearly values in [σ; 3σ], and super-linarly all >= 3σ
"""
@inline function adj_value_norm(x::Float64, μ::Float64, σ::Float64)::Float64
    dx = x - μ
    z = abs(dx) / σ
    if z < 1
        return x
    end

    if z < 3
        return μ + sign(dx) * (1 + (z - 1) / 4) * σ
    end

    return μ + sign(dx) * (sqrt(z) + 1.5 - sqrt(3)) * σ
end

function pairwise_gene_spatial_cor(genes::Vector{Int}, confidence::Vector{Float64}, adjacent_points::Array{Vector{Int}, 1}, adjacent_weights::Array{Vector{Float64}, 1};
        confidence_threshold::Float64=0.95)::Matrix{Float64}
    gene_cors = zeros(maximum(genes), maximum(genes))
    sum_weight_per_gene = zeros(maximum(genes))
    for gi in 1:length(genes)
        cur_adj_points = adjacent_points[gi]
        cur_adj_weights = adjacent_weights[gi]
        g2 = genes[gi]
        if confidence[gi] < confidence_threshold
            continue
        end

        for ai in 1:length(cur_adj_points)
            if confidence[cur_adj_points[ai]] < confidence_threshold
                continue
            end

            g1 = genes[cur_adj_points[ai]]
            cw = cur_adj_weights[ai]
            gene_cors[g2, g1] += cw
            sum_weight_per_gene[g1] += cw
            sum_weight_per_gene[g2] += cw
        end
    end

    for ci in 1:length(sum_weight_per_gene)
        for ri in 1:length(sum_weight_per_gene)
            gene_cors[ri, ci] /= fmax(sqrt(sum_weight_per_gene[ri] * sum_weight_per_gene[ci]), 0.1)
        end
    end

    return gene_cors
end

max_stars_repo_name: repushko/Baysor

max_issues_repo_name: repushko/Baysor

max_forks_repo_name: repushko/Baysor


text:
#=
exp_mult.jl
translated from matlab/C file mri_exp_mult_mex.c
originally by Sangwoo Lee, University of Michigan, Jun, 2004
based on the 2003 file int_tim_seg.m by Brad Sutton
2004-6-21 modified by JF
2020-02 by Daniel Wan
=#

export exp_mult

#using MIRT: exp_xform
using Test:@test

"""
    D = exp_mult(A, u, v ; warnboth)

Memory efficient and fast implementation of `D = A' * exp(-u * v^T)`
that is useful for B0-field-corrected MRI image reconstruction.

in:
* `A [N L]` matrix
* `u [N]` vector
* `v [M]` vector
* `warnboth` warn if both `u` and `v` are complex; default: true

out:
* `D [L M]` complex vector: `D = A' * exp(-u * v^T)`
`D_lm = sum_n A_nl^* exp(-u_n v_m)`
"""
function exp_mult(A, u::AbstractVector{<:Number}, v::AbstractVector{<:Number}
        ; warnboth::Bool = true)

    N, L = size(A)
    n = size(u,1)
    M = size(v,1)

    ~isreal(u) && ~isreal(v) && warnboth &&
        @warn "only one of u and v should be complex"

    T = promote_type(eltype(u), eltype(v), ComplexF32)
    D = zeros(T, L, M)

    for m = 1:M
        col = exp.(-u * v[m]) # [N] mth column of N × M matrix B = exp(-u * v^T)
        D[:,m] = A' * col
    end
    return D
end

"""
    exp_mult(:test)
self test
"""
function exp_mult(test::Symbol)
    test != :test && throw("bad test $test")

    L = 10
    N = 30
    M = 20
    A = randn(ComplexF32, N, L)
    ur = randn(N)
    ui = randn(N)
    vr = randn(M)
    vi = randn(M)
    u = ur + im * ui
    v = vr + im * vi

    d1 = A' * exp.(-ur * transpose(vr))
    @test d1 ≈ exp_mult(A, ur, vr)

    d1 = A' * exp.(-ur * transpose(v))
    @test d1 ≈ exp_mult(A, ur, v)

    d1 = A' * exp.(-u * transpose(vr))
    @test d1 ≈ exp_mult(A, u, vr)

    d1 = A' * exp.(-u * transpose(v))
    @test d1 ≈ exp_mult(A, u, v ; warnboth=false)

    d1 = A' * exp.(-ur * transpose(v))
    ur = reshape(ur, 1, :)
    v = reshape(v, 1, :)
    d2 = exp_xform(transpose(A'), ur, v, mode = :matrix)
    d2 = transpose(d2)
    @test d1 ≈ d2

    true
end

ext: jl


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie, Random
CairoMakie.activate!() # HIDE
let
    Random.seed!(123)
    n = 50
    x, y, color = rand(n), rand(n), rand(n)
    cmaps =[:cool,:viridis,:plasma,:inferno,:thermal,:leonardo,:winter,:spring,:ice]
    markers = [:+,:diamond,:star4,:rtriangle,:rect,:circle,:pentagon,:cross,:star5]

    function FigGridScatters()
        fig = Figure(resolution = (1200, 800))
        c = 1
        for i in 1:2, j in 1:2:5
            ax = Axis(fig[i, j],aspect = 1,xgridstyle=:dash,ygridstyle=:dash,
                                xtickalign=1, ytickalign=1)
            pnts = scatter!(x, y.^c, color = color, colormap=cmaps[c],
                markersize = 15, marker = markers[c], strokewidth=0)
            limits!(ax, -0.1,1.1,-0.1,1.1)
            ax.xticks = [0,1]
            ax.yticks = [0,1]
            ax.xticklabelsize = 20
            ax.yticklabelsize = 20
            cbar = Colorbar(fig, pnts,height = Relative(0.75), tickwidth = 2,
                            tickalign=1, width = 14, ticksize=14)
            fig[i, j+1] = cbar
            c+=1
        end
        fig
    end
    fig = FigGridScatters()
    #save("FigLineScatter.png", fig, px_per_unit = 2)
    save(joinpath(@__DIR__, "output", "ScattersCmaps.png"), fig, px_per_unit = 2) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie","Random"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
abstract type CompilerHint end
abstract type ProgramStructureHint <: CompilerHint end
abstract type AddressingHint <: CompilerHint end

include("static/kernel_hint.jl")
include("static/switch_hint.jl")
include("static/dynamic_address_hint.jl")


text:
# # Tutorial
#md #
#md # This tutorial is available in the Jupyter notebook format, togeter with other example notebooks, in the doc folder.
#md # To open Jupyter notebook in the correct folder simply type:
#md # ```julia
#md # using IJulia, Interact
#md # notebook(dir = Interact.notebookdir)
#md # ```
#md # in your Julia REPL. You can also view it online [here](https://github.com/JuliaGizmos/Interact.jl/blob/master/doc/notebooks/tutorial.ipynb).
#
# ## Installing everything
#
# To install Interact, simply type
# ```julia
# Pkg.add("Interact")
# ```

# in the REPL.
#
# The basic behavior is as follows: Interact provides a series of widgets. Each widget has an output that can be directly inspected or used to trigger some callbacks (i.e. run some code as soon as the widget changes value): the abstract supertype that gives this behavior is called `AbstractObservable`. Let's see this in practice.
#
# ## Displaying a widget
using Interact
ui = button()
display(ui)

# Note that `display` works in a [Jupyter notebook](https://github.com/JuliaLang/IJulia.jl) or in [Atom/Juno IDE](https://github.com/JunoLab/Juno.jl).
# Interact can also be deployed in Jupyter Lab, but that requires installing an extension first:
cd(Pkg.dir("WebIO", "assets"))
;jupyter labextension install webio
;jupyter labextension enable webio/jupyterlab_entry
# To deploy the app as a standalone Electron window, one would use [Blink.jl](https://github.com/JunoLab/Blink.jl):
using Blink
w = Window()
body!(w, ui);
# The app can also be served in a webpage via [Mux.jl](https://github.com/JuliaWeb/Mux.jl):
using Mux
WebIO.webio_serve(page("/", req -> ui), rand(8000:9000)) # serve on a random port
#
# ## Adding behavior
# The value of our button can be inspected using `getindex`:
ui[]
# In the case of a button, the observable represents the number of times it has been clicked: click on it and check the value again.
# For now however this button doesn't do anything. This can be changed by adding callbacks to it.
#
# To add some behavior to the widget we can use the `on` construct. `on` takes two arguments, a function and an `AbstractObservable`. As soon as the observable is changed, the function is called with the latest value.
on(println, ui)
# If you click again on the button you will see it printing the number of times it has been clicked so far.
#
# *Tip*: anonymous function are very useful in this programming paradigm. For example, if you want the button to say "Hello!" when pressed, you should use:
on(n -> println("Hello!"), ui)
#
# *Tip n. 2*: using the `[]` syntax you can also set the value of the widget:
ui[] = 33;

# ### Observables: the implementation of a widget's output
# The updatable container that only has the output of the widget but not the widget itself is a `Observable` and can be accessede using `observe(ui)`, though it should normally not be necessary to do so.
# To learn more about `Observables` and `AbstractObservable`, check out their documentation [here](https://juliagizmos.github.io/Observables.jl/latest/).

# ## What widgets are there?
#
# Once you have grasped this paradigm, you can play with any of the many widgets available:
filepicker() |> display # value is the path of selected file
textbox("Write here") |> display # value is the text typed in by the user
autocomplete(["Mary", "Jane", "Jack"]) |> display # as above, but you can autocomplete words
checkbox(label = "Check me!") |> display # value is a boolean describing whether it's ticked
toggle(label = "I have read and agreed") |> display # same as a checkbox but styled differently
slider(1:100, label = "To what extent?", value = 33) |> display # value is the number selected

# As well as the option widgets, that allow to choose among options:

dropdown(["a", "b", "c"]) |> display # value is option selected
togglebuttons(["a", "b", "c"]) |> display # value is option selected
radiobuttons(["a", "b", "c"]) |> display # value is option selected

# The option widgets can also take as input a dictionary (ordered dictionary is preferrable, to avoid items getting scrambled), in which case the label displays the key while the output stores the value:

s = dropdown(OrderedDict("a" => "Value 1", "b" => "Value 2"))
display(s)
#-
s[]
#
# ## Creating custom widgets
#
# Interact allows the creation of custom composite widgets starting from simpler ones.
# Let's say for example that we want to create a widget that has three sliders and a color
# that is updated to match the RGB value we gave with the sliders.

import Colors
using Plots

function mycolorpicker()
    r = slider(0:255, label = "red")
    g = slider(0:255, label = "green")
    b = slider(0:255, label = "blue")
    output = Interact.@map Colors.RGB(&r/255, &g/255, &b/255)
    plt = Interact.@map plot(sin, color = &output)
    wdg = Widget(["r" => r, "g" => g, "b" => b], output = output)
    @layout! wdg hbox(plt, vbox(:r, :g, :b)) ## custom layout: by default things are stacked vertically
end

# And now you can simply instantiate the widget with
mycolorpicker()
# Note the `&r` syntax: it means automatically update the widget as soon as the
# slider changes value. See [`Interact.@map`](@ref) for more details.
# If instead we wanted to only update the plot when a button is pressed we would do:
function mycolorpicker()
    r = slider(0:255, label = "red")
    g = slider(0:255, label = "green")
    b = slider(0:255, label = "blue")
    update = button("Update plot")
    output = Interact.@map (&update; Colors.RGB(r[]/255, g[]/255, b[]/255))
    plt = Interact.@map plot(sin, color = &output)
    wdg = Widget(["r" => r, "g" => g, "b" => b, "update" => update], output = output)
    @layout! wdg hbox(plt, vbox(:r, :g, :b, :update)) ## custom layout: by default things are stacked vertically
end

# ## A simpler approach for simpler cases
#
# While the approach sketched above works for all sorts of situations, there is a specific marcro to simplify it in some specific case. If you just want to update some result (maybe a plot) as a function of some parameters (discrete or continuous) simply write `@manipulate` before the `for` loop. Discrete parameters will be replaced by `togglebuttons` and continuous parameters by `slider`: the result will be updated as soon as you click on a button or move the slider:
#
width, height = 700, 300
colors = ["black", "gray", "silver", "maroon", "red", "olive", "yellow", "green", "lime", "teal", "aqua", "navy", "blue", "purple", "fuchsia"]
color(i) = colors[i%length(colors)+1]
ui = @manipulate for nsamples in 1:200,
        sample_step in slider(0.01:0.01:1.0, value=0.1, label="sample step"),
        phase in slider(0:0.1:2pi, value=0.0, label="phase"),
        radii in 0.1:0.1:60
    cxs_unscaled = [i*sample_step + phase for i in 1:nsamples]
    cys = sin.(cxs_unscaled) .* height/3 .+ height/2
    cxs = cxs_unscaled .* width/4pi
    dom"svg:svg[width=$width, height=$height]"(
        (dom"svg:circle[cx=$(cxs[i]), cy=$(cys[i]), r=$radii, fill=$(color(i))]"()
            for i in 1:nsamples)...
    )
end
# or, if you want a plot with some variables taking discrete values:
using Plots

x = y = 0:0.1:30

freqs = OrderedDict(zip(["pi/4", "π/2", "3π/4", "π"], [π/4, π/2, 3π/4, π]))

mp = @manipulate for freq1 in freqs, freq2 in slider(0.01:0.1:4π; label="freq2")
    y = @. sin(freq1*x) * sin(freq2*x)
    plot(x, y)
end

# ## Widget layout
#
# To create a full blown web-app, you should learn the layout tools that the CSS framework you are using provides. See for example the [columns](https://bulma.io/documentation/columns/) and [layout](https://bulma.io/documentation/layout/) section of the Bulma docs. You can use [WebIO](https://github.com/JuliaGizmos/WebIO.jl) to create from Julia the HTML required to create these layouts.
#
# However, this can be overwhelming at first (especially for users with no prior experience in web design). A simpler solution is [CSSUtil](https://github.com/JuliaGizmos/CSSUtil.jl), a package that provides some tools to create simple layouts.
loadbutton = filepicker()
hellobutton = button("Hello!")
goodbyebutton = button("Good bye!")
ui = vbox( # put things one on top of the other
    loadbutton,
    hbox( # put things one next to the other
        pad(1em, hellobutton), # to allow some white space around the widget
        pad(1em, goodbyebutton),
    )
)
display(ui)
#
# ## Update widgets as function of other widgets
#
# Sometimes the full structure of the GUI is not known in advance. For example, let's imagine we want to load a DataFrame and create a button per column. Not to make it completely trivial, as soon as a button is pressed, we want to plot a histogram of the corresponding column.
#
# *Important note*: this app needs to run in Blink, as the browser doesn't allow us to get access to the local path of a file.
#
# We start by adding a `filepicker` to choose the file, and only once we have a file we want to update the GUI. this can be done as follows:
loadbutton = filepicker()
columnbuttons = Observable{Any}(dom"div"())
# `columnbuttons` is the `div` object that will contain all the relevant buttons. it is an `Observable` as we want its value to change over time.
# To add behavior, we can use `map!`:
using CSV, DataFrames
data = Observable{Any}(DataFrame)
map!(CSV.read, data, loadbutton)
#
# Now as soon as a file is uploaded, the `Observable` `data` gets updated with the correct value. Now, as soon as `data` is updated, we want to update our buttons.
function makebuttons(df)
    buttons = button.(names(df))
    dom"div"(hbox(buttons))
end

map!(makebuttons, columnbuttons, data)
# We are almost done, we only need to add a callback to the buttons. The cleanest way is to do it during button initialization, meaning during our `makebuttons` step:
using Plots
plt = Observable{Any}(plot()) # the container for our plot
function makebuttons(df)
    buttons = button.(string.(names(df)))
    for (btn, name) in zip(buttons, names(df))
        map!(t -> histogram(df[name]), plt, btn)
    end
    dom"div"(hbox(buttons))
end
#
# To put it all together:
using CSV, DataFrames, Interact, Plots
loadbutton = filepicker()
columnbuttons = Observable{Any}(dom"div"())
data = Observable{Any}(DataFrame)
plt = Observable{Any}(plot())
map!(CSV.read, data, loadbutton)

function makebuttons(df)
    buttons = button.(string.(names(df)))
    for (btn, name) in zip(buttons, names(df))
        map!(t -> histogram(df[name]), plt, btn)
    end
    dom"div"(hbox(buttons))
end

map!(makebuttons, columnbuttons, data)

ui = dom"div"(loadbutton, columnbuttons, plt)
#
# And now to serve it in Blink:
using Blink
w = Window()
body!(w, ui)


text:
# # Advection-diffusion of tracer by cellular flow
#
#md # This example can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/literated/cellularflow.ipynb).
#
# An example demonstrating the advection-diffusion of a tracer by a cellular flow.
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add PassiveTracerFlows, Plots, Printf"
# ```

# ## Let's begin
# Let's load `PassiveTracerFlows.jl` and some other needed packages.
#
using PassiveTracerFlows, Plots, Printf

# ## Choosing a device: CPU or GPU

dev = CPU()     # Device (CPU/GPU)
nothing # hide

# ## Numerical parameters and time-stepping parameters

      n = 128            # 2D resolution = n²
stepper = "RK4"          # timestepper
     dt = 0.02           # timestep
 nsteps = 800            # total number of time-steps
 nsubs  = 25             # number of time-steps for intermediate logging/plotting (nsteps must be multiple of nsubs)
nothing # hide

# ## Numerical parameters and time-stepping parameters

L = 2π        # domain size
κ = 0.002     # diffusivity
nothing # hide

# ## Set up cellular flow
# We create a two-dimensional grid to construct the cellular flow. Our cellular flow is derived
# from a streamfunction ``ψ(x, y) = ψ₀ \cos(x) \cos(y)`` as ``(u, v) = (-∂_y ψ, ∂_x ψ)``.
grid = TwoDGrid(n, L)

ψ₀ = 0.2
mx, my = 1, 1

ψ = [ψ₀ * cos(mx * grid.x[i]) * cos(my * grid.y[j]) for i in 1:grid.nx, j in 1:grid.ny]

uvel(x, y) =  ψ₀ * my * cos(mx * x) * sin(my * y)
vvel(x, y) = -ψ₀ * mx * sin(mx * x) * cos(my * y)
nothing # hide

# ## Problem setup
# We initialize a `Problem` by providing a set of keyword arguments.
prob = TracerAdvectionDiffusion.Problem(dev; nx=n, Lx=L, κ=κ, steadyflow=true, u=uvel, v=vvel,
                                          dt=dt, stepper=stepper)
nothing # hide

# and define some shortcuts
sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
x, y = grid.x, grid.y
nothing # hide

# ## Setting initial conditions

# Our initial condition for the tracer ``c`` is a gaussian centered at ``(x, y) = (L_x/5, 0)``.

gaussian(x, y, σ) = exp(-(x^2 + y^2) / (2σ^2))

amplitude, spread = 0.5, 0.15
c₀ = [amplitude * gaussian(x[i] - 0.2 * grid.Lx, y[j], spread) for i=1:grid.nx, j=1:grid.ny]

TracerAdvectionDiffusion.set_c!(prob, c₀)
nothing # hide

# Let's plot the initial tracer concentration and streamlines. Note that when plotting, we decorate
# the variable to be plotted with `Array()` to make sure it is brought back on the CPU when
# `vars` live on the GPU.

function plot_output(prob)
  c = prob.vars.c

  p = heatmap(x, y, Array(vars.c'),
         aspectratio = 1,
              c = :balance,
         legend = :false,
           clim = (-0.2, 0.2),
          xlims = (-grid.Lx/2, grid.Lx/2),
          ylims = (-grid.Ly/2, grid.Ly/2),
         xticks = -3:3,
         yticks = -3:3,
         xlabel = "x",
         ylabel = "y",
          title = "initial tracer concentration (shading) + streamlines",
     framestyle = :box)

  contour!(p, x, y, Array(ψ'),
     levels=0.0125:0.025:0.2,
     lw=2, c=:black, ls=:solid, alpha=0.7)

  contour!(p, x, y, Array(ψ'),
     levels=-0.1875:0.025:-0.0125,
     lw=2, c=:black, ls=:dash, alpha=0.7)

  return p
end
nothing # hide

# ## Time-stepping the `Problem` forward

# We time-step the `Problem` forward in time.

startwalltime = time()

p = plot_output(prob)

anim = @animate for j = 0:round(Int, nsteps/nsubs)

 if j % (200 / nsubs) == 0
    log = @sprintf("step: %04d, t: %d, walltime: %.2f min",
                   clock.step, clock.t, (time()-startwalltime)/60)

    println(log)
  end

  p[1][1][:z] = Array(vars.c)
  p[1][:title] = "concentration, t=" * @sprintf("%.2f", clock.t)

  stepforward!(prob, nsubs)
  TracerAdvectionDiffusion.updatevars!(prob)
end

mp4(anim, "cellularflow.mp4", fps=12)

max_stars_repo_name: FourierFlows/PassiveTracerFlows

max_issues_repo_name: FourierFlows/PassiveTracerFlows

max_forks_repo_name: FourierFlows/PassiveTracerFlows


text:
# Generic Spiking Layer driving an HH layer

# verifies that identical stimulus presented to 2 disconnected cells produces exactly the same datas
tic();

using Base.Test
using NeuralNets

N_stim = 100 # input dim
T = 300
dt = NeuralNets.dt

net = buildnet()

spike_train = image_stim("../images/021108_spasticroomjpg_1_.jpg", int(sqrt(N_stim)), T, dt)
stim_spikes = cell(length(net.l))
for li=1:length(net.l)
	if li in 1:3
		stim_spikes[li] = spike_train[:,:,li]
	else
		stim_spikes[li] = nothing
	end
end
rs = RunSettings(T,stim_spikes)

out = runSpikingNet!(net,rs)
# save out the data
using HDF5, JLD
save("trial1.jld","results",out)

toc()

# PROFILING
# rs.T = dt # single iteration to compile it
# println("Initial Compilation using JIT...")
# out = runSpikingNet!(net,rs)
# println("Profiling Network Simulation...")
# rs.T = T
# Profile.clear()
# @profile out = runSpikingNet!(net,rs)
# using ProfileView
# #Profile.init(1_000_000, 1)
# ProfileView.view()


text:
using Documenter, Animations

makedocs(
    sitename="Animations.jl",
    format = Documenter.HTML(
            prettyurls = get(ENV, "CI", nothing) == "true"
        )
    )

deploydocs(
    repo = "github.com/jkrumbiegel/Animations.jl.git",
)


text:
function decoder!(item::Int64, buckets::Array{Int64,2}, keys::Array{Float64,2})
    @inbounds for i in 1:size(buckets, 1)
        if keys[i,item] < 0.0001
            buckets[i,item] = 1
        else
            buckets[i,item] = Int64(ceil(keys[i,item] * size(buckets, 1)))
        end
    end
end

function objective!(item::Int64, buckets::Array{Int64,2}, fitness::Array{Int64,1}, instance::OBOPDataset)
    fitness[item] = 0
    @inbounds for i in 1:size(instance.C, 1) - 1
        @inbounds for j in i + 1:size(instance.C, 1)
            if buckets[i,item] == buckets[j,item]
                fitness[item] += instance.C[i,j]
                fitness[item] += instance.C[j,i]
            else
                if buckets[i,item] < buckets[j,item]
                    fitness[item] += instance.C[i,j]
                    fitness[item] -= instance.C[j,i]
                else
                    fitness[item] -= instance.C[i,j]
                    fitness[item] += instance.C[j,i]
                end
            end
        end
    end
end

max_stars_repo_name: LuizHNLorena/AbrkgaObop

max_issues_repo_name: LuizHNLorena/AbrkgaObop

max_forks_repo_name: LuizHNLorena/AbrkgaObop


text:
function check_derivative(f::Function, g::Function, x::Number)
	auto_g = derivative(f)
	return maximum(abs(g(x) - auto_g(x)))
end

function check_gradient{T <: Number}(f::Function, g::Function, x::Vector{T})
	auto_g = gradient(f)
	return maximum(abs(g(x) - auto_g(x)))
end

function check_second_derivative(f::Function, h::Function, x::Number)
	auto_h = second_derivative(f)
	return maximum(abs(h(x) - auto_h(x)))
end

function check_hessian{T <: Number}(f::Function, h::Function, x::Vector{T})
	auto_h = hessian(f)
	return maximum(abs(h(x) - auto_h(x)))
end


text:
function MLJBase.fit(model::EvoTypes, verbosity::Int, X, y)
    Xmatrix = MLJBase.matrix(X)
    fitresult, cache = init_evotree(model, Xmatrix, y, verbosity = verbosity)
    grow_evotree!(fitresult, cache, verbosity = verbosity)
    report = nothing
    return fitresult, cache, report
end

function okay_to_continue(new, old)
    new.nrounds - old.nrounds >= 0 &&
    new.loss == old.loss &&
    new.λ == old.λ &&
    new.γ == old.γ &&
    new.max_depth  == old.max_depth &&
    new.min_weight == old.min_weight &&
    new.rowsample ==  old.rowsample &&
    new.colsample ==  old.colsample &&
    new.nbins ==  old.nbins &&
    new.α ==  old.α &&
    new.metric ==  old.metric
end

function MLJBase.update(model::EvoTypes, verbosity,
    fitresult, cache, X, y)

    if okay_to_continue(model, cache.params)
        grow_evotree!(fitresult, cache, verbosity = verbosity)
    else
        Xmatrix = MLJBase.matrix(X)
        fitresult, cache = init_evotree(model, Xmatrix, y, verbosity = verbosity)
        grow_evotree!(fitresult, cache, verbosity = verbosity)
    end
    report = nothing
    return fitresult, cache, report
end

function predict(model::EvoTypes, fitresult, Xnew)
    Xnew = MLJBase.matrix(Xnew)
    pred = predict(fitresult, Xnew)
    return pred
end

function predict(model::EvoTreeClassifier, fitresult, Xnew)
    Xnew = MLJBase.matrix(Xnew)
    pred = predict(fitresult, Xnew)
    return [MLJBase.UnivariateFinite(fitresult.levels, pred[i,:]) for i in 1:size(pred,1)]
end

function predict(model::EvoTreeCount, fitresult, Xnew)
    Xnew = MLJBase.matrix(Xnew)
    λ = predict(fitresult, Xnew)
    return [Distributions.Poisson(λᵢ) for λᵢ ∈ λ]
end

function predict(model::EvoTreeGaussian, fitresult, Xnew)
    Xnew = MLJBase.matrix(Xnew)
    pred = predict(fitresult, Xnew)
    return [Distributions.Normal(pred[i]...) for i in 1:size(pred,1)]
end

# MLJBase.predict_mean(model::Union{EvoTreeRegressor, EvoTreeCount, EvoTreeGaussian}, fitresult, Xnew) =
#     mean.(MLJBase.predict(model, fitresult, Xnew))
#
# MLJBase.predict_mode(model::Union{EvoTreeRegressor, EvoTreeClassifier,EvoTreeCount, EvoTreeGaussian}, fitresult, Xnew) =
#     mode.(MLJBase.predict(model, fitresult, Xnew))
#
# MLJBase.predict_median(model::Union{EvoTreeRegressor, EvoTreeCount, EvoTreeGaussian}, fitresult, Xnew) =
#     median.(MLJBase.predict(model, fitresult, Xnew))

# shared metadata
MLJBase.package_name(::Type{<:EvoTypes}) = "EvoTrees"
MLJBase.package_uuid(::Type{<:EvoTypes}) = "f6006082-12f8-11e9-0c9c-0d5d367ab1e5"
MLJBase.package_url(::Type{<:EvoTypes}) = "https://github.com/Evovest/EvoTrees.jl"
MLJBase.is_pure_julia(::Type{<:EvoTypes}) = true

MLJBase.load_path(::Type{<:EvoTreeRegressor}) = "EvoTrees.EvoTreeRegressor"
MLJBase.input_scitype(::Type{<:EvoTreeRegressor}) = MLJBase.Table(MLJBase.Continuous)
MLJBase.target_scitype(::Type{<:EvoTreeRegressor}) = AbstractVector{<:MLJBase.Continuous}

MLJBase.load_path(::Type{<:EvoTreeCount}) = "EvoTrees.EvoTreeCount"
MLJBase.input_scitype(::Type{<:EvoTreeCount}) = MLJBase.Table(MLJBase.Continuous)
MLJBase.target_scitype(::Type{<:EvoTreeCount}) = AbstractVector{<:MLJBase.Count}

MLJBase.load_path(::Type{<:EvoTreeClassifier}) = "EvoTrees.EvoTreeClassifier"
MLJBase.input_scitype(::Type{<:EvoTreeClassifier}) = MLJBase.Table(MLJBase.Continuous)
MLJBase.target_scitype(::Type{<:EvoTreeClassifier}) = AbstractVector{<:MLJBase.Finite}

MLJBase.load_path(::Type{<:EvoTreeGaussian}) = "EvoTrees.EvoTreeGaussian"
MLJBase.input_scitype(::Type{<:EvoTreeGaussian}) = MLJBase.Table(MLJBase.Continuous)
MLJBase.target_scitype(::Type{<:EvoTreeGaussian}) = AbstractVector{<:MLJBase.Continuous}

# function MLJ.clean!(model::EvoTreeRegressor)
#     warning = ""
#     if model.nrounds < 1
#         warning *= "Need nrounds ≥ 1. Resetting nrounds=1. "
#         model.nrounds = 1
#     end
#     if model.λ < 0
#         warning *= "Need λ ≥ 0. Resetting λ=0. "
#         model.λ = 0.0
#     end
#     if model.γ < 0
#         warning *= "Need γ ≥ 0. Resetting γ=0. "
#         model.γ = 0.0
#     end
#     if model.η <= 0
#         warning *= "Need η > 0. Resetting η=0.001. "
#         model.η = 0.001
#     end
#     if model.max_depth < 1
#         warning *= "Need max_depth ≥ 0. Resetting max_depth=0. "
#         model.max_depth = 1
#     end
#     if model.min_weight < 0
#         warning *= "Need min_weight ≥ 0. Resetting min_weight=0. "
#         model.min_weight = 0.0
#     end
#     if model.rowsample < 0
#         warning *= "Need rowsample ≥ 0. Resetting rowsample=0. "
#         model.rowsample = 0.0
#     end
#     if model.rowsample > 1
#         warning *= "Need rowsample <= 1. Resetting rowsample=1. "
#         model.rowsample = 1.0
#     end
#     if model.colsample < 0
#         warning *= "Need colsample ≥ 0. Resetting colsample=0. "
#         model.colsample = 0.0
#     end
#     if model.colsample > 1
#         warning *= "Need colsample <= 1. Resetting colsample=1. "
#         model.colsample = 1.0
#     end
#     if model.nbins > 250
#         warning *= "Need nbins <= 250. Resetting nbins=250. "
#         model.nbins = 250
#     end
#     return warning
# end


text:
# ---
# title: 143. Reorder List
# id: problem143
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Linked List
# link: <https://leetcode.com/problems/reorder-list/description/>
# hidden: true
# ---
#
# Given a singly linked list _L_ : _L_ 0-> _L_ 1->…-> _L_ _n_ -1-> _L_ n,
# reorder it to: _L_ 0-> _L_ _n_ -> _L_ 1-> _L_ _n_ -1-> _L_ 2-> _L_ _n_ -2->…
#
# You may **not** modify the values in the list's nodes, only nodes itself may
# be changed.
#
# **Example 1:**
#
#
#
#     Given 1->2->3->4, reorder it to 1->4->2->3.
#
# **Example 2:**
#
#
#
#     Given 1->2->3->4->5, reorder it to 1->5->2->4->3.
#
#
#
## @lc code=start
using LeetCode

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


text:
# module Bukdu

function render(::Type{HTML}, obj::Any)::Conn
    filtering(render,HTML,obj) do
        Conn(200, Dict("Content-Type"=>"text/html"), obj)
    end
end

include("html/tag.jl")


text:
# Methods and structs from this file are not meant to be part ot the public interface, so
# documentation is in the comments, no docstrings here.

###############################################################
# Check for LC_NUMERIC
# Possibly there are other options to handle this, e.g. replace strtod call in triangle
# See e.g.  https://github.com/JuliaLang/julia/issues/5928
# It might be the case that this problem occurs only in vendor distributions of Julia

function checklocale()
    endptr::Ptr{Ptr{UInt8}}=C_NULL
    if ccall((:strtod), Cdouble, (Cstring,Ptr{Ptr{UInt8}},), "0.5",endptr)!=0.5
        print(
            """

             Triangle uses  the C  library function strtod  to convert
             strings to double numbers.   The interpretation of double
             numbers  (whether  they  are  written with  ','  or  '.')
             depends on  the language  settings of your  computer.  In
             order  to ensure  portability  of  programs written  with
             Triangulate.jl, one should insist  on assuming '.' as the
             decimal point.

             In  the moment,  due to  whatever reason  (e.g. due  to a
             change  of  settings  performed  by  PyPlot),  the  wrong
             behavior has been detected.  The correct behaviour can be
             enforced by setting the environment variable 'LC_NUMERIC'
             to 'C' before starting Julia.

             See e.g. https://github.com/JuliaComputing/TextParse.jl/issues/30
         """)
        error("Missing or wrong value of LC_NUMERIC")
        return false
    end
    true
end

###############################################################
# Handling of triunsuitable callback

# Trivial default trinunsuitable function
function trivial_triunsuitable(org_x, org_y, dest_x, dest_y, apex_x, apex_y, area)
    return 0
end

# Global variable containing triunsuitable function
triunsuitable_func=trivial_triunsuitable

# triunsuitable function called from C (by triangulate(::CTriangulateIO)) if -u flag has been set
function jl_wrap_triunsuitable(org_x::Cdouble, org_y::Cdouble, dest_x::Cdouble, dest_y::Cdouble, apex_x::Cdouble, apex_y::Cdouble, area::Cdouble)::Cint
    return Cint(triunsuitable_func(org_x, org_y, dest_x, dest_y, apex_x, apex_y, area))
end

#
# Struct mapping Triangle's triangulateio
#
# We use a mutable struct because we need fields to be able  initialized
# in differerent ways.
#
mutable struct CTriangulateIO

    pointlist :: Ptr{Cdouble}
    pointattributelist :: Ptr{Cdouble}
    pointmarkerlist :: Ptr{Cint}
    numberofpoints :: Cint
    numberofpointattributes :: Cint

    trianglelist :: Ptr{Cint}
    triangleattributelist :: Ptr{Cdouble}
    trianglearealist :: Ptr{Cdouble}
    neighborlist :: Ptr{Cint}
    numberoftriangles :: Cint
    numberofcorners :: Cint
    numberoftriangleattributes :: Cint

    segmentlist :: Ptr{Cint}
    segmentmarkerlist :: Ptr{Cint}
    numberofsegments :: Cint

    holelist :: Ptr{Cdouble}
    numberofholes :: Cint

    regionlist :: Ptr{Cdouble}
    numberofregions :: Cint

    edgelist :: Ptr{Cint}
    edgemarkerlist :: Ptr{Cint}
    normlist :: Ptr{Cdouble}
    numberofedges :: Cint
end

#
# Constructor for `CTriangulateIO`. Initializes everything as `NULL`.
#
function CTriangulateIO()
    return CTriangulateIO(C_NULL, C_NULL, C_NULL, 0, 0,
                          C_NULL, C_NULL, C_NULL, C_NULL, 0, 0, 0,
                          C_NULL, C_NULL, 0,
                          C_NULL, 0,
                          C_NULL, 0,
                          C_NULL, C_NULL, C_NULL, 0)
end

#
# "Raw" triangulation call to Triangle library
#
function triangulate(triangle_switches::String,
                     ctio_in::CTriangulateIO,
                     ctio_out::CTriangulateIO,
                     vor_out::CTriangulateIO)

    # Check locale settings for decimal point
    checklocale()

    # Set unsuitable callback
    if occursin("u",triangle_switches)
        c_wrap_triunsuitable=@cfunction(jl_wrap_triunsuitable, Cint, (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,))
        ccall((:triunsuitable_callback,libtriangle),Cvoid,(Ptr{Cvoid},),c_wrap_triunsuitable)
    end

    # Call triangulate
    rc=ccall((:triangulate_catch_exit,libtriangle),
             Cint,
             ( Cstring,
               Ref{CTriangulateIO},
               Ref{CTriangulateIO},
               Ref{CTriangulateIO}
               ),
             triangle_switches,
             Ref(ctio_in),
             Ref(ctio_out),
             Ref(vor_out)
             )
    return rc
end


text:
using Markdown

mutable struct MapRayClassGrpNew{T} #<: Hecke.Map{T, Hecke.NfOrdIdlSet}
  header::Hecke.MapHeader
  modulus_fin::Hecke.NfOrdIdl
  modulus_inf::Vector{Hecke.InfPlc}

  function MapRayClassGrpNew{T}() where {T}
    return new{T}()
  end
end

#
# Modify the map of the class group so that the chosen representatives are coprime to m
#

function _coprime_ideal(C::GrpAbFinGen, mC::Map, m::NfOrdIdl)

  O=parent(m).order
  K=nf(O)
  L=NfOrdIdl[]
  for i=1:ngens(C)
    a=mC(C[i])
    if iscoprime(a,m)
      push!(L,a)
    else
      J=inv(a)
      s=K(rand(J.num,5))//J.den  # Is the bound acceptable?
      I=s*a
      simplify(I)
      I = numerator(I)
      while !iscoprime(I,m)
        s=K(rand(J.num,5))//J.den
        I=s*a
        simplify(I)
        I = numerator(I)
      end
      push!(L,I)
    end
  end

  function exp(a::GrpAbFinGenElem)
    I=ideal(O,1)
    for i=1:ngens(C)
      if Int(a.coeff[1,i])!= 0
        I=I*(L[i]^(Int(a.coeff[1,i])))
      end
    end
    return I
  end

  mp=Hecke.MapClassGrp{typeof(C)}()
  mp.header=Hecke.MapHeader(C,Hecke.NfOrdIdlSet(O),exp, mC.header.preimage)

  return mp

end

@doc Markdown.doc"""
***
    ray_class_group(m::NfOrdIdl, A::Vector{InfPlc}=[]) -> FinGenGrpAb, Map

> Compute the ray class group of the maximal order $L$ with respect to the modulus given by $m$ (the finite part) and the infinite primes of $A$
> and return an abstract group isomorphic to the ray class group with a map
> from the group to the ideals of $L$

"""
function ray_class_group_std(m::NfOrdIdl, primes::Vector{InfPlc}=InfPlc[])

#
# We compute the group using the sequence U -> (O/m)^* _> Cl^m -> Cl -> 1
# First of all, we compute all these groups with their own maps
#
  O=parent(m).order
  K=nf(O)

  C, mC = class_group(O)
  mC=_coprime_ideal(C,mC,m)
  U, mU = unit_group(O)
  M, pi= quo(O,m)
  G, mG=unit_group(M)

  p = [ x for x in primes if isreal(x) ]
  if !isempty(p)
    H,lH,eH=_infinite_primes(O,p,m)
    T=G
    G=direct_product(G,H)
  end

#
# We start to construct the relation matrix
#
  RG=rels(G)
  RC=rels(C)

  A=vcat(RC, MatrixSpace(FlintZZ, ngens(G)+ngens(U), ncols(RC))())
  B=vcat(MatrixSpace(FlintZZ, ngens(C), ncols(RG))(), RG)
  B=vcat(B, MatrixSpace(FlintZZ, ngens(U) , ncols(RG))())

#
# We compute the relation matrix given by the image of the map U -> (O/m)^*
#
  for i=1:ngens(U)
    u=mU(U[i])
    a=(mG\(pi(u))).coeff
    if !isempty(primes)
      a=hcat(a, (lH(K(u))).coeff)
    end
    for j=1:ngens(G)
      B[i+nrows(RC)+nrows(RG),j]=a[1,j]
    end
  end

#
# We compute the relation between generators of Cl and (O/m)^* in Cl^m
#

  for i=1: ngens(C)
    if order(C[i])!=1
      y=Hecke.principal_generator((mC(C[i]))^(Int(order(C[i]))))
      b=(mG\(pi(y))).coeff
      if primes != []
        b=hcat(b, (lH(K(y))).coeff)
      end
      for j=1: ngens(G)
        B[i,j]=-b[1,j]
      end
    end
  end

  R=hcat(A,B)
  X=abelian_group(R)

#
# Discrete logarithm
#

  function disclog(J::NfOrdIdl)

    !iscoprime(J,m) && error("The ideal is not coprime to the modulus")
    if isone(J)
      return X([0 for i=1:ngens(X)])
    else
      L=mC\J
      s=mC(L)
      I= J // s
      simplify(I)
      gamma=Hecke.principal_generator(I.num)
      y=((mG\(pi(gamma)))-(mG\(pi(O(I.den))))).coeff
      if primes!=[]
        z=lH(K(gamma))
        y=hcat(y, z.coeff)
      end
      return X(hcat(L.coeff,y))
    end
  end

#
# Exp map
#

  function expo(a::GrpAbFinGenElem)
    b=C([a.coeff[1,i] for i=1:ngens(C)])
    if isempty(primes)
      c=G([a.coeff[1,i] for i=ngens(C)+1:ngens(X)])
      return mC(b)*(pi\(mG(c)))
    else
      c=T([a.coeff[1,i] for i=ngens(C)+1:ngens(T)+ngens(C)])
      d=H([a.coeff[1,i] for i=ngens(T)+ngens(C)+1: ngens(X)])
      el=pi\(mG(c))
      # I need to modify $el$ so that it has the correct sign at the embeddings contained in primes
      vect=(lH(K(el))).coeff
      if vect==d.coeff
        return el*mC(b)
      else
        correction=O(1)
        for i=1:ngens(H)
          if d.coeff[1,i]==1
            correction=correction*eH(H[i])
          end
        end
        while vect!=d.coeff
          el=el+correction
          vect=(lH(K(el))).coeff
        end
        return el*mC(b)
      end
    end
  end

  mp=MapRayClassGrpNew{typeof(X)}()
  mp.header = Hecke.MapHeader(X, Hecke.NfOrdIdlSet(O), expo, disclog)
  mp.modulus_fin=m
  mp.modulus_inf=p

  return X, mp

end

########################################################
#
#  Ray Class Group - p-part
#
########################################################
function ray_class_group_p_part(p::Integer, m::NfOrdIdl, inf_plc::Vector{InfPlc}=InfPlc[])

  O=parent(m).order
  K=nf(O)
  Q,pi=quo(O,m)
  G, mG, lp = Hecke._mult_grp(Q,p)
  C, mC = class_group(O)
  Hecke._assure_princ_gen(mC)

  if p==2
    pr = [ x for x in inf_plc if isreal(x) ]
    if !isempty(pr)
      H,eH,lH=Hecke._infinite_primes(O,pr,m)
      T=G
      G =Hecke.direct_product(G,H)
    end
  end

  valclassp=Int(p^(valuation(order(C[ngens(C)]),p)))
  nonppartclass=div(order(C[ngens(C)]),valclassp)

  if valclassp==1 && order(G)==1
    return Hecke.empty_ray_class(m)
  end

  C, mC, vect = Hecke._class_group_mod_n(C,mC,valclassp)
  U, mU = unit_group_fac_elem(O)
  exp_class, Kel = Hecke._elements_to_coprime_ideal(C,mC,m)

  if order(G)==1
    return Hecke.class_as_ray_class(C,mC,exp_class,m)
  end

  n=ideal(O,1)
  for (q,vq) in lp
    n*=q^vq
  end

#
# We start to construct the relation matrix
#

  expo=exponent(G)
  inverse_d=gcdx(nonppartclass,expo)[2]

  R=zero_matrix(FlintZZ, ngens(C)+ngens(U)+ngens(G), ngens(C)+ngens(G))
  for i=1:ngens(C)
    R[i,i]=C.snf[i]
  end
  if issnf(G)
    for i=1:ngens(G)
      R[i+ngens(C),i+ngens(C)]=G.snf[i]
    end
  else
    for i=1:ngens(G)
      R[i+ngens(C),i+ngens(C)]=G.rels[i,i]
    end
  end

#
# We compute the relation matrix given by the image of the map U -> (O/m)^*
#

  @assert issnf(U)
  @vprint :RayFacElem 1 "Collecting elements to be evaluated; first, units \n"
  evals = Hecke.NfOrdQuoRingElem[]
  tobeeval = FacElem{nf_elem, AnticNumberField}[]
  if gcd(U.snf[1],p)!=1
    if U.snf[1]==2
      push!(evals,Q(-1))
    else
      push!(tobeeval, mU(U[1]))
    end
  else
    push!(evals, Q(1))
  end
  append!(tobeeval,[mU(U[i]) for i=2:ngens(U)])

  @vprint :RayFacElem 1 "then principal ideal generators \n"
  for i=1:ngens(C)
    @vtime :RayFacElem 1 push!(tobeeval, mC.princ_gens[i][2]*Kel[i])
  end

  @vprint :RayFacElem 1 "Time for elements evaluation: "
  ev,quots,idemps= Hecke.fac_elems_eval(O,Q,tobeeval,lp,expo)
  @vtime :RayFacElem 1 append!(evals, ev)
  @vprint :RayFacElem 1 "\n"

  for i=1:ngens(U)
    @vprint :RayFacElem 1 "Disclog of unit $i \n"
    a=(mG\(evals[i])).coeff
    if p==2 && !isempty(pr)
      if i==1
        a=hcat(a, matrix(FlintZZ,1,length(pr), [1 for i in pr]))
      else
        b=lH(tobeeval[length(tobeeval)-ngens(C)-ngens(U)+i])
        a=hcat(a, b.coeff)
      end
    end
    for j=1:ngens(G)
      R[i+ngens(G)+ngens(C),ngens(C)+j]=a[1,j]
    end
  end

#
# We compute the relation between generators of Cl and (O/m)^* in Cl^m
#

  for i=1: ngens(C)
    @vprint :RayFacElem 1 "Disclog of class group element $i \n"
    invn=gcdx(vect[i], C.snf[i])[2]
    a=((mG\(evals[i+ngens(U)]))*invn).coeff
    if p==2 && !isempty(pr)
      b=lH(tobeeval[length(tobeeval)-ngens(C)+i])
      a=hcat(a, b.coeff)
    end
    for j=1: ngens(G)
      R[i,ngens(C)+j]=-a[1,j]
    end
  end

  X=abelian_group(R)

#
# Discrete logarithm
#

  function disclog(J::FacElem)
    a= X([0 for i=1:ngens(X)])
    for (f,k) in J.fac
      a+=k*disclog(f)
    end
    return a
  end

  function disclog(J::NfOrdIdl)
    if isone(J)
      return X([0 for i=1:ngens(X)])
    else
      L=mC\J
      s=exp_class(L)
      I=J* inv(s)
      I=I^Int(nonppartclass)
      z=Hecke.principal_generator_fac_elem(I)
      el=Hecke._fac_elem_evaluation(O, Q, quots, idemps, z, lp, expo)
      y=((mG\pi(el))*inverse_d).coeff
      if p==2 && !isempty(pr)
        b=lH(z)
        y=hcat(y, b.coeff)
      end
      return X(hcat(L.coeff,y))
    end
  end

#
# Exp map
#

  function expon(a::GrpAbFinGenElem)
    b=C([a.coeff[1,i] for i=1:ngens(C)])
    if p!=2 || isempty(pr)
      c=G([a.coeff[1,i] for i=ngens(C)+1:ngens(X)])
      return exp_class(b)*ideal(O,pi\(mG(c)))
    else
      c=T([a.coeff[1,i] for i=ngens(C)+1:ngens(T)+ngens(C)])
      d=H([a.coeff[1,i] for i=ngens(T)+ngens(C)+1: ngens(X)])
      el=pi\(mG(c))
      # I need to modify $el$ so that it has the correct sign at the embeddings contained in primes
      vect=(lH(K(el))).coeff
      if vect==d.coeff
        return exp_class(b)*ideal(O,el)
      else
        correction=eH(d)
        while vect!=d.coeff
          el=el+correction
          vect=(lH(K(el))).coeff
        end
        return exp_class(b)*ideal(O,el)
      end
    end
  end

  mp=Hecke.MapRayClassGrpNew{typeof(X)}()
  mp.header = Hecke.MapHeader(X, FacElemMon(parent(m)) , expon, disclog)
  mp.modulus_fin=n
  mp.modulus_inf=inf_plc
  mp.fact_mod=lp
  mp.tame_mult_grp=mG.tame
  mp.wild_mult_grp=mG.wild
  mp.defining_modulus = (m, inf_plc)

  return X,mp
end


text:
using DashBootstrapComponents

badge = dbc_button(
    ["Notifications", dbc_badge("4", color="light", className="ml-1")],
    color="primary",
);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
using MetaGraphs: AbstractMetaGraph

## Convolution layers accepting AbstractMetaGraph

GCNConv(g::AbstractMetaGraph, ch::Pair{<:Integer,<:Integer}, σ=identity; kwargs...) =
    GCNConv(g.graph, ch, σ; kwargs...)

ChebConv(g::AbstractMetaGraph, ch::Pair{<:Integer,<:Integer}, k::Integer; kwargs...) =
    ChebConv(g.graph, ch, k; kwargs...)

GraphConv(g::AbstractMetaGraph, ch::Pair{<:Integer,<:Integer}, σ=identity, aggr=:add; kwargs...) =
    GraphConv(g.graph, ch, σ, aggr; kwargs...)

GATConv(g::AbstractMetaGraph, ch::Pair{<:Integer,<:Integer}; kwargs...) =
    GATConv(g.graph, ch; kwargs...)

GatedGraphConv(g::AbstractMetaGraph, out_ch::Integer, num_layers::Integer; kwargs...) =
    GatedGraphConv(g.graph, out_ch, num_layers; kwargs...)

EdgeConv(g::AbstractMetaGraph, nn; kwargs...) = EdgeConv(g.graph, nn; kwargs...)


text:
#InspectDR: Secondary annotation (Cairo layer)
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#
const ALIGN_MAP = Dict{Symbol, CAlignment}(
	:tl => ALIGN_TOP | ALIGN_LEFT,
	:tc => ALIGN_TOP | ALIGN_HCENTER,
	:tr => ALIGN_TOP | ALIGN_RIGHT,

	:cl => ALIGN_VCENTER | ALIGN_LEFT,
	:cc => ALIGN_VCENTER | ALIGN_HCENTER,
	:cr => ALIGN_VCENTER | ALIGN_RIGHT,

	:bl => ALIGN_BOTTOM | ALIGN_LEFT,
	:bc => ALIGN_BOTTOM | ALIGN_HCENTER,
	:br => ALIGN_BOTTOM | ALIGN_RIGHT,
)

#==Rendering text annotation
===============================================================================#
function render(canvas::PCanvas2D, a::TextAnnotation, graphinfo::Graph2DInfo)
	align = get(ALIGN_MAP, a.align, ALIGN_BOTTOM | ALIGN_LEFT)
	angle = deg2rad(a.angle)
	pt = map2dev(a.pos, canvas.xf, graphinfo.ixf, canvas.graphbb)
	render(canvas.ctx, a.text, pt, a.font, angle=angle, align=align)
	return
end

#==Rendering markers
===============================================================================#

function render(canvas::PCanvas2D, mkr::HVMarker, graphinfo::Graph2DInfo)
	ctx = canvas.ctx #WANTCONST
	graphbb = canvas.graphbb #WANTCONST
	if :none == mkr.line.style
		return
	end

	setlinestyle(ctx, LineStyle(mkr.line))
	pt = read2axis(mkr.pos, graphinfo.ixf)
	pt = map2dev(canvas.xf, pt)

	if mkr.vmarker
		drawline(ctx, Point2D(pt.x, graphbb.ymin), Point2D(pt.x, graphbb.ymax))
	end
	if mkr.hmarker
		drawline(ctx, Point2D(graphbb.xmin, pt.y), Point2D(graphbb.xmax, pt.y))
	end

	return
end

#==Rendering polyline for annotation
===============================================================================#
function render(canvas::PCanvas2D, a::PolylineAnnotation, graphinfo::Graph2DInfo)
	ctx = canvas.ctx #WANTCONST

	x = a.x; y = a.y
	setlinestyle(ctx, LineStyle(a.line))
	pt = map2dev(canvas.xf, Point2D(x[1], y[1]))
	Cairo.move_to(ctx, pt.x, pt.y)
	for i in 2:length(x)
		pt = map2dev(canvas.xf, Point2D(x[i], y[i]))
		Cairo.line_to(ctx, pt.x, pt.y)
	end
	if a.closepath
		Cairo.close_path(ctx)
		renderfill(ctx, a.fillcolor)
	end
	Cairo.stroke(ctx)

#	function xf(x, y)
#		pt = map2dev(canvas.xf, map2axis(Point2D(x,y), graphinfo.ixf))
#		return (x, y)
#	end
	return
end

#==Generic rendering algorithm for annotations
===============================================================================#
#TODO: warn undefined???
render(canvas::PCanvas2D, a::PlotAnnotation, graphinfo::Graph2DInfo) = nothing

#Basic render function for PlotAnnotation types:
function render(canvas::PCanvas2D, a::PlotAnnotation, graphinfo::Graph2DInfo, strip::Int)
	if 0 == a.strip || a.strip == strip
		render(canvas, a, graphinfo)
	end
	return
end

function render(canvas::PCanvas2D, alist::Vector{T}, graphinfo::Graph2DInfo, strip::Int) where T<:PlotAnnotation
	for a in alist
		render(canvas, a, graphinfo, strip)
	end
	return
end

#Last line


text:
using Blosc2
import Blosc2.Lib as lb
@testset "init destroy" begin
    lb.blosc_init()
    lb.blosc_destroy()
    @test 1 == 1
end

@testset "blosc_compress / decompress" begin
    lb.blosc_init()
    n = 10000
    data = rand(1:1000, n)
    buff = Vector{UInt8}(undef, sizeof(data) + lb.BLOSC_MAX_OVERHEAD)
    sz = lb.blosc_compress(9, lb.BLOSC_BITSHUFFLE, sizeof(Int64), sizeof(data), data, buff, sizeof(buff))
    @test sz > 0
    res = Vector{Int64}(undef, n + 500)
    sz = lb.blosc_decompress(buff, res, sizeof(res))
    @test sz == sizeof(data)
    @test data == res[1:div(sz, sizeof(Int64))]
    lb.blosc_destroy()
end

@testset "blosc_getitem" begin
    lb.blosc_init()
    n = 10000
    data = rand(1:1000, n)
    buff = Vector{UInt8}(undef, sizeof(data) + lb.BLOSC_MAX_OVERHEAD)
    sz = lb.blosc_compress(9, lb.BLOSC_BITSHUFFLE, sizeof(Int64), sizeof(data), data, buff, sizeof(buff))
    res = Int64[1, 2]
    for i = 1 : div(n, 2)
        sz = lb.blosc_getitem(buff, (i - 1) * 2, 2, res)
        @test sz == sizeof(res)
        @test res == data[(i - 1)*2 + 1:i*2]
    end
    lb.blosc_destroy()
end

@testset "blosc_get_nthreads / set_nthreads" begin
    lb.blosc_init()
    n = lb.blosc_get_nthreads()
    @test n == 1
    n = lb.blosc_set_nthreads(2)
    @test n == 1
    n = lb.blosc_get_nthreads()
    @test n == 2
    lb.blosc_destroy()
end

@testset "blosc_get_compressor / set_compressor" begin
    lb.blosc_init()
    r = lb.blosc_get_compressor()
    @test r == "blosclz"
    @test lb.blosc_set_compressor("lz4") >= 0
    r = lb.blosc_get_compressor()
    @test r == "lz4"
    lb.blosc_destroy()
end

@testset "set_delta / set_blocksize / free_resouces" begin
    lb.blosc_init()
    lb.blosc_set_delta(1)
    lb.blosc_set_blocksize(1000)
    r = lb.blosc_free_resources()
    @test r == 0
    lb.blosc_destroy()
end

@testset "blosc_cbufer_sizes / metainfo / versions / complib" begin
    lb.blosc_init()
    n = 10000
    data = rand(1:1000, n)
    buff = Vector{UInt8}(undef, sizeof(data) + lb.BLOSC_MAX_OVERHEAD)
    sz = lb.blosc_compress(9, lb.BLOSC_BITSHUFFLE, sizeof(Int64), sizeof(data), data, buff, sizeof(buff))
    sizes = lb.blosc_cbuffer_sizes(buff)
    @test sizes.nbytes == sizeof(data)
    @test sizes.cbytes == sz
    meta = lb.blosc_cbuffer_metainfo(buff)
    @test meta.typesize == 8
    versions = lb.blosc_cbuffer_versions(buff)
    @test 1 == 1
    complib = lb.blosc_cbuffer_complib(buff)
    @test complib == "LZ4"
    lb.blosc_destroy()
end

@testset "utilites" begin
    lb.blosc_init()
    @test lb.blosc_compcode_to_compname(0) == "blosclz"
    @test_throws ErrorException lb.blosc_compcode_to_compname(100)
    @test lb.blosc_compname_to_compcode("blosclz") == 0
    @test lb.blosc_compname_to_compcode("blosclz_rrrr") == -1
    @test length(lb.blosc_list_compressors()) > 1
    v = lb.blosc_get_version_string()
    @test !isempty(v)

    info = lb.blosc_get_complib_info("blosclz")
    @test info.code == 0
    @test info.complib == "BloscLZ"

    lb.blosc_destroy()
end


text:
using VT100
using Test

failed_tests = 0
for test in [
    "basic/hello", "basic/hello-nocr", "basic/long",
    "basic/color", "basic/boserase"
]
    em = ScreenEmulator(80,24)
    basepath = joinpath(dirname(@__FILE__),test)
    if isfile(string(basepath,".raw.in"))
        open(f->(parseall!(em,f)),string(basepath,".raw.in"))
    else
        @assert isfile(string(basepath,".in"))
        open(f->(parseall!(em,IOBuffer(unescape_string(read(f, String))))),
            string(basepath,".in"))
    end
    output = open(read,string(basepath,".out"))
    buf = IOBuffer()
    VT100.dump(buf,devnull,em)
    outbuf = take!(buf)
    if outbuf != output
        printstyled("Failed test $test\n"; color=:red)
        global failed_tests += 1
        println(output)
        println(outbuf)
        println(String(output))
        println(String(outbuf))
    end
end

exit(failed_tests)


text:
abstract type GeometricIntegratorAlgorithm <: DiffEqBase.DEAlgorithm end
struct GIEuler <: GeometricIntegratorAlgorithm end
struct GIMidpoint <: GeometricIntegratorAlgorithm end
struct GIHeun <: GeometricIntegratorAlgorithm end
struct GIKutta <: GeometricIntegratorAlgorithm end
struct GIERK4 <: GeometricIntegratorAlgorithm end
struct GIERK438 <: GeometricIntegratorAlgorithm end
struct GIImplicitEuler <: GeometricIntegratorAlgorithm end
struct GIImplicitMidpoint <: GeometricIntegratorAlgorithm end
struct GIRadIIA2 <: GeometricIntegratorAlgorithm end
struct GIRadIIA3 <: GeometricIntegratorAlgorithm end
struct GISRK3 <: GeometricIntegratorAlgorithm end
struct GIGLRK <: GeometricIntegratorAlgorithm
    s::Int
end
struct GISymplecticEulerA <: GeometricIntegratorAlgorithm end
struct GISymplecticEulerB <: GeometricIntegratorAlgorithm end
struct GILobattoIIIAIIIB2 <: GeometricIntegratorAlgorithm end
struct GILobattoIIIBIIIA2 <: GeometricIntegratorAlgorithm end


text:
import Base: show

"""
Construct a data structure describing a fluid-like weak scatterer with a deformed
cylindrical shape.  This shape is approximated by a series of `N` discrete segments,
described by the [x, y, z] coordinates, radii, and material properties at their
endpoints.

#### Parameters
- `r` : 3x`N` matrix describing the scatterer's centerline.  Each column is the
location in 3-D space of one of the centerline points (these should be arranged
in the proper order).
- `a` : Vector of radii, length `N`.
- `h`, `g` : Vectors of sound speed and density contrasts (i.e., the ratio
of sound speed or density inside the scatter to the same quantity in the
surrounding medium).
"""
:Scatterer
type Scatterer{T}
	r::Array{T, 2}
	a::Array{T, 1}
	h::Array{T, 1}
	g::Array{T, 1}
end

function Scatterer(r::Array, a::Array, h::Array, g::Array)
	T = Base.promote_eltype(r, a, h, g)
	r = convert(Matrix{T}, r)
	a = convert(Vector{T}, a)
	h = convert(Vector{T}, h)
	g = convert(Vector{T}, g)
	return Scatterer(r, a, h, g)
end

function Scatterer(r::Matrix, a::Vector, h::Real, g::Real)
	g1 = g * ones(length(a))
	h1 = h * ones(length(a))
	return Scatterer(r, a, h1, g1)
end

function show(io::IO, s::Scatterer)
	println("$(typeof(s)) with $(length(s.a)) segments")
	print("Length $(signif(length(s), 3))")
end

"""
Return the length of the scatterer (cartesian distance from one end to the other).
"""
length(s::Scatterer) = norm(s.r[:, 1] - s.r[:, end])

copy(s::Scatterer) = Scatterer(s.r, s.a, s.h, s.g)

"""
Scale the scatterer's size (overall or along a particular dimension) by a
constant factor.

#### Parameters
- `s` : Scatterer object.
- `scale` : Factor by which to grow/shrink the scatterer.
- `radius`, `x`, `y`, `z` : Optional factors, scaling the scatterer's radius
and along each dimension in space. All default to 1.0.

#### Returns
A rescaled scatterer.

#### Details
When making a scatterer larger, it is important to make sure it's body has enough
segments to accurately represent the shape at the frequencies of interest.
Specifically, the ratio L / (N λ), where L is the length of the animal, N is the
number of segments, and λ is the acoustic wavelength, should remain constant, which
may require interpolating new points between the existing ones. See Conti and
Demer (2006) or Calise and Skaret (2011) for details.
"""
function rescale(s::Scatterer; scale=1.0, radius=1.0, x=1.0, y=1.0, z=1.0)
	s = copy(s)
	M = diagm([x, y, z]) * scale
	s.r = M * s.r
	s.a = s.a * scale * radius
	return s
end

rescale(s::Scatterer, scale) = rescale(s, scale=scale)

"""
Resize a scatterer.  This is a convenience wrapper around `rescale`, for the
common situation where you want to change a scatterer to a specific length.
The scatterer's relative proportions are preserved.

#### Parameters
- `s` : Scatterer
- `len` : Desired length to which the scatterer should be scaled.

#### Returns
A resized scatterer.
"""
function resize(s::Scatterer, len)
	L0 = length(s)
	return rescale(s, len / L0)
end

"""
Resample a scatterer's measurement points by interpolating between them. Used
to change the resolution, for instance when increasing the scatterer's body
size or decreasing the acoustic wavelength.

#### Parameters
- `s` : Scatterer
- `n` : Number of body segments desired in the interpolated scatterer.

#### Returns
A Scatterer with a different number of body segments.
"""
function interpolate(s::Scatterer, n)
	# Length along scatterer's centerline
	along = [0, cumsum(vec(sqrt(sum(diff(s.r, 2).^2, 1))));]
	new_along = linspace(0, maximum(along), n)
	x = interpolate((along,), vec(s.r[1, :]), Gridded(Linear()))[new_along]
	y = interpolate((along,), vec(s.r[2, :]), Gridded(Linear()))[new_along]
	z = interpolate((along,), vec(s.r[3, :]), Gridded(Linear()))[new_along]
	a = interpolate((along,), s.a, Gridded(Linear()))[new_along]
	h = interpolate((along,), s.h, Gridded(Linear()))[new_along]
	g = interpolate((along,), s.g, Gridded(Linear()))[new_along]
	r = [x'; y'; z']
	return Scatterer(r, a, h, g)
end

"""
Rotate the scatterer in space, returning a rotated copy.

#### Parameters
- `roll` : Angle to roll the scatterer, in degrees. Defaults to 0.
- `tilt` : Angle to tilt the scatterer, in degrees. Defaults to 0.
- `yaw` : Angle to yaw the scatterer, in degrees. Defaults to 0.

#### Returns
A Scatterer with the same shape and properties, but a new orientation.

The roll, tilt, and yaw refer to rotations around the x, y, and z axes,
respectively. They are applied in that order.
"""
function rotate(s::Scatterer; roll=0.0, tilt=0.0, yaw=0.0)
	roll, tilt, yaw = deg2rad([roll, tilt, yaw])
	Rx = [1 0 0; 0 cos(roll) -sin(roll); 0 sin(roll) cos(roll)]
	Ry = [cos(tilt) 0 sin(tilt); 0 1 0; -sin(tilt) 0 cos(tilt)]
	Rz = [cos(yaw) -sin(yaw) 0; sin(yaw) cos(yaw) 0; 0 0 1]
	R = Rz * Ry * Rx
	s1 = copy(s)
	s1.r = R * s1.r
	return s1
end

function DWBAintegrand(s, rr, aa, gg, hh, k)
	r = vec(rr[:, 1] + s * diff(rr, 2))
	a = aa[1] + s * diff(aa)[1]
	g = gg[1] + s * diff(gg)[1]
	h = hh[1] + s * diff(hh)[1]
	alphatilt = acos(dot(k, r) / (norm(k) * norm(r)))
	betatilt = abs(alphatilt - pi/2.0)
	gammagamma = 1.0 / (g * h^2.0) + 1 / g - 2.0
	if abs(abs(betatilt) - pi/2.0) < 1e-10
		bessy = norm(k) * a / h
	else
		bessy = besselj1(2.0*norm(k) * a / h * cos(betatilt)) / cos(betatilt)
	end
	return (norm(k) / 4.0 * gammagamma * exp(2.0 * im * dot(k,r) / h) * a *
		bessy * norm(diff(rr, 2)))
end

scattering_function_param_docs = """
#### Parameters

- `s` : Scatterer object.
- `k` : Acoustic wavenumber vector.  Its magnitude is 2 * pi * f / c (where
f is the frequency and c is the sound speed) and it points in
the direction of propagation.  For downward-propagating sound waves,
it is [0, 0, -2pi * f / c].
- `phase_sd` : Standard deviation of the phase variability for each segment (in radians).
Defaults to 0.0, that is, an ordinary deterministic DWBA.  If > 0, the
return value will be stochastic (i.e., the SDWBA).
"""

"""
Calculate the complex-valued form function of a scatterer using the (S)DWBA.

$scattering_function_param_docs
"""
:form_function
function form_function(s::Scatterer, k::Vector, phase_sd=0.0)
	fbs = 0 + 0im
	m, n = size(s.r)

	for i in 1:(n-1)
		f(x) = DWBAintegrand(x, s.r[:, i:i+1], s.a[i:i+1],
			s.g[i:i+1], s.h[i:i+1], k)
		fbs += quadgk(f, 0.0, 1.0)[1] * exp(im * randn() * phase_sd)
	end
	return fbs
end

"""
Calculate the backscattering cross-section (sigma_bs) of a scatterer using the (S)DWBA.
This is the absolute square of the form function.

$scattering_function_param_docs
"""
function backscatter_xsection{T}(s::Scatterer{T}, k::Vector{T}, phase_sd=0.0)
	return abs2(form_function(s, k, phase_sd))
end

"""
Calculate the target strength (TS) of a scatterer using the (S)DWBA.  This is just
10 * log10(sigma_bs).

$scattering_function_param_docs
"""
function target_strength{T}(s::Scatterer{T}, k::Vector{T}, phase_sd=0.0)
	return 10 * log10(backscatter_xsection(s, k, phase_sd))
end

for func in [:form_function, :backscatter_xsection, :target_strength]
	eval(quote
		function ($func)(s::Scatterer, freq::Real, sound_speed::Real, phase_sd=0.0)
			k = [0.0, 0.0, -2pi * freq / sound_speed]
			return ($func)(s, k, phase_sd)
		end
	end)
end

"""
Calculate backscatter over a range of angles.

#### Parameters

- `s` : Scatterer object
- `angle1`, `angle2` : Endpoints of the angle range to calculate.
- `k` : Acoustic wavenumber vector
- `n` : Number of angles to calculate; defaults to 100

#### Returns

A dictionary containing elements "angles", "sigma_bs", and "TS",
each a length-n vector.
"""
function tilt_spectrum(s::Scatterer, angle1, angle2, k, n=100)
	angles = linspace(angle1, angle2, n)
	sigma = zeros(angles)
	for i in 1:n
		tilt = angles[i]
		sigma[i] = backscatter_xsection(rotate(s, tilt=tilt), k)
	end
	TS = 10 * log10(sigma)
	return Dict([("angles", angles), ("sigma_bs", sigma), ("TS", TS)])
end

"""
Calculate backscatter over a range of frequencies.  The insonifying sound comes
from above (i.e., traveling in the -z direction).

#### Parameters
- `s` : Scatterer object
- `freq1`, `freq2` : Endpoints of the angle range to calculate.
- `sound_speed` : Sound speed in the surrounding medium
- `n` : Number of frequencies to calculate; defaults to 100

Returns: A dictionary containing elements "freqs", "sigma_bs", and "TS",
	each a length-n vector.
"""
function freq_spectrum(s::Scatterer, freq1, freq2, sound_speed, n=100)
	freqs = linspace(freq1, freq2, n)
	sigma = zeros(freqs)
	for i in 1:n
		k = [0, 0, -2pi * freqs[i] / sound_speed]
		sigma[i] = backscatter_xsection(s, k)
	end
	TS = 10 * log10(sigma)
	return Dict([("freqs", freqs), ("sigma_bs", sigma), ("TS", TS)])
end

"""
Load a scatterer from a file on disk with comma-separated values.

#### Parameters
- `filename` : String.  Path to the datafile.  This should be a standard .csv file
with columns for the x, y, and z coordinates of the scatterer's centerline, as well
as the `a`, `h`, and `g` arguments to Scatterer().
- `columns` : Optional dictionary of column names. If the columns do not have the names
- `x`, `y`, `z`, `h`, and `g`, this must be provided.  The keys are the standard column
names and the values are the actual ones in the file.
"""
function from_csv(filename, columns=Dict([("x","x"),("y","y"),("z","z"),
		("a","a"), ("h","h"), ("g","g")]))
	data, header = readdlm(filename, ',', header=true)
	x = data[:, vec(header .== columns["x"])]
	y = data[:, vec(header .== columns["y"])]
	z = data[:, vec(header .== columns["z"])]
	a = vec(data[:, vec(header .== columns["a"])])
	h = vec(data[:, vec(header .== columns["h"])])
	g = vec(data[:, vec(header .== columns["g"])])
	r = [x y z]'
	return Scatterer(r, a, h, g)
end

"""
Save a scatterer's shape to a file on disk with comma-separated values.

#### Parameters
- `s` : Scatterer object to save.
- `filename` : Where to save it.

"""
function to_csv(s::Scatterer, filename)
	header = ["x" "y" "z" "a" "h" "g"]
	data = [s.r' s.a s.h s.g]
	writedlm(expanduser(filename), [header; data], ',')
end


text:
using ReverseDiff
using MNIST

type Dataset
  X::AbstractArray
  y::AbstractArray
end

type MNISTData
  train::Dataset
  dev::Dataset
  test::Dataset
end

function onehot{T}(y::Vector{T}, n)
  yy = zeros(T, n, length(y))
  for i=1:length(y)
    label = y[i]
    yy[Int(label)+1, i] = 1
  end
  return yy
end

function load_mnist()
  X, y = MNIST.traindata()
  X_train, y_train = X[:, 1:50000], y[1:50000]
  X_dev, y_dev = X[:, 50001:end], y[50001:end]
  X_test, y_test = MNIST.testdata()
  return MNISTData(
    Dataset(X_train/255, onehot(y_train, 10)),
    Dataset(X_dev/255, onehot(y_dev, 10)),
    Dataset(X_test/255, onehot(y_test, 10))
  )
end

function sgd_update(v, Δv, α)
  for i in eachindex(v)
    v[i] -= α * Δv[i]
  end
end

data = load_mnist()

const BATCH_SIZE = 100
const STEPS_PER_EPOCH = floor(size(data.train.X, 2) / BATCH_SIZE)
const EPOCHS = 20
α = 1e-3

# parameters
W = randn(Float64, 10, 784) * 0.01
b = zeros(Float64, 10)

inputs = (zeros(Float64, 784, BATCH_SIZE), zeros(Float64, 10, BATCH_SIZE), W, b)
outputs = map(similar, inputs)

ReverseDiff.@forward minus_log(x) = -log(x)
softmax(x) = (exp_x = exp.(x); exp_x ./ sum(exp_x, 1))
cross_entropy_loss(x, y) = mean(sum(y .* minus_log.(x), 1))

model(x, y, W, b) = cross_entropy_loss(softmax(W*x .+ b), y)
predict(x, W, b) = W*x .+ b
accuracy(y, ŷ) = mean(mapslices(indmax, y, 1) .== mapslices(indmax, ŷ, 1))

@time model! = ReverseDiff.compile_gradient(model, inputs)

for i=1:EPOCHS
  @time for j=1:STEPS_PER_EPOCH
    s = Int((j-1) * BATCH_SIZE + 1)
    e = Int(j * BATCH_SIZE)
    X_batch = @view data.train.X[:, s:e]
    y_batch = @view data.train.y[:, s:e]
    inputs = (X_batch, y_batch, W, b)
    model!(outputs, inputs)
    sgd_update(W, outputs[3], α)
    sgd_update(b, outputs[4], α)
  end

  dev_acc = accuracy(predict(data.dev.X, W, b), data.dev.y)
  println("Epoch $(i), Dev accuracy = $(dev_acc)")
end

test_acc = accuracy(predict(data.test.X, W, b), data.test.y)
println("Test accuracy = $(test_acc)")


text:
## PART 2

include("meansofbins.jl");

# Sample sizes
N1 = 1000;       # One thousand
N2 = 10000;      # Ten thousand
println("\nPart 2");
println("\n\tN1 = $N1\n\tN2 = $N2\n")

# Leave space for the means v̄1 and v̄2
v̄1 = zeros(int64(M/N1),5);
v̄2 = zeros(int64(M/N2),5);

# Save all those means
for a in [1:5]
    v̄1[:,a] = meansofbins(N1, M, v[:,a]);
    v̄2[:,a] = meansofbins(N2, M, v[:,a]);
end

# Declare standard devs of the sample means
σ̂1 = Float64[];
σ̂2 = Float64[];

# Calculate standard devs for sample means
for a in [1:5]
    push!(σ̂1,std(v̄1[:,a]));
    push!(σ̂2,std(v̄2[:,a]));
end

# Print the results
println("\t Standard deviations of the sample means are:\n");
println("\tσ̂1 = $σ̂1");
println("\tσ̂2 = $σ̂2");

# See if this is what we'd expect
println("\n\tWe expect (σ̂1./σ̂2).^2 = N2/N1 = 10:\n");
for i in [1:5]
    ratiosqr = (σ̂1[i]/σ̂2[i])^2
    println("\t(σ̂1[$i]/σ̂2[$i])^2 = $ratiosqr");
end
println("\n\t...which is pretty close.");

# Histogram for dataset 1
println("\n\tMaking histograms and placing in 'histograms' vector")
using Gadfly
histograms = Plot[];
for j in [1:5]
  push!(
    histograms,
    plot(
      layer(
          x = v̄1[:,j],
          color = ["N1 = 1,000" for i in 1:25],
          Geom.histogram(bincount=25, density=true),
          order=1
      ),
      layer(
          x = v̄2[:,j],
          color = ["N2 = 10,000" for i in 1:25],
          Geom.histogram(bincount=25, density=true),
          order=2,
          Theme(default_color=color("red"))
      ),
      Guide.xlabel("Sample mean"),
      Guide.ylabel("Relative Frequency"),
      Guide.title(join(["Sample Mean Histogram in Dataset",string(j)]," ")),
      Guide.colorkey("Legend title"),
    )
  )
end


text:
using ProgressMeter
using .Threads
using Lasso
using NPZ
using Statistics
using Random
using StatsBase
using LinearAlgebra

struct InvResult{T,S}
    i::Int
    λmin::T
    dn::T
    m::S
end

function fit_model(G,dv,i)
    dn = std(dv)
    dv ./= dn
    path = fit(LassoPath, G, dv, λminratio = 1e-5, standardize=false, intercept=false, maxncoef=size(G,2))
    λmin = path.λ[argmin(aicc(path))]
    m = coef(path; select=MinAICc())
    return InvResult(i, λmin, dn, m)
end

T = Float64

G1 = npzread("Kernels/G1.npy")
G2 = npzread("Kernels/G2.npy")
G3 = npzread("Kernels/G4.npy")

d1 = npzread("Testing/das_wvt_data.npy")
d2 = npzread("Testing/nodal_wvt_data_x.npy")
d3 = npzread("Testing/nodal_wvt_data_y.npy")
nfactor = npzread("Testing/nfactor.npy")

G = [nfactor*G1 nfactor*G2;
     G3 zeros(T, size(G3));
     zeros(T, size(G3)) G3]

Gn = std(G)

G ./= Gn

d = [d1; d2; d3]

s_size = [4, 8, 16, 32, 64]

p = Progress(length(s_size)*size(d,2), dt=0.5, barglyphs=BarGlyphs("[=> ]"), barlen=50, color=:yellow)

for B in s_size

    pvec = npzread("IncoherenceResults/res_inc_$B.npy")
    sensors = sortperm(pvec)[end-B+1:end]

    μl = vcat(1:size(G1, 1), size(G1, 1).+sensors, size(G1, 1).+sensors.+64)

    Gμ = G[μl,:]
    results = Array{Any,1}(undef, size(d,2))

    @threads for j = 1:size(d,2)
        results[j] = fit_model(Gμ, d[μl,j], j)
        next!(p)
    end

    sres = reduce(hcat, [x.m.*x.dn for x in results])

    npzwrite("Testing/LassoOutputInc/32B_$(B)_max.npy", Array{Float64}(sres))
end

max_stars_repo_name: jbmuir/DAS-Reconstruction

max_issues_repo_name: jbmuir/DAS-Reconstruction

max_forks_repo_name: jbmuir/DAS-Reconstruction


text:
# # Decaying barotropic QG beta-plane turbulence
#
#md # This example can be run online via [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/barotropicqg_betadecay.ipynb).
#md # Also, it can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/barotropicqg_betadecay.ipynb).
#
# An example of decaying barotropic quasi-geostrophic turbulence on a beta plane.

using FourierFlows, Plots, Printf, Random

using Statistics: mean
using FFTW: irfft

import GeophysicalFlows.BarotropicQG
import GeophysicalFlows.BarotropicQG: energy, enstrophy

# ## Choosing a device: CPU or GPU

dev = CPU()     # Device (CPU/GPU)
nothing # hide

# ## Numerical parameters and time-stepping parameters

      n = 128            # 2D resolution = n²
stepper = "FilteredRK4"  # timestepper
     dt = 0.05           # timestep
 nsteps = 1500           # total number of time-steps
 nsubs  = 10             # number of time-steps for intermediate logging/plotting (nsteps must be multiple of nsubs)
nothing # hide

# ## Physical parameters

L = 2π        # domain size
β = 10.0      # planetary PV gradient
μ = 0.0       # bottom drag
nothing # hide

# ## Problem setup
# We initialize a `Problem` by providing a set of keyword arguments. Not providing
# a viscosity coefficient `ν` leads to the module's default value: `ν=0`. In this
# example numerical instability due to accumulation of enstrophy at high wavenumbers
# is taken care with the `FilteredTimestepper` we picked.
prob = BarotropicQG.Problem(dev; nx=n, Lx=L, β=β, μ=μ, dt=dt, stepper=stepper)
nothing # hide

# and define some shortcuts
sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
x, y = grid.x, grid.y
nothing # hide

# ## Setting initial conditions

# Our initial condition consist of a flow that has power only at wavenumbers with
# ``8 < \frac{L}{2\pi} \sqrt{k_x^2 + k_y^2} < 10`` and initial energy ``E_0``:

E₀ = 0.1 # energy of initial condition

K = @. sqrt(grid.Krsq)                          # a 2D array with the total wavenumber
k = [grid.kr[i] for i=1:grid.nkr, j=1:grid.nl]  # a 2D array with the zonal wavenumber

Random.seed!(1234)
qih = randn(Complex{eltype(grid)}, size(sol))
@. qih = ifelse(K < 2  * 2π/L, 0, qih)
@. qih = ifelse(K > 10 * 2π/L, 0, qih)
@. qih = ifelse(k == 0 * 2π/L, 0, qih)   # no power at zonal wavenumber k=0 component
Ein = energy(qih, grid)           # compute energy of qi
qih *= sqrt(E₀ / Ein)             # normalize qi to have energy E₀
qi = irfft(qih, grid.nx)

BarotropicQG.set_zeta!(prob, qi)
nothing # hide

# Let's plot the initial vorticity field:

p1 = heatmap(x, y, vars.q',
         aspectratio = 1,
              c = :balance,
           clim = (-12, 12),
          xlims = (-grid.Lx/2, grid.Lx/2),
          ylims = (-grid.Ly/2, grid.Ly/2),
         xticks = -3:3,
         yticks = -3:3,
         xlabel = "x",
         ylabel = "y",
          title = "initial vorticity ζ=∂v/∂x-∂u/∂y",
     framestyle = :box)

p2 = contourf(x, y, vars.psi',
        aspectratio = 1,
             c = :viridis,
        levels = range(-0.65, stop=0.65, length=10),
          clim = (-0.65, 0.65),
         xlims = (-grid.Lx/2, grid.Lx/2),
         ylims = (-grid.Ly/2, grid.Ly/2),
        xticks = -3:3,
        yticks = -3:3,
        xlabel = "x",
        ylabel = "y",
         title = "initial streamfunction ψ",
    framestyle = :box)

l = @layout Plots.grid(1, 2)
p = plot(p1, p2, layout=l, size = (900, 800))

# ## Diagnostics

# Create Diagnostics -- `energy` and `enstrophy` functions are imported at the top.
E = Diagnostic(energy, prob; nsteps=nsteps)
Z = Diagnostic(enstrophy, prob; nsteps=nsteps)
diags = [E, Z] # A list of Diagnostics types passed to "stepforward!" will  be updated every timestep.
nothing # hide

# ## Output

# We choose folder for outputing `.jld2` files and snapshots (`.png` files).
filepath = "."
plotpath = "./plots_decayingbetaturb"
plotname = "snapshots"
filename = joinpath(filepath, "decayingbetaturb.jld2")
nothing # hide

# Do some basic file management,
if isfile(filename); rm(filename); end
if !isdir(plotpath); mkdir(plotpath); end
nothing # hide

# and then create Output.
get_sol(prob) = sol # extracts the Fourier-transformed solution
out = Output(prob, filename, (:sol, get_sol))
nothing # hide

# ## Visualizing the simulation

# We define a function that plots the vorticity and streamfunction and
# their corresponding zonal mean structure.

function plot_output(prob)
  ζ = prob.vars.zeta
  ψ = prob.vars.psi
  ζ̄ = mean(ζ, dims=1)'
  ū = mean(prob.vars.u, dims=1)'

  pζ = heatmap(x, y, ζ',
       aspectratio = 1,
            legend = false,
                 c = :balance,
              clim = (-12, 12),
             xlims = (-grid.Lx/2, grid.Lx/2),
             ylims = (-grid.Ly/2, grid.Ly/2),
            xticks = -3:3,
            yticks = -3:3,
            xlabel = "x",
            ylabel = "y",
             title = "vorticity ζ=∂v/∂x-∂u/∂y",
        framestyle = :box)

  pψ = contourf(x, y, ψ',
       aspectratio = 1,
            legend = false,
                 c = :viridis,
            levels = range(-0.65, stop=0.65, length=10),
              clim = (-0.65, 0.65),
             xlims = (-grid.Lx/2, grid.Lx/2),
             ylims = (-grid.Ly/2, grid.Ly/2),
            xticks = -3:3,
            yticks = -3:3,
            xlabel = "x",
            ylabel = "y",
             title = "streamfunction ψ",
        framestyle = :box)

  pζm = plot(ζ̄, y,
            legend = false,
         linewidth = 2,
             alpha = 0.7,
            yticks = -3:3,
             xlims = (-2.2, 2.2),
            xlabel = "zonal mean ζ",
            ylabel = "y")
  plot!(pζm, 0*y, y, linestyle=:dash, linecolor=:black)

  pum = plot(ū, y,
            legend = false,
         linewidth = 2,
             alpha = 0.7,
            yticks = -3:3,
             xlims = (-0.55, 0.55),
            xlabel = "zonal mean u",
            ylabel = "y")
  plot!(pum, 0*y, y, linestyle=:dash, linecolor=:black)

  l = @layout Plots.grid(2, 2)
  p = plot(pζ, pζm, pψ, pum, layout = l, size = (900, 800))

  return p
end
nothing # hide

# ## Time-stepping the `Problem` forward

# We time-step the `Problem` forward in time.

startwalltime = time()

p = plot_output(prob)

anim = @animate for j = 0:round(Int, nsteps/nsubs)

  if j % (500 / nsubs) == 0
    cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy])

    log = @sprintf("step: %04d, t: %d, cfl: %.2f, E: %.4f, Q: %.4f, walltime: %.2f min",
      clock.step, clock.t, cfl, E.data[E.i], Z.data[Z.i], (time()-startwalltime)/60)

    println(log)
  end

  p[1][1][:z] = vars.zeta
  p[1][:title] = "vorticity, t="*@sprintf("%.2f", clock.t)
  p[3][1][:z] = vars.psi
  p[2][1][:x] = mean(vars.zeta, dims=1)'
  p[4][1][:x] = mean(vars.u, dims=1)'

  stepforward!(prob, diags, nsubs)
  BarotropicQG.updatevars!(prob)

end

gif(anim, "barotropicqg_betadecay.gif", fps=8)

# ## Save

# Finally save the last snapshot.
savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), clock.step)
savefig(savename)


text:
"""
$(TYPEDEF)
Apex type for coordinate systems
"""
abstract type AbstractCoordinateSystem <: AbstractExtendableGridApexType end

"""
$(TYPEDSIGNATURES)

List possible coordinate systems. These describe the meaning of the grid coordinates.
"""
coordinatesystems()=AbstractTrees.print_tree(AbstractCoordinateSystem,5,indicate_truncation=false)

"""
$(TYPEDEF)
1D cartesion coordinate system (unknown `x`)
"""
abstract type Cartesian1D   <: AbstractCoordinateSystem end

"""
$(TYPEDEF)
2D cartesion coordinate system (unknowns `x,y`)
"""
abstract type Cartesian2D   <: AbstractCoordinateSystem end

"""
$(TYPEDEF)
2D cartesion coordinate system (unknowns `x,y,z`)
"""
abstract type Cartesian3D   <: AbstractCoordinateSystem end

"""
$(TYPEDEF)
3D cylindrical coordinate system (unknowns `r,ϕ,z`)
"""
abstract type Cylindrical3D <: AbstractCoordinateSystem end

"""
$(TYPEDEF)
2D cylindrical coordinate system (unknowns `r,z`)
"""
abstract type Cylindrical2D <: AbstractCoordinateSystem end

"""
$(TYPEDEF)
2D polar coordinate system (unknowns `r,ϕ`)
"""
abstract type Polar2D       <: AbstractCoordinateSystem end

"""
$(TYPEDEF)
1D polar coordinate system (unknown `r`)
"""
abstract type Polar1D       <: AbstractCoordinateSystem end

"""
$(TYPEDEF)
3D spheriacal coordinate system (unknowns `r,ϕ,θ`)
"""
abstract type Spherical3D   <: AbstractCoordinateSystem end

"""
$(TYPEDEF)
1D spheriacal coordinate system (unknown `r`)
"""
abstract type Spherical1D   <: AbstractCoordinateSystem end  #r (integral over ϕ,\theta)


text:
module PMC

#=
int max_clique(long long nedges, int *ei, int *ej, int index_offset,
                int outsize, int *clique);
=#

const libpmc = joinpath(dirname(@Base.__FILE__),"libpmc")

function pmc(A::SparseMatrixCSC)
 maxd = Int(maximum(sum(spones(A),1)))
 ei,ej = findnz(tril(A,1))

 ei = map(Int32,ei)
 ej = map(Int32,ej)

 offset = Cint(1)

 outsize = maxd
 output = zeros(Int32,maxd)

 clique_size = ccall(
    (:max_clique, libpmc), Cint,
    (Clonglong, Ptr{Cint}, Ptr{Cint}, Cint, Cint, Ptr{Cint}),
    length(ei), ei, ej, offset, outsize, output)

  return map(Int64, output[1:clique_size])
end

end

#@show PMC.pmc(sprandn(10000,10000,10/10000))

max_stars_repo_name: veeresht/pmc

max_issues_repo_name: veeresht/pmc

max_forks_repo_name: veeresht/pmc


text:
#==============================================================================#
# AWSCredentials.jl
#
# Load AWS Credentials from:
# - EC2 Instance Profile,
# - Environment variables, or
# - ~/.aws/credentials file.
#
# Copyright Sam O'Connor 2014 - All rights reserved
#==============================================================================#

export AWSCredentials,
       localhost_is_lambda,
       localhost_is_ec2,
       aws_user_arn,
       aws_account_number

type AWSCredentials
    access_key_id::String
    secret_key::String
    token::String
    user_arn::String
    account_number::String

    function AWSCredentials(access_key_id, secret_key,
                            token="", user_arn="", account_number="")
        new(access_key_id, secret_key, token, user_arn, account_number)
    end
end

function Base.show(io::IO,c::AWSCredentials)
    println(io, string(c.user_arn,
                       " (",
                       c.account_number,
                       c.access_key_id,
                       c.secret_key == "" ? "" : ", $(c.secret_key[1:3])...",
                       c.token == "" ? "" : ", $(c.token[1:3])..."),
                       ")")
end

# Discover AWS Credentials for local host.

function AWSCredentials()

    if localhost_is_ec2()

        return ec2_instance_credentials()

    elseif haskey(ENV, "AWS_ACCESS_KEY_ID")

        return env_instance_credentials()

    elseif isfile(dot_aws_credentials_file())

        return dot_aws_credentials()
    end

    error("Can't find AWS credentials!")
end

# Is Julia running in an AWS Lambda sandbox?

localhost_is_lambda() = haskey(ENV, "LAMBDA_TASK_ROOT")

# Is Julia running on an EC2 virtual machine?

function localhost_is_ec2()

    if localhost_is_lambda()
        return false
    end

    @static if is_unix()
        return isfile("/sys/hypervisor/uuid") &&
               bytestring(readbytes("/sys/hypervisor/uuid",3)) == "ec2"
    end

    return false
end

# Get User ARN for "creds".

function aws_user_arn(aws::AWSConfig)

    creds = aws[:creds]

    if creds.user_arn == ""

        r = do_request(post_request(aws, "sts", "2011-06-15",
                                    Dict("Action" => "GetCallerIdentity",
                                         "ContentType" => "JSON")))
        creds.user_arn = r["Arn"]
        creds.account_number = r["Account"]
    end

    return creds.user_arn
end

# Get Account Number for "creds".

function aws_account_number(aws::AWSConfig)
    creds = aws[:creds]
    if creds.account_number == ""
        aws_user_arn(aws)
    end
    return creds.account_number
end

# Fetch EC2 meta-data for "key".
# http://docs.aws.amazon.com/AWSEC2/latest/UserGuide/AESDG-chapter-instancedata.html

function ec2_metadata(key)

    @assert localhost_is_ec2()

    bytestring(http_request("169.254.169.254", "/latest/meta-data/$key").data)
end

# Load Instance Profile Credentials for EC2 virtual machine.

using JSON

function ec2_instance_credentials()

    @assert localhost_is_ec2()

    info  = ec2_metadata("iam/info")
    info  = JSON.parse(info, dicttype=Dict{String,String})

    name  = ec2_metadata("iam/security-credentials/")
    creds = ec2_metadata("iam/security-credentials/$name")
    new_creds = JSON.parse(creds, dicttype=Dict{String,String})

    AWSCredentials(new_creds["AccessKeyId"],
                   new_creds["SecretAccessKey"],
                   new_creds["Token"],
                   info["InstanceProfileArn"])
end

# Load Credentials from environment variables (e.g. in Lambda sandbox)

function env_instance_credentials()

    AWSCredentials(ENV["AWS_ACCESS_KEY_ID"],
                   ENV["AWS_SECRET_ACCESS_KEY"],
                   get(ENV, "AWS_SESSION_TOKEN", ""),
                   get(ENV, "AWS_USER_ARN", ""))
end

# Load Credentials from AWS CLI ~/.aws/credentials file.

using IniFile

dot_aws_credentials_file() = get(ENV, "AWS_CONFIG_FILE",
                                 joinpath(homedir(), ".aws", "credentials"))

function dot_aws_credentials()

    @assert isfile(dot_aws_credentials_file())

    ini = read(Inifile(), dot_aws_credentials_file())

    profile = get(ENV, "AWS_DEFAULT_PROFILE", "default")

    AWSCredentials(get(ini, profile, "aws_access_key_id"),
                   get(ini, profile, "aws_secret_access_key"))
end

#==============================================================================#
# End of file.
#==============================================================================#


text:
aircraft = data_epema_rotor_sweep[1]
parameters = data_epema_rotor_sweep[2]

struct TestParams
    us
    vs
    wake_function
    wake_shape_functions
    axial_interpolations
    swirl_interpolations
    axial_multipliers
    swirl_multipliers
end

nsteps = length(parameters.us)
stepi = 5

testparams = TestParams(deepcopy(parameters.us), deepcopy(parameters.vs), Vector{Any}(nothing,nsteps), fill((Rtip, x) -> Rtip, length(aircraft.rotor_system.index)), fill((rs, us, r, Rtip) -> FM.linear(rs, us, r), length(aircraft.rotor_system.index)), fill((rs, vs, r, Rtip) -> FM.linear(rs, vs, r), length(aircraft.rotor_system.index)), fill((distance2plane, Rtip) -> 2, length(aircraft.rotor_system.index)), fill((distance2plane, Rtip) -> 1, length(aircraft.rotor_system.index)))

for i=1:nsteps
    AS.solve_rotor_wake(aircraft, testparams, nothing, nothing, nothing, i, nothing)
end

wakefun_aoa4 = testparams.wake_function[2]
r = aircraft.rotor_system.rotors[1].Rtip
# thetas = range(0, 2*pi, length=12)
# rs = range(r/10, 2*r, length=4)
# x = [r * cos(theta) for theta in thetas, r in rs]
# y = [r * cos(theta) for theta in thetas, r in rs]
# Vs = wakefun_aoa4.()
x = range(-r * 2.0, stop = r * 8.5, length=3)
y = range(-r * 2.5, stop = r * 2.5, length=35)
z = range(-r * 2.5, stop = r * 2.5, length=35)
lx = length(x)
ly = length(y)
lz = length(z)
xs = [x[i] for i in 1:lx, j in 1:ly, k in 1:lz]
ys = [y[j] for i in 1:lx, j in 1:ly, k in 1:lz]
zs = [z[k] for i in 1:lx, j in 1:ly, k in 1:lz]
Vs = wakefun_aoa4.([[x[i], y[j], z[k]] for i in 1:lx, j in 1:ly, k in 1:lz])
# Vs ./ AS.LA.norm.(Vs)
us = [Vs[i,j,k][1] for i in 1:lx, j in 1:ly, k in 1:lz]
vs = [Vs[i,j,k][2] for i in 1:lx, j in 1:ly, k in 1:lz]
ws = [Vs[i,j,k][3] for i in 1:lx, j in 1:ly, k in 1:lz]

fig = plt.figure("test_wakefunction")
fig.clear()
ax = fig.add_subplot(111, projection="3d")
ax.quiver3D(xs / r, ys / r, zs / r, us, vs, ws, length=0.01)#, normalize=true)

ax.set_xlabel(LS.L"x/r")
ax.set_ylabel(LS.L"y/r")
ax.set_zlabel(LS.L"z/r")

fig2d = plt.figure("test_wakefunction_2d")
fig2d.clear()
ax2 = fig2d.add_subplot(111)
ax2.quiver(ys[1,:,:] / r, zs[1,:,:] / r, (vs[1,:,:] .+ 1e-10), (ws[1,:,:] .+ 1e-10))#, length=0.01, normalize=true)
thetas = range(0,2*pi,length=100)
rotor_y = cos.(thetas)
rotor_z = sin.(thetas)
ax2.plot(rotor_y, rotor_z, color="blue")
ax2.set_xlabel(LS.L"x/r")
ax2.set_ylabel(LS.L"y/r")

rs = range(r/200, 1.5*r, length=100)
xs_3 = vcat(-reverse(rs), rs)
Vs_3 = wakefun_aoa4.([[0.0, r, 0.0] for r in xs_3])
us = [v[1] for v in Vs_3]
vs = [v[2] for v in Vs_3]
ws = [v[3] for v in Vs_3]

fig3 = plt.figure("test_wakefunction_us_vs")
fig3.clear()
ax_u = fig3.add_subplot(311)
ax_v = fig3.add_subplot(312)
ax_w = fig3.add_subplot(313)
ax_u.plot(xs_3 ./ r, us)
ax_u.set_ylabel(LS.L"u[m/s]")
ax_v.plot(xs_3 ./ r, vs)
ax_v.set_ylabel(LS.L"v[m/s]")
ax_w.plot(xs_3 ./ r, ws)
ax_w.set_ylabel(LS.L"w[m/s]")
ax_w.set_xlabel(LS.L"y/R, x=z=0")

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
#=
Implement broadcast functionality for Sol
=#

Base.similar(sol::Sol{N}) where {N} =
    Sol(similar(sol.x), similar(sol.α))

Broadcast.broadcasted(::typeof(identity), s::Sol{N}) where {N} = s

function Broadcast.materialize!(dest::Sol{N}, src::Sol{N}) where {N}
    dest.x .= src.x
    dest.α .= src.α
end

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
function reg_term(omega, B = size(omega,2); alpha=0.01)
  reg = 0.0
  for b in 1:B
    t = maximum(@views omega[:, b]) - omega[b,b]
    reg += t
  end
  return alpha*reg
end

function reg_term_turbo(omega, B = size(omega,2); alpha=0.01)
  reg = 0.0
  @turbo for b in 1:B
    m = -Inf
    for i in axes(omega,1)
      m = max(m, omega[i,b])
    end
    t = m - omega[b,b]
    reg += t
  end
  return alpha*reg
end

function awmean_lv(x::AbstractArray{T1}, σ::AbstractArray{T2}) where {T1<:Number,T2<:Number}
  n = length(x)
  T3 = promote_type(T1,T2)
  T = sizeof(T3) ≤ 4 ? Float32 : Float64
  sum_of_values = sum_of_weights = χ2 = zero(T)
  @turbo for i=1:n
    sum_of_values += x[i] / (σ[i]*σ[i])
    sum_of_weights += 1 / (σ[i]*σ[i])
  end
  wx = sum_of_values / sum_of_weights

  @turbo for i=1:n
    χ2 += (x[i] - wx) * (x[i] - wx) / (σ[i] * σ[i])
  end
  mswd = χ2 / (n-1)
  wσ = sqrt(one(T) / sum_of_weights)
  return wx, wσ, mswd
end
function awmean_simd(x::AbstractArray{T1}, σ::AbstractArray{T2}) where {T1<:Number,T2<:Number}
  n = length(x)
  T3 = promote_type(T1,T2)
  T = sizeof(T3) ≤ 4 ? Float32 : Float64
  sum_of_values = sum_of_weights = χ2 = zero(T)
  @fastmath @inbounds @simd for i=1:n
    sum_of_values += x[i] / (σ[i]*σ[i])
    sum_of_weights += 1 / (σ[i]*σ[i])
  end
  wx = sum_of_values / sum_of_weights

  @fastmath @inbounds @simd for i=1:n
    χ2 += (x[i] - wx) * (x[i] - wx) / (σ[i] * σ[i])
  end
  mswd = χ2 / (n-1)
  wσ = sqrt(one(T) / sum_of_weights)
  return wx, wσ, mswd
end

function test_awmean(::Type{T}) where {T}
    for n ∈ 2:100
      if T <: Integer
        x = view(rand(T(-50):T(100), n + 32), 17:n+16)
        σ = view(rand(T(1):T(10), n + 32), 17:n+16)
      else
        x = view((randn(T, n + 32) .+= T(2)), 17:n+16)
        σ = view(rand(T, n + 32), 17:n+16)
      end
      wx, wσ, mswd = awmean_simd(x, σ)
      @test iszero(@allocated((wxlv, wσlv, mswdlv) = awmean_lv(x, σ)))
      wxlv, wσlv, mswdlv = awmean_lv(x, σ)
      isfinite(wx)   && @test wx ≈ wxlv
      isfinite(wσ)   && @test wσ ≈ wσlv
      isfinite(mswd) && @test mswd ≈ mswdlv atol=eps(typeof(mswdlv)) rtol=sqrt(eps(typeof(mswdlv)))
    end
end

function logℒ_fast(α, β, t, c, x)
    eα = abs(α)
    n, k = size(x)

    (n == length(t) == length(c) && length(β) == k + 1) || throw(DimensionMismatch())
    s = zero(typeof(α))
    @inbounds for i in 1:n
      xb = sum(@inbounds(x[i, j] * β[j+1]) for j in 1:k) + β[1]
      ti = t[i]
      s += (1 - (c[i] == ti)) * (log(eα) + (eα - 1) * log(ti) + xb) - ti^eα * exp(xb)
    end
    return s
end
function logℒ_fast_turbo(α, β, t, c, x)
    eα = abs(α)
    n, k = size(x)

    (n == length(t) == length(c) && length(β) == k + 1) || throw(DimensionMismatch())
    s = zero(typeof(α))
    @turbo for i in 1:n
      xb0 = 0.0
      for j in 1:k
        xb0 += x[i,j] * β[j+1]
      end
      xb = xb0 + β[1]
      ti = t[i]
      s += (1 - (c[i] == ti)) * (log(eα) + (eα - 1) * log(ti) + xb) - ti^eα * exp(xb)
    end
    return s
end

function test_logℒ(n, k)
  t = rand(n)
  c = copy(t);
  b = rand(n) .> 0.5
  c[b] .= rand.();
  x = rand(n,k)
  α = 2.85
  β = rand(k+1)
  @test logℒ_fast(α, β, t, c, x) ≈ logℒ_fast_turbo(α, β, t, c, x)
end

function not_an_outer_reduct!(r, N::Int, x = 2.0, y= nothing) # there was a bug where this was classified as one
  @turbo for i ∈ eachindex(r)
    acc = y === nothing ? x : r[i]
    for n ∈ 1:N
      acc += 0
    end
    r[i] = acc
  end
  r
end
@testset "Outer Reductions" begin
  for T ∈ [Float32,Float64,Int32,Int64]
    test_awmean(T)
  end
  @test all(==(7.4), not_an_outer_reduct!(Vector{Float64}(undef, 5), 17, 7.4))
  for n ∈ 1:20, k ∈ 1:5
    test_logℒ(n,k)
  end
  omega = rand(87,87);
  @test reg_term(omega) ≈ reg_term_turbo(omega)

end