text:
# bias_variance-arwmh-pseduohuber
#
# This script compares the Stein discrepancy of sample sequences drawn
# from asymptotically biased approximate random walk Metropolis-Hastings
# (ARWMH) chains (http://arxiv.org/pdf/1304.5299v4.pdf) and asymptotically
# exact RWMH chains.  By plotting the Stein discrepancy results against
# a measure of computation required per sample (e.g., the number of
# datapoint likelihood evaluations needed), we can quantify a bias-variance
# trade-off underlying the approximate MCMC scheme.
#
# The target distribution is a Bayesian huber loss regression posterior
# under independent Gaussian priors based on Gelman's Radon dataset.
# The regression tries to model the relationship between a few exploratory
# variables and the levels of Radon present in houses.

include("experiment_utils.jl")

# Set random number generator seed
@parseintcli seed "s" "seed" 7
# set the epsilon param
@parsestringcli epsiloninput "e" "epsilon" "1e-5"
epsilon = float(epsiloninput)

using DataFrames

using SteinDiscrepancy: gsd
using SteinDistributions: SteinHuberLossRegressionGaussianPrior, gradlogdensity
using SteinSamplers: runapproxmh

# Select a solver for Stein discrepancy optimization problem
#solver = "clp"
solver = "gurobi"
# Setup the data
data_set = "radon"
rdat = readtable("src/experiments/radon/srrs2.dat")
cdat = readtable("src/experiments/radon/cty.dat")
rename!(cdat, :ctfips, :cntyfips)

dat = join(rdat, cdat, on = [:stfips, :cntyfips])
# subset only to Minnesota! (Gelman does this)
dat = dat[dat[:st] .== "MN",:]
# log radon
dat[:logradon] = log(dat[:activity])
dat[dat[:activity] .== 0, :logradon] = log(0.1)
# log uranium
dat[:loguranium] = log(dat[:Uppm])
dat[dat[:Uppm] .== 0, :loguranium] = log(0.1)
# make up X and y
X = Array(Float64, nrow(dat), 3)
X[:,1] = 1.0
X[:,2] = dat[:floor]
X[:,3] = dat[:loguranium]
y = convert(Array, dat[:logradon])
# Parameter dimension
d = size(X,2)
# the max number of samples to run in the solver
n = 1000
# Associate sample weights
q = ones(n, 1) ./ n
# name of prior + posterior
distname = "huberloss-gaussianprior"
# setup distribution
target = SteinHuberLossRegressionGaussianPrior(X, y)
# define gradlogp
function gradlogp(x::Array{Float64,1})
    gradlogdensity(target, x)
end
# intiial guess is OLS estimate
beta0 = inv(X' * X) * (X' * y)
# burnin likelihoods
burninlikelihoods = 1e5
# max number of likelihoods
maxlikelihoods = 1e7
# batchsize for different algorithms [ataset is ~1000 rows]
batchsize = 5
# Sample sizes at which Stein discrepancy will be computed
ns = vcat(10:10:min(400,n), 420:20:min(600,n), 650:50:min(1000,n))

@printf("Generating samples for epsilon=%f\n", epsilon)

(betas, numlikelihoods) =
    @setseed seed runapproxmh(target, beta0;
                              epsilon=epsilon,
                              numlikelihood=maxlikelihoods,
                              batchsize=batchsize,
                              proposalvariance=1e-4)

# Thin to n samples
firstn = find([(l >= burninlikelihoods) for l in numlikelihoods])[1]
lastn = length(numlikelihoods)
idx = round(Int, linspace(firstn, lastn, min(lastn - firstn + 1, n)))

# Keep track of yhats
yhat = X * betas[idx,:]'
scaledyhat = yhat ./ vec(maximum(X, 2))

for i in ns
    println("Beginning optimization for n=$(i), d=$(d), epsilon=$(epsilon)")
    # Compute Stein discrepancy for first i points
    betai = betas[idx[1:i],:]
    res = gsd(points=betai, gradlogdensity=gradlogp, solver=solver)
    println("\tn = $(i), objective = $(res.objectivevalue)")

    instance_data = Dict{Any, Any}(
        "n" => i,
        "distname" => distname,
        "d" => d,
        "seed" => seed,
        "epsilon" => epsilon,
        "batchsize" => batchsize,
        "X" => res.points,
        "q" => res.weights,
        "xnorm" => res.xnorm,
        "g" => res.g,
        "gradg" => res.gradg,
        "objectivevalue" => vec(res.objectivevalue),
        "scaled.yhat" => vec(mean(scaledyhat[:,1:i], 2)),
        "numlikelihood" => numlikelihoods[idx[i]],
        "edgetime" => res.edgetime,
        "solvetime" => res.solvetime,
    )

    save_json(
        instance_data;
        dir="bias_variance-arwmh-huberloss",
        distname=distname,
        dataset=data_set,
        n=i,
        burninlikelihoods=int(burninlikelihoods),
        epsilon=epsilon,
        batchsize=batchsize,
        d=d,
        seed=seed
    )
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    GraphMol, GraphMolReaction,
    graphmol, remapnodes, todict, tojson,
    SDFile, SMILES,
    getatom, getbond, hasbond,
    setatom!, setbond!, addatom!, addbond!,
    atomcount, bondcount

using JSON
import Unmarshal

struct GraphMol{A<:Atom,B<:Bond} <: OrderedGraph
    neighbormap::Vector{Dict{Int,Int}}
    edges::Vector{Tuple{Int,Int}}
    nodeattrs::Vector{A}
    edgeattrs::Vector{B}
    cache::Dict{Symbol,Any}
    attributes::Dict{Symbol,Any}
end

GraphMol{A, B}() where {A<:Atom,B<:Bond} = GraphMol{A,B}([], [], [], [], Dict(), Dict())

struct GraphMolReaction
    reactants::Vector{GraphMol}
    products::Vector{GraphMol}
end

GraphMolReaction() = GraphMolReaction(GraphMol[], GraphMol[])

"""
    graphmol() -> GraphMol

Generate empty `GraphMol`.
"""
graphmol(::Type{A}, ::Type{B}
    ) where {A<:Atom,B<:Bond} = GraphMol{A,B}([], [], [], [], Dict(), Dict())

"""
    graphmol(atoms::Vector{Atom}, bonds::Vector{Bond}) -> GraphMol

Generate `GraphMol` that has the given atom objects and edge objects.
"""
function graphmol(
        edges, atoms::Vector{A}, bonds::Vector{B}) where {A<:Atom,B<:Bond}
    nbrmap = [Dict{Int,Int}() for i in 1:length(atoms)]
    edges = collect(edges)
    for (i, (u, v)) in enumerate(edges)
        nbrmap[u][i] = v
        nbrmap[v][i] = u
    end
    return GraphMol(
        nbrmap, edges, atoms, bonds, Dict{Symbol,Any}(), Dict{Symbol,Any}())
end

"""
    graphmol(mol::GraphMol; clearcache=true) -> GraphMol

Copy `GraphMol`. Cached properties will be removed.
"""
function graphmol(graph::GraphMol{A,B}; clearcache=true
        ) where {A<:Atom,B<:Bond}
    newg = graphmol(A, B)
    for nbr in graph.neighbormap
        push!(newg.neighbormap, copy(nbr))
    end
    append!(newg.edges, graph.edges)
    append!(newg.nodeattrs, graph.nodeattrs)
    append!(newg.edgeattrs, graph.edgeattrs)
    clearcache || merge!(newg.cache, graph.cache)
    merge!(newg.attributes, graph.attributes)
    return newg
end

Graph.clone(mol::GraphMol) = graphmol(mol, clearcache=false)

"""
    graphmol(mol::SubgraphView{GraphMol}) -> GraphMol

Generate a new `GraphMol` from a substructure view.

Graph property caches and attributes are not inherited.
"""
function graphmol(view::SubgraphView)
    newg = graphmol(nodeattrtype(view), edgeattrtype(view))
    nkeys = sort(collect(nodeset(view)))
    ekeys = sort(collect(edgeset(view)))
    nmap = Dict{Int,Int}()
    for (i, n) in enumerate(nkeys)
        nmap[n] = i
        push!(newg.nodeattrs, nodeattr(view, n))
        push!(newg.neighbormap, Dict())
    end
    for (i, e) in enumerate(ekeys)
        (oldu, oldv) = getedge(view, e)
        u = nmap[oldu]
        v = nmap[oldv]
        push!(newg.edges, (u, v))
        push!(newg.edgeattrs, edgeattr(view, e))
        newg.neighbormap[u][i] = v
        newg.neighbormap[v][i] = u
    end
    return newg
end

function remapnodes(graph::GraphMol, mapping::Dict{Int,Int})
    newg = graphmol(graph)
    for i in 1:nodecount(graph)
        newg.nodeattrs[mapping[i]] = graph.nodeattrs[i]
        newg.neighbormap[mapping[i]] = Dict(
            k => mapping[v] for (k, v) in graph.neighbormap[i])
    end
    for i in 1:edgecount(graph)
        u, v = graph.edges[i]
        newg.edges[i] = (mapping[u], mapping[v])
        newg.edgeattrs[i] = graph.edgeattrs[i]
    end
    merge!(newg.attributes, graph.attributes)
    return newg
end

"""
    todict(graph::GraphMol) -> Dict{String,Any}

Convert molecule object into JSON compatible dictionary.
"""
todict(graph::GraphMol) = Dict(
    "edges" => graph.edges,
    "nodetype" => string(nodeattrtype(graph)),
    "edgetype" => string(edgeattrtype(graph)),
    "nodeattrs" => [todict(atom) for atom in graph.nodeattrs],
    "edgeattrs" => [todict(bond) for bond in graph.edgeattrs],
    "cache" => Dict(
        string(k) => (string(typeof(v)), v) for (k, v) in graph.cache),
    "attributes" => Dict(string(k) => v for (k, v) in graph.attributes)
)

tojson(graph::GraphMol) = JSON.json(todict(graph))

# Note: Unmarshal patch for Set
function Unmarshal.unmarshal(::Type{Set{E}},
    parsedJson::Union{Vector, AbstractArray},
    verbose::Bool=false, verboseLvl::Int=0) where E
    if verbose
        prettyPrint(verboseLvl, "Set{$E}")
        verboseLvl += 1
    end
    Set(Unmarshal.unmarshal(
        E, field, verbose, verboseLvl) for field in parsedJson)
end

