text:
module GRAPELinesearchAnalysis

using Printf
using QuantumControlBase.ConditionalThreads: @threadsif
using QuantumPropagators
using GRAPE
using Optim
using LBFGSB
using LinearAlgebra
using Plots
using Plots.PlotMeasures: px

# Given a search direction (e.g., the negative gradient) and a vector of step
# widths α, return a vector of J_T for going downhill in the direction of the
# gradient by a distance α, from `pulsevals0`. Plotting the resul is useful for
# getting a feeling for the gradient and the linesearch.
function explore_linesearch(search_direction, α_vals, pulsevals0, wrk)
    J_T_vals = zeros(length(α_vals))
    J_T_func = wrk.kwargs[:J_T]
    N_T = length(wrk.result.tlist) - 1
    N = length(wrk.objectives)
    Ψ = wrk.fw_states
    τ = wrk.result.tau_vals
    for (i, α) in enumerate(α_vals)
        pulsevals = pulsevals0 + α * search_direction
        @threadsif wrk.use_threads for k = 1:N
            copyto!(Ψ[k], wrk.objectives[k].initial_state)
            for n = 1:N_T  # `n` is the index for the time interval
                local (G, dt) = GRAPE._fw_gen(pulsevals, k, n, wrk)
                propstep!(Ψ[k], G, dt, wrk.fw_prop_wrk[k])
            end
            τ[k] = dot(wrk.objectives[k].target_state, Ψ[k])
        end
        J_T_vals[i] = J_T_func(Ψ, wrk.objectives; τ=τ)
    end
    return J_T_vals
end

function _get_linesearch_data(
    wrk,
    optimization_state::Optim.OptimizationState,
    optimizer_state::Optim.AbstractOptimizerState,
)
    # TODO: use wrk only (and compare with current implementation)
    α = optimizer_state.alpha # the stepwidth chosen by the linesearch
    lbfgs_direction = optimizer_state.s # the quasi-Newton update direction
    gradient_direction = -1 * optimizer_state.g_previous
    pulsevals_guess = optimizer_state.x_previous
    pulsevals_opt = optimizer_state.x

    return α, lbfgs_direction, gradient_direction, pulsevals_guess, pulsevals_opt

end

function _get_linesearch_data(
    wrk,
    iter_res::GRAPE.LBFGSB_Result,
    optimizer::LBFGSB.L_BFGS_B,
)
    # TODO: use wrk only (and compare with current implementation)
    lbfgs_direction = wrk.searchdirection
    gradient_direction = -1 * wrk.gradient
    if norm(lbfgs_direction) ≈ 0.0
        # LBFGSB hasn't yet built up Hessian (first iteration)
        # TODO: this should be fixed in GRAPE
        lbfgs_direction = gradient_direction
    end
    pulsevals_opt = iter_res.x
    pulsevals_guess = iter_res.x_previous
    Δϵ = pulsevals_opt - pulsevals_guess
    norm_Δϵ = norm(Δϵ)
    norm_ls = norm(lbfgs_direction)
    proj_ls = Δϵ ⋅ lbfgs_direction
    cosθ = proj_ls / (norm_Δϵ * norm_ls)
    @assert cosθ ≈ 1.0  # search direction and Δϵ are parallel
    α = proj_ls / norm_ls^2
    @assert α ≈ optimizer.dsave[14]
    @assert norm(Δϵ - α * lbfgs_direction) < 1e-12

    return α, lbfgs_direction, gradient_direction, pulsevals_guess, pulsevals_opt

end

function plot_linesearch(
    outdir;
    samples=100,
    verbose=(get(ENV, "GRAPE_LINESEARCH_ANALYSIS_VERBOSE", "0") == "1"),
    use_current_backend=(
        get(ENV, "GRAPE_LINESEARCH_ANALYSIS_USE_CURRENT_BACKED", "0") == "1"
    ),
    kwargs...
)

    defaults = Dict(
        :linewidth => 3,
        :marker_size => 3,
        :foreground_color_legend => nothing,
        :background_color_legend => RGBA(1, 1, 1, 0.8)
    )
    merge!(default, kwargs)

    mkpath(outdir)

    function _plot_linesearch(wrk, optimization_state, optimizer_state, iteration, args...)

        current_backend = backend()
        if !use_current_backend && !(current_backend isa Plots.GRBackend)
            gr()
        end

        (iteration == 0) && (return nothing)

        # TODO: linesearch data sources should be entirely in wrk, so we don't
        # need different methods.
        α, lbfgs_direction, gradient_direction, pulsevals_guess, pulsevals_opt =
            _get_linesearch_data(wrk, optimization_state, optimizer_state)
        α_vals = collect(range(0, 2α, length=samples))
        J_α_gradient = explore_linesearch(gradient_direction, α_vals, pulsevals_guess, wrk)
        J_α_lbfgs = explore_linesearch(lbfgs_direction, α_vals, pulsevals_guess, wrk)
        J_T = wrk.result.J_T
        J_T_guess = wrk.result.J_T_prev

        ax1 = plot(pulsevals_guess; label="guess", defaults...)
        plot!(
            ax1,
            pulsevals_opt;
            label="opt",
            xlabel="time step (control parameter index)",
            ylabel="control amplitude",
            defaults...
        )
        plot!(
            ax1,
            title="Iteration $iteration; J = $(@sprintf("%.2e", J_T_guess)) → $(@sprintf("%.2e", J_T))"
        )

        ax2 = plot(
            gradient_direction;
            label="gradient",
            xlabel="time step (control parameter index)",
            ylabel="gradient directon",
            defaults...
        )

        ax3 = plot(gradient_direction; linestyle=:dash, label="gradient", defaults...)
        plot!(
            ax3,
            lbfgs_direction;
            label="LBFGS direction",
            xlabel="time step (control parameter index)",
            ylabel="LBFGS direction",
            defaults...
        )

        ax4 = plot(
            α_vals,
            J_α_gradient;
            label="gradient",
            xlabel="step width α",
            ylabel="Functional",
            title="Gradient Linesarch",
            defaults...
        )
        scatter!(ax4, [0,], [J_T_guess,], label="")
        hline!(ax4, [J_T_guess], color="black", linewidth=0.5, label="")
        vline!(ax4, [α], color="black", linewidth=0.5, label="")

        ax5 = plot(
            α_vals,
            J_α_gradient;
            linestyle=:dash,
            label="gradient",
            title="LBFGS Linesearch α = $α",
            defaults...
        )
        plot!(
            ax5,
            α_vals,
            J_α_lbfgs;
            label="LBFGS",
            xlabel="step width α",
            ylabel="Functional",
            defaults...
        )
        scatter!(ax5, [0, α], [J_T_guess, J_T], label="")
        hline!(ax5, [J_T_guess], color="black", linewidth=0.5, label="")
        hline!(ax5, [J_T], color="black", linewidth=0.5, label="")
        vline!(ax5, [α], color="black", linewidth=0.5, label="")
        if J_T < 0.1
            plot!(ax5; yaxis=:log)
            ylim = ylims(ax5)
            ylims!(ax5, (ylim[1], 1))
        else
            ylim = ylims(ax5)
            ylims!(ax5, (ylim[1], 1))
        end

        fig = plot(ax1, ax2, ax3, ax4, ax5, layout=(5, 1),)
        plot!(fig, size=(700, 1500), left_margin=50px,)

        outfile = joinpath(outdir, @sprintf("linesearch-%03i.png", iteration))
        savefig(fig, outfile)
        if verbose
            @info "Written GRAPELinesearchAnalysis plot to $outfile"
        end

        if !use_current_backend && !(current_backend isa Plots.GRBackend)
            current_backend() # switch back
        end

        return nothing

    end

    return _plot_linesearch
end

end # module


text:
function majority(s::Array)
    # Create a counter dictionary Dict(element => number of times element appears in input array)
    counter = Dict()
    length_s = 0
    # Count the number of times each element in s appears
    for i in s
       length_s += 1
       if i in keys(counter)
           counter[i] += 1
       else
           counter[i] = 1
       end
    end
    # Loop through the items in the counter dictionary, return majority element
    for (key, value) in counter
       if value > length_s/2
           return key
       end
    end
    # If no majority element, return a message:
    return "No majority element"
end

# The example case from the prompt
println(majority([2,2,3,7,5,7,7,7,4,7,2,7,4,5,6,7,7,8,6,7,7,8,10,12,29,30,19,10,7,7,7,7,7,7,7,7,7]))
# And another case
println(majority([1,1,1,2,3,4,1,1,1]))
# And a case with no majority element
println(majority([9,9,9,9,8,8,8,8]))

ext: jl


text:
import Logging
@Logging.configure(level=DEBUG, output=STDOUT)


text:
"""
    NDAS{N,S} <: DAS
wraps multiple independent symmetries to work together.

# Example
```julia-repl
julia>NDAS(U1(), Z2())
NDAS{2,(U1, Z2)}()
```
"""
struct NDAS{N,S} <: DAS end
NDAS(s...) where S<:Tuple = NDAS{length(s),s}()

"""
    NDASCharges{N,T} <: DASCharges
holds combinations of multiple independent `DASCharges` that are grouped together.
Yields `NDASCharge{N,T}`s upon iteration.

# Example
```julia-repl
julia> a = NDASCharges(U1Charges(-1:1), ZNCharges{3}());
julia> foreach(println, a)
NDASCharge(U1Charge(-1), Z3Charge(0))
NDASCharge(U1Charge(0), Z3Charge(0))
NDASCharge(U1Charge(1), Z3Charge(0))
NDASCharge(U1Charge(-1), Z3Charge(1))
NDASCharge(U1Charge(0), Z3Charge(1))
NDASCharge(U1Charge(1), Z3Charge(1))
NDASCharge(U1Charge(-1), Z3Charge(2))
NDASCharge(U1Charge(0), Z3Charge(2))
NDASCharge(U1Charge(1), Z3Charge(2))
```
"""
struct NDASCharges{N,S} <: DASCharges
    v::S
    NDASCharges(v...) = new{length(v),typeof(v)}(v)
end
Base.show(io::IO, s::NDASCharges) = print(io,"NDASCharges", s.v)

⊕(a::T, b::T) where {T<:NDASCharges{N}} where N = NDASCharges(map(⊕, a.v, b.v)...)
Base.inv(a::T) where {T<:NDASCharges} = NDASCharges(inv.(a.v)...)
Base.length(a::NDASCharges) = prod(length.(a.v))
Base.intersect(a::NDASCharges{N}, b::NDASCharges{N}) where N = NDASCharges(map(intersect, a.v, b.v)...)
function Base.getindex(a::NDASCharges, i::Int)
    i <= length(a) || throw(BoundsError(a,i))
    inds = Tuple(CartesianIndices(map(x -> 1:length(x), a.v))[i])
    NDASCharge(map(getindex, a.v, inds)...)
end
chargestype(::NDAS{N,S}) where {N,S} =
    NDASCharges{N,Tuple{(chargestype(TT) for TT in S)...}}

"""
    NDASCharge{N} <: DASCharge
holds the charge of multiple independent symmetries grouped together.

# Example
```julia-repl
julia> a = NDASCharge(U1Charge(1), Z2Charge(1));
julia> a ⊕ a
NDASCharge(U1Charge(2), Z2Charge(0))
```
"""
struct NDASCharge{N,S} <: DASCharge
    ch::S
    NDASCharge(ch...) = new{length(ch),typeof(ch)}(ch)
    NDASCharge{N,S}(ch...) where {N,S} = new{N,S}(ch)
end
Base.show(io::IO, s::NDASCharge) = print(io,"NDASCharge", s.ch)

Base.eltype(::NDASCharges{N,S}) where {N,S} =
    NDASCharge{N,Tuple{(eltype(T) for T in S.types)...}}
Base.iszero(a::NDASCharge) = all(iszero,a.ch)
⊕(a::T, b::T) where {T<:NDASCharge} = NDASCharge(map(⊕, a.ch, b.ch)...)
Base.:-(a::T) where {T<:NDASCharge} = T(map(x -> -x, a.ch)...)
Base.isless(a::T, b::T) where {T<:NDASCharge{N}} where N = a.ch < b.ch
Base.zero(::Type{NDASCharge{N,T}}) where {N,T} = NDASCharge{N,T}(ntuple(i -> zero(T.types[i]), N)...)
Base.zero(::NDASCharge{N,T}) where {N,T} = NDASCharge{N,T}(ntuple(i -> zero(T.types[i]), N)...)
Base.:+(chs::T, ch::S) where {T<:NDASCharges, S<:NDASCharge} = NDASCharges(map(+, chs.v, ch.ch)...)
function chargeindex(ch::NDASCharge, chs::NDASCharges)
    ch in chs || throw(ArgumentError(string(ch, " not in ", chs)))
    linds = LinearIndices(map(x -> 1:x, length.(chs.v)))
    inds = map(chargeindex, ch.ch, chs.v)
    linds[inds...]
end
chargetype(::NDAS{N,S}) where {N,S} = NDASCharge{N,Tuple{(chargetype(TT) for TT in S)...}}


text:
# Euler's Totient function, φ(n) [sometimes called the phi function], is used
# to determine the number of numbers less than n which are relatively prime to
# n. For example, as 1, 2, 4, 5, 7, and 8, are all less than nine and relatively
# prime to nine, φ(9)=6.
#
# n   Relatively Prime  φ(n)   n/φ(n)
# --------------------------------------
# 2   1                 1      2
# 3   1,2               2      1.5
# 4   1,3               2      2
# 5   1,2,3,4           4      1.25
# 6   1,5               2      3
# 7   1,2,3,4,5,6       6      1.1666...
# 8   1,3,5,7           4      2
# 9   1,2,4,5,7,8       6      1.5
# 10  1,3,7,9           4      2.5
#
# It can be seen that n=6 produces a maximum n / φ(n) for n ≤ 10.
#
# Find the value of n ≤ 1,000,000 for which n / φ(n) is a maximum.

