text:
#################################################################################################
#                                                                                #
#                              L2 loss - Exact Assignment                                       #
#                                                                                #
#################################################################################################

"""
    function select_closest_center(data::AbstractMatrix, qd::QuantizerData, i::Int, j::Int)

Helper function for L2 partition assignment.\n
Selects the closest center for the given datapoint `i` and codebook index `j`
"""
function select_closest_center(data::AbstractMatrix, qd::QuantizerData, i::Int, j::Int)
    d_i1, d_i2 = (j-1)*qd.n_dims_center+1, j*qd.n_dims_center
    distances = zeros(qd.n_centers)
    for k in 1:qd.n_centers
        c_i = (j-1)*qd.n_centers*qd.n_dims_center + (k-1)*qd.n_dims_center + 1
        distances[k] = sqeuclidean(@view(data[d_i1:d_i2, i]),
                                   @view(qd.C[c_i:(c_i+qd.n_dims_center-1)]))
    end
    choice = argmin(distances)
    local_loss = distances[choice]
    return local_loss, choice
end

"""
    function reset_assignments(qd::QuantizerData, method::UpdateMethod)

Helper function for L2 partition assignment.\n
Resets all cluster/codebook assignments in `qd.I`
"""
function reset_assignments!(qd::QuantizerData, method::BMatrix)
    Threads.@threads for i in 1:qd.n_dp
        qd.I.B[i] = spzeros(Bool, (qd.n_dims_center * qd.n_codebooks * qd.n_centers),qd.n_dims)
    end
    fill!(qd.I.assignments, 0)
end

function reset_assignments!(qd::QuantizerData, method::InvertedIndex)
    qd.I.IVF[:] = [[Vector{Int}() for i in 1:qd.n_centers] for j in 1:qd.n_codebooks]
    fill!(qd.I.assignments, 0)
end

"""
    function assign_cluster(qd::QuantizerData, c_ij::Int, i::Int, j::Int, method::UpdateMethod)

Helper function for L2 partition assignment.\n
Updates the indexes in `qd.I` for the selected method for the given cluster `c_ij`, datapoint `i` and dimension unit (codebook) `j`.
"""
function assign_cluster!(qd::QuantizerData, c_ij::Integer, i::Int, j::Int, method::BMatrix)
    b_1, b_2 = (j-1)*qd.n_dims_center+1, ((qd.n_centers * (j-1) * qd.n_dims_center) + 1) + ((c_ij-1)*qd.n_dims_center)
    qd.I.B[i][b_2: (b_2+qd.n_dims_center-1),b_1:(b_1+qd.n_dims_center-1)] = sparse(I, qd.n_dims_center, qd.n_dims_center)
    qd.I.assignments[j,i] = c_ij
end
function assign_cluster!(qd::QuantizerData, c_ij::Integer, i::Int, j::Int, method::InvertedIndex)
    append!(qd.I.IVF[j][c_ij],i)
    qd.I.assignments[j,i] = c_ij
end

"""
    function append_with_datalock(qd::QuantizerData,data_locks::Array{Bool, 2}, i::Int, j::Int, c_ij::Int)

Helper function for the assignment_loop! function, specifically for the IVF MultiThreaded configuration to
avoid data race while appending.
"""
function append_with_datalock(qd::QuantizerData,data_locks::Array{Bool, 2}, i::Int, j::Int, c_ij::Int)
    task_completed = false
    while !task_completed
        if data_locks[j,c_ij]
            @inbounds data_locks[j,c_ij] = false
            @inbounds append!(qd.I.IVF[j][c_ij],i)
            @inbounds data_locks[j,c_ij] = true
            task_completed=true
        end
    end
end

"""
    function assignment_loop!(data::AbstractMatrix, qd::QuantizerData, method::UpdateMethod, threading::Threading)

Helper function for L2 partition assingment.\n
Loops through all datapoints and all dimension units (codebooks) and assigns the subdimensions to the L2 closest center.
"""
function assignment_loop!(data::AbstractMatrix, qd::QuantizerData, method::BMatrix, threading::MultiThreaded)
    Threads.@threads for i in 1:qd.n_dp
        for j in 1:qd.n_codebooks
            local_loss, c_ij = select_closest_center(data, qd, i, j)
            assign_cluster!(qd, c_ij, i, j, method)
        end
    end
    return nothing
end

function assignment_loop!(data::AbstractMatrix, qd::QuantizerData, method::InvertedIndex, threading::MultiThreaded)
    data_locks = Array{Bool, 2}(undef, qd.n_codebooks, qd.n_centers)
    fill!(data_locks, true)
    Threads.@threads for i in 1:qd.n_dp
        for j in 1:qd.n_codebooks
            _, c_ij = select_closest_center(data, qd, i, j)
            append_with_datalock(qd, data_locks, i, j, c_ij)
            qd.I.assignments[j,i] = c_ij
        end
    end
    return nothing
end

function assignment_loop!(data::AbstractMatrix, qd::QuantizerData, method::UpdateMethod, threading::SingleThreaded)
    total_loss = 0
    for i in 1:qd.n_dp
        for j in 1:qd.n_codebooks
            local_loss, c_ij = select_closest_center(data, qd, i, j)
            assign_cluster!(qd, c_ij, i, j, method)
            total_loss += local_loss
        end
    end
    return total_loss / qd.n_dp
end

"""
    function assignment_step!(data::AbstractMatrix, qd::QuantizerData, method::UpdateMethod, threading::Processing)

Assignment step for L2_loss: after having recomputed / initialized the codebooks,
we find for each dimension unit `u_j` per datapoint `x_i` the favorable (l2-closest) cluster assignment `c_ij`.
"""
function assignment_step!(data::AbstractMatrix, qd::QuantizerData, method::UpdateMethod, threading::Processing)
    reset_assignments!(qd, method)
    loss = assignment_loop!(data, qd, method, threading)
    return loss
end

#################################################################################################
#                                                                                #
#                          Anisotropic loss - Approximate Assignments                           #
#                                                                                #
#################################################################################################

"""
    function rebuild_Bmatrix(qd::QuantizerData)

Helper function of the Anisotropic Assignment Step.\nRebuilds the B matrices (identiy matrices)
from the given assignments.
"""
function rebuild_Bmatrix!(qd::QuantizerData)
    for i in 1:qd.n_dp
        qd.I.B[i] = spzeros(Bool, (qd.n_dims_center * qd.n_codebooks * qd.n_centers),qd.n_dims)
        for j in 1:qd.n_codebooks
            assign_cluster!(qd, qd.I.assignments[j,i], i, j, BMatrix())
        end
    end
end

"""
    function reconstruct_codeword_newcenter(qd::QuantizerData, codeword::AbstractArray, j::Int, k::Int)

Helper function for coordinate descent.\n
Updates the full-dimensional codeword at the cluster location given the dimension unit/codebook `j` and cluster `k`.
"""
function reconstruct_codeword_newcenter(qd::QuantizerData, codeword::AbstractArray, j::Int, k::Int)
    cw_i1, cw_i2 = (j-1)*qd.n_dims_center+1, j*qd.n_dims_center
    C_i1 = ((j-1) * qd.n_dims_center * qd.n_centers) + ((k-1) * qd.n_dims_center) + 1
    C_i2 = ((j-1) * qd.n_dims_center * qd.n_centers) + ((k) * qd.n_dims_center)
    codeword[cw_i1:cw_i2] = @view qd.C[C_i1:C_i2]
    return codeword
end

"""
    function coordinate_descent(qd::QuantizerData, codeword::AbstractArray, i::Int, j::Int, η::AnisotropicWeights)

Approximate assignment helper for non-orthogonal Anisotropic loss funciton. Selects the optimal center for a given
codebook, holding the rest of the codeword constant.
"""
function coordinate_descent(data::AbstractMatrix, qd::QuantizerData, codeword::AbstractArray, i::Int, j::Int, η::AnisotropicWeights)
    d_i1, d_i2 = (j-1)*qd.n_dims_center+1, j*qd.n_dims_center
    distances = zeros(qd.n_centers)
    for k in 1:qd.n_centers
        codeword = reconstruct_codeword_newcenter(qd, codeword, j, k)
        distances[k] = anisotropic_loss(@view(data[:, i]),
                                        codeword,
                                        η)
    end
    choice = argmin(distances)
    local_loss = distances[choice]
    return local_loss, choice
end

"""
    function approx_assignment_dp(data::AbstractMatrix, qd::QuantizerData, i::Integer, max_iter::Int)

Approximate anisotropic loss-based assignment on the data point level. Given a data set,
an index for the datapoint, the up-to-date QuantizerData and the weights for the Anisotropic loss function,
finds the optimal cluster assignment for all dimension units / codebooks and updates the assignments. Returns the loss for that dp.
"""
function approx_assignment_dp!(data::AbstractMatrix, qd::QuantizerData, i::Integer,η::AnisotropicWeights, max_iter::Int)
    count = 0
    dict_old = 0
    codeword = qd.I.B[i]'qd.C
    dp_loss = 0
    while (count < max_iter) && (dict_old != qd.I.assignments[:,i])
        dict_old = deepcopy(qd.I.assignments[:,i])
        for j in 1:qd.n_codebooks
            dp_loss, c_ij = coordinate_descent(data, qd, codeword, i, j, η)
            qd.I.assignments[j,i] = UInt16(c_ij)
            codeword = reconstruct_codeword_newcenter(qd, codeword, j, c_ij)
        end
        count +=1
    end
    return dp_loss
end

"""
`function assignment_step!(data::AbstractMatrix, qd::QuantizerData, η::AnisotropicWeights, max_iter::Int,
    threading::Threading)`

Assignment step for Anisotropic Loss: after having recomputed / initialized the codebooks,
we find for each `u_j ` (per codebook) per datapoint x_i the favorable (closest in terms of anisotropic loss) cluster assignment.
"""
function assignment_step!(data::AbstractMatrix, qd::QuantizerData, η::AnisotropicWeights, max_iter::Int,
                                        threading::SingleThreaded)
    total_loss = 0
    for i in 1:qd.n_dp
        total_loss += approx_assignment_dp!(data, qd, i, η, max_iter)
    end
    rebuild_Bmatrix!(qd)
    return total_loss / qd.n_dp
end

function assignment_step!(data::AbstractMatrix, qd::QuantizerData, η::AnisotropicWeights, max_iter::Int,
                                        threading::MultiThreaded)
    Threads.@threads for i in 1:qd.n_dp
        approx_assignment_dp!(data, qd, i, η, max_iter)
    end
    rebuild_Bmatrix!(qd)
    return
end


text:
using Pkg, InteractiveUtils
ENV["PYTHON"] = ""
Pkg.build("Conda")
Pkg.build("PyCall")
using Conda
Conda.add("pytorch"; channel = "pytorch")
pkg"precompile"

ext: jl

max_stars_repo_name: ericphanson/Attempt-diff-eigmin

max_issues_repo_name: ericphanson/Attempt-diff-eigmin

max_forks_repo_name: ericphanson/Attempt-diff-eigmin


text:
using InteractiveUtils: code_llvm
using Restacker: restack
using Test

"""
    llvm_ir(f, args) :: String

Get LLVM IR of `f(args...)` as a string.
"""
llvm_ir(f, args) = sprint(code_llvm, f, Base.typesof(args...))

nmatches(r, s) = count(_ -> true, eachmatch(r, s))

@noinline function filter_map_restack!(ys, xs)
    ys = restack(ys)
    @inbounds for i in eachindex(ys, xs)
        x = xs[i]
        if -0.5 < x < 0.5
            ys[i] = 2x
        end
    end
end

@testset "filter_map" begin
    ir = llvm_ir(filter_map_restack!, (Float64[], Float64[]))
    @debug "filter_map_restack!" LLVM_IR = Text(ir)
    @test nmatches(r"fmul <[0-9]+ x double>", ir) >= 4
    @test nmatches(r"fcmp [a-z]* <[0-9]+ x double>", ir) >= 4
end


text:
# indexnotation/product.jl
#
# A wrapper to store the contraction (product) of two indexed objects and
# evaluate lazily, i.e. evaluate upon calling `deindexify`.

immutable ProductOfIndexedObjects{IA,IB,CA,CB,OA,OB,TA,TB} <: AbstractIndexedObject
    A::IndexedObject{IA,CA,OA,TA}
    B::IndexedObject{IB,CB,OB,TB}
end

Base.conj(P::ProductOfIndexedObjects) = ProductOfIndexedObjects(conj(P.A),conj(P.B))

Base.eltype(P::ProductOfIndexedObjects) = promote_type(eltype(P.A),eltype(P.B))

*(P::ProductOfIndexedObjects, β::Number) = ProductOfIndexedObjects(β*P.A, P.B)
*(β::Number, P::ProductOfIndexedObjects) = *(P,β)
-(P::ProductOfIndexedObjects) = *(-1, P)

@generated function indices{IA,IB,CA,CB,OA,OB,TA,TB}(::Type{ProductOfIndexedObjects{IA,IB,CA,CB,OA,OB,TA,TB}})
    J = unique2([IA...,IB...])
    J = tuple(J...)
    meta = Expr(:meta, :inline)
    Expr(:block, meta, :($J))
