text:
export Operator
export bandinds, bandrange, linsolve, periodic
export dirichlet, neumann
export ldirichlet,rdirichlet,lneumann,rneumann
export ldiffbc,rdiffbc,diffbcs
export domainspace,rangespace

abstract Operator{T} #T is the entry type, Float64 or Complex{Float64}

Base.eltype{T}(::Operator{T}) = T
Base.eltype{T}(::Type{Operator{T}}) = T
Base.eltype{OT<:Operator}(::Type{OT}) = eltype(supertype(OT))

# default entry type
# we assume entries depend on both the domain and the basis
# realdomain case doesn't use

op_eltype(sp::Space) = promote_type(eltype(sp),prectype(domain(sp)))
op_eltype_realdomain(sp::Space) = promote_type(eltype(sp),real(prectype(domain(sp))))

 #Operators are immutable
Base.copy(A::Operator) = A

## We assume operators are T->T
rangespace(A::Operator) = error("Override rangespace for $(typeof(A))")
domainspace(A::Operator) = error("Override domainspace for $(typeof(A))")
domain(A::Operator) = domain(domainspace(A))

isconstspace(::) = false
## Functionals
isafunctional(A::Operator) = size(A,1)==1 && isconstspace(rangespace(A))

isbandedbelow(A::Operator) = isfinite(bandinds(A,1))
isbandedabove(A::Operator) = isfinite(bandinds(A,2))

isbanded(A::Operator) = isbandedbelow(A) && isbandedabove(A)
isbandedblockbanded(::) = false

israggedbelow(A::Operator) = isbandedbelow(A) || isbandedblockbanded(A)

macro functional(FF)
    quote
        Base.size(A::$FF,k::Integer) = k==1?1:∞
        ApproxFun.rangespace(::$FF) = ConstantSpace()
        ApproxFun.isafunctional(::$FF) = true
        function ApproxFun.defaultgetindex(f::$FF,k::Integer,j::Integer)
            @assert k==1
            f[j]
        end
    end
end

Base.size(A::Operator) = (size(A,1),size(A,2))
Base.size(A::Operator,k::Integer) = k==1?dimension(rangespace(A)):dimension(domainspace(A))

# used to compute "end" for last index
function Base.trailingsize(A::Operator, n::Integer)
    if n > 2
        1
    elseif n==2
        size(A,2)
    elseif isinf(size(A,2)) || isinf(size(A,1))
        ∞
    else
        size(A,1)*size(A,2)
    end
end

Base.ndims(::Operator) = 2

## bandrange and indexrange
bandwidth(A::Operator) = bandwidth(A,1) + bandwidth(A,2) + 1
bandwidth(A::Operator,k::Integer) = k==1?-bandinds(A,1):bandinds(A,2)
bandwidths(A::Operator) = (bandwidth(A,1),bandwidth(A,2))
# we are always banded by the size
bandinds(A::Operator) = (1-size(A,1),size(A,2)-1)
bandinds(A,k::Integer) = bandinds(A)[k]
bandrange(b::Operator) = UnitRange(bandinds(b)...)
function bandrangelength(B::Operator)
    bndinds=bandinds(B)
    bndinds[end]-bndinds[1]+1
end

function columninds(b::Operator,k::Integer)
    ret = bandinds(b)

    (ret[1]  + k < 1) ? (1,(ret[end] + k)) : (ret[1]+k,ret[2]+k)
end

## Strides
# lets us know if operators decouple the entries
# to split into sub problems
# A diagonal operator has essentially infinite stride
# which we represent by a factorial, so that
# the gcd with any number < 10 is the number
Base.stride(A::Operator) =
    isdiag(A)?factorial(10):1

Base.isdiag(A::Operator) = bandinds(A)==(0,0)

## Construct operators

include("SubOperator.jl")

#
# Base.sparse(B::Operator,n::Integer)=sparse(BandedMatrix(B,n))
# Base.sparse(B::Operator,n::Range,m::Range)=sparse(BandedMatrix(B,n,m))
# Base.sparse(B::Operator,n::Colon,m::Range)=sparse(BandedMatrix(B,n,m))
# Base.sparse(B::Operator,n::Range,m::Colon)=sparse(BandedMatrix(B,n,m))

## geteindex

Base.getindex(B::Operator,k,j) = defaultgetindex(B,k,j)
Base.getindex(B::Operator,k) = defaultgetindex(B,k)

## override getindex.

defaultgetindex(B::Operator,k::Integer) = error("Override [k] for $(typeof(B))")
defaultgetindex(B::Operator,k::Integer,j::Integer) = error("Override [k,j] for $(typeof(B))")

# Ranges

defaultgetindex(op::Operator,kr::Range) = eltype(op)[op[k] for k in kr]
defaultgetindex(B::Operator,k::Range,j::Range) = AbstractMatrix(view(B,k,j))

defaultgetindex(op::Operator,k::Integer,j::Range) = reshape(eltype(op)[op[k,j] for j in j],1,length(j))
defaultgetindex(op::Operator,k::Range,j::Integer) = eltype(op)[op[k,j] for k in k]

# Colon casdes
defaultgetindex(A::Operator,kr::Range,::Colon) = view(A,kr,:)
defaultgetindex(A::Operator,::Colon,jr::Range) = view(A,:,jr)
defaultgetindex(A::Operator,::Colon,::Colon) = A

defaultgetindex(A::Operator,kr::AbstractCount,jr::AbstractCount) = view(A,kr,jr)
defaultgetindex(B::Operator,k::AbstractCount,::Colon) = B[k,1:end]
defaultgetindex(B::Operator,::Colon,j::AbstractCount) = B[1:end,j]

# FiniteRange gives the nonzero entries in a row/column
immutable FiniteRange end

# default is to use bandwidth
# override for other shaped operators
default_colstart(A::Operator, i::Integer) = min(max(i-bandwidth(A,2), 1), size(A, 2))
default_colstop(A::Operator, i::Integer) = min(i+bandwidth(A,1), size(A, 1))
default_rowstart(A::Operator, i::Integer) = min(max(i-bandwidth(A,1), 1), size(A, 1))
default_rowstop(A::Operator, i::Integer) = min(i+bandwidth(A,2), size(A, 2))

for OP in (:colstart,:colstop,:rowstart,:rowstop)
    defOP = parse("default_"*string(OP))
    @eval begin
        $OP(A::Operator,i::Integer) = $defOP(A,i)
        $OP(A::Operator,i::Infinity{Bool}) = ∞
    end
end

function defaultgetindex(A::Operator,::Type{FiniteRange},::Type{FiniteRange})
    if isfinite(size(A,1)) && isfinite(size(A,2))
        A[1:size(A,1),1:size(A,2)]
    else
        error("Only exists for finite operators.")
    end
end
defaultgetindex(A::Operator,::Type{FiniteRange},jr) =
    A[1:colstop(A,maximum(jr)),jr]

defaultgetindex(A::Operator,kr,::Type{FiniteRange}) =
    A[kr,1:rowstop(A,maximum(kr))]

## Composition with a Fun, LowRankFun, and ProductFun

defaultgetindex{BT,S,T}(B::Operator{BT},f::Fun{S,T}) = B*Multiplication(domainspace(B),f)
defaultgetindex{BT,S,M,SS,T}(B::Operator{BT},f::LowRankFun{S,M,SS,T}) =
    mapreduce(i->f.A[i]*B[f.B[i]],+,1:rank(f))
defaultgetindex{BT,S,V,SS,T}(B::Operator{BT},f::ProductFun{S,V,SS,T}) =
    mapreduce(i->f.coefficients[i]*B[Fun([zeros(promote_type(BT,T),i-1);
                                            one(promote_type(BT,T))],f.space[2])],
                +,1:length(f.coefficients))

# Convenience for wrapper ops
unwrap_axpy!(α,P,A) = BLAS.axpy!(α,view(parent(P).op,P.indexes[1],P.indexes[2]),A)
iswrapper(::)=false

macro wrappergetindex(Wrap)
    ret = quote
        Base.getindex(OP::$Wrap,k::Integer...) =
            OP.op[k...]

        BLAS.axpy!{T,OP<:$Wrap}(α,P::ApproxFun.SubOperator{T,OP},A::AbstractMatrix) =
            ApproxFun.unwrap_axpy!(α,P,A)
    end

    for TYP in (:(BandedMatrices.BandedMatrix),:(ApproxFun.BandedBlockBandedMatrix),
                :Matrix,:AbstractMatrix,:Vector,:AbstractVector)
        ret = quote
            $ret

            Base.convert{T,OP<:$Wrap}(::Type{$TYP},P::ApproxFun.SubOperator{T,OP}) =
                $TYP(view(parent(P).op,P.indexes[1],P.indexes[2]))
        end
    end

    esc(ret)
end

macro wrapper(Wrap)
    ret = quote
        ApproxFun.@wrappergetindex($Wrap)

        ApproxFun.iswrapper(::$Wrap) = true
    end
    for func in (:(ApproxFun.rangespace),:(ApproxFun.domainspace),
                 :(ApproxFun.bandinds),:(ApproxFun.domain),:(Base.stride),
                 :(ApproxFun.isbandedblockbanded),:(ApproxFun.isconstop),
                 :(ApproxFun.israggedbelow))
        ret = quote
            $ret

            $func(D::$Wrap) = $func(D.op)
        end
    end
    for func in (:(ApproxFun.bandwidth),:(ApproxFun.colstart),:(ApproxFun.colstop),
                    :(ApproxFun.rowstart),:(ApproxFun.rowstop))
        ret = quote
            $ret

            $func(D::$Wrap,k::Integer) = $func(D.op,k)
        end
    end
    esc(ret)
end

## Standard Operators and linear algebra

include("linsolve.jl")

include("spacepromotion.jl")
include("banded/banded.jl")
include("general/general.jl")

include("functionals/functionals.jl")
include("almostbanded/almostbanded.jl")

include("systems.jl")

include("qrfact.jl")
include("nullspace.jl")

## Conversion

Base.zero{T<:Number}(::Type{Operator{T}}) = ZeroOperator(T)
Base.zero{O<:Operator}(::Type{O}) = ZeroOperator(eltype(O))

Base.eye(S::Space) = IdentityOperator(S)
Base.eye(S::Domain) = eye(Space(S))

Base.convert{T}(A::Type{Operator{T}},f::Fun) =
    norm(f.coefficients)==0?zero(A):convert(A,Multiplication(f))

## Promotion

Base.promote_rule{N<:Number}(::Type{N},::Type{Operator}) = Operator{N}
Base.promote_rule{N<:Number}(::Type{UniformScaling{N}},::Type{Operator}) =
    Operator{N}
Base.promote_rule{S,N<:Number}(::Type{Fun{S,N}},::Type{Operator}) = Operator{N}
Base.promote_rule{N<:Number,O<:Operator}(::Type{N},::Type{O}) =
    Operator{promote_type(N,eltype(O))}
Base.promote_rule{N<:Number,O<:Operator}(::Type{UniformScaling{N}},::Type{O}) =
    Operator{promote_type(N,eltype(O))}
Base.promote_rule{S,N<:Number,O<:Operator}(::Type{Fun{S,N}},::Type{O}) =
    Operator{promote_type(N,eltype(O))}

Base.promote_rule{BO1<:Operator,BO2<:Operator}(::Type{BO1},::Type{BO2}) =
    Operator{promote_type(eltype(BO1),eltype(BO2))}

## Wrapper

#TODO: Should cases that modify be included?
typealias WrapperOperator Union{SpaceOperator,MultiplicationWrapper,DerivativeWrapper,IntegralWrapper,
                                    ConversionWrapper,ConstantTimesOperator,TransposeOperator}

# The following support converting an Operator to a Matrix or BandedMatrix

## BLAS and matrix routines
# We assume that copy may be overriden

BLAS.axpy!(a,X::Operator,Y::AbstractMatrix) = BLAS.axpy!(a,AbstractMatrix(X),Y)

# this is for operators that implement copy via axpy!

bzeros(S::Operator) = bzeros(eltype(S),size(S,1),size(S,2),bandwidth(S,1),bandwidth(S,2))
Base.zeros(S::Operator) = zeros(eltype(S),size(S,1),size(S,2))
bbbzeros(S::Operator) = bbbzeros(eltype(S),blockbandwidth(S,1),blockbandwidth(S,2),
                    subblockbandwidth(S,1),subblockbandwidth(S,2),
            blocklengthrange(rangetensorizer(S),1:size(S,1)),
            blocklengthrange(domaintensorizer(S),1:size(S,2)))

rzeros(S::Operator) = rzeros(eltype(S),size(S,1),Int[colstop(S,j) for j=1:size(S,2)])

banded_convert_axpy!(S::Operator) =
    BLAS.axpy!(one(eltype(S)),S,bzeros(S))
matrix_convert_axpy!(S::Operator) =
    BLAS.axpy!(one(eltype(S)),S,zeros(S))
bandedblockbanded_convert_axpy!(S::Operator) =
        BLAS.axpy!(one(eltype(S)),S,bbbzeros(S))

function Base.convert(::Type{Matrix},S::Operator)
   if isinf(size(S,1)) || isinf(size(S,2))
       error("Cannot convert $S to a AbstractVector")
   end

   eltype(S)[S[k,j] for k=1:size(S,1),j=1:size(S,2)]
end

Base.convert(::Type{BandedMatrix},S::Operator) = default_bandedmatrix(S)
Base.convert(::Type{RaggedMatrix},S::Operator) = default_raggedmatrix(S)