using ProjectEulerSolutions

# Brute force solution, considering every value less than n.
function p069solution_bruteforce(n::Integer=1_000)::Integer

    max_val = 0.0
    max_index = 0
    for i in 2:n
        test_val = i / totient(i)
        if max_val < test_val
            max_val = test_val
            max_index = i
        end
    end

    return max_index
end

# Only consider values that are likely candidates, which are primodial numbers.
# Primodial numbers are the product of consecutive primes starting at 2, so
# 2, 6, 30, 210, etc.  We only consider primordial numbers up to that who's
# product more than n since they have maximal n / φ(n) values.
function p069solution_factors(n::Integer=1_000)::Integer

    primes = sieve_eratosthenes(30)
    index = 2
    p = primes[1]
    while p * primes[index] < n
        p *= primes[index]
        index += 1
    end
    return p
end

p069 = Problems.Problem(Dict("Brute force" => p069solution_bruteforce,
                             "Factors" => p069solution_factors))

Problems.benchmark(p069, 1_000_000)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
module TestIterators

using Compat.Iterators
using Compat: Compat, Iterators
using Test

# https://github.com/JuliaLang/julia/pull/33437
@testset "takewhile" begin
    @test collect(takewhile(<(4),1:10)) == [1,2,3]
    @test collect(takewhile(<(4),Iterators.countfrom(1))) == [1,2,3]
    @test collect(takewhile(<(4),5:10)) == []
    @test collect(takewhile(_->true,5:10)) == 5:10
    @test collect(takewhile(isodd,[1,1,2,3])) == [1,1]
    @test collect(takewhile(<(2), takewhile(<(3), [1,1,2,3]))) == [1,1]
end

# https://github.com/JuliaLang/julia/pull/33437
@testset "dropwhile" begin
    @test collect(dropwhile(<(4), 1:10)) == 4:10
    @test collect(dropwhile(<(4), 1:10)) isa Vector{Int}
    @test isempty(dropwhile(<(4), []))
    @test collect(dropwhile(_->false,1:3)) == 1:3
    @test isempty(dropwhile(_->true, 1:3))
    @test collect(dropwhile(isodd,[1,1,2,3])) == [2,3]
    @test collect(dropwhile(iseven,dropwhile(isodd,[1,1,2,3]))) == [3]
end

# https://github.com/JuliaLang/julia/pull/34352
@testset "Iterators.map" begin
    @test collect(Iterators.map(string, 1:3)::Base.Generator) == map(string, 1:3)
    @test collect(Iterators.map(tuple, 1:3, 4:6)::Base.Generator) == map(tuple, 1:3, 4:6)
end

@testset "Iterators.filter" begin
    # `Iterators.filter` already existed in Julia 1.0.  We just make
    # sure that it is available under `Compat.Iterators` namespace.
    @test Compat.Iterators.filter === Base.Iterators.filter
end

end  # module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule LevelDB_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("LevelDB")
JLLWrappers.@generate_main_file("LevelDB", UUID("e5add9ea-9f03-521a-be84-f8c1bde3cbe9"))
end  # module LevelDB_jll


text:
# Michelle Medina
# September 16, 2020

function derivBasis(l::Int, m::Int, i::Int, j::Int, x::T, y::T)::T where {T}
	return derivBasis(l, i, x)*derivBasis(m, j, y)
end

function derivEvaluate(H::HierarchicalBasis{T}, x::T, y::T)::T where {T}
	value = T(0)
	for l in 0:maxlevel(H), for m in 0:maxlevel(H)
		for i in 0:2^l, for j in 0:2^m
			value = value + H[l, m, i, j]*derivBasis(l, m, i, j, x, y)
		end
	end
	return value
end


text:
# ----------------------------------------------------------------------
# The repulsive core is built from shifted Buckingham potentials

@doc raw"""
`struct BuckPot` : Buckingham potential,
```math
   V(r) = e_0 + B  \exp\big( - A (r/r_i - 1)\big) \cdot \frac{r_i}{r}
```
"""
struct BuckPot{T} <: SimplePairPotential
   e0::T
   A::T
   ri::T
   B::T
end

@pot BuckPot

evaluate(V::BuckPot, r::Number) = V.e0 + V.B * exp( - V.A * (r/V.ri-1) ) * V.ri/r

evaluate_d(V::BuckPot, r::Number) = V.B * exp( - V.A * (r/V.ri-1) ) * V.ri * (
                            - V.A / V.ri / r  - 1/r^2  )

write_dict(V::BuckPot) = Dict(
         "__id__" => "ACE_BuckPot",
             "e0" => V.e0,
              "A" => V.A,
             "ri" => V.ri,
              "B" => V.B   )

read_dict(::Val{:SHIPs_BuckPot}, D::Dict) =
   read_dict(Val{:ACE_BuckPot}(), D)

read_dict(::Val{:ACE_BuckPot}, D::Dict) =
      BuckPot(D["e0"], D["A"], D["ri"], D["B"])

==(V1::BuckPot, V2::BuckPot) = _allfieldsequal(V1, V2)

# ----------------------------------------------------------------------

struct RepulsiveCore{T, TOUT, NZ} <: PairPotential
   Vout::TOUT
   Vin::SMatrix{NZ, NZ, BuckPot{T}}
end

@pot RepulsiveCore

==(V1::RepulsiveCore, V2::RepulsiveCore) = _allfieldsequal(V1, V2)

cutoff(V::RepulsiveCore) = cutoff(V.Vout)

alloc_temp(V::RepulsiveCore, N::Integer) = alloc_temp(V.Vout, N)
alloc_temp_d(V::RepulsiveCore, N::Integer) = alloc_temp_d(V.Vout, N)

function evaluate!(tmp, V::RepulsiveCore, r::Number, z, z0)
   iz, iz0 = z2i(V.Vout, z), z2i(V.Vout, z0)
   Vin = V.Vin[iz, iz0]
   if r > Vin.ri
      return evaluate!(tmp, V.Vout, r, z, z0)
   else
      return evaluate(Vin, r)
   end
end

function evaluate_d!(tmp, V::RepulsiveCore, r::Number, z, z0)
   iz, iz0 = z2i(V.Vout, z), z2i(V.Vout, z0)
   Vin = V.Vin[iz, iz0]
   if r > Vin.ri
      return evaluate_d!(tmp, V.Vout, r, z, z0)
   else
      return evaluate_d(Vin, r)
   end
end

function _simple_repulsive_core(Vout, ri, e0, verbose, z, z0)
   v = Vout(ri, z, z0)
   dv = @D Vout(ri, z, z0)
   if dv >= 0.0 && verbose
      @warn("The slope `Vout'(ri)` should be negative")
   end
   if dv > -1.0 && verbose
      @warn("""The slope `Vout'(ri) = $dv` may not be steep enough to attach a
               repulsive core. Proceed at your own risk.""")
   end
   if v-e0 <= 0.0 && verbose
      @warn("it is recommended that `Vout(ri) > 0`.")
   end
   if v-e0 <= 1.0 && verbose
      @warn("""Ideally the repulsive core should not be attached at small
               values of `Vout(ri) = $v`. Proceed at your own risk.""")
   end
   # e0 + B e^{-A (r/ri-1)} * ri/r
   #    => e0 + B = Vout(ri) => = Vout(ri) - e0 = v - e0
   # dv = - A*B/ri e^{-A (r/ri-1)} * ri/r - B*ri*e^{...} / r^2
   #    = - A/ri * (v - 1/ri * (v = - (1+A)/ri * (v-e0)
   #    => -(1+A)/ri * (v-e0) = dv
   #    => 1+A = - ri dv / (v-e0)
   #    => A = -1 - ri dv / (v-e0)
   A = -1 - ri * dv / (v-e0)
   B = v-e0
   Vin = BuckPot(e0, A, ri, B)
   if verbose
      @show ri
      @show Vout(ri), (@D Vout(ri))
      @show Vin(ri), (@D Vin(ri))
   end
   return Vin
end

function RepulsiveCore(Vout, ri::Number, e0=0.0; verbose=false)
   nz = numz(Vout)
   Vin = Matrix{Any}(undef, nz, nz)
   for i0 = 1:nz, i1 = 1:i0
      z0, z1 = i2z(Vout, i0), i2z(Vout, i1)
      Vin[i0, i1] = Vin[i1, i0] =
            _simple_repulsive_core(Vout, ri, e0, verbose, z1, z0)
   end
   # construct the piecewise potential
   return RepulsiveCore(Vout, Vin)
end

function RepulsiveCore(Vout, Vin::AbstractArray)
   nz = numz(Vout)
   return RepulsiveCore(Vout, SMatrix{nz, nz}(Vin...))
end

function RepulsiveCore(Vout, D::Dict; verbose=false)
   _sort_key(V, key) = tuple(sort([z2i.(Ref(V), AtomicNumber.(key))...])...)
   Di = Dict([_sort_key(Vout, key) => val for (key, val) in D]...)
   nz = numz(Vout)
   Vin = Matrix{Any}(undef, nz, nz)
   for i0 = 1:nz, i1 = 1:i0
      z0, z1 = i2z(Vout, i0), i2z(Vout, i1)
      ri, e0 = Di[(i1, i0)].ri, Di[(i1, i0)].e0
      Vin[i0, i1] = Vin[i1, i0] =
            _simple_repulsive_core(Vout, ri, e0, verbose, z1, z0)
   end
   # construct the piecewise potential
   return RepulsiveCore(Vout, Vin)
end

# ----------------------------------------------------
#  File IO

write_dict(V::RepulsiveCore) = Dict(
      "__id__" => "ACE_RepulsiveCore",
        "Vout" => write_dict(V.Vout),
         "Vin" => write_dict.(V.Vin[:]) )

function read_dict(::Val{:ACE_RepulsiveCore}, D::Dict)
   Vout = read_dict(D["Vout"])
   nz = numz(Vout)
   Vin = read_dict.(D["Vin"])
   return RepulsiveCore(Vout, reshape(Vin, (nz, nz)))
end


text:
push!(LOAD_PATH, "/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/main")
using DHC_2DUtils
using Profile
using BenchmarkTools
using FFTW
using HDF5
using Test
using SparseArrays
using Statistics
using Plots
using LinearAlgebra
using Measures
using Plots
using MLDatasets
using Images
theme(:juno)
using ImageFiltering

filter_hash = fink_filter_hash(1, 8, nx=256, pc=1, wd=2)

train_x, train_y = CIFAR10.traindata();
test_x, test_y   = CIFAR10.testdata();

using TestImages
using ImageMetadata
using ImageAxes
using ImageDraw
using ImageView
using ImageTransformations
using ImageCore, Colors

ind = 10
s = StackedView(transpose(train_x[:,:,1,ind]),transpose(train_x[:,:,2,ind]),transpose(train_x[:,:,3,ind]));
sc = colorview(RGB, s)

image_x = train_x[:,:,:,ind]
size(train_x[:,:,:,ind])

mean(image_x)

@benchmark mean(image_x, dims=(1,2))

mean_x = dropdims(mean(image_x, dims=(1,2)),dims=(1,2))
mean_x_shaped = mean(image_x, dims=(1,2))
norm_im = image_x .- mean_x_shaped

pwr_im = sum(norm_im .* norm_im,dims=(1,2))
norm_im./sqrt.(32*32 .*pwr_im)

import AbstractFFTs

@benchmark P0 = plan_fft(norm_im,(1,2))
@benchmark im_fd_0 = fft(norm_im,(1,2))
@benchmark im_fd_0 = P0*norm_im

S1M = filter_hash["S1_iso_mat"]
Mtest = [S1M 0I 0I; 0I S1M 0I; 0I 0I S1M]
blockdiag(S1M,S1M,S1M)

1+1

