text:
export graph_horner, graph_horner_degopt

"""
     (graph,crefs)=graph_horner(a; input=:A, B=:B, C=:C, scaling=1.0)

Generates the graph for the polynomial using Horner's method. More precisely, it
corresponds to the evaluation of the polynomial

    p(s) = a[1] + a[2]*(αs) + ... + a[n-1]*(αs)^(n-2) + a[n]*(αs)^(n-1),

as the recursion

    Cj=s*B{j+1}
    Bj=a[j]*I+α*Cj,

where `α=scaling`.

The kwargs `B` and `C` specifies the base-names of these intermediate variables.
"""
function graph_horner(a; input = :A, B = :B, C = :C, scaling = 1.0)

    # Initial setup
    n = length(a)
    cref = Vector{Tuple{Symbol,Int}}(undef, n)
    graph = Compgraph(eltype(a))
    out = Symbol("$(B)1")
    d = n - 1 # Polynomial degree

    if d == 0
        add_lincomb!(graph, out, a[1], :I, 0.0, :I)
        add_output!(graph, out)
        cref[1] = (out, 1)
        return (graph, cref)
    end

    for j = d:-1:1
        Bj = Symbol("$B$j")
        Bjp1 = Symbol("$B$(j+1)")
        Cj = Symbol("$C$j")
        if (j == d)
            add_lincomb!(graph, Bj, a[j], :I, a[j+1] * scaling, input)
            cref[j+1] = (Bj, 1)
            cref[j] = (Bj, 2)
        else
            #   Cj=X*B{j+1}
            #   Bj=a[j]*I+α*Cj
            add_mult!(graph, Cj, input, Bjp1)
            add_lincomb!(graph, Bj, a[j], :I, scaling, Cj)
            cref[j] = (Bj, 1)
        end
    end
    push!(graph.outputs, Symbol("$(B)1"))
    return (graph, cref)
end

"""
     (graph,crefs)=graph_horner_degopt(a; scaling=1.0, input=:A)

Generates a polynomial using Horner's evaluation scheme. The polynomial

     p(s) = a[1] + a[2]*(αs) + ... + a[n-1]*(αs)^(n-2) + a[n]*(αs)^(n-1),

is evaluated as

     p(s) = a[1] + (αs)*(a[2] + (αs)*(... + (αs)*(a[n-1] + a[n]*(αs))...)),

where α=`scaling`.
However, the function uses a call to [`graph_degopt`](@ref), resulting in more
degrees of freedom in `crefs`. See also [`graph_horner`](@ref).
"""
function graph_horner_degopt(a; scaling = 1.0, input = :A)
    n = length(a)
    T = eltype(a)

    if n == 1
        error("Does not implement degree-zero polynomial.")
    end
    if n == 2
        return graph_degopt([], [a[1], a[2] * scaling])
    end

    x = Vector{Tuple{Vector{T},Vector{T}}}(undef, n - 2)
    x[1] = ([a[n-1], a[n] * scaling], [zero(T), scaling])
    for i = 2:n-2
        x[i] = (
            vcat(a[n-i], zeros(T, i - 1), one(T)),
            vcat(zero(T), scaling, zeros(T, i - 1)),
        )
    end

    return graph_degopt(x, vcat(a[1], zeros(T, n - 2), one(T)), input = input)
end


text:
"""
    randstate(basis)

Calculate a random normalized ket state.
"""
function randstate(b::Basis)
    psi = Ket(b, rand(ComplexF64, length(b)))
    normalize!(psi)
    psi
end

"""
    randoperator(b1[, b2])

Calculate a random unnormalized dense operator.
"""
randoperator(b1::Basis, b2::Basis) = DenseOperator(b1, b2, rand(ComplexF64, length(b1), length(b2)))
randoperator(b::Basis) = randoperator(b, b)

"""
    thermalstate(H,T)

Thermal state ``exp(-H/T)/Tr[exp(-H/T)]``.
"""
function thermalstate(H::AbstractOperator,T::Real)
    return normalize(exp(-dense(H)/T))
end

"""
    coherentthermalstate(basis::FockBasis,H,T,alpha)

Coherent thermal state ``D(α)exp(-H/T)/Tr[exp(-H/T)]D^†(α)``.
"""
function coherentthermalstate(basis::B,H::AbstractOperator{B,B},T::Real,alpha::Number) where B<:FockBasis
    return displace(basis,alpha)*thermalstate(H,T)*dagger(displace(basis,alpha))
end

"""
    phase_average(rho)

Returns the phase-average of ``ρ`` containing only the diagonal elements.
"""
function phase_average(rho::DenseOperator)
    return DenseOperator(basis(rho),diagm(0 => diag(rho.data)))
end

"""
    passive_state(rho,IncreasingEigenenergies::Bool=true)

Passive state ``π`` of ``ρ``. IncreasingEigenenergies=true means that higher indices correspond to higher energies.
"""
function passive_state(rho::DenseOperator,IncreasingEigenenergies::Bool=true)
    return DenseOperator(basis(rho),diagm(0 => sort!(abs.(eigvals(rho.data)),rev=IncreasingEigenenergies)))
end

ext: jl


text:
using RoamesGeometry
VERSION < v"0.7" ? using Base.Test : using Test
using StaticArrays
using CoordinateTransformations
using Rotations

include("Line.jl")
include("LineString.jl")
include("Quadratic.jl")
include("Catenary.jl")

include("AbstractRegion.jl")
include("BoundingBox.jl")
include("TriangularPrism.jl")
include("Cylinder.jl")
include("Sphere.jl")
include("Triangle.jl")
include("Circle.jl")
include("Polygon.jl")

include("distances.jl")
include("wkt.jl")

if VERSION >= v"0.7.0"
	using Random
	using AcceleratedArrays
	using Colors
	using FixedPointNumbers
	using TypedTables

	include("GridIndex.jl")
	include("pointcloud_io.jl")
end


text:
@testset "Cell" begin

@testset "Constructors and macros" begin
    @test @inferred(Cell{:A}(1)) == Cell{:A, Int}(1)

    @test (@Cell A::Int64 = 1) == Cell{:A, Int}(1)
    @test (@Cell A = 1) == Cell{:A, Int}(1)
    @test @Cell(B::Float64 = 2.0) == Cell{:B, Float64}(2.0)
    @test @Cell(B = 2.0) == Cell{:B, Float64}(2.0)
end

@testset "Introspection" begin
    c1 = @Cell A::Int = 1

    @test name(c1) == :A
    @test name(typeof(c1)) == :A
    @test eltype(c1) == Int
    @test eltype(typeof(c1)) == Int
    @test show(c1) == nothing
    println()
    @test length(c1) == 1
    #@test length(typeof(c1)) == 1

    @test @inferred(rename(c1, Val{:A_new})) == @Cell(A_new = 1)
end

@testset "Accessing and interating" begin
    c1 = @Cell A::Int = 1

    @test copy(c1) == c1
    @test c1[] == c1.data
    @test c1[1] == c1.data
    @test c1[Val{:A}] == c1.data
    @test @inferred(first(c1)) == c1.data
    @test endof(c1) == c1.data
end

end


text:
module RetroSignalModel

export RtgMTK, scan_params, load_conditions, load_parameters, load_data
export optim_params, optim_params_threads, find_steady_states

include("common.jl")
include("models.jl")
include("params.jl")
include("steadystates.jl")

end


text:
using Sundials, Test, LinearAlgebra

# create mass matrix problems
function make_mm_probs(mm_A, ::Type{Val{iip}}) where {iip}
    # iip
    mm_b = vec(sum(mm_A; dims = 2))
    function mm_f(du, u, p, t)
        LinearAlgebra.mul!(du, mm_A, u)
        du .+= t * mm_b
        nothing
    end
    mm_g(du, u, p, t) = (@. du = u + t; nothing)

    # oop
    mm_f(u, p, t) = mm_A * (u .+ t)
    mm_g(u, p, t) = u .+ t

    mm_analytic(u0, p, t) = @. 2 * u0 * exp(t) - t - 1

    u0 = ones(3)
    tspan = (0.0, 1.0)

    prob = ODEProblem(
        ODEFunction{iip, true}(mm_f, analytic = mm_analytic, mass_matrix = mm_A),
        u0,
        tspan,
    )
    prob2 = ODEProblem(ODEFunction{iip, true}(mm_g, analytic = mm_analytic), u0, tspan)

    prob, prob2
end

mm_A = Float64[-2 1 4; 4 -2 1; 2 1 3]
prob, prob2 = make_mm_probs(mm_A, Val{true})

sol = solve(prob, ARKODE(), abstol = 1e-8, reltol = 1e-8)
sol2 = solve(prob2, ARKODE(), abstol = 1e-8, reltol = 1e-8)

@test norm(sol .- sol2) ≈ 0 atol = 1e-7


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
include(joinpath("..", "common.jl"))

name = "FFplay"

# The products that we will ensure are always built
products = [
    ExecutableProduct("ffplay", :ffplay),
]