function Base.convert(::Type{Vector},S::Operator)
    if size(S,2) ≠ 1  || isinf(size(S,1))
        error("Cannot convert $S to a AbstractVector")
    end

    eltype(S)[S[k] for k=1:size(S,1)]
end

Base.convert(::Type{AbstractMatrix},S::Operator) = Matrix(S)

function Base.convert(::Type{AbstractMatrix},S::SubOperator)
    if isinf(size(S,1)) || isinf(size(S,2))
        throw(BoundsError())
    end
    if isbanded(parent(S))
        BandedMatrix(S)
    elseif isbandedblockbanded(parent(S))
        BandedBlockBandedMatrix(S)
    elseif israggedbelow(parent(S))
        RaggedMatrix(S)
    else
        Matrix(S)
    end
end

Base.convert(::Type{AbstractVector},S::Operator) = Vector(S)


text:
function one_shot_gradient_estimator(sel::K,
                                     params::P,
                                     v_mod::Function,
                                     v_args::Tuple,
                                     mod::Function,
                                     args::Tuple;
                                     scale = 1.0) where {K <: ConstrainedSelection, P <: Parameters}
    _, cl = simulate(params, v_mod, v_args...)
    obs = merge(cl, sel)
    _, mlw = score(obs, params, mod, args...)
    lw = mlw - get_score(cl)
    gs = get_parameter_gradients(params, cl, nothing, lw * scale)
    return gs, lw, cl
end

function multi_shot_gradient_estimator(sel::K,
                                       params::P,
                                       v_mod::Function,
                                       v_args::Tuple,
                                       mod::Function,
                                       args::Tuple;
                                       num_samples::Int = 5000) where {K <: ConstrainedSelection, P <: Parameters}
    cs = Vector{CallSite}(undef, num_samples)
    lws = Vector{Float64}(undef, num_samples)
    Threads.@threads for i in 1:num_samples
        _, cs[i] = simulate(params, v_mod, v_args...)
        obs = merge(cs[i], sel)
        ret, mlw = score(obs, params, mod, args...)
        lws[i] = mlw - get_score(cs[i])
    end
    ltw = lse(lws)
    L = ltw - log(num_samples)
    nw = exp.(lws .- ltw)
    gs = Gradients()
    Threads.@threads for i in 1:num_samples
        ls = L - nw[i]
        accumulate_parameter_gradients!(params, gs, cs[i], nothing, ls)
    end
    return gs, L, cs, nw
end

# ------------ Automatic differentiation variational inference ------------ #

function automatic_differentiation_variational_inference(sel::K,
                                                         params::P,
                                                         v_mod::Function,
                                                         v_args::Tuple,
                                                         mod::Function,
                                                         args::Tuple;
                                                         opt = ADAM(),
                                                         iters = 1000,
                                                         gs_samples = 100) where {K <: ConstrainedSelection, P <: Parameters}
    cls = Vector{CallSite}(undef, gs_samples)
    elbows = Vector{Float64}(undef, iters)
    Threads.@threads for i in 1 : iters
        elbo_est = 0.0
        gs_est = Gradients()
        for s in 1 : gs_samples
            gs, lw, cl = one_shot_gradient_estimator(sel, params, v_mod, v_args, mod, args; scale = 1 / gs_samples)
            elbo_est += lw / gs_samples
            gs_est += gs
            cls[s] = cl
        end
        elbows[i] = elbo_est
        params = update_parameters(opt, params, gs_est)
    end
    params, elbows, cls
end

ext: jl


text:
using MagnitudeDistributions, Distributions, Random
using Test

tests = ["grclassic",
         "utils",
         "completeness",
         "likelihood"]

printstyled("Running tests:\n", color=:blue)

for t in tests
    @testset "Test $t" begin
        Random.seed!(10)
        include("$t.jl")
    end
end

# print method ambiguities
println("Potentially stale exports: ")
display(Test.detect_ambiguities(MagnitudeDistributions))
println()


text:
module debug

using ..Ahorn, Maple

config = Dict{String, Any}()

const defaultIgnoredTooltipAttrs = String["x", "y", "width", "height", "originX", "originY", "nodes"]

function setConfig(fn::String, buffertime::Int=0)
    global config = Ahorn.attemptLoadConfig(fn, buffertime)
end

function log(s::String, shouldPrint::Bool)
    if shouldPrint
        println(s)
    end
end

function log(s::String, key::String)
    if get(config, key, false)
        println(s)
    end
end

function reloadTools!()
    empty!(Ahorn.loadedTools)
    append!(Ahorn.loadedTools, joinpath.(Ahorn.abs"tools", readdir(Ahorn.abs"tools")))
    Ahorn.initExternalTools()

    Ahorn.loadModule.(Ahorn.loadedTools)
    Ahorn.loadExternalModules!(Ahorn.loadedModules, Ahorn.loadedTools, "tools")
    Ahorn.changeTool!(Ahorn.loadedTools[1])
    Ahorn.selectRow!(Ahorn.toolList, 1)

    return true
end

function reloadEntities!()
    dr = Ahorn.getDrawableRoom(Ahorn.loadedState.map, Ahorn.loadedState.room)

    empty!(Ahorn.loadedEntities)
    append!(Ahorn.loadedEntities, joinpath.(Ahorn.abs"entities", readdir(Ahorn.abs"entities")))
    Ahorn.initExternalEntities()

    Ahorn.loadModule.(Ahorn.loadedEntities)
    Ahorn.loadExternalModules!(Ahorn.loadedModules, Ahorn.loadedEntities, "entities")
    Ahorn.registerPlacements!(Ahorn.entityPlacements, Ahorn.loadedEntities)

    Ahorn.getLayerByName(dr.layers, "entities").redraw = true
    Ahorn.selectRow!(Ahorn.roomList, row -> row[1] == Ahorn.loadedState.roomName)

    Ahorn.fillPlacementCache!(Ahorn.entityPlacementsCache, Ahorn.entityPlacements)

    empty!(Ahorn.entityNameLookup)

    return true
end

function reloadEffects!()
    empty!(Ahorn.loadedEffects)
    append!(Ahorn.loadedEffects, joinpath.(Ahorn.abs"effects", readdir(Ahorn.abs"effects")))
    Ahorn.initExternalEffects()

    Ahorn.loadModule.(Ahorn.loadedEffects)
    Ahorn.loadExternalModules!(Ahorn.loadedModules, Ahorn.loadedEffects, "effects")
    Ahorn.registerPlacements!(Ahorn.effectPlacements, Ahorn.loadedEffects)

    return true
end

function reloadTriggers!()
    dr = Ahorn.getDrawableRoom(Ahorn.loadedState.map, Ahorn.loadedState.room)

    empty!(Ahorn.loadedTriggers)
    append!(Ahorn.loadedTriggers, joinpath.(Ahorn.abs"triggers", readdir(Ahorn.abs"triggers")))
    Ahorn.initExternalTriggers()

    Ahorn.loadModule.(Ahorn.loadedTriggers)
    Ahorn.loadExternalModules!(Ahorn.loadedModules, Ahorn.loadedTriggers, "triggers")
    Ahorn.registerPlacements!(Ahorn.triggerPlacements, Ahorn.loadedTriggers)

    Ahorn.getLayerByName(dr.layers, "triggers").redraw = true
    Ahorn.selectRow!(Ahorn.roomList, row -> row[1] == Ahorn.loadedState.roomName)

    Ahorn.fillPlacementCache!(Ahorn.triggerPlacementsCache, Ahorn.triggerPlacements)

    return true
end

function clearMapDrawingCache!(map::Maple.Map=Ahorn.loadedState.map)
    Ahorn.loadChangedExternalSprites!()

    Ahorn.deleteDrawableRoomCache(map)
    Ahorn.draw(Ahorn.canvas)

    return true
end

function forceDrawWholeMap!(map::Maple.Map=Ahorn.loadedState.map)
    Ahorn.loadChangedExternalSprites!()

    ctx = Ahorn.Gtk.getgc(Ahorn.canvas)

    for room in map.rooms
        dr = Ahorn.getDrawableRoom(map, room)

        Ahorn.drawRoom(ctx, Ahorn.camera, dr, alpha=1.0)
    end

    Ahorn.draw(Ahorn.canvas)

    return true
end

function checkTooltipCoverage(ignore::Array{String, 1}=defaultIgnoredTooltipAttrs)
    missing = Dict{String, Array{String, 1}}()

    for (cache, placements) in ((Ahorn.entityPlacementsCache, Ahorn.entityPlacements), (Ahorn.triggerPlacementsCache, Ahorn.triggerPlacements))
        for (name, placement) in placements
            target = Ahorn.getCachedPlacement!(cache, placements, name)
            key = isa(target, Maple.Entity) ? "entities" : "triggers"
            tooltips = get(Ahorn.langdata, ["placements", key, target.name, "tooltips"])

            if !haskey(missing, target.name)
                missing[target.name] = String[]
            end

            for (attr, value) in target.data
                if !(attr in missing[target.name]) && !haskey(tooltips, Symbol(attr)) && !(attr in ignore)
                    push!(missing[target.name], attr)
                end
            end
        end
    end

    println("-- Missing Tooltips --")

    for (name, attrs) in missing
        if !isempty(attrs)
            println(name)
        end

        for attr in attrs
            println("    $attr")
        end
    end
end

function reloadLangdata()
    return Ahorn.loadLangfile()
end

# Put on a method definition to time every call of it.
macro timemachine(expr)
    @assert expr.head in (:function, :(=))
    Expr(expr.head, expr.args[1], quote
        println("Method call: ", $(string(expr.args[1])))
        @time $(expr.args[2])
    end) |> esc
end

end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
# compile markdowns
Literate.markdown("examples/src/ex1_pricing_pe_nofeedback.jl", "docs/src/examples/ex1_pricing_nofeedback"; documenter=true, execute = true, codefence = "```julia" => "```")
Literate.markdown("examples/src/ex2_pricing_pe_with_feedback.jl", "docs/src/examples/ex2_pricing_wfeedback"; documenter=true, execute = true, codefence = "```julia" => "```")
Literate.markdown("examples/src/ex3_Mackowiak_Wiederholt_2009.jl", "docs/src/examples/ex3_mw2009"; documenter=true, execute = true, codefence = "```julia" => "```")
Literate.markdown("examples/src/ex4_Sims_2010.jl", "docs/src/examples/ex4_sims2010"; documenter=true, execute = true, codefence = "```julia" => "```")
Literate.markdown("examples/src/ex5_Mackowiak_Matejka_Wiederholt_2018.jl", "docs/src/examples/ex5_mmw2018"; documenter=true, execute = true, codefence = "```julia" => "```")
Literate.markdown("examples/src/ex6_Afrouzi_Yang_2020.jl", "docs/src/examples/ex6_ay2020"; documenter=true, execute = true, codefence = "```julia" => "```")