end

@generated function *(A::AbstractIndexedObject,B::AbstractIndexedObject)
    IA = indices(A)
    IB = indices(B)
    IC = symdiff(IA,IB)
    oindA, cindA, oindB, cindB, = contract_indices(IA,IB,IC)
    meta = Expr(:meta,:inline)
    if A <: IndexedObject{IA}
        argA = :A
    else
        IA = IA[vcat(oindA,cindA)]
        IA = tuple(IA...)
        indA = :(Indices{$IA}())
        argA = :(indexify(deindexify(A,$indA),$indA))
    end
    if B <: IndexedObject{IB}
        argB = :B
    else
        IB = IB[vcat(cindB,oindB)]
        IB = tuple(IB...)
        indB = :(Indices{$IB}())
        argB = :(indexify(deindexify(B,$indB),$indB))
    end
    meta = Expr(:meta,:inline)
    Expr(:block,meta,:(ProductOfIndexedObjects($argA,$argB)))
end

@generated function deindexify{IA,IB,CA,CB,IC}(P::ProductOfIndexedObjects{IA,IB,CA,CB}, I::Indices{IC}, T::Type = eltype(P))
    meta = Expr(:meta, :inline)
    oindA, cindA, oindB, cindB, indCinoAB = contract_indices(IA, IB, IC)
    indCinAB = vcat(oindA,length(IA)+oindB)[indCinoAB]
    conjA = Val{CA}
    conjB = Val{CB}
    quote
        $meta
        deindexify!(similar_from_indices(T, $indCinAB, P.A.object, P.B.object, $conjA, $conjB), P, I)
    end
end

@generated function deindexify!{IA,IB,CA,CB,IC}(dst, P::ProductOfIndexedObjects{IA,IB,CA,CB}, ::Indices{IC}, β=0)
    oindA, cindA, oindB, cindB, indCinoAB = contract_indices(IA, IB, IC)
    conjA = Val{CA}
    conjB = Val{CB}
    meta = Expr(:meta,:inline)
    quote
        $meta
        contract!(P.A.α*P.B.α,P.A.object,$conjA,P.B.object,$conjB,β,dst,$oindA,$cindA,$oindB,$cindB,$indCinoAB)
    end
end


text:
using SDPSymmetryReduction
using LinearAlgebra
using Test

@testset "SDPSymmetryReduction.jl" begin

    @test iszero(SDPSymmetryReduction.roundToZero(1e-10))

    M = rand(1:10,10,10)
    @test SDPSymmetryReduction.part(M).n == length(unique(M))

    P1 = Partition(3, [1 2 2; 2 3 3; 2 3 3])
    P2 = Partition(3, [1 1 2; 1 1 2; 1 1 3])
    P3 = Partition(6, [1 2 4; 2 3 5; 2 3 6])
    @test SDPSymmetryReduction.coarsestPart(P1,P2).P == P3.P

    @test SDPSymmetryReduction.part(SDPSymmetryReduction.rndPart(P1)).P == P1.P

    M = rand(10,10)
    @test SDPSymmetryReduction.roundMat(M) ≈ M atol = 1e-4

    A = rand(9,3)
    M = rand(3,3)
    @test maximum(abs.(A \ vec(SDPSymmetryReduction.projectAndRound(M, A; round = false)))) ≈ 0.0 atol=1e-10

    T = M - SDPSymmetryReduction.projectAndRound(M, A; round = false)
    @test all(isapprox.(A*(A\vec(T)) - vec(T), 0, atol=1e-8))

    q = 7
    PG2q = vcat([[0, 0, 1]],
    [[0, 1, b] for b = 0:q-1],
    [[1, a, b] for a = 0:q-1 for b = 0:q-1])
    Adj = [x' * y % q == 0 && x != y for x in PG2q, y in PG2q]

    N = length(PG2q)
    C = ones(N^2)
    A = vcat(vec(Adj)', vec(Matrix{Float64}(I, N, N))')
    b = [0.0, 1.0]

    @test admPartSubspace(C, A, b, true).n == 18

    @test sort(blockDiagonalize(admPartSubspace(C, A, b, true), true).blkSizes) == [2,2,2,2,3]

    @test SDPSymmetryReduction.unSymmetrize(P1).P == Partition(4, [1 3 3; 2 4 4; 2 4 4]).P

    # complex block diagonalization tests
    P = Partition(4,[1 2 3 2; 2 1 2 3; 3 2 1 2; 2 3 2 1])
    @test blockDiagonalize(P; complex = true).blkSizes == [1,1,1]

    function failsBlockDiagonalize()
        try
            blockDiagonalize(P; complex = false).blkSizes == [1,1,1]
        catch
            return true
        end
        return false
    end

    @test failsBlockDiagonalize()
end


text:
using DFControl
set_default_pseudo_dir(:pbesolrel,"/home/ponet/pseudos/pbesolrel/PSEUDOPOTENTIALS/")
set_default_pseudo_dir(:pbesol,"/home/ponet/pseudos/pbesol/PSEUDOPOTENTIALS/")
set_default_pseudo_dir(:sssp,"/home/ponet/pseudos/sssp/")
configure_default_pseudos()
set_default_input(:scf,read_qe_input(DFControl.assets_dir*"inputs/qe/scf.in",run_command="mpirun -np 24 ~/bin/pw.x"))
set_default_input(:bands,read_qe_input(DFControl.assets_dir*"inputs/qe/bands.in",run_command="mpirun -np 24 ~/bin/pw.x"))
set_default_input(:nscf,read_qe_input(DFControl.assets_dir*"inputs/qe/nscf.in",run_command="mpirun -np 24 ~/bin/pw.x"))
default_inputs[:scf]
atoms = get_data(HfO2,"scf",:atomic_positions)
change_atoms!(HfO2,atoms,pseudo_set=:sssp)
HfO2 = create_job("HfO2",phd_dir*"HfO2/NSOC",default_inputs[:scf],default_inputs[:bands],default_inputs[:projwfc],default_inputs[:nscf],server_dir = "HfO2/NSOC")
HfO2 = load_job(phd_dir*"HfO2/NSOC",server_dir="HfO2/NSOC")
pull_outputs(HfO2)
using Plots
bands = read_qe_bands_file(HfO2.local_dir*"/bands.out")
plot(read_qe_bands_file(HfO2.local_dir*"/bands.out"),fermi=read_fermi_from_qe_file(HfO2.local_dir*"scf.out"),ylims=[-2,2],legend=false)
remove_flags!(HfO2,:occupations,:smearing,:degauss)
submit_job(HfO2)
change_flow!(HfO2,[("nscf",false)])

using DFControl
hfo2    = load_server_job("HfO2/NSOC/",phd_dir*"HfO2/NSOC")
atoms = get_data(hfo2,"bands",:atomic_positions)
new_atoms = Dict{Symbol,Array{Point3D{Float32},1}}(:O=>[atoms[:O1],atoms[:O2]],:Hf=>[atoms[:Hf]])
atoms[:O] = [atoms[:O1]...,atoms[:O2]...]
change_atoms!(hfo2,atoms,pseudo_set=:pbesol,pseudo_fuzzy="paw")
pop!(atoms,:O2)
print_flow(hfo2)
default_pseudos[:Hf]
atoms
print_data(hfo2,:atomic_species)
outputs = pull_outputs(hfo2)
using Plots
plot(read_qe_bands_file(outputs[2]),fermi=-1.1936,ylims=[-5,10])
print_flags(hfo2)
add_flags!(hfo2,:system,:nbnd=>40)
print_data(hfo2,"scf")
remove_flags!(hfo2,:occupations,:smearing,:degauss)
submit_job(hfo2)
hfo2.header


text:
module Minesweeper
    using StatsBase
    import Base: show
    export Game, Cell, add_mines!, get_neighbors, reveal_zeros!, game_over
    export get_neighbor_indices, show, reveal, play, compute_score!
    export update_reveal!, select_cell!, initilize_cells, mine_count!, import_gui
    export flag_button!, update_flag_count!, flag_cell!
    export generate_gui, reveal_button!, start_new_game!
    include("Game.jl")
end


text:
# Test sample functions

using Distributions
using Base.Test

function est_p(x, K)
	h = zeros(Int, K)
	for xi in x
		h[xi] += 1
	end
	p = h / length(x)
end

#### sample with replacement

n = 10^5
x = sample([10,20,30], n)
@test isa(x, Vector{Int})
@test length(x) == n

h = [sum(x .== 10), sum(x .== 20), sum(x .== 30)]
@test sum(h) == n
ph = h / n
p0 = fill(1/3, 3)
@test_approx_eq_eps ph p0 0.02

#### sample without replacement

# case: K == 2

n = 10^5
x = zeros(Int, 2, n)
for i = 1:n
	v = sample(11:15, 2; replace=false)
	@assert v[1] != v[2]
	x[:,i] = v
end
@test min(x) == 11
@test max(x) == 15

x[:] -= 10  # brings x to 1:5
p0 = fill(1/5, 5)
@test_approx_eq_eps est_p(x, 5) p0 0.02

# case: K == 4 with moderate a (using Fisher-Yates)

n = 10^5
x = zeros(Int, 4, n)
for i = 1 : n
	v = sample(11:20, 4; replace=false)
	sv = sort(v)
	@assert sv[1] < sv[2] < sv[3] < sv[4]
	x[:,i] = v
end
@test min(x) == 11
@test max(x) == 20

x[:] -= 10
p0 = fill(0.1, 10)
@test_approx_eq_eps est_p(x, 10) p0 0.01

# case: K == 4 with very large a (using self-avoid)

n = 10^4
x = zeros(Int, 4, n)
a = 10^7 + 1
b = 2 * 10^7
for i = 1 : n
	v = sample(a:b, 4; replace=false)
	sv = sort(v)
	@assert sv[1] < sv[2] < sv[3] < sv[4]
	x[:,i] = v
end
@test min(x) >= a
@test max(x) <= b

#### weighted sampling

w = [2., 5., 3.]
n = 10^5
x = wsample([10,20,30], w, n)

h = [sum(x .== 10), sum(x .== 20), sum(x .== 30)]
@test sum(h) == n
p0 = w / sum(w)
ph = h / n
@test_approx_eq_eps ph p0 0.02


text:
#FIXME: PyObject not found when pyscf only
import PySCF: scf_rhf #, pyscf

# Workaround: repeat the imports.
using PyCall: pyimport
pyscf = pyimport("pyscf")

@testset "SCF RHF: H₂O" begin
   mol = pyscf.gto.Mole()
   mol.build(verbose = 0,
             basis = "sto-3g",
             atom = """8 0  0.    0
                       1 0 -0.757 0.587
                       1 0  0.757 0.587""",
            )

   mf = pyscf.scf.RHF(mol)

   rhf         = scf_rhf(mol,verbose=false)
   e_rhf_pyscf = mf.kernel()

   @test isapprox(e_rhf_pyscf, rhf.Emol, atol=1e-10)
end


text:
# This is an auto-generated file; do not edit

# Pre-hooks

# Macro to load a library
macro checked_lib(libname, path)
    path = string(path)
    ((VERSION >= v"0.4.0-dev+3844" ? Base.Libdl.dlopen_e : Base.dlopen_e)(path) == C_NULL) && error("Unable to load \n\n$libname ($path)\n\nPlease re-run Pkg.build(package), and restart Julia.")
    quote const $(esc(libname)) = $path end
end

include("constants.jl")

const csdp = dlpath

info("Using $csdp")

# Load-hooks


text:
function displacement(obj::StationDisplacement, arg0::BodiesElements, arg1::Frame, arg2::Vector3D)
    return jcall(obj, "displacement", Vector3D, (BodiesElements, Frame, Vector3D), arg0, arg1, arg2)
end


text:
using Documenter, PkgTemplates

makedocs(;
    modules=[PkgTemplates],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Package Generation" => "pages/package_generation.md",
        "Plugins" => "pages/plugins.md",
        "Plugin Development" => "pages/plugin_development.md",
        "Licenses" => "pages/licenses.md",
        "Index" => "pages/index.md",
    ],
    repo="https://github.com/invenia/PkgTemplates.jl/blob/{commit}{path}#L{line}",
    sitename="PkgTemplates.jl",
    authors="Chris de Graaf, Invenia Technical Computing Corporation",
    assets=[],
)

deploydocs(;
    repo="github.com/invenia/PkgTemplates.jl",
)


text:
xml_evfile1 = path*"/SampleFiles/XML/fdsnws-event_2017-01-12T03-18-55Z.xml"
xml_evfile2 = path*"/SampleFiles/XML/ISC_2011-tohoku-oki.xml"
xml_stfile  = path*"/SampleFiles/XML/fdsnws-station_2017-01-12T03-17-42Z.xml"
d1 = "2019-03-14T02:18:00"
d2 = "2019-03-14T02:28:00"
to = 30