function DHC_compute_RGB(image::Array{Float64}, filter_hash::Dict, filter_hash2::Dict=filter_hash;
    doS20::Bool=true, norm=true, iso=false, FFTthreads=2)
    # image        - input for WST
    # filter_hash  - filter hash from fink_filter_hash
    # filter_hash2 - filters for second order.  Default to same as first order.
    # doS2         - compute S2 coeffs
    # doS12        - compute S2 coeffs
    # doS20        - compute S2 coeffs
    # norm         - scale to mean zero, unit variance
    # iso          - sum over angles to obtain isotropic coeffs

    # Use 2 threads for FFT
    FFTW.set_num_threads(FFTthreads)

    # array sizes
    (Nx, Ny, Nc)  = size(image)
    if Nx != Ny error("Input image must be square") end
    (Nf, )    = size(filter_hash["filt_index"])
    if Nf == 0  error("filter hash corrupted") end
    @assert Nx==filter_hash["npix"] "Filter size should match npix"
    @assert Nx==filter_hash2["npix"] "Filter2 size should match npix"

    # allocate coeff arrays
    out_coeff = []
    S0  = zeros(Float64, Nc*2)
    S1  = zeros(Float64, Nc*Nf)

    if doS20 S20 = zeros(Float64, Nf, Nf) end  # real space correlatio
    # anyM2 = doS2 | doS12 | doS20
    anyrd = doS20 #| doS2             # compute real domain with iFFT

    # allocate image arrays for internal use
    mean_im = zeros(Float64,1,1,Nc)
    pwr_im = zeros(Float64,1,1,Nc)
    norm_im = zeros(Float64,Nx,Ny,Nc)
    im_fd_0 = zeros(ComplexF64, Nx, Ny, Nc)
    im_fd_0_sl = zeros(ComplexF64, Nx, Ny)

    if doS20
        Amat1 = zeros(Nx*Ny, Nf)
        Amat2 = zeros(Nx*Ny, Nf)
    end

    if anyrd im_rd_0_1  = Array{Float64, 4}(undef, Nx, Ny, Nf, Nc) end

    ## 0th Order
    mean_im = mean(image, dims=(1,2))
    S0[1:Nc]   = dropdims(mean_im,dims=(1,2))
    norm_im = image.-mean_im
    pwr_im = sum(norm_im .* norm_im,dims=(1,2))
    S0[1+Nc:end]   = dropdims(pwr_im,dims=(1,2))./(Nx*Ny)
    if norm
        norm_im ./= sqrt.(pwr_im)
    else
        norm_im = copy(image)
    end

    append!(out_coeff,S0[:])

    ## 1st Order
    im_fd_0 .= fft(norm_im,(1,2))  # total power=1.0

    # unpack filter_hash
    f_ind   = filter_hash["filt_index"]  # (J, L) array of filters represented as index value pairs
    f_val   = filter_hash["filt_value"]

    zarr = zeros(ComplexF64, Nx, Ny)  # temporary array to fill with zvals

    # make a FFTW "plan" for an array of the given size and type
    if anyrd
        P = plan_ifft(zarr) end  # P is an operator, P*im is ifft(im)

    ## Main 1st Order and Precompute 2nd Order
    for f = 1:Nf
        S1tot = 0.0
        f_i = f_ind[f]  # CartesianIndex list for filter
        f_v = f_val[f]  # Values for f_i
        # for (ind, val) in zip(f_i, f_v)   # this is slower!
        for chan = 1:Nc
            im_fd_0_sl .= im_fd_0[:,:,chan]
            if length(f_i) > 0
                for i = 1:length(f_i)
                    ind       = f_i[i]
                    zval      = f_v[i] * im_fd_0_sl[ind]
                    S1tot    += abs2(zval)
                    zarr[ind] = zval        # filter*image in Fourier domain
                end
                S1[f+(chan-1)*Nf] = S1tot/(Nx*Ny)  # image power
                if anyrd
                    im_rd_0_1[:,:,f,chan] .= abs2.(P*zarr) end
                zarr[f_i] .= 0
            end
        end
    end

    if iso
        S1M = filter_hash["S1_iso_mat"]
        M1 = blockdiag(S1M,S1M,S1M)
    end
    append!(out_coeff, iso ? M1*S1 : S1)

    # we stored the abs()^2, so take sqrt (this is faster to do all at once)
    if anyrd im_rd_0_1 .= sqrt.(im_rd_0_1) end

    # Real domain 2nd order
    if doS20
        for chan1 = 1:Nc
            for chan2 = 1:Nc
                Amat1 = reshape(im_rd_0_1[:,:,:,chan1], Nx*Ny, Nf)
                Amat2 = reshape(im_rd_0_1[:,:,:,chan2], Nx*Ny, Nf)
                S20  = Amat1' * Amat2
                append!(out_coeff, iso ? filter_hash["S2_iso_mat"]*S20[:] : S20[:])
            end
        end
    end

    return out_coeff
end

DHC_compute_RGB(train_x[:,:,:,ind],filter_hash)

using Interpolations
using DSP

function wind_2d_RGB(nx)
    dx   = nx/2-1
    filter = zeros(Float64, nx, nx)
    A = DSP.tukey(nx, 0.3)
    itp = extrapolate(interpolate(A,BSpline(Linear())),0)
    @inbounds for x = 1:nx
        sx = x-dx-1    # define sx,sy so that no fftshift() needed
        for y = 1:nx
            sy = y-dx-1
            r  = sqrt.((sx).^2 + (sy).^2) + nx/2
            filter[x,y] = itp(r)
        end
    end
    return reshape(filter,nx,nx,1)
end

wind_2d_RGB(256)

function apodizer(data,sp1,sp2,im_size)
    temp2d = data[sp1:sp1+im_size,sp2:sp2+im_size]
    datad_w = fweights(wind_2d(256));
    meanVal = mean(temp2d,datad_w)
    temp2d_a = (temp2d.-meanVal).*wind_2d(256).+meanVal
    return temp2d_a
end

function wind_2d(nx)
    dx   = nx/2-1
    filter = zeros(Float64, nx, nx)
    A = DSP.tukey(nx, 0.3)
    itp = extrapolate(interpolate(A,BSpline(Linear())),0)
    @inbounds for x = 1:nx
        sx = x-dx-1    # define sx,sy so that no fftshift() needed
        for y = 1:nx
            sy = y-dx-1
            r  = sqrt.((sx).^2 + (sy).^2) + nx/2
            filter[x,y] = itp(r)
        end
    end
    return filter
end

using StatsBase

temp = imresize(train_x[:,:,:,ind],(64,64,3))
mean(convert(Array{Float64,3},temp),dims=(1,2))
mean(temp,dims=(1,2))
convert(Array{Float64,3},temp)

test = zeros(64,64,3)
test[:,:,1] = wind_2d_RGB(64)
test[:,:,2] = wind_2d_RGB(64)
test[:,:,3] = wind_2d_RGB(64)

datad_w = fweights(wind_2d_RGB(64));
datad_w = fweights(test);
mean(temp,reshape(datad_w,64,64,3),dims=1)
mean(temp,datad_w)

size(temp)
size(reshape(datad_w,64,64,3))

function cifar_pad_RGB(im; θ=0.0)
    imbig = convert(Array{Float64,3},imresize(im,(64,64,3)))
    datad_w = fweights(wind_2d(64));
    mu_imbig = zeros(1,1,3)
    for chan = 1:3
        mu_imbig[chan] = mean(imbig[:,:,chan],datad_w)
    end
    imbig .-= mu_imbig
    imbig .*= wind_2d_RGB(64)
    impad = zeros(Float64,128,128,3)
    impad[96:-1:33,33:96,:] = imbig

    if θ != 0.0
        imrot = imrotate(impad, θ, axes(impad), Cubic(Throw(OnGrid())))
        imrot[findall(imrot .!= imrot)] .= 0.0
        return imrot .+ mu_imbig
    end

    return impad.+ mu_imbig
end

out_test = cifar_pad_RGB(train_x[:,:,:,ind])

function mnist_pad(im; θ=0.0)
    imbig = convert(Array{Float64,2},imresize(im,(64,64)))
    mu_imbig = mean(imbig)
    imbig .-= mu_imbig
    imbig .*= wind_2d(64)
    impad = zeros(Float64,128,128)
    impad[96:-1:33,33:96] = imbig

    if θ != 0.0
        imrot = imrotate(impad, θ, axes(impad), Cubic(Throw(OnGrid())))
        imrot[findall(imrot .!= imrot)] .= 0.0
        return imrot .+ mu_imbig
    end

    return impad.+ mu_imbig
end

d2_test = mnist_pad(train_x[:,:,1,ind])
@time d2_test = mnist_pad(train_x[:,:,1,ind])

s = StackedView(transpose(out_test[:,:,1]),transpose(out_test[:,:,2]),transpose(out_test[:,:,3]));
sc = colorview(RGB, s)

@time filter_hash = fink_filter_hash(1, 8, nx=128, pc=1, wd=2)

DHC_compute(d2_test,filter_hash)
DHC_compute_RGB(out_test,filter_hash)

@time DHC_compute(d2_test,filter_hash)

@time DHC_compute_RGB(out_test,filter_hash,doS20=true)

chan_img = channelview(train_x[:,:,:,ind])
cv = colorview(RGB, StackedView(train_x[:,:,1,ind],train_x[:,:,2,ind],train_x[:,:,3,ind]))
convert(RGB, HSL(270, 0.5, 0.5))
y_im = YCbCr.(cv)
channels = permutedims(channelview(float.(y_im)),(2,3,1))

cifar_DHC_out = h5read("../DHC/scratch_AKS/data/cifar10_RGB.h5", "main/train_data")

function transformMaker2(coeff, S1Mat, S2Mat; Nc=1)
    NS1 = size(S1Mat)[2]
    NS2 = size(S2Mat)[2]
    if Nc==1
        S0iso = coeff[:,1:2]
        S1iso = transpose(S1Mat*transpose(coeff[:,2+1:2+NS1]))
        S2iso = transpose(S2Mat*transpose(coeff[:,2+NS1+1:2+NS1+NS2]))
    else
        S0iso = coeff[:,1:2*Nc]
        S1MatChan = blockdiag(collect(Iterators.repeated(S1Mat,Nc))...)
        S2MatChan = blockdiag(collect(Iterators.repeated(S2Mat,Nc*Nc))...)
        S1iso = transpose(S1MatChan*transpose(coeff[:,2*Nc+1:2*Nc+Nc*NS1]))
        S2iso = transpose(S2MatChan*transpose(coeff[:,2*Nc+Nc*NS1+1:end]))
    end
    return hcat(S0iso,S1iso,S2iso)
end