function update_notebooks(content)
    content = replace(content, "# [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/afrouzi/DRIPs.jl/binder?filepath=examples) to run and modify the following code (no software is needed on the local machine).

" => "")
    return content
end

# compile jupyter notebooks
Literate.notebook("examples/src/ex1_pricing_pe_nofeedback.jl", "examples/notebooks/"; documenter = true, preprocess = update_notebooks)
Literate.notebook("examples/src/ex2_pricing_pe_with_feedback.jl", "examples/notebooks/"; documenter = true, preprocess = update_notebooks)
Literate.notebook("examples/src/ex3_Mackowiak_Wiederholt_2009.jl","examples/notebooks/"; documenter = true, preprocess = update_notebooks)
Literate.notebook("examples/src/ex4_Sims_2010.jl", "examples/notebooks/"; documenter = true, preprocess = update_notebooks)
Literate.notebook("examples/src/ex5_Mackowiak_Matejka_Wiederholt_2018.jl", "examples/notebooks/"; documenter = true, preprocess = update_notebooks)
Literate.notebook("examples/src/ex6_Afrouzi_Yang_2020.jl", "examples/notebooks"; documenter=true, preprocess = update_notebooks)


text:
function get_Ψ(S, Λ, Ω, tol, k=0)
    E = sum(Λ)
    Eᵣ = 0
    if k == 0
        k = 1
        while 1 - (Eᵣ + Λ[k])/E > tol
            Eᵣ += Λ[k]
            k+=1
        end
    end
    Ψ = S*Ω[:,1:k]
    for i in 1:k
        Ψ[:,i] ./= sqrt(abs.(Λ[i]))
    end
    return k, Ψ
end


text:
using LinearAlgebra

@testset "Core" begin

    points = [Point(1,1), Point(2,2), Point(3,3)]
    img = zeros(Gray{N0f8}, 3, 3)
    @inferred draw!(img, points)
    expected = Matrix(Diagonal([1,1,1]))
    @test all(img .== expected) == true

    img = zeros(Gray{N0f8}, 3, 3)
    @inferred draw!(img, points, Gray{N0f8}(0.8))
    expected = Matrix(Diagonal([0.8,0.8,0.8]))
    @test all(img .== expected) == true

    img = zeros(Gray{N0f8}, 3, 3)
    @inferred draw!(img, points, [Gray{N0f8}(0.8)])
    expected = Matrix(Diagonal([0.8,1,1]))
    @test all(img .== expected) == true

    points = [Point((1,1)), Point((2,2)), Point((3,3))]
    img = zeros(Gray{N0f8}, 3, 3)
    @inferred draw!(img, points)
    expected = Matrix(Diagonal([1,1,1]))
    @test all(img .== expected) == true

    ci = [CartesianIndex(1,2), CartesianIndex(1,3), CartesianIndex(3,1)]
    points = [Point(i) for i in ci]
    img = zeros(Gray{N0f8}, 3, 3)
    expected = copy(img)
    @inferred draw!(img, points)
    expected[1,2] = expected[1,3] = expected[3,1] = oneunit(Gray{N0f8})
    @test all(img .== expected) == true

    point = Point(4,1)
    img = zeros(Gray{N0f8}, 3, 3)
    @inferred draw!(img, point)
    expected = zeros(Gray{N0f8})
    @test all(img .== expected) == true
end


text:
using ITensors
using ITensorInfiniteMPS
using Test
using Random

#Ref time is 21.6s with negligible compilation time
@time @testset "vumpsmpo_ising" begin
  Random.seed!(1234)

  model = Model"ising"()
  model_kwargs = (J=1.0, h=1.2)

  function space_shifted(::Model"ising", q̃sz; conserve_qns=true)
    if conserve_qns
      return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1]
    else
      return [QN() => 2]
    end
  end

  # VUMPS arguments
  cutoff = 1e-8
  maxdim = 20
  tol = 1e-8
  maxiter = 50
  outer_iters = 3

  initstate(n) = "↑"

  # DMRG arguments
  Nfinite = 100
  sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true)
  Hfinite = MPO(model, sfinite; model_kwargs...)
  ψfinite = randomMPS(sfinite, initstate)
  sweeps = Sweeps(20)
  setmaxdim!(sweeps, 10)
  setcutoff!(sweeps, 1E-10)
  energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite, sweeps; outputlevel=0)
  Szs_finite = expect(ψfinite, "Sz")

  function energy(ψ, h, n)
    ϕ = ψ[n] * ψ[n + 1] * ψ[n + 2]
    return (noprime(ϕ * h) * dag(ϕ))[]
  end

  nfinite = Nfinite ÷ 2
  hnfinite = ITensor(model, sfinite[nfinite], sfinite[nfinite + 1]; model_kwargs...)
  orthogonalize!(ψfinite, nfinite)
  energy_finite = energy(ψfinite, hnfinite, nfinite)

  @testset "VUMPS/TDVP with: multisite_update_alg = $multisite_update_alg, conserve_qns = $conserve_qns, nsites = $nsites" for multisite_update_alg in
                                                                                [
      "sequential", "parallel"
    ],
    conserve_qns in [true, false],
    nsites in [1, 2],
    time_step in [-Inf]

    vumps_kwargs = (
      multisite_update_alg=multisite_update_alg,
      tol=tol,
      maxiter=maxiter,
      outputlevel=0,
      time_step=time_step,
    )
    subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)

    space_ = fill(space_shifted(model, 1; conserve_qns=conserve_qns), nsites)
    s = infsiteinds("S=1/2", nsites; space=space_)
    ψ = InfMPS(s, initstate)

    Hmpo = InfiniteMPOMatrix(model, s; model_kwargs...)
    # Alternate steps of running VUMPS and increasing the bond dimension
    ψ = tdvp(Hmpo, ψ; vumps_kwargs...)
    for _ in 1:outer_iters
      println("Subspace expansion")
      ψ = @time subspace_expansion(ψ, Hmpo; subspace_expansion_kwargs...)
      println("TDVP")
      ψ = @time tdvp(Hmpo, ψ; vumps_kwargs...)
    end

    @test norm(
      contract(ψ.AL[1:nsites]..., ψ.C[nsites]) - contract(ψ.C[0], ψ.AR[1:nsites]...)
    ) ≈ 0 atol = 1e-5
    #@test contract(ψ.AL[1:nsites]..., ψ.C[nsites]) ≈ contract(ψ.C[0], ψ.AR[1:nsites]...)

    H = InfiniteSum{MPO}(model, s; model_kwargs...)
    energy_infinite = expect(ψ, H)
    Szs_infinite = [expect(ψ, "Sz", n) for n in 1:nsites]

    @test energy_finite ≈ sum(energy_infinite) / nsites rtol = 1e-4
    @test Szs_finite[nfinite:(nfinite + nsites - 1)] ≈ Szs_infinite rtol = 1e-3
  end
end

@testset "vumpsmpo_extendedising" begin
  Random.seed!(1234)

  model = Model"ising_extended"()
  model_kwargs = (J=1.0, h=1.1, J₂=0.2)

  function space_shifted(::Model"ising_extended", q̃sz; conserve_qns=true)
    if conserve_qns
      return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1]
    else
      return [QN() => 2]
    end
  end

  # VUMPS arguments
  cutoff = 1e-8
  maxdim = 20
  tol = 1e-8
  maxiter = 20
  outer_iters = 4

  initstate(n) = "↑"

  # DMRG arguments
  Nfinite = 100
  sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true)
  Hfinite = MPO(model, sfinite; model_kwargs...)
  ψfinite = randomMPS(sfinite, initstate)
  sweeps = Sweeps(20)
  setmaxdim!(sweeps, 10)
  setcutoff!(sweeps, 1E-10)
  energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite, sweeps; outputlevel=0)
  Szs_finite = expect(ψfinite, "Sz")

  function energy(ψ, h, n)
    ϕ = ψ[n] * ψ[n + 1] * ψ[n + 2]
    return (noprime(ϕ * h) * dag(ϕ))[]
  end

  nfinite = Nfinite ÷ 2
  hnfinite = ITensor(
    model, sfinite[nfinite], sfinite[nfinite + 1], sfinite[nfinite + 2]; model_kwargs...
  )
  orthogonalize!(ψfinite, nfinite)
  energy_finite = energy(ψfinite, hnfinite, nfinite)

  for multisite_update_alg in ["sequential"],
    conserve_qns in [true],
    nsites in [1, 2],
    time_step in [-Inf]

    vumps_kwargs = (
      multisite_update_alg=multisite_update_alg,
      tol=tol,
      maxiter=maxiter,
      outputlevel=0,
      time_step=time_step,
    )
    subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)

    space_ = fill(space_shifted(model, 1; conserve_qns=conserve_qns), nsites)
    s = infsiteinds("S=1/2", nsites; space=space_)
    ψ = InfMPS(s, initstate)

    Hmpo = InfiniteMPOMatrix(model, s; model_kwargs...)
    # Alternate steps of running VUMPS and increasing the bond dimension
    ψ = tdvp(Hmpo, ψ; vumps_kwargs...)
    for _ in 1:outer_iters
      ψ = subspace_expansion(ψ, Hmpo; subspace_expansion_kwargs...)
      ψ = tdvp(Hmpo, ψ; vumps_kwargs...)
    end

    @test norm(
      contract(ψ.AL[1:nsites]..., ψ.C[nsites]) - contract(ψ.C[0], ψ.AR[1:nsites]...)
    ) ≈ 0 atol = 1e-5

    H = InfiniteSum{MPO}(model, s; model_kwargs...)
    energy_infinite = expect(ψ, H)
    Szs_infinite = [expect(ψ, "Sz", n) for n in 1:nsites]

    @test energy_finite ≈ sum(energy_infinite) / nsites rtol = 1e-4
    @test Szs_finite[nfinite:(nfinite + nsites - 1)] ≈ Szs_infinite rtol = 1e-3
  end
end

@testset "vumpsmpo_extendedising_translator" begin
  Random.seed!(1234)

  model = Model"ising_extended"()
  model_kwargs = (J=1.0, h=1.1, J₂=0.2)

  function space_shifted(::Model"ising_extended", q̃sz; conserve_qns=true)
    if conserve_qns
      return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1]
    else
      return [QN() => 2]
    end
  end

  # VUMPS arguments
  cutoff = 1e-8
  maxdim = 20
  tol = 1e-8
  maxiter = 20
  outer_iters = 4

  initstate(n) = "↑"

  # DMRG arguments
  Nfinite = 100
  sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true)
  Hfinite = MPO(model, sfinite; model_kwargs...)
  ψfinite = randomMPS(sfinite, initstate)
  sweeps = Sweeps(20)
  setmaxdim!(sweeps, 10)
  setcutoff!(sweeps, 1E-10)
  energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite, sweeps; outputlevel=0)
  Szs_finite = expect(ψfinite, "Sz")

  function energy(ψ, h, n)
    ϕ = ψ[n] * ψ[n + 1] * ψ[n + 2]
    return (noprime(ϕ * h) * dag(ϕ))[]
  end

  nfinite = Nfinite ÷ 2
  hnfinite = ITensor(
    model, sfinite[nfinite], sfinite[nfinite + 1], sfinite[nfinite + 2]; model_kwargs...
  )
  orthogonalize!(ψfinite, nfinite)
  energy_finite = energy(ψfinite, hnfinite, nfinite)

  temp_translatecell(i::Index, n::Integer) = ITensorInfiniteMPS.translatecelltags(i, n)

  for multisite_update_alg in ["sequential"],
    conserve_qns in [true],
    nsite in [1, 2, 3],
    time_step in [-Inf]

    vumps_kwargs = (
      multisite_update_alg=multisite_update_alg,
      tol=tol,
      maxiter=maxiter,
      outputlevel=0,
      time_step=time_step,
    )
    subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)

    space_ = fill(space_shifted(model, 1; conserve_qns=conserve_qns), nsite)
    s_bis = infsiteinds("S=1/2", nsite; space=space_)
    s = infsiteinds("S=1/2", nsite; space=space_, translator=temp_translatecell)
    ψ = InfMPS(s, initstate)

    Hmpo = InfiniteMPOMatrix(model, s; model_kwargs...)
    # Alternate steps of running VUMPS and increasing the bond dimension
    ψ = tdvp(Hmpo, ψ; vumps_kwargs...)
    for _ in 1:outer_iters
      ψ = subspace_expansion(ψ, Hmpo; subspace_expansion_kwargs...)
      ψ = tdvp(Hmpo, ψ; vumps_kwargs...)
    end

    @test norm(
      contract(ψ.AL[1:nsite]..., ψ.C[nsite]) - contract(ψ.C[0], ψ.AR[1:nsite]...)
    ) ≈ 0 atol = 1e-5

    H = InfiniteSum{MPO}(model, s; model_kwargs...)
    energy_infinite = expect(ψ, H)
    Szs_infinite = [expect(ψ, "Sz", n) for n in 1:nsite]

    @test energy_finite ≈ sum(energy_infinite) / nsite rtol = 1e-4
    @test Szs_finite[nfinite:(nfinite + nsite - 1)] ≈ Szs_infinite rtol = 1e-3

    #@test tags(s[nsite + 1]) == tags(s_bis[1 + 2nsite])
    @test ITensorInfiniteMPS.translator(ψ) == temp_translatecell
    @test ITensorInfiniteMPS.translator(s) == temp_translatecell
    @test ITensorInfiniteMPS.translator(Hmpo) == temp_translatecell
  end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Libdl, LibOSXUnwind_jll

@testset "LibOSXUnwind_jll" begin
    @static if Sys.isapple()
        @test dlsym(LibOSXUnwind_jll.libosxunwind_handle, :_Unwind_Backtrace; throw_error=false) !== nothing
    end
end

max_stars_repo_name: udesou/julia

max_issues_repo_name: udesou/julia

max_forks_repo_name: udesou/julia


text:
using GeometryTypes
using MeshAlgorithms
using Test
using StaticArrays

@testset "MeshAlgorithms.jl" begin
    # Write your own tests here.
    include("simplification.jl")
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 5c563230-df6e-4f65-a106-a707d2454fe7
using PlutoUI

# ╔═╡ 8bd774c7-1acf-4beb-bc99-4c303882d262
md"""
# Starving Neurons - Neural Networks and Gradient Starvation

Author: Peter Merseburger

## Introduction

Artificial intelligence (AI) and machine learning are disciplines whose methods have become ubiquitously applied in a broad range of scientific domains, enabling the implicit learning of features inherent to a dataset. A common AI tool for self-learning of data characteristics by examples are neural networks (NNs) which attempt to mimic biological neurons. NNs are broadly applicable, e.g. for tasks such as image and speech recognition, weather forecasting, stock market prediction, robotics and data mining among many others. But how do these artificial neural networks learn?
There is a myriad of neural network subtypes and topologies, but in this notebook, I want to focus on the multi-layer perceptron (MLP).

Basically speaking, a NN consists of multiple layers: an input layer, one or more hidden layers and an output layer (see fig. 1).

"""

# ╔═╡ 8f0cdbf1-0174-4d18-9072-2a8c5f015177
html"""
<img
	style="display: block;
           margin-left: auto;
           margin-right: auto;
           width: 50%;"
	src="https://upload.wikimedia.org/wikipedia/commons/thumb/e/e4/Artificial_neural_network.svg/538px-Artificial_neural_network.svg.png"
    alt="An image visualizing a basic NN topology with an input layer, a single hidden layer and an output layer">
</img>
"""

# ╔═╡ 6a79cfe1-12f9-4ae0-a917-14b4f3a1b31b
md"""
Figure 1: Basic NN/MLP topology with input, a single hidden and output layer.
CC-BY-SA 3.0. Source: [Wikimedia Commons](https://commons.wikimedia.org/wiki/File:Artificial_neural_network.svg)
"""

# ╔═╡ 512e6240-76e7-11ec-2511-07ffc4b74d27
md"""
The input layer encodes the data presented to the network and feeds it weighted into the first hidden layer; all other layers consist of artificial neurons with non-linear activation functions. MLPs are fully connected, meaning that every node of one layer connects to every node of the next layer and each connection has a weight assigned to it. The weights can be understood as the actual learned parameters being adjusted as they get passed from layer to layer. The MLP learns supervised, i.e. by presentation of labelled training data. Without going into too much detail, I want to briefly formalize the learning process by the example of a classification task, that is deciding to which discrete class a given input belongs. This is exactly the use case we will evaluate later on.
"""