function graphmol(data::Dict{String,Any})
    atomtype = eval(Meta.parse(data["nodetype"]))
    bondtype = eval(Meta.parse(data["edgetype"]))
    mol = graphmol(atomtype, bondtype)
    for a in data["nodeattrs"]
        push!(mol.nodeattrs, atomtype(a))
        push!(mol.neighbormap, Dict())
    end
    for (i, (u, v)) in enumerate(data["edges"])
        mol.neighbormap[u][i] = v
        mol.neighbormap[v][i] = u
        push!(mol.edgeattrs, bondtype(data["edgeattrs"][i]))
        push!(mol.edges, (u, v))
    end
    for (k, (vt, v)) in data["cache"]
        mol.cache[Symbol(k)] = Unmarshal.unmarshal(eval(Meta.parse(vt)), v)
    end
    for (k, v) in data["attributes"]
        mol.attributes[Symbol(k)] = v
    end
    return mol
end

graphmol(json::String) = graphmol(JSON.parse(json))

"""
    setcache!(graph::, key; kwargs...)

Set calculated property caches.
"""
function Graph.setcache!(mol::GraphMol, key; kwargs...)
    mol.cache[key] = getfield(MolecularGraph, key)(mol; kwargs...)
end

# Aliases

SDFile = GraphMol{SDFileAtom,SDFileBond}
SMILES = GraphMol{SmilesAtom,SmilesBond}

function Base.getindex(graph::GraphMol, sym::Symbol)
    if haskey(graph.cache, sym)
        return graph.cache[sym]
    else
        return eval(Expr(:call, sym, graph))
    end
end

function Base.getindex(view::SubgraphView, sym::Symbol)
    if haskey(view.graph.cache, sym)
        return view.graph.cache[sym]
    else
        return eval(Expr(:call, sym, view.graph))
    end
end

ext: jl


text:
using CovarianceMatrices
using Random
using GLM
using DataFrames

Random.seed!(1)

u = zeros(6000*50)
for j in 1:2999
    u[j+1] = 0.97*u[j] + randn()
end

df = DataFrame(y = randn(300*50) .+ u[15001:30000])

for j in Symbol.("x".*string.(collect(1:5)))
    df[!, j] = randn(300*50)
end
df[!, :cluster] = repeat(1:50, inner = [300])

df2 = DataFrame(y = sqrt(2).*randn(250))
for j in Symbol.("x".*string.(collect(1:5)))
    df2[!, j] .= randn(250)
end
df2[!, :cluster] .= repeat(1:50, inner = [5])

frm = @formula(y ~ x1 + x2 + x3 + x4 + x5)
lm1 = glm(frm, df, Normal(), IdentityLink())
lm2 = glm(frm, df2, Normal(),IdentityLink())

using Profile
Profile.clear()
@profile for j in 1:200; vcov(CRHC0(df[!,:cluster]), lm1); end
Juno.profiler()
using LinearAlgebra
k0 = CRHC0(df[!,:cluster])
cache = CovarianceMatrices.install_cache(k0, lm1)

using BenchmarkTools
@btime CovarianceMatrices.__vcov(k0, lm1, cache, Matrix, Cholesky, CovarianceMatrices.dofadjustment(k0, cache))

k2 = CRHC2(df[!,:cluster])
cache = CovarianceMatrices.install_cache(k2, lm1)
@btime CovarianceMatrices.__vcov(k2, lm1, cache, Matrix, Cholesky, CovarianceMatrices.dofadjustment(k2, cache))


text:
module DFMs

include("HelperFunctions.jl")
include("EstimateDFM.jl")
include("Update.jl")

end


text:
@testset "ReducedSpaceEvaluators ($case)" for case in ["case9.m", "case30.m"]
    if has_cuda_gpu()
        ITERATORS = zip([CPU(), CUDADevice()], [Array, CuArray])
    else
        ITERATORS = zip([CPU()], [Array])
    end
    datafile = joinpath(INSTANCES_DIR, case)
    @testset "Constructor" for (device, M) in ITERATORS
        powerflow_solver = NewtonRaphson(; tol=1e-11)
        nlp = ExaPF.ReducedSpaceEvaluator(datafile; device=device, powerflow_solver=powerflow_solver)
        # Test that arguments are passed as expected in the constructor:
        @test nlp.powerflow_solver == powerflow_solver
    end

    pf = PowerSystem.PowerNetwork(datafile)
    @testset "Test API on $device" for (device, M) in ITERATORS
        polar = PolarForm(pf, device)
        x0 = ExaPF.initial(polar, State())
        u0 = ExaPF.initial(polar, Control())

        constraints = Function[
            ExaPF.voltage_magnitude_constraints,
            ExaPF.active_power_constraints,
            ExaPF.reactive_power_constraints,
            ExaPF.flow_constraints,
        ]
        nlp = ExaPF.ReducedSpaceEvaluator(polar; constraints=constraints)
        # Test printing
        println(devnull, nlp)

        # Test consistence
        n_state = get(polar, NumberOfState())
        @test length(nlp.λ) == n_state
        n = ExaPF.n_variables(nlp)
        m = ExaPF.n_constraints(nlp)
        @test n == length(u0)
        @test isless(nlp.u_min, nlp.u_max)
        @test isless(nlp.g_min, nlp.g_max)
        @test length(nlp.g_min) == m

        # Test API
        @test isa(get(nlp, ExaPF.Constraints()), Array{Function})
        @test get(nlp, State()) == x0
        buffer = get(nlp, ExaPF.PhysicalState())
        @test isa(buffer, ExaPF.AbstractBuffer)
        @test ExaPF.initial(nlp) == u0

        vm, va, pg, qg = buffer.vmag, buffer.vang, buffer.pinj, buffer.qinj

        u = u0
        # Update nlp to stay on manifold
        ExaPF.update!(nlp, u)
        # Compute objective
        c = ExaPF.objective(nlp, u)
        @test isa(c, Real)
        # Compute gradient of objective
        g = similar(u)
        fill!(g, 0)
        ExaPF.gradient!(nlp, g, u)
        function reduced_cost(u_)
            ExaPF.update!(nlp, u_)
            return ExaPF.objective(nlp, u_)
        end
        grad_fd = FiniteDiff.finite_difference_gradient(reduced_cost, u)
        @test isapprox(grad_fd, g, rtol=1e-6)

        ExaPF.update!(nlp, u)
        hv = similar(u) ; fill!(hv, 0)
        w = similar(u) ; fill!(w, 0)
        w[1] = 1.0
        ExaPF.hessprod!(nlp, hv, u, w)
        H = similar(u, n, n) ; fill!(H, 0)
        ExaPF.hessian!(nlp, H, u)
        # Is Hessian vector product relevant?
        @test H * w == hv
        # Is Hessian correct?
        hess_fd = FiniteDiff.finite_difference_hessian(reduced_cost, u)
        # Take attribute data as hess_fd is of type Symmetric
        @test H ≈ hess_fd.data rtol=1e-6

        # Constraint
        ## Evaluation of the constraints
        cons = similar(nlp.g_min)
        fill!(cons, 0)
        ExaPF.constraint!(nlp, cons, u)

        ## Evaluation of the transpose-Jacobian product
        v = similar(cons) ; fill!(v, 0)
        fill!(g, 0)
        ExaPF.jtprod!(nlp, g, u, v)
        @test iszero(g)
        fill!(v, 1) ; fill!(g, 0)
        ExaPF.jtprod!(nlp, g, u, v)

        ## Evaluation of the Jacobian (only on CPU)
        if isa(device, CPU)
            jac = M{Float64, 2}(undef, m, n)
            ExaPF.jacobian!(nlp, jac, u)
            # Test transpose Jacobian vector product
            @test isapprox(g, transpose(jac) * v)
            # Test Jacobian vector product
            S = typeof(u)
            jv = ExaPF.xzeros(S, m)
            v = ExaPF.xzeros(S, n)
            ExaPF.jprod!(nlp, jv, u, v)
            @test jac * v == jv
        end

        # Utils
        inf_pr1 = ExaPF.primal_infeasibility(nlp, u)
        inf_pr2 = ExaPF.primal_infeasibility!(nlp, cons, u)
        @test inf_pr1 == inf_pr2

        # test reset!
        ExaPF.reset!(nlp)
        @test get(nlp, State()) == x0
        @test iszero(nlp.λ)

        # setters
        nbus = get(nlp, PS.NumberOfBuses())
        loads = similar(u, nbus) ; fill!(loads, 1)
        ExaPF.setvalues!(nlp, PS.ActiveLoad(), loads)
        ExaPF.setvalues!(nlp, PS.ReactiveLoad(), loads)
    end
end


text:
# Apply c-posterior to autoregressive model of variable order.
# Using analytical expression for evidence, assuming known variance.
module Run
using PyPlot
draw() = (pause(0.001); get_current_fig_manager()[:window][:raise_]())

# Settings
ns = [10 .^ [2:4]]
alphas = [Inf,500] # robustification param (max sample size)
K = 20 # maximum value of k to use

# Data generation parameters
s0 = 1.0 # std dev of noise
sm = 0.5 # scale of misspecification
a0 = [0.25,0.25,-0.25,0.25] # autoregression coefficients
k0 = length(a0)
f_noise(t) = randn()*s0 + sm*sin(t) # Time-varying noise

# Model parameters
s = s0
sa = 1.0 # std dev of prior on a's

# Code for autoregressive model
include("AR.jl")

srand(0) # Reset RNG
x_all = generate(maximum(ns),a0,f_noise) # Generate data
for (i_a,alpha) in enumerate(alphas)
    for (i_n,n) in enumerate(ns)
        zeta = (1/n)/(1/n + 1/alpha) # robustification param
        x = x_all[1:n]

        # Compute marginal likelihood for each k
        log_m = zeros(K+1)
        for k = 0:K
            log_m[k+1] = log_marginal(x,k,s,sa,zeta)
        end

        # Plot
        figure(i_a*10+i_n,figsize=(8,2.5)); clf(); hold(true)
        subplots_adjust(top=.8,bottom=0.25,left=.13)
        plot(0:K,log_m,"ko-",markersize=8,linewidth=2)
        xlim(0,K)
        xticks(0:K)
        mx,mn = maximum(log_m),minimum(log_m)
        pad = (mx-mn)/20
        ylim(mn-pad,mx+pad)
        ticklabel_format(axis="y",style="sci",scilimits=(-4,4))
        if alpha<Inf; title("\$\\mathrm{Coarsened},\\,n=$n\$",fontsize=17)
        else; title("\$\\mathrm{Standard},\\,n=$n\$",fontsize=17)
        end
        xlabel("\$k \\,\\,\\mathrm{(\\#\\,of\\,coefficients)}\$",fontsize=16)
        if n==ns[1]
            ylabel("\$\\log\\,p(\\mathrm{data}|k)\$",fontsize=15)
        end
        draw()
        savefig("AR-alpha=$alpha-n=$n.png",dpi=150)
    end
