text:
import Base: @deprecate

@deprecate DataFrame!(args...; kwargs...) DataFrame(args...; copycols=false, kwargs...)

# TODO: remove these definitions in year 2021
by(args...; kwargs...) = throw(ArgumentError("by function was removed from DataFrames.jl. " *
                                             "Use the `combine(groupby(...), ...)` or `combine(f, groupby(...))` instead."))

aggregate(args...; kwargs...) = throw(ArgumentError("aggregate function was removed from DataFrames.jl. " *
                                                    "Use the `combine` function instead."))

export categorical, categorical!
function CategoricalArrays.categorical(df::AbstractDataFrame,
                                       cols::Union{ColumnIndex, MultiColumnIndex};
                                       compress::Union{Bool, Nothing}=nothing)
    if compress === nothing
        compress = false
        categoricalstr = "categorical"
    else
        categoricalstr = "(x -> categorical(x, compress=$compress))"
    end
    if cols isa AbstractVector{<:Union{AbstractString, Symbol}}
        Base.depwarn("`categorical(df, cols)` is deprecated. " *
                     "Use `transform(df, cols .=> $categoricalstr, renamecols=false)` instead.",
                     :categorical)
        return transform(df, cols .=> (x -> categorical(x, compress=compress)), renamecols=false)
    elseif cols isa Union{AbstractString, Symbol}
        Base.depwarn("`categorical(df, cols)` is deprecated. " *
                     "Use `transform(df, cols => $categoricalstr, renamecols=false)` instead.",
                     :categorical)
        return transform(df, cols => (x -> categorical(x, compress=compress)), renamecols=false)
    else
        Base.depwarn("`categorical(df, cols)` is deprecated. " *
                     "Use `transform(df, names(df, cols) .=> $categoricalstr, renamecols=false)` instead.",
                     :categorical)
        return transform(df, names(df, cols) .=> (x -> categorical(x, compress=compress)), renamecols=false)
    end
end

function CategoricalArrays.categorical(df::AbstractDataFrame,
                                       cols::Union{Type, Nothing}=nothing;
                                       compress::Bool=false)
    if compress === nothing
        compress = false
        categoricalstr = "categorical"
    else
        categoricalstr = "categorical(x, compress=$compress)"
    end
    if cols === nothing
        cols = Union{AbstractString, Missing}
        Base.depwarn("`categorical(df)` is deprecated. " *
                     "Use `transform(df, names(df, $cols) .=> $categoricalstr, renamecols=false)` instead.",
                     :categorical)
    else
        Base.depwarn("`categorical(df, T)` is deprecated. " *
                     "Use transform(df, names(df, T) .=> $categoricalstr, renamecols=false)` instead.",
                     :categorical)
    end
    return transform(df, names(df, cols) .=> (x -> categorical(x, compress=compress)), renamecols=false)
end

function categorical!(df::DataFrame, cols::Union{ColumnIndex, MultiColumnIndex};
                      compress::Union{Bool, Nothing}=nothing)
    if compress === nothing
        compress = false
        categoricalstr = "categorical"
    else
        categoricalstr = "(x -> categorical(x, compress=$compress))"
    end
    if cols isa AbstractVector{<:Union{AbstractString, Symbol}}
        Base.depwarn("`categorical!(df, cols)` is deprecated. " *
                     "Use `transform!(df, cols .=> $categoricalstr, renamecols=false)` instead.",
                     :categorical!)
        return transform!(df, cols .=> (x -> categorical(x, compress=compress)), renamecols=false)
    elseif cols isa Union{AbstractString, Symbol}
        Base.depwarn("`categorical!(df, cols)` is deprecated. " *
                     "Use `transform!(df, cols => $categoricalstr, renamecols=false)` instead.",
                     :categorical!)
        return transform!(df, cols => (x -> categorical(x, compress=compress)), renamecols=false)
    else
        Base.depwarn("`categorical!(df, cols)` is deprecated. " *
                     "Use `transform!(df, names(df, cols) .=> $categoricalstr, renamecols=false)` instead.",
                     :categorical!)
        return transform!(df, names(df, cols) .=> (x -> categorical(x, compress=compress)), renamecols=false)
    end
end

function categorical!(df::DataFrame, cols::Union{Type, Nothing}=nothing;
                      compress::Bool=false)
    if compress === nothing
        compress = false
        categoricalstr = "categorical"
    else
        categoricalstr = "(x -> categorical(x, compress=$compress))"
    end
    if cols === nothing
        cols = Union{AbstractString, Missing}
        Base.depwarn("`categorical!(df)` is deprecated. " *
                     "Use `transform!(df, names(df, $cols) .=> $categoricalstr, renamecols=false)` instead.",
                     :categorical!)
    else
        Base.depwarn("`categorical!(df, T)` is deprecated. " *
                     "Use `transform!(df, names(df, T) .=> $categoricalstr, renamecols=false)` instead.",
                     :categorical!)
    end
    return transform!(df, names(df, cols) .=> (x -> categorical(x, compress=compress)), renamecols=false)
end

@deprecate DataFrame(pairs::NTuple{N, Pair}; makeunique::Bool=false,
          copycols::Bool=true) where {N} DataFrame(pairs..., makeunique=makeunique, copycols=copycols)
@deprecate DataFrame(columns::NTuple{N, AbstractVector}, cnames::NTuple{N, Symbol}; makeunique::Bool=false,
          copycols::Bool=true) where {N} DataFrame(collect(columns), collect(cnames);
              makeunique=makeunique, copycols=copycols)
@deprecate DataFrame(columns::NTuple{N, AbstractVector}, cnames::NTuple{N, AbstractString}; makeunique::Bool=false,
                     copycols::Bool=true) where {N} DataFrame(collect(columns), [Symbol(c) for c in cnames];
                                                              makeunique=makeunique, copycols=copycols)
@deprecate DataFrame(columns::NTuple{N, AbstractVector};
                     copycols::Bool=true) where {N} DataFrame(collect(columns),
                                                              Symbol.(:x, 1:length(columns)), copycols=copycols)

# this deprecation is very important, becuase without it users will
# get strange results with old code as described in https://github.com/JuliaData/Tables.jl/issues/208
@deprecate DataFrame(columns::AbstractVector{<:AbstractVector}; makeunique::Bool=false,
                     copycols::Bool=true) DataFrame(columns, :auto, copycols=copycols)

@deprecate DataFrame(columns::AbstractMatrix) DataFrame(columns, :auto)

function DataFrame(column_eltypes::AbstractVector{T}, cnames::AbstractVector{Symbol},
                   nrows::Integer=0; makeunique::Bool=false)::DataFrame where T<:Type
    Base.depwarn("`DataFrame` constructor with passed eltypes is deprecated. " *
                 "Pass explicitly created columns to a `DataFrame` constructor instead.",
                     :DataFrame)
    columns = AbstractVector[elty >: Missing ?
                             fill!(Tables.allocatecolumn(elty, nrows), missing) :
                             Tables.allocatecolumn(elty, nrows)
                             for elty in column_eltypes]
    return DataFrame(columns, Index(convert(Vector{Symbol}, cnames),
                     makeunique=makeunique), copycols=false)
end

DataFrame(column_eltypes::AbstractVector{<:Type},
          cnames::AbstractVector{<:AbstractString},
          nrows::Integer=0; makeunique::Bool=false) =
    DataFrame(column_eltypes, Symbol.(cnames), nrows; makeunique=makeunique)

import Base: convert
@deprecate convert(::Type{DataFrame}, A::AbstractMatrix) DataFrame(Tables.table(A, header=Symbol.(:x, axes(A, 2))))


text:
# This code is based on the gridap hyperelasticity demo: https://gridap.github.io/Tutorials/dev/pages/t005_hyperelasticity/
# Here I expanded it to 3D and added Makie based model visualisation.

# Note this code currently requires: ] add Makie@0.15.2 GLMakie@0.4.6

using Gridap
using Gridap.Visualization
using Gridap.ReferenceFEs
using Gridap.Geometry
using FileIO
using LineSearches: BackTracking
using GLMakie, GeometryBasics
using Colors, ColorSchemes

# Geometry and BC parameters
sample_dim = [1,1,1] #Sample dimensions
numElem    = [5,5,5] #Number of elements in each direction
disp_max   = 0.3 #Maximum displacement
disp_inc   = disp_max/10 #Desired displacement increment per step
degree     = 2 #Mesh order

# Material parameters
const λ = 100.0
const μ = 1.0

# Deformation Gradient
F(∇u) = one(∇u) + ∇u'

#Jacobian = volume ratio
J(F) = sqrt(det(C(F)))