# ╔═╡ 93bc80a8-673d-457b-96d8-e1bb0b09682e
md"""
### Supervised Learning – Learning from Examples

Supervised learning tries to learn a function which maps input to output values based on known input-output pairs (i.e. the pair ($x$, $y$); a correct classification $y$ for a given input $x$). The observed outputs $y$ originate from an unknown function $f(x) = y$. Thus, the goal of the learning process is to approximate $f(x)$ by a hypothesis $h(x) \approx f(x)$ while minimizing the error of the approximation, referred to as loss. More formally, the loss is defined as the utility lost by substituting the correct classification $f(x) = y$ by the hypothesis $h(x) = \hat{y}$:

$L(x, y, ŷ) = \text{Utility(result of using}\ y\ \text{given}\ x\\)
- \text{Utility(result of using}\ ŷ\ \text{given}\ x\\)$

Commonly applied loss functions are the mean absolute error and mean squared error loss. Especially popular for classification tasks is the cross-entropy loss, also called logarithmic loss:

$L(y, ŷ) = -\sum_{i}y_i\cdot \log(ŷ_i)$

The learning process maximizes the expected utility through choosing the optimal hypothesis by minimizing the loss function. In practice, this is estimated by the empirical loss observed on the set $X$ of $N$ training examples.

$\text{EmpiricalLoss}(h) = \frac{1}{N}\cdot\sum_{(x, y) \in X} L(y, h(x))$

But what about the weights I’ve introduced earlier as the “actual learned parameters”? Keep in mind that a hypothesis $h$ is a function and the weights simply represent the function coefficients. For the linear case that would mean $h(x) = w_1 x + w_0$.
In the multivariate case there is no closed-form solution for the minimum loss. Therefore, we are confronted with an optimization problem which can be solved by gradient descent. The weight update with learning rate $\alpha$ and the vector of weights $\bf{w}$ looks like this:

$w_i = w_i - \alpha\cdot \frac{\partial}{\partial w_i}\cdot \text{Loss}(\bf{w})$

"""

# ╔═╡ 66a3b5cb-a465-42fb-a8c8-404d574312d2
md"""
### Regularization

> "Entities should not be multiplied beyond necessity." – William of Ockham
This quote, known as Ockham’s razor and usually paraphrased as “simple solutions are better solutions” captures the core idea of regularization. The aim of regularization is to find the hypothesis which minimizes the total cost, defined as the sum of empirical loss and its complexity:

$\text{Cost}(h) = \text{EmpiricalLoss}(h) + \lambda\cdot \text{Complexity}(h)$

Regularization penalizes complex hypotheses (mind the added regularization term in the formula $\lambda\cdot \text{Complexity}(h)$) and favors more regular functions. It limits the values that the weights can take during the learning process. Therefore, regularization avoids overfitting of the data.
The estimated best hypothesis $ĥ^*$ from the hypothesis space $\mathcal{H}$ now becomes:

$ĥ^* = \arg\min_{h \in \mathcal{H}} \text{Cost}(h)$

### Optimization
Given the major workings of supervised learning, it should have become clear that learning is a mathematical optimization problem. A classic way to achieve the minimization of the cost is the application of gradient descent (GD). However, the path the gradient take on the surface of the solution space can limit its ability to find the global optimum. To conquer that, the steps can be modified by additional optimization methods. Common choices for that purpose are the addition of momentum, resulting in stochastic gradient descent (SGD); an adaptive learning rate, different forms of parameter initialization as well as normalization techniques.
"""

# ╔═╡ 4a14831d-7736-4af5-875f-05ae6e3cf7c4
md"""
### Side Effects of Gradient-based Learning
Despite its popularity, gradient descent learning in NNs can come with some unwanted side effects such as the vanishing gradient problem, which condenses down to the gradient being too small to effectively update a weight and thus preventing learning. This phenomenon is complemented in the other direction by exploding gradients.

Another not yet published but currently in preprint described observation is the effect of [**_Gradient Starvation_ (GS)**](https://arxiv.org/abs/2011.09468). Pezeshki et al. coined and formalized this phenomenon of over-parameterized NNs in their work. According to them, _Gradient Starvation_ is a predisposition of NNs learned by cross-entropy, which occurs when the loss minimization is driven by only a feature subset regardless of the presence of other meaningful features. Thus, the NN becomes biased towards superficial features. Beyond demonstrating the effects of GS, Pezeshki et al. introduce **_Spectral Decoupling_ (SD)** regularization as a means to counteract gradient starvation.
"""

# ╔═╡ 60905936-9373-44de-96e6-cb5d6dd72d33
md"""
### Spectral Decoupling
Proceeding from the popular "ridge-regularized cross entropy":

$\mathcal{L}(\boldsymbol{θ}) = \boldsymbol{1} \cdot \log(1 + \exp(-\boldsymbol{Y}\boldsymbol{ŷ})) + \frac{λ}{2} ||\boldsymbol{θ}||^2$

with $\boldsymbol{Y}$ the diagonal matrix of $y$ and $\boldsymbol{θ}$ the concatenation of all layer's weights, SD is a variation where the L2 regularization term is replaced by a penalty on the raw predictions $\boldsymbol{ŷ}$:

$\mathcal{L}(\boldsymbol{θ}) = \boldsymbol{1} \cdot \log(1 + \exp(-\boldsymbol{Y}\boldsymbol{ŷ})) + \frac{λ}{2} ||\boldsymbol{ŷ}||^2$

This, as stated by the Pezeshki et al., uncouples the learning of multiple features. According to their work, a feature trained close to its optimum inhibits (starves) the training of other features which is thus prevented by SD.
"""

# ╔═╡ 00d4161d-6f8f-42d6-8f13-7d5aa7e70ab3
md"""
### Objective

My project and thereby this notebook aims to reproduce and demonstrate the effects of GS on a simple 2D classification task. Furthermore, the proposed SD regularization method is implemented and compared to other methods. On top of that, spectral decoupling is not only applied in the context of SGD but all other investigated optimization approaches.
"""

# ╔═╡ ac02597f-a7ca-4ede-ab8a-c0739c60d0d5
md"""
### Waning or Waxing Moon? A 2D Classification Experiment
A modified version of the [Moons dataset](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_moons.html) was used with two different topologies as exemplified in figure 1.
The first topology allows to draw a line by a very small margin between the two moons. Therefore, it represents the fully linearly separable case (fig. 1, left). The second one applies a smaller offset to interleave both moons which only allows to separate them by learning a curved classification boundary (fig. 1, right). The task of the NN is to discriminate if a given data point belongs to the waning (blue) or waxing moon (red).
"""

# ╔═╡ 2f766a9c-de6e-4f95-aa73-aacb486aaf61
md"""
$(Resource("https://raw.githubusercontent.com/justinsane1337/GradientStarvation/master/plots/figure_1_edited.png"))
Figure 1: Two moons topologies; one linear separable (offset Δ1.0), the other separable by a curved boundary (offset Δ0.5).
"""

# ╔═╡ 7f7ed965-fa08-4d6f-8908-dfc8876bdd02
md"""
The reasoning of the researchers was that in the linear separable case the NN is not encouraged to learn a curved decision boundary because the loss becomes negligible by only discriminating the two moons along one axis, therefore learning only one feature and neglecting the other. In contrast to that, the interleaved moons can only be separated by learning both features which better accommodates the curvatures of the data structure in general.
They observed that a different choice of regularization and optimization methods like weight decay (WD) and _Adaptive Moment Estimation_ (ADAM), among others, are not able to stimulate the NN to learn a curved decision boundary.
Moreover, Pezeshki et al. explained that circumstance by cross-entropy loss learning which starves the gradients from the other feature, only reinforcing the superficial one.

#### Repercussions of Gradient Starvation
Why does the shape of the decision boundary even matter as long as it discriminates the data points confidently? As stated by the authors, GS results in a very small distance between the data points and the boundary which translates into a lack of robustness when generalizing to new data. Contrarily, GS could also have a positive effect in terms of preventing overfitting by not learning non-dominant features. The authors made the comparison to the aformentioned Ockham's razor: GS leads to simpler decision boundaries and hence NNs that generalize better, thus GS being a form of implicit regularization.

#### Implementation
The implementation is based on [Flux](https://fluxml.ai/), a ML framework for Julia. The NN consists of two hidden layers with 500 nodes each. Rectified Linear Unit (ReLu) is used as an activation function and learning is facilitated by cross-entropy loss for 1000 epochs on a data set with 300 observations (150 blue and 150 red points). In contrast to Pezeshki et al., I've opted to train in batches of 50 (6 batches in total) since this is a common practice for training NNs.
All networks are trained with gradient descent but with different regularization and optimization approaches. This yields the following different setups:
- Gradient descent (GD; no regularization/optimization)
- Gradient descent with momentum (stochastic gradient descent (SGD); no regularization)
- Gradient descent with weight decay (WD; L2 regularization)
- Gradient descent with _Adaptive Moment Estimation_ (ADAM; no regularization)

Different sets of hyperparameters were evaluated while tinkering with the NNs. However, since the main objective is to reproduce the research conducted by Pezeshki et al., their set of hyperparameters were used:
- learning rate = 0.01
- weight decay coefficient = 0.01
- spectral decoupling coefficient = 0.003
"""

# ╔═╡ 88ddf822-622b-4d46-a0a4-a7d749c8bc48
md"""
#### Results
"""

# ╔═╡ 884bc572-9daa-40ad-8d76-56ec58728f63
md"""
##### Linear separable case
$(Resource("https://raw.githubusercontent.com/justinsane1337/GradientStarvation/master/plots/GD_%CE%941.0.png"))
Figure 2: GD, moon offset = Δ1.0. The decision boundary is almost linear and the loss continuously falls while training. Accuracy reaches 100%.

$(Resource("https://raw.githubusercontent.com/justinsane1337/GradientStarvation/master/plots/SGD_%CE%941.0.png"))
Figure 3: SGD, moon offset = Δ1.0. The decision boundary is similar to the one learned with GD but equally distant from both classes. Loss declines quicker and reaches an overall lower level compared to GD. Accuracy reaches 100% after less training.

$(Resource("https://raw.githubusercontent.com/justinsane1337/GradientStarvation/master/plots/WD_%CE%941.0.png"))
Figure 4: WD, moon offset = Δ1.0, WD coefficient = 0.01. The boundary is similar to GD without WD.

$(Resource("https://raw.githubusercontent.com/justinsane1337/GradientStarvation/master/plots/ADAM_%CE%941.0.png"))
Figure 5: ADAM, moon offset = Δ1.0. The decision boundary is slightly curved around the data points closest to the boundary. Loss converges to zero almost immediately while full accuracy is attained.
"""

# ╔═╡ ac453412-0f71-45ee-aefb-0cc519467c6c
md"""
As we can see from figures 2 to 5, the decision boundaries of all networks exhibit some specificities. The GD network learned an almost linear decision boundary which appears to be closer to the blue-labeled moon. The boundary looks similar with SGD, but its distance to both moons is equal. Only the network with ADAM optimization learns a slight curvature around the points closest to the boundary. Furthermore, the ADAM classification learned the fastest as can be seen from the rapid decrease in loss and increase in accuracy. At the same time, the contour lines all closely condense around the decision boundary which indicates that the ADAM network is most confident in discriminating both classes.
"""

# ╔═╡ 5e282e1a-62fc-4911-9b7e-f2631eeb2669
md"""
##### Linear inseparable case

$(Resource("https://raw.githubusercontent.com/justinsane1337/GradientStarvation/master/plots/GD_%CE%940.5.png"))
Figure 6: GD, moon offset = Δ0.5. The network fails to learn a decision boundary to separate the two classes within the 1000 training iterations.

$(Resource("https://raw.githubusercontent.com/justinsane1337/GradientStarvation/master/plots/SGD_%CE%940.5.png"))
Figure 7: SGD, moon offset = Δ0.5. The network learns a curved decision boundary which is very close to the interleaving data points. Loss becomes negligible during learning.

$(Resource("https://raw.githubusercontent.com/justinsane1337/GradientStarvation/master/plots/WD_%CE%940.5.png"))
Figure 8: WD, moon offset = Δ0.5. The network fails to learn a decision boundary to separate the two classes. Again, similar result compared to GD.

$(Resource("https://raw.githubusercontent.com/justinsane1337/GradientStarvation/master/plots/ADAM_%CE%940.5.png"))
Figure 9: ADAM, moon offset = Δ0.5. The network learns a curved decision boundary. Loss fluctuates especially during the early iterations of training.
"""

# ╔═╡ f7bae71f-885a-404d-985b-92d612e1f3d6
md"""
SGD and ADAM are able to learn a curved decision boundary (fig. 7 and 9) while GD and WD fail to do so (fig. 6 and 8). Moreover, the decision boundary learned by ADAM better resembles the data structure. The contours of the ADAM network show a very strong confidence on the networks classification. This contrasts the results from the GS paper, where ADAM was not able to learn a decision boundary as curved and well-fitting the data structure as shown in fig. 9.
"""