# fdsn_hdr
# @test fdsn_hdr("http://http://service.ncedc.org/fdsnws/station/1/") == ["Host" => "service.ncedc.org", "User-Agent" => "curl/7.60.0", "Accept" => "*/*"]
# @test fdsn_hdr("http://http://service.ncedc.org/fdsnws/dataselect/1/") == ["Host" => "service.ncedc.org", "User-Agent" => "curl/7.60.0", "Accept" => "*/*"]
# @test fdsn_hdr("http://http://service.ncedc.org/fdsnws/event/1/") == ["Host" => "service.ncedc.org", "User-Agent" => "curl/7.60.0", "Accept" => "*/*"]
# @test fdsn_hdr("http://service.scedc.caltech.edu/fdsnws/dataselect/1/")  == ["User-Agent" => "Julia"]
# @test fdsn_hdr("http://service.scedc.caltech.edu/fdsnws/station/1/")  == ["User-Agent" => "Julia"]
# @test fdsn_hdr("http://service.scedc.caltech.edu/fdsnws/event/1/")  == ["User-Agent" => "Julia"]
# @test fdsn_hdr("http://www.cnn.com.") == webhdr

# From FDSN the response code is 200
url = "http://service.iris.edu/fdsnws/dataselect/1/query?format=miniseed&net=UW&sta=XNXNX&loc=99&cha=QQQ&start="*d1*"&end="*d2*"&szsrecs=true"
req_info_str = datareq_summ("FDSNWS data", "UW.XNXNX.99.QQQ", d1, d2)
(req, parsable) = get_http_req(url, req_info_str, to, status_exception=false)
@test typeof(req) == Array{UInt8,1}
@test startswith(String(req), "HTTP.Messages.Response")
@test parsable == false
(req,parsable) = get_http_req(url, req_info_str, to, status_exception=true)
@test typeof(req) == Array{UInt8,1}
@test startswith(String(req), "HTTP.Messages.Response")
@test parsable == false

# From IRIS the response code is 400 and we can get an error
url = "http://service.iris.edu/irisws/timeseries/1/query?net=DE&sta=NENA&loc=99&cha=LUFTBALLOONS&start="*d1*"&end="*d2*"&scale=AUTO&format=miniseed"
req_info_str = datareq_summ("IRISWS data", "DE.NENA.99.LUFTBALLOONS", d1, d2)
(req,parsable) = get_http_req(url, req_info_str, to, status_exception=false)
@test typeof(req) == Array{UInt8,1}
@test startswith(String(req), "HTTP.Messages.Response")
@test parsable == false
(req, parsable) = get_http_req(url, req_info_str, to, status_exception=true)
@test typeof(req) == Array{UInt8,1}
@test startswith(String(req), "HTTP.Messages.Response")
@test parsable == false
(req, parsable) = get_http_post(url, NOOF, to, status_exception=false)
@test typeof(req) == Array{UInt8,1}
@test startswith(String(req), "HTTP.Messages.Response")
@test parsable == false
(req, parsable) = get_http_post(url, NOOF, to, status_exception=true)
@test typeof(req) == Array{UInt8,1}
@test startswith(String(req), "HTTP.Messages.Response")
@test parsable == false

printstyled("  FDSN XML\n", color=:light_green)

id_err = "error in Station ID creation!"
unit_err = "units don't match instrument code!"
true_id = String["3337497", "3279407", "2844986", "2559759", "2092067", "1916079", "2413"]
true_ot = DateTime("2011-03-11T05:46:23.200")
true_loc = Float64[2.2376 38.2963; 93.0144 142.498; 26.3 19.7]
true_mag = Float32[8.6, 9.1, 8.8, 8.5, 8.6, 9.0, 8.5]
true_msc = String["MW", "MW", "MW", "MW", "MW", "MW", ""]
r1 = PZResp(f0 = 0.02f0, p = ComplexF32[-981.0+1009.0im, -981.0-1009.0im, -3290.0+1263.0im, -3290.0-1263.0im])
resp_a0!(r1)
r2 = PZResp(Complex{Float32}.([   0.0+0.0im       -0.037-0.037im
                                  0.0+0.0im       -0.037+0.037im
                                  -15.15+0.0im    -15.64+0.0im
                                  -176.6+0.0im    -97.34-400.7im
                                  -463.1-430.5im  -97.34+400.7im
                                  -463.1+430.5im  -374.8+0.0im
                                  0.0+0.0im       -520.3+0.0im
                                  0.0+0.0im       -10530.0-10050.0im
                                  0.0+0.0im       -10530.0+10050.0im
                                  0.0+0.0im       -13300.0+0.0im
                                  0.0+0.0im       -255.097+0.0im ]),rev=true)
r2.z = r2.z[1:6]
r2.f0 = 0.02f0
resp_a0!(r2)

printstyled("    QuakeML test 1\n", color=:light_green)
(EC,RR) = read_qml(xml_evfile1)
Nev = length(EC)
@test Nev == length(true_id)
for i = 1:Nev
  @test EC[i].id == true_id[i]
  @test EC[i].mag.val == true_mag[i]
  @test EC[i].mag.scale == true_msc[i]
end

@test EC[2].ot==true_ot
for i = 1:2
  @test ≈(EC[i].loc.lat, true_loc[1,i])
  @test ≈(EC[i].loc.lon, true_loc[2,i])
  @test ≈(EC[i].loc.dep, true_loc[3,i])
end

printstyled("    QuakeML test 2\n", color=:light_green)
H, R = read_qml(xml_evfile2)
H = H[1]
R = R[1]

# Check basic headers
@test H.typ == "earthquake"
@test H.id == "16461282"

# Check that the correct magnitude is retained
@test H.mag.val ≥ 9.0f0
@test H.mag.scale == "MW"

# Check H.loc
@test H.loc.lat == 38.2963
@test H.loc.lon == 142.498
@test H.loc.dep == 19.7152
@test H.loc.se == 2.1567
@test H.loc.nst == 2643
@test H.loc.src == "smi:ISC/origid=602227159,ISC"

# Check source params
@test R.id == "600002952"
@test R.m0 == 5.312e22
@test R.mt == [1.73e22, -2.81e21, -1.45e22, 2.12e22, 4.55e22, -6.57e21]
@test R.dm == [6.0e19, 5.0e19, 5.0e19, 6.8e20, 6.5e20, 4.0e19]
@test R.pax == [295.0 115.0 205.0; 55.0 35.0 0.0; 5.305e22 -5.319e22 1.4e20]
@test R.planes == [25.0 203.0; 80.0 10.0; 90.0 88.0]
@test R.st.dur == 70.0
@test R.misc["methodID"] == "smi:ISC/methodID=Best_double_couple"
@test R.misc["pax_desc"] == "azimuth, plunge, length"
@test R.misc["planes_desc"] == "strike, dip, rake"
@test R.misc["derivedOriginID"] == "600126955"

printstyled("    station XML\n", color=:light_green)
io = open(xml_stfile, "r")
xsta = read(io, String)
close(io)

web_chanspec()
web_chanspec()
web_chanspec()

# Test tracking of changes in SeisData structure S ==========================
# Does S change? Let's see.
printstyled("  Tracking with track_on!, track_off!\n", color=:light_green)

S = SeisData()
track_on!(S)  # should do nothing but no error
@test track_off!(S) == nothing # should return nothing for empty struct

S = SeisData(3)
@test track_off!(S) == [true,true,true]
# @test_throws ErrorException("Tracking not enabled!") track_off!(S)

# Now replace S with randSeisData
S = randSeisData(3)
track_on!(S)
@test haskey(S.misc[1], "track")
push!(S, randSeisChannel())
u = track_off!(S)
@test (u == [false, false, false, true])
@test haskey(S.misc[1], "track") == false

# Now turn tracking on again and move things around
track_on!(S)
@test haskey(S.misc[1], "track")
Ch1 = deepcopy(S[1])
Ch3 = deepcopy(S[3])
S[3] = deepcopy(Ch1)
S[1] = deepcopy(Ch3)
@test haskey(S.misc[3], "track")
@test haskey(S.misc[1], "track") == false
@test haskey(S.misc[2], "track") == false
append!(S.x[1], rand(Float64, 1024))        # Should flag channel 1 as updated
S.id[2] = reverse(S.id[2])                  # Should flag channel 2 as updated
u = track_off!(S)
@test (u == [true, true, false, false])
@test haskey(S.misc[3], "track") == false


text:
"""
# PSFModels

Statistical models for constructing point-spread functions (PSFs).

## Models

The following models are currently implemented
* [`gaussian`](@ref)/[`normal`](@ref)
* [`airydisk`](@ref)
* [`moffat`](@ref)

## Parameters

In general, the PSFs have a position, a full-width at half-maximum (FWHM) measure, and an amplitude. The position follows a 1-based pixel coordinate system, where `(1, 1)` represents the *center* of the bottom left pixel. This matches the indexing style of Julia as well as DS9, IRAF, SourceExtractor, and WCS. The FWHM is a consistent scale parameter for the models. That means a `gaussian` with a FWHM of 5 will be visually similar to an `airydisk` with a FWHM of 5. All models support a scalar (isotropic) FWHM and a FWHM for each axis (diagonal), as well as arbitrarily rotating the PSF.

!!! warning "Pixel-convention"

    The pixel convention adopted here is that the bottom-left pixel *center* is `(1, 1)`. The column-major memory layout of julia puts the `x` axis as the rows of a matrix and the `y` axis as the columns. In other words, the axes unpack like
    ```julia
    xs, ys = axes(image)
    ```

## Usage

### Evaluating models

Directly evaluating the functions is the most straightforward way to use this package

```jldoctest model
julia> gaussian(0, 0; x=0, y=0, fwhm=3)
1.0

julia> gaussian(BigFloat, 0, 0; x=0, y=0, fwhm=3, amp=0.1, bkg=1)
1.100000000000000088817841970012523233890533447265625
```

We also provide "curried" versions of the functions, which allow you to specify the parameters and evaluate the PSF later

```jldoctest model
julia> model = gaussian(x=0, y=0, fwhm=3);

julia> model(0, 0)
1.0
```

If we want to collect the model into a dense matrix, simply iterate over indices

```jldoctest model
julia> inds = CartesianIndices((-2:2, -2:2));

julia> model.(inds) # broadcasting
5×5 Matrix{Float64}:
 0.0850494  0.214311  0.291632  0.214311  0.0850494
 0.214311   0.54003   0.734867  0.54003   0.214311
 0.291632   0.734867  1.0       0.734867  0.291632
 0.214311   0.54003   0.734867  0.54003   0.214311
 0.0850494  0.214311  0.291632  0.214311  0.0850494
```

This makes it very easy to evaluate the PSF on the same axes as an image (array)

```jldoctest model
julia> img = randn(5, 5);

julia> model.(CartesianIndices(img))
5×5 Matrix{Float64}:
 0.54003      0.214311     0.0459292    0.00531559   0.000332224
 0.214311     0.0850494    0.018227     0.00210949   0.000131843
 0.0459292    0.018227     0.00390625   0.000452087  2.82555e-5
 0.00531559   0.00210949   0.000452087  5.2322e-5    3.27013e-6
 0.000332224  0.000131843  2.82555e-5   3.27013e-6   2.04383e-7
```

this is trivially expanded to fit "stamps" in images

```jldoctest model
julia> big_img = randn(1000, 1000);

julia> stamp_inds = (750:830, 400:485);

julia> stamp = @view big_img[stamp_inds...];

julia> stamp_model = model.(CartesianIndices(stamp_inds));
```

or we can create a loss function for fitting PSFs without allocating any memory. We are simply iterating over the image array!

```jldoctest model
julia> using Statistics

julia> mse = mean(I -> (big_img[I] - model(I))^2, CartesianIndices(stamp_inds));
```

### Fitting data

There exists a simple, yet powerful, API for fitting data with [`PSFModels.fit`](@ref).

```julia
# `fit` is not exported to avoid namespace clashes
using PSFModels: fit

data = # load data
stamp_inds = # optionally choose indices to "cutout"

# use an isotropic Gaussian
params, synthpsf = fit(gaussian, (x=12, y=13, fwhm=3.2, amp=0.1),
                       data, stamp_inds)
# elliptical, rotated Gaussian
params, synthpsf = fit(gaussian, (x=12, y=13, fwhm=(3.2, 3.2), amp=0.1, theta=0)
                       data, stamp_inds)
# obscured Airy disk
params, synthpsf = fit(airydisk, (x=12, y=13, fwhm=3.2, amp=0.1, ratio=0.3),
                       data, stamp_inds)
# bivariate Moffat with arbitrary alpha
params, synthpsf = fit(moffat, (x=12, y=13, fwhm=(3.2, 3.2), amp=0.1, alpha=1),
                       data, stamp_inds)
```

### Plotting

Finally, we provide plotting recipes (`psfplot`/`psfplot!`) from [RecipesBase.jl](https://github.com/JuliaPlots/RecipesBase.jl), which can be seen in use in the [API/Reference](@ref) section.

```julia
using Plots

model = gaussian(x=0, y=0, fwhm=(8, 10), theta=12)
psfplot(model, -30:30, -30:30, colorbar_scale=:log10)
```
"""
module PSFModels