#Green-Lagrange strain
#E(F) = 0.5*( F'*F - one(F) ) #Green-Lagrange strain
dE(∇du,∇u) = 0.5*( ∇du⋅F(∇u) + (∇du⋅F(∇u))' )

# Right Cauchy-green deformation tensor
C(F) = (F')⋅F

# Hyperelastic constitutive law for the Neo hookean material
function S(∇u)
  Cinv = inv(C(F(∇u))) #Inverse of C i.e. B
  μ*(one(∇u)-Cinv) + λ*log(J(F(∇u)))*Cinv
end

function dS(∇du,∇u)
  Cinv = inv(C(F(∇u)))
  _dE = dE(∇du,∇u)
  λ*(Cinv⊙_dE)*Cinv + 2*(μ-λ*log(J(F(∇u))))*Cinv⋅_dE⋅(Cinv')
end

# Cauchy stress tensor
σ(∇u) = (1.0/J(F(∇u)))*F(∇u)⋅S(∇u)⋅(F(∇u))'

# Model
domain = (0,sample_dim[1],0,sample_dim[2],0,sample_dim[3])
partition = (numElem[1],numElem[2],numElem[3])
model = CartesianDiscreteModel(domain,partition)

# Define new boundaries
labels = get_face_labeling(model)
add_tag_from_tags!(labels,"diri_0",[1,3,5,7,13,15,17,19,25])
add_tag_from_tags!(labels,"diri_1",[2,4,6,8,14,16,18,20,26])

# Setup integration
Ω = Triangulation(model)
dΩ = Measure(Ω,degree)

# Weak form
res(u,v) = ∫( (dE∘(∇(v),∇(u))) ⊙ (S∘∇(u)) )*dΩ

jac_mat(u,du,v) =  ∫( (dE∘(∇(v),∇(u))) ⊙ (dS∘(∇(du),∇(u))) )*dΩ

jac_geo(u,du,v) = ∫( ∇(v) ⊙ ( (S∘∇(u))⋅∇(du) ) )*dΩ

jac(u,du,v) = jac_mat(u,du,v) + jac_geo(u,du,v)

# Construct the FEspace
reffe = ReferenceFE(lagrangian,VectorValue{3,Float64},1)
V = TestFESpace(model,reffe,conformity=:H1,dirichlet_tags = ["diri_0", "diri_1"])

# Setup non-linear solver
nls = NLSolver(show_trace=true,method=:newton,linesearch=BackTracking())
solver = FESolver(nls)

function run(x0,disp_x,step,nsteps,cache)

  g0 = VectorValue(0.0,0.0,0.0)
  g1 = VectorValue(disp_x,0.0,0.0)
  U = TrialFESpace(V,[g0,g1])

  #FE problem
  op = FEOperator(res,jac,U,V)

  println("\n+++ Solving for disp_x $disp_x in step $step of $nsteps +++\n")

  uh = FEFunction(U,x0)

  uh, cache = solve!(uh,solver,op,cache)

  return get_free_dof_values(uh), uh, cache

end

function runs(disp_max,disp_inc)

 nsteps = ceil(Int,abs(disp_max)/disp_inc)

 x0 = zeros(Float64,num_free_dofs(V))
 nodalDisplacements = Vector{Vector{VectorValue{3, Float64}}}(undef,nsteps+1)

 cache = nothing
 for step in 1:nsteps
  disp_x = step * disp_max / nsteps
  x0, uh, cache = run(x0,disp_x,step,nsteps,cache)

  vd = visualization_data(Ω,"",cellfields=["u"=>uh])
  nodalDisplacements[step+1] = vd[1].nodaldata["u"]
 end

 nodalDisplacements[1]=nodalDisplacements[2].*0 #Add zeros for initial state

 return nodalDisplacements
end

function nodesToPointset(V)
  #Convert gridap coordinate type to Makie Point3f type
  P=Vector{GeometryBasics.Point{3, Float32}}(undef,size(V,1))
  for q=1:1:size(V,1)
    P[q]=convert(GeometryBasics.Point,convert(Tuple{Float64, Float64, Float64},V[q]))
  end
  return P
end

function convertToFacePointSet(Ω)
  #TO DO: Implement other element types, hex->quads only shown here

  #Get gridap element and node descriptions
  # E=Ω.cell_node_ids[:] #Elements
  # V=Ω.node_coords[:] #Nodes/Vertices
  vd=visualization_data(Ω,"");
  grid = vd[1].grid
  E = get_cell_node_ids(grid)
  V = get_node_coordinates(grid)

  #Get faces and convert to QuadFace type
  F=Vector{QuadFace{Int64}}(undef,size(E,1)*6)
  for q=1:1:size(E,1)
    F[q]=convert(QuadFace{Int64},E[q][[1,2,4,3],1])             #top
    F[q+size(E,1)*1]=convert(QuadFace{Int64},E[q][[5,6,8,7],1]) #bottom
    F[q+size(E,1)*2]=convert(QuadFace{Int64},E[q][[1,2,6,5],1]) #side 1
    F[q+size(E,1)*3]=convert(QuadFace{Int64},E[q][[4,3,7,8],1]) #side 2
    F[q+size(E,1)*4]=convert(QuadFace{Int64},E[q][[2,4,8,6],1]) #front
    F[q+size(E,1)*5]=convert(QuadFace{Int64},E[q][[3,1,5,7],1]) #back
  end

  #Create face type labels
  faceTypeLabel=[ones(Int64,size(E,1))*1;
                 ones(Int64,size(E,1))*2;
                 ones(Int64,size(E,1))*3;
                 ones(Int64,size(E,1))*4;
                 ones(Int64,size(E,1))*5;
                 ones(Int64,size(E,1))*6;]

  P=nodesToPointset(V)

  return P,F, faceTypeLabel
end

#Do the work!
nodalDisplacements = runs(disp_max,disp_inc)

#Create makie compatible face and point set
pointSet,faceSet,faceTypeLabel=convertToFacePointSet(Ω)

function getCoordStep(Ω,nodalDisplacement)
  vd = visualization_data(Ω,"")
  grid = vd[1].grid
  V = get_node_coordinates(grid)
  V2= V + nodalDisplacement
  pointSet2=nodesToPointset(V2)
  return pointSet2
end

function getMagnitude(U)
  M=zeros(size(U,1))
  for q=1:1:size(U,1)
    M[q]=sqrt(U[q][1]^2 + U[q][2]^2 + U[q][3]^2)
  end
  return M
end

pointSet=getCoordStep(Ω,nodalDisplacements[1])

#Gather face and point set as GeometryBasics mesh
M=GeometryBasics.Mesh(pointSet,faceSet)

nodalColor=getMagnitude(nodalDisplacements[1])

#Visualize mesh
fig = Figure()

sl_step = Slider(fig[2, 1], range = 1:1:size(nodalDisplacements,1), startvalue = size(nodalDisplacements,1))

nodalColor = lift(sl_step.value) do stepIndex
  getMagnitude(nodalDisplacements[stepIndex])
end

M = lift(sl_step.value) do stepIndex
  GeometryBasics.Mesh(getCoordStep(Ω,nodalDisplacements[stepIndex]),faceSet)
end

titleString = lift(sl_step.value) do stepIndex
  "Step: "*string(stepIndex-1)
end

ax=Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = titleString)
hp=poly!(M, strokewidth=1,shading=false,color=nodalColor, transparency=false, overdraw=false,
colormap = (RGB(255.0, 215.0, 0.0)/255,RGB(0.0, 87.0, 183.0)/255),colorrange=(0,disp_max))
Colorbar(fig[1, 2],hp.plots[1],label = "Displacement magnitude [mm]")
fig

# ax=Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = titleString)
# hp=poly!(M, strokewidth=3,shading=true,color=nodalColor, transparency=false, overdraw=false
# ,colormap = Reverse(:Spectral),colorrange=(0,0.8))
# Colorbar(fig[1, 2],hp.plots[1],label = "Displacement magnitude [mm]")
# fig

ext: jl

max_stars_repo_name: Kevin-Mattheus-Moerman/JuliaAdventures

max_issues_repo_name: Kevin-Mattheus-Moerman/JuliaAdventures

max_forks_repo_name: Kevin-Mattheus-Moerman/JuliaAdventures


text:
# # [Basic Epidemiology Models](@id epidemiology_example)
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/examples/covid/epidemiology.ipynb)

using AlgebraicPetri.Epidemiology

using Petri
using StochasticDiffEq
using Plots

using Catlab.Theories
using Catlab.CategoricalAlgebra.ShapeDiagrams
using Catlab.Graphics

display_wd(ex) = to_graphviz(ex, orientation=LeftToRight, labels=true);

# #### SIR Model:

# define model

sir = transmission ⋅ recovery

# get resulting petri net and visualize model

p_sir = decoration(F_epi(sir));
display_wd(sir)
#-
Graph(p_sir)

# define initial states and transition rates, then
# create, solve, and visualize ODE problem

u0 = [10.0, 1, 0];
p = [0.4, 0.4];

prob,cb = SDEProblem(p_sir,u0,(0.0,7.5),p);
sol = solve(prob,SRA1(),callback=cb)

plot(sol)

# #### SEIR Model:

# define model

sei = exposure ⋅ (illness ⊗ id(I)) ⋅ ∇(I)

seir = sei ⋅ recovery

# get resulting petri net and visualize model

p_seir = decoration(F_epi(seir));

display_wd(seir)
#-
Graph(p_seir)

# define initial states and transition rates, then
# create, solve, and visualize ODE problem

u0 = [10.0, 1, 0, 0];
p = [0.9, 0.2, 0.5];

prob,cb = SDEProblem(p_seir,u0,(0.0,15.0),p);
sol = solve(prob,SRA1(),callback=cb)

plot(sol)

# #### SEIRD Model:

# define model

seird = sei ⋅ Δ(I) ⋅ (death ⊗ recovery)

# get resulting petri net and visualize model

p_seird = decoration(F_epi(seird));

display_wd(seird)
#-
Graph(p_seird)

# define initial states and transition rates, then
# create, solve, and visualize ODE problem

u0 = [10.0, 1, 0, 0, 0];
p = [0.9, 0.2, 0.5, 0.1];

prob,cb = SDEProblem(p_seird,u0,(0.0,15.0),p);
sol = solve(prob,SRA1(),callback=cb)

plot(sol)


text:
struct ShaderSource
    code::Vector{UInt8}
    language::ShaderLanguage
    stage::ShaderStageFlag
    entry_points::Vector{Symbol}
end

function pad_shader_code!(code::Vector{UInt8})
    size = cld(length(code), 4)
    rem = size * 4 - length(code)
    if rem ≠ 0
        resize!(code, size * 4)
        code[end - rem + 1:end] .= 0
    end
    @assert length(code) % 4 == 0
    code
end

struct Shader
    source::ShaderSource
    shader_module::ShaderModule
    entry_point::Symbol
    descriptor_infos::Vector{DescriptorInfo}
end

struct ShaderCache
    device::Device
    compiled::Dictionary{String,ShaderSource}
    shaders::Dictionary{ShaderSource,Shader}
end

Base.hash(source::ShaderSource, h::UInt64) = objectid(source.code) + h

ShaderCache(device) = ShaderCache(device, Dictionary{String,ShaderSource}(), Dictionary{ShaderSource,Shader}())

function find_source!(cache::ShaderCache, spec::ShaderSpecification)
    file = string(spec.source_file)
    if haskey(cache.compiled, file)
        source = cache.compiled[file]
        if spec.entry_point in source.entry_points
            return source
        end
    else
        source = ShaderSource(spec)
        if source.language ≠ SPIR_V
            source = compile(source)
        end
        insert!(cache.compiled, file, source)
        source
    end
end

function ShaderSource(spec::ShaderSpecification)
    ShaderSource(read(spec.source_file), spec.language, spec.stage, [spec.entry_point])
end

function find_shader!(cache::ShaderCache, source::ShaderSource, entry_point::Symbol)
    if source.language ≠ SPIR_V
        source = compile(source)
    end
    find_shader!(cache, source, ShaderSpecification("", false, entry_point, SPIR_V))
end

find_shader!(cache::ShaderCache, spec::ShaderSpecification) = find_shader!(cache, find_source!(cache, spec), spec)

function find_shader!(cache::ShaderCache, source::ShaderSource, spec::ShaderSpecification)
    if haskey(cache.shaders, source)
        cache.shaders[source]
    else
        shader_module = ShaderModule(cache.device, source)
        ir = IR(SPIRV.Module(IOBuffer(source.code)))
        infos = descriptor_infos(ir)
        shader = Shader(source, shader_module, spec.entry_point, infos)
        insert!(cache.shaders, source, shader)
        shader
    end
end

function descriptor_infos(ir::IR)
    filter!(!isnothing, map(pairs(ir.global_vars)) do (id, var)
        decs = var.decorations
        # only select variables with a descriptor set and binding assigned
        haskey(decs, SPIRV.DecorationDescriptorSet) && haskey(decs, SPIRV.DecorationBinding) || return nothing

        type = @match s = var.storage_class begin
            &SPIRV.StorageClassStorageBuffer => DESCRIPTOR_TYPE_STORAGE_BUFFER
            &SPIRV.StorageClassImage => @match var.type begin
                ::ImageType => DESCRIPTOR_TYPE_STORAGE_IMAGE
            end
            &SPIRV.StorageClassUniformConstant => @match var.type begin
                ::SamplerType => DESCRIPTOR_TYPE_SAMPLER
                ::SampledImageType => DESCRIPTOR_TYPE_SAMPLED_IMAGE
                ::ImageType => DESCRIPTOR_TYPE_STORAGE_IMAGE
                if haskey(decs, SPIRV.DecorationBlock) end => DESCRIPTOR_TYPE_UNIFORM_BUFFER
                _ => error(var)
            end
            &SPIRV.StorageClassUniform => @match var.type begin
                if haskey(decs, SPIRV.DecorationBufferBlock) end => DESCRIPTOR_TYPE_STORAGE_BUFFER
                _ => DESCRIPTOR_TYPE_UNIFORM_BUFFER
            end
            _ => error("Could not map variable $var to a Vulkan descriptor type.")
        end

        DescriptorInfo(type, decs[SPIRV.DecorationDescriptorSet][], decs[SPIRV.DecorationBinding][])
    end)
end


text:
module LinAlg

importall Base
import Base: USE_BLAS64, size, copy, copy_transpose!, power_by_squaring, print_matrix, transpose!

export
# Modules
    LAPACK,
    BLAS,

# Types
    SymTridiagonal,
    Tridiagonal,
    Bidiagonal,
    Woodbury,
    Factorization,
    BunchKaufman,
    Cholesky,
    CholeskyPivoted,
    Eigen,
    GeneralizedEigen,
    GeneralizedSVD,
    GeneralizedSchur,
    Hessenberg,
    LU,
    LUTridiagonal,
    LDLt,
    QR,
    QRPivoted,
    Schur,
    SVD,
    Hermitian,
    Symmetric,
    Triangular,
    Diagonal,
    UniformScaling,

# Functions
    axpy!,
    bkfact,
    bkfact!,
    chol,
    cholfact,
    cholfact!,
    cholpfact,
    cholpfact!,
    cond,
    condskeel,
    copy!,
    cross,
    ctranspose,
    det,
    diag,
    diagind,
    diagm,
    diff,
    dot,
    eig,
    eigfact,
    eigfact!,
    eigmax,
    eigmin,
    eigs,
    eigvals,
    eigvecs,
    expm,
    sqrtm,
    eye,
    factorize,
    givens,
    gradient,
    hessfact,
    hessfact!,
    ishermitian,
    isposdef,
    isposdef!,
    issym,
    istril,
    istriu,
    kron,
    ldltfact!,
    ldltfact,
    linreg,
    logdet,
    lu,
    lufact,
    lufact!,
    norm,
    null,
    peakflops,
    pinv,
    qr,
    qrfact!,
    qrfact,
    qrp,
    qrpfact!,
    qrpfact,
    rank,
    rref,
    scale,
    scale!,
    schur,
    schurfact!,
    schurfact,
    solve,
    svd,
    svdfact!,
    svdfact,
    svdvals!,
    svdvals,
    trace,
    transpose,
    tril,
    triu,
    tril!,
    triu!,
    vecnorm,

# Operators
    \,
    /,
    A_ldiv_B!,
    A_ldiv_Bc,
    A_ldiv_Bt,
    A_mul_B,
    A_mul_B!,
    A_mul_Bc,
    A_mul_Bc!,
    A_mul_Bt,
    A_mul_Bt!,
    A_rdiv_Bc,
    A_rdiv_Bt,
    Ac_ldiv_B,
    Ac_ldiv_Bc,
    Ac_mul_b_RFP,
    Ac_mul_B,
    Ac_mul_B!,
    Ac_mul_Bc,
    Ac_mul_Bc!,
    Ac_rdiv_B,
    Ac_rdiv_Bc,
    At_ldiv_B,
    At_ldiv_Bt,
    At_mul_B,
    At_mul_B!,
    At_mul_Bt,
    At_mul_Bt!,
    At_rdiv_B,
    At_rdiv_Bt,

# Constants
    I

typealias BlasFloat Union(Float64,Float32,Complex128,Complex64)
typealias BlasReal Union(Float64,Float32)
typealias BlasComplex Union(Complex128,Complex64)
typealias BlasChar Char

if USE_BLAS64
    typealias BlasInt Int64
    blas_int(x) = int64(x)
else
    typealias BlasInt Int32
    blas_int(x) = int32(x)
end

#Check that stride of matrix/vector is 1
function chkstride1(A::StridedVecOrMat...)
    for a in A
        stride(a,1)== 1 || error("Matrix does not have contiguous columns")
    end
end

#Check that matrix is square
function chksquare(A...)
    sizes=Int[]
    for a in A
        size(a,1)==size(a,2) || throw(DimensionMismatch("Matrix is not square: dimensions are $(size(a))"))
        push!(sizes, size(a,1))
    end
    length(A)==1 ? sizes[1] : sizes
end

#Check that upper/lower (for special matrices) is correctly specified
macro chkuplo()
   :((uplo=='U' || uplo=='L') || throw(ArgumentError("""invalid uplo = $uplo

Valid choices are 'U' (upper) or 'L' (lower).""")))
end

include("linalg/exceptions.jl")
include("linalg/generic.jl")

include("linalg/blas.jl")
include("linalg/matmul.jl")
include("linalg/lapack.jl")

include("linalg/dense.jl")
include("linalg/tridiag.jl")
include("linalg/factorization.jl")
include("linalg/lu.jl")

include("linalg/bunchkaufman.jl")
include("linalg/triangular.jl")
include("linalg/symmetric.jl")
include("linalg/woodbury.jl")
include("linalg/diagonal.jl")
include("linalg/bidiag.jl")
include("linalg/uniformscaling.jl")
include("linalg/rectfullpacked.jl")
include("linalg/givens.jl")
include("linalg/special.jl")
include("linalg/bitarray.jl")
include("linalg/ldlt.jl")

include("linalg/sparse.jl")
include("linalg/umfpack.jl")
include("linalg/cholmod.jl")

include("linalg/arpack.jl")
include("linalg/arnoldi.jl")

function __init__()
    Base.check_blas()
    if Base.blas_vendor() == :mkl
        ccall((:MKL_Set_Interface_Layer, Base.libblas_name), Void, (Cint,), USE_BLAS64 ? 1 : 0)
    end
end

end # module LinAlg

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
using Test

using Beamforming
using Seis

# Test array responses are calculated for a zero-slowness set of plane waves
# at 3, 4, 5, 6 and 7 Hz for three stations of the ARCES array in Norway
# and can be compared to Ruigrok et al., J Seismol, 2017.
# Data are written little-endian
@testset "Array response" begin
    @testset "$T" for T in (Float32, Float64)
        # Geographic stations
        ara1 = Station{T}(net="NO", sta="ARA1", lon=25.5071, lat=69.5363)
        ara2 = Station{T}(net="NO", sta="ARA2", lon=25.5078, lat=69.5338)
        arb2 = Station{T}(net="NO", sta="ARB2", lon=25.5134, lat=69.5357)

        stas = [ara1, ara2, arb2]

        # Geographic traces
        n = length(stas)
        t = [Trace{T}(0, 1, 0) for _ in 1:n]
        t.sta = stas

        s_max = 1
        ds = 0.01
        f1 = 3
        f2 = 7
        df = 1

        x, y, z, lon, lat = Beamforming.array_geometry(t)

        # Cartesian stations
        cstas = [Station{T,Seis.Cartesian{T}}() for _ in 1:n]
        cstas.net = stas.net
        cstas.sta = stas.sta
        cstas.x = x.*1000
        cstas.y = y.*1000
        cstas.z = z.*1000

        # Cartesian traces
        ct = [CartTrace{T}(0, 1, 0) for _ in 1:n]
        ct.sta = cstas

        @testset "Conventional ARF" begin
            true_power = ltoh.(reinterpret(Float64,
                read(joinpath(@__DIR__, "data", "arf_conventional_$(f1)_$(f2)_$(df).bin"))))
            true_power = reshape(true_power, Int(sqrt(length(true_power))), :)
            arf = array_response(t, s_max, ds, f1, f2, df)
            @test arf.power ≈ T.(true_power) rtol=5e-3
            @test arf isa Beamforming.ArrayResponse{T}

            @test array_response(stas, s_max, ds, f1, f2, df) == arf

            # Cartesian coordinates
            carf = array_response(ct, s_max, ds, f1, f2, df)
            @testset "$f" for f in fieldnames(typeof(carf))
                @test getfield(carf, f) == getfield(arf, f) ||
                    ≈(getfield(carf, f), getfield(arf, f), rtol=1e-3)
            end
            csarf = array_response(cstas, s_max, ds, f1, f2, df)
            @test carf == csarf

            @test arf.sx == -s_max:ds:s_max
            @test arf.sy == -s_max:ds:s_max
            @test arf.freqs == f1:df:f2
            @test arf.x == x
            @test arf.y == y
        end

        @testset "Cross-correlation ARF" begin
            true_power = ltoh.(reinterpret(Float64,
                read(joinpath(@__DIR__, "data", "arf_crosscorrelation_$(f1)_$(f2)_$(df).bin"))))
            true_power = reshape(true_power, Int(sqrt(length(true_power))), :)
            arf = crosscorrelation_array_response(t, s_max, ds, f1, f2, df)
            @test arf.power ≈ T.(true_power) rtol=5e-3
            @test arf isa Beamforming.ArrayResponseCrosscorrelation{T}

            @test crosscorrelation_array_response(stas, s_max, ds, f1, f2, df) == arf
            @test crosscorrelation_array_response(ct, s_max, ds, f1, f2, df) == arf
            @test crosscorrelation_array_response(cstas, s_max, ds, f1, f2, df) == arf

            α, Δ = Beamforming._azimuth_distance_ccbf(x, y)
            @test all(arf.sx .== T.(-s_max:ds:s_max))
            @test all(arf.sy .== T.(-s_max:ds:s_max))
            @test all(arf.freqs .== f1:df:f2)
            @test arf.azimuth == α
            @test arf.distance == Δ
        end
    end

    @testset "_azimuth_distance_ccbf" begin
        x = [0, 1]
        y = [1, 2]
        α = [45, 45-180]
        Δ = [√2, √2]
        @test Beamforming._azimuth_distance_ccbf(x, y) == (α, Δ)
    end
end

max_stars_repo_name: anowacki/Beamforming

max_issues_repo_name: anowacki/Beamforming

max_forks_repo_name: anowacki/Beamforming


text:
using MAT
using SparseArrays
using LinearAlgebra
include("../src/recursive_spectral.jl")
mat = matread("../data/emailEUcore.mat")
A = mat["A"]            # Largest connected component of email EU
ASim = mat["Asim"]      # Undirected version of same graph: ASim = spones(A+A')
m = sum(A)
label = vec(round.(Int64,mat["labels"]))

include("../src/Motif_Adj_Mat.jl")
Ata = ExtractTriMat(ASim)
Acc, p = largest_component(Ata)
core_inds = findall(x->x==true,p)
Ata = Ata[p,p]      # Just extracting triangles
ASim = ASim[p,p]    # Symmetric version of the graph
label = label[p]
n = size(ASim,1)

numpoints = 20
mins = 10 .^range(log10(5),log10(500),length=numpoints)
times = zeros(20,2)

SpecTri = zeros(n,numpoints)
SpecSim = zeros(n,numpoints)

for i = 1:length(mins)
    minsize = mins[i]

    # Run recursive spectral on the original graph
    start = time()
    @show minsize
    sets = recursive_spectral(ASim,minsize)
    c = set2clusvec(sets,n)
    times[i,1] = time()-start
    SpecSim[:,i] .= c

    # Run it on the motif adjacency matrix
    start = time()
    sets = recursive_spectral(Ata,minsize)
    c = set2clusvec(sets,n)
    times[i,2] = time()-start
    SpecTri[:,i] .= c

end

matwrite("Output/RecursiveSpectral_Email.mat", Dict("times"=>times,
"SpecTri"=>SpecTri, "SpecSim"=>SpecSim))

max_stars_repo_name: nveldt/ParamCC

max_issues_repo_name: nveldt/ParamCC

max_forks_repo_name: nveldt/ParamCC


text:
import Luxor: Point

isinside(p, a, b) = (b.x - a.x) * (p.y - a.y) > (b.y - a.y) * (p.x - a.x)

function intersection(a, b, s, f)
    dc = [a.x - b.x, a.y - b.y]
    dp = [s.x - f.x, s.y - f.y]
    n1 = a.x * b.y - a.y * b.x
    n2 = s.x * f.y - s.y * f.x
    n3 = 1.0 / (dc[1] * dp[2] - dc[2] * dp[1])
    Point((n1 * dp[1] - n2 * dc[1]) * n3, (n1 * dp[2] - n2 * dc[2]) * n3)
end

function clipSH(spoly, cpoly)
    outarr = spoly
    q = cpoly[end]
    for p in cpoly
        # @show p
        inarr = outarr
        outarr = Point[]
        # DEBUG if inarr is empty, I assume the two poly does not overlap
        # if isempty(inarr) return Point[] end
        s = inarr[end]
        for vtx in inarr
            if isinside(vtx, q, p)
                if !isinside(s, q, p)
                    push!(outarr, intersection(q, p, s, vtx))
                end
                push!(outarr, vtx)
            elseif isinside(s, q, p)
                push!(outarr, intersection(q, p, s, vtx))
            end
            s = vtx
        end
        q = p
    end
    outarr
end

function isoverlap_old(spoly, cpoly)
    if isempty(clipSH(spoly, cpoly)) return false
    else return true
    end
end

function test()
    isoverlap([Point(50, 150), Point(200, 50), Point(350, 150), Point(350, 300),
               Point(250, 300), Point(200, 250), Point(150, 350), Point(100, 250),
               Point(100, 200)],
              [Point(100, 100), Point(300, 100), Point(300, 300), Point(100, 300)])
    clipSH([Point(50, 150), Point(200, 50), Point(350, 150), Point(350, 300),
            Point(250, 300), Point(200, 250), Point(150, 350), Point(100, 250),
            Point(100, 200)],
           [Point(100, 100), Point(300, 100), Point(300, 300), Point(100, 300)])
    isoverlap([Point(100, 100), Point(300, 100), Point(300, 300), Point(100, 300)],
              [Point(10, 10), Point(20,10), Point(20,20), Point(10, 20)])

    isoverlap([Point(0,0), Point(0,100), Point(100,100), Point(100,0)],
              [Point(200,0), Point(200,100), Point(300,100), Point(300,0)])
    isoverlap([Point(0,0), Point(0,100), Point(100,100), Point(100,0)],
              [Point(50,-50), Point(50,50), Point(150,50), Point(150,-50)])
    clipSH([Point(0,0), Point(0,100), Point(100,100), Point(100,0)],
           [Point(50,-50), Point(50,50), Point(150,50), Point(150,-50)])
    clipSH([Point(0,0), Point(0,100), Point(100,100), Point(100,0)],
           [Point(-50,30), Point(-50,50), Point(150,50)])

    clipSH([Point(100,150), Point(200,250), Point(300,200)],
           [Point(150,150), Point(150,200),
            Point(200,200), Point(200,150)])
end

# and there's also PolygonClipping.jl
# https://github.com/JuliaGeometry/PolygonClipping.jl
#
# Clipper.jl https://github.com/JuliaGeometry/Clipper.jl
#
# which is a wrapper for Clipper library
# http://www.angusj.com/delphi/clipper.php
#
# some algorithms
# https://en.wikipedia.org/wiki/Sutherland%E2%80%93Hodgman_algorithm
#
# But this snippet does not work (shown above):
# https://rosettacode.org/wiki/Sutherland-Hodgman_polygon_clipping#Julia
using Clipper
function test()
    path1 = Vector{IntPoint}()
    push!(path1, IntPoint(0, 0))
    push!(path1, IntPoint(0, 1))
    push!(path1, IntPoint(1, 1))
    push!(path1, IntPoint(1, 0))

    path2 = Vector{IntPoint}()
    push!(path2, IntPoint(1, 0))
    push!(path2, IntPoint(1, 1))
    push!(path2, IntPoint(2, 1))
    push!(path2, IntPoint(2, 0))

    c = Clip()
    add_path!(c, path1, PolyTypeSubject, true)
    add_path!(c, path2, PolyTypeSubject, true)

    result, polys = execute(c, ClipTypeUnion, PolyFillTypeEvenOdd, PolyFillTypeEvenOdd)
    # @test result == true
    # @test polys[1][1] == Clipper.IntPoint(0, 0)
    # @test polys[1][2] == Clipper.IntPoint(2, 0)
    # @test polys[1][3] == Clipper.IntPoint(2, 1)
    # @test polys[1][4] == Clipper.IntPoint(0, 1)
end

function isoverlap_clip(poly0, polys...)
    path1 = Vector{IntPoint}()
    for p in poly0
        push!(path1, p)
    end
    c = Clip()
    add_path!(c, path1, PolyTypeSubject, true)

    for poly in polys
        path2 = Vector{IntPoint}()
        for p in poly
            push!(path2, p)
        end
        add_path!(c, path2, PolyTypeClip, true)
    end

    result, polys = execute(c,
                            # ClipTypeUnion,
                            # ClipTypeDifference,
                            ClipTypeIntersection,
                            PolyFillTypeEvenOdd,
                            PolyFillTypeEvenOdd)
    # @show polys
    # FIXME isempty(polys) does not work
    return length(polys) > 0
end

function tuple2intpoint(t)
    map(t) do (x,y)
        IntPoint(round(x*1e6),round(y*1e6))
    end
end

function isoverlap(poly0, polys...)
    isoverlap_clip(tuple2intpoint(poly0), map(tuple2intpoint, polys)...)
end

function test()
    isoverlap([(0,0), (0, 1), (1, 1), (1, 0)],
              [(1, 0), (1, 1), (2, 1), (2, 0)])
    isoverlap([(0,0), (0, 1), (1, 1), (1, 0)],
              [(1, 0), (1, 1), (2, 1), (2, 0)])

    # true, overlap
    isoverlap([(0,0), (0,100), (100,100), (100,0)],
              [(50,-50), (50,50), (150,50), (150,-50)])
    isoverlap([(0,0), (0,100), (100,100), (100,0)],
              [(-50,30), (-50,50), (150,50)])
    isoverlap([(100,150), (200,250), (300,200)],
              [(150,150), (150,200),
               (200,200), (200,150)])
    # false, non-overlapping
    isoverlap([(0,0), (0,100), (100,100), (100,0)],
              [(200,0), (200,100), (300,100), (300,0)])
end

max_stars_repo_name: lihebi/bhdl

max_issues_repo_name: lihebi/bhdl

max_forks_repo_name: lihebi/bhdl


text:
# Note, this requires a special branch of Zygote to count ops
using Zygote, Flux

function build_model(num_layers::Int, size::Int)
    layers = Any[]
    for layer_idx in 1:num_layers
        push!(layers, LSTM(size, size))
    end
    return mapleaves(Flux.data, Chain(layers...))
end

for num_layers in 1:4
	model = build_model(num_layers, 4)

	seq_len = 4
	x = randn(4, 4)
	test(m, x) = Zygote.gradient(m, x) do m, x
		x_t = x
		for idx in 1:seq_len
			x_t = m(x_t)
		end
		return sum(x_t)
	end
	test(model, x)
	Zygote.reset_num_ops!()
	test(model, x)
	@info "ops for layers", num_layers, Zygote.get_num_ops()
end

max_stars_repo_name: MikeInnes/zygote-paper

max_issues_repo_name: MikeInnes/zygote-paper

max_forks_repo_name: MikeInnes/zygote-paper


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    RegularGridData(data)
    RegularGridData(data, origin, spacing)

Regularly spaced `data` georeferenced with `origin` and `spacing`.
The `data` argument is a dictionary mapping variable names to Julia
arrays with the actual data.

`NaN` or `missing` values in the Julia arrays are interpreted as
non-valid. They can be used to mask the variables on the grid.

## Examples

Given `poro` and `perm` two 2-dimensional Julia arrays containing
values of porosity and permeability, the following code can be used
to georeference the data:

```julia
julia> data = OrderedDict(:porosity => poro, :permeability => perm)
julia> RegularGridData(data, (0.,0.,0.), (1.,1.,1.))
```

Alternatively, one can omit `origin` and `spacing` for default
values of zeros and ones:

```julia
julia> RegularGridData{Float64}(data)
```

See also: [`RegularGrid`](@ref)
"""
struct RegularGridData{T,N} <: AbstractData{T,N}
  data::OrderedDict{Symbol,<:AbstractArray}
  domain::RegularGrid{T,N}

  function RegularGridData{T,N}(data, domain) where {N,T}
    new(data, domain)
  end
end

function RegularGridData(data::OrderedDict{Symbol,<:AbstractArray},
                         origin::NTuple{N,T}, spacing::NTuple{N,T}) where {N,T}
  sizes = [size(array) for array in values(data)]
  @assert length(unique(sizes)) == 1 "data dimensions must be the same for all variables"
  @assert length(sizes[1]) == N "inconsistent number of dimensions for given origin/spacing"
  RegularGridData{T,length(origin)}(data, RegularGrid(sizes[1], origin, spacing))
end

RegularGridData{T}(data::OrderedDict{Symbol,<:AbstractArray{<:Any,N}}) where {N,T} =
  RegularGridData(data, ntuple(i->zero(T), N), ntuple(i->one(T), N))

Base.size(geodata::RegularGridData) = size(geodata.domain)
origin(geodata::RegularGridData) = origin(geodata.domain)
spacing(geodata::RegularGridData) = spacing(geodata.domain)

function Base.getindex(geodata::RegularGridData, var::Symbol)
  vals = [getindex(geodata, ind, var) for ind in 1:npoints(geodata)]
  reshape(vals, size(geodata))
end

function Base.getindex(geodata::RegularGridData,
                       icoords::Vararg{Int,N}) where {N}
  lin  = LinearIndices(size(geodata.domain))
  ind  = lin[icoords...]
  vars = [var for (var,V) in variables(geodata)]
  vals = [geodata.data[var][ind] for var in vars]
  NamedTuple{tuple(vars...)}(vals)
end

# ------------
# IO methods
# ------------
function Base.show(io::IO, geodata::RegularGridData{T,N}) where {N,T}
  dims = join(size(geodata), "×")
  print(io, "$dims RegularGridData{$T,$N}")
end

function Base.show(io::IO, ::MIME"text/plain", geodata::RegularGridData{T,N}) where {N,T}
  println(io, geodata)
  println(io, "  origin:  ", Tuple(origin(geodata)))
  println(io, "  spacing: ", Tuple(spacing(geodata)))
  println(io, "  variables")
  varlines = ["    └─$var ($(eltype(array)))" for (var, array) in geodata.data]
  print(io, join(varlines, "\n"))
end


text:
const GITHUB_REPO_ISSUES_URL = "https://github.com/felipenoris/BLPData.jl/issues"

Base.showerror(io::IO, err::BLPResultErrException{SecurityErr}) = print(io, err.err.security, " error: ", err.err.message)
Base.showerror(io::IO, err::BLPUnknownException) = print(io, "$(err.cause). Please, send a bug report to $GITHUB_REPO_ISSUES_URL with the full log of this error.")
Base.showerror(io::IO, err::BLPTimeoutException) = print(io, "Timeout: exceeded timeout limit of $(err.timeout_milliseconds) milliseconds.")

function BLPUnknownException(cause::AbstractString)
    BLPException(BLPUnknownErr(String(cause)))
end

function FieldErr(security::AbstractString, element::Element{false, BLPAPI_DATATYPE_SEQUENCE})
    @assert has_name(element, "fieldExceptions")
    error_info = element["errorInfo"]
    return FieldErr(
            String(security),
            Symbol(get_element_value(element["fieldId"])),
            get_element_value(error_info["source"]),
            get_element_value(error_info["code"]),
            get_element_value(error_info["category"]),
            get_element_value(error_info["subcategory"]),
            get_element_value(error_info["message"])
        )
end

get_result(r::NamedTupleResult) = r.result
get_result(r::FieldDataVecResult) = r.field_data_vec
unwrap(::Unwrap, r::BLPResultErr) = throw(BLPResultErrException(r))
unwrap(::Unwrap, r::BLPResultOk) = get_result(r)
unwrap(::NoUnwrap, r::BLPResult) = r

function parse_reference_data_response_into!(
            element::Element{false, BLPAPI_DATATYPE_CHOICE},
            result::Dict{String, Any},
            securities::Vector{T1},
            fields::Vector{T2},
            error_handling::ErrorHandling
        ) where {T1<:AbstractString, T2<:AbstractString}

    if !has_name(element, "ReferenceDataResponse")
        throw(BLPUnknownException("Expected response element with name `ReferenceDataResponse`. Got `$(get_name(element))`."))
    end

    response_element = get_choice(element)

    if has_name(response_element, "responseError")
        throw(BLPResponseException("Got responseError. \n$response_element"))
    end

    if !has_name(response_element, "securityData")
        throw(BLPUnknownException("Expected response_element with name `securityData`. Got `$(get_name(response_element))`."))
    end

    parse_security_data_into!(response_element, result, securities, fields, error_handling)
end

function parse_historical_data_response_into!(
            element::Element{false, BLPAPI_DATATYPE_CHOICE},
            result::Dict{String, Any},
            security::AbstractString,
            fields::Vector{T},
            error_handling::ErrorHandling
        ) where {T<:AbstractString}

    if !has_name(element, "HistoricalDataResponse")
        throw(BLPUnknownException("Expected response element with name `HistoricalDataResponse`. Got `$(get_name(element))`."))
    end

    response_element = get_choice(element)

    if has_name(response_element, "responseError")
        throw(BLPResponseException("Got responseError. \n$response_element"))
    end

    if !has_name(response_element, "securityData")
        throw(BLPUnknownException("Expected response_element with name `securityData`. Got `$(get_name(response_element))`."))
    end

    parse_security_data_into!(response_element, result, String(security), fields, error_handling)
end

function parse_security_data_into!(
            security_data_vec::Element{true, BLPAPI_DATATYPE_SEQUENCE},
            result::Dict{String, Any},
            securities::Vector{T1},
            fields::Vector{T2},
            error_handling::ErrorHandling
        ) where {T1<:AbstractString, T2<:AbstractString}

    for security_data_element in get_element_value(security_data_vec)
        security = get_element_value(security_data_element["security"])

        if !(security ∈ securities)
            throw(BLPUnknownException("Security $security not in the requested securities list $securities."))
        end

        parse_security_data_into!(security_data_element, result, security, fields, error_handling)
    end

    nothing
end

function parse_security_data_into!(security_data_element::Element{false, BLPAPI_DATATYPE_SEQUENCE}, result::Dict{String, Any}, security::String,
            fields::Vector{T2},
            error_handling::ErrorHandling
        ) where {T1<:AbstractString, T2<:AbstractString}

    @assert Symbol(get_name(security_data_element)) == :securityData
    @assert isa(security_data_element, Element{false, BLPAPI_DATATYPE_SEQUENCE})

    if haskey(result, security)
        throw(BLPUnknownException("Got repeated response for security $security."))
    end

    if haskey(security_data_element, "securityError")
        result[security] = unwrap(error_handling, SecurityErr(security, security_data_element["securityError"]))
    else
        result[security] = parse_field_data(security_data_element["fieldData"], security_data_element["fieldExceptions"], security, fields, error_handling)
    end

    nothing
end

function parse_field_data(field_data::Element{false, BLPAPI_DATATYPE_SEQUENCE}, field_exceptions::Element{true, BLPAPI_DATATYPE_SEQUENCE}, security::String, fields::Vector{T}, error_handling::E) where {T<:AbstractString, E<:ErrorHandling}
    data = parse_field_exceptions(field_exceptions, security, error_handling)

    for child_element in each_child_element(field_data)
        data[Symbol(get_name(child_element))] = get_element_value(child_element)
    end

    return unwrap(error_handling, NamedTupleResult(fields, data))
end

function parse_field_data(field_data_vec::Element{true, BLPAPI_DATATYPE_SEQUENCE}, field_exceptions::Element{true, BLPAPI_DATATYPE_SEQUENCE}, security::String, fields::Vector{T}, error_handling::E) where {T<:AbstractString, E<:ErrorHandling}

    field_exceptions_dict = parse_field_exceptions(field_exceptions, security, error_handling)

    result = Vector()
    data_buffer = Dict{Symbol, Any}()

    for field_data in get_element_value(field_data_vec)
        empty!(data_buffer)

        for child_element in each_child_element(field_data)
            data_buffer[Symbol(get_name(child_element))] = get_element_value(child_element)
        end

        push!(result, unwrap(error_handling, NamedTupleResult(fields, data_buffer)))
    end

    return unwrap(error_handling, FieldDataVecResult(field_exceptions_dict, result))
end

function parse_field_exceptions(field_exceptions_vec::Element{true, BLPAPI_DATATYPE_SEQUENCE}, security::AbstractString, error_handling::ErrorHandling) where {T}
    result = Dict{Symbol, Any}()

    for field_exception in get_element_value(field_exceptions_vec)
        field_sym = Symbol(get_element_value(field_exception["fieldId"]))
        result[field_sym] = unwrap(error_handling, FieldErr(security, field_exception))
    end

    return result
end

function NamedTupleResult(fields::Vector{T}, data::Dict{Symbol, A}) where {T<:AbstractString, A}
    tuple_keys = Symbol.(fields)
    tuple_values = Vector(undef, length(tuple_keys))
    for (i, k) in enumerate(tuple_keys)
        if haskey(data, k)
            tuple_values[i] = data[k]
        else
            tuple_values[i] = missing
        end
    end

    # trick based on the docstring for NamedTuple
    nt = (; zip(tuple_keys, tuple_values)...)
    return NamedTupleResult(nt)
end

function SecurityErr(security::AbstractString, err_element::Element{false, BLPAPI_DATATYPE_SEQUENCE})
    @assert has_name(err_element, "securityError")

    return SecurityErr(
            security,
            get_element_value(err_element["source"]),
            get_element_value(err_element["code"]),
            get_element_value(err_element["category"]),
            get_element_value(err_element["subcategory"]),
            get_element_value(err_element["message"]),
        )
end

max_stars_repo_licenses: OpenSSL MIT

max_issues_repo_licenses: OpenSSL MIT

max_forks_repo_licenses: OpenSSL MIT


text:
# Utilities
using LinearAlgebra: checksquare, opnorm

include("SymplecticForm.jl")

# calculate the 'nonSymplecticity".
nonSymplecticity(A::AbstractMatrix) = opnorm(A'*Ω*A - Ω) / opnorm(A)

# check if a matrix is symplectic, returning its non-symplecticity
function checkSymplectic(A::AbstractMatrix)
    A <: Real || throw(ArgumentError("matrix is not real-valued"))
    checksquare(A) % 2 == 0 || throw(DimensionMismatch("matrix is not even-dimensional"))
    return nonSymplecticity(A)
end
function checkSymplecticMatrix(A...)
    nonSymplecticities = []
    for a in A
        a <: Real || throw(ArgumentError("matrix is not real-valued"))
        checksquare(a) % 2 == 0 || throw(DimensionMismatch("matrix is not even-dimensional"))
        push!(nonSymplecticities, nonSymplecticity(a))
    end
    return nonSymplecticities
end

using LinearAlgebra

# Define the AbstractSymplecticMatrix data type
abstract type AbstractSymplecticMatrix{T} <: AbstractArray{T, 2} end
Base.size(A::AbstractSymplecticMatrix) = error("size not implemented")
Base.getindex(A::AbstractSymplecticMatrix, i::Int) = error("getindex not implemented")
Base.getindex(A::AbstractSymplecticMatrix, I::Vararg{Int, N}) where {N} = error("getindex not implemented")
Base.setindex!(A::AbstractSymplecticMatrix, v, i::Int) = error("setindex! not implemented")
Base.setindex!(A::AbstractSymplecticMatrix, v, I::Vararg{Int, N}) where {N} = error("setindex! not implemented")
# get a Q-quadrature from columns
getQCol(A::AbstractSymplecticMatrix, i::Int) = error("getQCol not implemented")
# get a P-quadrature from columns
getPCol(A::AbstractSymplecticMatrix, i::Int) = error("getPCol not implemented")
# get a Q-quadrature from rows
getQRow(A::AbstractSymplecticMatrix, i::Int) = error("getQCol not implemented")
# get a P-quadrature from rows
getPRow(A::AbstractSymplecticMatrix, i::Int) = error("getPCol not implemented")

Base.transpose(S::AbstractSymplecticMatrix) = error("transpose not implemented")
Base.adjoint(S::AbstractSymplecticMatrix) = error("adjoin not implemented")

Base.:-(S::AbstractSymplecticMatrix) = error("- not implemented")
Base.:+(S::AbstractSymplecticMatrix) = error("+ not implemented")
Base.:*(S1::AbstractSymplecticMatrix, S2::AbstractSymplecticMatrix, Ss::Vararg{AbstractSymplecticMatrix}) = error("* not implemented")
Base.inv(S::AbstractSymplecticMatrix) = error("inv not implemented")
Base.:^(S::AbstractSymplecticMatrix, p::Int) = error("^ not implemented")

# Define the SymplecticMatrix data type
struct SymplecticMatrix{T} <: AbstractSymplecticMatrix{T}
    matrix:: AbstractArray{T, 2}
    function SymplecticMatrix{T}(matrix:: AbstractArray{T, 2}) where T
        if !(T<:Real) return error("not real") end
        if !(LinearAlgebraUtilities.isSymplectic(matrix)) return error("not symplectic") end
        return new{T}(AbstractMatrix{T}(matrix))
    end
end

Base.size(S::SymplecticMatrix) = size(S.matrix)
Base.getindex(S::SymplecticMatrix, i::Int) = getindex(S.matrix, i)
Base.getindex(S::SymplecticMatrix, I::Vararg{Int, N}) where {N} = getindex(S.matrix, I)
getQCol(S::SymplecticMatrix, i::Int) = S.matrix[:, 2*i-1]
getPCol(S::SymplecticMatrix, i::Int) = S.matrix[:, 2*i]
getQRow(S::SymplecticMatrix, i::Int) = S.matrix[2*i-1, :]
getPRow(S::SymplecticMatrix, i::Int) = S.matrix[2*i, :]

SymplecticMatrix(M::AbstractArray{T, 2}) where T = SymplecticMatrix{T}(M)
SymplecticMatrix(S::SymplecticMatrix) = S

Base.transpose(S::SymplecticMatrix)::SymplecticMatrix = transpose( S.matrix )
Base.adjoint(S::SymplecticMatrix)::SymplecticMatrix = adjoint( S.matrix )

Base.:-(S::SymplecticMatrix)::SymplecticMatrix = Base.:-( S.matrix )
Base.:+(S::SymplecticMatrix)::SymplecticMatrix = Base.:+( S.matrix )
Base.:*(S1::SymplecticMatrix, S2::SymplecticMatrix, Ss::Vararg{SymplecticMatrix})::SymplecticMatrix = *(S1.matrix, S2.matrix,  [s.matrix for s in Ss]... )
Base.inv(S::SymplecticMatrix)::SymplecticMatrix = - Ω(S) * S' * Ω(S)
Base.:^(S::SymplecticMatrix, p::Int)::SymplecticMatrix = p < 0 ? inv(S).matrix ^ (-p)  : S.matrix^p

# Unary Minus and Unary Plus
⊗ = Base.kron

# Direct Sum
dsum(M1::AbstractMatrix, M2::AbstractMatrix)::AbstractMatrix = LinearAlgebraUtilities.dsum(M1, M2)
dsum(S1::SymplecticMatrix, S2::SymplecticMatrix, Ss::Vararg{SymplecticMatrix})::SymplecticMatrix = dsum(S1.matrix, S2.matrix,  [s.matrix for s in Ss]... )
⊕ = dsum

Base.Matrix(S::SymplecticMatrix)::Array{T, 2} = S.matrix
Base.convert(::Type{SymplecticMatrix}, S::SymplecticMatrix) = S
Base.convert(::Type{SymplecticMatrix}, M::AbstractArray{T, 2}) where T = SymplecticMatrix(M)
Base.convert(::Type{AbstractArray}, S::SymplecticMatrix) = Matrix(S)

# Omega(n::Int)::SymplecticMatrix = LinearAlgebraUtilities.Omega(n)
# Omega(S::SymplecticMatrix)::SymplecticMatrix = LinearAlgebraUtilities.Omega(S.matrix)
# Ω = Omega

isGeneric(S::SymplecticMatrix)::Bool = LinearAlgebraUtilities.isGeneric(S.matrix)

# Override functions in Base
# Base.Matrix(S::SymplecticMatrix)::AbstractMatrix = S.matrix
# Base.convert(::Type{SymplecticMatrix}, matrix::AbstractMatrix) = SymplecticMatrix(matrix)
# Base.convert(::Type{SymplecticMatrix}, S::SymplecticMatrix) = S
# Base.convert(::Type{AbstractMatrix}, S::SymplecticMatrix) = Matrix(S)
# Base.size(S::SymplecticMatrix, dim::Integer) = size(S.matrix, dim)
# Base.size(S::SymplecticMatrix) = size(S.matrix)
## Import !!!
# function Base.getindex(S::SymplecticMatrix{T}, i::Integer, j::Integer)::T where {T}
#     return S.matrix[i, j]
# end


text:
println("\nTesting the Land ODE functions...");
@testset "Land --- Land model" begin
    # initialize parameters
    @test true;
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
module Lasso
using Compat
module Util
    # Extract fields from object into function locals
    # See https://github.com/JuliaLang/julia/issues/9755
    macro extractfields(from, fields...)
        esc(Expr(:block, [:($(fields[i]) = $(from).$(fields[i])) for i = 1:length(fields)]...))
    end
    export @extractfields
end

import Base.LinAlg.BlasReal
include("FusedLasso.jl")
include("TrendFiltering.jl")

using Reexport, StatsBase, .Util
@reexport using GLM, Distributions, .FusedLassoMod, .TrendFiltering
using GLM.FPVector, GLM.wrkwt!
export LassoPath, fit, fit!, coef

## HELPERS FOR SPARSE COEFFICIENTS

immutable SparseCoefficients{T} <: AbstractVector{T}
    coef::Vector{T}              # Individual coefficient values
    coef2predictor::Vector{Int}  # Mapping from indices in coef to indices in original X
    predictor2coef::Vector{Int}  # Mapping from indices in original X to indices in coef

    SparseCoefficients(n::Int) = new(T[], Int[], zeros(Int, n))
end

function Base.A_mul_B!{T}(out::Vector, X::Matrix, coef::SparseCoefficients{T})
    fill!(out, zero(eltype(out)))
    @inbounds for icoef = 1:nnz(coef)
        ipred = coef.coef2predictor[icoef]
        c = coef.coef[icoef]
        @simd for i = 1:size(X, 1)
            out[i] += c*X[i, ipred]
        end
    end
    out
end

function Base.A_mul_B!{T}(out::Vector, X::SparseMatrixCSC, coef::SparseCoefficients{T})
    @extractfields X colptr rowval nzval
    fill!(out, zero(eltype(out)))
    @inbounds for icoef = 1:nnz(coef)
        ipred = coef.coef2predictor[icoef]
        c = coef.coef[icoef]
        @simd for i = colptr[ipred]:colptr[ipred+1]-1
            out[rowval[i]] += c*nzval[i]
        end
    end
    out
end

function Base.dot{T}(x::Vector{T}, coef::SparseCoefficients{T})
    v = 0.0
    @inbounds @simd for icoef = 1:nnz(coef)
        v += x[coef.coef2predictor[icoef]]*coef.coef[icoef]
    end
    v
end

Base.size(x::SparseCoefficients) = (length(x.predictor2coef),)
Base.nnz(x::SparseCoefficients) = length(x.coef)
Base.getindex{T}(x::SparseCoefficients{T}, ipred::Int) =
    x.predictor2coef[ipred] == 0 ? zero(T) : x.coef[x.predictor2coef[ipred]]

function Base.setindex!{T}(A::Matrix{T}, coef::SparseCoefficients, rg::UnitRange{Int}, i::Int)
    A[:, i] = zero(T)
    for icoef = 1:nnz(coef)
        A[rg[coef.coef2predictor[icoef]], i] = coef.coef[icoef]
    end
    A
end

function Base.copy!(x::SparseCoefficients, y::SparseCoefficients)
    length(x) == length(y) || throw(DimensionMismatch())
    n = length(y.coef)
    resize!(x.coef, n)
    resize!(x.coef2predictor, n)
    copy!(x.coef, y.coef)
    copy!(x.coef2predictor, y.coef2predictor)
    copy!(x.predictor2coef, y.predictor2coef)
    x
end

# Add a new coefficient to x, returning its index in x.coef
function addcoef!{T}(x::SparseCoefficients{T}, ipred::Int)
    push!(x.coef, zero(T))
    push!(x.coef2predictor, ipred)
    coefindex = nnz(x)
    x.predictor2coef[ipred] = coefindex
end

# Add newcoef to column i of coefs
function addcoefs!(coefs::SparseMatrixCSC, newcoef::SparseCoefficients, i::Int)
    n = nnz(coefs)
    nzval = coefs.nzval
    rowval = coefs.rowval
    resize!(nzval, n+nnz(newcoef))
    resize!(rowval, n+nnz(newcoef))
    @inbounds for ipred = 1:length(newcoef.predictor2coef)
        icoef = newcoef.predictor2coef[ipred]
        if icoef != 0
            cval = newcoef.coef[icoef]
            if cval != 0
                n += 1
                nzval[n] = cval
                rowval[n] = ipred
            end
        end
    end
    resize!(nzval, n)
    resize!(rowval, n)
    coefs.colptr[i+1:end] = n+1
end

## COEFFICIENT ITERATION IN SEQUENTIAL OR RANDOM ORDER

if VERSION >= v"0.4-dev+1915"
    # Julia 0.4 has a nice interface that lets us do random coefficient
    # iteration quickly.
    immutable RandomCoefficientIterator
        rng::MersenneTwister
        rg::Base.Random.RangeGeneratorInt{Int,@compat UInt}
        coeforder::Vector{Int}
    end
    const RANDOMIZE_DEFAULT = true

    RandomCoefficientIterator() =
        RandomCoefficientIterator(MersenneTwister(1337), Base.Random.RangeGenerator(1:2), Int[])
else
    immutable RandomCoefficientIterator
        RandomCoefficientIterator() = error("randomization not supported on Julia 0.3")
    end
    const RANDOMIZE_DEFAULT = false
end

typealias CoefficientIterator @compat Union{UnitRange{Int},RandomCoefficientIterator}

# Iterate over coefficients in random order
function Base.start(x::RandomCoefficientIterator)
    if !isempty(x.coeforder)
        @inbounds for i = length(x.coeforder):-1:2
            j = rand(x.rng, x.rg)
            x.coeforder[i], x.coeforder[j] = x.coeforder[j], x.coeforder[i]
        end
    end
    return 1
end
Base.next(x::RandomCoefficientIterator, i) = (x.coeforder[i], i += 1)
Base.done(x::RandomCoefficientIterator, i) = i > length(x.coeforder)

# Add an additional coefficient and return a new CoefficientIterator
function addcoef(x::RandomCoefficientIterator, icoef::Int)
    push!(x.coeforder, icoef)
    RandomCoefficientIterator(x.rng, Base.Random.RangeGenerator(1:length(x.coeforder)), x.coeforder)
end
addcoef(x::UnitRange{Int}, icoef::Int) = 1:length(x)+1

## LASSO PATH

type LassoPath{S<:@compat(Union{LinearModel,GeneralizedLinearModel}),T} <: RegressionModel
    m::S
    nulldev::T                    # null deviance
    nullb0::T                     # intercept of null model, if one was fit
    λ::Vector{T}                  # shrinkage parameters
    autoλ::Bool                   # whether λ is automatically determined
    Xnorm::Vector{T}              # original squared norms of columns of X before standardization
    pct_dev::Vector{T}            # percent deviance explained by each model
    coefs::SparseMatrixCSC{T,Int} # model coefficients
    b0::Vector{T}                 # model intercepts
    niter::Int                    # number of coordinate descent iterations

    LassoPath(m, nulldev::T, nullb0::T, λ::Vector{T}, autoλ::Bool, Xnorm::Vector{T}) =
        new(m, nulldev, nullb0, λ, autoλ, Xnorm)
end

function Base.show(io::IO, path::LassoPath)
    prefix = isa(path.m, GeneralizedLinearModel) ? string(typeof(path.m.rr.d).name.name, " ") : ""
    println(io, prefix*"Lasso Solution Path ($(size(path.coefs, 2)) solutions for $(size(path.coefs, 1)) predictors in $(path.niter) iterations):")

    coefs = path.coefs
    ncoefs = zeros(Int, size(coefs, 2))
    for i = 1:size(coefs, 2)-1
        ncoefs[i] = coefs.colptr[i+1] - coefs.colptr[i]
    end
    ncoefs[end] = nnz(coefs) - coefs.colptr[size(coefs, 2)] + 1
    Base.showarray(io, [path.λ path.pct_dev ncoefs]; header=false)
end

StatsBase.coef(path::LassoPath) = path.coefs

## MODEL CONSTRUCTION

# Controls early stopping criteria with automatic λ
const MIN_DEV_FRAC_DIFF = 1e-5
const MAX_DEV_FRAC = 0.999

# Compute automatic λ values based on X'y and λminratio
function computeλ(Xy, λminratio, α, nλ)
    λmax = abs(Xy[1])
    for i = 1:length(Xy)
        x = abs(Xy[i])
        λmax = ifelse(x > λmax, x, λmax)
    end
    λmax /= α
    logλmax = log(λmax)
    λ = exp(linspace(logλmax, logλmax + log(λminratio), nλ))
end

function StatsBase.fit{T<:AbstractFloat,V<:FPVector}(::Type{LassoPath},
                                                     X::AbstractMatrix{T}, y::V, d::UnivariateDistribution=Normal(),
                                                     l::Link=canonicallink(d);
                                                     wts::@compat(Union{FPVector,Void})=ones(T, length(y)),
                                                     offset::V=similar(y, 0),
                                                     α::Number=one(eltype(y)), nλ::Int=100,
                                                     λminratio::Number=ifelse(size(X, 1) < size(X, 2), 0.01, 1e-4),
                                                     λ::@compat(Union{Vector,Void})=nothing, standardize::Bool=true,
                                                     intercept::Bool=true,
                                                     naivealgorithm::Bool=(!isa(d, Normal) || !isa(l, IdentityLink) || size(X, 2) > 5*size(X, 1)),
                                                     dofit::Bool=true,
                                                     irls_tol::Real=1e-7, randomize::Bool=RANDOMIZE_DEFAULT,
                                                     maxncoef::Int=min(size(X, 2), 2*size(X, 1)), fitargs...)
    size(X, 1) == size(y, 1) || DimensionMismatch("number of rows in X and y must match")
    n = length(y)
    length(wts) == n || error("length(wts) = $(length(wts)) should be 0 or $n")

    # Standardize predictors if requested
    if standardize
        Xnorm = vec(full(std(X, 1, corrected=false)))
        for i = 1:length(Xnorm)
            @inbounds Xnorm[i] = 1/Xnorm[i]
        end
        X = scale(X, Xnorm)
    else
        Xnorm = T[]
    end

    # Lasso initialization
    α = convert(T, α)
    λminratio = convert(T, λminratio)
    coefitr = randomize ? RandomCoefficientIterator() : (1:0)
    cd = naivealgorithm ? NaiveCoordinateDescent{T,intercept,typeof(X),typeof(coefitr)}(X, α, maxncoef, 1e-7, coefitr) :
                          CovarianceCoordinateDescent{T,intercept,typeof(X),typeof(coefitr)}(X, α, maxncoef, 1e-7, coefitr)

    # GLM response initialization
    autoλ = λ == nothing
    wts .*= convert(T, 1/sum(wts))
    off = convert(Vector{T}, offset)

    if isa(d, Normal) && isa(l, IdentityLink)
        # Special no-IRLS case
        mu = isempty(offset) ? y : y + off
        nullb0 = intercept ? mean(mu, weights(wts)) : zero(T)
        nulldev = 0.0
        @simd for i = 1:length(mu)
            @inbounds nulldev += abs2(mu[i] - nullb0)*wts[i]
        end

        if autoλ
            # Find max λ
            if intercept
                muscratch = Array(T, length(mu))
                @simd for i = 1:length(mu)
                    @inbounds muscratch[i] = (mu[i] - nullb0)*wts[i]
                end
            else
                muscratch = mu.*wts
            end
            Xy = X'muscratch
            λ = computeλ(Xy, λminratio, α, nλ)
        else
            λ = convert(Vector{T}, λ)
        end

        # First y is just a placeholder here
        model = LinearModel(LmResp{typeof(y)}(mu, off, wts, y), cd)
    else
        # Fit to find null deviance
        # Maybe we should reuse this GlmResp object?
        nullmodel = fit(GeneralizedLinearModel, ones(T, n, ifelse(intercept, 1, 0)), y, d, l;
                        wts=wts, offset=offset, convTol=irls_tol)
        nulldev = deviance(nullmodel)

        if autoλ
            # Find max λ
            Xy = X'*broadcast!(*, nullmodel.rr.wrkresid, nullmodel.rr.wrkresid, nullmodel.rr.wrkwts)
            λ = computeλ(Xy, λminratio, α, nλ)
            nullb0 = intercept ? coef(nullmodel)[1] : zero(T)
        else
            λ = convert(Vector{T}, λ)
            nullb0 = zero(T)
        end

        eta = GLM.initialeta!(d, l, similar(y), y, wts, off)
        rr = GlmResp{typeof(y),typeof(d),typeof(l)}(y, d, l, eta, similar(eta), offset, wts)
        model = GeneralizedLinearModel(rr, cd, false)
    end

    # Fit path
    path = LassoPath{typeof(model),T}(model, nulldev, nullb0, λ, autoλ, Xnorm)
    dofit && fit!(path; irls_tol=irls_tol, fitargs...)
    path
end

include("coordinate_descent.jl")

end


text:
using TuringModels

Turing.setadbackend(:reverse_diff)
#nb Turing.turnprogress(false)

# Turing model

@model m8_3(y) = begin
    α ~ Normal(1, 10)
    σ ~ Truncated(Cauchy(0, 1), 0, Inf)

    for i ∈ 1:length(y)
        y[i] ~ Normal(α, σ)
    end
end

y = [-1,1]

# Sample

chns = sample(m8_3(y), NUTS(0.95), 1000);

# Results rethinking

m83rethinking = "
      mean   sd  5.5% 94.5% n_eff Rhat
alpha 0.09 1.63 -2.13  2.39   959    1
sigma 2.04 2.05  0.68  4.83  1090    1
";

# Describe the posterior samples

describe(chns)

# End of `08/m8.3t.jl`


text:
###############################
#
#        SECOND ORDER
#
###############################

# Thanks to https://web.media.mit.edu/~crtaylor/calculator.html
# For different stencil sizes:
# c1 = [1, -2, 1]
# c2 = [-1, 16, -30, 16, -1] / 12
# c3 = [2, -27, 270, -490, 270, -27, 2] / 180
# c4 = [-9, 128, -1008, 8064, -14350, 8064, …] / 5040
# etc

function fd_diff2_ord2!(dx2_U, U, dx⁻¹)
    local c1 = SA[1.0, -2.0, 1.0]
    for i = 2 : length(U)-1
        dx2_U[i] = dot(c1, U[i-1 : i+1])
    end
    dx2_U[1] = dx2_U[2]
    dx2_U[end] = dx2_U[end-1]

    @. dx2_U *= dx⁻¹ ^ 2

    nothing
end

function fd_diff2_ord4!(dx2_U, U, dx⁻¹)
    local c1 = SA[1.0, -2.0, 1.0]
    local c2 = SA[-1.0, 16.0, -30.0, 16.0, -1.0] * 0.08333333333333333333
    for i in 3 : length(U)-2
        dx2_U[i] = dot(c2, U[i-2 : i+2])
    end

    dx2_U[2] = dot(c1, U[1:3])
    dx2_U[end-1] = dot(c1, U[end-2 : end])

    dx2_U[1] = dx2_U[2]
    dx2_U[end] = dx2_U[end-1]

    @. dx2_U *= dx⁻¹ ^ 2

    nothing
end

function fd_diff2_ord6!(dx2_U, U, dx⁻¹)
    local c1 = SA[1.0, -2.0, 1.0]
    local c2 = SA[-1.0, 16.0, -30.0, 16.0, -1.0] * 0.08333333333333333333
    local c3 = SA[2.0, -27.0, 270.0, -490.0, 270.0, -27.0, 2.0]
    c3 *= 0.00555555555555555555
    for i in 4 : length(U)-3
        dx2_U[i] = dot(c3, U[i-3 : i+3])
    end

    dx2_U[3] = dot(c2, U[1:5])
    dx2_U[end-2] = dot(c2, U[end-4 : end])

    dx2_U[2] = dot(c1, U[1:3])
    dx2_U[end-1] = dot(c1, U[end-2 : end])

    dx2_U[1] = 2.0*dx2_U[2] - dx2_U[3]
    dx2_U[end] = 2.0*dx2_U[end-1] - dx2_U[end-2]

    @. dx2_U *= dx⁻¹ ^ 2

    nothing
end

max_stars_repo_name: tremelow/StiffKinetic

max_issues_repo_name: tremelow/StiffKinetic

max_forks_repo_name: tremelow/StiffKinetic


text:
@symbol_func function K_BR(cur_reactor::AbstractReactor)
  cur_nu_n = cur_reactor.nu_n

  cur_nu_T = cur_reactor.nu_T

  cur_K = 0.1056

  cur_K *= ( 1 + cur_nu_n ) ^ 2

  cur_K *= sqrt( 1 + cur_nu_T )

  cur_K /= 1 + 2 * cur_nu_n + cur_nu_T / 2

  cur_K *= cur_reactor.Z_eff

  cur_K *= cur_reactor.epsilon ^ 2

  cur_K *= cur_reactor.kappa_95

  cur_K *= g(cur_reactor)

  cur_K
end


text:
"""
    Conv(kernel, in=>out, σ; stride = 1, pading = 0, set_w = "Xavier")
This is the traditional convolution layer. `kernel` is a tuple of integers that specifies the kernel size, it must have one or two elements. And, `in` and `out` specifies number of input and out channels.

The input data must have a dimensions WHCB(weight, width, channel, batch). If you want to use a data which has a dimentions WHC, you must be add a dimentioins of B.

`stride` and `padding` are single integers or tuple(stride is tuple of 2 elements, padding is tuple of 2 elements), and if you specifies `same` to padding, we adjust sizes of input and return a matrix which has the same sizes.
`set_w` is `Xavier` or `He`, it decide a method to create a first parameter. This parameter is the same as `Dense()`.

# Example
```
julia> C = Conv((2, 2), 2=>2, relu)
Convolution(k:(2, 2), IO:2 => 2, σ:relu)

julia> C(rand(10, 10, 2, 5)) |> size
(9, 9, 2, 5)
```

!!! warning
    When you specidies `same` to `padding`, in some cases, it will be returned one size smaller.
    Because of its expression.
```

julia> C = Conv((2, 2), 2=>2, relu, padding = "same")
Convolution(k:(2, 2), IO:2 => 2, σ:relu

julia> C(rand(10, 10, 2, 5)) |> size
(9, 9, 2, 5)
```
"""
struct Conv <: Param
    σ
    w
    W
    b
    stride::Tuple
    padding::Union{Tuple, String}
    shape::IOShape
end

function Conv(k::NTuple{N, Int}, io::Pair{<:Integer, <:Integer}, σ; stride = 1, padding = 0, set_w = "Xavier") where {N}
    weight = conv_w(k, io, set_w)
    M = io[2]
    w = reshape(weight, M, div(length(weight), M))
    n = ndims(weight)
    padding = convert(2N, padding)
    stride = convert(N, stride)
    b = create_bias(weight, size(weight, n))
    return Conv(σ, w, weight, b, stride, padding, IOShape())
end

function Base.show(io::IO, C::Conv)
    k = size(C.W)[1:end-2]
    i, o = size(C.W)[end-1:end]
    σ = C.σ
    print("Convolution(k:$k, IO:$i => $o, σ:$σ)")
end

function (C::Conv)(x::AbstractArray)
    IC = Im2Col(x, size(C.W)[1:2], C.stride, C.padding)
    if size(C.w)[2] != size(IC.x)[1]
        kc = size(C.w)[end-1]
        ic = IC.c
        throw(DimensionMismatch("the number of channels must match between the kernel and the input!(now kernel: $kc, input: $ic)"))
    end
    σ, b = C.σ, C.b
    M = size(C.W)[end]
    x = reshape(σ.(C.w * IC.x .+ b), IC.Oh, IC.Ow, M, IC.b)
    return x
end

function (C::Conv)(x::AbstractArray, record::Dict, i)
    Ih, Iw = size(x)[1:2]
    IC = Im2Col(x, size(C.W)[1:2], C.stride, C.padding)
    C.shape(Ih, Iw, IC.Oh, IC.Ow)
    if size(C.w)[2] != size(IC.x)[1]
        kc = size(C.w)[end-1]
        ic = IC.c
        throw(DimensionMismatch("the number of channels must match between the kernel and the input!(now kernel: $kc, input: $ic)"))
    end
    σ, b = C.σ, C.b
    M = size(C.W)[end]
    z = C.w * IC.x
    _z = reshape(z, IC.Oh, IC.Ow, M, IC.b)
    a = σ.(z .+ b)
    a = reshape(a, IC.Oh, IC.Ow, M, IC.b)
    record[i+1] = (_z, a)
    return a
end

function (C::Conv)(Δ, z, back::Bool)
    z = Im2Col(z, size(C.W)[1:2], C.stride, C.padding).x
    w, σ = C.w, C.σ
    Δ = w' * Δ .* σ.(z, true)
    Fh, Fw,  = size(C.W)
    Δ = Col2Im(Δ, size(C.W)[1:2], C.shape, C.stride, C.padding).x
    #I have to use Col2Im to Δ.
    return Δ
end


text:
@testset "runtestRepeat" begin
    println("Testing runtestRepeat")
    include("../src/bfgaRepeat.jl")
    using .bfgaRepeat

    include("../src/Bf.jl")
    using .BfInterpreter

    @time model = test_serial()
    #@time model = test_parallel()
    prems = model.population[length(model.population)]
    println("Generation : $(model.params.currentGeneration)")
    println("& fitness : $(model.specific_fitness.fitness(prems, model.instructionsSet))")
    #bfcode :: Array{Float64,1}  = prems.dna
    println("code :  $(prems.program)") #$(join(genesToBfInstructions(bfcode), "") ) ")

    #bft = bfType(prems.dna)
    input1 = "One Two One Two, repeat after me"
    input2 = "The RMA had been infected"
    n1 , n2 = length(input1), length(input2)

    output1, _ = execute(prems.program, input1, model.instructionsSet)
    output2, _ = execute(prems.program, input2, model.instructionsSet)

    @test length(output1) >= n1
    @test collect(input1[1:n1]) == output1[1:n1]
    @test length(output2) >= n2
    @test collect(input2[1:n2]) == output2[1:n2]

    # do the test_parallel after ...
    #test_parallel()
end


text:
using BinDeps
using BinDeps: MakeTargets

basedir = dirname(@__FILE__)
prefix = joinpath(basedir, "usr")

cmake_version = v"3.7.2"
base_url = "https://cmake.org/files/v$(cmake_version.major).$(cmake_version.minor)"
@static if Sys.iswindows()
    binary_name = "cmake.exe"
else
    binary_name = "cmake"
end

function install_binaries(file_base, file_ext, binary_dir)
    filename = "$(file_base).$(file_ext)"
    url = "$(base_url)/$(filename)"
    binary_path = joinpath(basedir, "downloads", file_base, binary_dir)

    @static if Sys.iswindows()
        install_step = () -> begin
            for dir in readdir(dirname(binary_path))
                cp(joinpath(dirname(binary_path), dir),
                   joinpath(prefix, dir);
                   force=true)
            end
        end
    else
        install_step = () -> begin
            for file in readdir(binary_path)
                symlink(joinpath(binary_path, file),
                        joinpath(prefix, "bin", file))
            end
        end
    end

    function test_step()
        try
            run(`$(joinpath(prefix, "bin", binary_name)) --version`)
        catch e
            error("""
Running the precompiled cmake binary failed with the error
$(e)
To build from source instead, run:
    julia> ENV["CMAKE_JL_BUILD_FROM_SOURCE"] = 1
    julia> Pkg.build("CMake")
""")
        end
    end
    (@build_steps begin
        FileRule(joinpath(prefix, "bin", binary_name),
            (@build_steps begin
                FileDownloader(url, joinpath(basedir, "downloads", filename))
                FileUnpacker(joinpath(basedir, "downloads", filename),
                             joinpath(basedir, "downloads"),
                             "")
                CreateDirectory(joinpath(prefix, "bin"))
                install_step
                test_step
            end))
    end)
end

function install_from_source(file_base, file_ext)
    filename = "$(file_base).$(file_ext)"
    url = "$(base_url)/$(filename)"

    (@build_steps begin
        FileRule(joinpath(prefix, "bin", binary_name),
            (@build_steps begin
                FileDownloader(url, joinpath(basedir, "downloads", filename))
                CreateDirectory(joinpath(basedir, "src"))
                FileUnpacker(joinpath(basedir, "downloads", filename),
                             joinpath(basedir, "src"),
                             "")
                begin
                    ChangeDirectory(joinpath(basedir, "src", file_base))
                    `./configure --prefix=$(prefix)`
                    MakeTargets()
                    MakeTargets("install")
                end
            end))
    end)
end

force_source_build = lowercase(get(ENV, "CMAKE_JL_BUILD_FROM_SOURCE", "")) in ["1", "true"]

process = @static if Sys.islinux()
    if Sys.ARCH == :x86_64 && !force_source_build
        install_binaries(
            "cmake-$(cmake_version)-Linux-x86_64",
            "tar.gz",
            "bin")
    else
        install_from_source("cmake-$(cmake_version)", "tar.gz")
    end
elseif Sys.isapple()
    if !force_source_build
        install_binaries(
            "cmake-$(cmake_version)-Darwin-x86_64",
            "tar.gz",
            joinpath("CMake.app", "Contents", "bin"))
    else
        install_from_source("cmake-$(cmake_version)", "tar.gz")
    end
elseif Sys.iswindows()
    if sizeof(Int) == 8
        install_binaries(
            "cmake-$(cmake_version)-win64-x64",
            "zip",
            "bin")
    elseif sizeof(Int) == 4
        install_binaries(
            "cmake-$(cmake_version)-win32-x86",
            "zip",
            "bin")
    else
        error("Only 32- or 64-bit architectures are supported")
    end
else
    error("Sorry, I couldn't recognize your operating system.")
end

run(process)

open(joinpath(dirname(@__FILE__), "deps.jl"), "w") do f
    write(f, """
const cmake = "$(escape_string(joinpath(prefix, "bin", binary_name)))"
""")

end


text:
if isdefined(OffsetArrays, :OffsetMatrix)
    # This alias is introduced in OffsetArrays v1.1.3
    OffsetMatrix = OffsetArrays.OffsetMatrix
else
    const OffsetMatrix{T} = OffsetArrays.OffsetArray{T,2}
end

if VERSION < v"1.1"
    # newly added in https://github.com/JuliaLang/julia/pull/30630
    function require_one_based_indexing(A...)
        return !Base.has_offset_axes(A...) || throw(
            ArgumentError(
                "offset arrays are not supported but got an array with index other than 1"
            ),
        )
    end
else
    using Base: require_one_based_indexing
end


text:
"""
Created on Sun Apr 26 2020

@author: Yoann Pradat

    CentraleSupelec
    MICS laboratory
    9 rue Juliot Curie, Gif-Sur-Yvette, 91190 France

Code for the functions cophcorr, dispersion and silhouette.
"""

"""
    _build_lists_ancestors!(list_nodes, lists_ancestors, node, list_ancestors, merges)

Helper recursive function for `_get_lists_ancestors`. Recursively go through a binary tree to get the list of ancestors
of each node.
"""
function _build_lists_ancestors!(list_nodes::Vector{Int64}, lists_ancestors::Vector{Vector{Int64}}, node::Int64, list_ancestors::Vector{Int64}, merges::Array{Int64,2})
    if node < 0
        push!(list_nodes, node)
        push!(lists_ancestors, list_ancestors)

        nothing
    else
        push!(list_ancestors, node)

        node_left = merges[node,1]
        node_right = merges[node,2]

        #### recursive search
        _build_lists_ancestors!(list_nodes, lists_ancestors, node_left, list_ancestors[:], merges)
        _build_lists_ancestors!(list_nodes, lists_ancestors, node_right, list_ancestors[:], merges)
    end
end

"""
    _get_lists_ancestors(merges)

Helper function for `_cophenetic_distances`. Retrieve the list of ancestors of each node in a binary tree.
"""
function _get_lists_ancestors(merges::Array{Int64, 2})
    n = size(merges, 1)
    list_nodes      = Vector{Int64}()
    lists_ancestors = Vector{Vector{Int64}}()

    node_left = merges[n,1]
    node_right = merges[n,2]
    list_ancestors  = Vector{Int64}()

    _build_lists_ancestors!(list_nodes, lists_ancestors, node_left, list_ancestors[:], merges)
    _build_lists_ancestors!(list_nodes, lists_ancestors, node_right, list_ancestors[:], merges)

    list_nodes, lists_ancestors
end

"""
    _get_lowest_common_ancestor(ancestors_a, ancestors_b)

Find lowest common element from two vectors ordered in decreasing order.
"""
function _get_lowest_common_ancestor(ancestors_a::Vector{Int64}, ancestors_b::Vector{Int64})
    n_a = length(ancestors_a)
    n_b = length(ancestors_b)

    if n_a == 0 || n_b == 0
        #### -1 to mean the tree's root node
        return -1
    end

    #### iterate backwards
    for i = n_a:-1:1
        for j = n_b:-1:1
            if ancestors_b[j] > ancestors_a[i]
                break
            elseif ancestors_b[j] == ancestors_a[i]
                return ancestors_a[i]
            end
        end
    end

    #### no common ancestor was found
    return -1
end

"""
    _compute_cophenetic_distances(merges, heights)

Compute matrix of cophenetic distances from the array of merges and the array of heights as produced in output of the hclust function.
"""
function _compute_cophenetic_distances(merges::Array{Int64, 2}, heights::Array{Float64,1})
    #### get the list of ancestors for each node
    list_nodes, lists_ancestors = _get_lists_ancestors(merges)

    #### init matrix of distances
    n = length(list_nodes)
    D = zeros(Float64, n, n)

    for i = 2:n
        for j=1:i-1
            node_i = -list_nodes[i]
            node_j = -list_nodes[j]

            lca = _get_lowest_common_ancestor(lists_ancestors[i], lists_ancestors[j])

            if lca == -1
                D[node_i, node_j] = heights[n-1]
            else
                D[node_i, node_j] = heights[lca]
            end
        end
    end

    D + transpose(D)
end

"""
    cons_cophenetic(cons, hc)

Computes the cophenetic correlation coefficient from the consensus matrix and hierarchical clustering result.
"""
function cons_cophenetic(cons::Matrix{T}, hc::Hclust{T}) where T <: Real
    #### matrix of distances
    D_con = 1 .- cons

    #### compute matrix of cophenetic distances
    D_cop = _compute_cophenetic_distances(hc.merges, hc.heights)

    #### compute cophenetic correlation
    D_con_lower_part = D_con[[x > y ? true : false for x in 1:size(D_con,1), y in 1:size(D_con,2)]]
    D_cop_lower_part = D_cop[[x > y ? true : false for x in 1:size(D_cop,1), y in 1:size(D_cop,2)]]

    #### Pearson correlation
    cor(D_con_lower_part, D_cop_lower_part)
end

"""
    cons_dispersion(cons)

Compute the dispersion coefficient for the consensus matrix.
"""
function cons_dispersion(cons::Matrix{T}) where T <: Real
    mean(4 * (cons .- 0.5) .^ 2)
end

"""
    cons_avg_silhouette(cons, hc, rank)

Compute the average silhouette width from the hierarchical clustering of the consensus matrix cut with the number
of clusters equal to the rank.
"""
function cons_avg_silhouette(cons::Matrix{T}, hc::Hclust{T}, rank::Integer) where T <: Real
    idxs = cutree(hc, k=rank)
    silh = silhouettes(idxs, 1 .- cons)
    silh_avg = 0

    for lab in 1:rank
        silh_avg += mean(silh[idxs .== lab])
    end

    silh_avg/rank
end

"""
    matr_avg_col_sparseness(X)

Compute the average sparseness of the columns of the input matrix using Hoyer's definition of sparseness.
"""
function matr_avg_col_sparseness(X::Matrix{T}) where T <: Real
    mean(mapslices(x -> (sqrt(length(x)) - norm(x, 1)/norm(x,2))/(sqrt(length(x))-1), X; dims=1))
end

"""
    matr_avg_row_sparseness(X)

Compute the average sparseness of the rows of the input matrix using Hoyer's definition of sparseness.
"""
function matr_avg_row_sparseness(X::Matrix{T}) where T <: Real
    mean(mapslices(x -> (sqrt(length(x)) - norm(x, 1)/norm(x,2))/(sqrt(length(x))-1), X; dims=2))
end


text:
# Implement parallel external sort for 4-byte integers
#
# For more information https://en.wikipedia.org/wiki/External_sorting

@everywhere include("es.jl")
@everywhere include("phelper.jl")

"""
External sort.
"""
function parallel_external_sort(number_of_machines,
                                max_memory,
                                input_filename="sample.bin",
                                output_filename="sorted-sample.bin")
    files_on_machines = Array(Array{UTF8String}, number_of_machines)
    for i in 1:number_of_machines
        files_on_machines[i] = Array(UTF8String, 0)
    end

    # 1. Get the number of integers because this allow to send more data
    #    over network at each call.
    number_of_integers = 0
    input_file = open(input_filename, "r")
    while ~ eof(input_file)
        number_of_integers += 1
        read(input_file, Int32)
    end
    close(input_file)

    # 2. Send data over network.
    input_file = open(input_filename, "r")
    @sync while ~ eof(input_file)
        for machine in 1:number_of_machines
            temporary_vector = Array(Int32, 0)
            # The loop is needed because otherwise it can raise
            # "EOFError: read end of file"
            for i in 1:max_memory
                if eof(input_file)
                    break
                else
                    push!(temporary_vector, read(input_file, Int32))
                end
            end

            # Array without elements will raise
            # UndefRefError: access to undefined reference
            if length(temporary_vector) == 0
                break
            end

            push!(files_on_machines[machine],
                  string("machine-", machine, "-", length(files_on_machines[machine]) + 1 ,"-sample.bin"))
            remotecall(machine,
                       write_chunck_to_file,
                       temporary_vector,
                       files_on_machines[machine][end],
                       true,
                       false)
        end
    end
    close(input_file)

    # 3. Sort files on remote machine
    #    using the single machine version.
    @sync for machine in 1:number_of_machines
        remotecall(machine,
                        external_sort_phase2,
                        max_memory,
                        files_on_machines[machine],
                        string("machine-", machine, "-sorted-sample.bin"))
    end

    # 4. Merge files locally.
    output_file = open(output_filename, "w")
    size_of_sorted_chunk = floor(Integer, max_memory / (number_of_machines + 1))
    input_remote_files_position = zeros(Int32, number_of_machines)  # We use -1 when the file ended.
    input_buffer = Array(Array{Int32}, number_of_machines)

    for machine in 1:number_of_machines
        # If using shared array we can't change the data.
        shared_array, input_remote_files_position[machine] = remotecall_fetch(
                machine,
                parallel_get_new_chunk,
                string("machine-", machine, "-sorted-sample.bin"),
                input_remote_files_position[machine],
                size_of_sorted_chunk)
        input_buffer[machine] = copy(shared_array)
    end

    while ~ input_buffer_is_empty(input_buffer)
        selected_machine = get_first_non_empty_input_buffer(input_buffer)
        for machine in selected_machine + 1:number_of_machines
            if length(input_buffer[machine]) > 0 && input_buffer[machine][1] < input_buffer[selected_machine][1]
                selected_machine = machine
            end
        end
        write(output_file, shift!(input_buffer[selected_machine]))

        if length(input_buffer[selected_machine]) == 0
            shared_array, input_remote_files_position[selected_machine] = remotecall_fetch(
                    selected_machine,
                    parallel_get_new_chunk,
                    string("machine-", selected_machine, "-sorted-sample.bin"),
                    input_remote_files_position[selected_machine],
                    size_of_sorted_chunk)
            input_buffer[selected_machine] = copy(shared_array)
        end
    end
    close(output_file)

    return
end

max_stars_repo_name: rgaiacs/external-sort

max_issues_repo_name: rgaiacs/external-sort

max_forks_repo_name: rgaiacs/external-sort


text:
module C
const NAMES = [
    "S"
    "kE"
    "kM"
    "kCD"
    "kCDS"
    "kR"
    "kRE"
    "kb"
    "KS"
    "kCE"
    "dM"
    "dE"
    "dCD"
    "dCE"
    "dR"
    "dRP"
    "dRE"
    "kP1"
    "kP2"
    "kDP"
    "KM"
    "KE"
    "KCD"
    "KCE"
    "KRP"
]

#name2idx
for (idx,name) in enumerate(NAMES)
    eval(Meta.parse("const $name = $idx"))
end

const NUM = length(NAMES)
end  # module

max_stars_repo_name: himoto/Bifurcation

max_issues_repo_name: himoto/Bifurcation

max_forks_repo_name: himoto/Bifurcation


text:
using Hilbert
using Test

# Example from https://www.mathworks.com/help/signal/ref/hilbert.html
@testset "Real vector" begin
    signal = [1 2 3 4]
    analytical_signal = hilbert(signal)
    @test analytical_signal == [1+1im 2-1im 3-1im 4+im]
end

@testset "Complex numbers" begin
    complex_signal = Complex.([1 2 3 4])
    # verify that the warning is thrown
    @test_logs (:warn, "Using real part, ignoring complex part") hilbert(complex_signal)
    # verify that the output is correct assuming only the real part is used
    analytical_signal = hilbert(complex_signal)
    @test analytical_signal == [1+1im 2-1im 3-1im 4+im]
end


text:
using DataFrames
using CSV

cd("/Users/kiristern/Documents/GitHub/SDM/")

predators = CSV.read("data/preds/originals/Procyon_lotor.csv", delim="\t")
predators = select(predators, [:species, :decimalLatitude, :decimalLongitude, :year])

# select only ursus arctos from imported all predators csv file
#df = predators[predators.species .== "Ursus arctos", :]

p1 = CSV.read("data/prey_sp/Bufo_bufo.csv")
p1 = select(p1,[:genus, :decimalLatitude, :decimalLongitude, :year])
rename!(p1, :genus => :species)

p2 = CSV.read("data/prey_sp/Caiman_yacare.csv", delim="\t")
p2 = select(p2,[:species, :decimalLatitude, :decimalLongitude, :year])

p3 = CSV.read("data/prey_sp/Emberiza_citrinella.csv", delim="\t")
p3 = select(p3,[:species, :decimalLatitude, :decimalLongitude, :year])

p4 = CSV.read("data/prey_sp/Larus_canus.csv", delim="\t")
p4 = select(p4,[:species, :decimalLatitude, :decimalLongitude, :year])

#create one dataframe with predators/prey dataframes
pred_prey = vcat(predators, p2)

pred_prey = select(pred_prey, [:species, :decimalLatitude, :decimalLongitude, :year])

pred_prey = dropmissing!(pred_prey, [:decimalLatitude, :decimalLongitude, :year])

CSV.write("data/pred_prey/Panthera_onca.csv", pred_prey)

show(by(pred_prey, :species, :species => length), allcols=true)

max_stars_repo_name: kiristern/SDM

max_issues_repo_name: kiristern/SDM

max_forks_repo_name: kiristern/SDM


text:
using RocksDB
using Test

# write your own tests here
db = RocksDB.open_db("/tmp/test.db", true)
a = Array{Int32}(1000)
for i in 1:1000
    a[i] = rand(Int32)
    RocksDB.db_put(db, string(i), a[i])
end

# Check if values are correct
for i in 1:1000
    val = RocksDB.db_get(db, string(i))
    @test val == a[i]
end