# ╔═╡ 2bfb5f10-cb87-41a4-9dca-f3f0b1157c6f
md"""
#### Spectral Decoupling

$(Resource("https://raw.githubusercontent.com/justinsane1337/GradientStarvation/master/plots/GD_%CE%940.5_%2B_SD.png"))
Figure 10: GD with SD, moon offset = Δ0.5. No curved decision boundary is learned but the loss also did not converge.

$(Resource("https://raw.githubusercontent.com/justinsane1337/GradientStarvation/master/plots/SGD_%CE%940.5_%2B_SD.png"))
Figure 11: SGD with SD, moon offset = Δ0.5. The network learns a curved decision boundary which greatly resembles the moon data structure and shows invaginations at the data points of both moons closest to the boundary. Loss converges to ~15%.

$(Resource("https://raw.githubusercontent.com/justinsane1337/GradientStarvation/master/plots/WD_%CE%940.5_%2B_SD.png"))
Figure 12: WD with SD, moon offset = Δ0.5. NN fails to separate the two data classes. Similar to GD with SD, loss could still be lowered with further training.

$(Resource("https://raw.githubusercontent.com/justinsane1337/GradientStarvation/master/plots/ADAM_%CE%940.5_%2B_SD.png"))
Figure 13: ADAM with SD, moon offset = Δ0.5. The network learns a curved decision boundary which greatly resembles the data structure. Loss heavily oscillates and steadily increases with longer training. Interestingly, all data points of each respective class belong to the same contour.
"""

# ╔═╡ 7460f609-c32b-4704-a361-71c98bb626cd
md"""
SD, originally applied as a regularization by Pezeshki et al. only with SGD, is able to yield curved decision boundaries which greatly resemble the data structure of the two moons in combination with SGD and ADAM optimization. While SGD learns a more angular boundary with low loss, ADAM learns a smoother one but exhibits a high loss with strong oscillations. Interestingly, ADAM learns to group all data points of a class within the same contour (fig. 13). GD and WD do not profit from SD regularization and don't learn a curved boundary. However, their loss profiles suggest improvement with further learning.
"""

# ╔═╡ 87cc8172-d5cc-49ff-9a0e-fa734e1010c1
md"""
#### Additional experiments
Given that the loss did not converge for GD + SD after 1000 epochs, the network was trained for 10000 iterations.

$(Resource("https://raw.githubusercontent.com/justinsane1337/GradientStarvation/master/plots/GD_%CE%940.5_%2B_SD_10000.png"))
Figure 14: GD with SD, moon offset = Δ0.5. NN was trained 10 times longer (10000 iterations). The learned curved decision boundary resembles the result obtained with SGD + SD.

To conquer the overfitting observed with ADAM + SD lower learning rates were evaluated.

$(Resource("https://raw.githubusercontent.com/justinsane1337/GradientStarvation/master/plots/1e-3_ADAM_%CE%940.5_%2B_SD.png"))
Figure 15: ADAM with SD on the moons with offset = Δ0.5 and a learning rate of 0.001 still suffers from severe overfitting. Nonetheless, oscillations in loss are dampened.

$(Resource("https://raw.githubusercontent.com/justinsane1337/GradientStarvation/master/plots/1e-4_ADAM_%CE%940.5_%2B_SD.png"))
Figure 16: ADAM with SD on the moons with offset = Δ0.5 and a learning rate of 0.0001 shows lower loss compared to 0.001 but a decision boundary with a smaller margin.

The additional experiments show that GD with SD performs equally well compared to SGD with SD but needs much more training for that.
The overfitting problem with ADAM + SD can be successfully conquered with lower learning rates. However, the shape of the decision boundary suffers with learning rates as low as 0.0001. From all three ADAM + SD experiments, the one learning at a rate of 0.001 performed best.
"""

# ╔═╡ 2c740919-1fd4-45c4-91d1-ef3d3e081139
md"""
## Discussion
Pezeshki et al. described gradient starvation as a learning proclivity that can exhibit both negative and positive effects on neural network learning. Considering the 2D classification task discussed in this notebook, the effect of learning a decision boundary with a small margin, thus supposedly hardly robust, could be exemplified for stochastic gradient descent. Combining SGD with spectral decoupling regularization was greatly able to enlarge the margin, therefore yielding a more robust network that not only focussed on a single dominant feature. On the other hand, the potential negative effect of overfitting by precisely learning less dominant features as well became evident through higher losses (SGD: 8.61% without vs. 13.37% with SD; ADAM: 2.00% without vs. 68.16%).
Contrarily to the results from Pezeshki et al., the here used neural network was perfectly able to learn curved decision boundaries featuring a large margin between the data points using ADAM without SD. While Pezeshki et al. trained their ADAM network with rates of 0.0001 and 0.001, I trained all networks under the same conditions, hence both SGD and ADAM with a learning rate of 0.01. Despite the higher learning rate, ADAM showed no sign of overfitting without SD. However, with SD the overfitting was extreme but ADAM coupled with SD regularization was also a combination of methods not assessed by the researchers.
Regarding the same matter, my approach to apply SD not only to SGD but all other optimizers simultaneously offered limited insight. First of all, the combination of WD (L2 regularization) and SD was expected to not work well since both methods are somewhat conflicting. ADAM performed well with SD but suffered from overfitting which can be conquered by a lower learning rate as demonstrated. Nonetheless, the network with purportedly lower loss learns a boundary with a smaller margin.

All in all, the results of Pezeshki et al. could be reproduced with the exception of the network using ADAM, which was able to learn decision boundaries with similar characteristics as the ones learned by SGD with SD.

Finally, it must be mentioned that Pezeshki et al. not only exemplified GS on the rather synthetical moon classification task but showed GS ramifications and the potential of SD as a regularization technique on more complex, well-established datasets like CIFAR, colored MNIST and CelebA which was sadly out of scope of this project.
"""

# ╔═╡ 21854d7b-13ce-477c-bc61-1cb15536e24f
md"""
### Conclusion
Gradient starvation is another relevant side effect of gradient descent-learned neural networks that may be worth the attention in the optimization process of NNs. However, if GS is a blessing or a curse depends on the use case. Starved models, i.e. models concentrating on superficial features during learning, can be less prone to overfitting while adressing GS promises to increase robustness.
"""

# ╔═╡ c95db58c-9d47-4ef0-9ede-234d7a693f41
md"""
### References
Pezeshki, M., Kaba, S. O., Bengio, Y., Courville, A., Precup, D., & Lajoie, G. (2020). Gradient starvation: A learning proclivity in neural networks. arXiv preprint arXiv:2011.09468.

Pizurica, A. (2020). E016330: Artificial Intelligence – Learning from examples [PDF course notes]. Ufora.

Pizurica, A. (2020). E016330: Artificial Intelligence - Neural networks [PDF course notes]. Ufora.

Russell, S., & Norvig, P. (2002). Artificial intelligence: a modern approach.

Kingma, D. P., & Ba, J. (2014). Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980.
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.32"
"""

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

julia_version = "1.7.0"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

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

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

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

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

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

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

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

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

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "0b5cfbb704034b5b4c1869e36634438a047df065"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.1"

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

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "ae6145ca68947569058866e443df69587acc1806"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.32"

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

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

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

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

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

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

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

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

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

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

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

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

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

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

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

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

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

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

# ╔═╡ Cell order:
# ╟─5c563230-df6e-4f65-a106-a707d2454fe7
# ╟─8bd774c7-1acf-4beb-bc99-4c303882d262
# ╟─8f0cdbf1-0174-4d18-9072-2a8c5f015177
# ╟─6a79cfe1-12f9-4ae0-a917-14b4f3a1b31b
# ╟─512e6240-76e7-11ec-2511-07ffc4b74d27
# ╟─93bc80a8-673d-457b-96d8-e1bb0b09682e
# ╟─66a3b5cb-a465-42fb-a8c8-404d574312d2
# ╟─4a14831d-7736-4af5-875f-05ae6e3cf7c4
# ╟─60905936-9373-44de-96e6-cb5d6dd72d33
# ╟─00d4161d-6f8f-42d6-8f13-7d5aa7e70ab3
# ╟─ac02597f-a7ca-4ede-ab8a-c0739c60d0d5
# ╟─2f766a9c-de6e-4f95-aa73-aacb486aaf61
# ╟─7f7ed965-fa08-4d6f-8908-dfc8876bdd02
# ╟─88ddf822-622b-4d46-a0a4-a7d749c8bc48
# ╟─884bc572-9daa-40ad-8d76-56ec58728f63
# ╟─ac453412-0f71-45ee-aefb-0cc519467c6c
# ╟─5e282e1a-62fc-4911-9b7e-f2631eeb2669
# ╟─f7bae71f-885a-404d-985b-92d612e1f3d6
# ╟─2bfb5f10-cb87-41a4-9dca-f3f0b1157c6f
# ╟─7460f609-c32b-4704-a361-71c98bb626cd
# ╟─87cc8172-d5cc-49ff-9a0e-fa734e1010c1
# ╟─2c740919-1fd4-45c4-91d1-ef3d3e081139
# ╟─21854d7b-13ce-477c-bc61-1cb15536e24f
# ╟─c95db58c-9d47-4ef0-9ede-234d7a693f41
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
function fcsc_simple{T}(X::Array{T, 3}, dictSize::Int, filterSize::Int, optParams::FCSCOpts{T})
# Input:
#   X: HxWxN stack of images (padded already)
#
# Alg. parameters:
#   mu_max: 1e5
#   tau: [1.01 1.5] larger means faster convergence but sacrifices primal feasibility,
#        smaller means slower convergence

  const mu_max = 1e5;
  const tau = 1.01;

  const Dinit    = optParams.Dinit;
  const beta::T  = optParams.beta;
        mu_s::T  = optParams.mu_s;
        mu_t::T  = optParams.mu_t;
  const max_iter = optParams.max_iter;
  const tol      = optParams.tol;
  const DEBUG    = optParams.DEBUG;
  const FIXED_D  = optParams.FIXED_D;

  (convH, convW, numImgs) = size(X);
  convLen = convH*convW;
  assert(!FIXED_D || !isempty(Dinit), "Dinit must be provided when FIXED_D is enabled.");

  # Initialization
  mask = ones(T, convH, convW, dictSize);
  mask[1:filterSize, 1:filterSize, :] = 0;
  mask = find(mask);

  if !FIXED_D;
    if isempty(Dinit);
      Dinit = reshape(mod([convert(T, 1):convH*convW*dictSize], 13), (convH, convW, dictSize));
    else;
      Dinit = padarray(Dinit, (convH-size(Dinit, 1), convW-size(Dinit, 2)), convert(T, 0), "post");
    end
    d = Dinit;
    s = copy(d);
    s[mask] = 0;
    normalize!(s);
  else;
    d = padarray(Dinit, (convH-size(Dinit, 1), convW-size(Dinit, 2)), convert(T, 0), "post");
    s = d;
  end
  lambda_s = zeros(T, convH, convW, dictSize);

  # Transform images to the frequency domain
  x_fft = reshape(fft2(X), convLen, numImgs);
  z = zeros(T, convLen*dictSize, numImgs);
  t = zeros(T, convLen*dictSize, numImgs);
  lambda_t = zeros(T, convLen*dictSize, numImgs);

  @printf("iter=000");
  if DEBUG;
    debugData = {"x"=>X,
                 "d"=>d, "s"=>s,
                 "lambda_s"=>lambda_s,
                 "mu_s"=>mu_s,
                 "z"=>z, "t"=>t,
                 "lambda_t"=>lambda_t,
                 "mu_t"=>mu_t,
                 "beta"=>beta,
                 "function"=>compute_lagrangian_dictonly};
    (prev_obj,rerr,spar) = compute_lagrangian(debugData);
    @printf(": obj=%.10f (rerr=%.10f, spar=%.10f)\n", prev_obj, rerr, spar);
  end

  d_fft = reshape(fft2(d), convLen, dictSize);
  D_fft = zeros(Complex{T}, convLen, convLen*dictSize);
  for k=1:dictSize;
    D_fft[:, (k-1)*convLen+[1:convLen]] = diagm(d_fft[:, k]);
  end
  DtD_fft = D_fft'*D_fft;
  for j=1:numImgs;
    #
    # subproblem z
    Dtx_fft = D_fft'*x_fft[:, j];

    # Initialize z to pure reconstruction solution
    z_fft = ( DtD_fft + mu_t.*eye(convLen*dictSize) ) \ Dtx_fft;
    z[:, j] = vec(real(ifft2(reshape(z_fft, convH, convW, dictSize))));

    #
    # subproblem t
    t[:, j] = shrink(z[:, j], beta/mu_t);
  end

  iter = 0;
  converged = false;
  while ~converged;
    iter += 1;

    @printf("\b\b\b\b\b\b\b\biter=%03d", iter);
    if DEBUG;
      debugData["z"] = z; debugData["t"] = t;
      debugData["lambda_t"] = lambda_t;
      debugData["mu_t"] = mu_t;
      debugData["d"] = d; debugData["s"] = s;
      debugData["lambda_s"] = lambda_s;
      debugData["mu_s"] = mu_s;
      (prev_obj,rerr,spar) = compute_lagrangian(debugData);
      @printf(": obj=%.10f (rerr=%.10f, spar=%.10f)\n", prev_obj, rerr, spar);
    end

    #
    # CODING
    #
    d_fft = reshape(fft2(d), convLen, dictSize);
    D_fft = zeros(Complex{T}, convLen, convLen*dictSize);
    for k=1:dictSize;
      D_fft[:, (k-1)*convLen+[1:convLen]] = diagm(d_fft[:, k]);
    end
    DtD_fft = D_fft'*D_fft;
    for j=1:numImgs;
      #
      # subproblem z
      Dtx_fft = D_fft'*x_fft[:, j];
      t_fft = vec(fft2(reshape(t[:, j], convH, convW, dictSize)));
      lambda_t_fft = vec(fft2(reshape(lambda_t[:, j], convH, convW, dictSize)));

      z_fft = ( DtD_fft + mu_t.*eye(convLen*dictSize) ) \ ( Dtx_fft + mu_t.*t_fft - lambda_t_fft );
      z[:, j] = vec(real(ifft2(reshape(z_fft, convH, convW, dictSize))));

      #
      # subproblem t
      t[:, j] = shrink(z[:, j]+lambda_t[:, j]./mu_t, beta/mu_t);
    end
    if DEBUG;
      debugData["z"] = z;
      debugData["t"] = t;
      (obj,rerr,spar) = compute_lagrangian(debugData);
      is_mono = (prev_obj>=obj) ? 'O' : 'X';
      @printf("          C: %c obj=%.10f (rerr=%.10f, spar=%.10f)\n", is_mono, obj, rerr, spar);
      prev_obj = obj;
    end

    #
    # LEARNING
    #
    if !FIXED_D;
      #
      # subproblem d
      Z_fft = zeros(Complex{T}, convLen, convLen*dictSize);
      ZtZ_fft = zeros(Complex{T}, convLen*dictSize, convLen*dictSize);
      Ztx_fft = zeros(Complex{T}, convLen*dictSize);
      for j=1:numImgs;
        z_fft = reshape(fft2(reshape(z[:, j], convH, convW, dictSize)), convLen, dictSize);
        for k=1:dictSize;
          Z_fft[:, (k-1)*convLen+[1:convLen]] += diagm(z_fft[:, k]);
        end

        ZtZ_fft += Z_fft'*Z_fft;
        Ztx_fft += Z_fft'*x_fft[:, j];
      end
      s_fft = vec(fft2(s));
      lambda_s_fft = vec(fft2(lambda_s));
      d_fft = ( ZtZ_fft + mu_s.*eye(convLen*dictSize) ) \ ( Ztx_fft + mu_s.*s_fft - lambda_s_fft );
      d = real(ifft2(reshape(d_fft, convH, convW, dictSize)));
      if DEBUG;
        debugData["d"] = d;
        (obj,rerr,spar) = compute_lagrangian(debugData);
        is_mono = (prev_obj>=obj) ? 'O' : 'X';
        @printf("          D: %c obj=%.10f (rerr=%.10f, spar=%.10f)\n", is_mono, obj, rerr, spar);
        prev_obj = obj;
      end

      #
      # subproblem s
      s = d + lambda_s/mu_s;
      s[mask] = 0;
      normalize!(s);
      if DEBUG;
        debugData["s"] = s;
        (obj,rerr,spar) = compute_lagrangian(debugData);
        is_mono = (prev_obj>=obj) ? 'O' : 'X';
        @printf("          S: %c obj=%.10f (rerr=%.10f, spar=%.10f)\n", is_mono, obj, rerr, spar);
        prev_obj = obj;
      end
    end

    #
    # update lambda multipliers
    for j=1:numImgs;
      for i=1:convLen*dictSize;
        lambda_t[i, j] += mu_t*( z[i, j]-t[i, j] );
      end
    end
    for i=1:convLen*dictSize;
      lambda_s[i] += mu_s*( d[i]-s[i] );
    end
    if DEBUG;