transformMaker2(cifar_DHC_out',filter_hash["S1_iso_mat"],filter_hash["S2_iso_mat"],Nc=3)

Nc = 3
NS1 = size(filter_hash["S1_iso_mat"])[2]
cifar_DHC_out[2*Nc+Nc*NS1+1:end,:]

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
module Geometry

export rotation_matrix, vector_components, line_coordinates, circular_curve, triangle_area, triangle_centroid

#counterclockwise, right hand rule
function rotation_matrix(θ)

    R = [cos(θ) -sin(θ)
         sin(θ) cos(θ)]

end

function vector_components(ΔL, θ)

    num_lines = length(ΔL)

    Δxy = zeros(Float64, (num_lines, 2))

    for i = 1:num_lines

        Δxy[i,:] = [ΔL[i] * cos(deg2rad(θ[i]))  ΔL[i] * sin(deg2rad(θ[i]))]

    end

    return Δxy

end

function line_coordinates(Δxy, closed_or_open)

    xcoords = []
    ycoords = []

    num_lines = size(Δxy)[1]

    for i = 1:num_lines

        if i==1
            xcoords = [0.0, Δxy[i,1]]
            ycoords = [0.0, Δxy[i,2]]
        else

            if ((closed_or_open == 0) & (i == num_lines))
                #avoid last node for closed section
            else
                xcoords = [xcoords; xcoords[end] .+ Δxy[i,1]]
                ycoords = [ycoords; ycoords[end] .+ Δxy[i,2]]
            end
        end

    end

    return xcoords, ycoords

end

function circular_curve(radius, interior_angle, xy_PI, PI_unit_normal, n)

    #follow this nomenclature for a horizontal curve
    #https://www.cpp.edu/~hturner/ce220/circular_curves.pdf

    #interior angle
    θ = deg2rad(interior_angle)

    #angle over which the arc sweeps
    Δ = π - θ

    #define orthogonal distance from curve to PI
    E = radius*(1/(cos(Δ/2))-1)

    #define distance along tangent from PI to BC and EC
    T = radius*tan(Δ/2)

    #calculate vector that runs tangent to curve through BC
    #rotate PI_normal vector

    # if concave_up_or_down == 1
        BC_rotation = Δ + (π/2 - Δ/2)
    # elseif concave_up_or_down == 0
        # BC_rotation = -(Δ + (π/2 - Δ/2))
    # end
    R_BC = rotation_matrix(BC_rotation)
    BC_unit_tangent = -R_BC*PI_unit_normal  #- sign is to flip direction of vector, face towards direction of curve calculation

    #calculate vector that runs tangent to curve through EC
    #rotate PI_normal vector
    #if concave_up_or_down == 1
        EC_rotation = -(Δ + (π/2 - Δ/2))
    #elseif concave_up_or_down == 0
    #    EC_rotation = Δ + (π/2 - Δ/2)
    #end
    R_EC = rotation_matrix(EC_rotation)
    EC_unit_tangent = R_EC*PI_unit_normal

    #define x-y coordinates of BC
    xy_BC = xy_PI - T * BC_unit_tangent

    #define x-y coordinates of EC
    xy_EC = xy_PI + T * EC_unit_tangent

    #define x-y coordinates at center of circle
    xy_o = xy_PI - (E+radius) * PI_unit_normal

    #define unit radius vector going through BC
    #this is the unit vector normal to the tangent vector at BC
    radius_unit_vector_BC = [-BC_unit_tangent[2], BC_unit_tangent[1]]

    #rotate unit radius vector
    Δ_sweep = range(0.0, -Δ, length=n+1)
    xy_curve = zeros(Float64, (n+1, 2))
    xy_curve[1,:] = xy_BC
    xy_curve[end,:] = xy_EC

    for i = 2:n
        R = rotation_matrix(Δ_sweep[i])
        xy_curve[i,:] = xy_o + radius*(R*radius_unit_vector_BC)
    end

    return xy_curve, Δ, E, T, xy_o, xy_BC, xy_EC, BC_unit_tangent, EC_unit_tangent, radius_unit_vector_BC

end

# https://keisan.casio.com/has10/SpecExec.cgi?path=05000000.Mathematics%252F01000500.Plane%2520geometry%252F10010300.Area%2520of%2520a%2520triangle%2520with%2520three%2520points%252Fdefault.xml&charset=utf-8
function triangle_area(x1, y1, x2, y2, x3, y3)

    A = abs((x1*y2 + x2*y3 + x3*y1 - y1*x2 - y2*x3 - y3*x1)/2)

end

#https://www.mathopenref.com/coordcentroid.html
function triangle_centroid(x1, y1, x2, y2, x3, y3)

    cx = (x1 + x2 + x3)/3
    cy = (y1 + y2 + y3)/3

    return cx, cy

end

end #module


text:
export AbstractMeganetElement, nTheta, nFeatIn, nFeatOut,nDataOut, splitWeights, initTheta
export JYTmv, JYmv, Jthetamv, JthetaTmv, Jmv, JTmv
export linearizeY, linearizeTheta
export apply
export getJOp, getJYOp, getJthetaOp

"""
 here go all functions that are shared among elements of a Meganet
 (i.e., transformations, layers, networks, blocks,...)

 An AbstractMeganetElement does the following

  Y_k+1 = apply(theta,Y_k)

 where 'apply' can be everything from a single affine transformation to a
 ResNet block. All these operations have a similar structure, e.g.,
 provide derivatives w.r.t. theta and Y_k, ...

 Example: Consider a Neural network consisting of two layers.
 T1  = dense([12,8]);
 T2  = dense([24,12])
 net = NN({T1, T2});

 Calling apply(net,theta,Y) results in a nested evaluation

 apply(net,theta,Y) = apply(T2, theta2, apply(T1, theta1, Y));

 This example shows that each element of the network needs the following
 functions

  split - partition the input parameters into parameters of elements
          describing this object (in our case theta -> theta1, theta2
  apply - evaluate the action (e.g., forward propagation, filtering, ..)
          in many cases this involves calling 'apply' for other objects
          (e.g., for different layers, kernels,...)
  Jthetamv  - compute the action of the Jacobian w.r.t theta on a vector
  JthetaTmv - compute the action of the transpose(Jacobian) w.r.t theta on a vector
  JYmv      - compute the action of the Jacobian w.r.t Y on a vector
  JYTmv     - compute the action of the transpose(Jacobian) w.r.t Y on a vector

 In addition, elements of this class also need to provide the folowing
 methods

  nTheta   - return the number of parameters, numel(theta) for this
              element (may have to ask lower-level elements for this)
  nFeatIn  - number of input features
  nFeatOut - number of output features
  initTheta - initialize parameters

"""
abstract type AbstractMeganetElement{T <: Number} end

"""
return number of trainable weights
"""
nTheta(this::AbstractMeganetElement) = error("An AbstractMeganetElement must provide method nTheta")

"""
returns the precision of the weights and features
"""
getPrecision{T}(this::AbstractMeganetElement{T}) = T

"""
number of input features
"""
nFeatIn(this::AbstractMeganetElement) = error("An AbstractMeganetElement must provide method nFeatIn")

"""
number of data going to classifier
"""
nDataOut(this::AbstractMeganetElement) = nFeatOut(this)

"""
number of output features
"""
nFeatOut(this::AbstractMeganetElement) = error("An AbstractMeganetElement must provide method nFeatOut")

"""
split weights into parts
"""
splitWeights(this::AbstractMeganetElement,theta::Array) = error("An AbstractMeganetElement must provide method splitWeights")

"""
initialize weights
"""
initTheta(this::AbstractMeganetElement) = randn(nTheta(this))

"""
transforms feature vector Y using the weights theta
"""
apply(this::AbstractMeganetElement,theta::Array,Y::Array,doDerivative) = error("An AbstractMeganetElement must provide method apply ")

"""
computes dY = transpose(J_Y(theta,Y))\*W

Input:

   W     - vector or matrix
   theta - current theta
   Y     - current Y
   tmp   - intermediates used in derivative computations
           (e.g., hidden features, activations, derivatives,
           ... )

Output:

  dY     - directional derivative, numel(dY)==numel(Y)
"""
function JYTmv(this::AbstractMeganetElement{T},Wdata::Array{T},W::Array{T},theta::Array{T},Y::Array{T},tmp=nothing) where {T <: Number}
    return JTmv(this,Wdata,W,theta,Y,tmp)[2]
end

"""
computes dZ = J_Y(theta,Y)\*dY

Input:

   dY    - perturbation in Y
   theta - current theta
   Y     - current Y
   tmp   - intermediates used in derivative computations
           (e.g., hidden features, activations, derivatives,
           ... )

Output:

  dZ     - directional derivative, numel(dZ)==numel(Z)
"""
function JYmv(this::AbstractMeganetElement{T},dY::Array{T},theta::Array{T},Y::Array{T},tmp=nothing) where {T <: Number}
    return Jmv(this,(T)[],dY,theta,Y,tmp);
end

"""
linearization with respect to Y, i.e.,

Z(theta,Y+dY) approx Z(theta,Y) + J\*dY

Input:

   theta - current theta
   Y     - current Y
   tmp   - intermediates used in derivative computations
           (e.g., hidden features, activations, derivatives,
           ... )

Output:

  Z     - current output features
  J     - Jacobian, LinearOperator
"""
function linearizeY(this::AbstractMeganetElement{T},theta::Array{T},Y::Array{T}) where {T <: Number}
    Z,~,tmp  = apply(this,theta,Y)
    J        = getJYOp(this,theta,Y,tmp)
    return Z,J
end

"""
J = getJYOp(this,theta,Y,tmp)

constructs Jacobian w.r.t. Y around current (theta,Y)

Input:

   W     - vector or matrix
   theta - current theta
   Y     - current Y
   tmp   - intermediates used in derivative computations
           (e.g., hidden features, activations, derivatives,
           ... )

Output:

  dY     - directional derivative, numel(dY)==numel(Y)
"""
function getJYOp(this::AbstractMeganetElement{T},theta::Array{T},Y::Array{T},tmp=nothing) where {T <: Number}
    nex    = div(length(Y),nFeatIn(this))
    m      = nex*nDataOut(this)
    n      = length(Y);
    Amv    = x -> JYmv(this,x,theta,Y,tmp)
    ATmv   = x -> JYTmv(this,x,(T)[],theta,Y,tmp)
    return LinearOperator(m,n,false,false,Amv,ATmv,ATmv);
end

"""
dZ = abstractMeganetElement.Jthetamv(this,W,theta,Y,tmp)

computes dZ = J_theta(theta,Y)\*dtheta

Input:

   dtheta- perturbation in theta
   theta - current theta
   Y     - current Y
   tmp   - intermediates used in derivative computations
           (e.g., hidden features, activations, derivatives,
           ... )

Output:

  dZ     - directional derivative, numel(dZ)==numel(Z)
"""
function Jthetamv(this::AbstractMeganetElement{T},dtheta::Array{T},theta::Array{T},Y::Array{T},tmp) where {T <: Number}
    return Jmv(this,dtheta,0*Y,theta,Y,tmp);
end

"""
computes dtheta = transpose(J_theta(theta,Y))\*W

Input:

   W     - vector or matrix
   theta - current theta
   Y     - current Y
   tmp   - intermediates used in derivative computations
           (e.g., hidden features, activations, derivatives,
           ... )

Output:

  dtheta - directional derivative, numel(dtheta)==numel(theta)
"""
function JthetaTmv(this::AbstractMeganetElement{T},Wdata::Array{T},W::Array{T},theta::Array{T},Y::Array{T},tmp=nothing) where {T <: Number}
     # error("There is a never ending recursive call over here... Need to fix this");
	 return JTmv(this,Wdata,W,theta,Y,tmp)[1]
end

"""
J = abstractMeganetElement.getJthetaOp(this,theta,Y,tmp)

constructs Jacobian w.r.t. theta as LinearOperator

Input:

   theta - current theta
   Y     - current Y
   tmp   - intermediates used in derivative computations
           (e.g., hidden features, activations, derivatives,
           ... )

Output:

  J     - Jacobian, LinearOperator
"""
function getJthetaOp(this::AbstractMeganetElement{T},theta::Array{T},Y::Array{T},tmp=nothing) where {T <: Number}
    nex    = div(length(Y),nFeatIn(this))
    m      = nex*nDataOut(this)
    n      = length(theta)
    Amv    = x -> Jthetamv(this,x,theta,Y,tmp)
    ATmv   = x -> JthetaTmv(this,x,(T)[],theta,Y,tmp)
    return LinearOperator(m,n,false,false,Amv,ATmv,ATmv)
end

"""
linearization with respect to theta, i.e.,

Z(theta+dth,Y) \approx Z(theta,Y) + J*dth

Input:

   theta - current theta
   Y     - current Y

Output:

  Z     - output features
  J     - Jacobian, LinearOperator
"""
function linearizeTheta(this::AbstractMeganetElement{T},theta::Array{T},Y::Array{T}) where {T <: Number}
    Z,Z2,tmp = apply(this,theta,Y)
    J       = getJthetaOp(this,theta,Y,tmp)
    return Z,J
end

"""
computes dZ = J_theta(theta,Y)\*dtheta + J_Y(theta,Y)\*dY

Input:

   dtheta- perturbation in theta
   dY    - perturbation of input features Y
   theta - current theta
   Y     - current Y
   tmp   - intermediates used in derivative computations
           (e.g., hidden features, activations, derivatives,
           ... )

Output:

  dZ     - directional derivative, numel(dZ)==numel(Z)
"""
function Jmv(this::AbstractMeganetElement{T},dtheta::Array{T},dY::Array{T},theta::Array{T},Y::Array{T},tmp=nothing) where {T <: Number}

    if isempty(dtheta) || norm(vec(dtheta))==0
        dZdata = (T)[]
        dZ     = (T)[]
    else
        dZdata,dZ = Jthetamv(this,dtheta,theta,Y,tmp)
    end

    if !isempty(dY) && norm(vec(dY))>0
        dZdt,dZt = JYmv(this,dY,theta,Y,tmp);
		if !isempty(dZdata)
        	dZdata +=  dZdt;
		else
			dZdata  = dZdt;
		end
		if !isempty(dZ)
        	dZ     +=  dZt;
		else
			dZ = dZt;
		end
    end
    return dZdata, dZ
end

"""
computes [dtheta;dY] = [J_theta(theta,Y)'; J_Y(theta)']\*Z

Input:

   Z     - perturbation of output
   theta - current theta
   Y     - current Y
   tmp   - intermediates used in derivative computations
           (e.g., hidden features, activations, derivatives,
           ... )
   doDerivative - vector with two elements for theta and Y
   derivative. Important only when nargout==1. default=[1,0];

"""
function JTmv(this::AbstractMeganetElement{T},Wdata::Array{T},W::Array{T},theta::Array{T},Y::Array{T},tmp=nothing) where {T <: Number}
    dtheta = JthetaTmv(this,Wdata,W,theta,Y,tmp)
    dY     = JYTmv(this,Wdata,W,theta,Y,tmp)
    return vec(dtheta), vec(dY)
end

"""
constructs Jacobian J(theta,Y) such that

Z(theta+dth,Y+dY) \approx Z(theta,Y) + J*[dth; dY]

Input:

   theta - current theta
   Y     - current Y
   tmp   - intermediates used in derivative computations
           (e.g., hidden features, activations, derivatives,
           ... )

Output:

  J     - Jacobian, LinearOperator
"""
function getJOp(this::AbstractMeganetElement{T},theta::Array{T},Y::Array{T},tmp=nothing) where {T <: Number}
    nex    = div(length(Y),nFeatIn(this))
    m      = nex*nDataOut(this)
    nth    = length(theta)
    nY     = length(Y)
    Amv    = x -> Jmv(this,x[1:nth],x[nth+1:end],theta,Y,tmp)
    ATmv   = x -> JTmv(this,x,(T)[],theta,Y,tmp)
    J      = LinearOperator(m,nth+nY,false,false,Amv,ATmv,ATmv)
    return J
end


text:
using BitSAD
using Test
using Statistics: mean
using LinearAlgebra: norm

include("bitstream.jl")
include("sbitstream.jl")
# include("dbitstream.jl")


text:
using Base.Broadcast: broadcastable, Broadcasted
using Swizzles.ValArrays
using Swizzles.ArrayifiedArrays
using Swizzles

foo(x) = x

@testset "ValArrays" begin
    A0 = ValArray(0)
    A1 = ValArray(1)

    @test A1.+A1 == 2

    @test A1.*A1 == 1

    B = rand(3,3)
    B .*= A0
    @test B == zeros(3,3)

    B = rand(3,3)
    B .= A0
    @test B == zeros(3,3)

    @test lift_vals(1) == A1
    @test typeof(lift_vals(1)) == typeof(A1)

    B = reshape(1:27, 3, 3, 3)
    C = reshape(1:9, 3, 3)

    @test Sum(2, 3).(B .* (A1 .+ C .+ A0) .+ A1) == Sum(2, 3).(B .* (1 .+ C .+ 0) .+ 1)
    @test typeof(Sum(2, 3).(B .* (A1 .+ C .+ A0))) == typeof(Sum(2, 3).(B .* (1 .+ C .+ 0)))

    @test typeof(Delay().(Sum(2, 3).(A0, B .* (A1 .+ C .+ A0) .+ A1))) ==
        typeof(lift_vals(Delay().(Sum(2, 3).(B .* (1 .+ C .+ 0) .+ 1))))

end


text:
@testset "Mass" begin

  #
  # Examples from arguments checks: indefinite and rank 0 cases
  #

  Qx, x = PolynomialRing(FlintQQ, "x")
  f = x^2 - 2
  K, a = NumberField(f, "a", cached = false)
  Kt, t = PolynomialRing(K, "t")
  g = t^2 - a*t + 1
  E, b = NumberField(g, "b", cached = false)
  D = matrix(E, 3, 3, [-1, 0, 0, 0, 1, 0, 0, 0, 1])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [2, 0, 0]), map(E, [a, 0, 0]), map(E, [b + 1, 0, 0]), map(E, [a*b + a, 0, 0]), map(E, [0, 10, 0]), map(E, [0, 10*a, 0]), map(E, [0, 2*b + 6*a + 10, 0]), map(E, [0, a*b + 5*a + 6, 0]), map(E, [0, 5, 5]), map(E, [0, 5*a, 5*a]), map(E, [0, b + 3*a, b + 3*a]), map(E, [0, a*b + 6, a*b + 6])]
  L = hermitian_lattice(E, gens, gram = D)
  @test !isdefinite(L)
  @test_throws ArgumentError mass(L)

  D = matrix(E, 0, 0, [])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[]
  L = hermitian_lattice(E, gens, gram = D)
  @test isdefinite(L)
  @test rank(L) == 0
  @test @inferred mass(L) == 1

  #
  # Bunch of examples
  #

  K, a = rationals_as_number_field()
  Kt, t = PolynomialRing(K, "t")
  g = t^2 + 5
  E, b = NumberField(g, "b", cached = false)
  D = matrix(E, 3, 3, [1, 0, 0, 0, 1, 0, 0, 0, 2])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [b + 8, b + 9, 0]), map(E, [-25*b + 66, -51//2*b + 171//2, -5//2*b + 1]), map(E, [104*b + 150, 132*b + 145, 5//2*b + 35//2]), map(E, [529*b - 47, 1243//2*b - 437//2, 28*b + 95//2])]
  L = hermitian_lattice(E, gens, gram = D)
  @test isdefinite(L)
  m = @inferred mass(L)
  @test m == sum([inv(fmpq(automorphism_group_order(LL))) for LL in genus_representatives(L)])

  K, a = rationals_as_number_field()
  Kt, t = PolynomialRing(K, "t")
  g = t^2 + 2
  E, b = NumberField(g, "b", cached = false)
  D = matrix(E, 4, 4, [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [-11*b - 18, 5*b - 8, 0, 0]), map(E, [-763//4*b - 373, 471//4*b - 279//2, -3//2*b + 5, 0]), map(E, [-205589//4*b - 45927, 54261//4*b - 71509//2, -13*b + 1018, -1//2*b + 2]), map(E, [-267023*b + 283211, -206429//2*b - 168273, 14671//4*b + 2083, 25//4*b + 15//2])]
  L = hermitian_lattice(E, gens, gram = D)
  m = @inferred mass(L)
  @test m == 1//512

  D = matrix(E, 4, 4, [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [14*b - 18, 8*b - 10, 0, 0]), map(E, [85//2*b - 493//2, 25*b - 277//2, 1, 0]), map(E, [-399//4*b - 2843//2, -205//4*b - 803, 11//4*b + 8, -1//4*b + 3//2]), map(E, [13929//2*b - 6826, 3955*b - 3803, -121//4*b + 61, -35//4*b + 7//2])]
  L = hermitian_lattice(E, gens, gram = D)
  m = @inferred mass(L)
  @test m == 1//512

  D = matrix(E, 4, 4, [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [-9*b - 38, 3*b + 8, 0, 0]), map(E, [-1//2*b - 5//2, -1//2, 1, 0]), map(E, [-716*b - 173, 329//2*b, -1//4*b + 1//2, -1//4*b]), map(E, [59*b - 84, -11*b + 22, 0, 0])]
  L = hermitian_lattice(E, gens, gram = D)
  m = @inferred mass(L)
  @test m == sum([inv(fmpq(automorphism_group_order(LL))) for LL in genus_representatives(L)])

  D = matrix(E, 4, 4, [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [-6*b - 38, 4*b + 10, 0, 0]), map(E, [1321//4*b - 2323//2, -91//4*b + 369, -1//2*b, 0]), map(E, [-2089//2*b + 5256, -23*b - 1610, 5//4*b - 5//2, 5//4*b]), map(E, [6376*b + 13976, -2645*b - 3174, -3//2*b - 9, 5//2*b - 4])]
  L = hermitian_lattice(E, gens, gram = D)
  m = @inferred mass(L)
  @test m == sum([inv(fmpq(automorphism_group_order(LL))) for LL in genus_representatives(L)])

  D = matrix(E, 3, 3, [1, 0, 0, 0, 1, 0, 0, 0, 1])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [-2*b-2, b+6, 0]), map(E, [0, 1, 1]), map(E, [b-6, -6*b+6, 0])]
  L = hermitian_lattice(E, gens, gram = D)
  m = @inferred mass(L)
  @test m == 1//32

  g = t^2 + 1
  E, b = NumberField(g, "b", cached = false)
  D = matrix(E, 4, 4, [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [8*b + 14, 4*b + 10, 0, 0]), map(E, [57*b + 443//2, 17*b + 303//2, -3//2*b + 3//2, 0]), map(E, [5//2*b - 7//2, b - 2, -3//2*b + 1, 3//2*b]), map(E, [1//2*b - 3, -3//2, -1//2*b + 1, b - 1//2])]
  L = hermitian_lattice(E, gens, gram = D)
  m = @inferred mass(L)
  @test m == sum([inv(fmpq(automorphism_group_order(LL))) for LL in genus_representatives(L)])

  D = matrix(E, 4, 4, [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 3])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [1, 1, 0, 0]), map(E, [-3*b + 5, 0, 8, 0]), map(E, [-40*b - 24, 0, -129//2*b + 2, -3//2*b + 1]), map(E, [-312*b - 208, 0, -1037//2*b - 17//2, -25//2*b + 15//2])]
  L = hermitian_lattice(E, gens, gram = D)
  m = @inferred mass(L)
  @test m == sum([inv(fmpq(automorphism_group_order(LL))) for LL in genus_representatives(L)])

  D = matrix(E, 4, 4, [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 3])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [-b - 2, 1, 0, 0]), map(E, [-9*b - 9, 0, -8, 0]), map(E, [-30*b - 24, 0, -9//2*b - 31//2, -1//2*b - 3//2]), map(E, [-48*b - 96, 0, 13//2*b - 89//2, 1//2*b - 9//2])]
  L = hermitian_lattice(E, gens, gram = D)
  @test isdefinite(L)
  m = @inferred mass(L)
  @test m == sum([inv(fmpq(automorphism_group_order(LL))) for LL in genus_representatives(L)])

  D = matrix(E, 6, 6, [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [-11*b + 1, -3*b + 1, 0, 0, 0, 0]), map(E, [-42*b + 27//2, -4*b + 17//2, 3//2*b + 13//2, 0, 0, 0]), map(E, [-2074*b + 4314, 370*b + 952, 624*b + 368, 7*b - 1, 0, 0]), map(E, [-1//2, -1//2, 0, 1//2*b - 1//2, 1, 0]), map(E, [-64401//2*b - 379225, -127585//2*b - 50257, -57600*b, -440*b + 352, -1//2*b - 1, -1//2*b - 1]), map(E, [178435*b + 414553, 90613*b + 32617, 64800*b - 21600, 363*b - 561, b + 1, b + 1])]
  L = hermitian_lattice(E, gens, gram = D)
  m = @inferred mass(L)
  @test m == 1//147456

  Qx, x = PolynomialRing(FlintQQ, "x")
  f = x^2-3
  K, a = NumberField(f, "a", cached = false)
  Kt, t = PolynomialRing(K, "t")
  g = t^2+1
  E, b = NumberField(g, "b", cached = false)
  D = matrix(E, 3, 3, [1, 0, 0, 0, 1, 0, 0, 0, 1])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [(-1//2*a-1//2)*b+(1//2*a+5//2), -1, 0]), map(E, [(a-9)*b+ 15*a+3, 0, -a*b+ 8*a+13]), map(E, [(-1610*a+2990)*b + 2070*a-2990, 0, (-68*a+393)*b + 19*a+462])]
  L = hermitian_lattice(E, gens, gram = D)
  m = @inferred mass(L)
  @test m == 1//108

  Qx, x = PolynomialRing(FlintQQ, "x")
  f = x^4-x^3-4*x^2+4*x+1
  K, a = NumberField(f, "a", cached = false)
  Kt, t = PolynomialRing(K, "t")
  g = t^2+(a^3 - 1*a^2 - 4*a + 5)
  E, b = NumberField(g, "b", cached = false)
  D = matrix(E, 3, 3, [1, 0, 0, 0, 1, 0, 0, 0, 1])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [1, 0, 0]), map(E, [0, 1, 0]), map(E, [0, 0, 1])]
  L = hermitian_lattice(E, gens, gram = D)
  m = @inferred mass(L)
  @test m == 1//10125

end


text:
using Combinatorics
using Test

include("numbers.jl")
include("factorials.jl")
include("combinations.jl")
include("permutations.jl")
include("partitions.jl")
include("multinomials.jl")
include("youngdiagrams.jl")


text:
# see results in README.md

using BenchmarkTools
using PEG # this package
using StringParserPEG # another similar package
using JSON # hand-written JSON parser for comparison

# PEG grammar
commasPEG(head, tail) = [head, (t[2] for t ∈ tail)...]
bracketsPEG(left, middle, right) = (isempty(middle) ? [] : middle[1])
@rule object =
  r"{"p & (pair & (r","p & pair)[*] > commasPEG)[:?] & r"}"p > Dict∘bracketsPEG
@rule pair = string & r":"p & value > (k,c,v) -> (k => v)
@rule array =
  r"\["p & (value & (r","p & value)[*] > commasPEG)[:?] & r"]"p > bracketsPEG
@rule value =
  string, number, object, array,
  r"true|false"w |> Meta.parse, r"null"w |> x->nothing
@rule string = r"\"([^\"\\]|\\[\"\\/bfnrt]|\\u[0-9A-Fa-f]{4})*\""p |> s->Meta.parse(replace(s, r"\$"=>"\\\$"))
@rule number = r"-?(0|[1-9]\d*)(\.\d+)?([Ee][+-]?\d+)?"p |> Meta.parse
@rule startPEG = r"\s*" & value > (w,v)->v

# StringParserPEG grammar
commasSPP(r,v,f,l,c) = (length(c) == 1 ? [c[1]] : [c[1], (t.children[2] for t ∈ c[2].children)...])
bracketsSPP(r,v,f,l,c) = (isempty(c[2].children) ? [] : c[2].children[1])
braces = Dict∘bracketsSPP
parsevalue(r,v,f,l,c) = Meta.parse(v)
parsestr(r,v,f,l,c) = Meta.parse(replace(v, r"\$"=>"\\\$"))
grammarSPP = Grammar("""
  object => (r({\\s*)r & ?((pair & *(r(,\\s*)r & pair)) {Main.commasSPP}) & r(}\\s*)r) {Main.braces}
  pair => (string & r(:\\s*)r & value) {(r,v,f,l,c)->(c[1] => c[3])}
  array => (r(\\[\\s*)r & ?((value & *(r(,\\s*)r & value)) {Main.commasSPP}) & r(]\\s*)r) {Main.bracketsSPP}
  value => string | number | object | array | r((true|false)\\b\\s*)r {Main.parsevalue} | r(null\\b\\s*)r {(r,v,f,l,c)->:null}
  string => r(\\"([^\\"\\\\]|\\\\[\\"\\\\/bfnrt]|\\\\u[0-9A-Fa-f]{4})*\\"\\s*)r {Main.parsestr}
  number => r(-?(0|[1-9]\\d*)(\\.\\d+)?([Ee][+-]?\\d+)?\\s*)r {Main.parsevalue}
  start => (r(\\s*)r & value) {(r,v,f,l,c)->c[1]}
""")

# load JSON string from file from
# https://github.com/Chevrotain/chevrotain/blob/gh-pages/performance/samples/10K_json.js
json_str = open(x->read(x,String), "10K.json")

println("JSON:")
display(@benchmark JSON.parse(json_str) seconds=60)
println("PEG:")
display(@benchmark parse_whole(startPEG, json_str) seconds=60)
println("StringParserPEG:")
display(@benchmark parse(grammarSPP, json_str) seconds=60)


text:
using Dates
using Unitful: unit, ustrip

import Base: *, convert, round
import Statistics: mean

# Handle working with Missings
const UnitfulMissing = Union{<:Unitful.Quantity, Missing}

# Helper functions
round(x::AbstractArray{<:Quantity}, r::Int) = round(ustrip(x), r) * unit(eltype(x))
round(x::T, r::Int) where {T<:Quantity} = round(ustrip(x), r) * unit(x)

mean(x::AbstractArray{<:Quantity}, r::Int) = mean(ustrip(x), r) * unit(eltype(x))

"""
    asqtype(x::T) where {T<:Unitful.Units} -> Type

A helper function to convert from the "dimensions" of a Unitful quantity to the "quantities",
as they are treated separately.
"""
asqtype(x::T) where {T<:Unitful.Units} = typeof(1.0*x)

"""
    fustrip(x::Array{T}) where {T<:Any} = Array

Operation to strip the units an Parametric Array{T} Type. Needed for operating on DataFrames?
"""
fustrip(x::Array{T}) where {T<:Any} = map(t -> ustrip(t), x)

# Handle working with `Period`s
*(x::Unitful.Units, y::Period) = *(y, x)
*(x::Period, y::Unitful.Units) = convert(y, x)
function convert(a::Unitful.Units, x::Period)
    sec = Dates.value(Dates.Second(x))
    uconvert(a, (sec)u"s")
end

# Methods to drop
# Exist to test that (offsets)u"hr" should work the same way
dt2umin(t::AbstractArray{Dates.Minute}) = Dates.value.(t).*u"minute"

max_stars_repo_name: invenia/PowerSystemsUnits

max_issues_repo_name: invenia/PowerSystemsUnits

max_forks_repo_name: invenia/PowerSystemsUnits


text:
using DataFrames
df = DataFrame(a = [1,2], b = [2.3, 4.5])

using PlotlyJS, VegaLite, ElectronDisplay

electrondisplay(df |> @vlplot(:point, x=:a, y=:b))
electrondisplay(plot(scatter(df, x=:a, y=:b, mode="markers")))

if !Sys.iswindows()
    @info "compiling Blink, TableView"
    using Blink, TableView
    w = Window()
    body!(w, showtable(df))
else
    @warn "compiling Blink does not work on Windows currently, only including TableView"
    using TableView
    t = showtable(df)
end

using Query
q = @from i in df begin
    @where i.a > 1
    @select {i.a, b=i.b^2}
    @collect
    end


text:
export User

"""
A Discord user.
More details [here](https://discordapp.com/developers/docs/resources/user#user-object).
"""
struct User
    id::Snowflake
    # The User inside of a Presence only needs its ID set.
    username::Optional{String}
    discriminator::Optional{String}
    avatar::OptionalNullable{String}
    bot::Optional{Bool}
    mfa_enabled::Optional{Bool}
    locale::Optional{String}
    verified::Optional{Bool}
    email::OptionalNullable{String}  # Not supposed to be nullable.
end
@boilerplate User :constructors :docs :lower :merge :mock


text:
using LinearAlgebra

import Base: size, inv, getindex, setindex, copy, display
import LinearAlgebra: ldiv!, det, logdet, isposdef, diag

#####
##### Incremental Cholesky Decomposition
#####

@doc raw"""
"""
mutable struct IncrementalCholesky{T} <: Factorization{T}
    n::Int
    info::Int
    R::Matrix{T}
end

function incremental_cholesky!(A::AbstractMatrix, n; check=true)
    size(A, 1) == size(A, 2) || throw(ErrorException("A is not square"))
    @views chol = cholesky!(Symmetric(A[1:n, 1:n]); check=check)
    return IncrementalCholesky(n, chol.info, A)
end

function incremental_cholesky(T, capacity)
    return incremental_cholesky!(Array{T}(undef, capacity, capacity), 0)
end

valid_update(R::IncrementalCholesky, k) = 0 <= R.n + k <= size(R.R, 1)

# extend is deprecated in favor of view_next
extend(R::IncrementalCholesky, k) = view_next(R, k)
function view_next(R::IncrementalCholesky, k)
    valid_update(R, k) || throw(ErrorException("Size cannot exceed capacity"))
    return @views R.R[1:R.n, R.n+1:R.n+k], R.R[R.n+1:R.n+k, R.n+1:R.n+k]
end

# update! is depracated in favor of compute_next!
update!(R::IncrementalCholesky, k) = compute_next!(R, k)
function compute_next!(R::IncrementalCholesky, k; check=true)
    valid_update(R, k) || throw(ErrorException("Size cannot exceed capacity"))
    @views A12, A2 = R.R[1:R.n, R.n+1:R.n+k], R.R[R.n+1:R.n+k, R.n+1:R.n+k]

    ldiv!(get_chol(R).L, A12)
    mul!(A2, A12', A12, -1, 1)

    chol = cholesky!(Symmetric(A2); check=check)
    R.info = chol.info
    R.n += k
    return nothing
end

# remove! is depracated in favor of remove_last!
remove!(R::IncrementalCholesky, k) = remove_last!(R, k)
function remove_last!(R::IncrementalCholesky, k)
    valid_update(R, -k) || throw(ErrorException("Size cannot be negative"))
    R.n -= k
    return nothing
end

cap(R::IncrementalCholesky) = size(R.R, 1)

#error handling slow in julia
get_chol(R::IncrementalCholesky) = Cholesky(view(R.R, 1:R.n, 1:R.n), :U, R.info) #info indicates if matrix is still positive definite

size(R::IncrementalCholesky, args...) = size(get_chol(R), args...)

ldiv!(R::IncrementalCholesky, args...) = ldiv!(get_chol(R), args...)

inv(R::IncrementalCholesky) = inv(get_chol(R))

det(R::IncrementalCholesky) = det(get_chol(R))

logdet(R::IncrementalCholesky) = logdet(get_chol(R))

isposdef(R::IncrementalCholesky) = R.info == 0

#lastindex(R::IncrementalCholesky{Float64}, t::Int64) = size(R)[1] #for doing R[2:end], for example

#####
##### Growable Data Array
#####

mutable struct DataArray{T,N} <: AbstractArray{T,N}
    n::Int
    A::Array{T,N}
end

data_array(A::Array, n) = DataArray(n, A)
data_array(T, capacity, dims...) = DataArray(0, Array{T}(undef, dims..., capacity))

function array_view(A::DataArray)
    return view(A.A, ntuple(_ -> Colon(), ndims(A) - 1)..., 1:A.n)
end

cap(A::DataArray) = size(A.A)[end]

size(A::DataArray, args...) = size(array_view(A), args...)

getindex(A::DataArray, args...) = getindex(array_view(A), args...)

setindex!(A::DataArray, args...) = setindex!(array_view(A), args...)

valid_update(A::DataArray, k) = 0 <= A.n + k <= size(A.A)[end]

# extend is deprecated in favor of view_next
extend(A::DataArray, k) = view_next(A, k)
function view_next(A::DataArray, k)
    valid_update(A, k) || throw(ErrorException("Size cannot exceed capacity"))
    return view(A.A, :, A.n+1:A.n+k)
end

# update! is depracated in favor of compute_next!
update!(A::DataArray, k) = compute_next!(A, k)
function compute_next!(A::DataArray, k)
    valid_update(A, k) || throw(ErrorException("Size cannot exceed capacity"))
    A.n += k
    return nothing
end

# remove! is depracated in favor of remove_last!
remove!(A::DataArray, k) = remove_last!(A, k)
function remove_last!(A::DataArray, k)
    valid_update(A, -k) || throw(ErrorException("Size cannot be negative"))
    A.n -= k
    return nothing
end

const DataMatrix{T} = DataArray{T, 2}
const DataVector{T} = DataArray{T, 1}

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
# LP programming using MATLAB-like construction
#
#   maximize 1000 x + 350 y
#
#       s.t. x >= 30, y >= 0
#            x - 1.5y >= 0  (i.e. -x + 1.5 y <= 0)
#            12 x + 8 y <= 1000
#            1000 x + 300 y <= 70000
#
#   solution: (59.0909, 36.3636)
#   objv: 71818.1818
#

using Gurobi, Test

@testset "LP 02" begin

	env = Gurobi.Env()
	setparam!(env, "OutputFlag", 0)

	model = gurobi_model(env;
		name="lp_02",
		sense=:maximize,
		f = [1000., 350.],
		A = [-1. 1.5; 12. 8.; 1000. 300.],
		b = [0., 1000., 70000.],
		lb = [0., 30.])

	optimize(model)

	@test isapprox(get_solution(model), [59.0909, 36.3636], atol=1e-4)
	@test isapprox(get_objval(model), 71818.1818, atol=1e-4)
end


text:
# mean square loss: f(x,y)=(sum_{i=1:I} (x_i-y_i)^2)/I
# One way to do this is meansquare( x-y). However, this would define a new now x-y. This would be fast, but require storing x-y.
# As an alternative, we'll recompute here x-y, rather than storing it. This will be slower, but saves on storage.

FmeanSquareLoss(x,y)=(FmeanSquare(axpy(-1.0,y,x))[1],nothing)
function FmeanSquareLoss_inplace(value,auxvalue,x,y)
        FmeanSquare_inplace(value,auxvalue,axpy(-1.0,y,x))
end

function DmeanSquareLoss(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::Array,y::Array)
    if derivativeIDX==1
        tmp=copy(grad_c)  # how to avoid this tmp?
        scale!(2./length(x),tmp)
        axpy!(tmp[1],axpy(-1.0,y,x),grad_n)
    elseif derivativeIDX==2
        tmp=copy(grad_c)
        scale!(2./length(x),tmp)
        axpy!(tmp[1],axpy(-1.0,x,y),grad_n)
    end
end

if PROC=="GPU"

    function FmeanSquareLoss(x::CudaArray,y::CudaArray)
        tmp=copy(x)
        axpy!(-1.0,y,tmp)
        return (FmeanSquare(tmp)[1],nothing)
    end

    function FmeanSquareLoss_inplace(value,auxvalue,x::CudaArray,y::CudaArray)
        #copy!(value,x)
        #axpy!(-1.0,y,value);FmeanSquare_inplace(value,auxvalue,value)
        # why doesn't this work??

        tmp=copy(x)
        axpy!(-1.0,y,tmp);FmeanSquare_inplace(value,auxvalue,tmp)
        free(tmp)
    end

    function DmeanSquareLoss(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::CudaArray,y::CudaArray)
        if derivativeIDX==1
            DmeanSquareLoss!(grad_c,x,y,grad_n)
        elseif derivativeIDX==2
            DmeanSquareLoss!(grad_c,y,x,grad_n)
        end
    end
end

Derivative[FmeanSquareLoss]=DmeanSquareLoss
Inplace[FmeanSquareLoss]=FmeanSquareLoss_inplace

meanSquareLoss(nx::ADnode,ny::ADnode)=ADnode(FmeanSquareLoss,[nx ny])
export  meanSquareLoss

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
using .Default_params

function create_or_open_group(file::HDF5.File, group_name::String, subgroup_name::String = "None")

    if !haskey(file, group_name)
        group = create_group(file, group_name)
    else
        group = open_group(file, group_name)
    end

    if subgroup_name != "None"
        if !haskey(group, subgroup_name)
            subgroup = create_group(group, subgroup_name)
        else
            subgroup = open_group(group, subgroup_name)
        end
    else
        subgroup = nothing
    end

    return group, subgroup

end

function write_timers(lengths::Vector{Int}, size::Int, chars::AbstractVector{UInt8}, params::FilterParameters)

    write_timers(lengths, size, chars, params.output_filename)

end

function write_timers(lengths::Vector{Int}, size::Int, chars::AbstractVector{UInt8}, filename::String)

    group_name = "timer"

    h5open(filename, "cw") do file

        if !haskey(file, group_name)
            group = create_group(file, group_name)
        else
            group = open_group(file, group_name)
        end

        sum_lengths = cumsum(lengths)
        for i in 1:size
            timer_string = String(chars[1 + (i > 1 ? sum_lengths[i - 1] : 0):sum_lengths[i]])
            dataset_name = "rank" * string(i-1)

            if !haskey(group, dataset_name)
                group[dataset_name] = timer_string
            else
                @warn "Write failed, dataset " * group_name * "/" * dataset_name *  " already exists in " * file.filename * "!"
            end
        end
    end
end


text:
using Agents
using Agents.CA1D

# 0. Define the rule
# rules = Dict("111"=>"0", "110"=>"0", "101"=>"0", "100"=>"1", "011"=>"1", "010"=>"1", "001"=>"1", "000"=>"0")  # rule 30
rules = Dict("111"=>"0", "110"=>"0", "101"=>"0", "100"=>"1", "011"=>"0", "010"=>"1", "001"=>"1", "000"=>"0")  # rule 22

# 1. Build the model
model = CA1D.build_model(rules=rules, ncols=101)  # creates a model where all columns are "0"
model.agents[50].status="1"

# 2. Run the model, collect data, and visualize it
runs = 100
CA1D.ca_run(model, runs)


text:
macro run_mixedmode(expr)
    quote
        thunk()=$(esc(expr))
        Cassette.recurse(MixModeDebugCtx(), thunk)
    end
end


text:
function summary(data::ModelSelection.ModelSelectionData, filename::String; resultnum::Union{Nothing, Int}=nothing)
    return summary(data, filename=filename, resultnum=resultnum)
end

function summary(data::ModelSelection.ModelSelectionData; filename::Union{Nothing, String}=nothing, resultnum::Union{Nothing, Int}=nothing)
    addextras(data, :summary, filename, nothing)
    if size(data.results, 1) > 0
        if resultnum != nothing
            return summary(data, data.results[resultnum], filename=filename)
        else
            str = ""
            for i in 1:size(data.results, 1)
                str = string(str, summary(data, data.results[i]))
            end
            writefile(str, filename)
            return str
        end
    end
    return ""
end

function summary(data::ModelSelection.ModelSelectionData, result::ModelSelection.AllSubsetRegression.AllSubsetRegressionResult, filename::String)
    return summary(data, result, filename=filename)
end

function summary(data::ModelSelection.ModelSelectionData, result::ModelSelection.AllSubsetRegression.AllSubsetRegressionResult; filename::Union{Nothing, String}=nothing)
    str = ModelSelection.AllSubsetRegression.to_string(data, result)
    writefile(str, filename)
    return str
end

function summary(data::ModelSelection.ModelSelectionData, result::ModelSelection.CrossValidation.CrossValidationResult, filename::String)
    return summary(data, result, filename=filename)
end

function summary(data::ModelSelection.ModelSelectionData, result::ModelSelection.CrossValidation.CrossValidationResult; filename::Union{Nothing, String}=nothing)
    str = ModelSelection.CrossValidation.to_string(data, result)
    writefile(str, filename)
    return str
end

function writefile(str, filename=nothing)
    if filename != nothing
        file = open(filename, "w")
        write(file, sum)
        close(file)
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "open62541"
version = v"1.2.2"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/open62541/open62541.git",
              "ecf5a703785877a8719a0cda863a98455f7d5d12"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/open62541/
atomic_patch -p1 ../patches/0001-fix-core-Explicit-cast-to-avoid-compiler-warning-420.patch
atomic_patch -p1 ../patches/0002-refactor-pubsub-Fix-check-macros-and-slightly-clean-.patch
mkdir build && cd build/
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DUA_MULTITHREADING=100 \
    -DUA_ENABLE_SUBSCRIPTIONS=ON \
    -DUA_ENABLE_METHODCALLS=ON \
    -DUA_ENABLE_PARSING=ON \
    -DUA_ENABLE_NODEMANAGEMENT=ON \
    -DUA_ENABLE_AMALGAMATION=ON \
    -DUA_ENABLE_IMMUTABLE_NODES=ON \
    -DBUILD_SHARED_LIBS=ON \
    ..
make -j${nproc}
make install
install_license ../LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

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

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
#using LinearAlgebra
#using DelimitedFiles
#using Statistics
#using Plots
#using DelimitedFiles
#using Random
#using PyCall
#using SparseArrays

#plt = pyimport("matplotlib.pyplot")

#include("Geometry_twisted.jl")

### for triangular lattice, use pi-flux hopping!!!

function get_valley_operator_twisted_bilayer(R,r,mode,theta,r0,layer_factor=1)  # d: inter-layer distance, lambda: parameter, mode="uniform_hop" or "pi-flux"
    t0=now()
    N=size(R,1)

    I,J,value=ones(Int,N*13),ones(Int,N*13),zeros(Complex,N*13)
    n=1
    for i=1:N
        for j=1:N
            dR=R[i]-(R[j]+r)
            z=dR[3]
            L=norm(dR)

            if abs(z)<0.01 && L<2
                t=get_valley_operator_hopping(R[i],R[j]+r,mode,theta,r0,layer_factor)
                if abs(t)>0
                    I[n]=i
                    J[n]=j
                    value[n]=t
                    n+=1
                end
            end
        end
    end
    P=sparse(I,J,value,N,N)
    P=dropzeros(P)
    t1=now()
    println("time to construct P: ",t1-t0)
    return P
end

function get_valley_operator_hopping(R,R2,mode,theta,r0,layer_factor=1)
    factor=1
    U=rotation_along_z_3D(-(theta))
    if R[3]==1
        R=U*(R-r0)
        R2=U*(R2-r0)
        factor=layer_factor
    end
    x=R[1]-R2[1]
    y=R[2]-R2[2]
    if mode=="uniform_hop"    #valley operator for honeycomb
        hop=0
        if round(x/sqrt(3),digits=3)==1 && round(y,digits=3)==0
            hop=im
        elseif round(x/sqrt(3),digits=3)==-1 && round(y,digits=3)==0
            hop=-im
        elseif round(x*2/sqrt(3),digits=3)==-1 && round(y,digits=3)==3/2
            hop=im
        elseif round(x*2/sqrt(3),digits=3)==1 && round(y,digits=3)==-3/2
            hop=-im
        elseif round(x*2/sqrt(3),digits=3)==-1 && round(y,digits=3)==-3/2
            hop=im
        elseif round(x*2/sqrt(3),digits=3)==1 && round(y,digits=3)==3/2
            hop=-im
        end
        return hop*layer_factor/(3*sqrt(3))
    elseif mode=="pi_flux"
        hop=0
        if round(x,digits=3)==1 && round(y,digits=3)==0
            hop=-im
        elseif round(x,digits=3)==-1 && round(y,digits=3)==0
            hop=im
        elseif round(x,digits=3)==0 && round(y/sqrt(3),digits=3)==-1
            hop=-im
        elseif round(x,digits=3)==0 && round(y/sqrt(3),digits=3)==1
            hop=im
        end
        return hop*factor/4
    else
        println("invalid mode for valley operator")
    end
end

function get_valley_operator_inter_twisted_bilayer(g::geometry_twisted,mode,layer_factor=1)
    R=g.sites
    theta=g.twist_angle
    r0=g.inter_alignment

    P_intra=get_valley_operator_twisted_bilayer(R,[0,0,0],mode,theta,r0,layer_factor)
    P_tx=get_valley_operator_twisted_bilayer(R,g.inter_vector.x,mode,theta,r0,layer_factor)
    P_ty=get_valley_operator_twisted_bilayer(R,g.inter_vector.y,mode,theta,r0,layer_factor)
    P_txy=get_valley_operator_twisted_bilayer(R,g.inter_vector.xy,mode,theta,r0,layer_factor)
    P=get_hamiltonian(P_intra,P_tx,P_ty,P_txy,g)
    return P
end

function get_Pk(P::get_hamiltonian,k::Array{Float64,1})
    P0=exp(im*dot(k,P.geometry.inter_vector.x))*P.tx+exp(im*dot(k,P.geometry.inter_vector.y))*P.ty+exp(im*dot(k,P.geometry.inter_vector.xy))*P.txy
    Pk=P.intra+P0+adjoint(P0)
    return Pk
end

function get_Pk_ribbon(P::get_hamiltonian,k::Array{Float64,1})   # here k represents the phase difference
    @assert k[1]*k[2]==0 "please input 1D momenta"
    k[1]==0 ? P0=exp(im*k[2])*P.ty : P0=exp(im*k[1])*P.tx

    Pk=P.intra+P0+adjoint(P0)
    return Pk
end

function plot_P(P,R)
    plt.figure(figsize=(6,6),dpi=80)
    N=size(R,1)
    for i=1:N
        for j=i:N
            if imag(P[i,j])>0.1
                plt.plot((R[i][1],R[j][1]), (R[i][2],R[j][2]),color="red",linewidth=norm(P[i,j])*1)
            elseif imag(P[i,j])<-0.1
                plt.plot((R[i][1],R[j][1]), (R[i][2],R[j][2]),color="blue",linewidth=norm(P[i,j])*1)
            end
        end
    end
    plt.axis([-5,5,-5,5])
    plt.xlabel("")
    plt.ylabel("")
    plt.show()
end

function get_layer_operator_twisted_bilayer(R)  # d: inter-layer distance, lambda: parameter, mode="uniform_hop" or "pi-flux"
    N=size(R,1)

    L=spzeros(Complex,N,N)
    for i=1:N
        if R[i][3]==0
            L[i,i]=1
        else
            L[i,i]=-1
        end
    end
    return L
end


text:
const data_root_dir = joinpath(homedir(), "Data", "project_data")
const train_data_full = joinpath(data_root_dir, "dep-wiki-1based")
const wordlist = joinpath(data_root_dir, "map-words")
const rellist = joinpath(data_root_dir, "map-dep-rel")
const wordmaps_full = joinpath(data_root_dir, "maps-word.jld")
const train_data_reduced = joinpath(data_root_dir,
                             "pairs-wiki-1based-global_index-30k")

const project_root = joinpath(data_root_dir, "SpSkipGram")
const train_data = joinpath(project_root, "pairs-wiki-1based-20k.idx")
const wordmaps = joinpath(project_root, "maps-word.jld")
const init_model = joinpath(project_root, "rw2vec-model-20k.jld")

using ..SpWord2Vec: BACKEND
using ..SpWord2Vec: DEVICE
using ..SpWord2Vec: vec_lib


text:
# A simple relation that uses hash-based lookup
struct IsEqualRelation{names, T <: Tuple} <: Relation{names}
    values::T
end

Base.isequal(::Label{name}, value) where {name} = IsEqualRelation{(name,)}((value,))

Base.iteratorsize(::Type{<:IsEqualRelation}) = Base.HasLength()
Base.length(::IsEqualRelation) = 1
Base.start(::IsEqualRelation) = false
Base.next(r::IsEqualRelation{names}, i) = (Row{names}(r.values), true)
Base.done(::IsEqualRelation, i) = i
Base.first(r::IsEqualRelation) = Row{names}(r.values)

# Hash-based lookup table
struct UniqueHashIndex{names, T} <: UniqueIndex{names}
    hashmap::Dict{T, Int}
end

function UniqueHashIndex{hashnames}(t::Table{names}) where {hashnames, names}
    HashIndex{hashnames}(t::Table[labels(hashnames)...])
end

function UniqueHashIndex{names}(t::Table{names}) where {names}
    T = eltype(t, labels(names))
    hashmap = Dict{T, Int}()
    for i = 1:length(t)
        @inbounds row = t[i]
        hashmap[row] = i
    end
    return UniqueHashIndex{names, T}(hashmap)
end

function makeindex(t::Table, UniqueHashIndex{names})

end


text:
module TAPAS

using CSV
using DataFrames
using StatsBase
using Dates
using Printf
using Plots

struct Node
    v::Int64                                # Node value
    T::Vector{Int64}                        # Predecessor nodes (tail node values)
    H::Vector{Int64}                        # Successor nodes (head node values)
end
# NOTE: To avoid circular reference predecessor and sucessor arrays store node indices and not the node itself

struct Origin
    n::Node                                 # Origin node
    k::Int64                                # Origin node index in vector G.O; k = findfirst(x -> (x == o), O)
    S::Vector{Node}                         # Vector of destination nodes
    Q::Vector{Float64}                      # Vector of flows
    Lᵖ::Vector{Int64}                       # Predecessor label
end

mutable struct Arc
    t::Node                                 # Tail node (From)
    h::Node                                 # Head node (To)
    V::Float64                              # Volume capacity of the arc
    d::Float64                              # Length of the arc
    tₒ::Float64                             # Free flow travel time on the arc
    α::Float64                              # BPR parameter
    β::Float64                              # BPR parameter
    xʳ::Vector{Float64}                     # Origin based arc flow
    c::Float64                              # Arc cost
    c′::Float64                             # Arc cost derivative
    ϕ::Bool                                 # Assignment boolean (UE => ϕ = false; SO => ϕ = true)
end

struct PAS
    e₁::Vector{Arc}                         # First segment
    e₂::Vector{Arc}                         # Second segment
    o::Origin                               # Assosciated origin
end

struct Graph
    name::String                            # Network name
    N::Vector{Node}                         # Vector of nodes
    A::Vector{Arc}                          # Vector of arcs
    K::Dict{Tuple{Int64, Int64}, Int64}     # Collection of arc indices mapped to tuple of arc head and tail node value (one to one mapping - (i,j) => findfirst(x -> (x == (i,j)), [(a.t.v, a.h.v) for a in A]))
    O::Vector{Origin}                       # Vector of origins
end

function Base.:∈(pₒ::PAS, P::Vector{PAS})
    for p in P if pₒ.e₁ == p.e₁ && pₒ.e₂ == p.e₂ return true end end
    return false
end
Base.:∉(pₒ::PAS, P::Vector{PAS}) = !(Base.:∈(pₒ, P))

include("build.jl")
include("func.jl")
include("TA.jl")

export itapas

end


text:
###################################################
## complex virtual city, with main city and suburbs
###################################################

"""
Elaborate network, depending on parameters
- width = main city width
- nSub  = number of city's suburb (half the city width)
Represents whole city, with commuting effects
"""
function urbanNetwork(width::Int=8; distance::Float64=200.)
    subwidth = div(width,2)

    function coordToLoc(i::Int, j::Int, c::Int)
        if c ==0
            return j + (i-1)*width
        else
            return width^2 + (c-1)*subwidth^2 + j + (i-1)*subwidth
        end
    end
    function coordToLoc(k::Int, c::Int)
        if c ==0
            return k
        else
            return width^2 + (c-1)*subwidth^2 + k
        end
    end

    #square cities:
    urban = squareNetwork(width, distance=distance)
    suburbs = Network[squareNetwork(subwidth, distance=distance) for i in 1:8]
    resize!(urban.nodes, width^2 + 8*subwidth^2)

    # 4  3  2
    # 5  0  1
    # 6  7  8

    #Positionning the adjacent suburbs:
    for (i,a,b) in [(1,1,0), (3,0,1), (5,-1,0), (7,0,-1)]
        offset = rand()*width*0.55
        for (k,n) in enumerate(suburbs[i].nodes)
            #offset
            x =  n.x + a*distance*((width-1)/2 + (subwidth+1)/2 + offset)
            y =  n.y + b*distance*((width-1)/2 + (subwidth+1)/2 + offset)
            urban.nodes[coordToLoc(k,i)] =  Node(x,y)
        end
        for ((l,m),r) in suburbs[i].roads
            n1,n2 = coordToLoc(l,i), coordToLoc(m,i)
            urban.roads[n1,n2] = Road(n1,n2, r.distance, r.roadType)
        end
    end

    #Positionning the diagonal suburbs:
    for (i,a,b) in [(2,1,1), (4,-1,1), (6,-1, -1), (8,1,-1)]
        offset = rand()*width*0.3
        for (k,n) in enumerate(suburbs[i].nodes)
            #rotation and offset
            x =  (n.x - n.y)/sqrt(2) + a*distance*((width-1)/2 + offset + (subwidth+1)/(2*sqrt(2)))
            y =  (n.y + n.x)/sqrt(2) + b*distance*((width-1)/2 + offset + (subwidth+1)/(2*sqrt(2)))
            urban.nodes[coordToLoc(k,i)] =  Node(x,y)
        end
        for ((l,m),r) in suburbs[i].roads
            n1, n2 = coordToLoc(l,i), coordToLoc(m,i)
            urban.roads[n1,n2] = Road(n1, n2, r.distance, r.roadType)
        end
    end

    #Adding connecting roads
    connections = [ #city1, i1, j1, city2, i2, j2, type
        #perimeter
        (1, 1              , div(subwidth,2), 2, subwidth       , div(subwidth,2), 2),
        (2, 1              , div(subwidth,2), 3, div(subwidth,2), subwidth       , 2),
        (3, div(subwidth,2), 1              , 4, div(subwidth,2), subwidth       , 2),
        (4, div(subwidth,2), 1              , 5, 1              , div(subwidth,2), 2),
        (5, subwidth       , div(subwidth,2), 6, 1              , div(subwidth,2), 2),
        (6, subwidth       , div(subwidth,2), 7, div(subwidth,2), 1              , 2),
        (7, div(subwidth,2), subwidth       , 8, div(subwidth,2), 1              , 2),
        (8, div(subwidth,2), subwidth       , 1, subwidth       , div(subwidth,2), 2),
        # link to city
        (1, rand(1:subwidth), 1               , 0, rand(1:width), width        , 2),
        (2, rand(1:subwidth), 1               , 0, 1            , width        , 2),
        (3, subwidth        , rand(1:subwidth), 0, 1            , rand(1:width), 2),
        (4, subwidth        , rand(1:subwidth), 0, 1            , 1            , 2),
        (5, rand(1:subwidth), subwidth        , 0, rand(1:width), 1            , 2),
        (6, rand(1:subwidth), subwidth        , 0, width        , 1            , 2),
        (7, 1               , rand(1:subwidth), 0, width        , rand(1:width), 2),
        (8, 1               , rand(1:subwidth), 0, width        , width        , 2)
    ]
    for (c1,a1,b1,c2,a2,b2,t) in connections
        # if rand() <= 0.7
            n1, n2 = coordToLoc(a1,b1,c1), coordToLoc(a2,b2,c2)
            d = distanceCoord(urban.nodes[n1],urban.nodes[n2])
            urban.roads[n1,n2] = Road(n1,n2,d,t)
            urban.roads[n2,n1] = Road(n2,n1,d,t)
        # end
    end

    #Add highway around downtown
    for i in 1:(width-1), j in [1,width]
        n1, n2 = coordToLoc(j,i,0), coordToLoc(j,i+1,0)
        d = distanceCoord(urban.nodes[n1],urban.nodes[n2])
        urban.roads[n1,n2] = Road(n1,n2,d,1)
        urban.roads[n2,n1] = Road(n2,n1,d,1)

        n1, n2 = coordToLoc(i,j,0), coordToLoc(i+1,j,0)
        d = distanceCoord(urban.nodes[n1],urban.nodes[n2])
        urban.roads[n1,n2] = Road(n1,n2,d,1)
        urban.roads[n2,n1] = Road(n2,n1,d,1)
    end

    return Network(urban.nodes,urban.roads)
end

max_stars_repo_name: andrew-zhm/RoutingNetworks

max_issues_repo_name: andrew-zhm/RoutingNetworks

max_forks_repo_name: andrew-zhm/RoutingNetworks


text:
@testset "ece.jl" begin
    @testset "Trivial tests" begin
        ece = ECE(UniformBinning(10))

        # categorical distributions
        for predictions in ([[0, 1], [1, 0]], ColVecs([0 1; 1 0]), RowVecs([0 1; 1 0]))
            @test iszero(@inferred(ece(predictions, [2, 1])))
        end
        for predictions in (
            [[0, 1], [0.5, 0.5], [0.5, 0.5], [1, 0]],
            ColVecs([0 0.5 0.5 1; 1 0.5 0.5 0]),
            RowVecs([0 1; 0.5 0.5; 0.5 0.5; 1 0]),
        )
            @test iszero(@inferred(ece(predictions, [2, 2, 1, 1])))
        end

        # probabilities
        for predictions in ([0, 1], [0.0, 1.0])
            @test iszero(@inferred(ece(predictions, [false, true])))
        end
        @test iszero(@inferred(ece([0, 0.5, 0.5, 1], [false, false, true, true])))
    end

    @testset "Uniform binning: Basic properties" begin
        ece = ECE(UniformBinning(10))
        estimates = Vector{Float64}(undef, 1_000)

        # categorical distributions
        for nclasses in (2, 10, 100)
            dist = Dirichlet(nclasses, 1.0)

            predictions = [Vector{Float64}(undef, nclasses) for _ in 1:20]
            targets = Vector{Int}(undef, 20)

            for i in 1:length(estimates)
                rand!.(Ref(dist), predictions)
                targets .= rand.(Categorical.(predictions))

                estimates[i] = ece(predictions, targets)
            end

            @test all(x -> zero(x) < x < one(x), estimates)
        end

        # probabilities
        predictions = Vector{Float64}(undef, 20)
        targets = Vector{Bool}(undef, 20)
        for i in 1:length(estimates)
            rand!(predictions)
            map!(targets, predictions) do p
                rand() < p
            end
            estimates[i] = ece(predictions, targets)
        end
        @test all(x -> zero(x) < x < one(x), estimates)
    end

    @testset "Median variance binning: Basic properties" begin
        ece = ECE(MedianVarianceBinning(10))
        estimates = Vector{Float64}(undef, 1_000)

        for nclasses in (2, 10, 100)
            dist = Dirichlet(nclasses, 1.0)

            predictions = [Vector{Float64}(undef, nclasses) for _ in 1:20]
            targets = Vector{Int}(undef, 20)

            for i in 1:length(estimates)
                rand!.(Ref(dist), predictions)
                targets .= rand.(Categorical.(predictions))

                estimates[i] = ece(predictions, targets)
            end

            @test all(x -> zero(x) < x < one(x), estimates)
        end
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random, StatsBase, DataFrames
	Random.seed!(1)
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 2.1"

# ╔═╡ e4591176-1eca-11eb-1e04-15ae5da7dbb0

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
begin
	df = DataFrame()
	for n in 4:7
		N, faces = 10^n, 1:6

		numSol = sum([iseven(i+j) for i in faces, j in faces]) / length(faces)^2
		mcEst  = sum([iseven(rand(faces) + rand(faces)) for i in 1:N]) / N

		append!(df, DataFrame(:n => n, :N => N, :numsol => numSol, :mcEst => mcEst))
	end
	df
end

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 2.1"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═e4591176-1eca-11eb-1e04-15ae5da7dbb0
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3


text:
using Pkg
pkg"activate Burner"
# using Revise
using Burner
using Images
using VideoIO
using BSON

width = 512
height = 288

demodisp = Burner.Display(ones(3, height, width), 0.0015)
video = VideoIO.openvideo("oblivion.mp4")
imageshow = Float64.(Images.channelview(imresize(load("descargar-images-full-hd-1080p.jpg"), (height, width))))

Burner.burnVideo(demodisp, video, 1.0/30.0)
imagepresented = colorview(RGB, Burner.present(demodisp, imageshow))
imagepresented
bson("results.bson", Dict(:display => demodisp, :demoimage => imagepresented))

max_stars_repo_name: thekineticeffect/burninsim

max_issues_repo_name: thekineticeffect/burninsim

max_forks_repo_name: thekineticeffect/burninsim


text:
function tchanelSolve(
    p0,
    T,
    IRScale,
    UVScale,
    massfun,
    Imlambdaflowfun,
    Relambdaflowfun,
    ϵ = 0.1;
    u0,
    progress = true,
    atol = 1e-14,
    rtol = 1e-14,
    adaptive = true,
    dtmax = 0.01,
    kwargs...,
)
    function flow_t(du, u, p, k)
        #display(plot(x->-h2(x)[2],k,Λ))
        #println(k, " m2=", -u[2], " λ=", u[1])
        du[1] = Imlambdaflowfun(p, k, massfun(k), T, u[1], u[2], ϵ)
        du[2] = Relambdaflowfun(p, k, massfun(k), T, u[1], u[2], ϵ)
    end
    alg = Tsit5()
    kspan = (UVScale, IRScale)
    prob = ODEProblem(flow_t, u0, kspan, p0)
    solve(
        prob,
        alg,
        progress = progress,
        atol = atol,
        rtol = rtol,
        adaptive = adaptive,
        dtmax = dtmax,
        kwargs...,
    )
end

function tchanelSolveParallel(
    T,
    IRScale,
    UVScale,
    massfun,
    Imlambdaflowfun,
    Relambdaflowfun,
    numberOfParameters,
    ϵ = 0.1;
    u0,
    atol = 1e-14,
    rtol = 1e-14,
    adaptive = true,
    dtmax = 0.01,
    kwargs...,
)
    function flow_t(du, u, p, k)
        #display(plot(x->-h2(x)[2],k,Λ))
        #println(k, " m2=", -u[2], " λ=", u[1])
        du[1] = Imlambdaflowfun(p, k, massfun(k), T, u[1], u[2], ϵ)
        du[2] = Relambdaflowfun(p, k, massfun(k), T, u[1], u[2], ϵ)
    end
    alg = Tsit5()
    kspan = (UVScale, IRScale)
    p = 1.0
    prob = ODEProblem(flow_t, u0, kspan, p)
    parameterList = range(1.0, stop = 500.0, length = numberOfParameters)
    # function parameterChange(prob, i, repeat)
    #     prob.p = parameterList[i]
    #     prob
    # end
    parameterChange = (prob, i, repeat) -> remake(prob, p = parameterList[i])
    ensembleProb = EnsembleProblem(prob, prob_func = parameterChange)
    solve(
        ensembleProb,
        alg,
        EnsembleThreads(),
        trajectories = numberOfParameters,
        atol = atol,
        rtol = rtol,
        adaptive = adaptive,
        dtmax = dtmax,
        kwargs...,
    )
end

function tchanelSolveTwoPointParallel(
    T,
    IRScale,
    UVScale,
    massfun,
    lamdafun,
    Imlambdaflowfun,
    Relambdaflowfun,
    numberOfParameters,
    ϵ = 0.1;
    u0,
    atol = 1e-14,
    rtol = 1e-14,
    adaptive = true,
    dtmax = 0.01,
    kwargs...,
)
    function flow_t(du, u, p, k)
        #display(plot(x->-h2(x)[2],k,Λ))
        #println(k, " m2=", -u[2], " λ=", u[1])
        du[1] = Imlambdaflowfun(p, k, massfun(k), lamdafun(k), T, u[1], u[2], ϵ)
        du[2] = Relambdaflowfun(p, k, massfun(k), lamdafun(k), T, u[1], u[2], ϵ)
    end
    alg = Tsit5()
    kspan = (UVScale, IRScale)
    p = 0.0
    prob = ODEProblem(flow_t, u0, kspan, p)
    parameterList = range(0, stop = 400.0, length = numberOfParameters)
    # function parameterChange(prob, i, repeat)
    #     prob.p = parameterList[i]
    #     prob
    # end
    parameterChange = (prob, i, repeat) -> remake(prob, p = parameterList[i])
    function imrefun(sol, i)
        p0 = parameterList[i]
        imgamm2 =
            -hcubature(
                x -> twopointflow2(p0, x[1], massfun(x[1]), T, y -> sol(y)[1]),
                [1.0],
                [800.0],
                atol = 1e-6,
                rtol = 1e-6,
            )[1]
        regamm2 =
            hcubature(
                x -> twopointflow2(p0, x[1], massfun(x[1]), T, y -> sol(y)[2]),
                [1.0],
                [800.0],
                atol = 1e-6,
                rtol = 1e-6,
            )[1] + p0^2 - massfun(800.0)
        return ([imgamm2, regamm2], false)
    end
    ensembleProb = EnsembleProblem(
        prob,
        prob_func = parameterChange,
        output_func = imrefun,
    )
    solve(
        ensembleProb,
        alg,
        EnsembleThreads(),
        trajectories = numberOfParameters,
        atol = atol,
        rtol = rtol,
        adaptive = adaptive,
        dtmax = dtmax,
        kwargs...,
    )
end

function tchanelSolveTwoPointParallel3(
    p0::T,
    Temper::T,
    krang::RGscale,
    massfun::Function,
    lamdafun::Function,
    Imlambdaflowfun::Function,
    Relambdaflowfun::Function,
    numberOfParameters::Int,
    ϵ = 0.1;
    u0::AbstractArray{T},
    config::ODEConfig
) where {T<:Number}
    function flow_t(du, u, p, k)
        #display(plot(x->-h2(x)[2],k,Λ))
        #println(k, " m2=", -u[2], " λ=", u[1])
        du[1] =
            Imlambdaflowfun(p0, p, k, massfun(k), lamdafun(k), Temper, u[1], u[2], ϵ)
        du[2] =
            Relambdaflowfun(p0, p, k, massfun(k), lamdafun(k), Temper, u[1], u[2], ϵ)
    end
    alg = Tsit5()
    p = krang.IR
    prob = ODEProblem(flow_t, u0, krang.kspan, p)
    kList = range(krang.IR, stop = krang.UV, length = numberOfParameters)
    parameterChange =
        (prob, i, repeat) -> remake(
            prob,
            tspan = (krang.UV, kList[i]),
            p = Epi(kList[i], massfun(kList[i])),
        )
    ensembleProb = EnsembleProblem(
        prob,
        prob_func = parameterChange,
        output_func = (sol, i) -> (sol[end], false),
    )
    soltemp = solve(
        ensembleProb,
        alg,
        EnsembleSerial(),
        trajectories = numberOfParameters,
        atol = config.atol,
        rtol = config.rtol,
        adaptive = config.adaptive,
        dense = config.dense,
        save_on = config.save_on,
        save_start = config.save_start,
        save_end = config.save_end,
    )
    imdata = hcat(soltemp...)'[:, 1]
    redata = hcat(soltemp...)'[:, 2]
    imfun = Spline1D(kList, imdata)
    refun = Spline1D(kList, redata)
    imgamm2 =
        -hcubature(
            x -> twopointflow2(p0, x[1], massfun(x[1]), Temper, imfun),
            [krang.IR],
            [krang.UV],
            atol = 1e-6,
            rtol = 1e-6,
            initdiv = 100,
        )[1]
    regamm2 =
        hcubature(
            x -> twopointflow2(p0, x[1], massfun(x[1]), Temper, refun),
            [krang.IR],
            [krang.UV],
            atol = 1e-6,
            rtol = 1e-6,
            initdiv = 100,
        )[1] + p0^2 - massfun(krang.UV)
    return [imgamm2, regamm2]
end

function tchanelSolveFourPointParallel(
    k0::T,
    Temper::T,
    krang::RGscale,
    massfun::Union{Function,Spline1D},
    lamdafun::Union{Function,Spline1D},
    Imlambdaflowfun::Function,
    Relambdaflowfun::Function,
    numberOfParameters::Int,
    ϵ = 0.1;
    u0::AbstractArray{T},
    config::ODEConfig
) where {T<:Number}
    function flow_t(du, u, p, k)
        du[1] =
            Imlambdaflowfun(p, Epi(k0, massfun(k0)), k, massfun(k), lamdafun(k), Temper, u[1], u[2], ϵ)
        du[2] =
            Relambdaflowfun(p, Epi(k0, massfun(k0)), k, massfun(k), lamdafun(k), Temper, u[1], u[2], ϵ)
    end
    alg = Vern9()
    p = 0.0#p0 initial
    prob = ODEProblem(flow_t, u0, (krang.UV,k0), p)
    p0List = range(0.0, stop = 400.0, length = numberOfParameters)
    parameterChange =
        (prob, i, repeat) -> remake(
            prob,
            p = p0List[i],
        )
    ensembleProb = EnsembleProblem(
        prob,
        prob_func = parameterChange,
        output_func = (sol, i) -> (sol[end], false),
    )
    soltemp = solve(
        ensembleProb,
        alg,
        EnsembleSplitThreads(),
        trajectories = numberOfParameters,
        dtmax=config.dtmax,
        abstol = config.atol,
        reltol = config.rtol,
        calck = false,
        # dt=0.0000001,
        dtmin=1e-15,
        adaptive = config.adaptive,
        dense = config.dense,
        save_on = config.save_on,
        save_start = config.save_start,
        save_end = config.save_end,
        progress = config.progress
    )
    imdata = hcat(soltemp...)'[:, 1]
    redata = hcat(soltemp...)'[:, 2]
    return [p0List,imdata, redata]
end