end

figure(3,figsize=(8,3.2)); clf(); hold(true)
subplots_adjust(bottom=0.2)
x = generate(200,a0,f_noise)
plot(x)
title("\$\\mathrm{AR($k0)\\,data\\,with\\,time-varying\\,noise}\$",fontsize=17)
xlabel("\$t \\, \\mathrm{(time)}\$",fontsize=16)
ylabel("\$x_t\$",fontsize=18)
ylim(-4,4)
draw()
savefig("AR-data.png",dpi=150)

end # module

max_stars_repo_name: jwmi/CoarsenedPosterior

max_issues_repo_name: jwmi/CoarsenedPosterior

max_forks_repo_name: jwmi/CoarsenedPosterior


text:
function max_error(g_true, g_approx)
    n = length(g_true)
    @assert length(g_approx) == n
    err = 0.0
    for i in 1:n
        abs_err = g_true[i] - g_approx[i]
        err = max(err, abs(abs_err))
    end
    return err
end


text:
# See also fastmath_able.jl for where rules are defined simple base functions
# that also have FastMath versions.

@scalar_rule one(x) zero(x)
@scalar_rule zero(x) zero(x)
@scalar_rule transpose(x) One()

# `adjoint`

frule((_, Δz), ::typeof(adjoint), z::Number) = (z', Δz')