#      debugData["lambda_t"] = lambda_t;
#      debugData["lambda_s"] = lambda_s;
#      (obj,rerr,spar) = compute_lagrangian(debugData);
#      is_mono = (prev_obj>=obj) ? 'O' : 'X';
#      @printf("          L: %c obj=%.10f (rerr=%.10f, spar=%.10f)\n", is_mono, obj, rerr, spar);
#      prev_obj = obj;
    end

    #
    # update penalties
    mu_t = min(tau*mu_t, mu_max);
    mu_s = min(tau*mu_s, mu_max);
    if DEBUG;
#      debugData["mu_t"] = mu_t;
#      debugData["mu_s"] = mu_s;
#      (obj,rerr,spar) = compute_lagrangian(debugData);
#      is_mono = (prev_obj>=obj) ? 'O' : 'X';
#      @printf("          M: %c obj=%.10f (rerr=%.10f, spar=%.10f)\n", is_mono, obj, rerr, spar);
#      prev_obj = obj;
    end

    learnConverged = FIXED_D || mean((d-s).^2)<=tol;
    codeConverged = true;
    for j=1:numImgs;
      if mu_t>0 && mean((z[:, j]-t[:, j]).^2)>tol;
        codeConverged = false;
        break
      end
    end
    if iter>=max_iter || (learnConverged && codeConverged);
      converged = true;
    end
  end
  s = s[1:filterSize, 1:filterSize, :];

  return (s, z);
end

max_stars_repo_name: bkong/FCSC


text:
#=
countDfs:
- Julia version:
- Author: Dan
- Date: 2021-11-28
=#

push!(LOAD_PATH, "/mnt/g/Python Programme/BALL.jl/src" )
using BALL
using BenchmarkTools
using StatProfilerHTML

pdb_filepath = "1AKE.pdb"

comp = System(pdb_filepath)
dfs = DataFrameSystem(comp)

@profilehtml ramachandranangles(dfs)

elapsed = @elapsed ramachandranangles(dfs)

println(elapsed)
@show @benchmark ramachandranangles(dfs)


text:
using Distributed
using SharedArrays

@everywhere import FIRRemez
@everywhere import PyPlot
@everywhere import Random
@everywhere import Printf

import FileIO
import Printf

import VisualizationTools

using FFTW
#import JLD
import BSON

fig_num = 1
PyPlot.close("all")

Random.seed!(25)

fig_num = 1

output_folder = "./output"
#L = 1000 # change this.
L = 200 # change this.

# try 1
#𝑟0 = 0.00625*2 # relative central transition band freq, in rads, but /π.

𝑟0 = 0.00625*3 # try 3

#κ = 𝑟0*0.2
κ = 𝑟0*0.5 # try 2, 3
passband_ω = BigFloat(𝑟0*π-κ) # stopband_ω*0.9
stopband_ω = BigFloat(𝑟0*π+κ)
passband_ω_Float64 = convert(Float64, passband_ω)
stopband_ω_Float64 = convert(Float64, stopband_ω)

f = xx->FIRRemez.cosinetransitionlowpassfunction(xx,passband_ω,stopband_ω)
save_name_tag = "lowpass"

println("passband is ", passband_ω)
println("stopband is ", stopband_ω)

filter_type_num = 1 # symmetric, odd length.

savefig_flag = false
plot_output_flag = false

opt_params_file_name = "../default/config_opt.txt" # location of the configuration file.

# omit opt_params_file_name is using default configuration file provided by the package.
h_BigFloat, X_BigFloat, fig_num = FIRRemez.frontendfilterdesign(filter_type_num, f,
                            L, fig_num, savefig_flag, plot_output_flag,
                            opt_params_file_name)

h = convert(Vector{Float64},h_BigFloat)
X = convert(Vector{Float64}, X_BigFloat)

Printf.@printf("Order of Chebyshev is %d, length of filter is %d, FIR type %d\n", L, length(h), filter_type_num)

band_info_string = Printf.@sprintf("passband_%f_stopband_%f", passband_ω_Float64, stopband_ω_Float64)
save_name = Printf.@sprintf("type%d_%s_L%d_%s.bson",
    filter_type_num, save_name_tag, L, band_info_string)

passband = passband_ω_Float64
stopband = stopband_ω_Float64
BSON.bson(joinpath(output_folder, save_name),
    h = h, X = X, passband = passband, stopband = stopband)
#JLD.save(joinpath(output_folder, save_name), "h" = h, "X" = X, "passband" = passband, "stopband" = stopband)
#
save_name = Printf.@sprintf("type%d_%s_L%d_%s.xlsx",
    filter_type_num, save_name_tag, L, band_info_string)
XLSX.openxlsx(joinpath(output_folder, save_name), mode="w") do xf
    sheet = xf[1]
    XLSX.rename!(sheet, "filter")

    sheet["A1", dim = 1] = ["passband"; passband]
    sheet["B1", dim = 1] = ["passband"; passband]
    sheet["C1", dim = 1] = ["Chebyshev node positions"; X]
    sheet["D1", dim = 1] = ["coefficients"; h]
end

println("Plotting")
@time fig_num = FIRRemez.plotmagnitudersp(h, fig_num, "filter's magnitude response")

# this was the output on threadripper.
# Reference set no longer changing.
# 6738.074528 seconds (1.68 G allocations: 80.051 GiB, 0.43% gc time)
# Completed mini-max algorithm. Status: Reference set no longer changing.
# Start to package filter solution.
# 1561.200300 seconds (32.29 G allocations: 1.643 TiB, 20.47% gc time)
# conversion discrepany between barycenter and Chebyshev is 8.833208284570889182364488789150169650254777059595584072429122220208126304367404e-19
# 1561.382912 seconds (32.29 G allocations: 1.644 TiB, 20.47% gc time)
# Finished packaging filter solution.
# Order of Chebyshev is 1000, length of filter is 2003, FIR type 1
# Plotting
#   1.885682 seconds (80.37 k allocations: 3.159 MiB)

max_stars_repo_name: RoyCCWang/FIRRemez

max_issues_repo_name: RoyCCWang/FIRRemez

max_forks_repo_name: RoyCCWang/FIRRemez


text:
module ULID

using Random
using Dates

export ulid

# Crockford's base-32 encoding
const _ENCODING = "0123456789ABCDEFGHJKMNPQRSTVWXYZ"
const _BASE = 32

@inline prng() = rand(RandomDevice(), UInt16) / 0xFFFF

function encodetime(now, n)
    # now should be in milliseconds-since-epoch
    s = Vector{Char}(undef, n)
    @inbounds for i = n:-1:1
        m = now % _BASE
        s[i] = _ENCODING[m+1]
        now = (now - m) ÷ _BASE
    end
    String(s)
end

function encoderandom(n)
    s = Vector{Char}(undef, n)
    @inbounds for i = n:-1:1
        r = floor(Int, _BASE * prng()) + 1
        s[i] = _ENCODING[r]
    end
    String(s)
end

"""
    ulid()

Generate a Universally Unique Lexicographically Sortable Identifier
(ULID) as a string.
"""
ulid() = encodetime(trunc(Int, datetime2unix(Dates.now()) * 1000), 10) * encoderandom(16)

end # module


text:
struct YouTube <: ShortCode
    id::String
    seekto::Int32
end

function Base.show(io::IO, ::MIME"text/plain", video::YouTube)
    video_as_text = "https://www.youtube.com/watch?v=$(video.id)&start=$(video.seekto)"
    print(io, video_as_text)
end

function Base.show(io::IO, ::MIME"text/html", video::YouTube)
    print(io, youtube(video.id, video.seekto))
end

YouTube(id, seektomin, seektosec) = YouTube(id, seektomin*60 + seektosec)
YouTube(id) = YouTube(id, 0)
"""
    Embed youtube video id that seeks seekto seconds into the video and pauses there to
    make it possible to show a particular still from the video by default. Note that
    the youtube API disallows hiding the annoying "More videos" overlay.
"""
@memoize function youtube(id, seekto)
    vid = string(uuid1()) # Assign each video unique id to allow multiple instances.
    htm = """

"""
        return htm
end

"""
    Embed youtube video by id, seek to seektomin minutes and seektosek second in
"""
function youtube(id, seektomin, seektosec)
    youtube(id, seektomin * 60 + seektosec)
end

"""
    Embed youtube video by id, uses default oembed code with reasonable size.
"""
function youtube(id)
    url = "https://youtube.com/oembed?url=http://www.youtube.com/watch?v=$id&format=json&maxwidth=600&maxheight=500"
    response = HTTP.get(url)
    json = JSON3.read(String(response.body))
    return HTML(json[:html])
end

struct Flickr <: ShortCode
    id::Integer
end

Flickr(url::String) = Flickr(parse(Int,split(url,"/")[6])) # This should be more robust

function Base.show(io::IO, ::MIME"text/plain", img::Flickr)
    json = fetch_flickr(img.id)
    flickr_as_text = "$(json[:author_name]): $(json[:title])\n$(json[:web_page])"
    print(io, flickr_as_text)
end

function Base.show(io::IO, ::MIME"text/html", img::Flickr)
    json = fetch_flickr(img.id)
    print(io, json[:html])
end

@memoize function fetch_flickr(id)
    url = "https://www.flickr.com/services/oembed/?format=json&url=http%3A//www.flickr.com/photos/bees/$id"
    response = HTTP.get(url)
    json = JSON3.read(String(response.body))
end


text:
@testset "Covariance kernels" begin
    @test expkern isa AbstractCovarianceKernel
    @test matkern isa AbstractCovarianceKernel
    @test ar1kern isa AbstractCovarianceKernel

    @test expΣ isa PDMat
    @test matΣ isa PDMat
    @test ar1Σ isa PDMat

    n = length(Ω)
    @test size(expΣ) == (n, n)
    @test size(matΣ) == (n, n)
    @test all(diag(expΣ.mat) .== σ²)
    @test all(diag(matΣ.mat) .== σ²)
end


text:
include("AbstractCausalAnalysis.jl")
include("CausalAnalysis.jl")
include("summarise.jl")
include("recipes.jl")

export AbstractCausalAnalysis, CausalAnalysis, causality, summarise


text:
# # [SIR example with age structure](@id sir_age)
# Sean L. Wu (@slwu89), 2021-1-22