using CoordinateTransformations
using KeywordCalls
using Optim
using Rotations
using SpecialFunctions
using StaticArrays
using Statistics

export gaussian, normal, airydisk, moffat

const BivariateLike = Union{<:Tuple,<:AbstractVector}

function rotate_point(dx, dy, theta)
    # generate rotation matrix
    # (theta is degrees CCW from x-axis)
    R = RotMatrix{2}(-deg2rad(theta))
    # rotate points
    return R * SA[dx, dy]
end

include("gaussian.jl")
include("airy.jl")
include("moffat.jl")
include("plotting.jl")

# codegen for common functionality
# if you add a new model, make sure it gets added here
for model in (:gaussian, :airydisk, :moffat)
    @eval begin
        $model(px, py; kwargs...) = $model(Float64, px, py; kwargs...)
        $model(point::BivariateLike; kwargs...) = $model(point...; kwargs...)
        $model(T, point::BivariateLike; kwargs...) = $model(T, point...; kwargs...)
        $model(idx::CartesianIndex; kwargs...) = $model(idx.I; kwargs...)
        $model(T, idx::CartesianIndex; kwargs...) = $model(T, idx.I; kwargs...)
        $model(; kwargs...) = (point...) -> $model(_curried_point(point...); kwargs...)
        $model(::Type{T}; kwargs...) where {T} = (point...) -> $model(T, _curried_point(point...); kwargs...)
    end
end

_curried_point(P::BivariateLike) = P
_curried_point(point...) = Tuple(point)

include("fitting.jl")

end # module PSFModels


text:
# Copyright 2021 The CLVR Authors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

# File containing helper functions for algorithms.

""" Compute the nonzero rows in A_T for each block."""
function compute_nzrows_for_blocks(A_T::SparseMatrixCSC, blocksize::Int)

    _, n = size(A_T)

    blocks = Array{UnitRange{Int}}([])
    C = Array{Vector{Int}}([])
    len_b = n ÷ blocksize
    for i = 1:len_b
        if i == len_b
            push!(blocks, (1 + (i-1) * blocksize): n)
        else
            push!(blocks, (1 + (i-1) * blocksize): (i * blocksize))
        end

        row_set = Set{Int}()
        for j in blocks[i]
            loc = A_T.colptr[j]:(A_T.colptr[j+1]-1)
            union!(row_set, A_T.rowval[loc])
        end
        row_vec = collect(row_set)
        push!(C, row_vec)
    end

    blocks, C
end

"""Export results into CSV formatted file."""
function exportresultstoCSV(results::Results, outputfile::String)

    CSV.write(
        outputfile, (
            iterations = results.iterations,
            times = results.times,
            fvaluegaps = results.fvaluegaps,
            metricLPs = results.metricLPs
        )
    )
end

"""Computing a common metric for LP. See Eqn (20) in Applegate et al 2020."""
function compute_fvaluegap_metricLP(x_out::Vector, y_out::Vector, problem::StandardLinearProgram)

    A_T, b, c = problem.A_T, problem.b, problem.c

    norm1 = norm(max.(-x_out, 0))
    norm2 = norm(max.((x_out' * A_T)' - b, 0))
    norm3 = norm(max.(-(x_out' * A_T)' + b, 0))  # TODO: Can combine with norm2
    norm4 = norm(max.(-A_T * y_out - c, 0))
    norm5 = norm(c' * x_out + b' * y_out)

    norm5, sqrt(norm1^2 + norm2^2 + norm3^2 + norm4^2 + norm5^2)
end

"""Compute number of nonzero elements of each row in a sparse column matrix."""
function nnz_per_row(A_T::SparseMatrixCSC)

    ret = zeros(Int, size(A_T)[1])
    rows, _, values = findnz(A_T)
    for (r, v) in zip(rows, values)
        if !isapprox(v, 0.)
            ret[r] += 1
        end
    end
    return ret
end

max_stars_repo_name: ericlincc/Efficient-GLP

max_issues_repo_name: ericlincc/Efficient-GLP

max_forks_repo_name: ericlincc/Efficient-GLP


text:
"""
    blocklossfn(predblock, yblock)

Construct a loss function that compares a batch of model outputs
`ŷs` and encoded targets `ys` and returns a scalar loss.

For example for `block = OneHotLabel(classes)` (i.e. an encoded
`Label(classes)`), we have
`blocklossfn(block, block) == Flux.Losses.logitcrossentropy`.
"""
function blocklossfn end

blocklossfn(predblock) = blocklossfn(predblock, predblock)


text:
#! format: off
@deprecate imfill(img::AbstractArray{Bool}, interval::Tuple{Real,Real}, dims::Union{Dims, AbstractVector{Int}}) imfill(img, interval; dims=dims)

@deprecate dilate!(img; kwargs...) dilate!(img, copy(img); kwargs...)
@deprecate erode!(img; kwargs...) erode!(img, copy(img); kwargs...)

@deprecate opening!(img; kwargs...) opening!(img, copy(img), similar(img); kwargs...)
@deprecate closing!(img; kwargs...) closing!(img, copy(img), similar(img); kwargs...)

Base.@deprecate_binding morphogradient mgradient
Base.@deprecate_binding morpholaplace mlaplacian

@deprecate extremefilt!(A, select; kwargs...) extreme_filter!(select, A, copy(A); kwargs...)

import .FeatureTransform: feature_transform
@deprecate feature_transform(img, weights; kwargs...) feature_transform(img; weights=weights, kwargs...)


text:
__precompile__()

module SortedVectors

import Base: getindex, length, push!, isempty,
        pop!, popfirst!

export SortedVector

import ..StrategyBase:filter_elements!
using ..StrategyBase

struct SortedVector{T, F<:Function} <: Strategy

    data::Vector{T}
    by::F

    function SortedVector(data::Vector{T}, by::F) where {T,F}
        new{T,F}(sort(data, by = by), by)
    end
end

SortedVector(data::Vector{T}) where {T} = SortedVector(data, identity)

function show(io::IO, v::SortedVector)
    print(io, "SortedVector($(v.data))")
end

getindex(v::SortedVector, i::Int) = v.data[i]
length(v::SortedVector) = length(v.data)

function push!(v::SortedVector{T}, x::T) where {T}
    i = searchsortedfirst(v.data, x, by=v.by)
    insert!(v.data, i, x)
    return v
end

isempty(v::SortedVector) = isempty(v.data)

pop!(v::SortedVector) = pop!(v.data)

popfirst!(v::SortedVector) = popfirst!(v.data)

function filter_elements!(v::SortedVector{T}, x::T) where {T}
    cutoff = searchsortedfirst(v.data, x, by=v.by)
    if cutoff >= length(v.data)
        return v
    end
    resize!(v.data, cutoff)
    return v
end

end


text:
using CausalityToolsBase
using StateSpaceReconstruction

# You'd use an
pts = invariantize([rand(3) for i = 1:200])
spts = [SVector{3, Float64}(pt) for pt in pts]
mpts = [MVector{3, Float64}(pt) for pt in pts]
D = Dataset(pts)

# Initialising test
@test VisitationFrequency().b isa Number
@test VisitationFrequency(b = 2).b isa Number
@test VisitationFrequency(b = 10).b isa Number
@test TransferOperatorGrid().b isa Number
@test TransferOperatorGrid(b = 10).b isa Number

vars = TEVars([1], [2], [3])
@test transferentropy(pts, vars, RectangularBinning(0.2), VisitationFrequency()) >= 0
@test transferentropy(spts, vars, RectangularBinning(0.2), VisitationFrequency()) >= 0
@test transferentropy(mpts, vars,  RectangularBinning(0.2), VisitationFrequency()) >= 0
@test transferentropy(D, vars, RectangularBinning(0.2), VisitationFrequency()) >= 0

# Different logarithms should give different TE values
te_b2 = transferentropy(pts, vars, RectangularBinning(0.2), VisitationFrequency(b = 2))
te_b10 = transferentropy(pts, vars, RectangularBinning(0.2), VisitationFrequency(b = 10))
@test te_b2 != te_b10

@test transferentropy(pts, vars, RectangularBinning(0.2), TransferOperatorGrid()) >= 0
@test transferentropy(spts, vars, RectangularBinning(0.2), TransferOperatorGrid()) >= 0
@test transferentropy(mpts, vars, RectangularBinning(0.2), TransferOperatorGrid()) >= 0
@test transferentropy(D, vars, RectangularBinning(0.2), TransferOperatorGrid()) >= 0

@test NearestNeighbourMI(b = 10).b isa Number

pts = invariantize([rand(3) for i = 1:15])
spts = [SVector{3, Float64}(pt) for pt in pts]
mpts = [MVector{3, Float64}(pt) for pt in pts]
D = Dataset(pts)

te_D = transferentropy(D, vars, NearestNeighbourMI(b = 2))
te_pts = transferentropy(pts, vars, NearestNeighbourMI(b = 2))
te_spts = transferentropy(spts, vars, NearestNeighbourMI(b = 2))
te_mpts = transferentropy(mpts, vars, NearestNeighbourMI(b = 2))

@test te_D isa Number
@test te_pts isa Number
@test te_spts isa Number
@test te_mpts isa Number

te_D ≈ te_pts ≈ te_spts ≈ te_mpts

# Compute invariant measure over a triangulation using approximate and exact
# simplex intersections. This is relatively slow.
μapprox = invariantmeasure(pts, TriangulationBinning(), ApproximateIntersection())
μexact = invariantmeasure(pts, TriangulationBinning(), ExactIntersection())

# Compute transfer entropy from the invariant measure over multiple
# bin sizes. This is fast, because the measure has been precomputed.
@test transferentropy(μapprox, vars, RectangularBinning(0.2)) >= 0
@test transferentropy(μexact, vars, RectangularBinning(0.2)) >= 0

#@test transferentropy_transferoperator_grid(pts, 0.3, vars), transferentropy(pts,  vars, RectangularBinning(0.3),TransferOperatorGrid())
#@test transferentropy_visitfreq(pts, 0.3, vars), transferentropy(pts, vars, RectangularBinning(0.3), VisitationFrequency())


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ 36d2a4e8-f81f-11eb-3c33-bb29b29df2fb
using Viznet,Compose,GeometryBasics,LinearAlgebra

# ╔═╡ 6b17ade8-171f-4e56-9632-97691f9f5fbd
# define you atoms/lattice
begin
	# OBC here
	irange,jrange,krange = 0.05:0.2:0.9,0.05:0.2:0.9,0.3:0.4:0.9
	vertices = [(i,j,k) for i = irange for j = jrange for k = krange]
	edges_h = [((i,j,k),(i+step(irange),j,k)) for i=irange[1:end-1] for j=jrange  for k = krange]
	edges_v = [((i,j,k),(i,j+step(jrange),k)) for i=irange for j=jrange[1:end-1]  for k = krange]
	edges_z = [((i,j,k),(i,j,k+step(krange))) for i=irange for j=jrange for k = krange[1:end-1]]
	edges = vcat(edges_h,edges_v,edges_z)
	graph = (vertices,edges)
end

# ╔═╡ aef35114-62b0-406c-9dc8-ecfb2a36b24b
# set the line brush and node brush
begin
	linebrush = bondstyle(:default, stroke("grey"))
	nodebrush = nodestyle(:default, stroke("lightblue"),fill("lightblue"),r=0.011)
end;

# ╔═╡ 027510e1-7b89-410b-a8b7-79cb1fc7ffe3
function PerspectiveMap3D(kx::Real,ky::Real)
	pmap(p) = (p[1]+kx*p[3],p[2]+ky*p[3])
	return pmap
end

# ╔═╡ ca796e3b-6ffc-4dbc-9e67-163c0dd9a1aa
pfv = PerspectiveMap3D(0.19,0.15)

# ╔═╡ 8d981c85-325f-4062-a9f4-cae7b795beec
pfe(line) = (pfv(line[1]),pfv(line[2]))

# ╔═╡ ace074e5-5227-48ed-9bf9-975d0d10b170
begin
	vertices2D = map(pfv,vertices)
	edges2D = map(pfe,edges)
	graph2D = (vertices2D,edges2D)
end

# ╔═╡ c0b772d5-8345-4dc4-8801-45d4e614206d
# draw something on the canvas
canvas() do
    for i in vertices2D
        nodebrush >> i
    end
    for (i, j) in edges2D
        linebrush >> (i,j)
    end
end

# ╔═╡ 0e253811-3d8d-4c3e-a9d7-0e7e071ce7a1
function GeoGraph(graph2D)
	points = Point.(graph2D[1])
	lines = map(x->Line(map(Point,x)...),graph2D[2])
	return (points,lines)
end

# ╔═╡ 9f33023b-cb6a-4271-a2c8-5406e8cd18fb
function intersects2D(geograph2D;eps=1E-10)
	points = geograph2D[1]
	lines = geograph2D[2]
	nedges = length(lines)
	buffer = []
	for i = 1:nedges-1
		for j = i+1:nedges
			isinter,point = intersects(lines[i],lines[j])
			if isinter
				tp = (lines[i][1],lines[i][2],lines[j][1],lines[j][2])
				if norm(tp[1]-point)>eps && norm(tp[2]-point)>eps && norm(tp[3]-point)>eps && norm(tp[4]-point)>eps
					push!(buffer,point)
				end
			end
		end
	end
	return buffer
end

# ╔═╡ 6a74d1df-6e26-4fbe-bec7-242292910180
fgraph = GeoGraph(graph2D)

# ╔═╡ dc216f7c-ee56-4510-aa05-891af745da03
gates = intersects2D(fgraph)

# ╔═╡ 4fd816a7-1c9e-46d5-ba80-62914ff54998
canvas() do
    for i in vertices2D
        nodebrush >> i
    end
    for (i, j) in edges2D
        linebrush >> (i,j)
    end
	gatesbrush = nodestyle(:square, stroke("black");r=0.005)
	for i in gates
		gatesbrush >> Tuple(i)
	end
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Compose = "a81c6b42-2e10-5240-aca2-a61377ecd94b"
GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Viznet = "52a3aca4-6234-47fd-b74a-806bdf78ede9"

[compat]
Compose = "~0.9.2"
GeometryBasics = "~0.4.1"
Viznet = "~0.3.3"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.1"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "344f143fa0ec67e47917848795ab19c6a455f32c"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.32.0"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[Compose]]
deps = ["Base64", "Colors", "DataStructures", "Dates", "IterTools", "JSON", "LinearAlgebra", "Measures", "Printf", "Random", "Requires", "Statistics", "UUIDs"]
git-tree-sha1 = "c6461fc7c35a4bb8d00905df7adafcff1fe3a6bc"
uuid = "a81c6b42-2e10-5240-aca2-a61377ecd94b"
version = "0.9.2"

[[DataAPI]]
git-tree-sha1 = "ee400abb2298bd13bfc3df1c412ed228061a2385"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.7.0"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "4437b64df1e0adccc3e5d1adbc3ac741095e4677"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.9"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Dierckx]]
deps = ["Dierckx_jll"]
git-tree-sha1 = "5fefbe52e9a6e55b8f87cb89352d469bd3a3a090"
uuid = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
version = "0.5.1"