function rrule(::typeof(adjoint), z::Number)
    adjoint_pullback(ΔΩ) = (NO_FIELDS, ΔΩ')
    return (z', adjoint_pullback)
end

# `real`

@scalar_rule real(x::Real) One()

frule((_, Δz), ::typeof(real), z::Number) = (real(z), real(Δz))

function rrule(::typeof(real), z::Number)
    # add zero(z) to embed the real number in the same number type as z
    real_pullback(ΔΩ) = (NO_FIELDS, real(ΔΩ) + zero(z))
    return (real(z), real_pullback)
end

# `imag`

@scalar_rule imag(x::Real) Zero()

frule((_, Δz), ::typeof(imag), z::Complex) = (imag(z), imag(Δz))

function rrule(::typeof(imag), z::Complex)
    imag_pullback(ΔΩ) = (NO_FIELDS, real(ΔΩ) * im)
    return (imag(z), imag_pullback)
end

# `Complex`

frule((_, Δz), ::Type{T}, z::Number) where {T<:Complex} = (T(z), Complex(Δz))
function frule((_, Δx, Δy), ::Type{T}, x::Number, y::Number) where {T<:Complex}
    return (T(x, y), Complex(Δx, Δy))
end

function rrule(::Type{T}, z::Complex) where {T<:Complex}
    Complex_pullback(ΔΩ) = (NO_FIELDS, Complex(ΔΩ))
    return (T(z), Complex_pullback)
end
function rrule(::Type{T}, x::Real) where {T<:Complex}
    Complex_pullback(ΔΩ) = (NO_FIELDS, real(ΔΩ))
    return (T(x), Complex_pullback)
end
function rrule(::Type{T}, x::Number, y::Number) where {T<:Complex}
    Complex_pullback(ΔΩ) = (NO_FIELDS, real(ΔΩ), imag(ΔΩ))
    return (T(x, y), Complex_pullback)
end

# `hypot`

@scalar_rule hypot(x::Real) sign(x)

function frule((_, Δz), ::typeof(hypot), z::Complex)
    Ω = hypot(z)
    ∂Ω = _realconjtimes(z, Δz) / ifelse(iszero(Ω), one(Ω), Ω)
    return Ω, ∂Ω
end

function rrule(::typeof(hypot), z::Complex)
    Ω = hypot(z)
    function hypot_pullback(ΔΩ)
        return (NO_FIELDS, (real(ΔΩ) / ifelse(iszero(Ω), one(Ω), Ω)) * z)
    end
    return (Ω, hypot_pullback)
end

@scalar_rule fma(x, y, z) (y, x, One())
@scalar_rule muladd(x, y, z) (y, x, One())
@scalar_rule rem2pi(x, r::RoundingMode) (One(), DoesNotExist())
@scalar_rule(
    mod(x, y),
    @setup((u, nan) = promote(x / y, NaN16), isint = isinteger(x / y)),
    (ifelse(isint, nan, one(u)), ifelse(isint, nan, -floor(u))),
)

@scalar_rule deg2rad(x) π / oftype(x, 180)
@scalar_rule rad2deg(x) oftype(x, 180) / π

@scalar_rule(ldexp(x, y), (2^y, DoesNotExist()))

# Can't multiply though sqrt in acosh because of negative complex case for x
@scalar_rule acosh(x) inv(sqrt(x - 1) * sqrt(x + 1))
@scalar_rule acoth(x) inv(1 - x ^ 2)
@scalar_rule acsch(x) -(inv(x ^ 2 * sqrt(1 + x ^ -2)))
@scalar_rule acsch(x::Real) -(inv(abs(x) * sqrt(1 + x ^ 2)))
@scalar_rule asech(x) -(inv(x * sqrt(1 - x ^ 2)))
@scalar_rule asinh(x) inv(sqrt(x ^ 2 + 1))
@scalar_rule atanh(x) inv(1 - x ^ 2)

@scalar_rule acosd(x) (-(oftype(x, 180)) / π) / sqrt(1 - x ^ 2)
@scalar_rule acotd(x) (-(oftype(x, 180)) / π) / (1 + x ^ 2)
@scalar_rule acscd(x) ((-(oftype(x, 180)) / π) / x ^ 2) / sqrt(1 - x ^ -2)
@scalar_rule acscd(x::Real) ((-(oftype(x, 180)) / π) / abs(x)) / sqrt(x ^ 2 - 1)
@scalar_rule asecd(x) ((oftype(x, 180) / π) / x ^ 2) / sqrt(1 - x ^ -2)
@scalar_rule asecd(x::Real) ((oftype(x, 180) / π) / abs(x)) / sqrt(x ^ 2 - 1)
@scalar_rule asind(x) (oftype(x, 180) / π) / sqrt(1 - x ^ 2)
@scalar_rule atand(x) (oftype(x, 180) / π) / (1 + x ^ 2)

@scalar_rule cot(x) -((1 + Ω ^ 2))
@scalar_rule coth(x) -(csch(x) ^ 2)
@scalar_rule cotd(x) -(π / oftype(x, 180)) * (1 + Ω ^ 2)
@scalar_rule csc(x) -Ω * cot(x)
@scalar_rule cscd(x) -(π / oftype(x, 180)) * Ω * cotd(x)
@scalar_rule csch(x) -(coth(x)) * Ω
@scalar_rule sec(x) Ω * tan(x)
@scalar_rule secd(x) (π / oftype(x, 180)) * Ω * tand(x)
@scalar_rule sech(x) -(tanh(x)) * Ω

@scalar_rule acot(x) -(inv(1 + x ^ 2))
@scalar_rule acsc(x) -(inv(x ^ 2 * sqrt(1 - x ^ -2)))
@scalar_rule acsc(x::Real) -(inv(abs(x) * sqrt(x ^ 2 - 1)))
@scalar_rule asec(x) inv(x ^ 2 * sqrt(1 - x ^ -2))
@scalar_rule asec(x::Real) inv(abs(x) * sqrt(x ^ 2 - 1))

@scalar_rule cosd(x) -(π / oftype(x, 180)) * sind(x)
@scalar_rule cospi(x) -π * sinpi(x)
@scalar_rule sind(x) (π / oftype(x, 180)) * cosd(x)
@scalar_rule sinpi(x) π * cospi(x)
@scalar_rule tand(x) (π / oftype(x, 180)) * (1 + Ω ^ 2)

@scalar_rule(
    clamp(x, low, high),
    @setup(
        islow = x < low,
        ishigh = high < x,
    ),
    (!(islow | ishigh), islow, ishigh),
)
@scalar_rule x \ y (-((y / x) / x), inv(x))

function frule((_, ẏ), ::typeof(identity), x)
    return (x, ẏ)
end

function rrule(::typeof(identity), x)
    function identity_pullback(ȳ)
        return (NO_FIELDS, ȳ)
    end
    return (x, identity_pullback)
end

#####
##### `evalpoly`
#####

if VERSION ≥ v"1.4"
    function frule((_, Δx, Δp), ::typeof(evalpoly), x, p)
        N = length(p)
        @inbounds y = p[N]
        Δy = Δp[N]
        @inbounds for i in (N - 1):-1:1
            Δy = muladd(Δx, y, muladd(x, Δy, Δp[i]))
            y = muladd(x, y, p[i])
        end
        return y, Δy
    end

    function rrule(::typeof(evalpoly), x, p)
        y, ys = _evalpoly_intermediates(x, p)
        function evalpoly_pullback(Δy)
            ∂x, ∂p = _evalpoly_back(x, p, ys, Δy)
            return NO_FIELDS, ∂x, ∂p
        end
        return y, evalpoly_pullback
    end

    # evalpoly but storing intermediates
    function _evalpoly_intermediates(x, p::Tuple)
        return if @generated
            N = length(p.parameters)
            exs = []
            vars = []
            ex = :(p[$N])
            for i in 1:(N - 1)
                yi = Symbol("y", i)
                push!(vars, yi)
                push!(exs, :($yi = $ex))
                ex = :(muladd(x, $yi, p[$(N - i)]))
            end
            push!(exs, :(y = $ex))
            Expr(:block, exs..., :(y, ($(vars...),)))
        else
            _evalpoly_intermediates_fallback(x, p)
        end
    end
    function _evalpoly_intermediates_fallback(x, p::Tuple)
        N = length(p)
        y = p[N]
        ys = (y, ntuple(N - 2) do i
            return y = muladd(x, y, p[N - i])
        end...)
        y = muladd(x, y, p[1])
        return y, ys
    end
    function _evalpoly_intermediates(x, p)
        N = length(p)
        @inbounds yn = one(x) * p[N]
        ys = similar(p, typeof(yn), N - 1)
        @inbounds ys[1] = yn
        @inbounds for i in 2:(N - 1)
            ys[i] = muladd(x, ys[i - 1], p[N - i + 1])
        end
        @inbounds y = muladd(x, ys[N - 1], p[1])
        return y, ys
    end

    # TODO: Handle following cases
    #     1) x is a UniformScaling, pᵢ is a matrix
    #     2) x is a matrix, pᵢ is a UniformScaling
    @inline _evalpoly_backx(x, yi, ∂yi) = ∂yi * yi'
    @inline _evalpoly_backx(x, yi, ∂x, ∂yi) = muladd(∂yi, yi', ∂x)
    @inline _evalpoly_backx(x::Number, yi, ∂yi) = conj(dot(∂yi, yi))
    @inline _evalpoly_backx(x::Number, yi, ∂x, ∂yi) = _evalpoly_backx(x, yi, ∂yi) + ∂x

    @inline _evalpoly_backp(pi, ∂yi) = ∂yi

    function _evalpoly_back(x, p::Tuple, ys, Δy)
        return if @generated
            exs = []
            vars = []
            N = length(p.parameters)
            for i in 2:(N - 1)
                ∂pi = Symbol("∂p", i)
                push!(vars, ∂pi)
                push!(exs, :(∂x = _evalpoly_backx(x, ys[$(N - i)], ∂x, ∂yi)))
                push!(exs, :($∂pi = _evalpoly_backp(p[$i], ∂yi)))
                push!(exs, :(∂yi = x′ * ∂yi))
            end
            push!(vars, :(_evalpoly_backp(p[$N], ∂yi))) # ∂pN
            Expr(
                :block,
                :(x′ = x'),
                :(∂yi = Δy),
                :(∂p1 = _evalpoly_backp(p[1], ∂yi)),
                :(∂x = _evalpoly_backx(x, ys[$(N - 1)], ∂yi)),
                :(∂yi = x′ * ∂yi),
                exs...,
                :(∂p = (∂p1, $(vars...))),
                :(∂x, Composite{typeof(p),typeof(∂p)}(∂p)),
            )
        else
            _evalpoly_back_fallback(x, p, ys, Δy)
        end
    end
    function _evalpoly_back_fallback(x, p::Tuple, ys, Δy)
        x′ = x'
        ∂yi = Δy
        N = length(p)
        ∂p1 = _evalpoly_backp(p[1], ∂yi)
        ∂x = _evalpoly_backx(x, ys[N - 1], ∂yi)
        ∂yi = x′ * ∂yi
        ∂p = (
            ∂p1,
            ntuple(N - 2) do i
                ∂x = _evalpoly_backx(x, ys[N-i-1], ∂x, ∂yi)
                ∂pi = _evalpoly_backp(p[i+1], ∂yi)
                ∂yi = x′ * ∂yi
                return ∂pi
            end...,
            _evalpoly_backp(p[N], ∂yi), # ∂pN
        )
        return ∂x, Composite{typeof(p),typeof(∂p)}(∂p)
    end
    function _evalpoly_back(x, p, ys, Δy)
        x′ = x'
        ∂yi = one(x′) * Δy
        N = length(p)
        @inbounds ∂p1 = _evalpoly_backp(p[1], ∂yi)
        ∂p = similar(p, typeof(∂p1))
        @inbounds begin
            ∂x = _evalpoly_backx(x, ys[N - 1], ∂yi)
            ∂yi = x′ * ∂yi
            ∂p[1] = ∂p1
            for i in 2:(N - 1)
                ∂x = _evalpoly_backx(x, ys[N - i], ∂x, ∂yi)
                ∂p[i] = _evalpoly_backp(p[i], ∂yi)
                ∂yi = x′ * ∂yi
            end
            ∂p[N] = _evalpoly_backp(p[N], ∂yi)
        end
        return ∂x, ∂p
    end
end


text:
module TestScitypes

using MLJ
using Test

S = scitype(ConstantClassifier())
@test S().prediction_type == :probabilistic
U = scitype(FeatureSelector())
@test U().input_scitype == MLJ.Table

@test scitype(OneHotEncoder()) ==
    MLJ.UnsupervisedScitype{Table,Table}

@test scitype(ConstantRegressor()) ==
    MLJ.SupervisedScitype{Table,
                          AbstractVector{Continuous},
                          :probabilistic}

end
true


text:
"""
    LassoResult

Result of lasso estimation for a grid of values of regularization
parameter. Contains:

| Atributes | Type               | Brief Description                                            |
|:--------- |:------------------ |:------------------------------------------------------------ |
| mdl       | `NarmaxRegressors` | Regressors                                                   |
| λ         | `Matrix{Float64}`  | Sequence of regularization parameters λ                      |
| β         | `Matrix{Float64}`  | Parameter vectors arranged columnwise, for the sequence of λ |
| df        | `Matrix{Int}`      | Number of nonzero terms for the sequence of λ                |
| time      | `Float64`          | Execution time                                               |
| total_iter| `Int`              | Total number of internal iterations                          |
"""
struct LassoResult
    mdl::NarmaxRegressors
    β::Matrix{Float64}
    λ::Vector{Float64}
    df::Vector{Int}
    time::Float64
    total_iter::Int
end

# Update noise terms
function update_noise_terms!(
        X::Vector{Vector{Float64}}, i::Int, N::Int, r::Vector{Float64},
        noise_terms_start::Int, noise_lags::Vector{Int})
    if i >= noise_terms_start
        lag!(X[i], r, noise_lags[i-noise_terms_start+1], N)
    end
end

# Update basis matrix
function update_basis_matrix!{T<:Basis}(
    X::Vector{Vector{Float64}}, i::Int, N::Int, r::Vector{Float64},
    noise_components_start::Int,
    noise_terms_start::Int,
    noise_lags::Vector{Int},
    basis::Vector{T},
    Y::Vector{Vector{Float64}},
    p::Int)
    if i >= noise_components_start
        evaluate!(X[i], Y, basis[i], N, p, update_noise_terms!,
                  (r, noise_terms_start, noise_lags))
    end
end

## Narmax Lasso function
"""
    narmax_lasso(y, u, mdl::NarmaxRegressors[; use_glmnet]) -> r::LassoResult

Given input and output signals `u` and `y`, and a structure `mdl` specifying
the regressors through a structure [`NarmaxRegressors`](@ref). Find the solution
of the lasso problem:
```math
\\min_\\beta ||\\mathbf{e}||^2 + \\lambda * \\sum_{i} |\\beta_i|
```
for a grid of values of the regularization parameter ``\\lambda``. Return
a [`LassoResult`](@ref) structure containing the obtained values. The option
`use_glmnet=true` changes the internal solver to GLMNet (only available
when no noise term is present).
"""
function narmax_lasso(y::Vector{Float64}, u::Vector{Float64},
                      mdl::NarmaxRegressors; use_glmnet=false,
                      kwargs...)
    # Lag matrix
    v = copy(y)
    Y = lag_matrix(mdl, y, u, v)
    Y_views, N, n = matrix2subviews(Y)
    # Generate regressor matrix
    X = regressor_matrix(mdl, y, u, v)
    # Pathwise coordinate descent
    tic()
    if use_glmnet
        if isempty(mdl.vterms) && glmnet_loaded
            path = GLMNet.glmnet(X, y, intercept=false,
                          standardize=false, tol=1e-8; kwargs...)
            β = Matrix{Float64}(path.betas)
            λ = path.lambda
            df_matrix = [sum(β[:, i].==0) for i = 1:length(λ)]
            total_iter = path.npasses
        elseif ~glmnet_loaded
            throw(ArgumentError(string("GlMNet not installed. Run:\n",
                                "\tjulia> Pkg.add(\"GLMNet\")\n",
                                "\tjulia> Pkg.build(\"GLMNet\")")))
        else
            throw(ArgumentError(string("GLMNet can only be used when there ",
                                 "is no noise terms present")))
        end
    else
        β, λ, df_matrix, total_iter = pathwise_coordinate_optimization(
                X, y, update_basis_matrix!,
                (mdl.noise_components_starts, mdl.noise_terms_starts, mdl.vterms,
                 mdl.basis, Y_views, n); kwargs...)
    end
    time = toq()
    return LassoResult(mdl, β, λ, df_matrix,
                       time, total_iter)
end


text:
"""
    AbstractAudioGenerator

Interface for streaming audio generators.  Concrete types must implement
`reset!`, `has_samples`, and `next_samples!`.
"""
abstract type AbstractAudioGenerator end

"""
    NullAudioGenerator

An audio generator containing no audio samples.
"""
struct NullAudioGenerator <: AbstractAudioGenerator end
reset!(gen::NullAudioGenerator) = nothing
has_samples(gen::NullAudioGenerator, sample_rate::Float64) = false
function next_samples!(gen::NullAudioGenerator, max_n::Int,
                       sample_rate::Float64)
    @assert(false, "NullAudioGenerator has no samples")
end


text:
function stderr_SparKLIE(Ψx, Ψy, θ, ω)
    nx = size(Ψx, 2)
    ny = size(Ψy, 2)
    supp = findall(!iszero, ω)
    S = (cov(transpose(Ψx[supp, :]), corrected=false) ./ nx) .+ (cov(rhat(θ, Ψy) .* transpose(Ψy[supp, :]), corrected=false) ./ ny)
    σ2 = 0.
    for k in 1:length(supp)
        for l in 1:length(supp)
            σ2 += S[k,l] * ω[supp][k] * ω[supp][l]
        end
    end
    sqrt(σ2)
end


text:
module FrameConnection

include("structdefinitions.jl")
include("isolatesmd.jl")
include("precluster.jl")
include("organizeclusters.jl")
include("computeclusterinfo.jl")
include("estimateparams.jl")
include("estimatedensities.jl")
include("create_costmatrix.jl")
include("solveLAP.jl")
include("linkclusters.jl")
include("smoothen_connectID.jl")
include("combinelocalizations.jl")
include("connectlocalizations.jl")
include("frameconnect.jl")

end


text:
#raw kRPC constructs
function SendRawProto(outp::IO, obj::Any)
    iob = PipeBuffer()
    writeproto(iob, obj)
    bytes = readavailable(iob)
    varib = PipeBuffer()
    ProtoBuf.write_varint(varib, length(bytes))
    write(outp, readavailable(varib), bytes)
    flush(outp)
end

const MSB = 0x80
const MASK7 = 0x7f
function RecvRawProto(inp::IO)
    flush(inp)
    iob = PipeBuffer()
    res = zero(Int64)
    n = 0
    byte = UInt8(MSB)
    while (byte & MSB) != 0
        flush(inp)
        byte = @sync read(inp, UInt8)
        res |= (convert(Int64, byte & MASK7) << (7*n))
        n += 1
    end
    write(iob, read(inp, res))
    return iob
end

function SendBiMessage(conn::kRPCConnection, req::kRPC.krpc.schema.Request)
    iob = PipeBuffer()
    SendRawProto(conn.conn, req)
    res = readproto(RecvRawProto(conn.conn), krpc.schema.Response())

    if has_field(res, :error)
        throw(make_error(res.error))
    end
    return res
end

function GetInfo(conn::kRPCConnection)
    rv = SendBiMessage(conn,
        krpc.schema.Request(calls=[krpc.schema.ProcedureCall(service="KRPC", procedure="GetStatus")]))
    if ProtoBuf.has_field(rv, :error)
        throw(make_error(rv.error))
    end
    if ProtoBuf.has_field(rv.results[1], :error)
        throw(make_error(rv.results[1].error))
    end
    iob = PipeBuffer()
    write(iob, rv.results[1].value)
    return readproto(iob, krpc.schema.Status())
end

function GetServices(conn::kRPCConnection)
    rv = SendBiMessage(conn,
        krpc.schema.Request(calls=[krpc.schema.ProcedureCall(service="KRPC", procedure="GetServices")]))
    if ProtoBuf.has_field(rv, :error)
        throw(make_error(rv.error))
    end
    if ProtoBuf.has_field(rv.results[1], :error)
        throw(make_error(rv.results[1].error))
    end
    iob = PipeBuffer()
    write(iob, rv.results[1].value)
    return readproto(iob, krpc.schema.Services())
end


text:
# define and test barycentric velocity interpolation, to be used
# for computing seepage force in particle simulations

using LinearAlgebra: dot

"""
computes velocity of single particle without using cell list

efficient in that it allocates no new memory, but inefficient in that
it loops over all elements (no cell list given)
"""
function compute_seepage_force_nocl!(sfX::Vector{T},sfY::Vector{T},mesh,u::Vector{T},v::Vector{T},pList,data) where T<:Real
    fill!(sfX,zero(T))
    fill!(sfY,zero(T))

    for ti=1:data.n_particles
        elFound = false
        elInd = 0
        data.extremePoint.y = pList[ti].pos.y

        # find which element the particle belongs to
        for el=1:length(mesh.cm)
            # construct polygon
            for ti=1:4
                data.polygon[ti].x = mesh.xy[mesh.cm[el].NodeList[ti]].x
                data.polygon[ti].y = mesh.xy[mesh.cm[el].NodeList[ti]].y
            end

            if StokesParticles.is_inside(data.polygon,4,pList[ti].pos;extreme=data.extremePoint)
                elInd = el
                elFound = true
                break
            end
        end

        # find weights for barycentric
        computeBaryWeights!(data.w,data.polygon,pList[ti].pos,data.a,data.b,data.c,data.d)

        if elFound
            # interpolate values
            for tj=1:4
                data.uEl[tj] = u[mesh.cm[elInd].NodeList[tj]]
                data.vEl[tj] = v[mesh.cm[elInd].NodeList[tj]]
            end

            sfX[ti] = dot(data.uEl, data.w)
            sfY[ti] = dot(data.vEl, data.w)
        else
            sfX[ti] = 0.0
            sfY[ti] = 0.0
        end
    end

    nothing
end

function compute_single_seepage_force_nocl(mesh,u::Vector{T},v::Vector{T},particle,data) where T<:Real
    sfX = 0.0
    sfY = 0.0
    elFound = false
    elInd = 0
    data.extremePoint.y = particle.pos.y

    # find which element the particle belongs to
    for el=1:length(mesh.cm)
        # construct polygon
        for ti=1:4
            data.polygon[ti].x = mesh.xy[mesh.cm[el].NodeList[ti]].x
            data.polygon[ti].y = mesh.xy[mesh.cm[el].NodeList[ti]].y
        end

        if StokesParticles.is_inside(data.polygon,4,particle.pos;extreme=data.extremePoint)
            elInd = el
            elFound = true
            break
        end
    end

    # find weights for barycentric
    computeBaryWeights!(data.w,data.polygon,particle.pos,data.a,data.b,data.c,data.d)

    if elFound
        # interpolate values
        for tj=1:4
            data.uEl[tj] = u[mesh.cm[elInd].NodeList[tj]]
            data.vEl[tj] = v[mesh.cm[elInd].NodeList[tj]]
        end

        sfX = dot(data.uEl, data.w)
        sfY = dot(data.vEl, data.w)
    end

    return sfX, sfY
end

"""
efficient velocity interpolation.

Given nodelist of particles (via point positions), interpolates velocity from FEM using
cell list.

for max efficiency, velocities are interpolated according to cell ordering, NOT particle
index ordering

passes in a vector for velocity interpolations (expected u and v are same size!!)
"""
function compute_seepage_force!(sfX::Vector{T},sfY::Vector{T},mesh,u::Vector{T},v::Vector{T},pList,data,cl,femap) where T<:Real
    fill!(sfX,zero(T))
    fill!(sfY,zero(T))

    # loop over cells
    for cellInd = 1:length(cl.cells)
        for nodeInd in cl.cells[cellInd].particleIDList
            foundNode = false

            data.extremePoint.y = pList[nodeInd].pos.y

            # check whether cell belongs to element INSIDE cell
            for elInd in femap[cellInd]
                # construct polygon -- FEM element
                for ti=1:4
                    data.polygon[ti].x = mesh.xy[mesh.cm[elInd].NodeList[ti]].x
                    data.polygon[ti].y = mesh.xy[mesh.cm[elInd].NodeList[ti]].y
                end

                if StokesParticles.is_inside(data.polygon,4,pList[nodeInd].pos;extreme=data.extremePoint)
                    foundNode = true
                    # find weights for barycentric
                    computeBaryWeights!(data.w,data.polygon,pList[nodeInd].pos,data.a,data.b,data.c,data.d)

                    # interpolate values
                    for ti=1:4
                        data.uEl[ti] = u[mesh.cm[elInd].NodeList[ti]]
                        data.vEl[ti] = v[mesh.cm[elInd].NodeList[ti]]
                    end
                    data.sfX[nodeInd] = dot(data.uEl, data.w)
                    data.sfY[nodeInd] = dot(data.vEl, data.w)
                    break
                end
            end

            if !(foundNode)
                sfX[nodeInd],sfY[nodeInd] = compute_single_seepage_force_nocl(mesh,u,v,pList[nodeInd],data)
            end
        end
    end

    nothing
end

max_stars_repo_name: pseastham/eFEMpart

max_issues_repo_name: pseastham/eFEMpart


text:
# utilities
# commented out are imported functions that haven't been thoroughly tested

# Migrate these to MaskRCNN when done testing
using StatsBase
using StatsBase: shuffle
using Base.Iterators
using Statistics

const _mcnurl = "http://www.vlfeat.org/matconvnet/models"
const _mcndir = joinpath(@__DIR__, "imagenet")

function matconvnet(name, pass = false)
    global _mcncache
    if !@isdefined(_mcncache); _mcncache=Dict(); end
    if !haskey(_mcncache, name)
        matfile = "$name.mat"
        @info("Loading $matfile...")
        path = joinpath(_mcndir,matfile)
        if !isfile(path)
            if !pass
                println("Should I download $matfile? (y if Yes)")
                readline()[1] == 'y' || error(:ok)
                isdir(_mcndir) || mkpath(_mcndir)
                download("$_mcnurl/$matfile",path)
            else
                isdir(_mcndir) || mkpath(_mcndir)
                download("$_mcnurl/$matfile",path)
            end
        end
        _mcncache[name] = matread(path)
    end
    return _mcncache[name]
end

function imgdata(img, averageImage)
    global _imgcache
    if !@isdefined(_imgcache); _imgcache = Dict(); end
    if !haskey(_imgcache,img)
        if occursin("://",img)
            @info("Downloading $img")
            a0 = load(download(img))
        else
            a0 = load(img)
        end
        new_size = ntuple(i->div(size(a0,i)*224,minimum(size(a0))),2)
        a1 = Images.imresize(a0, new_size)
        i1 = div(size(a1,1)-224,2)
        j1 = div(size(a1,2)-224,2)
        b1 = a1[i1+1:i1+224,j1+1:j1+224]
        # ad-hoc solution for Mac-OS image
        macfix = convert(Array{FixedPointNumbers.Normed{UInt8,8},3}, channelview(b1))
        c1 = permutedims(macfix, (3,2,1))
        d1 = convert(Array{Float32}, c1)
        e1 = reshape(d1[:,:,1:3], (224,224,3,1))
        f1 = (255 * e1 .- averageImage)
        g1 = permutedims(f1, [2,1,3,4])
        _imgcache[img] = g1
    end
    return _imgcache[img]
end

# Resize image
function resize_image(img; min_dim = false, max_dim = false, padding = false)
    h, w = size(img)[1:2]
    window = (0.,0.,h,w)
    scale = 1.f0

    if min_dim > 0
	scale = max(1, min_dim / min(h, w))
    end

    if max_dim > 0
	image_max = max(h, w)
	if round(image_max * scale) > max_dim
	    scale = max_dim / image_max
	end
    end

    if scale != 1.
	img = imresize(img, (round(Int, h*scale), round(Int, w*scale), 3))
    end

    resized_img = zeros(max_dim, max_dim, 3)
    h, w = size(img)[1:2]
    top_pad = max(1, round(Int, (max_dim - h) / 2))
    bottom_pad = max_dim - h - top_pad
    left_pad = max(1, round(Int, (max_dim - w) / 2))
    right_pad = max_dim - w - left_pad
    resized_img[top_pad:(h+top_pad-1), left_pad:(w+left_pad-1), 1:3] .= img
    padding = [(top_pad, bottom_pad), (left_pad, right_pad), (0, 0)]
    window = (top_pad, left_pad, h + top_pad, w + left_pad)
    return resized_img, window, scale, padding
end

# Computes the IoU of a box
function compute_iou(box, boxes, box_area, boxes_area)
    #=
    Calculates the IoU of the given box with the array of the given boxes

    box: 1D vector [y1 x1 y2 x2]
    boxes: [boxes_count, (y1, x1, y2, x2)]
    box_area: float, the area of the box
    boxes_area:  array of the length of boxes_count

    Note: Areas are passed to avoid duplicate calculations
    =#
    y1 = max.(box[1], boxes[:,1])
    y2 = min.(box[3], boxes[:,3])
    x1 = max.(box[2], boxes[:,2])
    x2 = min.(box[4], boxes[:,4])
    intersection = max.(x2 .- x1, 0) .* max.(y2 .- y1, 0)
    union = box_area .+ boxes_area .- intersection
    iou = intersection ./ union
    return iou
end

# Computes overlaps of two sets of bboxes
function compute_overlaps_bbox(boxes1, boxes2)
    #=
    Computes the IoU overlaps between two sets of boxes

    boxes1, boxes2: [N, (y1, x1, y2, x2)]
    =#
    areas1 = (boxes1[:,3] .- boxes1[:,1]) .* (boxes1[:, 4] .- boxes1[:, 2])
    areas2 = (boxes2[:,3] .- boxes2[:,2]) .* (boxes2[:, 4] .- boxes2[:, 2])

    overlaps = zeros(size(boxes1,1), size(boxes2,1))
    for i = 1:size(overlaps, 2)
	box2 = boxes2[i,:]
	overlaps[:,i] = compute_iou(box2, boxes1, areas2[i], areas1)
    end
    return overlaps
end

# Computes the refinement needed to transform box to gt_box
function box_refinement_graph(box, gt_box)
    # box = Float32.(box)
    # gt_box = Float32.(gt_box)
    height = box[:, 3] .- box[:, 1]
    width = box[:, 4] .- box[:, 2]
    # width = replace(x -> x < 1f-3 ? 1.f0 : x, width)
    # height = replace(x -> x < 1f-3 ? 1.f0 : x, height)
    width = clamp.(width, 1f-3, 1.f0) # remove
    height = clamp.(height, 1f-3, 1.f0) # remove
    center_y = box[:, 1] .+ 0.5f0 .* height
    center_x = box[:, 2] .+ 0.5f0 .* width
    # return center_x
    gt_height = gt_box[:, 3] .- gt_box[:, 1]
    gt_width = gt_box[:, 4] .- gt_box[:, 2]
    gt_center_y = gt_box[:, 1] .+ 0.5f0 .* gt_height
    gt_center_x = gt_box[:, 2] .+ 0.5f0 .* gt_width

    dy = (gt_center_y .- center_y) ./ height
    dx = (gt_center_x .- center_x) ./ width
    # return dy
    dh = log.(abs.(gt_height ./ height))
    dw = log.(abs.(gt_width ./ width))

    return hcat(dy, dx, dh, dw)
end

#=
function apply_box_deltas(boxes, deltas)
    @show boxes[1,:]
    @show deltas[1,:]
    @show findall(isnan, boxes)
    @show findall(isnan, deltas)
    @show typeof(deltas)
    # bs = copy(boxes)

    heights = boxes[:,3] .- boxes[:,1]
    widths = boxes[:,4] .- boxes[:,2]
    centers_y = boxes[:,1] .+ (0.5f0 .* heights)
    centers_x = boxes[:,2] .+ (0.5f0 .* widths)

    centers_y = centers_y .+ (deltas[:,1] .* heights)
    centers_x = centers_x .+ (deltas[:,2] .* widths)
    # hs = deepcopy(heights)
    # @show findall(x -> isinf(x) || isnan(x), deltas[:,3])
    # @show findall(x -> isinf(x) || isnan(x), deltas[:,4])
    # @show findall(x -> isinf(x) || isnan(x), heights)
    # @show findall(x -> isinf(x) || isnan(x), widths)
    # @show findall(x -> isinf(x) || isnan(x), centers_x)
    # @show findall(x -> isinf(x) || isnan(x), centers_y)

    # @show maximum(deltas[:,3])
    # @show maximum(deltas[:,4])
    # @show maximum(heights)
    # @show maximum(widths)
    # @show maximum(centers_x)
    # @show maximum(centers_y)
    # @show std(deltas[:,3])
    # @show std(deltas[:,4])

    heights = heights .* exp.(deltas[:,3])
    widths = widths .* exp.(deltas[:,4])
    # global gh = heights
    # global gw = widths
    # error()

    infinds = findall(x -> isinf(x) || isnan(x), heights)
    @show deltas[:,3][infinds]
    @warn "something"
    # @show hs[findall(x -> isinf(x) || isnan(x), heights)]
    @show exp.(deltas[:,3])[findall(x -> isinf(x) || isnan(x), heights)]
    # @show findall(x -> isinf(x) || isnan(x), widths)

    # t = findall(isinf, heights)
    # v = findall(isinf, deltas[:,3])
    # @show v
    # # push!(t, 1)
    # @show exp.(deltas[:,3])[t]
    # @show hs[t]
    # @show hs[t] .* exp.(deltas[t,3])
    # p = findall(isinf, heights)
    # @show all(v.==p)
    # heights = replace(x -> isinf(x) ? 0.0f0 : x, heights)
    # widths = replace(x -> isinf(x) ? 0.0f0 : x, widths)
    y1s = centers_y .- (0.50f0 .* heights)
    x1s = centers_x .- (0.5f0 .* widths)
    y2s = y1s .+ heights
    x2s = x1s .+ widths
    Flux.stack([y1s, x1s, y2s, x2s], 2)
    # m = hcat(y1s, x1s, y2s, x2s)
    # ff = findall(isinf, heights)
    # if length(ff) > 0
    # 	@show size(y1s)
    # 	@show size(y2s)
    # 	t = findall(isinf, heights)
    # 	@show y1s[t]
    # 	@show heights[t]
    # 	@warn "nans in s"
    # 	@show "!!!!!!!!!!!!!!!"
    # end
end
=#
#=
function crop(boxes, window)
    y1 = clamp.(boxes[:,1], window[1], window[3])
    x1 = clamp.(boxes[:,2], window[2], window[4])
    y2 = clamp.(boxes[:,3], window[1], window[3])
    x2 = clamp.(boxes[:,4], window[2], window[4])
    hcat(y1,x1,y2,x2)
end
=#
function generate_anchors
end

function generate_pyramid_anchors
end

function add_mask
end

function add_bbox
end

function add_class
end

function resize_mask
end

function unmold_mask
end

# mask is the shape of an image wtf
# mask = rand(28,28, 10) => 10 masks of size 28 x 28 (to be applied to every channel)
function extract_bboxes(mask)
    nth = last(size(mask))
    boxes = zeros(Integer, nth, 4)
    for i =1:nth
	m = mask[:,:,i]
	cluster = findall(!iszero, m)
	if length(cluster) > 0
	    Is = map(x -> [x.I[1], x.I[2]], cluster) |> x -> hcat(x...)'
	    x1, x2 = extrema(Is[:,1])
	    y1, y2 = extrema(Is[:,2])
	else
	    x1 ,x2, y1, y2 = 0, 0, 0, 0
	end

	boxes[i,:] = [y1, x1, y2, x2]
    end
    return boxes
end

function upsample(x, scale_factor::Tuple)
    return repeat(x, inner = scale_factor)
end

function clip_boxes(boxes, window)
    y1 = max.(min.(boxes[:,1], window[:,3]), window[:,1])
    x1 = max.(min.(boxes[:,2], window[:,4]), window[:,2])
    y2 = max.(min.(boxes[:,3], window[:,3]), window[:,1])
    x2 = max.(min.(boxes[:,4], window[:,4]), window[:,2])
    return hcat(y1,x1,y2,x2)
end

#=
generate all the boxes that need to be refined and chosen and set
scales = rand(n)
ratios = rand(m)
shape = (28,28)
feature_stride = 0.2 # stride compared to the image to align the features with the image
anchor_stride = 2 # generate an anchor every `2` steps

pytorch expected [... ...
... ...]
the input shape remains this one ^

but we operate on it transposed [. .
. .

. .
. .]

For reference:
	xx, yy = np.meshgrid(l1, l2)
	xx == collect(flatten(l1'))' .* ones(length(l1), length(l2)))
	yy == collect(flatten(l2)) .* ones(length(l1), length(l2)))

	could also use `repeat`, but would allocate off the butt

	use `broadcasting` here

	xx == repeat(l1', length(l2))
=#
#=
function generate_anchors(scales, ratios, shape, feature_stride, anchor_stride)
   ls, lr = length(scales), length(ratios)
   sr = size(ratios)

   scales = collect(flatten(scales'))
   ratios = sqrt.(flatten(ratios'))

   heights = collect(flatten(scales / ratios))
   widths = collect(flatten(scales * ratios'))
   lw, lh = length(widths), length(heights)

   shifts_y = collect((0:anchor_stride:shape[1]) .* feature_stride)[1:end-1]
   shifts_x = collect((0:anchor_stride:shape[2]) .* feature_stride)[1:end-1]
   ly = length(shifts_y)
   lx = length(shifts_x)

   shifts_x = repeat(collect(shifts_x), ly)    #|> x -> reshape(x, ly,ly)
   shifts_y = repeat(collect(shifts_y), lx)  #|> x -> reshape(x, lx,lx)

   box_widths = repeat(widths', length(shifts_x))
   box_centers_x = repeat(shifts_x, 1, lw) # |> x -> reshape(x, length(widths), length(shifts_x))
   # @show lw, lh
   box_heights = repeat(heights', length(shifts_y))
   box_centers_y = repeat(shifts_y[1:ly], inner = (ly, lh)) # |> x -> reshape(x, length(widths), length(shifts_x))

   box_centers = Flux.stack([box_centers_y, box_centers_x], 3)
   global gbc = box_centers
   box_centers = permutedims(box_centers, (2,1,3)) |> x -> reshape(x, :, 2)

   box_sizes = Flux.stack([box_heights, box_widths], 3)
   global gbs = box_sizes
   box_sizes = permutedims(box_sizes, (2,1,3)) |> x -> reshape(x, :, 2)

   cat(box_centers .- (0.5f0 .* box_sizes), box_centers .+ (0.5f0 .* box_sizes), dims = 2)
   # return box_centers
end
=#
#=
function generate_anchors2(scales, ratios, shape, feature_stride, anchor_stride)
    ls = length(scales)
    lr = length(ratios)
    s = repeat(scales, lr)
    r = repeat(ratios, inner = ls)
    h = s ./ sqrt.(r)
    w = s .* sqrt.(r)

    shifts_y = collect((0:anchor_stride:shape[1]) * feature_stride)[1:end-1]
    shifts_x = collect((0:anchor_stride:shape[2]) * feature_stride)[1:end-1]

    lx = length(shifts_x)
    ly = length(shifts_y)
    shifts_x = repeat(shifts_x', ly)
    shifts_y = repeat(shifts_y, outer = (1,lx))

    lsx = length(shifts_x)
    lsy = length(shifts_y)
    lw = length(w)
    lh = length(h)
    box_widths = repeat(w', lsx)
    box_centers_x = repeat(shifts_x[1,:], outer = (size(shifts_x, 1),lw))

    box_heights = repeat(h', lsy)
    box_centers_y = repeat(shifts_y[:,1], inner = (size(shifts_y, 1),lh))

    box_centers = Flux.stack([box_centers_y, box_centers_x], 3)

    Flux.stack([box_centers_y, box_centers_x], 3)
end
=#

#=
function generate_pyramid_anchors(scales, ratios, feature_shapes,
				  feature_strides, anchor_stride)
    anchors = []
    # scales = transpose(scales)
    for i = 1:size(scales, 1)
	push!(anchors, generate_anchors(scales[i,:], ratios, feature_shapes[i,:],
					feature_strides[i], anchor_stride))
    end

    gpu(cat(anchors..., dims = 1))
end
=#
#=
inputs => (boxes, feature_maps...)
boxes => rand(10,4) (TODO: => rand(4,10,10)
feature_maps = similar(batch)
batch => rand(299, 299, 3, 10)
pool_size => [7,7]
=#
#=
function bbox_overlaps2(boxes1, boxes2)
    boxes1_repeat = (size(boxes2, 1), 1, 1)
    boxes2_repeat = (size(boxes1, 1), 1, 1)
    boxes1 = repeat(boxes1, outer = boxes1_repeat)
    boxes2 = repeat(boxes2, outer = boxes2_repeat)
    b1_y1, b2_y1 = boxes1[:, 1], boxes2[:, 1]
    b1_x1, b2_x1 = boxes1[:, 2], boxes2[:, 2]
    b1_y2, b2_y2 = boxes1[:, 3], boxes2[:, 3]
    b1_x2, b2_x2 = boxes1[:, 4], boxes2[:, 4]
    y1 = max.(b1_y1, b2_y1)
    x1 = max.(b1_x1, b2_x1)
    y2 = min.(b1_y2, b2_y2)
    x2 = min.(b1_x2, b2_x2)
    zs = zeros(size(y1, 1)) |> gpu
    intersection = max.(x2 .- x1 .+ 1, zs) .* max.(y2 .- y1 .+ 1, zs)
    b1_area = (b1_y2 .- b1_y1 .+ 1) .* (b1_x2 .- b1_x1 .+ 1)
    b2_area = (b2_y2 .- b2_y1 .+ 1) .* (b2_x2 .- b2_x1 .+ 1)
    unions = b1_area .+ b2_area .- intersection
    iou = intersection ./ unions
    iou = reshape(iou, boxes2_repeat[1], boxes1_repeat[1])
    return iou
end
=#

function clip_to_window(window, boxes)
    boxes[:, 1] = clamp.(boxes[:, 1], window[1], window[3])
    boxes[:, 2] = clamp.(boxes[:, 2], window[2], window[4])
    boxes[:, 3] = clamp.(boxes[:, 3], window[1], window[3])
    boxes[:, 4] = clamp.(boxes[:, 4], window[2], window[4])
    return boxes
end

#=
function clip_to_window2(window, boxes)
    @show window
    b1 = clamp.(boxes[:, 1], window[1], window[3])
    b2 = clamp.(boxes[:, 2], window[2], window[4])
    b3 = clamp.(boxes[:, 3], window[1], window[3])
    b4 = clamp.(boxes[:, 4], window[2], window[4])

    return hcat(b1, b2, b3, b4)
end
=#
#=
function refine_detections(rois, probs, deltas, window, config = nothing)
    # global grois = rois
    # global gprobs = probs
    # global gdeltas = deltas
    # global gwindow = window
    @show "in refine"
    @show size(probs)
    # gradient dropped here ?
    _, class_ids = findmax(Tracker.data(probs), dims = 1)
    @show size(deltas)
    @show size(probs)
    @show size(class_ids)
    # @show class_ids
    @warn "find max"
    idx = 1:length(class_ids)
    class_scores = vec(probs[class_ids])
    @show "maybe here"
    class_ids = vec(map(x -> x.I[1], class_ids))
    # @show idx
    # @show class_ids
    z = zip(idx, class_ids)
    # @show class_ids
    # @show idx
    # @show size(deltas)
    # @show size(rois)
    # @show size(probs)
    # deltas_specific = [deltas[x[1],:,x[2]] for x in z]
    # @show size(deltas_specific[1])
    # deltas_specific = reduce(hcat, deltas_specific)'
    # global gprobs = probs
    # global gclass_ids = class_ids
    # global gidx = idx
    # global gdeltas = deltas
    # deltas_specific = deltas[class_ids,:,idx]
    std_dev = [0.1 0.1 0.2 0.2] |> gpu
    @show "guess what"
    # deltas_specific .* std_dev
    deltas_specific = []
    # @show size(deltas_specific)
    for (m,k) in zip(idx, class_ids)
        push!(deltas_specific, deltas[k,:,m])
    end
    deltas_specific = cat(deltas_specific..., dims = 2)
    # @show size(deltas_specific' .* std_dev)
    # @show size(transpose(rois))
    refined_rois = apply_box_deltas(transpose(rois), transpose(deltas_specific) .* std_dev)
    @warn "applied deltas"
    height, width = 1024, 1024
    scale = [height width height width] |> gpu
    refined_rois = refined_rois .* scale
    refined_rois = clip_to_window2(window, refined_rois)
    @show typeof(refined_rois)
    refined_rois = round.(refined_rois)
    # @show refined_rois

    keep_bool = class_ids .> 0
    @show sum(keep_bool)
    DETECTION_MIN_CONFIDENCE = .5f0
    @show minimum(class_scores), mean(class_scores), maximum(class_scores)
    @show sum(class_scores .> DETECTION_MIN_CONFIDENCE)
    cs = copy(class_scores) |> cpu
    keep_bool = keep_bool .& (cs .>= DETECTION_MIN_CONFIDENCE)
    # @show sum(keep_bool), size(keep_bool)
    keep = findall(!iszero, keep_bool)
    @show keep
    pre_nms_class_ids = class_ids[keep]
    pre_nms_scores = class_scores[keep]
    # idxs = map(x -> x.I[1], findall(keep))
    pre_nms_rois = refined_rois[keep, :]

    unis = unique(pre_nms_class_ids) |> sort
    @warn "going to loop"
    nms_keep = Int[]
    for (i,class_id) in enumerate(unis)
	ixs = findall(pre_nms_class_ids .== class_id)

	ix_rois = pre_nms_rois[ixs, :]
	@show typeof(ix_rois)
	ix_scores = pre_nms_scores[ixs]
	order = sortperm(ix_scores, rev = true)
	ix_scores = ix_scores[order]
	ix_rois = ix_rois[order, :]
	# @show ix_rois
	DETECTION_NMS_THRESHOLD = 0.3f0
	class_keep = nms2(hcat(ix_rois .* 1.0f0, ix_scores), DETECTION_NMS_THRESHOLD)
	class_keep = keep[ixs[order[class_keep]]]
	if i==1
            push!(nms_keep, class_keep...)
        else
            push!(nms_keep, unique(class_keep)...)
        end
    end
    nms_keep = sort(nms_keep)
    keep = sort(intersect(keep, nms_keep))
    DETECTION_MAX_INSTANCES = 100
    roi_count = DETECTION_MAX_INSTANCES
	ends = min(roi_count, length(keep))
    top_ids = sortperm(class_scores[keep], rev = true)[1:ends]
    keep = keep[top_ids]

    @show size(refined_rois)
    @show size(class_ids)
    @show size(class_scores)

    # global grefined_rois = refined_rois
    # global gclass_ids = class_ids
    # global gclass_scores = class_scores
    # global gkeep = keep
    hcat(refined_rois[keep, :] .* 1.0f0, cu(class_ids[keep]) .* 1.0f0, class_scores[keep])

end
=#
#=
Loss Functions
=#

#=
function compute_rpn_class_loss(rpn_match, rpn_class_logits; labels = 1:80)
    anchor_class = Int.(rpn_match .== 1)
    indices = findall(!iszero, rpn_match .!= 0)

    rpn_class_logits = rpn_class_logits[:, indices, :]
    rpn_class_logits = dropdims(rpn_class_logits, dims = ndims(rpn_class_logits))
    anchor_class = anchor_class[indices]

    anchor_class = Flux.onehotbatch(anchor_class, 0:1) |> gpu
    Flux.logitcrossentropy(rpn_class_logits, anchor_class)
end
=#

function smooth_l1_loss(y, fx; δ = 1)
    α = abs(y - fx)
    abs(α) <= δ && return 0.5f0 * α ^ 2
    δ * α - (0.5f0 * δ ^ 2)
end

#=
huber_loss(y, ŷ; kwargs...) = smooth_l1_loss(y, ŷ, kwargs...)
=#
#=
function compute_rpn_bbox_loss(target_bbox, rpn_match, rpn_bbox)
    inds = findall(rpn_match .== 1)

    rpn_bbox = rpn_bbox[:, inds, :]
    rpn_bbox = dropdims(rpn_bbox, dims = ndims(rpn_bbox))

    target_bbox = target_bbox[1:size(rpn_bbox,2), :]
    target_bbox = transpose(target_bbox)

    # smooth L1 loss
    mean(smooth_l1_loss.(target_bbox, rpn_bbox))
end
=#
#=
function compute_mrcnn_class_loss(target_class_ids, pred_class_logits; labels = 0:80)
    if length(target_class_ids) > 0
	target_class_ids = Int.(target_class_ids)
	y = Flux.onehotbatch(target_class_ids, labels) |> gpu
	return Flux.logitcrossentropy(pred_class_logits, y)
    else
	return param(0.0f0)
    end
end
=#
#=
function compute_mrcnn_bbox_loss(target_deltas, target_class_ids, pred_bbox; labels = 0:80)
    if length(target_class_ids) > 0
	target_class_ids = map(y -> findall(x -> x == y, labels)[1], target_class_ids)
	positive_roi_ix = findall(target_class_ids .> 0)

	positive_roi_class_ids = target_class_ids[positive_roi_ix]
	target_deltas = target_deltas[positive_roi_ix, :]
	target_deltas = transpose(target_deltas)
	# pred_bbox = pred_bbox[positive_roi_class_ids, :, positive_roi_ix]
	bb = []
	for i = 1:length(positive_roi_ix)
	    a = pred_bbox[positive_roi_class_ids[i], :, i]
	    push!(bb, a)
	end
	pred_bbox = reduce(hcat, bb)

	mean(smooth_l1_loss.(pred_bbox, target_deltas))
    else
	return param(0.0f0)
    end
end
=#
#=
function compute_mrcnn_mask_loss(target_mask, target_class_ids, mrcnn_mask; labels = 0:80)
    if length(target_class_ids) > 0
	target_class_ids = map(y -> findall(x -> x == y, labels)[1], target_class_ids)
	positive_ix = findall(target_class_ids .> 0.f0)
	positive_class_ids = target_class_ids[positive_ix]

	y_true = target_mask[:,:,positive_ix]
	bb = []
	for i = 1:length(positive_ix)
	    a = mrcnn_mask[:,:,positive_class_ids[i], i]
	    push!(bb, a)
	end
	y_pred = cat(bb..., dims = 3)
		# y_pred = cpu(y_pred)
	# y_true = cpu(y_true)

	mean(bce(y_pred, y_true))
    else
	return param(0.0f0)
    end
end
=#
##################
# Data Generator #
##################

function mold_image(image, config = Nothing)
    MEAN_PIXEL = [123.7f0 / 255.0f0, 116.8f0 / 255.0f0, 103.9f0 / 255.0f0]
    # MEAN_PIXEL = [123.7f0, 116.8f0, 103.9]
    image[:,:,1] .-= MEAN_PIXEL[1]
    image[:,:,2] .-= MEAN_PIXEL[2]
    image[:,:,3] .-= MEAN_PIXEL[3]
    return image
end

function compose_image_meta(image_id, image_shape, window, active_class_ids)
    return ([image_id], image_shape, window, active_class_ids)
end

parse_image_meta(image_meta) = return image_meta

function mold_inputs(images, config = Nothing)
    IMAGE_MIN_DIM = 800 # just for testing
    IMAGE_MAX_DIM = 1024 # 128
    IMAGE_PADDING = true # just for testing
    NUM_CLASSES = 81 # just for testing - coco
    molded_images = []
    image_metas = []
    windows = []
    for image in images
	molded_image, window, scale, padding = resize_image(
	    image,
            min_dim = IMAGE_MIN_DIM,
            max_dim = IMAGE_MAX_DIM,
            padding = IMAGE_PADDING)

	molded_image = mold_image(molded_image)
	image_meta = compose_image_meta(0, size(image), window, zeros(Int32, NUM_CLASSES))
	push!(molded_images, molded_image)
	push!(windows, window)
	push!(image_metas, image_meta)
    end

    molded_images = cat(molded_images..., dims = 4)
    windows = reduce(hcat, windows)
    image_metas = reduce(hcat, image_metas)
    return molded_images, image_metas, windows
end


text:
export BrickletIndustrialDualACRelayValue
struct BrickletIndustrialDualACRelayValue
    channel0::Bool
    channel1::Bool
end

export BrickletIndustrialDualACRelayMonoflop
struct BrickletIndustrialDualACRelayMonoflop
    value::Bool
    time::Integer
    time_remaining::Integer
end

export BrickletIndustrialDualACRelaySPITFPErrorCount
struct BrickletIndustrialDualACRelaySPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletIndustrialDualACRelayIdentity
struct BrickletIndustrialDualACRelayIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletIndustrialDualACRelay
"""
Two relays to switch AC devices
"""
mutable struct BrickletIndustrialDualACRelay <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletIndustrialDualACRelay(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_industrial_dual_ac_relay")
        deviceInternal = package.BrickletIndustrialDualACRelay(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export set_value
"""
    $(SIGNATURES)

Sets the state of the relays, *true* means on and *false* means off.
For example: (true, false) turns relay 0 on and relay 1 off.

If you just want to set one of the relays and don't know the current state
of the other relay, you can get the state with :func:`Get Value` or you
can use :func:`Set Selected Value`.

All running monoflop timers will be aborted if this function is called.
"""
function set_value(device::BrickletIndustrialDualACRelay, channel0, channel1)
    device.deviceInternal.set_value(channel0, channel1)
end

export get_value
"""
    $(SIGNATURES)

Returns the state of the relays, *true* means on and *false* means off.
"""
function get_value(device::BrickletIndustrialDualACRelay)
    return device.deviceInternal.get_value()
end

export set_channel_led_config
"""
    $(SIGNATURES)

Each channel has a corresponding LED. You can turn the LED off, on or show a
heartbeat. You can also set the LED to "Channel Status". In this mode the
LED is on if the channel is high and off otherwise.
"""
function set_channel_led_config(device::BrickletIndustrialDualACRelay, channel, config)
    device.deviceInternal.set_channel_led_config(channel, config)
end

export get_channel_led_config
"""
    $(SIGNATURES)

Returns the channel LED configuration as set by :func:`Set Channel LED Config`
"""
function get_channel_led_config(device::BrickletIndustrialDualACRelay, channel)
    return device.deviceInternal.get_channel_led_config(channel)
end

export set_monoflop
"""
    $(SIGNATURES)

The first parameter can be 0 or 1 (relay 0 or relay 1). The second parameter
is the desired state of the relay (*true* means on and *false* means off).
The third parameter indicates the time that the relay should hold
the state.

If this function is called with the parameters (1, true, 1500):
Relay 1 will turn on and in 1.5s it will turn off again.

A monoflop can be used as a failsafe mechanism. For example: Lets assume you
have a RS485 bus and a Industrial Dual AC Relay Bricklet connected to one of the
slave stacks. You can now call this function every second, with a time parameter
of two seconds. The relay will be on all the time. If now the RS485
connection is lost, the relay will turn off in at most two seconds.
"""
function set_monoflop(device::BrickletIndustrialDualACRelay, channel, value, time)
    device.deviceInternal.set_monoflop(channel, value, time)
end

export get_monoflop
"""
    $(SIGNATURES)

Returns (for the given relay) the current state and the time as set by
:func:`Set Monoflop` as well as the remaining time until the state flips.

If the timer is not running currently, the remaining time will be returned
as 0.
"""
function get_monoflop(device::BrickletIndustrialDualACRelay, channel)
    return device.deviceInternal.get_monoflop(channel)
end

export set_selected_value
"""
    $(SIGNATURES)

Sets the state of the selected relay, *true* means on and *false*
means off.

A running monoflop timer for the selected relay will be aborted if this function
is called.

The other relay remains untouched.
"""
function set_selected_value(device::BrickletIndustrialDualACRelay, channel, value)
    device.deviceInternal.set_selected_value(channel, value)
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletIndustrialDualACRelay)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletIndustrialDualACRelay, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletIndustrialDualACRelay)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletIndustrialDualACRelay, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletIndustrialDualACRelay, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletIndustrialDualACRelay, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletIndustrialDualACRelay)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletIndustrialDualACRelay)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletIndustrialDualACRelay)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletIndustrialDualACRelay, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletIndustrialDualACRelay)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletIndustrialDualACRelay)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletIndustrialDualACRelay, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
using LinearAlgebra, StaticArrays
I4 = SMatrix{4,4}(I)

# Define your noisy linear system
# dF(t)/dt = B*F(T) + β + σ*Noise
B = SMatrix{4,4,Float64}([-0.1 0 1.0 0 ; 0 -0.1 0 1.0; 0 0 -1 8.0; 0 0 -8.0 -1])
σ = SMatrix{4,4}(Diagonal([0.001, 0.001, 1.0, 1.0]))
β = SVector(0.0, 0.0, 0.0, 0.0)
plin = (B, β, σ)

import MitosisStochasticDiffEq as MSDE
using Mitosis

# define SDE function
f(u,p,t) = p[1]*u + p[2]
g(u,p,t) = diag(p[3])

# time span
tstart = 0.0
tend = 5.0
dt = 0.001
trange = tstart:dt:tend

# Start point
u0 = SVector(0.0, 1.0, 1.0, 0.0)

# End point
uT = u0 # loop

sdekernel = MSDE.SDEKernel(f,g,trange,plin)
kerneltilde = MSDE.SDEKernel(Mitosis.AffineMap(B, β), Mitosis.ConstantMap(σ), trange, plin)

WG = WGaussian{(:μ,:Σ,:c)}(uT, SMatrix{4,4}(Diagonal([0.0001, 0.0001, 0.0001, 0.0001].^(-2.0))), 0.0) # move back from endpoint plus a bit of uncertainty
WG = (;logscale = 0.0, μ=uT, Σ=SMatrix{4,4}(Diagonal([0.0000001, 0.0000001, 0.0000001, 0.0000001].^(2.0))))

message, solend = MSDE.backwardfilter(kerneltilde, WG)

ts, u, uend, noise, ll = MSDE.forwardguiding(sdekernel, message, (u0, 0.0),
   Z=nothing; inplace=false, save_noise=true)

using Makie
pl = lines(getindex.(u, 1), getindex.(u, 2))
#lines!(pl, getindex.(u, 3), getindex.(u, 4), color=:blue)
#scatter!(pl, [u[1][1], u[1][3], u[end][1], u[end][3]], [u[1][2], u[1][4], u[end][2], u[end][4]]  )
#save("curve.png", pl)
pl


text:
export tanhActivation, tanhActivation!

"""
 hyperbolic tan activation A = tanh(Y)

 Input:

   Y - array of features

 Optional Input:

   doDerivative - flag for computing derivative, set via varargin
                  Ex: tanhActivation(Y,true);

 Output:

  A  - activation
  dA - derivatives
"""
function tanhActivation(Y::Array{T,2},doDerivative::Bool=false) where {T <: Number}

    A = tanh.(Y)
    dA = zeros(A)
    if doDerivative
         dA .= one(T) .- A.^2
    end
    return A, dA
end

function tanhActivation!(A::Array{T,2},dA=[],doDerivative::Bool=false) where {T <: Number}

    A .= tanh.(A)
    if doDerivative
        if isempty(dA)
            dA = one(T) .- A.^2
        else
            dA .= one(T) .- A.^2
        end
    end
    return A, dA
end


text:
@testset "I474" begin
    @test isempty(messageset(lintstr("""
    function f(x)
        d = Dict{Symbol,Int}(:a=>1, :b=>2)
        for i in d
            @show i
        end
        return x
    end
    """)))

    @test messageset(lintstr("""
    function f(x)
        d = Dict{Symbol,Int}(:a=>1, :b=>2)
        for (k,) in d
            @show k
        end
        return x
    end
    """)) == Set([:I474])

    @test isempty(messageset(lintstr("""
    function f(a::Array{Int,1})
        for i in enumerate(a)
            println(i)
        end
    end
    """)))

    @test messageset(lintstr("""
    function f(x)
        d = [(1,2,3), (4,5,6)]
        for (i, j) in d
            @show i j
        end
        return x
    end
    """)) == Set([:I474])
end