using Individual.Sampling
using Individual.SchemaBase

using Catlab.Present, Catlab.CSetDataStructures, Catlab.Theories, Catlab.CategoricalAlgebra, Catlab.Graphics, Catlab.Graphs
using Plots, GraphViz
using LinearAlgebra
using Random

# ## Introduction
# This tutorial shows how to extend the existing abstract ACSet types to incorporate an additional attribute for modeling
# age structure. The model is an age-structured SIR model where the force of infection experienced by an individual is the
# sum of their effective contacts with each age group. The contact matrix ``C_{ij}`` gives the number of daily contacts a person in
# age group ``i`` has with persons in age group ``j``.

# ## Schema

# We need to define a new schema inheriting from `TheoryIBM`. Please note that when we use `@acset_type` to define the type
# for the ACSet data structure for our model, we need to set `:statelabel` as a unique index. This is because the simulation
# framework assumes that there is a one to one mapping between states and their unique labels. Please look at the type definition
# for `IBM` and `SchedulingIBM` in the source files to see what indices are expected, and the [ACSet documentation in Catlab.jl](https://algebraicjulia.github.io/Catlab.jl/stable/apis/categorical_algebra/#Acsets)
# for more information about these data types.

@present TheoryAgeIBM <: TheoryIBM begin
    Age::AttrType
    age::Attr(Person, Age)
end

@abstract_acset_type AbstractAgeIBM <: AbstractIBM
@acset_type AgeIBM(TheoryAgeIBM, index = [:state, :state_update, :age], unique_index = [:statelabel]) <: AbstractAgeIBM

# The schema looks like this. Note that there is an additional attribute in the schema, for age of each individual.

to_graphviz(TheoryAgeIBM)

# ## Parameters

# To parameterize the model from ``R_{0}`` and ``\gamma``, we follow the method from ["Social contact patterns and control strategies for influenza in the elderly"](http://www.sherrytowers.com/towers_feng_2012.pdf).
# where ``R_{0}`` is ``\beta / \gamma`` multiplied by the largest eigenvalue of the reciprocal (symmetric) contact matrix. Note that the
# interpretation of ``\beta`` is slightly different in the age-structured model; here it is only the probability of infection given infectious contact,
# in the other examples it was the rate of contact multiplied by the probability of infection given infectious contact.

# We use the estimated contact matrix for Taiwan from ["Projecting social contact matrices in 152 countries using contact surveys and demographic data"](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005697), given in `contact_unnorm`.
# There are 16 age groups, in 5 year bins, up to age 75, at which all remaining persons are included in the last bin.
# Because the contact matrix is directly estimated from survey data, it is not reciprocal (i.e. the number of contacts sent from
# age group ``i`` to ``j`` must equal the number from ``j`` to ``i``). We use a simple pairwise correction to generate a reciprocal contact matrix.
# For more information on reciprocity in parameterizing contact matrices, please see the article ["Projecting social contact matrices to different demographic structures"](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1006638).

# The vector `pop_TW_actual` gives the census population in 16 bins, which we scale to our population `N`.

N = 1000
I0 = 8
S0 = N - I0
Δt = 0.1
tmax = 100
steps = Int(tmax/Δt)
γ = 1/10 # recovery rate
R0 = 2.5

pop_TW_actual = [1011137, 1041749, 975803, 1213008, 1502279, 1618075, 1621625, 1926961, 2032452, 1755391, 1806638, 1852580, 1684376, 1416638, 903130, 1454933]
pop_TW = N * (pop_TW_actual / sum(pop_TW_actual))
pop_TW = Int.(round.(pop_TW))
N_ages = Float64.(pop_TW)

contact_unnorm = [1.13133995414878	0.557101193757009	0.277403588569526	0.166283484217312	0.243637116538954	0.501429856885401	0.773456505042611	0.664734678057495	0.373030323264845	0.207403718009079	0.219768245286902	0.183713318566528	0.100598824431586	0.0723635363224441	0.0451987854057666	0.0264190015627899
0.522116359073187	5.13993582380529	0.999881750990576	0.256872404884024	0.152965888814777	0.440037571740546	0.81718457664922	0.901816339582795	0.776770796127256	0.324978689513501	0.203946781656598	0.170308924697168	0.114410112258284	0.0766043987718478	0.0377257124824009	0.0309979523804595
0.233720556183555	1.92985936809443	11.0692392115508	0.983950046684279	0.331933259964229	0.36113787699594	0.563341231085396	0.903032808933502	1.17790738474796	0.651636521136264	0.342842258726613	0.173034004951767	0.085317189948624	0.0792891389004607	0.0627981424580935	0.0524757101169723
0.117682987086876	0.345910173819922	3.30103123371644	10.0791939893023	1.50916450884044	0.791285884658043	0.640497694236322	0.824953404104752	1.08583985409548	1.09515652926247	0.535792916941561	0.197364288066501	0.0713874830715285	0.0516466911367817	0.0315463617042823	0.0232625811949096
0.184283846973784	0.180619839511097	0.291480297355471	2.42727014100382	4.36737726042213	2.04934997728767	1.39723781506946	1.20032801881538	1.01070904475343	1.23524796447696	0.800832127889001	0.39902019518293	0.0944429030797382	0.0416382822716035	0.0470826453808244	0.0394951197797532
0.56156192817194	0.342992537769556	0.200122810772951	0.916989232128664	2.3980129020946	4.30977975618343	2.22681401915761	1.59261033191491	1.30238635661337	1.05207634128304	1.00949814315475	0.536928086612407	0.168547286689823	0.0624501078266029	0.0317257423796732	0.0269431049949397
0.817794867047499	1.1441520682733	0.841196707142586	0.551555851206977	1.22865834990464	2.17367993548303	3.72979423098863	2.14929466728635	1.55306730764128	1.11703404262428	0.859302034849221	0.601599257124208	0.256942068421412	0.106203837431624	0.0555526815130332	0.0511357782774989
0.665898107120277	1.15314199919716	0.9740043186808	0.697233734859189	0.825555277703876	1.54165452016297	1.99243882433986	3.25956559702894	2.06558139725526	1.22850364764608	0.791670106907652	0.438605150114473	0.244794735730732	0.148777023280774	0.0869686510606626	0.0349673148163165
0.344192143343282	0.807995410469339	1.14415043227526	1.22438695056131	0.988030619038175	1.30424464223099	1.72028688625099	1.87814450829015	2.90470388452872	1.49758830410338	0.934324050254837	0.330600903622231	0.182537686520201	0.114580140863156	0.0813846406012246	0.0400784886603155
0.235791530653026	0.585194161111923	0.776189668181761	1.67596626847397	0.98479644076933	1.01459833535067	1.25482711892082	1.36074838781177	1.44439037412333	2.06915205205562	0.978401494457824	0.39999986541145	0.147110687316517	0.0835852658511084	0.0775170890988546	0.0761516911425131
0.247247471139332	0.643310594857866	1.05811260803268	1.31212557944165	1.09555967723574	1.35073741525286	1.16660463514035	1.06882688561953	1.37118566862	1.49192984016939	1.73933382380217	0.740193371293642	0.243542236335681	0.109266882350377	0.0798121160310869	0.0833033422190585
0.442684358021657	0.702087727411373	0.662177052932204	0.816859526244507	0.738982978496004	1.2358668685183	1.27567127460055	0.896236373088294	0.950326986928672	0.797183567958301	1.02749599715298	1.51681655407295	0.451707898396268	0.195099577202936	0.0961333966552319	0.0769572529007728
0.327638342195265	0.328285016941782	0.238165901349415	0.354821937305308	0.320966945158153	0.529029787102833	0.677350792373831	0.650949037046765	0.501900491233428	0.407240319726415	0.395920768113556	0.541100083728198	0.836769891969799	0.267847588746602	0.145638827229407	0.059691205799689
0.193568331575816	0.313138348784034	0.2640962452464	0.160094097408356	0.204553272242237	0.325505767060564	0.516407578574901	0.473341415992634	0.438186035505819	0.281001067157269	0.296302451590512	0.375946282497057	0.331749859023468	0.679792406234649	0.173000119847818	0.0702778253123571
0.0963863876721027	0.282393513167451	0.295599962015119	0.30895124164141	0.1384465024105	0.257953175229946	0.27963826844241	0.436956466193058	0.509042820432503	0.403187215044375	0.316140357668137	0.277393058839257	0.385660305905714	0.370945733661385	0.606083091791185	0.185944796718132
0.188903590501605	0.259894019899239	0.37803587788717	0.31873977192663	0.125868047138266	0.160125548829692	0.269195296463481	0.305611946209338	0.348244256184042	0.390220623327373	0.393237304756271	0.240097073151863	0.144076910885202	0.19701354194507	0.184345002762545	0.325239463933643];

contact = (contact_unnorm + transpose(contact_unnorm) .* ( pop_TW * transpose(1 ./ pop_TW))) ./ 2

C = contact
for i = 1:16
    for j = 1:16
        C[i,j] = contact[i,j]*N_ages[i]/N_ages[j]
    end
end

β = R0 * γ / max(eigvals(C)...);

# The initial states are the same as the other tutorials, 1000 individuals, 8 of whom are intially infected.

initial_states = fill("S", N)
initial_states[rand(1:N, I0)] .= "I"
state_labels = ["S", "I", "R"];

# ## Model object

# The `AgeIBM` schema we made needs two type parameters.
# The first `String` is for the set of state labels, inherited from `IBM`.
# The second is for the age attribute.
# Because age is discretized into bins to match the survey data used to parameterize the contact matrix, the attribute type is an integer.

const SIR = AgeIBM{String, Int64}()
initialize_states(SIR, initial_states, state_labels);

# We need to sample the age bins for each person, such that the sizes of each age bin are correct.

ages = vcat([fill(i, pop_TW[i]) for i in 1:16]...)
shuffle!(ages)
set_subpart!(SIR, 1:N, :age, ages);

# ## Processes

# The force of infection on a person in age class ``i`` is computed according to:
# ```math
#  \lambda_{i} = \beta \sum\limits_{j} C_{ij} \left( \frac{I_{j}}{N_{j}} \right)
# ```
# We use a helper function `tabulate_ages` to calculate the infectious population sizes ``I_{j}``.

function tabulate_ages(ages_vec)
    [Float64(sum(ages_vec .== i)) for i = 1:16]
end

function infection_process(t::Int)
    I = get_index_state(SIR, "I")
    I_ages = subpart(SIR, I, :age)
    I_ages = tabulate_ages(I_ages)

    λ = β * sum(C * diagm(I_ages ./ N_ages), dims = 2)

    S = get_index_state(SIR, "S")
    S_ages = subpart(SIR, S, :age)

    λ = λ[S_ages]

    S = bernoulli_sample(S, λ, Δt)
    queue_state_update(SIR, S, "I")
end

# The recovery process is the same as in the [basic SIR tutorial](@ref sir_basic).

function recovery_process(t::Int)
    I = get_index_state(SIR, "I")
    I = bernoulli_sample(I, γ, Δt)
    queue_state_update(SIR, I, "R")
end

# ## Simulation

# We use `render_process` to create a rendering (output) process and
# a matrix giving state counts by time step.
# Then we pass the necessary processes to `simulation_loop` to draw a trajectory and plot the results.

state_out, render_process = render_states(SIR, steps)

simulation_loop(SIR, [infection_process, recovery_process, render_process], steps)

plot(
    (1:steps) * Δt,
    state_out,
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number"
)

# We can also plot the size of the epidemic in each age bin, using a stacked bar chart.

ever_infected = vcat(incident(SIR, [2, 3], :state)...)
ever_infected_ages = subpart(SIR, ever_infected, :age)
ever_infected_ages = [sum(ever_infected_ages .== i) for i = 1:16]

bar(pop_TW, label = false, xlabel = "Age bins", ylabel = "Number")
bar!(ever_infected_ages, label = false)


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

JLLWrappers.@generate_main_file_header("RTKLib")
JLLWrappers.@generate_main_file("RTKLib", UUID("d19eca65-acdf-50fd-a11a-2de6aab1dc1f"))
end  # module RTKLib_jll


text:
using Pkg
Pkg.rm("AzureClusterlessHPC")
Pkg.develop(url="https://github.com/microsoft/AzureClusterlessHPC.jl.git")
Pkg.rm("JUDI4Cloud")
Pkg.develop(url="https://github.com/slimgroup/JUDI4Cloud.jl.git")


text:
function mul_analytic!(τ::Vector{<:AbstractFloat},
                       dx1::AbstractFloat,
                       ie1::Integer,
                       xs1::Integer)
    for i1 = 1:ie1
        @inbounds τ[i1] *= sqrt(sum(((i1-xs1)*dx1)^2))
    end
end

function mul_analytic!(τ::Matrix{<:AbstractFloat},
                       dx1::AbstractFloat, dx2::AbstractFloat,
                       ie1::Integer, ie2::Integer,
                       xs1::Integer, xs2::Integer)
    for i2 = 1:ie2
        for i1 = 1:ie1
            @inbounds τ[i1,i2] *= sqrt(sum(((i1-xs1)*dx1)^2+((i2-xs2)*dx2)^2))
        end
    end
end