# Dependencies that must be installed before this package can be built
# TODO: Theora once it's available
dependencies = [
    BuildDependency("nv_codec_headers_jll"),
    BuildDependency("Xorg_xorgproto_jll"),
    BuildDependency("LibVPX_jll"), # We use the static archive
    Dependency(PackageSpec(name="FFMPEG_jll", version=version)),
    Dependency("SDL2_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script(; ffplay=true), platforms, products, dependencies; preferred_gcc_version=preferred_gcc_version)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using TuringModels, StatsFuns

Turing.setadbackend(:reversediff);

d = CSV.read(joinpath(@__DIR__, "..", "..", "data", "UCBadmit.csv"), delim=';');
size(d) # Should be 12x5

@model m_pois(admit, reject) = begin
   α₁ ~ Normal(0,100)
   α₂ ~ Normal(0,100)

   for i ∈ 1:length(admit)
       λₐ = exp(α₁)
       λᵣ = exp(α₂)
       admit[i] ~ Poisson(λₐ)
       reject[i] ~ Poisson(λᵣ)
   end
end;

chns = sample(m_pois(d[:, :admit], d[:, :reject]), Turing.NUTS(0.65), 1000);

# Rethinking/CmdStan result

m_10_yyt_result = "
    mean   sd 5.5% 94.5% n_eff Rhat
 a1 4.99 0.02 4.95  5.02  2201    1
 a2 5.44 0.02 5.41  5.47  2468    1
";

# Describe the draws

describe(chns)

# End of 10/m10.yyt.jl


text:
export bissect
function bissect_Cub(h :: AbstractNLPModel,
                     nlpstop :: AbstractStopping;
                     t1 = h.meta.x0, t2 = 100.0,
                     verbose = false)

  (tₐ, tᵦ) = trouve_intervalle(h, t1, 2.5)

  γ = 0.8; t = tᵦ; tₚ = tₐ; tqnp = tₐ
  hₖ = [0.0]; hkm1 = [0.0]; gkm1 = [0.0]; hₚ = [0.0]
  iter = 0

  hₖ   = obj(h, t);    gₖ   = grad(h, t)
  hkm1 = obj(h, tqnp); gkm1 = grad(h, tqnp)

  OK = update_and_start!(nlpstop, x = t, gx = gₖ, g0 = copy(gₖ))

  verbose &&
    @printf(" iter        tₚ        t         dN         ")
  verbose &&
    @printf("gₖ          gplus        \n")
  verbose &&
    @printf(" %7.2e %7.2e %7.2e  %7.2e  %7.2e  %7.2e\n",
             iter, tₚ[1], t[1], 0.0, gₖ[1],  0.0)

  while !OK
    s = t .- tqnp; y = gₖ .- gkm1; α = -s
    z = gₖ .+ gkm1 .+ 3 .* (hₖ .- hkm1) ./ α
    discr = z .^ 2 .- gₖ .* gkm1; denom = gₖ .+ gkm1 + 2 .* z
    if true in (discr .> 0.0) && true in (abs.(denom) .> eps(Float64))
      #si on peut on utilise l'interpolation cubique
      w = sqrt.(discr)
      dN = -s .* (gₖ .+ z .+ sign.(α) .* w) ./ (denom)
    else #on se rabat sur une étape de sécante
      dN = -gₖ .* s ./ y
    end

    if true in ((tₚ .- t) .* dN .> 0.0) && true in (dN ./ (tₚ .- t) .< γ)
      tplus = t .+ dN
      hplus = obj(h, tplus)
      gplus = grad(h, tplus)
      # verbose && println("N")
    else
      tplus = (t + tₚ) / 2
      hplus = obj(h, tplus)
      gplus = grad(h, tplus)
      # verbose && println("B")
    end

    if true in (t .> tₚ)
      if true in (gplus .< 0.0)
        tₚ = t; tqnp = t; t = tplus
      else
        tqnp = t; t = tplus
      end
    else
      if true in (gplus .> 0.0)
        tₚ = t; tqnp = t; t = tplus
      else
        tqnp = t; t = tplus
      end
    end

    # update values
    hkm1 = hₖ; gkm1 = gₖ; hₖ = hplus; gₖ = gplus
    OK = update_and_stop!(nlpstop, x = t, fx = hₖ, gx = gₖ)
    iter = iter + 1
    verbose &&
        @printf(" %7.2e %7.2e %7.2e  %7.2e  %7.2e  %7.2e\n",
                iter, tₚ[1], t[1], dN[1], gₖ[1], gplus[1])
  end

  optimal = nlpstop.meta.optimal
  return optimal, nlpstop
end # function


text:
# DWARF support
#=
function read(io::IO,file::ELFFile,h::ELFSectionHeader,::Type{DWARF.ARTable})
    seek(io,h.sh_offset)
    ret = DWARF.ARTable(Array(DWARF.ARTableSet,0))
    while position(io) < h.sh_offset + h.sh_size
        push!(ret.sets,read(io,DWARF.ARTableSet,f.endianness))
    end
    ret
end

function read(io::IO,file::ELFFile,h::ELFSectionHeader,::Type{DWARF.PUBTable})
    seek(io,h.sh_offset)
    ret = DWARF.PUBTable(Array(DWARF.PUBTableSet,0))
    while position(io) < h.sh_offset + h.sh_size
        push!(ret.sets,read(io,DWARF.PUBTableSet,f.endianness))
    end
    ret
end

function read(io::IO,f::ELFFile,h::ELFSectionHeader,::Type{DWARF.AbbrevTableSet})
    seek(io,h.sh_offset)
    read(io,AbbrevTableSet,f.endianness)
end

function read(io::IO,f::ELFFile,h::ELFSectionHeader,s::DWARF.PUBTableSet,::Type{DWARF.DWARFCUHeader})
    seek(io,h.sh_offset+s.header.debug_info_offset)
    read(io,DWARF.DWARFCUHeader,f.endianness)
end

function read(io::IO,f::ELFFile,debug_info::ELFSectionHeader,debug_abbrev::ELFSectionHeader,
    s::DWARF.PUBTableSet,e::DWARF.PUBTableEntry,header::DWARF.DWARFCUHeader,::Type{DWARF.DIE})
    ats = read(io,f,debug_abbrev,header,DWARF.AbbrevTableSet)
    seek(io,debug_info.sh_offset+s.header.debug_info_offset+e.offset)
    read(io,header,ats,DWARF.DIE)
end

function read(io::IO,f::ELFFile,h::ELFSectionHeader,s::DWARF.DWARFCUHeader,::Type{DWARF.AbbrevTableSet})
    seek(io,h.sh_offset+s.debug_abbrev_offset)
    read(io,AbbrevTableSet,f.endianness)
end

function debugsections(io::IO,f::ELFFile)
    snames = names(io,f,f.sheaders)
    sections = Dict{String,ELFSectionHeader}()
    for i in 1:length(snames)
        # Remove leading "."
        ind = findfirst(DEBUG_SECTIONS,snames[i][2:end])
        if ind != 0
            sections[DEBUG_SECTIONS[ind]] = f.sheaders[ind]
        end
    end
    sections
end

function read(io::IO,f::ELFFile,debug_info::ELFSectionHeader,debug_abbrev::ELFSectionHeader,
    s::DWARF.PUBTableSet,e::DWARF.PUBTableEntry,header::DWARF.DWARFCUHeader,::Type{DWARF.DIETree})
    ats = read(io,f,debug_abbrev,header,DWARF.AbbrevTableSet)
    seek(io,debug_info.sh_offset+s.header.debug_info_offset+e.offset)
    ret = DIETree(Array(DWARF.DIETreeNode,0))
    read(io,header,ats,ret,DWARF.DIETreeNode,f.endianness)
    ret
end
=#

struct dl_phdr_info
    dlpi_addr::UInt64
    dlpi_name::Ptr{UInt8}
    dlpi_phdr::Ptr{Cvoid}
    dlpi_phnum::UInt16
end

function callback(info::Ptr{dl_phdr_info},size::Csize_t, data::Ptr{Cvoid})
    push!(unsafe_pointer_to_objref(data),unsafe_load(info))
    convert(Cint,0)
end

function loaded_libraries()
    x = Array(dl_phdr_info,0)
    ccall(:dl_iterate_phdr, Cint, (Ptr{Cvoid}, Any), cfunction(callback, Cint, (Ptr{dl_phdr_info},Csize_t,Ptr{Cvoid})), x)
    x
end

## DWARF Support
function debugsections(h::ELFHandle)
    sects = collect(Sections(h))
    strt = strtab(h)
    snames = map(s->sectionname(s.header;strtab=strt),sects)
    sections = Dict{String,SectionRef}()
    for i in 1:length(snames)
        # remove leading "."
        ind = findfirst(ObjFileBase.DEBUG_SECTIONS,string(snames[i])[2:end])
        if ind != 0
            sections[ObjFileBase.DEBUG_SECTIONS[ind]] = sects[i]
        end
    end
    ObjFileBase.DebugSections(h,sections)
end


text:
# # Polarizability by linear response
#
# We compute the polarizability of a Helium atom. The polarizability
# is defined as the change in dipole moment
# ```math
# \mu = \int r ρ(r) dr
# ```
# with respect to a small uniform electric field ``E = -x``.
#
# We compute this in two ways: first by finite differences (applying a
# finite electric field), then by linear response. Note that DFTK is
# not really adapted to isolated atoms because it uses periodic
# boundary conditions. Nevertheless we can simply embed the Helium
# atom in a large enough box (although this is computationally wasteful).
#
# As in other tests, this is not fully converged, convergence
# parameters were simply selected for fast execution on CI,

using DFTK
using LinearAlgebra

a = 10.
lattice = a * I(3)  # cube of ``a`` bohrs
He = ElementPsp(:He, psp=load_psp("hgh/lda/He-q2"))
atoms = [He => [[1/2; 1/2; 1/2]]]  # Helium at the center of the box

kgrid = [1, 1, 1]  # no kpoint sampling for an isolated system
Ecut = 30
tol = 1e-8

## dipole moment of a given density (assuming the current geometry)
function dipole(ρ)
    basis = ρ.basis
    rr = [a * (r[1] - 1/2) for r in r_vectors(basis)]
    d = sum(rr .* ρ.real) * basis.model.unit_cell_volume / prod(basis.fft_size)
end;

# ## Polarizability by finite differences
# We first compute the polarizability by finite differences.
# First compute the dipole moment at rest:
model = model_LDA(lattice, atoms; symmetries=false)
basis = PlaneWaveBasis(model, Ecut; kgrid=kgrid)
res = self_consistent_field(basis, tol=tol)
μref = dipole(res.ρ)

# Then in a small uniform field:
ε = .01
model_ε = model_LDA(lattice, atoms; extra_terms=[ExternalFromReal(r -> -ε * (r[1] - a/2))],
                    symmetries=false)
basis_ε = PlaneWaveBasis(model_ε, Ecut; kgrid=kgrid)
res_ε = self_consistent_field(basis_ε, tol=tol)
με = dipole(res_ε.ρ)

#-
polarizability = (με - μref) / ε

println("Reference dipole:  $μref")
println("Displaced dipole:  $με")
println("Polarizability :   $polarizability")

# The result on more converged grids is very close to published results.
# For example [DOI 10.1039/C8CP03569E](https://doi.org/10.1039/C8CP03569E)
# quotes **1.65** with LSDA and **1.38** with CCSD(T).

# ## Polarizability by linear response
# Now we use linear response to compute this analytically; we refer to standard
# textbooks for the formalism. In the following, ``\chi_0`` is the
# independent-particle polarizability, and ``K`` the
# Hartree-exchange-correlation kernel. We denote with ``dV_{\rm ext}`` an external
# perturbing potential (like in this case the uniform electric field). Then:
# ```math
# d\rho = \chi_0 dV = \chi_0 (dV_{\rm ext} + K d\rho),
# ```
# which implies
# ```math
# d\rho = (1-\chi_0 K)^-1 \chi_0 dV_{\rm ext}.
# ```
# From this we identify the polarizability operator to be ``\chi = (1-\chi_0 K)^-1 \chi_0``.
# Numerically, we apply ``\chi`` to ``dV = -x`` by solving a linear equation
# (the Dyson equation) iteratively.

using KrylovKit

## KrylovKit cannot deal with the density as a 3D array, so we need `vec` to vectorize
## and `devec` to "unvectorize"
devec(arr) = from_real(basis, reshape(arr, size(res.ρ.real)))

## Apply (1- χ0 K) to a vectorized dρ
function dielectric_operator(dρ)
    dρ = devec(dρ)
    dv = apply_kernel(basis, dρ; ρ=res.ρ)
    χ0dv = apply_χ0(res.ham, res.ψ, res.εF, res.eigenvalues, dv...)[1]
    vec((dρ - χ0dv).real)
end

## dVext is the potential from a uniform field interacting with the dielectric dipole
## of the density.
dVext = from_real(basis, [-a * (r[1] - 1/2) for r in r_vectors(basis)])

## Apply χ0 once to get non-interacting dipole
dρ_nointeract = apply_χ0(res.ham, res.ψ, res.εF, res.eigenvalues, dVext)[1]

## Solve Dyson equation to get interacting dipole
dρ = devec(linsolve(dielectric_operator, vec(dρ_nointeract.real), verbosity=3)[1])

println("Non-interacting polarizability: $(dipole(dρ_nointeract))")
println("Interacting polarizability:     $(dipole(dρ))")

# As expected, the interacting polarizability matches the finite difference
# result. The non-interacting polarizability is higher.


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests related to the LaTeX backend.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Default" begin
    include("latex_backend/default.jl")
end

@testset "Formatters" begin
    include("latex_backend/formatters.jl")
end

@testset "Highlighters" begin
    include("latex_backend/highlighters.jl")
end

@testset "Issues" begin
    include("latex_backend/issues.jl")
end

@testset "Row names" begin
    include("latex_backend/row_names.jl")
end


text:
using Test
@time @testset "Meshing" begin include("test_meshing.jl") end
@time @testset "Basics" begin include("test_basics.jl") end
@time @testset "Miscellaneous" begin include("test_miscellaneous.jl") end
@time @testset "Matrix multiplication" begin include("test_matrix_multiplication.jl") end
true


text:
module Kwant

export Builder,
TranslationalSymmetry,
smatrix,
lattice,
builder

abstract type AbstractKwantObject end

include("load_kwant.jl")
include("Lattice.jl")
include("Builder.jl")

using PyPlot
using .builder

function Base.getproperty(k::AbstractKwantObject, name::Symbol)
    if name==:o
        p = getfield(k,:o)
    else
        p = getproperty(k.o,name)
    end
    if name==:finalized
        q = _FiniteSystem(p)
    elseif name==:neighbors
        q = _Neighbors(p)
    elseif name==:submatrix
        q(args...;kwargs...) = submatrix(pycall(p,PyObject,args...,kwargs...))
    else
        q=p
    end
    return q
end

(p::PyObject)(k::AbstractKwantObject,args...;kwargs...) = pycall(p,PyObject,k.o,args...;kwargs...)

struct _FiniteSystem o::PyObject end
struct _Neighbors o::PyObject end
(fs::_FiniteSystem)() = pycall(fs.o,PyObject)
(nn::_Neighbors)() = pycall(nn.o,PyObject)

PyPlot.plot(b::AbstractKwantObject; kwargs...) = kwant.plot(b.o;kwargs...)

TranslationalSymmetry(x) = kwant.TranslationalSymmetry(x)
wave_function(args...;kwargs...) = kwant.wave_function(args...;kwargs...)

struct smatrix <: AbstractKwantObject
    o::PyObject
    smatrix(args...; kwargs...) = new(pycall(kwant.smatrix,PyObject,args...,kwargs...))
end
struct submatrix{T,N} <: AbstractArray{T,N}
    o::PyObject
    sbm::Array{T,N}
    function submatrix(o::PyObject)
        sbm = convert(Array,o)
        return new{eltype(sbm),ndims(sbm)}(o,sbm)
    end
end
Base.show(io::IO,sbm::submatrix) = show(io,sbm.sbm)
function Base.getproperty(k::submatrix, name::Symbol)
    if name==:o
        p = getfield(k,:o)
    elseif name==:sbm
        p = getfield(k,:sbm)
    else
        p = getproperty(k.o,name)
    end
    return p
end

funs = [:+,:-,:conj,:size,:length,:getindex]
for fun ∈ funs
    @eval Base.$fun(sbm::submatrix,args...;kwargs...) = $fun(sbm.sbm,args...;kwargs...)
    @eval Base.$fun(sbm1::submatrix,sbm2::submatrix,args...;kwargs...) = $fun(sbm1.sbm,sbm2.sbm,args...;kwargs...)
    @eval Base.$fun(x::Any,sbm::submatrix,args...;kwargs...) = $fun(x,sbm.sbm,args...;kwargs...)
end

end # module


text:
# Prevent scalars
CUDA.allowscalar(false)

"""
 Gradient via trce estimation
"""
function grad_ev(X::AbstractArray{Float32, 4}, Y::AbstractArray{Float32, 4},
                 n::Integer, nw::Integer, stride::NTuple=(1, 1))

    nx, ny, nchi, batchsize = size(X)
    nxs, nys, ncho, batchsize = size(Y)
    dW = similar(X, nw*nw, nchi, ncho)
    fill!(dW, 0f0)

    # Get diagonals offsets
    offsets = vcat([((-div(nw,2):div(nw,2)) .+ i*nx) for i=-div(nw,2):div(nw,2)]...)

    # LR product temporaries
    Re =  similar(X, batchsize, n)
    LRe = similar(X, nx*ny*nchi, n)
    LRees = similar(X, nchi, ncho, n)
    e = similar(X, nx*ny*ncho, n)
    disprand!(e)

    # reshape X and Y
    Xloc = reshape(X, :, batchsize)
    Yloc = reshape(Y, :, batchsize)

    # Probing
    LR_probe_batched!(Xloc, Yloc, dW, Re, LRe, e, offsets, n, (nx,ny), LRees, stride)

    scal!(dW, n)
    return reshape(dW, nw, nw, nchi, ncho)[end:-1:1, end:-1:1, :, :]
end

function grad_ev(seeding::Union{UInt64, Tuple{GPUArrays.RNG,UInt64}}, eX::AbstractArray{Float32, 2},
		         Y::AbstractArray{Float32, 4}, w::AbstractArray{Float32}, stride::NTuple=(1, 1))
    n = _params[:p_size]
    nw, _, nchi, ncho = size(w)
    nx, ny = size(Y)[1:2]
    batchsize, probe_bsize = size(eX)
    dW = similar(w, nw*nw, nchi, ncho)
    fill!(dW, 0f0)

    # Get diagonals offsets
    offsets = vcat([((-div(nw,2):div(nw,2)) .+ i*nx*stride[1]) for i=-div(nw,2):div(nw,2)]...)

    # LR product temporaries
    LRe = similar(Y, nx*ny*ncho*prod(stride), probe_bsize)
    LRees = similar(Y, nchi, ncho, probe_bsize)
    e = similar(Y, nx*ny*nchi*prod(stride), probe_bsize)
    disprand!(e, seeding)

    # reshape X and Y
    Yloc = dilate(Y, nx*stride[1], stride)

    # Probing
    LR_probe_batched!(Yloc, eX, dW, LRe, e, offsets, probe_bsize, (nx, ny), LRees, stride)

    scal!(dW, probe_bsize)
    any(isnan.(dW)) && scal!(disprand!(dW), 1f3)
    return reshape(dW, nw, nw, nchi, ncho)[end:-1:1, end:-1:1, :, :]
end

"""
LR_probe_batched!(L::AbstractArray{Float32, 2}, R::AbstractArray{Float32, 2},
                           dW::AbstractArray{Float32, 3}, Re::AbstractArray{Float32, 2}, LRe::AbstractArray{Float32, 2},
                           es::AbstractArray{Float32, 3}, e::AbstractArray{Float32, 2},
                           offsets::AbstractArray{<:Integer, 1}, batch::Integer, nn::Integer,
                           LRees::AbstractArray{Float32, 3})
"""
function LR_probe_batched!(L::AbstractArray{Float32, 2}, Re::AbstractArray{Float32, 2},
                           dW::AbstractArray{Float32, 3}, LRe::AbstractArray{Float32, 2},
                           e::AbstractArray{Float32, 2}, offsets::AbstractArray{<:Integer, 1},
                           batch::Integer, N::NTuple, LRees::AbstractArray{Float32, 3},
                           stride::NTuple=(1, 1))
    se = maximum(offsets)
    nn = prod(N)
    eend = nn-se
    # L*R'*e
    dispgemm!('N', 'N', 1f0, L, Re, 0f0, LRe)
    # e'*L*R'*e
    Lree = view(reshape(LRe, nn, :, batch), se+1:eend, :, :)
    ee = reshape(e, nn, :, batch)
    @inbounds for i=1:length(offsets)
        # Reshape single probe as `nn x nci` and do Mat-vec with es dW for all input channels and sum
        ev = view(ee, se+1+offsets[i]:eend+offsets[i], :, :)
        Bgemm!('T', 'N', 1f0, ev, Lree, 0f0, LRees)
        cumsum!(LRees, LRees, dims=3)
        @views broadcast!(+, dW[i, :, :], dW[i, :, :], LRees[:, :, batch])
    end
end

function LR_probe_batched!(L::AbstractArray{Float32, 2}, R::AbstractArray{Float32, 2},
                           dW::AbstractArray{Float32, 3}, Re::AbstractArray{Float32, 2}, LRe::AbstractArray{Float32, 2},
                           e::AbstractArray{Float32, 2}, offsets::AbstractArray{<:Integer, 1},
                           batch::Integer, N::NTuple, LRees::AbstractArray{Float32, 3},
                           stride::NTuple=(1, 1))
    se = maximum(offsets)
    nn = prod(N)
    eend = nn-se
    # Apply stride
    e, e_s = strided_e(e, N..., stride)
    # R'*e
    dispgemm!('T', 'N', 1f0, R, e_s, 0f0, Re)
    # L*R'*e
    dispgemm!('N', 'N', 1f0, L, Re, 0f0, LRe)
    # e'*L*R'*e
    Lree = view(reshape(LRe, nn, :, batch), se+1:eend, :, :)
    ee = reshape(e, nn, :, batch)
    @inbounds for i=1:length(offsets)
        # Reshape single probe as `nn x nci` and do Mat-vec with es dW for all input channels and sum
        ev = view(ee, se+1-offsets[i]:eend-offsets[i], :, :)
        Bgemm!('T', 'N', 1f0, Lree, ev, 0f0, LRees)
        cumsum!(LRees, LRees, dims=3)
        @views broadcast!(+, dW[i, :, :], dW[i, :, :], LRees[:, :, batch])
    end
end

# rand
disprand!(e::Array{Float32}, seed::UInt64) = (Random.seed!(seed);disprand!(e))
disprand!(e::AbstractArray{Float32}) = (rand!(e);broadcast!(-, e, e, .5f0))
disprand!(e::CuArray{Float32}, seeding::Tuple{GPUArrays.RNG, UInt64}) = (Random.seed!(seeding...);randn!(seeding[1], e);)

# Utilities
function probe_X(X::AbstractArray{Float32, 4})
    seeding = make_rng(X)
    e = probe_vec(seeding, X)
    eX = similar(X, size(X, 4), _params[:p_size])
    dispgemm!('T', 'N', 1f0, reshape(X, :,  size(X, 4)), e, 0f0, eX)
    return seeding, eX
end

function probe_X(X::AbstractArray{Float32, 4}, e::AbstractArray{Float32})
    eX = similar(X, size(X, 4), _params[:p_size])
    dispgemm!('T', 'N', 1f0, reshape(X, :,  size(X, 4)), e, 0f0, eX)
    return eX
end

function probe_vec(seeding, X::AbstractArray{Float32, 4})
    e = similar(X, prod(size(X)[1:3]), _params[:p_size])
    disprand!(e, seeding)
    return e
end

@inline make_rng(::CuArray) = (GPUArrays.default_rng(CuArray), rand(UInt64))
@inline make_rng(::Array) = rand(UInt64)

@inline function strided_e(e::AbstractArray{Float32, 2}, nx::Integer, ny::Integer, stride::NTuple=(1, 1))
    stride == (1, 1) && (return e, e)
    inds = vec(sum.(product(0:stride[1]:nx-1, 1:nx*stride[2]:size(e, 1))))
    e_s = e[inds, :]
    return e, e_s
end

@inline function dilate(Y::AbstractArray{Float32, 4}, nx::Integer, stride::NTuple=(1, 1))
    Y = reshape(Y, :, size(Y, 4))
    stride == (1, 1) && (return Y)
    Yloc = similar(Y, size(Y, 1)*prod(stride), size(Y, 2))
    fill!(Yloc, 0f0)
    inds = vec(sum.(product(0:stride[1]:nx-1, 1:nx*stride[2]:size(Yloc, 1))))
    Yloc[inds, :] .= Y[:, :]
    return Yloc
end

max_stars_repo_name: sheevy/XConv

max_issues_repo_name: sheevy/XConv

max_forks_repo_name: sheevy/XConv


text:
using KernelAbstractions: @index, @kernel, Event
using KernelAbstractions.Extras.LoopInfo: @unroll
using Oceananigans.Utils
using Oceananigans.AbstractOperations: Δz
using Oceananigans.BoundaryConditions
using Oceananigans.Operators

# Evolution Kernels
#=
∂t(η) = -∇⋅U
∂t(U) = - gH∇η + f
=#

@kernel function split_explicit_free_surface_substep_kernel_1!(grid, Δτ, η, U, V, Gᵁ, Gⱽ, g, Hᶠᶜ, Hᶜᶠ)
    i, j = @index(Global, NTuple)
    # ∂τ(U) = - ∇η + G
    @inbounds U[i, j, 1] +=  Δτ * (-g * Hᶠᶜ[i, j] * ∂xᶠᶜᶜ(i, j, 1, grid, η) + Gᵁ[i, j, 1])
    @inbounds V[i, j, 1] +=  Δτ * (-g * Hᶜᶠ[i, j] * ∂yᶜᶠᶜ(i, j, 1, grid, η) + Gⱽ[i, j, 1])
end

@kernel function split_explicit_free_surface_substep_kernel_2!(grid, Δτ, η, U, V, η̅, U̅, V̅, velocity_weight, free_surface_weight)
    i, j = @index(Global, NTuple)
    # ∂τ(η) = - ∇⋅U
    @inbounds η[i, j, 1] -=  Δτ * div_xyᶜᶜᶜ(i, j, 1, grid, U, V)
    # time-averaging
    @inbounds U̅[i, j, 1] +=  velocity_weight * U[i, j, 1]
    @inbounds V̅[i, j, 1] +=  velocity_weight * V[i, j, 1]
    @inbounds η̅[i, j, 1] +=  free_surface_weight * η[i, j, 1]
end

function split_explicit_free_surface_substep!(η, state, auxiliary, settings, arch, grid, g, Δτ, substep_index)
    # unpack state quantities, parameters and forcing terms
    U, V, η̅, U̅, V̅     = state.U, state.V, state.η̅, state.U̅, state.V̅
    Gᵁ, Gⱽ, Hᶠᶜ, Hᶜᶠ  = auxiliary.Gᵁ, auxiliary.Gⱽ, auxiliary.Hᶠᶜ, auxiliary.Hᶜᶠ

    vel_weight = settings.velocity_weights[substep_index]
    η_weight   = settings.free_surface_weights[substep_index]

    fill_halo_regions!(η, arch)

    event = launch!(arch, grid, :xy, split_explicit_free_surface_substep_kernel_1!,
            grid, Δτ, η, U, V, Gᵁ, Gⱽ, g, Hᶠᶜ, Hᶜᶠ,
            dependencies=Event(device(arch)))

    wait(device(arch), event)

    # U, V has been updated thus need to refill halo
    fill_halo_regions!(U, arch)
    fill_halo_regions!(V, arch)

    event = launch!(arch, grid, :xy, split_explicit_free_surface_substep_kernel_2!,
            grid, Δτ, η, U, V, η̅, U̅, V̅, vel_weight, η_weight,
            dependencies=Event(device(arch)))

    wait(device(arch), event)

end

# Barotropic Model Kernels
# u_Δz = u * Δz

@kernel function barotropic_mode_kernel!(U, V, u, v, grid)
    i, j = @index(Global, NTuple)

    # hand unroll first loop
    @inbounds U[i, j, 1] = Δzᶠᶜᶜ(i, j, 1, grid) * u[i, j, 1]
    @inbounds V[i, j, 1] = Δzᶜᶠᶜ(i, j, 1, grid) * v[i, j, 1]

    @unroll for k in 2:grid.Nz
        @inbounds U[i, j, 1] += Δzᶠᶜᶜ(i, j, k, grid) * u[i, j, k]
        @inbounds V[i, j, 1] += Δzᶜᶠᶜ(i, j, k, grid) * v[i, j, k]
    end
end

# may need to do Val(Nk) since it may not be known at compile
function barotropic_mode!(U, V, grid, u, v)
    sum!(U, u * Δz)
    sum!(V, v * Δz)

    arch = architecture(grid)
    fill_halo_regions!(U, arch)
    fill_halo_regions!(V, arch)
end

function set_average_to_zero!(free_surface_state)
    fill!(free_surface_state.η̅, 0.0)
    fill!(free_surface_state.U̅, 0.0)
    fill!(free_surface_state.V̅, 0.0)
end

@kernel function barotropic_split_explicit_corrector_kernel!(u, v, U̅, V̅, U, V, Hᶠᶜ, Hᶜᶠ)
    i, j, k = @index(Global, NTuple)
    @inbounds begin
        u[i, j, k] = u[i, j, k] + (-U[i, j] + U̅[i, j]) / Hᶠᶜ[i, j]
        v[i, j, k] = v[i, j, k] + (-V[i, j] + V̅[i, j]) / Hᶜᶠ[i, j]
    end
end

# may need to do Val(Nk) since it may not be known at compile. Also figure out where to put H
function barotropic_split_explicit_corrector!(u, v, free_surface, grid)
    sefs = free_surface.state
    U, V, U̅, V̅ = sefs.U, sefs.V, sefs.U̅, sefs.V̅
    Hᶠᶜ, Hᶜᶠ = free_surface.auxiliary.Hᶠᶜ, free_surface.auxiliary.Hᶜᶠ
    arch = architecture(grid)

    # take out "bad" barotropic mode,
    # !!!! reusing U and V for this storage since last timestep doesn't matter
    barotropic_mode!(U, V, grid, u, v)
    # add in "good" barotropic mode

    event = launch!(arch, grid, :xyz, barotropic_split_explicit_corrector_kernel!,
        u, v, U̅, V̅, U, V, Hᶠᶜ, Hᶜᶠ,
        dependencies = Event(device(arch)))

    wait(device(arch), event)
end

@inline calc_ab2_tendencies(Gⁿ, G⁻, χ) = (convert(eltype(Gⁿ), (1.5)) + χ) * Gⁿ - (convert(eltype(Gⁿ), (0.5)) + χ) * G⁻

"""
Explicitly step forward η in substeps.
"""
ab2_step_free_surface!(free_surface::SplitExplicitFreeSurface, model, Δt, χ, velocities_update) =
    split_explicit_free_surface_step!(free_surface, model, Δt, χ, velocities_update)

function split_explicit_free_surface_step!(free_surface::SplitExplicitFreeSurface, model, Δt, χ, velocities_update)

    grid = model.grid
    arch = architecture(grid)

    # we start the time integration of η from the average ηⁿ
    η = free_surface.η
    state = free_surface.state
    auxiliary = free_surface.auxiliary
    settings = free_surface.settings
    g = free_surface.gravitational_acceleration

    U, V = (state.U, state.V)
    Δτ = 2 * Δt / settings.substeps  # we evolve for two times the Δt

    u, v, _ = model.velocities # this is u⋆

    Gu = calc_ab2_tendencies(model.timestepper.Gⁿ.u, model.timestepper.G⁻.u, χ)
    Gv = calc_ab2_tendencies(model.timestepper.Gⁿ.v, model.timestepper.G⁻.v, χ)

    # reset free surface averages
    set_average_to_zero!(state)

    # Wait for predictor velocity update step to complete and mask it if immersed boundary.
    wait(device(arch), velocities_update)
    masking_events = Tuple(mask_immersed_field!(q) for q in model.velocities)
    wait(device(arch), MultiEvent(masking_events))

    # Compute barotropic mode of tendency fields
    barotropic_mode!(auxiliary.Gᵁ, auxiliary.Gⱽ, grid, Gu, Gv)

    # Solve for the free surface at tⁿ⁺¹
    start_time = time_ns()

    for substep in 1:settings.substeps
        split_explicit_free_surface_substep!(η, state, auxiliary, settings, arch, grid, g, Δτ, substep)
    end

    # Reset eta for the next timestep
    # this is the only way in which η̅ is used: as a smoother for the
    # substepped η field
    set!(η, free_surface.state.η̅)

    @debug "Split explicit step solve took $(prettytime((time_ns() - start_time) * 1e-9))."

    fill_halo_regions!(η, arch)

    return NoneEvent()
end


text:
module Parametron

# types
export
    Variable,
    LinearTerm,
    QuadraticTerm,
    AffineFunction,
    QuadraticFunction,
    Model,
    Parameter

# enum values
export
    Minimize,
    Maximize

# functions
export
    setobjective!,
    add_nonnegative_constraint!,
    add_nonpositive_constraint!,
    add_zero_constraint!,
    solve!,
    value,
    objectivevalue,
    terminationstatus,
    primalstatus,
    dualstatus,
    findallocs

# macros
export
    @expression,
    @constraint,
    @objective

using LinearAlgebra
using DocStringExtensions

include("FunctionWrappersQuickFix.jl")
using .FunctionWrappersQuickFix: FunctionWrapper

import MathOptInterface
import MacroTools: @capture, postwalk

const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities

@enum Sense Minimize Maximize

include("util.jl")
include("functions.jl")
include("parameter.jl")

using .Functions

include("lazyexpression.jl")
include("moi_interop.jl")
include("model.jl")
include("mockmodel.jl")
include("debug.jl")

end # module


text:
function solve(Jac0, Qleft, Qx, Qy)
	# step756: whenever outputs and/or inputs are densities,
	# pre/postmultiply the jacobians from step 2 by an appropriate
	# matrix so things integrate to 1
	Qright = cat([1,2],Qx',Qy',Qx',Qy')
	Jac1 = Qleft*Jac0*Qright
	#println("finished step 3")

	# apply generalized Schur decomposition a la Klein (2000)
	n = size(Jac1,1)
	A =  Jac1[:,1:n]
	B = -Jac1[:,n+1:2*n]
	QZ        = schur(A,B)
	eigs      = QZ.β ./ QZ.α
	eigselect = abs.(eigs) .< 1 # returns false for NaN gen. eigenvalues
	                           # which is correct here bc they are > 1
	QZ = ordschur(QZ,eigselect)
	# note: because we did (A,B), beta gives eigs corresponding to B

	NK = size(Qx,1) # number of predetermined vars. NOTE: this is after
					# removing elements corresponding to densities
	nk = sum(convert(Array{Int64,1},eigselect)) # number of stable eigs
	println("NK = ",NK)
	println("nk = ",nk)
	if nk>NK
	    warn("equilibrium is locally indeterminate")
	elseif nk<NK
	    warn("no local equilibrium exists")
	end

	U = QZ[:Z]'
	T = QZ[:T]
	S = QZ[:S]
	U11 = U[1:NK,1:NK]
	U12 = U[1:NK,NK+1:end]
	U21 = U[NK+1:end,1:NK]
	U22 = U[NK+1:end,NK+1:end]
	S11 = S[1:NK,1:NK]
	T11 = T[1:NK,1:NK]

	gx_coef = - U22'*pinv(U22*U22')*U21
	# this is the minimum norm solution to U21 + U22*gx = 0
	# in principle one could calculate gx as -pinv(U22)*U21
	# however, the minimum norm solution is more numerically stable
	S11invT11 = S11\T11;
	Ustuff = (U11 + U12*gx_coef);
	invUstuff = Ustuff\eye(NK);
	invterm = pinv(eye(NK)+gx_coef'*gx_coef);
	hx_coef = invterm*Ustuff'*S11invT11*Ustuff;
	# again, in principle there are other ways to calculate hx
	# but this one is the most numerically stable
	# check: hx and S11invT11 should have same eigenvalues
	(eigst,valst) = eig(S11invT11);
	(eighx,valhx) = eig(hx_coef);
	if abs.(maximum(abs.(eighx))-maximum(abs.(eigst)))>1e-4
		warn("max abs eigenvalue of S11invT11 and hx are different!")
	end

	# next, want to represent policy functions in terms of meaningful things

	gx_fval = Qy'*gx_coef*Qx
	hx_fval = Qx'*hx_coef*Qx

	return gx_fval, hx_fval, gx_coef, hx_coef
end


text:
function get_element_type_dict(::AbaqusMeshFormat)

    quad = (Cell{2,4,4}, ("CPE4", "CPS4", "CPE4R", "CPS4R"))
    quad2 = (Cell{2,8,4}, ("CPS8", "CPS8R", "CPE8", "CPE8R"))
    tria = (Cell{2,3,3}, ("CPE3", "CPS3"))
    tria2 = (Cell{2,6,3}, ("CPE6", "CPS6", "CPE6M", "CPS6M"))
    tetra = (Cell{3,4,4}, ("C3D4",))
    tetra2 = (Cell{3,10,4}, ("C3D10",))
    hexa = (Cell{3,8,6}, ("C3D8","C3D8R"))
    hexa2 = (Cell{3,20,6}, ("C3D20","C3D20R"))

    dict = Dict{String,DataType}()
    for types in (quad, tria, quad2, tria2, tetra, tetra2, hexa, hexa2)
        merge!(dict, Dict(code=>types[1] for code in types[2]))
    end
    return dict
end

# Default creator for Ferrite cell types from Abaqus elements
create_cell(::Type{CellType}, node_numbers, ::AbaqusMeshFormat) where{CellType<:Ferrite.AbstractCell} = CellType(ntuple(j->node_numbers[j], length(node_numbers)))


text:
module Biomechanics

using Reexport, DSP, Interpolations, Statistics, StatsBase, LinearAlgebra, NearestNeighbors,
    Distances, LoopVectorization

# utils.jl
export intervals

#reductions.jl
export circmean, circmeand, circstd, circstdd, mean_std_range, avgextrema, intervalextrema

# transformations.jl
export continuousphase, crpensemble, demean, demean!, detrend, detrend!, centraldiff
export ForwardBackwardPad

# timenormalize.jl
export timenormalize, timestoindices, normtime, limitcycle

# gaitevents.jl
export matchevents, toindices, totimes

# gaitstability.jl
export xcom, margin_of_stability, phase_coordination_index

# helperfuncs.jl
export calcresiduals, optfc

include("utils.jl")
include("reductions.jl")
include("transformations.jl")
include("timenormalize.jl")

include("stepmetrics.jl")
@reexport using .SpatiotemporalMetrics

include("gaitevents.jl")
include("gaitstability.jl")
include("helperfuncs.jl")

end # module


text:
using Test
import DynamicPolynomials: @polyvar
# using Revise
using PolynomialBasis

PB = PolynomialBasis

@polyvar x

function allequal(v1,v2,tol)
    return allequal(v1,v2;tol=tol)
end

function allequal(v1,v2;tol=1e3eps())
    np = length(v1)
    f = length(v2) == np
    return f && all([isapprox(v1[i],v2[i],atol=tol) for i = 1:np])
end

funcs = [2.0+0.0x]
deleteat!(funcs,1)
@test_throws AssertionError PB.LagrangeBasis(funcs,zeros(0))

funcs = [2.0+0.0x]
@test_throws AssertionError PB.LagrangeBasis(funcs,[0.0])

funcs = [1.0+0.0x]
@test_throws AssertionError PB.LagrangeBasis(funcs,[0.0,1.0])

@polyvar y
funcs = [1.0+x+y,-x-y]
@test_throws AssertionError PB.LagrangeBasis(funcs,[-1.0,1.0])

funcs = [x^2,1.0+0.0x,-x^2]
@test_throws AssertionError PB.LagrangeBasis(funcs,[-1.0,0.0,1.0])

@test_throws AssertionError PB.LagrangeBasis(-1)
@test_throws MethodError PB.LagrangeBasis(2.5)
@test_throws AssertionError PB.LagrangeBasis(1,start=2)
@test_throws AssertionError PB.LagrangeBasis(2,stop=-2.0)

basis = PB.LagrangeBasis(0)
@test allequal(basis.points,[0.0])
@test allequal(basis(-1.0),[1.0])
@test allequal(basis(1.0),[1.0])
@test allequal(PB.derivative(basis,-1.0),[0.0])
@test allequal(PB.derivative(basis,1.0),[0.0])

basis = PB.LagrangeBasis(0,start=0.0)
@test allequal(basis.points,[0.5])
@test allequal(basis(2.0),[1.0])
@test allequal(PB.derivative(basis,1.0),[0.0])

basis = PB.LagrangeBasis(1)
@test allequal(basis.points,[-1.0,1.0])
@test allequal(basis(-1.0),[1.0,0.0])
@test allequal(basis(1.0),[0.0,1.0])
@test allequal(PB.derivative(basis,-1.0),[-0.5,0.5])
@test allequal(PB.derivative(basis,1.0),[-0.5,0.5])
@test allequal(PB.gradient(basis,-1.0),[-0.5,0.5])
@test allequal(PB.gradient(basis,1.0),[-0.5,0.5])

basis = PB.LagrangeBasis(1,start=0.0)
@test allequal(basis.points,[0.0,1.0])
@test allequal(basis(0.0),[1.0,0.0])
@test allequal(basis(1.0),[0.0,1.0])
@test allequal(PB.derivative(basis,0.0),[-1.0,1.0])
@test allequal(PB.derivative(basis,1.0),[-1.0,1.0])
@test allequal(PB.gradient(basis,-1.0),[-1.0,1.0])
@test allequal(PB.gradient(basis,1.0),[-1.0,1.0])

basis = PB.LagrangeBasis(2)
@test allequal(basis.points,[-1.0,0.0,1.0])
@test allequal(basis(-1.0),[1.0,0.0,0.0])
@test allequal(basis(0.0),[0.0,1.0,0.0])
@test allequal(basis(1.0),[0.0,0.0,1.0])
@test allequal(PB.derivative(basis,-1.0),[-1.5,+2,-0.5])
@test allequal(PB.derivative(basis,0.0),[-0.5,0.0,0.5])
@test allequal(PB.derivative(basis,1.0),[0.5,-2,1.5])
@test allequal(PB.derivative(basis,-1.0),PB.gradient(basis,-1.0))
@test allequal(PB.derivative(basis,0.0),PB.gradient(basis,0.0))
@test allequal(PB.derivative(basis,1.0),PB.gradient(basis,1.0))

basis = PB.LagrangeBasis(2,stop=0.0)
@test allequal(basis.points,[-1.0,-0.5,0.0])
@test allequal(basis(-1.0),[1.0,0.0,0.0])
@test allequal(basis(-0.5),[0.0,1.0,0.0])
@test allequal(basis(0.0),[0.0,0.0,1.0])
@test allequal(PB.derivative(basis,-1.0),[-3.,4,-1])
@test allequal(PB.derivative(basis,-0.5),[-1.0,0,1])
@test allequal(PB.derivative(basis,0.0),[1.0,-4,3])

@test PB.number_of_basis_functions(basis) == 3
@test typeof(basis) == PB.LagrangeBasis{3,typeof(1.0)}
@test supertype(PB.LagrangeBasis) == PB.AbstractBasis{1}

basis = PB.LagrangeBasis(1)
@test allequal(PB.hessian(basis,-1),[0.,0.])
@test allequal(PB.hessian(basis,+1),[0.,0.])

basis = PB.LagrangeBasis(2)
@test allequal(PB.hessian(basis,-1),[1,-2,1])
@test allequal(PB.hessian(basis,0),[1,-2,1])
@test allequal(PB.hessian(basis,1),[1,-2,1])

basis = PB.LagrangeBasis(3)
coeffs = [x^3 for x in basis.points]
h = PB.hessian(basis,-1)
@test allequal(coeffs'*h,[-6.0])
h = PB.hessian(basis,0)
@test allequal(coeffs'*h,[0.0])
h = PB.hessian(basis,+1)
@test allequal(coeffs'*h,[+6.0])

basis = PB.LagrangeBasis(4)
coeffs = [x^4 for x in basis.points]
@test all([allequal(coeffs'*PB.hessian(basis,p),[12p^2],100eps()) for p in basis.points])


text:
# For the purpose of infinite stream tests
include("../infst.jl")

function testptb()
    f = open("../data/ptb.train.txt")
    t_vocab = create_vocab("../data/ptb_3k_vocab")
    sdict = Dict{Int64, Array{Any, 1}}();
    batchsize = 10
    index_to_word = Array(AbstractString, length(t_vocab))
    for (k, v) in t_vocab; index_to_word[v] = k;end;
    res = nextbatch(f, sdict, t_vocab, batchsize)
    result = Any[]

    while (res != nothing)
        res = nextbatch(f, sdict, t_vocab, batchsize)
        (res != nothing) && (push!(result, res))
    end
    isinteractive() && (return result) || println("$(length(result)) many sequences in ptb")

end
!isinteractive() && testptb()

max_stars_repo_name: kirnap/wordsubs

max_issues_repo_name: kirnap/wordsubs

max_forks_repo_name: kirnap/wordsubs


text:
"""
    _broadcast_to(A, shape)
Apply broadcast to make `size(x.values)` as `shape`.
# Examples
```julia-repl
julia> _broadcast_to([1, 2, 3], (3, 3))
3×3 Matrix{Float64}:
 1.0  1.0  1.0
 2.0  2.0  2.0
 3.0  3.0  3.0

julia> _broadcast_to(1, (2, 2))
2×2 Matrix{Float64}:
 1.0  1.0
 1.0  1.0
```
"""
function _broadcast_to(x::R, shape) where R <: Real
    return fill(x, shape)
end

function _broadcast_to(A::AbstractArray{R}, shape) where R <: Real
    return zeros(R, shape) .+ A
end
mutable struct BroadcastTo <: DiffableFunction
    grad_field :: GradField
    shape
    x_shape
    BroadcastTo(shape, x_shape) = new(GradField(), shape, x_shape)
end

function forward(f::BroadcastTo, x)
    y = _broadcast_to(x, f.shape)
    return y
end

function backward(f::BroadcastTo, gy)
    gx = sum_to(gy, f.x_shape)
    return gx
end

"""
    broadcast_to(x::Variable, shape)

Apply broadcast to make `size(x.values)` as `shape`.

# Examples
```julia-repl
julia> x = Variable([1, 2, 3])
name: nothing
values: [1, 2, 3]
creator: User-Defined(nothing)

julia> JITrench.broadcast_to(x, (3, 2))
name: nothing
values: [1.0 1.0; 2.0 2.0; 3.0 3.0]
creator: JITrench.BroadcastTo
```
"""
function broadcast_to(x::Variable, shape)
    if size(x.values) == shape
        return x
    else
        return BroadcastTo(shape, size(x.values))(x)
    end
end


text:
# This file was generated, do not modify it. # hide
@pipeline SimplePipe(hot = OneHotEncoder(),
                     clf = MultinomialClassifier()) is_probabilistic=true
mach = machine(SimplePipe(), Xtrain, ytrain)
res = evaluate!(mach; resampling=Holdout(fraction_train=0.9),
                measure=cross_entropy)
round(res.measurement[1], sigdigits=3)


text:
using VDPTag
using Plots
using POMDPToolbox
using Reel
using ProgressMeter
using PmapProgressMeter
using MCTS
using Plots

@everywhere begin
    using POMDPs
    using POMDPToolbox
    using MCTS
    using VDPTag
    using PmapProgressMeter
end

n_points = 10
mrew = Array(Float64, n_points)
ns = Array(Float64, n_points)
for j = 1:n_points

    N = 1000
    # s_rewards = SharedArray(Float64, N)
    prog = Progress(N, desc="Simulating...")
    rewards = pmap(prog, 1:N, fill(j,N)) do i, j
    # @showprogress "Simulating..." for i in 1:N
        mdp = VDPTagMDP()
        seed = 27
        n = floor(Int, 10^(j/2))
        solver = DPWSolver(depth=20,
                           exploration_constant=40.0,
                           rng=MersenneTwister(seed),
                           estimate_value=RolloutEstimator(ToNextML(mdp)),
                           n_iterations=n,
                           k_action=8.0,
                           alpha_action=1/20,
                           k_state=4.0,
                           alpha_state=1/20)

        policy = solve(solver, mdp)

        hr = HistoryRecorder(max_steps=100, rng=MersenneTwister(i))
        hist = simulate(hr, mdp, policy)
        discounted_reward(hist)
    end
    @show ns[j] = 10^(j/2)
    @show mrew[j] = mean(rewards)
end

unicodeplots()
plot(ns, mrew)
gui()


text:
using Gadfly

set_default_plot_size(6.6inch, 6.6inch)

# contour_color_none, contour_function, contour_layers, contour_matrix
xs = 1:10.
ys = xs
zs = Float64[x^2*log(y) for x in xs, y in ys]

f(x,y) = x*exp(-(x-(round(Int, x)))^2-y^2)
M = Array{Float64}(undef, 200, 200)
for (i, x) in enumerate(range(-8., stop=8., length=200)), (j, y) in enumerate(range(-2., stop=2., length=200))
    M[i, j] = f(x, y)
end

p1 = plot(x=xs, y=ys, z=zs, Geom.contour, color=[colorant"black"])
p2 = plot((x,y) -> x*exp(-(x-(round(Int, x)))^2-y^2), -8., 8, -2., 2)
p3 = plot(layer((x,y) -> x*exp(-(x-(round(Int, x)))^2-y^2), -8., 8, -2., 2),
     layer((x, y) -> sqrt(hypot(x, y)), -2, 2, -1, 1),
     layer(x=[3], y=[1.5], Geom.point))
p4 = plot(z=M, Geom.contour, xmin=[-8.], xmax=[8], ymin=[-2.], ymax=[2])

gridstack([p1 p2; p3 p4])


text:
using ContinuumArrays, BlockArrays, InfiniteArrays, Test
import ContinuumArrays: PiecewiseBasis, VcatBasis, HvcatBasis, arguments, ApplyLayout, checkpoints, UnionDomain

@testset "ConcatBasis" begin
    @testset "hcat" begin
        x = Inclusion(-1..1)
        H = [ones(x) x]
        @test H[:,1:2] isa SubQuasiArray
        Q  =QuasiFill(2,(x,Base.oneto(∞)))
        H∞ = [ones(x) Q]
        @test H∞[:,1:∞] isa SubQuasiArray
        @test (H∞')[1:∞,:] isa SubQuasiArray
    end
    @testset "PiecewiseBasis" begin
        S1 = LinearSpline(0:1)
        S2 = LinearSpline(2:3)
        S = PiecewiseBasis(S1, S2)

        @test S == S == copy(S)
        @test checkpoints(S) == union(checkpoints(S1), checkpoints(S2))

        @test S[0.5,1:4] == [S1[0.5,1:2]; zeros(2)]
        @test S[2.5,1:4] == [zeros(2); S2[2.5,1:2]]
        @test_throws BoundsError S[1.5,2]
        @test_throws BoundsError S[0.5,5]

        D = Derivative(axes(S,1))
        D1 = Derivative(axes(S1,1))
        D2 = Derivative(axes(S2,1))
        @test (D*S)[0.5,1:4] == [(D1 * S1)[0.5,1:2]; zeros(2)]
        @test (D*S)[2.5,1:4] == [zeros(2); (D2 * S2)[2.5,1:2]]

        @test_throws BoundsError (D*S)[1.5,2]
        @test_throws BoundsError (D*S)[0.5,5]

        @test_throws DimensionMismatch D1*S

        @testset "Vec case" begin
            Sv = PiecewiseBasis([S1,S2])
            @test axes(Sv,2) isa BlockedUnitRange
            @test Sv[0.5,1:4] == S[0.5,1:4]
        end

        @testset "UnionDomain with point checkpoints" begin
            @test 0 ∈ checkpoints(UnionDomain(0, 1..2))
        end
    end

    @testset "VcatBasis" begin
        S1 = LinearSpline(0:1)
        S2 = LinearSpline(0:0.5:1)
        S = VcatBasis(S1, S2)

        @test size(S,2) == 5
        @test axes(S,1) == axes(S1,1) == axes(S2,1)
        @test blockaxes(S) == (Block.(1:1), Block.(1:2))

        @test S == S

        @test S[0.1,1:5] == [vcat.(S1[0.1,:],0); vcat.(0, S2[0.1,:])]
        @test_throws BoundsError S[1.1,1]
        @test_throws BoundsError S[0.1,6]

        @test permutedims(S)[1:5,0.1] == S[0.1,1:5]

        D = Derivative(axes(S,1))
        @test (D*S)[0.1,1:5] == [vcat.((D*S1)[0.1,:],0); vcat.(0, (D*S2)[0.1,:])]

        @test_throws BoundsError (D*S)[1.5,2]
        @test_throws BoundsError (D*S)[0.5,6]

        H = VcatBasis(HeavisideSpline(S1.points), HeavisideSpline(S2.points))
        @test H \ (D*S) == [-1 1 0 0 0; 0 0 -2 2 0; 0 0 0 -2 2]
    end

    @testset "HvcatBasis" begin
        S1 = LinearSpline(0:1)
        S2 = LinearSpline(0:0.5:1)
        S = HvcatBasis(2, S1, S2, S2, S1)
        D = Derivative(axes(S,1))

        @test S == S

        @test S[0.1, 1] == [S1[0.1,1] 0; 0 0]
        @test S[0.1,Block(1)] == [[S1[0.1,1] 0; 0 0], [S1[0.1,2] 0; 0 0]]
        @test S[0.1,Block(1)[1]] == [S1[0.1,1] 0; 0 0]
        @test S[0.1,getindex.(Block(1),1:2)] == [[S1[0.1,1] 0; 0 0], [S1[0.1,2] 0; 0 0]]
        D = Derivative(axes(S,1))
        @test_broken (D*S)[0.1,1] # throws error

        v = view(S, :, Block.(2:3))
        @test v[0.1,1] == S[0.1,3]
        @test blockisequal(axes(arguments(ApplyLayout{typeof(*)}(), v)[2],1), axes(S,2))
    end
end


text:
using Test
using BenchmarkTools

# Include function that generates the initial point
include("./initial-point.jl")

# Include function that computes the minimum distance
include("./mindist.jl")

# Function that computes the minimum distance, with a naive algorithm
function mindist_naive(p,cutoff)
  np = length(p)
  d = +Inf
  for ip in 1:np-1
    for jp in ip+1:np
      d = min(d,sqrt( (p[ip][1]-p[jp][1])^2 +
                      (p[ip][2]-p[jp][2])^2) )
    end
  end
  return d
end

# Compute the cell to which each particle belongs
function celllist(p,side,cutoff)
  n = round(Int64,side/cutoff)
  cellparticles = [ Int64[] for i in 1:n, j in 1:n ]
  for i in 1:length(p)
    icell = trunc(Int64,p[i][1]/cutoff) + 1
    jcell = trunc(Int64,p[i][2]/cutoff) + 1
    push!(cellparticles[icell,jcell],i)
  end
  return cellparticles
end

# Positions of 10,000 particles, with side 100. and cutoff 10.
# initial point guarantees that no atom pair is closer than 0.9Å
side = 100.
cutoff = 10.
p = initial_point(10_000,side,0.9)

# Compute the cell of each particle
cellparticles = celllist(p,side,cutoff)
println(" Time for particle classification: ")
@btime celllist($p,$side,$cutoff)

println(" Using mindist: ")
@btime mindist($p,$cellparticles,$cutoff)

println(" Using mindist_naive: ")
@btime mindist_naive($p,$cutoff)

# Check result
dmin_naive = mindist_naive(p,cutoff)
dmin = mindist(p,cellparticles,cutoff)
@test dmin ≈ dmin_naive


text:
noconic_options = (; solve_subproblems = false, solve_relaxation = false)

insts = OrderedDict()

insts["test"] = [
    ((3, 4, true),)
    ((3, 4, false),)
    ((3, 4, true), noconic_options)
    ((3, 4, false), noconic_options)
    ((5, 8, true),)
    ((5, 8, false),)
]

function matrixdecomposition_insts(use_nat::Bool, n_max::Int)
    return [((15, n, use_nat),) for n in vcat(50, 50:50:n_max)]
end

insts["nat"] = matrixdecomposition_insts(true, 700)
insts["ext"] = matrixdecomposition_insts(false, 400)

return (MatrixDecomposition, insts)


text:
using DataFrames
using KernelDensity

KWArgT = Union{Dict, NamedTuple, Nothing}

update_args(args::Union{Dict, NamedTuple}, update::Nothing) = args
update_args(args::Union{Dict, NamedTuple}, update::Union{Dict, NamedTuple}) =
    merge([Dict{Symbol, Any}(zip(keys(d), values(d))) for d in (args, update)]...)

estimate_density_kde(coords::Matrix{Float64}, points::Matrix{Float64}, bandwidth::T where T <: Real)::Vector{Float64} =
    InterpKDE(kde((coords[1,:], coords[2,:]), bandwidth=(Float64(bandwidth), Float64(bandwidth)))).itp.(points[1,:], points[2,:])

function val_range(arr::AT where AT <: AbstractArray{<:Real})
    if length(arr) == 0
        return (nothing, nothing)
    end

    min_val, max_val = arr[1], arr[1]
    for v in arr
        min_val = fmin(min_val, v)
        max_val = fmax(max_val, v)
    end

    return min_val, max_val
end

count_array(values::VT where VT<: AbstractVector{<:Integer}, args...; max_value::Union{<:Integer, Nothing}=nothing, kwargs...) =
    count_array!(zeros(Int, max_value !== nothing ? max_value : maximum(values)), values, args...; erase_counts=false, kwargs...)

function count_array!(counts::VT1 where VT1 <: AbstractVector{<:Integer}, values::VT2 where VT2 <: AbstractVector{<:Integer}; drop_zero::Bool=false, erase_counts::Bool=true)
    if erase_counts
        counts .= 0
    end

    has_zero = false
    for v in values
        if v == 0
            has_zero = true
            continue
        end

        counts[v] += 1
    end

    if !drop_zero && has_zero
        @warn "Array has zero values. It was ignored."
    end

    return counts
end

function count_array!(counts::VT1 where VT1 <: AbstractVector{RT}, values::VT2 where VT2 <: AbstractVector{<:Integer}, weights::VT3 where VT3 <: AbstractVector{RT};
        drop_zero::Bool=false, erase_counts::Bool=true) where RT<:Real
    if erase_counts
        counts .= 0
    end

    has_zero = false
    for i in 1:length(values)
        v = values[i]
        if v == 0
            has_zero = true
            continue
        end

        counts[v] += weights[i]
    end

    if !drop_zero && has_zero
        @warn "Array has zero values. It was ignored."
    end

    return counts
end

function prob_array(values::Union{Array{Int, 1}, SubArray{Int,1}}; max_value::Union{Int, Nothing}=nothing, smooth::Float64=0.0)
    if max_value === nothing
        max_value = maximum(values)
    end

    sum_value = length(values) + max_value * smooth
    counts = fill(smooth / sum_value, max_value)
    for v in values
        counts[v] += 1.0 / sum_value
    end

    return counts
end

function prob_array!(counts::Union{Array{Float64, 1}, SubArray{Float64,1}}, values::Array{Int, 1}; smooth::Float64=0.0)
    sum_value = length(values) + length(counts) * smooth
    counts .= smooth / sum_value
    for v in values
        counts[v] += 1.0 / sum_value
    end

    return counts
end

function split(vector::T where T <: AbstractVector; n_parts::Int)
    offset = ceil(Int, length(vector) / n_parts)
    return [vector[n:min(n + offset - 1, length(vector))] for n in 1:offset:length(vector)]
end

function split(array::T where T <: AbstractVector{TV}, factor::T2 where T2 <: AbstractVector{<:Integer}; max_factor::Union{Int, Nothing}=nothing, drop_zero::Bool=false)::Array{Vector{TV}, 1} where TV
    @assert length(array) == length(factor)
    if max_factor === nothing
        max_factor = maximum(factor)
    end

    splitted = [TV[] for i in 1:max_factor]
    for i in 1:length(array)
        if drop_zero && factor[i] == 0
            continue
        end

        push!(splitted[factor[i]], array[i])
    end

    return splitted
end

split(array::UnitRange{Int64}, factor::Array{Int64,1}; kwargs...) = split(collect(array), factor; kwargs...)

split_ids(factor::Array{Int, 1}; kwargs...) = split(1:length(factor), factor; kwargs...)

split(df::DataFrame, factor::Symbol; kwargs...) = split(df, Array(df[!, factor]); kwargs...)
split(df::DataFrame, factor::Vector{Int}; kwargs...) = [df[ids, :] for ids in split(1:size(df, 1), factor; kwargs...)]

@inline @fastmath function fmax(v1::T, v2::T) where T <: Real
    v1 > v2 ? v1 : v2
end

@inline @fastmath function fmin(v1::T, v2::T) where T <: Real
    v1 < v2 ? v1 : v2
end

function estimate_difference_l0(m1::Matrix{Float64}, m2::Matrix{Float64}; col_weights::Union{Nothing, Vector{Float64}}=nothing, change_threshold::Float64=1e-7)::Tuple{Float64, Float64}
    max_diff = 0.0
    n_changed = 0
    if !all(size(m1) .== size(m2))
        error("Matrices must be of the same size")
    end

    @inbounds for ci in 1:size(m1, 2)
        c_max = 0.0
        for ri in 1:size(m1, 1)
            c_max = fmax(abs(m1[ri, ci] - m2[ri, ci]), c_max)
        end

        if col_weights !== nothing
            c_max *= col_weights[ci]
        end

        if c_max > change_threshold
            n_changed += 1
        end

        max_diff = fmax(c_max, max_diff)
    end

    return max_diff, n_changed / size(m1, 2)
end

### Statistics

@inline function fsample(w::Vector{Float64})::Int
    n = length(w)
    if n == 0
        error("Empty vector for sampling")
    end

    t = rand(Random.GLOBAL_RNG) * sum(w)
    i = 1
    cw = w[1]
    while cw < t && i < length(w)
        i += 1
        @inbounds cw += w[i]
    end
    return i
end

@inline @inbounds fsample(arr::Vector{Int}, w::Vector{Float64})::Int = arr[fsample(w)]

function wmean(values::Vector{Float64}, weights::T where T <: AbstractVector{Float64}; non_zero_ids::Union{UnitRange{Int}, Vector{Int}}=1:length(values))
    s, ws = 0.0, 0.0
    for i in non_zero_ids
        s += values[i] * weights[i]
        ws += weights[i]
    end

    return s / ws
end

function wmean_std(values::Vector{Float64}, weights::T where T <: AbstractVector{Float64}; non_zero_ids::Union{UnitRange{Int}, Vector{Int}}=1:length(values))
    m = wmean(values, weights; non_zero_ids=non_zero_ids)
    s, ws = 0.0, 0.0
    for i in non_zero_ids
        dv = (values[i] - m)
        s += dv * dv * weights[i]
        ws += weights[i]
    end

    return m, sqrt(s / ws)
end

function estimate_hist(vec::Vector{<:Real}, weights=FrequencyWeights(ones(length(vec)));
        ext_cols::NamedTuple=NamedTuple(), rel_width::Float64=0.9, normalize::Union{Bool, Symbol}=false, center=true, bins=nothing, kwargs...)
    hf = (bins === nothing) ? fit(Histogram, vec, weights; kwargs...) : fit(Histogram, vec, weights, bins; kwargs...)
    diffs = rel_width * diff(hf.edges[1])[1]
    df = DataFrame(:s => hf.edges[1][1:end-1], :e => hf.edges[1][1:end-1] .+ diffs, :h => hf.weights)
    if center
        df.s .-= 0.5 * diffs
        df.e .-= 0.5 * diffs
    end
    for (k,v) in pairs(ext_cols)
        df[!,k] .= v
    end

    if isa(normalize, Bool)
        normalize = normalize ? :density : :none
    end

    if normalize == :density
        df[!, :h] = df.h ./ sum(df.h) ./ hf.edges[1].step.hi
    elseif normalize == :frac
        df[!, :h] = df.h ./ sum(df.h)
    elseif normalize != :none
        error("Unknown normalize")
    end

    return df
end

max_stars_repo_name: repushko/Baysor

max_issues_repo_name: repushko/Baysor

max_forks_repo_name: repushko/Baysor


text:
using DataDeps

register(DataDep(
  "glass-identification",
  "http://archive.ics.uci.edu/ml/datasets/Glass+Identification",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/glass/glass.data",
  "dd67373f4baf2807345df02cbfef2093d342e61ad0d82a4fb79af43ef8ce449d",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=11,
        feature_indices=2:10,
        id_indices=1,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


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

"""
    RegressionTask(x, y)

A regression task consists of finding a function `f` such that `y ~ f(x)`
for all training examples `(x,y)` with `y` a continuous variable.
"""
struct RegressionTask{N} <: AbstractLearningTask
  features::NTuple{N,Symbol}
  label::Symbol
end

RegressionTask(x::Symbol, y::Symbol) = RegressionTask{1}((x,), y)

RegressionTask(x::AbstractVector{Symbol}, y::Symbol) =
  RegressionTask{length(x)}(Tuple(x), y)

"""
    features(task)

Return the features of the regression `task`.
"""
features(task::RegressionTask) = task.features

"""
    label(task)

Return the label of the regression `task`.
"""
label(task::RegressionTask) = task.label

issupervised(task::RegressionTask) = true

# ------------
# IO methods
# ------------
function Base.show(io::IO, task::RegressionTask)
  x = task.features
  y = task.label
  lhs = length(x) > 1 ? "("*join(x, ", ")*")" : "$(x[1])"
  print(io, "Regression $lhs → $y")
end


text:
# This file contains some real first-rate type piracy. Sorry in advance.

LinearAlgebra.UpperTriangular(D::Diagonal) = D


text:
facts("FinancialID constructors works") do

    context("Strings work") do
        @fact Ticker("AAPL").id              => "AAPL"
        @fact CUSIP("008000AA7").id         =>  "008000AA7"
        @fact ReutersID("IBM.N").id          => "IBM.N"
        @fact BloombergID("627065740034").id => "627065740034"
    end

    context("non-strings fail") do
        @fact_throws Ticker(123).id
        @fact_throws CUSIP(123).id
        @fact_throws ReutersID(123).id
        @fact_throws BloombergID(123).id
    end
end

facts("Currency constructor works") do

    context("Strings work") do
        @fact Currency("AED").origin => "AED"
    end

    context("non-strings fail") do
        @fact_throws Currency(123).origin
    end
end


text:
module HWconstrained

greet() = print("Hello World!")

	using JuMP, NLopt, DataFrames, Ipopt
	using LinearAlgebra

	export data, table_NLopt, table_JuMP

	function data(a=0.5)

		return Dict("a"=>a,"na"=>na,"nc"=>nc,"ns"=>ns,"nss"=>nss,"e"=>e,"p"=>p,"z"=>z,"pi"=>pi)
	end

	function max_JuMP(a=0.5)

		return Dict("obj"=>objective_value(m),"c"=>value(c),"omegas"=>[value(omega[i]) for i in 1:length(omega)])
	end

	function table_JuMP()
		d = DataFrame(a=[0.5;1.0;5.0],c = zeros(3),omega1=zeros(3),omega2=zeros(3),omega3=zeros(3),fval=zeros(3))
		for i in 1:nrow(d)
			xx = max_JuMP(d[i,:a])
			d[i,:c] = xx["c"]
			d[i,:omega1] = xx["omegas"][1]
			d[i,:omega2] = xx["omegas"][2]
			d[i,:omega3] = xx["omegas"][3]
			d[i,:fval] = xx["obj"]
		end
		return d
	end

	function obj(x::Vector,grad::Vector,data::Dict)

	end

	function constr(x::Vector,grad::Vector,data::Dict)

	end

	function max_NLopt(a=0.5)

	end

	function table_NLopt()
		d = DataFrame(a=[0.5;1.0;5.0],c = zeros(3),omega1=zeros(3),omega2=zeros(3),omega3=zeros(3),fval=zeros(3))
		for i in 1:nrow(d)
			xx = max_NLopt(d[i,:a])
			for j in 2:ncol(d)-1
				d[i,j] = xx[2][j-1]
			end
			d[i,end] = xx[1]
		end
		return d
	end

end # module


text:
@enum GradientDescentStoppingCrit begin
    unknown
    max_iter_reached
    small_gradient
    small_step_size
    small_function_value_delta
end

struct GradientDescentStatistics
    x_sol
    y_sol
    iter
    stopping_crit
end

function Base.show(io::IO, stat::GradientDescentStatistics)
    print("Min search with $(stat.iter)-iterations\n")
    println("Cause of termination: $(stat.stopping_crit)")
    println("f: $(stat.y_sol)")
    print("x: ")
    show(
            stdout,
            MIME("text/plain"),
            stat.x_sol )
    print("\n")

end

struct ToleranceOptions
    step_size_tol::Float64
    f_abs_tol::Float64
    f_rel_tol::Float64
    g_abs_tol::Float64
    max_iter::Int
end

function ToleranceOptions(;
    step_size_tol = 1e-15,
    f_abs_tol = 1e-15,
    f_rel_tol = 1e-15,
    g_abs_tol = 1e-15,
    max_iter = 100)

    ToleranceOptions(
        step_size_tol,
        f_abs_tol,
        f_rel_tol,
        g_abs_tol,
        max_iter)
end

mutable struct OptiSearchSpacePath{T}
    x::Matrix{T}
end
function OptiSearchSpacePath(x0)
    obj = OptiSearchSpacePath{ eltype(x0) }( zeros( eltype(x0), length(x0), 1 ) )
    obj.x[:, 1] = x0
    return obj
end

abstract type DescentMethod
end

function gradient_descent_template(
    m::T,
    init_fcn!,
    step_direction_fcn!,
    update_fcn!,
    fcn,
    x0;
    tol_options::ToleranceOptions       = ToleranceOptions(),
    search_options::LineSearchOptions   = LineSearchOptions(),
    gradient_fcn                        = numder.gradient_fw,
    line_search_fcn                     = line_search,
    log_path                            = false) where {T <: DescentMethod}

    if log_path
        opti_search_path = OptiSearchSpacePath(x0)
    end

    mach_eps = tol_options.f_abs_tol * 1e-5

    init_fcn!(m, x0)
    m.f = fcn(m.x)

    stopping_crit = unknown::GradientDescentStoppingCrit
    x_sol = m.x
    y_sol = m.f

    iter = 0
    for outer iter = 1:tol_options.max_iter

        # m.g = gradient_fcn(fcn, m.x)

        # step direction algorithm
        step_direction_fcn!(m, fcn, gradient_fcn)

        # stopping criterion based on the gradient norm
        norm_g = norm(m.g)
        if norm_g <= tol_options.g_abs_tol
            x_sol = m.x
            y_sol = m.f
            stopping_crit = small_gradient::GradientDescentStoppingCrit
            break
        end

        # line search for step size
        step_size = line_search_fcn(fcn, m.x, m.step_direction, m.f, search_options)

        # stopping criterion based on the step size
        abs_step_size = norm( step_size * m.step_direction )
        if abs_step_size <= tol_options.step_size_tol
            x_sol = m.x
            y_sol = m.f
            stopping_crit = small_step_size::GradientDescentStoppingCrit
            break
        end

        # next step
        x1 = m.x + step_size * m.step_direction
        f1 = fcn(x1)

        # stopping criteria for function value
        if abs(f1 - m.f) <= tol_options.f_abs_tol || abs(f1 - m.f) / ( abs(m.f) == 0 ? mach_eps : abs(m.f)) < tol_options.f_rel_tol
            x_sol = x1
            y_sol = f1
            stopping_crit = small_function_value_delta::GradientDescentStoppingCrit
            break
        end

        # update for the next iteration
        update_fcn!(m, x1, f1)
        if log_path
            # push!(opti_search_path.x, x1)
            opti_search_path.x = hcat(opti_search_path.x, x1)
        end
    end

    # last step: iterations run out
    if iter == tol_options.max_iter
        x_sol = m.x
        y_sol = m.f
        stopping_crit = max_iter_reached::GradientDescentStoppingCrit
    end

    stat = GradientDescentStatistics(x_sol, y_sol, iter, stopping_crit)
    if log_path
        return (stat, opti_search_path)
    else
        return stat
    end

end

mutable struct SteepestGradientDescentMethod <: DescentMethod
    x
    f

    g
    step_direction
end
function SteepestGradientDescentMethod()
    empty_array = zeros(Float64, 0)
    empty_scalar = zero(Float64)
    SteepestGradientDescentMethod( empty_array, empty_scalar, empty_array, empty_array )
end

function init!(m::SteepestGradientDescentMethod, x0)
    m.g = zeros(eltype(x0), size(x0))
    m.x = x0
end

function step_direction_fcn_steepestgrad!(m, fcn, gradient_fcn)
    m.g = gradient_fcn(fcn, m.x)
    m.step_direction = -m.g ./ norm(m.g)
end

function update_fcn_steepestgrad!(m, x1, f1)
    m.f                   = f1
    m.x                   = x1
end

function steepest_gradient_descent(
    fcn,
    x0;
    tol_options::ToleranceOptions       = ToleranceOptions(),
    search_options::LineSearchOptions   = LineSearchOptions(),
    gradient_fcn                        = numder.gradient_fw,
    line_search_fcn                     = line_search,
    log_path                            = false)

    m                   = SteepestGradientDescentMethod();
    init_fcn            = init!
    step_direction_fcn  = step_direction_fcn_steepestgrad!
    update_fcn          = update_fcn_steepestgrad!

    stat = gradient_descent_template(
        m,
        init_fcn,
        step_direction_fcn,
        update_fcn,
        fcn,
        x0;
        tol_options     = tol_options,
        search_options  = search_options,
        gradient_fcn    = gradient_fcn,
        line_search_fcn = line_search_fcn,
        log_path        = log_path)

    return stat

end

mutable struct ConjugateGradientDescentMethod <: DescentMethod
    x                   # current position
    f                   # current function value

    g                   # gradient
    g_prev              # previous gradient
    step_direction      # step direction
    step_direction_prev # previous step direction

    beta                # prev step direction effect
end
function ConjugateGradientDescentMethod()
    empty_array = zeros(Float64, 0)
    empty_scalar = zero(Float64)
    ConjugateGradientDescentMethod( empty_array, empty_scalar, empty_array, empty_array, empty_array, empty_array, empty_scalar )
end

function init!(m::ConjugateGradientDescentMethod, x0)
    m.g_prev              = zeros(eltype(x0), size(x0))
    m.step_direction_prev = zeros(eltype(x0), size(x0))
    m.beta                = zero(eltype(x0))

    m.x = x0
end

function step_direction_fcn_cg_pr!(m, fcn, gradient_fcn)
    # conjugate gradient
    # polak-ribiere method

    # step direction algo
    m.g = gradient_fcn(fcn, m.x)
    g_prev_normsq = sum(m.g_prev .* m.g_prev)
    if g_prev_normsq <= 1e-15
        m.beta = 0
    else
        m.beta = sum(m.g .* (m.g - m.g_prev)) / g_prev_normsq
        m.beta = max(m.beta, 0)
    end
    m.step_direction = -m.g + m.beta * m.step_direction_prev

    return m.step_direction
end

function update_fcn_cg_pr!(m, x1, f1)
    m.f                   = f1
    m.x                   = x1
    m.g_prev              = m.g
    m.step_direction_prev = m.step_direction
end

function cg_pr(
    fcn,
    x0;
    tol_options::ToleranceOptions = ToleranceOptions(),
    search_options::LineSearchOptions = LineSearchOptions(),
    gradient_fcn    = numder.gradient_fw,
    line_search_fcn = line_search,
    log_path        = false)

    m                   = ConjugateGradientDescentMethod();
    init_fcn            = init!
    step_direction_fcn  = step_direction_fcn_cg_pr!
    update_fcn          = update_fcn_cg_pr!

    stat = gradient_descent_template(
        m,
        init_fcn,
        step_direction_fcn,
        update_fcn,
        fcn,
        x0;
        tol_options     = tol_options,
        search_options  = search_options,
        gradient_fcn    = gradient_fcn,
        line_search_fcn = line_search_fcn,
        log_path        = log_path)

    return stat

end

function steepest_gradient_descent_original(
    fcn,
    x0;
    tol_options::ToleranceOptions = ToleranceOptions(),
    search_options::LineSearchOptions = LineSearchOptions())

    # max_iter = 100, step_size_tol = 1e-15, f_abs_tol = 1e-15, f_rel_tol = 1e-15, g_abs_tol = 1e-16)

    max_iter      = tol_options.max_iter
    step_size_tol = tol_options.step_size_tol
    f_abs_tol     = tol_options.f_abs_tol
    f_rel_tol     = tol_options.f_rel_tol
    g_abs_tol     = tol_options.g_abs_tol

    mach_eps = f_abs_tol * 1e-5

    f0 = fcn(x0)
    iter = 0
    for outer iter = 1:max_iter
        g = numder.gradient_fw(fcn, x0)

        norm_g = norm(g)
        if norm_g <= g_abs_tol
            x_sol = x0
            y_sol = f0
            return GradientDescentStatistics(x_sol, y_sol, iter, small_gradient::GradientDescentStoppingCrit)
        end

        # step direction algo
        step_direction = -g ./ norm_g

        step_size = line_search(fcn, x0, step_direction, f0, search_options)
        if abs(step_size) <= step_size_tol
            x_sol = x0
            y_sol = f0
            return GradientDescentStatistics(x_sol, y_sol, iter, small_step_size::GradientDescentStoppingCrit)
        end

        x1 = x0 + step_size * step_direction
        f1 = fcn(x1)
        if abs(f1 - f0) <= f_abs_tol || abs(f1 - f0) / ( abs(f0) == 0 ? mach_eps : abs(f0)) < f_rel_tol
            x_sol = x1
            y_sol = f1
            return GradientDescentStatistics(x_sol, y_sol, iter, small_function_value_delta::GradientDescentStoppingCrit)
        end

        # update for the next iteration
        f0 = f1
        x0 = x1
    end
    if iter == max_iter
        return GradientDescentStatistics(x0, f0, iter, max_iter_reached::GradientDescentStoppingCrit)
    end

end

include("accelerated_grad_descent.jl")

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
"""
    erdos_renyi(n::Int, p::Real, G=Graph; seed=-1)
    erdos_renyi(n::Int, m::Int, G=Graph; seed=-1)

Creates an [Erdős–Rényi](http://en.wikipedia.org/wiki/Erdős–Rényi_model)
random graph of type `G` with `n` vertices.
Edges are added between pairs of vertices with probability `p` in the first method.
In the second method `m` edges are randomly chosen insted.

Undirected graphs are created by default. Directed graphs can be created
passing a directed graph type as last argument (e.g. `DiGraph`)

Note also that Erdős–Rényi graphs may be generated quickly using `erdos_renyi(n, ne)`
or the  `Graph(nv, ne)` constructor, which randomly select `ne` edges among all the potential
edges.
"""
function erdos_renyi(n::Int, p::Real, ::Type{G} = Graph;
        seed::Int=-1) where G<:AGraphOrDiGraph
    m = is_directed(G) ? n*(n-1) : div(n*(n-1),2)
    seed > 0 && seed_dsfmt(seed)
    nedg = randbinomial(m, p) # sadly StatsBase doesn't support non-global RNG
    return erdos_renyi(n, nedg, G; seed=seed)
end

function erdos_renyi(n::Int, m::Int, ::Type{G} = Graph;
        seed::Int = -1) where G<:AGraphOrDiGraph
    maxe = is_directed(G) ? n * (n-1) : div(n * (n-1), 2)
    @assert(m <= maxe, "Maximum number of edges for this generator is $maxe")
    m > 2/3 * maxe && return complement(erdos_renyi(n, maxe-m, G; seed=seed))

    rng = getRNG(seed)
    g = G(n)
    while ne(g) < m
        source = rand(rng, 1:n)
        dest = rand(rng, 1:n)
        source != dest && add_edge!(g, source, dest)
    end
    return g
end

"""
    watts_strogatz(n, k, β, G=Graph; seed=-1)

Creates a [Watts-Strogatz](https://en.wikipedia.org/wiki/Watts_and_Strogatz_model)
small model random graph with `n` vertices, each with degree `k`. Edges are
randomized per the model based on probability `β`.

Undirected graphs are created by default. Directed graphs can be created
passing a directed graph type as last argument (e.g. `DiGraph`).
"""
function watts_strogatz(n::Int, k::Int, β::Real, ::Type{G} = Graph;
        seed::Int = -1) where G<:AGraphOrDiGraph
    @assert k < n/2
    g = G(n)
    rng = getRNG(seed)
    for s in 1:n
        for i in 1:(floor(Int, k/2))
            target = ((s + i - 1) % n) + 1
            if rand(rng) > β && !has_edge(g, s, target) # TODO: optimize this based on return of add_edge!
                add_edge!(g, s, target)
            else
                while true
                    d = target
                    while d == target
                        d = rand(rng, 1:n-1)
                        if s < d
                            d += 1
                        end
                    end
                    if s != d
                        add_edge!(g, s, d)[1] && break
                    end
                end
            end
        end
    end
    return g
end

function _suitable(edges::Set{Edge{T}}, potential_edges::Dict{T, T}) where T
    isempty(potential_edges) && return true
    list = keys(potential_edges)
    for s1 in list, s2 in list
        s1 >= s2 && continue
        (Edge{T}(s1, s2) ∉ edges) && return true
    end
    return false
end

_try_creation_rrg(n::T, k::T, rng::AbstractRNG) where {T<:Integer} = _try_creation_rrg(n, fill(k,n), rng)

function _try_creation_rrg(n::T, k::Vector{T}, rng::AbstractRNG) where T<:Integer
    E = Edge{T}
    edges = Set{E}()
    m = 0
    stubs = zeros(T, sum(k))
    for i=1:n
        for j = 1:k[i]
            m += 1
            stubs[m] = i
        end
    end
    # stubs = vcat([fill(i, k[i]) for i=1:n]...) # slower

    while !isempty(stubs)
        potential_edges =  Dict{T,T}()
        shuffle!(rng, stubs)
        for i in 1:2:length(stubs)
            s1, s2 = stubs[i:i+1]
            if (s1 > s2)
                s1, s2 = s2, s1
            end
            e = E(s1, s2)
            if s1 != s2 && ∉(e, edges)
                push!(edges, e)
            else
                potential_edges[s1] = get(potential_edges, s1, 0) + 1
                potential_edges[s2] = get(potential_edges, s2, 0) + 1
            end
        end

        if !_suitable(edges, potential_edges)
            return Set{E}()
        end

        stubs = Vector{T}()
        for (e, ct) in potential_edges
            append!(stubs, fill(e, ct))
        end
    end
    return edges
end

"""
    barabasi_albert(n, k, G=Graph; seed=-1)
    barabasi_albert(n, n0, k, G=Graph; seed=-1)

Creates a random graph of type `G` with `n` vertices according to [Barabási–Albert model](https://en.wikipedia.org/wiki/Barab%C3%A1si%E2%80%93Albert_model).
It is grown by adding new vertices to an initial graph with `n0` vertices (`n0=k` if not specified).
Each new vertex is attached with `k` edges to `k` different vertices already present in the system by preferential attachment.
The initial graph is empty by default.

Undirected graphs are created by default. Directed graphs can be created
passing a directed graph type as last argument (e.g. `DiGraph`).

See also [`barabasi_albert!`](@ref) for growing a given graph.
"""
barabasi_albert(n::Int, k::Int, ::Type{G}=Graph; keyargs...) where {G<:AGraphOrDiGraph} =
    barabasi_albert(n, k, k, G; keyargs...)

function barabasi_albert(n::Int, n0::Int, k::Int, ::Type{G}=Graph;
        seed::Int = -1) where G<:AGraphOrDiGraph
    g = G(n0)
    barabasi_albert!(g, n, k; seed = seed)
    return g
end

"""
    barabasi_albert!(g, n::Int, k::Int; seed::Int = -1)

Grows the graph `g` according to
[Barabási–Albert](https://en.wikipedia.org/wiki/Barab%C3%A1si%E2%80%93Albert_model)
process into a graph with `n` vertices. At each step a new vertex is attached by
preferential attachment to `k` different vertices already present in the graph.

See also [`barabasi_albert`](@ref).
"""
function barabasi_albert!(g::AGraphOrDiGraph, n::Int, k::Int; seed::Int=-1)
    n0 = nv(g)
    1 <= k <= n0 <= n ||
        throw(ArgumentError("Barabási-Albert model requires 1 <= k <= nv(g) <= n"))
    n0 == n && return g

    # seed random number generator
    seed > 0 && srand(seed)

    add_vertices!(g, n - n0)

    # if initial graph doesn't contain any edges
    # expand it by one vertex and add k edges from this additional node
    if ne(g) == 0
        # expand initial graph
        n0 += 1

        # add edges to k existing nodes
        for target in sample!(collect(1:n0-1), k)
            add_edge!(g, n0, target)
        end
    end

    # vector of weighted nodes (each node is repeated once for each adjacent edge)
    weightedNodes = Vector{Int}(undef, 2*(n-n0)*k + 2*ne(g))

    # initialize vector of weighted nodes
    offset = 0
    for e in edges(g)
        weightedNodes[offset+=1] = src(e)
        weightedNodes[offset+=1] = dst(e)
    end

    # array to record if a node is picked
    picked = fill(false, n)

    # vector of targets
    targets = Vector{Int}(undef, k)

    for source in n0+1:n
        # choose k targets from the existing nodes
        # pick uniformly from weightedNodes (preferential attachement)
        i = 0
        while i < k
            target = weightedNodes[rand(1:offset)]
            if !picked[target]
                targets[i+=1] = target
                picked[target] = true
            end
        end

        # add edges to k targets
        for target in targets
            add_edge!(g, source, target)

            weightedNodes[offset+=1] = source
            weightedNodes[offset+=1] = target
            picked[target] = false
        end
    end

    return g
end

"""
    static_fitness_model(m, fitness, G=Graph; seed=-1)
    static_fitness_model(m, fitness_out, fitness_in, G=DiGraph; seed=-1)

Generates a random graph with `length(fitness)` nodes and `m` edges,
in which the probability of the existence of edge `(i, j)` is proportional
to `fitness[i]*fitness[j]`. Time complexity is O(|V| + |E| log |E|).

In and out fitness have to be supplied for generating directed graphs.

Reference:

* Goh K-I, Kahng B, Kim D: Universal behaviour of load distribution
in scale-free networks. Phys Rev Lett 87(27):278701, 2001.
"""
function static_fitness_model(m::Int, fitness::Vector{T},
        ::Type{G}=Graph; seed::Int=-1) where {T<:Real, G<:AGraph}
    @assert(m >= 0, "invalid number of edges")
    n = length(fitness)
    m == 0 && return G(n)
    nodes = 0
    for f in fitness
        # sanity check for the fitness
        f < zero(T) && error("fitness scores must be non-negative")
        f > zero(T) && (nodes += 1)
    end
    # avoid getting into an infinite loop when too many edges are requested
    max_no_of_edges = div(nodes*(nodes-1), 2)
    @assert(m <= max_no_of_edges, "too many edges requested")
    # calculate the cumulative fitness scores
    cum_fitness = cumsum(fitness)
    g = G(n)
    _create_static_fitness_graph!(g, m, cum_fitness, cum_fitness, seed)
    return g
end

function static_fitness_model(m::Int, fitness_out::Vector{T},
        fitness_in::Vector{S}, ::Type{G}=DiGraph; seed::Int=-1) where {T<:Real,S<:Real, G<:ADiGraph}
    @assert(m >= 0, "invalid number of edges")
    n = length(fitness_out)
    @assert(length(fitness_in) == n, "fitness_in must have the same size as fitness_out")
    m == 0 && return G(n)
    # avoid getting into an infinite loop when too many edges are requested
    outnodes = innodes = nodes = 0
    @inbounds for i=1:n
        # sanity check for the fitness
        (fitness_out[i] < zero(T) || fitness_in[i] < zero(S)) && error("fitness scores must be non-negative")
        fitness_out[i] > zero(T) && (outnodes += 1)
        fitness_in[i] > zero(S) && (innodes += 1)
        (fitness_out[i] > zero(T) && fitness_in[i] > zero(S)) && (nodes += 1)
    end
    max_no_of_edges = outnodes*innodes - nodes
    @assert(m <= max_no_of_edges, "too many edges requested")
    # calculate the cumulative fitness scores
    cum_fitness_out = cumsum(fitness_out)
    cum_fitness_in = cumsum(fitness_in)
    g = G(n)
    _create_static_fitness_graph!(g, m, cum_fitness_out, cum_fitness_in, seed)
    return g
end

function _create_static_fitness_graph!(g::AGraphOrDiGraph, m::Int, cum_fitness_out::Vector{T}, cum_fitness_in::Vector{S}, seed::Int) where {T<:Real,S<:Real}
    rng = getRNG(seed)
    max_out = cum_fitness_out[end]
    max_in = cum_fitness_in[end]
    while m > 0
        source = searchsortedfirst(cum_fitness_out, rand(rng)*max_out)
        target = searchsortedfirst(cum_fitness_in, rand(rng)*max_in)
        # skip if loop edge
        (source == target) && continue
        edge = Edge(source, target)
        # is there already an edge? If so, try again
        add_edge!(g, edge)[1] || continue
        m -= 1
    end
end

"""
    function static_scale_free(n, m, α, G=Graph;
            seed=-1, finite_size_correction=true)

Generates a random graph with `n` vertices, `m` edges and expected power-law
degree distribution with exponent `α`. `finite_size_correction` determines
whether to use the finite size correction proposed by Cho et al.
This generator calls internally the `static_fitness_model function`.
Time complexity is O(|V| + |E| log |E|).

    function static_scale_free(n, m, α_out, α_in, G=DiGraph;
            seed=-1, finite_size_correction=true)

Generates a random digraph

References:

* Goh K-I, Kahng B, Kim D: Universal behaviour of load distribution in scale-free networks. Phys Rev Lett 87(27):278701, 2001.

* Chung F and Lu L: Connected components in a random graph with given degree sequences. Annals of Combinatorics 6, 125-145, 2002.

* Cho YS, Kim JS, Park J, Kahng B, Kim D: Percolation transitions in scale-free networks under the Achlioptas process. Phys Rev Lett 103:135702, 2009.
"""
function static_scale_free(n::Int, m::Int, α::Float64, ::Type{G} = Graph;
        seed::Int=-1, finite_size_correction::Bool=true) where G<:AGraph
    @assert(n >= 0, "Invalid number of nodes")
    @assert(α >= 2, "out-degree exponent must be >= 2")
    fitness = _construct_fitness(n, α, finite_size_correction)
    static_fitness_model(m, fitness, G, seed=seed)
end

function static_scale_free(n::Int, m::Int, α_out::Float64, α_in::Float64, ::Type{G} = DiGraph;
        seed::Int=-1, finite_size_correction::Bool=true) where G<:ADiGraph
    @assert(n >= 0, "Invalid number of nodes")
    @assert(α_out >= 2, "out-degree exponent must be >= 2")
    @assert(α_in >= 2, "in-degree exponent must be >= 2")
    # construct the fitness
    fitness_out = _construct_fitness(n, α_out, finite_size_correction)
    fitness_in = _construct_fitness(n, α_in, finite_size_correction)
    # eliminate correlation
    shuffle!(fitness_in)
    static_fitness_model(m, fitness_out, fitness_in, G, seed=seed)
end

function _construct_fitness(n::Int, α::Float64, finite_size_correction::Bool)
    α = -1/(α-1)
    fitness = zeros(n)
    j = float(n)
    if finite_size_correction && α < -0.5
        # See the Cho et al paper, first page first column + footnote 7
        j += n^(1+1/2α) * (10sqrt(2)*(1+α)) ^ (-1/α) - 1
    end
    j = max(j, n)
    @inbounds for i=1:n
        fitness[i] = j ^ α
        j -= 1
    end
    return fitness
end

"""
    random_regular_graph(n::Int, k::Int; seed=-1)

Creates a random undirected
[regular graph](https://en.wikipedia.org/wiki/Regular_graph) with `n` vertices,
each with degree `k`.

For undirected graphs, allocates an array of `nk` `Int`s, and takes
approximately ``nk^2`` time. For ``k > n/2``, generates a graph of degree
`n-k-1` and returns its complement.
"""
function random_regular_graph(n::Int, k::Int, ::Type{G}=Graph;
        seed::Int=-1) where G<:AGraph
    @assert(iseven(n*k), "n * k must be even")
    @assert(0 <= k < n, "the 0 <= k < n inequality must be satisfied")
    if k == 0
        return G(n)
    end
    if (k > n/2) && iseven(n * (n-k-1))
        return complement(random_regular_graph(n, n-k-1, G, seed=seed))
    end

    rng = getRNG(seed)

    T = vertextype(G)
    edges = _try_creation_rrg(T(n), T(k), rng)
    while isempty(edges)
        edges = _try_creation_rrg(T(n), T(k), rng)
    end
    g = G(n)
    for edge in edges
        add_edge!(g, edge)
    end

    return g
end

"""
    random_configuration_model(n::Int, k::Vector{Int}; seed=-1, check_graphical=false)

Creates a random undirected graph according to the [configuraton model]
(http://tuvalu.santafe.edu/~aaronc/courses/5352/fall2013/csci5352_2013_L11.pdf).
It contains `n` vertices, the vertex `i` having degree `k[i]`.

Defining `c = mean(k)`, it allocates an array of `nc` `Int`s, and takes
approximately ``nc^2`` time.

If `check_graphical=true` makes sure that `k` is a graphical sequence (see `is_graphical`).
"""
function random_configuration_model(n::Int, k::Vector{Int}, ::Type{G}=Graph;
        seed::Int=-1, check_graphical::Bool=false) where G<:AGraph
    @assert(n == length(k), "a degree sequence of length n has to be provided")
    m = sum(k)
    @assert(iseven(m), "sum(k) must be even")
    @assert(all(0 .<= k .< n), "the 0 <= k[i] < n inequality must be satisfied")
    if check_graphical
        is_graphical(k) || error("Degree sequence non graphical")
    end
    rng = getRNG(seed)

    edges = _try_creation_rrg(n, k, rng)
    while m > 0 && isempty(edges)
        edges = _try_creation_rrg(n, k, rng)
    end

    g = G(n)
    for edge in edges
        add_edge!(g, edge)
    end

    return g
end

"""
    random_regular_digraph(n::Int, k::Int; dir::Symbol=:out, seed=-1)

Creates a random directed
[regular graph](https://en.wikipedia.org/wiki/Regular_graph) with `n` vertices,
each with degree `k`. The degree (in or out) can be
specified using `dir=:in` or `dir=:out`. The default is `dir=:out`.

For directed graphs, allocates an ``n \times n`` sparse matrix of boolean as an
adjacency matrix and uses that to generate the directed graph.
"""
function random_regular_digraph(n::Int, k::Int, ::Type{G}=DiGraph;
        dir::Symbol=:out, seed::Int=-1) where G<:ADiGraph
    #TODO remove the function sample from StatsBase for one allowing the use
    # of a local rng
    @assert(0 <= k < n, "the 0 <= k < n inequality must be satisfied")

    if k == 0
        return G(n)
    end
    if (k > n/2) && iseven(n * (n-k-1))
        return complement(random_regular_digraph(n, n-k-1, G, dir=dir, seed=seed))
    end
    rng = getRNG(seed)
    cs = collect(2:n)
    i = 1
    I = Vector{Int}(undef, n*k)
    J = Vector{Int}(undef, n*k)
    V = fill(true, n*k)
    for r in 1:n
        l = (r-1)*k+1 : r*k
        I[l] .= r
        J[l] .= sample!(rng, cs, k, exclude = r)
    end

    if dir == :out
        return G(sparse(I, J, V, n, n))
    else
        return G(sparse(I, J, V, n, n)')
    end
end

"""
    stochastic_block_model(c::Matrix{Float64}, n::Vector{Int}; seed::Int = -1)
    stochastic_block_model(cin::Float64, coff::Float64, n::Vector{Int}; seed::Int = -1)

Returns a Graph generated according to the Stochastic Block Model (SBM).

`c[a,b]` : Mean number of neighbors of a vertex in block `a` belonging to block `b`.
           Only the upper triangular part is considered, since the lower traingular is
           determined by ``c[b,a] = c[a,b] * n[a]/n[b]``.
`n[a]` : Number of vertices in block `a`

The second form samples from a SBM with `c[a,a]=cin`, and `c[a,b]=coff`.

For a dynamic version of the SBM see the `StochasticBlockModel` type and
related functions.
"""
function stochastic_block_model(c::Matrix{T}, n::Vector{Int},
        ::Type{G}=Graph; seed::Int = -1) where {T<:Real, G<:AGraph}
    @assert size(c,1) == length(n)
    @assert size(c,2) == length(n)
    # init dsfmt generator without altering GLOBAL_RNG
    seed > 0 && seed_dsfmt(seed)
    rng =  getRNG(seed)

    N = sum(n)
    K = length(n)
    nedg = zeros(Int,K, K)
    g = G(N)
    cum = [sum(n[1:a]) for a=0:K]
    for a=1:K
        ra = cum[a]+1:cum[a+1]
        for b=a:K
            @assert a==b ? c[a,b] <= n[b]-1 : c[a,b] <= n[b]   "Mean degree cannot be greater than available neighbors in the block."

            m = a==b ? div(n[a]*(n[a]-1),2) : n[a]*n[b]
            p = a==b ? n[a]*c[a,b] / (2m) : n[a]*c[a,b]/m
            nedg = randbinomial(m, p)
            rb = cum[b]+1:cum[b+1]
            i=0
            while i < nedg
                source = rand(rng, ra)
                dest = rand(rng, rb)
                if source != dest
                    if add_edge!(g, source, dest)[1]
                        i += 1
                    end
                end
            end
        end
    end
    return g
end

function stochastic_block_model(cint::T, cext::T, n::Vector{Int},
        ::Type{G}=Graph; seed::Int=-1) where {T<:Real, G<:AGraph}
    K = length(n)
    c = [ifelse(a==b, cint, cext) for a=1:K,b=1:K]
    stochastic_block_model(c, n, G, seed=seed)
end


text:
using Test, ExtendableGrids, GridVisualize
import PyPlot

@test true
if !Sys.isapple()
    plotting=joinpath(@__DIR__,"..","examples","plotting.jl")
    include(plotting)
    include("../docs/makeplots.jl")
    @test makeplots(mktempdir())
end


text:
module OPccmt
#
# Define generic functions used in cell-to-cell mapping of dynamic
# systems.
#

export nearestCells, surroundingCells, xFromCell, assignPhi

using SparseArrays

# ===================================================================
function nearestCells(x,xb)
#
# Return the nearest cell indices, for the cells surrounding each
# entry in a vector of states.  Also the weights for a linear
# interpolation.
#
# Version:        Changes:
# --------        -------------
# 26.03.2020      Original code.
# 27.07.2021      Streamlined the code.
#
# Version:        Verification:
# --------        -------------
# 26.03.2020      Verified that it reproduces the original Octave
#                 code, used in a variety of example problems.
# 27.07.2021      Verified that it reproduces the old nearestCells
#                 output.  Checked using @btime, about 20% faster.
#
# Inputs:
# -------
# x               : Ndim-by-Nvecs array containing the vectors.
# xb              : Ndim-by-3 array containing xmin, dx, xmax for each
#                   dimension; that is, the possible values.
#
# Outputs:
# --------
# ic              : nearest cell ids, 2^Ndim-by-Nvecs.
# wc              : weights associated with each of the ic cells.

Ndim = size(x,1)
Nvecs = size(x,2)
Nx = round.(Int, (xb[:,3] - xb[:,1])./xb[:,2]) .+ 1

vx  = (x .- xb[:,1])./xb[:,2]
ixx = min.(max.(vx .+ 1,1),Nx)
ixL = min.(max.(floor.(Int,vx) .+ 1,1),Nx.-1)
wL  = 1 .- (ixx - ixL)

ic = zeros(Int,2^Ndim,Nvecs)
wc = zeros(2^Ndim,Nvecs)

for jj = 1:2^Ndim

   b = string(jj-1, base=2, pad=Ndim)
   ix = deepcopy(ixL)
   wx = deepcopy(wL)

   # The values that are 1 in binary, make "high".
   for kk in 1:Ndim
      hi = parse(Bool, b[kk])
      ix[kk,:] = ix[kk,:] .+ hi
      if (hi)
         wx[kk,:] = 1 .- wL[kk,:]
      end
   end

   wc[jj,:] = prod(wx, dims=1)

   for idim = 1:Ndim

      if (idim == Ndim)
         ic[jj,:] = ic[jj,:] + ix[idim,:]
      else
         ic[jj,:] = ic[jj,:] + prod(Nx[idim+1:Ndim]).*(ix[idim,:].-1)
      end

   end
end

ic = min.(max.(ic,1),prod(Nx))

return ic, wc

end # nearestCells

# ===================================================================
function surroundingCells(x,xb)
#
# Return the cells surrounding (adjacent to) the nearest cell.
#
# Version:        Changes:
# --------        -------------
# 27.07.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 27.07.2021      Checked a simple 2D problem, including out-of-bounds.
#
# Inputs:
# -------
# x               : Ndim-by-Nvecs array containing the vectors.
# xb              : Ndim-by-3 array containing xmin, dx, xmax for each
#                   dimension; that is, the possible values.
#
# Outputs:
# --------
# ic              : Surrounding cell ids, 3^Ndim-by-Nvecs.  In the case
#                   where the x coordinate is out of bounds, there may
#                   be repeated entries, but they are truncated to
#                   valid indices.

Ndim = size(x,1)
Nvecs = size(x,2)
Nx = round.(Int, (xb[:,3] - xb[:,1])./xb[:,2]) .+ 1
Ncel = 3^Ndim   # Ndim box of surrounding cells, including the center.

# First find the nearest cell.
jc, wc = nearestCells(x,xb)
jnear = zeros(Int,Nvecs)
for jj in 1:Nvecs
   imax = argmax(wc[:,jj])
   jnear[jj] = jc[imax,jj]
end

# Next identify the lower bound among the surrounding cells.
xnear = xFromCell(jnear,xb)
xlb = xnear .- 0.1*xb[:,2]  # Move a fraction of dx towards the LB.

# Find the lower-bound index on the surrounding cells.
vx  = (xlb .- xb[:,1])./xb[:,2]
ixx = min.(max.(vx .+ 1,1),Nx)
ixL = min.(max.(floor.(Int,vx) .+ 1,1),Nx.-1)

ic = zeros(Int, Ncel, Nvecs)

for jj in 1:Ncel

   # Base 3 number: 0 = "low", 1 = "mid", 2 = "high".
   b = string(jj-1, base=3, pad=Ndim)

   # Increment the values of ix according to the present cell ID.
   ix = deepcopy(ixL)
   for kk in 1:Ndim
      ix[kk,:] = ix[kk,:] .+ parse(Int, b[kk])
   end

   for kk = 1:Ndim

      if (kk == Ndim)
         ic[jj,:] = ic[jj,:] + ix[kk,:]
      else
         ic[jj,:] = ic[jj,:] + prod(Nx[kk+1:Ndim]).*(ix[kk,:].-1)
      end

   end

end

ic = min.(max.(ic,1),prod(Nx))

end  # surroundingCells

# ===================================================================
function xFromCell(ic,xb)
#
# Return the cell states for a given cell ID.
#
# Version:        Changes:
# --------        -------------
# 28.03.2020      Original code.
#
# Version:        Verification:
# --------        -------------
# 28.03.2020      Verified that it reproduces the original Octave
#                 code, used in a variety of example problems.
#
# Inputs:
# -------
# ic              : nearest cell id.
# xb              : N-by-3 array containing xmin, dx, xmax for each x;
#                   that is, the possible values.
#
# Outputs:
# --------
# x               : the vector of N states.

icr = deepcopy(ic)
Ndim = size(xb,1)
Nvec = size(ic,1)
Nx = round.(Int, (xb[:,3] - xb[:,1])./xb[:,2]) .+ 1
ix = zeros(Int,Ndim,Nvec)
for idim = 1:Ndim

   if (idim == Ndim)
      ix[idim,:] = icr
   else
      pr = prod(Nx[idim+1:Ndim])
      n = floor.(Int,icr./pr)
      rem = mod.(icr,pr)
      ind = (rem .== 0)
      n[ind] = n[ind] .- 1
      ix[idim,:] = n .+ 1
      icr = icr - pr.*n
   end

end

x = xb[:,1] .+ (ix .- 1).*xb[:,2]

return x

end # xFromCell

# ===================================================================
function assignPhi(jc,wc,Np)

Nj1 = size(jc,1)
Nnz = Nj1*size(jc,2)
ir = reshape(jc,Nnz)
ind = range(0,Nnz-Nj1+1,step=Nj1)
ic = zeros(Int,Nnz)
for irep = 1:Nj1
   ic[irep.+ind] = range(1,size(jc,2),step=1)
end
ss = reshape(wc,Nnz)
Phi = sparse(ir,ic,ss,Np,size(jc,2))
return Phi

end # assignPhi

end # module

max_stars_repo_name: SINTEF-Energy-Wind/OutliersPrimer

max_issues_repo_name: SINTEF-Energy-Wind/OutliersPrimer

max_forks_repo_name: SINTEF-Energy-Wind/OutliersPrimer