[[Dierckx_jll]]
deps = ["CompilerSupportLibraries_jll", "Libdl", "Pkg"]
git-tree-sha1 = "a580560f526f6fc6973e8bad2b036514a4e3b013"
uuid = "cd4c43a9-7502-52ba-aa6d-59fb2a88580b"
version = "0.0.1+0"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[EarCut_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "92d8f9f208637e8d2d28c664051a00569c01493d"
uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
version = "2.1.5+1"

[[FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[GeometryBasics]]
deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
git-tree-sha1 = "58bcdf5ebc057b085e58d95c138725628dd7453c"
uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
version = "0.4.1"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[IterTools]]
git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.3.0"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Measures]]
git-tree-sha1 = "e498ddeee6f9fdb4551ce855a46f54dbd900245f"
uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
version = "0.3.1"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "477bf42b4d1496b454c10cce46645bb5b8a0cf2c"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.2"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3240808c6d463ac46f1c1cd7638375cd22abbccb"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.12"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StructArrays]]
deps = ["Adapt", "DataAPI", "StaticArrays", "Tables"]
git-tree-sha1 = "000e168f5cc9aded17b6999a560b7c11dda69095"
uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
version = "0.6.0"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "d0c690d37c73aeb5ca063056283fde5585a41710"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.5.0"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Viznet]]
deps = ["Compose", "Dierckx"]
git-tree-sha1 = "7a022ae6ac8b153d47617ed8c196ce60645689f1"
uuid = "52a3aca4-6234-47fd-b74a-806bdf78ede9"
version = "0.3.3"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═36d2a4e8-f81f-11eb-3c33-bb29b29df2fb
# ╠═6b17ade8-171f-4e56-9632-97691f9f5fbd
# ╠═aef35114-62b0-406c-9dc8-ecfb2a36b24b
# ╠═027510e1-7b89-410b-a8b7-79cb1fc7ffe3
# ╠═ca796e3b-6ffc-4dbc-9e67-163c0dd9a1aa
# ╠═8d981c85-325f-4062-a9f4-cae7b795beec
# ╠═ace074e5-5227-48ed-9bf9-975d0d10b170
# ╠═c0b772d5-8345-4dc4-8801-45d4e614206d
# ╠═0e253811-3d8d-4c3e-a9d7-0e7e071ce7a1
# ╠═9f33023b-cb6a-4271-a2c8-5406e8cd18fb
# ╠═6a74d1df-6e26-4fbe-bec7-242292910180
# ╠═dc216f7c-ee56-4510-aa05-891af745da03
# ╠═4fd816a7-1c9e-46d5-ba80-62914ff54998
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: qiyang-ustc/FTensorKit

max_issues_repo_name: qiyang-ustc/FTensorKit

max_forks_repo_name: qiyang-ustc/FTensorKit


text:
module Yodel

using LightXML # LightXML is required to parse the XML file

export YodelEngine,Route,getRoute,getRouteWithController,isRoute,isRouteWithController

type YodelEngine
  routes::Array # Routes of Yodel Object

  # The YodelEngine constructor parses the XML file
  function YodelEngine(xmlPath::ASCIIString)
    doc = parse_file(xmlPath)
    xroot = root(doc)
    routes = []
    xRoutes = get_elements_by_tagname(xroot,"route")
    for xRoute in xRoutes
      if typeof(attribute(xRoute,"variables")) == Void
        variables = []
      else
        variables = separateAttributes(attribute(xRoute,"variables"))
      end

      tempRoute = Route(attribute(xRoute,"url"),variables,attribute(xRoute,"controller"))

      push!(routes,tempRoute)
    end

    new(routes)
  end
end

# Object Route which represent each route
type Route
  url::ASCIIString
  variables::Array
  controller::ASCIIString
  function Route(url::ASCIIString,variables::Array,controller::ASCIIString)
    new(url,variables,controller)
  end
end

function separateAttributes(attributesString) # Parses the variable attribute of the Route XML element to an array
  attributesString = string(attributesString)
  attributes = []

  attributesMatches = eachmatch(r",",attributesString)
  lastIndex = 0

  for att in attributesMatches
    push!(attributes,attributesString[lastIndex+1:att.offset])
    lastIndex = att.offset
  end

  push!(attributes,attributesString[lastIndex+1:end])
  return attributes
end

# Returns whether a function has a route for a specific URL or not
function isRoute(ydl::YodelEngine,url::ASCIIString)
  for route in ydl.routes
    if ismatch(Regex(route.url),url)
      return true
    end
  end
  return false
end

# Get a route that matches URL
function getRoute(ydl::YodelEngine,url::ASCIIString)
  for route in ydl.routes
    if ismatch(Regex(route.url),url)
      return route
    end
  end
  return Route("/",[],"Default")
end

# Get a route with a given controller
function isRouteWithController(ydl::YodelEngine,controller::ASCIIString)
  for route in ydl.routes
    if controller == route.controller
      return true
    end
  end
  return false
end
# Get a route with a given controller
function getRouteWithController(ydl::YodelEngine,controller::ASCIIString)
  for route in ydl.routes
    if controller == route.controller
      return route
    end
  end
  return Route("/",[],"Default")
end

# Extract variable from the url of a given route
function getVariable(url::ASCIIString,route::Route,variable::ASCIIString)
  #TODO
end

end # end of Yodel Module

max_stars_repo_name: PhoenixMachina/Yodel

max_issues_repo_name: PhoenixMachina/Yodel

max_forks_repo_name: PhoenixMachina/Yodel


text:
using myexample
using Test

myexample.f(2)

@testset "myexample.jl" begin
    # Write your own tests here.
    @test myexample.f(2) == 4
    @test myexample.f(0) == 0
    @test myexample.f(1) == 2
end


text:
using Luna

radius = 125e-6 # HCF core radius
flength = 1.5 # HCF length

gas = :He
pressure = 0.8 # gas pressure at the start of the HCF in bar

λ0 = 800e-9 # central wavelength of the pump pulse
τfwhm = 7.5e-15 # FWHM duration of the pump pulse
energy = 275e-6 # energy in the pump pulse

# giving the pressure argument as (pressure, 0) defines a decreasing gradient to vacuum
vuv = prop_capillary(radius, flength, gas, (pressure, 0); λ0, τfwhm, energy,
                     modes=4, trange=400e-15, λlims=(90e-9, 4e-6))

Plotting.prop_2D(vuv, :λ; modes=:sum, trange=(-20e-15, 20e-15), λrange=(100e-9, 1000e-9),
                 dBmin=-30)
# plot the total time-dependent power at the end of the propagation:
Plotting.time_1D(vuv; modes=:sum)
# plot the VUV dispersive wave at the exit:
Plotting.time_1D(vuv; modes=:sum, bandpass=(100e-9, 140e-9))
# show the evolution of the VUV pulse during the propagation:
Plotting.prop_2D(vuv; λrange=(100e-9, 200e-9), trange=(0, 20e-15), modes=:sum, bandpass=(100e-9, 140e-9))


text:
mutable struct LineSampling
    lines::Integer
    points::Vector{<:Real}
    direction::NamedTuple

    function LineSampling(
        lines::Integer,
        points::Vector{<:Real}=collect(0.5:0.5:5),
        direction::NamedTuple=NamedTuple(),
    )
        return new(lines, points, direction)
    end
end

function sample(inputs::Array{<:UQInput}, sim::LineSampling)
    random_inputs = filter(i -> isa(i, RandomUQInput), inputs)
    deterministic_inputs = filter(i -> isa(i, DeterministicUQInput), inputs)

    n_rv = count_rvs(random_inputs)
    n_samples = length(sim.points) * sim.lines
    rv_names = names(random_inputs)

    α = map(n -> sim.direction[n], rv_names)
    α /= norm(α)

    θ = rand(Normal(), n_rv, sim.lines)

    θ = θ - α * (α' * θ)
    θ = repeat(θ; outer=[length(sim.points), 1])

    θ = θ[:] + repeat(α * sim.points'; outer=[1, sim.lines])[:]

    samples = transpose(reshape(θ, n_rv, n_samples))
    samples = DataFrame(rv_names .=> eachcol(samples))

    if !isempty(deterministic_inputs)
        samples = hcat(samples, sample(deterministic_inputs, n_samples))
    end

    to_physical_space!(inputs, samples)

    return samples
end

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
#=
This module houses the built-in feature values for a variety of non-contextual atom features and also some convenience functions for constructing them easily.
=#
module ElementFeatureUtils

using DataFrames
using Flux: onecold
using ...ChemistryFeaturization.Codec: build_onehot_vec
using ...ChemistryFeaturization.Data: atom_data_df, feature_info

# export things
export default_nbins, oom_threshold_log
export atom_data_df, avail_feature_names
export categorical_feature_names, categorical_feature_vals, continuous_feature_names
export default_log, fea_minmax, default_categorical
export get_bins, get_param_vec
export onehot_lookup_encoder, onecold_decoder, encodable_elements

# default number of bins for continuous features, if unspecified
const default_nbins = 10
# if values of a feature span more than this many orders of magnitude, log-space it by default (open to better names for this...)
const oom_threshold_log = 2

# read in features...
const categorical_feature_names = feature_info["categorical"]
const categorical_feature_vals = Dict(
    fea => sort(collect(Set(skipmissing(atom_data_df[:, fea])))) for
    fea in categorical_feature_names
)
# but I want blocks to be in my order
categorical_feature_vals["Block"] = ["s", "p", "d", "f"]
const continuous_feature_names = feature_info["continuous"]
const avail_feature_names =
    cat(categorical_feature_names, continuous_feature_names; dims = 1)

"Compute the minimum and maximum possible values of a feature."
function fea_minmax(feature_name::String, lookup_table::DataFrame = atom_data_df)
    @assert feature_name in names(lookup_table) "Feature $feature_name isn't in the lookup table!"
    return [
        f(skipmissing(lookup_table[:, Symbol(feature_name)])) for f in [minimum, maximum]
    ]
end

"""
    default_log(feature_name, lookup_table = atom_data_df; threshold = oom_threshold_log)

Determine whether a continuous-valued feature should have logarithmically spaced bins.

Operates by finding the minimum and maximum values the feature can take on and comparing their ratio to a specified order-of-magnitude threshold that defaults to a package constant if not provided.
"""
function default_log(
    feature_name::String,
    lookup_table::DataFrame = atom_data_df;
    threshold::Real = oom_threshold_log,
)
    min_val, max_val = fea_minmax(feature_name, lookup_table)
    local log
    if typeof(min_val) <: Number
        signs = sign.([min_val, max_val])
        same_sign = all(x -> x == signs[1], signs)
        if same_sign
            oom_arg = sign(min_val) < 0 ? min_val / max_val : max_val / min_val
            oom = log10(oom_arg)
            log = oom > threshold
        else
            log = false
        end
    else
        log = false
    end
    return log
end

"""
    default_categorical(feature_name, lookup_table = atom_data_df)

Determine if a feature should be treated as categorical or continuous-valued.

If the value type is not a number, always returns true. If it is, checks whether it is in the built-in list of categorical features.

TODO: possibly add behavior where it will default to categorical if there is below some threshold number of discrete values?
"""
function default_categorical(feature_name::String, lookup_table::DataFrame = atom_data_df)
    local categorical
    if feature_name in avail_feature_names
        if feature_name in categorical_feature_names
            categorical = true
        else
            categorical = false
        end
    else
        feature_type = eltype(skipmissing(lookup_table[:, Symbol(feature_name)]))
        if feature_type <: Number
            categorical = false
        else
            categorical = true
        end
    end
    return categorical
end

"Little helper function to check that the logspace/categorical vector/boolean is appropriate and convert it to a vector as needed."
function get_param_vec(vec, num_features::Integer; pad_val = false)
    if !(typeof(vec) <: Vector)
        output_vec = [vec for i = 1:num_features]
    elseif length(vec) == num_features # specified properly
        output_vec = vec
    elseif length(vec) < num_features
        @info "Parameter vector too short. Padding end with $pad_val."
        output_vec = vcat(vec, [pad_val for i = 1:num_features-size(vec, 1)])
    elseif size(vec, 1) > num_features
        @info "Parameter vector too long. Cutting off at appropriate length."
        output_vec = vec[1:num_features]
    end
    return output_vec
end

"Helper function for encoder and decoder...(nbins is ignored for categorical=true)"
function get_bins(
    feature_name::String,
    lookup_table::DataFrame = atom_data_df;
    nbins::Integer = default_nbins,
    logspaced::Bool = default_log(feature_name, lookup_table),
    categorical::Bool = default_categorical(feature_name, lookup_table),
)
    local bins, min_val, max_val

    if categorical
        if feature_name in categorical_feature_names
            bins = categorical_feature_vals[feature_name]
        else
            bins = sort(unique(skipmissing(lookup_table[:, Symbol(feature_name)])))
        end
    else
        min_val, max_val = fea_minmax(feature_name, lookup_table)

        if isapprox(min_val, max_val)
            @warn "It looks like the minimum and maximum possible values of $feature_name are approximately equal. This could cause numerical issues with binning, and also this feature is likely uninformative. Perhaps reconsider if it needs to be included?"
        end

        if logspaced
            @assert all(x -> sign(x) == sign(min_val), [min_val, max_val]) "I don't know how to do a logarithmically spaced feature whose value can be zero! :("
            if sign(min_val) > 0
                bins = 10 .^ range(log10(min_val), log10(max_val), length = nbins + 1)
            else
                bins =
                    -1 .* (
                        10 .^
                        range(log10(abs(min_val)), log10(abs(max_val)), length = nbins + 1)
                    )
            end
        else
            bins = range(min_val, max_val, length = nbins + 1)
        end
    end
    return bins
end

"Function for encoding onehot-style vectors based on values in a lookup table. Intended to be used as an `encode_f` for `ElementFeatureDescriptor` objects. See source code of convenience constructor for `ElementFeatureDescriptor` for more details."
function onehot_lookup_encoder(
    el::String,
    feature_name::String,
    lookup_table::DataFrame = atom_data_df;
    nbins::Integer = default_nbins,
    logspaced::Bool = default_log(feature_name, lookup_table),
    categorical::Bool = default_categorical(feature_name, lookup_table),
)
    colnames = names(lookup_table)
    @assert (feature_name in colnames) && ("Symbol" in colnames) "Your lookup table must have a column called :Symbol and one with the same name as your feature to be usable!"

    feature_vals = lookup_table[:, [:Symbol, Symbol(feature_name)]]

    @assert el in feature_vals.Symbol "Element $el is not in the database! :("

    bins = get_bins(
        feature_name,
        lookup_table;
        nbins = nbins,
        logspaced = logspaced,
        categorical = categorical,
    )

    # pull value of feature for this element
    val = getproperty(feature_vals[feature_vals.Symbol.==el, :][1, :], Symbol(feature_name))
    build_onehot_vec(val, bins, categorical)
end

"Function for decoding onehot-style vectors based on values in a lookup table. Intended to be used as an `decode_f` for `ElementFeatureDescriptor` objects. See source code of convenience constructor for `ElementFeatureDescriptor` for more details."
function onecold_decoder(
    encoded,
    feature_name::String,
    lookup_table::DataFrame = atom_data_df;
    nbins::Integer = default_nbins,
    logspaced::Bool = default_log(feature_name, lookup_table),
    categorical::Bool = default_categorical(feature_name, lookup_table),
)
    colnames = names(lookup_table)
    @assert feature_name in colnames && "Symbol" in colnames "Your lookup table must have a column called :Symbol and one with the same name as your feature to be usable!"

    bins = get_bins(
        feature_name,
        lookup_table;
        nbins = nbins,
        logspaced = logspaced,
        categorical = categorical,
    )
    local decoded
    if categorical # return value
        decoded = onecold(encoded, bins)
    else # return bounds (TODO: should this be a tuple or a vector..? I like tuple for distinguishing from encoded vectors, but it doesn't play so nice with broadcasting...)
        decoded = (onecold(encoded, bins[1:end-1]), onecold(encoded, bins[2:end]))
    end
    return decoded
end

end


text:
using BenchmarkTools
#
#	tool to replace params with values in an expression
#

expr_replace!(ex, s, v) = ex
expr_replace!(ex::Symbol, s, v) = s == ex ? v : ex

function expr_replace!(ex::Expr, s, v)
    for i=1:length(ex.args)
        ex.args[i] = expr_replace!(ex.args[i], s, v)
    end
    return ex
end

expr_replace(ex, s, v) = expr_replace!(copy(ex), s, v)

function expr_replace(ex, nt)
    e = copy(ex)
    k=keys(nt)
    v=values(nt)
    for i in 1:length(nt)
        expr_replace!(e, k[i], v[i])
    end
    return e
end

# see this section of docs https://docs.julialang.org/en/v1/base/arrays/#Broadcast-and-vectorization-1

t = collect(0.0:0.1:10000)
params = [1.0, 1.0, 1.0, 1.0, 1.0]
nt=(G₀=1., G₁=1., η₁=1., G₂=1., η₂=1.)

#
#	Vectorised function
#

function G_vec(t::Vector{T}, params::Vector{T}) where T<:Real
   G₀, G₁, η₁, G₂, η₂ = params
   G = G₀ .+ G₁*exp.(-t*G₁/η₁) .+ G₂*exp.(-t*G₂/η₂)
end

println("test vectorised version of function")
@btime G_vec($t, $params)
println("\n")

#
#	Broadcasting functions
#

function G_single(t::T, params::Array{T,1}) where T<:Real
   G₀, G₁, η₁, G₂, η₂ = params
   G = G₀ + G₁*exp(-t*G₁/η₁) + G₂*exp(-t*G₂/η₂)
end

println("test broadcasting using . notation")
@btime G_single.($t, $(params,))
println("\n")

println("test broadcast function explicitly (should be the same as . notation)")
@btime broadcast($G_single, $t, $(params,))
println("\n")

#
#	Array comprehension
#

function G_arraycomp(t::Vector{T}, params::Vector{T}) where T<:Real
    [G_single(ti, params) for ti in t]
end

println("test using array comprehension")
@btime G_arraycomp($t, $params)
println("\n")

#
#	parameter substitution
#

e = :( G₀ + G₁*exp(-t*G₁/η₁) + G₂*exp(-t*G₂/η₂) )

function G_arrayexp(ta::Vector{T}, nt) where T<:Real
	e2 = expr_replace(e, nt)
    @eval f(t::T) where T<:Real = begin  $e2   end
    #@eval [f(t) for t in $ta]
    @eval broadcast(f,$ta)
end

println("test using local function from expression with parameter replacement")
@btime G_arrayexp($t, $nt)
println("\n")

println("test using array comprehension in external function from expression with parameter replacement")

function paramsubst(e,nt)
	e2 = expr_replace(e, nt)
    @eval f(t::T) where T<:Real = begin  $e2   end
    @eval g(ta::Vector{T}) where T<:Real = begin [f(t) for t in ta] end
    #g = ta -> [f(t) for t in ta]
    return g
end

@btime paramsubst(e,nt)
g=paramsubst(e,nt)

@btime g($t)
println("\n")

println("test using array broadcast in external function from expression with parameter replacement")
    @eval g2(ta::Vector{T}) where T<:Real = begin broadcast(f, ta) end
    #g2 = t -> broadcast(f, t)

@btime g2($t)
println("\n")

println("test using array comprehension with external lambda function")
    gl = ta -> G_arraycomp(  ta, params)

@btime gl($t)
println("\n")

println("test using array braodcast with external lambda function")
    gl2 = ta -> broadcast(G_single, ta, (params,))

@btime gl2($t)
println("\n")

#
#	for loop
#

function G_loop(t::Vector{T}, params::Vector{T}) where T<:Real

    out = Array{T}(undef, length(t))

    for (i,ti) in enumerate(t)
        out[i] = G_single(ti, params)
    end

    out
end

println("test using for loop")
@btime G_loop($t, $params)
println("\n")

println("Benchmarks complete")

max_stars_repo_path: 20190321_alex follow up 2.jl

max_stars_repo_name: JuliaRheology/RheoBenchAndTest

max_issues_repo_path: 20190321_alex follow up 2.jl

max_issues_repo_name: JuliaRheology/RheoBenchAndTest

max_forks_repo_path: 20190321_alex follow up 2.jl

max_forks_repo_name: JuliaRheology/RheoBenchAndTest


text:
using PortAudio
import Base.Threads.@spawn

export keyboard

function getchar()
    ret = ccall(:jl_tty_set_mode, Int32, (Ptr{Cvoid},Int32), stdin.handle, true)
    ret == 0 || error("unable to switch to raw mode")
    c = read(stdin, Char)
    ccall(:jl_tty_set_mode, Int32, (Ptr{Cvoid},Int32), stdin.handle, false)
    return c
end

function keyboard(instruments::AbstractArray{Instrument})
    @assert length(instruments) > 0

    instr = 1
    instrument() = instruments[instr]

    isrunning = true
    audio_task = @spawn begin
        buf = zeros(Float64, 128)
        stream = PortAudioStream(0, 1; latency = 0.009, blocksize=16)
        try
            while isrunning
                write(stream, next!(buf, instrument()))
            end
        finally
            close(stream)
            println("audio stream is closed")
        end
    end

    while true
        c = getchar()

        tuning = instrument().tuning
        dict = keymap(tuning)
        if haskey(dict, c)
            tone = dict[c]
            resound!(instrument(), 440.0 * interval(tuning, tone))
        end
        val = UInt(c)
        if val == 3 || val == 27
            break
        end

        if c == '+' || c == '_'
            N = length(instruments)
            instr += ifelse(
                        c=='+',
                        ifelse(instr < N, 1, -N+1),
                        ifelse(instr > 1, -1, N-1))

            println("Now using instrument $instr: $(instrument()) with $tuning")

            @show instr
        end
    end

    isrunning = false
    wait(audio_task)
end


text:
using EliminateGraphs
using Test

@testset "mis1" begin
    @test EliminateGraphs.minx(x->x^2, [2,-1, 3]) == -1

    graph = EliminateGraph([0 1 0 0 0;
                            1 0 1 1 0;
                            0 1 0 1 0;
                            0 1 1 0 1;
                            0 0 0 1 0])

    @test mis1(graph) == 3
end


text:
# function gen_opcsp(r, positions, covariance, distance_limit=1.0, start=1, stop=-1; rng=nothing)
#     d = rand!(rng, Array(Float64, length(r)), MVN(zeros(length(r)), covariance))
#     return OPCSP(r, d, positions, covariance, distance_limit, start, stop)
# end

function gen_op(;distance_limit=1.5,
                 n_nodes=10,
                 rng=MersenneTwister())

    positions = Array(Vector{Float64}, n_nodes)
    for i in 1:n_nodes
        positions[i] = rand(rng, 2).-0.5
    end

    return SimpleOP(5.0*rand(rng, n_nodes) + 5.0,
                    positions,
                    distance_limit,
                    1,
                    n_nodes)
end

function gen_two_cluster_problem(;distance_limit=6.0,
                                  n_nodes=10,
                                  p=0.4, # probability of connection
                                  rng=MersenneTwister(),
                                  noise=1.0
                                  )
    positions = Array(Vector{Float64}, n_nodes)
    n_1 = trunc(Int,n_nodes/2)
    for i in 1:n_1
        positions[i] = rand(rng, 2)
        # positions[i] = sample(MVN([-0.5,0], eye(2)), rng=rng)
    end
    for i in n_1+1:n_nodes
        positions[i] = -rand(rng,2)
        # positions[i] = sample(MVN([0.5,0], eye(2)), rng=rng)
    end

    # generate correlations with erdos-renyi connectivity
    covariance = zeros(n_nodes, n_nodes)
    for i in 1:n_nodes
        for j in i+1:n_nodes
            if rand(rng) < p
                sgn = rand(rng)
                if sgn >= 0.5
                    rho = 0.99
                else
                    rho = -0.99
                end
                # rho = 2.0*rand(rng)-1.0
                # stdev = randn(rng,2)
                stdev = rand(rng,2)
                block = [     stdev[1]^2 rho*prod(stdev);
                         rho*prod(stdev)      stdev[2]^2]
                covariance[[i,j], [i,j]] += block
            end
        end
    end
    for i in 1:n_nodes
        if covariance[i,i] == 0.0
            covariance[i,i] = randn(rng)^2
        end
    end
    covariance = covariance.*(noise^2/(det(covariance)^(1/n_nodes)))

    @assert isposdef(covariance)
    @assert norm(positions[1]-positions[n_nodes]) <= distance_limit

    r = 5.0*rand(rng, n_nodes) + 5
    return OPCSP(r, positions, covariance, distance_limit, 1, n_nodes)
end

function gen_problem(;distance_limit=sqrt(2),
                      n_nodes=10,
                      p=0.4, # probability of connection
                      rng=MersenneTwister(),
                      noise=1.0
                      )

    positions = Array(Vector{Float64}, n_nodes)
    for i in 1:n_nodes
        positions[i] = rand(rng, 2)
    end

    # generate correlations with erdos-renyi connectivity
    covariance = zeros(n_nodes, n_nodes)
    for i in 1:n_nodes
        for j in i+1:n_nodes
            if rand(rng) < p
                sgn = rand(rng)
                if sgn >= 0.5
                    rho = 0.99
                else
                    rho = -0.99
                end
                # rho = 2.0*rand(rng)-1.0
                # stdev = randn(rng,2)
                stdev = rand(rng,2)
                block = [     stdev[1]^2 rho*prod(stdev);
                         rho*prod(stdev)      stdev[2]^2]
                covariance[[i,j], [i,j]] += block
            end
        end
    end
    for i in 1:n_nodes
        if covariance[i,i] == 0.0
            covariance[i,i] = randn(rng)^2
        end
    end
    covariance = covariance.*(noise^2/(det(covariance)^(1/n_nodes)))

    @assert isposdef(covariance)
    @assert norm(positions[1]-positions[n_nodes]) <= distance_limit

    # r = 5.0*rand(rng, n_nodes) + 5
    # r = 10.0*rand(rng, n_nodes)
    r = 10*ones(n_nodes)
    return OPCSP(r, positions, covariance, distance_limit, 1, n_nodes)
end

function gen_highly_connected(;distance_limit=sqrt(2),
                      n_nodes=10,
                      n_connected=2,
                      p=0.5, # probability of connection from a connected node to another node
                      rbar=10.0,
                      noise_2norm=1.0,
                      rng=MersenneTwister()
                      )

    positions = Array(Vector{Float64}, n_nodes)
    for i in 1:n_nodes
        positions[i] = rand(rng, 2)
    end

    covariance = zeros(n_nodes, n_nodes)
    open = IntSet(2:n_nodes-1)
    for i in 1:n_connected
        j = rand(rng, collect(open))
        delete!(open, j)
        others = IntSet(1:n_nodes)
        delete!(others, j)
        for k in others
            if rand(rng) < p
                sgn = rand(rng)
                if sgn >= 0.5
                    rho = 0.99
                else
                    rho = -0.99
                end
                # rho = 2.0*rand(rng)-1.0
                # stdev = randn(rng,2)
                stdev = rand(rng,2)
                block = [     stdev[1]^2 rho*prod(stdev);
                         rho*prod(stdev)      stdev[2]^2]
                covariance[[j,k], [j,k]] += block
            end
        end
    end
    for i in 1:n_nodes
        if covariance[i,i] == 0.0
            covariance[i,i] = 0.01^2
        end
    end

    # normalize so that noise sqrt(trace(covariance)) = noise
    alpha = noise_2norm^2/trace(covariance)
    covariance .*= alpha

    @assert abs(sqrt(trace(covariance))-noise_2norm) < 1e-5
    @assert isposdef(covariance + 1e-5*eye(covariance))
    @assert norm(positions[1]-positions[n_nodes]) <= distance_limit

    # r = 5.0*rand(rng, n_nodes) + 5
    # r = 10.0*rand(rng, n_nodes)
    r = rbar*ones(n_nodes)
    return OPCSP(r, positions, covariance, distance_limit, 1, n_nodes)
end

function gen_informative(;distance_limit=sqrt(2),
                      n_nodes=10,
                      n_connected=n_nodes-2,
                      rbar=10.0,
                      noise_2norm=1.0,
                      rng=MersenneTwister()
                      )

    positions = Array(Vector{Float64}, n_nodes)
    for i in 1:n_nodes
        positions[i] = rand(rng, 2)
    end

    covariance = zeros(n_nodes, n_nodes)
    open = IntSet(2:n_nodes-1)
    for i in 1:n_connected
        j = rand(rng, collect(open))
        delete!(open, j)
        others = IntSet(1:n_nodes)
        delete!(others, j)
        first = rand(rng, collect(others))
        delete!(others, first)
        second = rand(rng, collect(others))
        stdev = rand(rng,3)
        block = [         2.0*stdev[1]^2 .99*prod(stdev[1:2]) -.99*prod(stdev[[1,3]]);
                    .99*prod(stdev[1:2])           stdev[2]^2                       0;
                 -.99*prod(stdev[[1,3]])                    0               stdev[3]^2 ]
        covariance[[j,first,second], [j,first,second]] += block
    end
    for i in 1:n_nodes
        if covariance[i,i] == 0.0
            covariance[i,i] = 0.01^2
        end
    end

    # normalize so that noise sqrt(trace(covariance)) = noise
    alpha = noise_2norm^2/trace(covariance)
    covariance .*= alpha

    @assert abs(sqrt(trace(covariance))-noise_2norm) < 1e-5
    @assert isposdef(covariance + 1e-5*eye(covariance))
    @assert norm(positions[1]-positions[n_nodes]) <= distance_limit

    # r = 5.0*rand(rng, n_nodes) + 5
    # r = 10.0*rand(rng, n_nodes)
    r = rbar*ones(n_nodes)
    return OPCSP(r, positions, covariance, distance_limit, 1, n_nodes)
end


text:
using FountainCodes
using SparseArrays
using DelimitedFiles
using Random, StatsBase, Distributions

function load_matrix(filename="./examples/H0.txt")
    M = readdlm(filename)
    Is = M[:, 2] .+ 1
    Js = M[:, 1] .+ 1
    Vs = ones(GF256, length(Is))
    # Vs = FountainCodes.rand_nonzero(Random.GLOBAL_RNG, GF256, length(Is))
    # Vs = randn(length(Is))
    println("Loaded matrix $filename with $(eltype(Vs)) entries")

    # set the first entry of each constraint to 1
    A = sparse(Is, Js, Vs, 1536, 512)
    vals = nonzeros(A)
    for j in 1:size(A, 2)
        i = 1
        k = nzrange(A, j)[i]
        vals[k] = one(eltype(vals))
    end

    A
end

"""

Write a `AbstractSparseMatrixCSC` to file, with each line corresponding to a non-zero entry.
"""
function write_matrix(A::SparseArrays.AbstractSparseMatrixCSC, filename="R10.txt"; sep="\t")
    Is, Js, Vs = findnz(A)
    open(filename, "w") do io
        for i in 1:length(Is)
            println(io, "$(Is[i]-1)$(sep)$(Js[i]-1)\n")
        end
    end
    return
end

function simulate_ldpc(Ht, enc, erasurerate, erasures)
    @assert erasures == "fixed" || erasures == "random"
    enc .= zero(eltype(enc))
    nenc, nconstraints = size(Ht)
    Xs = collect(1:nenc)
    if erasures == "fixed"
        shuffle!(Xs)
        nerasures = round(Int, erasurerate*nenc)
        Xs_rec = sort(Xs[1:nenc-nerasures])
        Xs_era = sort(Xs[end-nerasures+1:end])
    elseif erasures == "random"
        mask = rand(nenc) .< erasurerate
        Xs_rec = Xs[.!mask]
        Xs_era = Xs[mask]
    end
    Ht_rec = Ht[Xs_rec, :] # constraints corresponding to received symbols
    Ht_era = Ht[Xs_era, :] # constraints corresponding to erased symbols
    enc_rec = enc[Xs_rec]
    enc_era = enc[Xs_era]

    # attempt to create the decoder
    nrowops, ninactivated = 0, 0
    decoding_success = false
    decoder = Decoder(Ht_era)

    # attempt decoding
    try
        rhs = Ht_rec'*enc_rec
        enc_era_dec = decode(Ht_era, rhs; decoder)
        # if eltype(enc_era_dec) == Float64
        #     # @assert
        #     if !(enc_era_dec ≈ enc_era)
        #         println(enc_era_dec)
        #     end
        # else
        #     @assert enc_era_dec == enc_era
        # end
        decoding_success = true
    catch e
        if e isa RankDeficiencyException
            # if eltype(enc_era) == Float64
            #     @assert rank(Matrix(Ht_era)) < length(enc_era)
            # end
        else
            rethrow()
        end
    end
    nrowops = decoder.num_rowops
    ninactivated = decoder.num_inactivated
    nrowops, ninactivated, decoding_success
end

function main_ldpc(;erasurerate, minsamples=10000, maxsamples=100000, minfails=100, erasures)
    Ht = load_matrix()
    H = sparse(Ht')
    nconstraints, nenc = size(H)
    enc = zeros(eltype(H), nenc)

    # run simulations
    samples = Vector{Any}()
    nfails = 0
    while length(samples) < minsamples || (length(samples) < maxsamples && (nfails < minfails || length(samples) - nfails < minfails))
        sample = simulate_ldpc(Ht, enc, erasurerate, erasures)
        nfails += 1 - sample[3]
        push!(samples, sample)
    end

    # compute and return means
    nops = [sample[1] for sample in samples]
    nina = [sample[2] for sample in samples]
    fails = 1 - mean([sample[3] for sample in samples])
    length(samples), mean(nops), mean(nina), fails
end

function simulate_raptor(code, int, nenc, erasurerate, erasures)
    @assert erasures == "fixed" || erasures == "random"

    # randomly selected nenc unique encodes symbol identifiers, corresponding to the nenc encoded
    # symbols received by the decoder (the upper limit is chosen arbitrarily)
    Xs = collect(0:nenc-1)
    if erasures == "fixed"
        shuffle!(Xs)
        nerasures = round(Int, erasurerate * nenc)
        Xs = sort!(Xs[1:nenc-nerasures])
    elseif erasures == "random"
        mask = rand(nenc) .< erasurerate
        Xs = Xs[.!mask]
    end

    # compute the encoded symbols
    # G = generator_matrix(code, Xs)
    # enc = G'*int

    # decode the intermediate symbols from the encoded symbols
    A = constraint_matrix(code, Xs)
    pre = zeros(GF256, size(A, 2))
    # pre[end-length(enc)+1:end] .= enc

    # attempt to create the decoder
    nrowops, ninactivated = 0, 0
    decoding_success = false
    decoder = Decoder(A)

    # attempt decoding
    try

        # permanent inactivations
        if code isa RQ
            for cpi in (code.L-code.P+1):code.L
                FountainCodes.mark_inactive!(decoder, A, cpi)
            end
        end

        # decode
        dec_int = decode(A, pre; decoder)
        decoding_success = true
    catch e
        if e isa RankDeficiencyException
        else
            rethrow()
        end
    end

    # only count non-permanent inactivations
    ninactivated = decoder.num_inactivated
    if code isa RQ
        ninactivated -= code.P
    end
    nrowops = decoder.num_rowops
    nrowops, ninactivated, decoding_success
end

function main_raptor(;nsrc=1024, rate=2/3, erasurerate, minsamples=10000, maxsamples=100000, minfails=100, style, erasures)
    @assert style == "R10" || style == "RQ"

    code = style == "R10" ? R10(nsrc) : RQ(nsrc)
    nsrc = style == "R10" ? nsrc : code.Kp

    nenc = round(Int, nsrc / rate)
    src = rand(GF256, nsrc)

    # compute the intermediate symbols
    A = constraint_matrix(code, 0:nsrc-1)
    pre = zeros(GF256, size(A, 2))
    pre[end-nsrc+1:end] .= src
    int = decode(A, pre)

    # run simulations
    samples = Vector{Any}()
    nfails = 0
    while length(samples) < minsamples || (length(samples) < maxsamples && (nfails < minfails || length(samples) - nfails < minfails))
        sample = simulate_raptor(code, int, nenc, erasurerate, erasures)
        nfails += 1 - sample[3]
        push!(samples, sample)
    end

    # compute and return means
    nops = [sample[1] for sample in samples]
    nina = [sample[2] for sample in samples]
    fails = 1 - mean([sample[3] for sample in samples])
    length(samples), mean(nops), mean(nina), fails
end

function simulate_lt(code, nenc, erasurerate, erasures)
    @assert erasures == "fixed" || erasures == "random"

    # randomly selected nenc unique encodes symbol identifiers, corresponding to the nenc encoded
    # symbols received by the decoder (the upper limit is chosen arbitrarily)
    Xs = collect(0:nenc-1)
    if erasures == "fixed"
        shuffle!(Xs)
        nerasures = round(Int, erasurerate * nenc)
        Xs = sort!(Xs[1:nenc-nerasures])
    elseif erasures == "random"
        mask = rand(nenc) .< erasurerate
        Xs = Xs[.!mask]
    end

    # compute the encoded symbols
    G = generator_matrix(code, Xs)
    enc = zeros(eltype(G), length(Xs))

    # attempt to create the decoder
    nrowops, ninactivated = 0, 0
    decoding_success = false
    decoder = Decoder(G)

    # attempt decoding
    try
        # decode
        decode(G, enc; decoder)
        decoding_success = true
    catch e
        if e isa RankDeficiencyException
        else
            rethrow()
        end
    end
    nrowops = decoder.num_rowops
    ninactivated = decoder.num_inactivated
    nrowops, ninactivated, decoding_success
end

function main_lt(;nsrc=1024, rate=2/3, erasurerate, minsamples=10000, maxsamples=100000, minfails=100, erasures)
    dd = Soliton(nsrc, 2, 1e-6, 1e-3)
    code = LT{GF256}(nsrc, dd) # w. generator matrix coefficients of type GF256
    nenc = round(Int, nsrc / rate)

    # run simulations
    samples = Vector{Any}()
    nfails = 0
    while length(samples) < minsamples || (length(samples) < maxsamples && (nfails < minfails || length(samples) - nfails < minfails))
        sample = simulate_lt(code, nenc, erasurerate, erasures)
        nfails += 1 - sample[3]
        push!(samples, sample)
    end

    # compute and return means
    nops = [sample[1] for sample in samples]
    nina = [sample[2] for sample in samples]
    fails = 1 - mean([sample[3] for sample in samples])
    length(samples), mean(nops), mean(nina), fails
end

function main(erasures="fixed")
    erasurerates = 0.28:0.01:0.33

    # LT simulations
    nsamples_lt = zeros(0)
    rowops_lt = zeros(0)
    inactivations_lt = zeros(0)
    fails_lt = zeros(0)
    for erasurerate in erasurerates
        println("Simulating lt at $erasurerate")
        nsamples, rowops, inactivations, fails = main_lt(;erasurerate, erasures)
        push!(nsamples_lt, nsamples)
        push!(rowops_lt, rowops)
        push!(inactivations_lt, inactivations)
        push!(fails_lt, fails)
    end
    println("LT")
    println(nsamples_lt)
    println(rowops_lt)
    println(inactivations_lt)
    println(fails_lt)

    # R10 simulations
    println("Simulating R10")
    nsamples_r10 = zeros(0)
    rowops_r10 = zeros(0)
    inactivations_r10 = zeros(0)
    fails_r10 = zeros(0)
    for erasurerate in erasurerates
        println("Simulating R10 at $erasurerate")
        nsamples, rowops, inactivations, fails = main_raptor(;erasurerate, style="R10", erasures)
        push!(nsamples_r10, nsamples)
        push!(rowops_r10, rowops)
        push!(inactivations_r10, inactivations)
        push!(fails_r10, fails)
    end
    println("R10")
    println(nsamples_r10)
    println(rowops_r10)
    println(inactivations_r10)
    println(fails_r10)

    # RQ simulations
    println("Simulating RQ")
    nsamples_rq = zeros(0)
    rowops_rq = zeros(0)
    inactivations_rq = zeros(0)
    fails_rq = zeros(0)
    for erasurerate in erasurerates
        println("Simulating RQ at $erasurerate")
        nsamples, rowops, inactivations, fails = main_raptor(;erasurerate, style="RQ", erasures)
        push!(nsamples_rq, nsamples)
        push!(rowops_rq, rowops)
        push!(inactivations_rq, inactivations)
        push!(fails_rq, fails)
    end
    println("RQ")
    println(nsamples_rq)
    println(rowops_rq)
    println(inactivations_rq)
    println(fails_rq)

    # LDPC simulations
    nsamples_ldpc = zeros(0)
    rowops_ldpc = zeros(0)
    inactivations_ldpc = zeros(0)
    fails_ldpc = zeros(0)
    for erasurerate in erasurerates
        println("Simulating LDPC at $erasurerate")
        nsamples, rowops, inactivations, fails = main_ldpc(;erasurerate, erasures)
        push!(nsamples_ldpc, nsamples)
        push!(rowops_ldpc, rowops)
        push!(inactivations_ldpc, inactivations)
        push!(fails_ldpc, fails)
    end
    println("LDPC")
    println(nsamples_ldpc)
    println(rowops_ldpc)
    println(inactivations_ldpc)
    println(fails_ldpc)
    return
end

max_stars_repo_name: severinson/RaptorCodes

max_issues_repo_name: severinson/RaptorCodes

max_forks_repo_name: severinson/RaptorCodes


text:
struct QLayout <: MemoryLayout end

MemoryLayout(::Type{<:AbstractQ}) = QLayout()

transposelayout(::QLayout) = QLayout()

copy(M::Lmul{QLayout}) = copyto!(similar(M), M)

function copyto!(dest::AbstractArray{T}, M::Lmul{QLayout}) where T
    A,B = M.A,M.B
    if size(dest,1) == size(B,1)
        copyto!(dest, B)
    else
        copyto!(view(dest,1:size(B,1),:), B)
        zero!(@view(dest[size(B,1)+1:end,:]))
    end
    materialize!(Lmul(A,dest))
end

function copyto!(dest::AbstractArray, M::Ldiv{QLayout})
    A,B = M.A,M.B
    copyto!(dest, B)
    materialize!(Ldiv(A,dest))
end

materialize!(M::Ldiv{QLayout}) = materialize!(Lmul(M.A',M.B))

_qr(layout, axes, A; kwds...) = Base.invoke(qr, Tuple{AbstractMatrix{eltype(A)}}, A; kwds...)
_qr(layout, axes, A, pivot::P; kwds...) where P = Base.invoke(qr, Tuple{AbstractMatrix{eltype(A)},P}, A, pivot; kwds...)
_lu(layout, axes, A; kwds...) = Base.invoke(lu, Tuple{AbstractMatrix{eltype(A)}}, A; kwds...)
_lu(layout, axes, A, pivot::P; kwds...) where P = Base.invoke(lu, Tuple{AbstractMatrix{eltype(A)},P}, A, pivot; kwds...)
_qr!(layout, axes, A, args...; kwds...) = error("Overload _qr!(::$(typeof(layout)), axes, A)")
_lu!(layout, axes, A, args...; kwds...) = error("Overload _lu!(::$(typeof(layout)), axes, A)")
_factorize(layout, axes, A) = Base.invoke(factorize, Tuple{AbstractMatrix{eltype(A)}}, A)

macro _layoutfactorizations(Typ)
    esc(quote
        LinearAlgebra.qr(A::$Typ, args...; kwds...) = ArrayLayouts._qr(ArrayLayouts.MemoryLayout(typeof(A)), axes(A), A, args...; kwds...)
        LinearAlgebra.qr!(A::$Typ, args...; kwds...) = ArrayLayouts._qr!(ArrayLayouts.MemoryLayout(typeof(A)), axes(A), A, args...; kwds...)
        LinearAlgebra.lu(A::$Typ, pivot::Union{Val{false}, Val{true}}; kwds...) = ArrayLayouts._lu(ArrayLayouts.MemoryLayout(typeof(A)), axes(A), A, pivot; kwds...)
        LinearAlgebra.lu(A::$Typ{T}; kwds...) where T = ArrayLayouts._lu(ArrayLayouts.MemoryLayout(typeof(A)), axes(A), A; kwds...)
        LinearAlgebra.lu!(A::$Typ, args...; kwds...) = ArrayLayouts._lu!(ArrayLayouts.MemoryLayout(typeof(A)), axes(A), A, args...; kwds...)
        LinearAlgebra.factorize(A::$Typ) = ArrayLayouts._factorize(ArrayLayouts.MemoryLayout(typeof(A)), axes(A), A)
    end)
end

macro layoutfactorizations(Typ)
    esc(quote
        ArrayLayouts.@_layoutfactorizations $Typ
        ArrayLayouts.@_layoutfactorizations SubArray{<:Any,2,<:$Typ}
    end)
end


text:
module WeakGradTests
   using Gridap
   using GridapGeosciences
   import Gridap.Fields: ∇, divergence
   using Plots
   using Test

   include("ConvergenceAnalysisTools.jl")

   function ω(xyz)
     θϕ = xyz2θϕ(xyz)
     θ,ϕ = θϕ
     cos(θ)*cos(ϕ)^2
   end

   function ωθϕ(θϕ)
     θ,ϕ = θϕ
     cos(θ)*cos(ϕ)^2
   end

   # spherical surface gradient computed analytically
   function gradω(xyz)
    θϕr = xyz2θϕr(xyz)
    θ,ϕ,r = θϕr
    dθ = -sin(θ)*cos(ϕ)
    dϕ = -2*cos(θ)*cos(ϕ)*sin(ϕ)
    spherical_to_cartesian_matrix(θϕr)⋅VectorValue(dθ,dϕ,0)
  end

  #  function divergence(::typeof(ωθϕ))
  #   divergence_unit_sphere(ωθϕ)
  #  end

   function compute_error_weak_grad(model,order,degree)
     # Setup geometry
     Ω=Triangulation(model)
     dΩ=Measure(Ω,degree)
     dω=Measure(Ω,degree,ReferenceDomain())

     # Setup H(div) spaces
     reffe_rt = ReferenceFE(raviart_thomas, Float64, order)
     V = FESpace(model, reffe_rt ; conformity=:HDiv)
     U = TrialFESpace(V)

     # L2 spaces
     reffe_lgn = ReferenceFE(lagrangian, Float64, order)
     S = FESpace(model, reffe_lgn; conformity=:L2)
     R = TrialFESpace(S)

     # Project the analytical field ω onto the L2 space
     a1(u,v) = ∫(v*u)dΩ
     b1(v)   = ∫(v*ω)dΩ
     op      = AffineFEOperator(a1,b1,S,R)
     wh      = solve(op)

     # Compute weak grad of wh
     a2(u,v) = ∫(v⋅u)dΩ
     b2(v)   = ∫((-wh)*DIV(v))dω
     op      = AffineFEOperator(a2,b2,U,V)
     gradwh  = solve(op)

     e    = gradwh-gradω
     dive = divergence(gradwh)-laplacian_unit_sphere(ωθϕ)∘xyz2θϕ
     sqrt(sum(∫(e⋅e+dive*dive)dΩ))
   end
   #model=CubedSphereDiscreteModel(24)
   #e,gradwh=compute_error_weak_grad(model,1,4)
   #writevtk(Triangulation(model),"RT124x24analyticalmap",cellfields=["gradwh"=>gradwh])

   @time ahs0,ak0errors,as0=convergence_study(compute_error_weak_grad,generate_n_values(2),0,4)
   @test as0 ≈ 0.8346320885900106

   @time ahs1,ak1errors,as1=convergence_study(compute_error_weak_grad,generate_n_values(2,n_max=50),1,8)
   @test as1 ≈ 1.1034326200306834

  #  plotd=plot([ahs0,ahs1],[ak0errors,ak1errors],
  #  xaxis=:log, yaxis=:log,
  #  label=["k=0 analytical map" "k=1 analytical map"],
  #  shape=:auto,
  #  xlabel="h",ylabel="H(div) error norm", legend=:bottomright)

  #  savefig(plotd,"L2_error_weak_div_perp.png")

end


text:
using Documenter
using Literate
using RadiativeTransfer
using RadiativeTransfer.Absorption
using RadiativeTransfer.Scattering
using RadiativeTransfer.vSmartMOM
using Plots

ENV["GKSwstype"] = "nul"

function build()

    tutorials = ["Tutorial_Absorption.jl", "Tutorial_Scattering.jl"] # ,
    tutorials_paths = [joinpath(@__DIR__, "src", "pages", "tutorials", tutorial) for tutorial in tutorials]

    for tutorial in tutorials_paths
        Literate.markdown(tutorial, joinpath(@__DIR__, "src", "pages", "tutorials"))
    end

    tutorials_md = [joinpath("pages", "tutorials", tutorial[1:end-3]) * ".md" for tutorial in tutorials]

    pages = Any[
        "Home"                  => "index.md",
        "RadiativeTransfer"     => "pages/RadiativeTransfer.md",
        "Absorption"            => "pages/Absorption.md",
        "Scattering"            => "pages/Scattering.md",
        "Tutorials"             => tutorials_md
    ]

    mathengine = MathJax(Dict(
        :TeX => Dict(
            :equationNumbers => Dict(:autoNumber => "AMS"),
            :Macros => Dict(),
        ),
    ))

    # The format will make other pages in parallel with the index page
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        mathengine = mathengine,
        collapselevel = 1,
    )

    # This way it shows warnings of functions that have not been documented
    makedocs(
        sitename = "Radiative Transfer",
        format = format,
        clean = false,
        modules = [RadiativeTransfer],
        pages = pages,
    )

end

build()

deploydocs(
    repo = "github.com/RupeshJey/RadiativeTransfer.jl.git",
    target = "build",
    push_preview = true,
)