function mul_analytic!(τ::Array{<:AbstractFloat,3},
                       dx1::AbstractFloat, dx2::AbstractFloat, dx3::AbstractFloat,
                       ie1::Integer, ie2::Integer, ie3::Integer,
                       xs1::Integer, xs2::Integer, xs3::Integer)
    for i3 = 1:ie3
        for i2 = 1:ie2
            for i1 = 1:ie1
                @inbounds τ[i1,i2,i3] *= sqrt(sum(((i1-xs1)*dx1)^2+((i2-xs2)*dx2)^2+((i3-xs3)*dx3)^2))
            end
        end
    end
end

function grad_analytic(τ0::Vector{R},
                       dx1::AbstractFloat,
                       ie1::Integer,
                       xs1::Integer) where {R <: AbstractFloat}
    ∇τ0 = [ones(R, size(τ0))]
    for i1 = 1:ie1
        @inbounds ∇τ0[1][i1] = dx1*(i1-xs1)/τ0[i1]
    end
    #set terms at source
    ∇τ0[1][xs1] = 1
    ∇τ0
end

function grad_analytic(τ0::Matrix{R},
                       dx1::AbstractFloat, dx2::AbstractFloat,
                       ie1::Integer, ie2::Integer,
                       xs1::Integer, xs2::Integer) where {R <: AbstractFloat}
    ∇τ0 = [ones(R, size(τ0)), ones(R, size(τ0))]
    for i2 = 1:ie2
        for i1 = 1:ie1
            @inbounds ∇τ0[1][i1,i2] = dx1*(i1-xs1)/τ0[i1,i2]
            @inbounds ∇τ0[2][i1,i2] = dx2*(i2-xs2)/τ0[i1,i2]
        end
    end
    ∇τ0[1][xs1, xs2] = 1/sqrt(2)
    ∇τ0[2][xs1, xs2] = 1/sqrt(2)
    ∇τ0
end

function grad_analytic(τ0::Array{R,3},
                       dx1::AbstractFloat, dx2::AbstractFloat, dx3::AbstractFloat,
                       ie1::Integer, ie2::Integer, ie3::Integer,
                       xs1::Integer, xs2::Integer, xs3::Integer) where {R <: AbstractFloat}
    ∇τ0 = [ones(R, size(τ0)), ones(R, size(τ0)), ones(R, size(τ0))]
    for i3 = 1:ie3
        for i2 = 1:ie2
            for i1 = 1:ie1
                @inbounds ∇τ0[1][i1,i2,i3] = dx1*(i1-xs1)/τ0[i1,i2,i3]
                @inbounds ∇τ0[2][i1,i2,i3] = dx2*(i2-xs2)/τ0[i1,i2,i3]
                @inbounds ∇τ0[3][i1,i2,i3] = dx3*(i3-xs3)/τ0[i1,i2,i3]
            end
        end
    end
    ∇τ0[1][xs1, xs2, xs3] = 1/sqrt(3)
    ∇τ0[2][xs1, xs2, xs3] = 1/sqrt(3)
    ∇τ0[3][xs1, xs2, xs3] = 1/sqrt(3)
    ∇τ0
end


text:
using MPI
using Plots
MPI.Init()

comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
size = MPI.Comm_size(comm)

dst = mod(rank+1, size)
src = mod(rank-1, size)

bw  = Array{Float64}(undef, 12)
print("$rank: Generating Data\n")
for n in 1:12
    N = 2^n
    send_msg = Array{Int8}(undef,N)
    recv_msg = Array{Int8}(undef,N)
    fill!(send_msg, Int8(rank))
    rreq = MPI.Irecv!(recv_msg, src, src+32, comm)
    # print("$rank: Sending   $rank -> $dst\n")
    t  = @elapsed (MPI.Send(send_msg, dst, rank+32, comm))
    stats = MPI.Waitall!([rreq])
    # print("$rank: Received $src -> $rank\n")
    val = (n/t)*10e-6
    bw[n] = val
end
print("$rank: Data Generated, Starting plot\n")
# print("$rank: $bw")

msize = [2^n for n in 1:12]

plot(msize, bw, xlabel="Message Size(bytes)", ylabel="Bandwith(MB/s)", title="Process $rank", size=(600, 600), label=false)

png("D:\\JuliaDev\\Scripts\\Plots\\mpass$rank")

MPI.Barrier(comm)


text:
#Examples:

"""
    cbioproc_example1()

A 3D example without land mask or MeshArrays.
"""
function cbioproc_example1()
    cbioproc_example1("")
end

function cbioproc_example1(dirIn::String)
    in=fill(1.,(90,1170,50))
    SPM,lon,lat=read_SPM(dirIn)
    siz=size(lon)
    out=MatrixInterp(in,SPM,siz)
    #
    outLoop,outName=loop_exampleB(1:13,SPM,siz)
    tmp=transpose(outLoop[1,1])
    heatmap(vec(lon[:,1]),vec(lat[1,:]),tmp)
    return outLoop,outName
end

"""
    cbioproc_example2()

A 2D example with MeshArrays + mask or loop.
"""
function cbioproc_example2()
    dirIn=""
    cbioproc_example2(dirIn)
end

function cbioproc_example2(dirIn::String)
    #
    GCMGridSpec()
    GCMGridLoad()
    msk2d=mask(view(MeshArrays.hFacC,:,:,1),NaN,0)
    SPM,lon,lat=read_SPM(dirIn)
    siz=size(lon)
    #
    in=read_bin(dirIn*"diags/state_2d_set1.0000000732.data",Float32)
    in=convert2gcmfaces(in[:,:,1].*msk2d);
    outMsk=MatrixInterp(in,SPM,siz)
    heatmap(vec(lon[:,1]),vec(lat[1,:]),transpose(outMsk))
    #
    outLoop=loop_exampleA(1:13,SPM,siz);
    heatmap(vec(lon[:,1]),vec(lat[1,:]),transpose(outLoop[1]))
    #
    return outLoop,outMsk
end

"""
    loop_exampleA(indx,SPM,siz)

Applies `MatrixInterp` in a loop over a subset of model output files (`filList[indx]`)
"""
function loop_exampleA(indx,SPM,siz)
   dirIn="diags/"
   filIn="state_2d_set1"
   tmp1=readdir(dirIn)
   tmp1=filter(x -> occursin(filIn,x),tmp1)
   filList=filter(x -> occursin(".data",x),tmp1)
   filList=filList[indx]
   nf=length(filList)
   #show(filList)

   out=Array{Any,1}(undef,nf)
   for ff=1:nf
      fil=dirIn*filList[ff]
      fid = open(fil);
      #println(fil)
      fld = Array{Float32,2}(undef,(90,1170));
      #fld = Array{prec,1}(undef,(n1*n2*n3));
      read!(fid,fld);
      fld = hton.(fld);
      #println(size(fld))
      out[ff]=MatrixInterp(fld,SPM,siz)
   end

   return out
end

"""
    loop_exampleB(indx,SPM,siz)

Applies `MatrixInterp` in a loop over a subset of model output files (`filList[indx]`)
"""
function loop_exampleB(indx,SPM,siz)
   dirIn="diags/"
   filIn="state_2d_set1"
   tmp1=readdir(dirIn)
   tmp1=filter(x -> occursin(filIn,x),tmp1)
   filList=filter(x -> occursin(".data",x),tmp1)
   filList=filList[indx]
   nf=length(filList)
   #show(filList)

   MetaFile=dirIn*filList[1]
   MetaFile=MetaFile[1:end-5]*".meta"
   MetaFile=MetaFileRead(MetaFile)
   nv=length(MetaFile["fldList"])
   nd=MetaFile["nDims"]
   dims=Int.(MetaFile["dimList"][:,1])
   prec=MetaFile["dataprec"]

   out=Array{Any,2}(undef,(nf,nv))
   tmp = Array{prec,nd}(undef,Tuple(dims));
   for ff=1:nf
      fil=dirIn*filList[ff]
      fid = open(fil);
      #println(fil)
      #tmp = Array{prec,1}(undef,(n1*n2*n3));
      for vv=1:nv
         read!(fid,tmp);
         tmp = hton.(tmp);
         out[ff,vv]=MatrixInterp(tmp,SPM,siz)
      end
   end

   return out,MetaFile["fldList"]
end


text:
"""
    struct LaguerreBasis{P} <: AbstractMultipleOrthogonalBasis{P}
        polynomials::Vector{P}
    end

Orthogonal polynomial with respect to the univariate weight function ``w(x) = \\exp(-x)`` over the interval ``[0, \\infty]``.
"""
struct LaguerreBasis{P} <: AbstractMultipleOrthogonalBasis{P}
    polynomials::Vector{P}
end

polynomial_type(::Type{<:LaguerreBasis}, V::Type) = MP.polynomialtype(V, Float64)

even_odd_separated(::Type{<:LaguerreBasis}) = false

reccurence_first_coef(::Type{<:LaguerreBasis}, degree) = -1
reccurence_second_coef(::Type{<:LaguerreBasis}, degree) = (2degree - 1)
reccurence_third_coef(::Type{<:LaguerreBasis}, degree) = -(degree - 1)
reccurence_deno_coef(::Type{<:LaguerreBasis}, degree) = degree

degree_one_univariate_polynomial(::Type{<:LaguerreBasis}, variable::MP.AbstractVariable) = MA.@rewrite(1 - variable)

max_stars_repo_name: JuliaAlgebra/PolynomialBases

max_issues_repo_name: JuliaAlgebra/PolynomialBases

max_forks_repo_name: JuliaAlgebra/PolynomialBases


text:
using ProjectMoana
using ColorTypes

world = World()
add_shape(world, Sphere(translation(Vec(2., 0., 0.)) * scaling(Vec(0.1, 0.1, 0.1))))

image = HdrImage(3, 3)
camera = PerspectiveCamera()
tracer = ImageTracer(image, camera)
renderer = OnOff_Renderer(world)

fire_all_rays(tracer, OnOff, renderer)

@testset "Renderer: OnOff renderer" begin
    for i ∈ 1:3
        for j ∈ 1:3
            image.pixels[get_pixel(image, i, j)] == RGB(1.f0,1.f0,1.f0)
        end
    end

    @test image.pixels[get_pixel(image,1, 1)]≈RGB(0.f0,0.f0,0.f0)
    @test image.pixels[get_pixel(image,2, 1)]≈RGB(0.f0,0.f0,0.f0)
    @test image.pixels[get_pixel(image,3, 1)]≈RGB(0.f0,0.f0,0.f0)

    @test image.pixels[get_pixel(image,1, 2)]≈RGB(0.f0,0.f0,0.f0)
    @test image.pixels[get_pixel(image,2, 2)]≈RGB(1.f0,1.f0,1.f0)
    @test image.pixels[get_pixel(image,3, 2)]≈RGB(0.f0,0.f0,0.f0)

    @test image.pixels[get_pixel(image,1, 3)]≈RGB(0.f0,0.f0,0.f0)
    @test image.pixels[get_pixel(image,2, 3)]≈RGB(0.f0,0.f0,0.f0)
    @test image.pixels[get_pixel(image,3, 3)]≈RGB(0.f0,0.f0,0.f0)
end

world = World()
sphere_color = RGB(1.f0, 2.0f0, 3.0f0)
sphere =Sphere(translation(Vec(2, 0, 0)) * scaling(Vec(0.2, 0.2, 0.2)), Material(DiffuseBRDF(UniformPigment(sphere_color))))
add_shape(world, sphere)

image = HdrImage(3, 3)
camera = OrthogonalCamera()
tracer = ImageTracer(image, camera)
renderer = Flat_Renderer(world)

fire_all_rays(tracer, Flat, renderer)

@testset "Renderer: Flat renderer" begin
    @test image.pixels[get_pixel(image,1, 1)]≈RGB(0.f0,0.f0,0.f0)
    @test image.pixels[get_pixel(image,2, 1)]≈RGB(0.f0,0.f0,0.f0)
    @test image.pixels[get_pixel(image,3, 1)]≈RGB(0.f0,0.f0,0.f0)

    @test image.pixels[get_pixel(image,1, 2)]≈RGB(0.f0,0.f0,0.f0)
    @test image.pixels[get_pixel(image,2, 2)]≈sphere_color
    @test image.pixels[get_pixel(image,3, 2)]≈RGB(0.f0,0.f0,0.f0)

    @test image.pixels[get_pixel(image,1, 3)]≈RGB(0.f0,0.f0,0.f0)
    @test image.pixels[get_pixel(image,2, 3)]≈RGB(0.f0,0.f0,0.f0)
    @test image.pixels[get_pixel(image,3, 3)]≈RGB(0.f0,0.f0,0.f0)
end

@testset "Renderer: Furnace test" begin
    pcg = PCG()

    for i in 1:10^3
        world = World()

        emitted_radiance = Float32(pcg_randf(pcg) * 0.8) #per assicurare la convergenza scelgo un numero casuale non troppo vicino a 1
        reflectance = Float32(pcg_randf(pcg) * 0.8)
        enclosure_material = Material( DiffuseBRDF(UniformPigment(RGB(1.f0, 1.f0, 1.f0) * reflectance)), UniformPigment(RGB(1.f0, 1.f0, 1.f0) * emitted_radiance))

        add_shape(world, Sphere(Transformation(), enclosure_material))

        path_tracer = PathTracer_Renderer(world, pcg=pcg, num_of_rays=1, max_depth=1000, russian_roulette_limit=101)

        ray = Ray(Point(0, 0, 0), Vec(1, 0, 0))
        color = PathTracer(ray, path_tracer)

        expected = emitted_radiance / (1.f0 - reflectance)
        @test isapprox(expected, color.r, atol=5e-2)
        @test isapprox(expected, color.g, atol=5e-2)
        @test isapprox(expected, color.b, atol=5e-2)
    end

end