text:
# Behavior cloning where the policy is approximated by a DNN

"""
    normalize(X::AbstractMatrix) -> (Matrix, UnitRangeTransform)

Min-max normalize the input matrix `X`. The normalized version and the associated `UnitRangeTransform`
object are returned.
"""
function normalize(X::AbstractMatrix)
    urt = StatsBase.fit(UnitRangeTransform, X; dims = 2)
    return StatsBase.transform(urt, X), urt
end

function shuffle_split_obs(X, y; at = 0.9)
    return splitobs(shuffleobs((X, y)); at)
end

struct Agent{T<:Flux.Chain}
    urt::UnitRangeTransform # min-max scaling
    model::T
end

function (agent::Agent)(X)
    Xn = StatsBase.transform(agent.urt, X)
    return agent.model(Xn)
end

"""
    eval_loss(loader, model, l; device=cpu) -> Float64

Compute the total loss of data in `loader` using the DNN `model` assessed by the `loss` function.
The aggregator of `loss` should be `mean`. The input data in `loader` should be normalized already.
"""
function eval_loss(loader, model, loss; device = Flux.cpu)::Float64
    n = 0
    tl = 0.0
    model = model |> device
    for (x, y) in loader
        x, y = x |> device, y |> device
        ŷ = model(x)
        tl += loss(ŷ, reshape(y, size(ŷ)...)) * length(y)
        n += length(y)
    end
    return tl / n
end

"""
    eval_loss(loader, agent::Agent, loss; device = Flux.cpu) -> Float64

Compute the total loss of data in `loader` using the DNN `agent` assessed by the `loss` function.
The aggregator of `loss` should be `mean`. The input data in `loader` will get first normalized by
the transform in `agent`.
"""
function eval_loss(loader, agent::Agent, loss; device = Flux.cpu)::Float64
    urt = agent.urt
    norm_loader = [(StatsBase.transform(urt, X), y) for (X, y) in loader]
    return eval_loss(norm_loader, agent.model, loss; device)
end

"""

Train a neural network `model`. The weights of the `model` are updated in place. The loss history of
the training set and the validation set during training is returned. At last, the `model` has parameters
corresponding to the minimum validation loss.

Note that the input `X` should have already been normalized properly.
`loss(ŷ, y)`
"""
function trainNN!(model, X, y, loss; val_ratio = 0.1,
    η = 1e-2, batchsize = 256, epochs = 500, device = Flux.gpu,
    slow_learning_threshold = 10, early_stopping_threshold = 20,
    report_freq = 5)

    (X_train, Y_train), (X_val, Y_val) = shuffle_split_obs(X, y; at = 1 - val_ratio)
    loader_train = Flux.Data.DataLoader((X_train, Y_train); batchsize, shuffle = true)
    loader_val = Flux.Data.DataLoader((X_val, Y_val); batchsize, shuffle = false)

    m = model |> device
    ps = Flux.params(m)
    opt = Flux.Optimise.ADAM(η)

    best_loss_val = Inf
    best_ps = nothing
    last_improvement = 0
    loss_history = zeros(epochs, 2)

    loss_train = eval_loss(loader_train, m, loss; device)
    loss_val = eval_loss(loader_val, m, loss; device)
    # loss_history[e, :] .= (loss_train, loss_val)
    report_freq > 0 && @info 0 loss_train loss_val

    for e = 1:epochs
        for (x, y) in loader_train  # mini_batch
            x, y = device(x), device(y)
            gs = gradient(ps) do
                ŷ = m(x)
                loss(ŷ, reshape(y, size(ŷ)...))
            end
            Flux.update!(opt, ps, gs)
        end

        # calculate all at once since the dataset is small
        loss_train = eval_loss(loader_train, m, loss; device)
        loss_val = eval_loss(loader_val, m, loss; device)
        loss_history[e, :] .= (loss_train, loss_val)
        if report_freq > 0 && e % report_freq == 0
            @info e loss_train loss_val
        end

        if loss_val < best_loss_val
            last_improvement = e
            best_loss_val = loss_val
            # @save model_file model = cpu(m)
            best_ps = deepcopy(m |> cpu |> params)
        end
        # learning rate schedule
        if e - last_improvement >= slow_learning_threshold && η > 1e-5
            η /= 1.5
            @info "No improvement in $(slow_learning_threshold) epochs. Drop learning rate to" η
            last_improvement = e
        end
        if e - last_improvement >= early_stopping_threshold
            @info "No improvement in $(early_stopping_threshold) epochs. Early stopped." best_loss_val
            loss_history = loss_history[1:e, :]
            break
        end
    end

    # write the optimized parameters into `model`
    @info "The minimum validation loss is $best_loss_val"
    Flux.loadparams!(model, best_ps)
    return loss_history
end

"""
    build_SU_net(tf=relu) -> Chain

State: ``[t, ρ, P^{PV}]``. Action: ``z^{su} \\in {0, 1}``
"""
build_SU_net(in; tf = relu) = Chain(
    Dense(in, 200, tf),
    Dropout(0.1),
    Dense(200, 100, tf),
    Dropout(0.3),
    Dense(100, 50, tf),
    Dropout(0.5),
    Dense(50, 1, sigmoid)
)

build_SU_net2(in; tf = relu) = Chain(
    Dense(in, 100, tf),
    Dense(100, 50, tf),
    # Dense(100, 50, tf),
    Dense(50, 1, sigmoid)
)

"""
    weighted_binarycrossentropy(ŷ, y; w1=1.0)

Add weight to Flux `binarycrossentropy` loss. The weight of class 0 is fixed to 1.0.
"""
function weighted_binarycrossentropy(ŷ, y; w1 = 1.0)
    # based on the source of Flux.binarycrossentropy
    ϵ = Flux.epseltype(ŷ)
    mean(@.(-w1 * Flux.Losses.xlogy(y, ŷ + ϵ) - Flux.Losses.xlogy(1 - y, 1 - ŷ + ϵ)))
end

function train_SU_agent(X::AbstractMatrix, y::AbstractVector; tf = relu, w1 = NaN, kwargs...)::Tuple{Agent,Matrix}
    Xn, urt = normalize(X)

    if isnothing(w1)
        loss = (ŷ, y) -> Flux.binary_focal_loss(ŷ, y; γ = 3.5)
        @show "gamma"
    else
        if isnan(w1)
            # imbalanced -> reweight. Fix w0 = 1.
            cm = countmap(y)
            w1 = cm[0] / cm[1]
        end
        loss = (ŷ, y) -> weighted_binarycrossentropy(ŷ, y; w1)
    end
    model = build_SU_net2(size(X)[1]; tf)
    hist = trainNN!(model, Xn, y, loss; kwargs...)
    return (Agent(urt, model), hist)
end

function eval_accuracy(loader, agent::Agent; device = cpu)
    n_acc = n_all = 0
    model = agent.model |> device
    for (x, y) in loader
        x, y = StatsBase.transform(agent.urt, x) |> device, y |> device
        ŷ = model(x) .> 0.5
        n_acc += sum(y .== reshape(ŷ, size(y)...))
        n_all += length(y)
        @show n_acc n_all
    end
    return n_acc / n_all
end

max_stars_repo_name: ShuhuaGao/HybridHEMS

max_issues_repo_name: ShuhuaGao/HybridHEMS

max_forks_repo_name: ShuhuaGao/HybridHEMS


text:
export AbstractEnvWrapper

abstract type AbstractEnvWrapper <: AbstractEnv end

Base.nameof(env::AbstractEnvWrapper) = "$(nameof(env.env)) |> $(nameof(typeof(env)))"

Base.getindex(env::AbstractEnvWrapper) = env.env

(env::AbstractEnvWrapper)(args...; kwargs...) = env.env(args...; kwargs...)

for f in vcat(RLBase.ENV_API, RLBase.MULTI_AGENT_ENV_API)
    @eval RLBase.$f(x::AbstractEnvWrapper, args...; kwargs...) = $f(x[], args...; kwargs...)
end

# avoid ambiguous
RLBase.state(env::AbstractEnvWrapper, ss::RLBase.AbstractStateStyle, p) =
    state(env[], ss, p)
RLBase.state(env::AbstractEnvWrapper, ss::RLBase.AbstractStateStyle) = state(env[], ss)
RLBase.state_space(env::AbstractEnvWrapper, ss::RLBase.AbstractStateStyle) =
    state_space(env[], ss)
RLBase.state_space(env::AbstractEnvWrapper, ss::RLBase.AbstractStateStyle, p) =
    state_space(env[], ss, p)

include("ActionTransformedEnv.jl")
include("DefaultStateStyle.jl")
include("MaxTimeoutEnv.jl")
include("RewardOverriddenEnv.jl")
include("StateCachedEnv.jl")
include("StateTransformedEnv.jl")
include("StochasticEnv.jl")

ext: jl


text:
# Julia wrapper for header: /usr/local/ffmpeg/2.4/include/libavutil/xtea.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_xtea_init,
    av_xtea_crypt

function av_xtea_init(ctx,key)
    ccall((:av_xtea_init,libavutil),Void,(Ptr{AVXTEA},Ptr{Uint8}),ctx,key)
end

function av_xtea_crypt(ctx,dst,src,count::Integer,iv,decrypt::Integer)
    ccall((:av_xtea_crypt,libavutil),Void,(Ptr{AVXTEA},Ptr{Uint8},Ptr{Uint8},Cint,Ptr{Uint8},Cint),ctx,dst,src,count,iv,decrypt)
end


text:
# This is very simple implementation of go-style slices. Nodes in the B-tree
# have fixed sized arrays that change when elements are inserted or deleted.
# We would like for those to behave as though they are resized without actually
# doing any reallocation.

mutable struct Slice{T, N} <: AbstractVector{T}
    data::Vector{T}
    n::Int # Number of stored elements

    function Slice{T,N}() where {T,N}
        return new{T,N}(Vector{T}(undef, N), 0)
    end
end

function Base.length(s::Slice)
    return s.n
end

function Base.size(s::Slice)
    return (s.n,)
end

@inline function Base.getindex(s::Slice{T, N}, i::Integer) where {T, N}
    @boundscheck checkbounds(s, i)
    @inbounds return s.data[i]
end

@inline function Base.setindex!(s::Slice{T, N}, value, i::Integer) where {T, N}
    @boundscheck checkbounds(s, i)
    @inbounds s.data[i] = value
    return s
end

function Base.push!(s::Slice{T, N}, value) where {T, N}
    if s.n < N
        s.n += 1
        @inbounds s.data[s.n] = value
    else
        throw(BoundsError("Slice expanded past fixed size."))
    end
end

function Base.pop!(s::Slice{T, N}) where {T, N}
    if s.n > 0
        @inbounds x = s.data[s.n]
        s.n -= 1
        return x
    else
        throw(BoundsError())
    end
end

function Base.insert!(s::Slice{T, N}, i::Integer, value) where {T, N}
    if s.n < N && 1 <= i <= s.n + 1
        # TODO: This should work but won't. Fix this in Julia.
        #copy!(s.data, i+1, s.data, i, s.n - i + 1)
        if isbitstype(T)
            unsafe_copyto!(pointer(s.data, i+1), pointer(s.data, i), s.n - i + 1)
        else
            for k in 0:(s.n - i)
                @inbounds s.data[s.n+1-k] = s.data[s.n-k]
            end
        end
        @inbounds s.data[i] = value
        s.n += 1
        return s
    else
        throw(BoundsError())
    end
end

function Base.resize!(s::Slice{T, N}, n::Integer) where {T, N}
    if 1 <= n <= N
        s.n = n
        return s
    else
        throw(BoundsError())
    end
end

function Base.splice!(s::Slice{T, N}, i::Integer) where {T, N}
    if 1 <= i <= s.n
        @inbounds x = s.data[i]
        for j in i:s.n-1
            @inbounds s.data[j] = s.data[j + 1]
        end
        s.n -= 1
        return x
    else
        throw(BoundsError())
    end
end

function Base.searchsortedfirst(s::Slice, x)
    return searchsortedfirst(s.data, x, 1, s.n, Base.Order.Forward)
end

function slice_insert!(xs::Vector{T}, count::Integer, i::Integer, value::T) where T
    if count < length(xs) && 1 <= i <= count + 1
        if isbitstype(T)
            unsafe_copyto!(pointer(xs, i+1), pointer(xs, i), count - i + 1)
        else
            for k in 0:(count - i)
                @inbounds xs[count+1-k] = xs[count-k]
            end
        end
        @inbounds xs[i] = value
        count += 1
        return count
    else
        throw(BoundsError())
    end
end

function slice_splice!(xs::Vector{T}, count::Integer, i::Integer) where T
    if 1 <= i <= count
        @inbounds x = xs[i]
        for j in i:count-1
            @inbounds xs[j] = xs[j + 1]
        end
        count -= 1
        return x, count
    else
        throw(BoundsError())
    end
end


text:
module TestSorting

using Test
using RLEVectors

@testset begin

    # issorted
    @test issorted(RLEVector([3, 4, 2], [3, 6, 9])) == false
    @test issorted(RLEVector([2, 3, 4], [3, 6, 9])) == true
    @test issorted(RLEVector(["bob", "fred", "joe"], [3, 6, 9])) == true

    # sort, sort!
    x = RLEVector([2, 4, 3, 5], [4, 7, 9, 10])
    @test sort(x) == RLEVector([2, 3, 4, 5], [4, 6, 9, 10])
    @test sort!(x) == RLEVector([2, 3, 4, 5], [4, 6, 9, 10])

    # sortperm
    x = RLEVector([2, 4, 3, 5], [4, 7, 9, 10])
    @test sortperm(x) == RLEVector([1, 3, 2, 4], [4, 6, 9, 10])

    # reverse
    x = RLEVector([2, 4, 3, 5], [4, 7, 9, 10])
    @test reverse(x) == RLEVector([5, 3, 4, 2], [1, 3, 6, 10])
    @test reverse!(x) == RLEVector([5, 3, 4, 2], [1, 3, 6, 10])

    # permute runs
    x = RLEVector([2, 4, 3, 5], [4, 7, 9, 10])
    y = RLEVector([4, 5, 3, 2], [3, 4, 6, 10])
    @test permute_runs(x, [2, 4, 3, 1]) == y

end # testset

end # module


text:
using TankDemo
using Base.Test

tank = Tank(0.0, 0.0, [Tank(100, 0.1, nothing), Tank(50, 0.2, nothing)])
fill_tank!(tank)

@test tank.level == 100*0.1 + 50*0.2
@test tank.tank_up[1].level == 100 - 10
@test tank.tank_up[2].level == 50 - 10

tank = Tank(100.0, 0.05, nothing)       # tank with no upstream tanks
fill_tank!(tank)

@test tank.level == 100                 # nothing should change


text:
## Result - data structure interaction

# These functions allow to pass investment decisions between two problems: the investment decision
# results of the first problem are copied into the data structure of the second problem; an
# appropriate constraint could be necessary in the model of the second problem to read the data
# prepared by these functions.

function add_ne_branch_indicator!(mn_data::Dict{String,Any}, result::Dict{String,Any})
    for (n, data_nw) in mn_data["nw"]
        sol_nw = result["solution"]["nw"][n]
        for (b, data_branch) in data_nw["ne_branch"]
            if data_branch["br_status"] == 1
                data_branch["built"] = sol_nw["ne_branch"][b]["built"]
            end
        end
    end
end

function add_ne_storage_indicator!(mn_data::Dict{String,Any}, result::Dict{String,Any})
    for (n, data_nw) in mn_data["nw"]
        sol_nw = result["solution"]["nw"][n]
        for (s, data_storage) in data_nw["ne_storage"]
            data_storage["isbuilt"] = sol_nw["ne_storage"][s]["isbuilt"]
        end
    end
end

function add_td_coupling_power_active!(mn_data::Dict{String,Any}, result::Dict{String,Any})
    for (n, data_nw) in mn_data["nw"]
        p = result["solution"]["nw"][n]["td_coupling"]["p"]
        d_gen_id = _FP.dim_prop(mn_data, parse(Int,n), :sub_nw, "d_gen")
        d_gen = data_nw["gen"]["$d_gen_id"] = deepcopy(data_nw["gen"]["$d_gen_id"]) # Gen data is shared among nws originally.
        d_gen["pmax"] = p
        d_gen["pmin"] = p
    end
end

function add_max_cost!(mn_data::Dict{String,Any}, max_cost::Float64)
    mn_data["max_cost"] = max_cost
end


text:
module Agents

using Requires
using Distributed
using DataStructures
using LightGraphs
using DataFrames
using Random
import ProgressMeter

import Base.length
import LinearAlgebra

# Core structures of Agents.jl
include("core/agents.jl")
include("core/model.jl")
include("core/space_interaction_API.jl")

# Existing spaces
include("spaces/nothing.jl")
include("spaces/graph.jl")
include("spaces/grid.jl")
include("spaces/discrete.jl")
include("spaces/continuous.jl")
include("spaces/openstreetmap.jl")
include("spaces/utilities.jl")

# Stepping and data collection functionality
include("simulations/step.jl")
include("simulations/collect.jl")
include("simulations/paramscan.jl")
include("simulations/sample.jl")
include("simulations/ensemblerun.jl")

# Other features that exist in submodules
include("submodules/pathfinding/all_pathfinders.jl")
include("submodules/schedulers.jl")

function __init__()
    # Plot recipes
    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
        include("visualization/plot-recipes.jl")
    end
    # Workaround for Documenter.jl, so we don't need to include
    # heavy dependencies to build documentation
    @require Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" begin
        include("visualization/plot-recipes.jl")
    end
end

# Deprecations, that will be removed in future versions
include("deprecated.jl")

# Predefined models
include("models/Models.jl")
export Models

# Update message:
display_update = true
version_number = "4.2"
update_name = "update_v$(version_number)"

if display_update
    if !isfile(joinpath(@__DIR__, update_name))
        printstyled(
            stdout,
            """
            \nUpdate message: Agents v$(version_number)
            Please see the changelog online. Some key features:

            * Self-contained features of Agents.jl will from now own exist in their own
              submodules. In the future, more features will be in submodules like this.
            * Full support for pathfinding, using the A* algorithm, in `GridSpace`
            * Several names have been reworked for more clarity (with deprecations)
            * New function `ensemblerun!` which replaces using `replicates` in `run!`
            * New documentation page "Performance Tips"
            * Parallelization of `paramscan`
            """;
            color = :light_magenta,
        )
        touch(joinpath(@__DIR__, update_name))
    end
end

end # module


text:
function inverseCantorPairing(z)
	w = floor((sqrt(8*z+1)-1)/2);
	t = w*(w+1)/2;
	y = z - t;
	x = w - y;
	return Int64[x ,y];
end

max_stars_repo_name: driesverhees/dobble

max_issues_repo_name: driesverhees/dobble

max_forks_repo_name: driesverhees/dobble


text:
# Pointers with address space information

#
# Address spaces
#

export AS, addrspace

abstract type AddressSpace end

module AS

import ..AddressSpace

struct Generic  <: AddressSpace end
struct Global   <: AddressSpace end
struct Shared   <: AddressSpace end
struct Constant <: AddressSpace end
struct Local    <: AddressSpace end

end

#
# Device pointer
#

"""
    DevicePtr{T,A}

A memory address that refers to data of type `T` that is accessible from the GPU. It is the
on-device counterpart of `CUDAdrv.CuPtr`, additionally keeping track of the address space
`A` where the data resides (shared, global, constant, etc). This information is used to
provide optimized implementations of operations such as `unsafe_load` and `unsafe_store!.`
"""
DevicePtr

if sizeof(Ptr{Cvoid}) == 8
    primitive type DevicePtr{T,A} <: Ref{T} 64 end
else
    primitive type DevicePtr{T,A} <: Ref{T} 32 end
end

# constructors
DevicePtr{T,A}(x::Union{Int,UInt,CuPtr,DevicePtr}) where {T,A<:AddressSpace} = Base.bitcast(DevicePtr{T,A}, x)
DevicePtr{T}(ptr::CuPtr{T}) where {T} = DevicePtr{T,AS.Generic}(ptr)
DevicePtr(ptr::CuPtr{T}) where {T} = DevicePtr{T,AS.Generic}(ptr)

## getters

Base.eltype(::Type{<:DevicePtr{T}}) where {T} = T

addrspace(x::DevicePtr) = addrspace(typeof(x))
addrspace(::Type{DevicePtr{T,A}}) where {T,A} = A

## conversions

# to and from integers
## pointer to integer
Base.convert(::Type{T}, x::DevicePtr) where {T<:Integer} = T(UInt(x))
## integer to pointer
Base.convert(::Type{DevicePtr{T,A}}, x::Union{Int,UInt}) where {T,A<:AddressSpace} = DevicePtr{T,A}(x)
Int(x::DevicePtr)  = Base.bitcast(Int, x)
UInt(x::DevicePtr) = Base.bitcast(UInt, x)

# between host and device pointers
Base.convert(::Type{CuPtr{T}},  p::DevicePtr)  where {T}                 = Base.bitcast(CuPtr{T}, p)
Base.convert(::Type{DevicePtr{T,A}}, p::CuPtr) where {T,A<:AddressSpace} = Base.bitcast(DevicePtr{T,A}, p)
Base.convert(::Type{DevicePtr{T}}, p::CuPtr)   where {T}                 = Base.bitcast(DevicePtr{T,AS.Generic}, p)

# between CPU pointers, for the purpose of working with `ccall`
Base.unsafe_convert(::Type{Ptr{T}}, x::DevicePtr{T}) where {T} = reinterpret(Ptr{T}, x)

# between device pointers
Base.convert(::Type{<:DevicePtr}, p::DevicePtr)                         = throw(ArgumentError("cannot convert between incompatible device pointer types"))
Base.convert(::Type{DevicePtr{T,A}}, p::DevicePtr{T,A})   where {T,A}   = p
Base.unsafe_convert(::Type{DevicePtr{T,A}}, p::DevicePtr) where {T,A}   = Base.bitcast(DevicePtr{T,A}, p)
## identical addrspaces
Base.convert(::Type{DevicePtr{T,A}}, p::DevicePtr{U,A}) where {T,U,A} = Base.unsafe_convert(DevicePtr{T,A}, p)
## convert to & from generic
Base.convert(::Type{DevicePtr{T,AS.Generic}}, p::DevicePtr)               where {T}     = Base.unsafe_convert(DevicePtr{T,AS.Generic}, p)
Base.convert(::Type{DevicePtr{T,A}}, p::DevicePtr{U,AS.Generic})          where {T,U,A} = Base.unsafe_convert(DevicePtr{T,A}, p)
Base.convert(::Type{DevicePtr{T,AS.Generic}}, p::DevicePtr{T,AS.Generic}) where {T}     = p  # avoid ambiguities
## unspecified, preserve source addrspace
Base.convert(::Type{DevicePtr{T}}, p::DevicePtr{U,A}) where {T,U,A} = Base.unsafe_convert(DevicePtr{T,A}, p)

## limited pointer arithmetic & comparison

isequal(x::DevicePtr, y::DevicePtr) = (x === y) && addrspace(x) == addrspace(y)
isless(x::DevicePtr{T,A}, y::DevicePtr{T,A}) where {T,A<:AddressSpace} = x < y

Base.:(==)(x::DevicePtr, y::DevicePtr) = UInt(x) == UInt(y) && addrspace(x) == addrspace(y)
Base.:(<)(x::DevicePtr,  y::DevicePtr) = UInt(x) < UInt(y)
Base.:(-)(x::DevicePtr,  y::DevicePtr) = UInt(x) - UInt(y)

Base.:(+)(x::DevicePtr, y::Integer) = oftype(x, Base.add_ptr(UInt(x), (y % UInt) % UInt))
Base.:(-)(x::DevicePtr, y::Integer) = oftype(x, Base.sub_ptr(UInt(x), (y % UInt) % UInt))
Base.:(+)(x::Integer, y::DevicePtr) = y + x

## memory operations

Base.convert(::Type{Int}, ::Type{AS.Generic})  = 0
Base.convert(::Type{Int}, ::Type{AS.Global})   = 1
Base.convert(::Type{Int}, ::Type{AS.Shared})   = 3
Base.convert(::Type{Int}, ::Type{AS.Constant}) = 4
Base.convert(::Type{Int}, ::Type{AS.Local})    = 5

function tbaa_make_child(name::String, constant::Bool=false; ctx::LLVM.Context=JuliaContext())
    tbaa_root = MDNode([MDString("ptxtbaa", ctx)], ctx)
    tbaa_struct_type =
        MDNode([MDString("ptxtbaa_$name", ctx),
                tbaa_root,
                LLVM.ConstantInt(0, ctx)], ctx)
    tbaa_access_tag =
        MDNode([tbaa_struct_type,
                tbaa_struct_type,
                LLVM.ConstantInt(0, ctx),
                LLVM.ConstantInt(constant ? 1 : 0, ctx)], ctx)

    return tbaa_access_tag
end

tbaa_addrspace(as::Type{<:AddressSpace}) = tbaa_make_child(lowercase(String(as.name.name)))

@generated function Base.unsafe_load(p::DevicePtr{T,A}, i::Integer=1,
                                     ::Val{align}=Val(1)) where {T,A,align}
    eltyp = convert(LLVMType, T)

    T_int = convert(LLVMType, Int)
    T_ptr = convert(LLVMType, DevicePtr{T,A})

    T_actual_ptr = LLVM.PointerType(eltyp)

    # create a function
    param_types = [T_ptr, T_int]
    llvm_f, _ = create_function(eltyp, param_types)

    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        ptr = inttoptr!(builder, parameters(llvm_f)[1], T_actual_ptr)

        ptr = gep!(builder, ptr, [parameters(llvm_f)[2]])
        ptr_with_as = addrspacecast!(builder, ptr, LLVM.PointerType(eltyp, convert(Int, A)))
        ld = load!(builder, ptr_with_as)

        if A != AS.Generic
            metadata(ld)[LLVM.MD_tbaa] = tbaa_addrspace(A)
        end
        alignment!(ld, align)

        ret!(builder, ld)
    end

    call_function(llvm_f, T, Tuple{DevicePtr{T,A}, Int}, :((p, Int(i-one(i)))))
end

@generated function Base.unsafe_store!(p::DevicePtr{T,A}, x, i::Integer=1,
                                       ::Val{align}=Val(1)) where {T,A,align}
    eltyp = convert(LLVMType, T)

    T_int = convert(LLVMType, Int)
    T_ptr = convert(LLVMType, DevicePtr{T,A})

    T_actual_ptr = LLVM.PointerType(eltyp)

    # create a function
    param_types = [T_ptr, eltyp, T_int]
    llvm_f, _ = create_function(LLVM.VoidType(JuliaContext()), param_types)

    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        ptr = inttoptr!(builder, parameters(llvm_f)[1], T_actual_ptr)

        ptr = gep!(builder, ptr, [parameters(llvm_f)[3]])
        ptr_with_as = addrspacecast!(builder, ptr, LLVM.PointerType(eltyp, convert(Int, A)))
        val = parameters(llvm_f)[2]
        st = store!(builder, val, ptr_with_as)

        if A != AS.Generic
            metadata(st)[LLVM.MD_tbaa] = tbaa_addrspace(A)
        end
        alignment!(st, align)

        ret!(builder)
    end

    call_function(llvm_f, Cvoid, Tuple{DevicePtr{T,A}, T, Int},
                  :((p, convert(T,x), Int(i-one(i)))))
end

## loading through the texture cache

export unsafe_cached_load

# NOTE: CUDA 8.0 supports more caching modifiers, but those aren't supported by LLVM yet

# TODO: this functionality should throw <sm_32

# operand types supported by llvm.nvvm.ldg.global
const CachedLoadOperands = Union{UInt8, UInt16, UInt32, UInt64,
                                 Int8, Int16, Int32, Int64,
                                 Float32, Float64}

# containing DevicePtr types
const CachedLoadPointers = Union{Tuple(DevicePtr{T,AS.Global}
                                 for T in Base.uniontypes(CachedLoadOperands))...}

@generated function unsafe_cached_load(p::DevicePtr{T,AS.Global}, i::Integer=1,
                                       ::Val{align}=Val(1)) where
                                      {T<:CachedLoadOperands,align}
    # NOTE: we can't `ccall(..., llvmcall)`, because
    #       1) Julia passes pointer arguments as plain integers
    #       2) we need to addrspacecast the pointer argument

    eltyp = convert(LLVMType, T)

    T_int = convert(LLVMType, Int)
    T_int32 = LLVM.Int32Type(JuliaContext())
    T_ptr = convert(LLVMType, DevicePtr{T,AS.Global})

    T_actual_ptr = LLVM.PointerType(eltyp)
    T_actual_ptr_as = LLVM.PointerType(eltyp, convert(Int, AS.Global))

    # create a function
    param_types = [T_ptr, T_int]
    llvm_f, _ = create_function(eltyp, param_types)

    # create the intrinsic
    intrinsic_name = let
        class = if isa(eltyp, LLVM.IntegerType)
            :i
        elseif isa(eltyp, LLVM.FloatingPointType)
            :f
        else
            error("Cannot handle $eltyp argument to unsafe_cached_load")
        end
        width = sizeof(T)*8
        typ = Symbol(class, width)
        "llvm.nvvm.ldg.global.$class.$typ.p1$typ"
    end
    mod = LLVM.parent(llvm_f)
    intrinsic_typ = LLVM.FunctionType(eltyp, [T_actual_ptr_as, T_int32])
    intrinsic = LLVM.Function(mod, intrinsic_name, intrinsic_typ)

    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        ptr = inttoptr!(builder, parameters(llvm_f)[1], T_actual_ptr)

        ptr = gep!(builder, ptr, [parameters(llvm_f)[2]])
        ptr_with_as = addrspacecast!(builder, ptr, T_actual_ptr_as)
        ld = call!(builder, intrinsic,
                   [ptr_with_as, ConstantInt(Int32(align), JuliaContext())])

        metadata(ld)[LLVM.MD_tbaa] = tbaa_addrspace(AS.Global)

        ret!(builder, ld)
    end

    call_function(llvm_f, T, Tuple{DevicePtr{T,AS.Global}, Int}, :((p, Int(i-one(i)))))
end

@inline unsafe_cached_load(p::DevicePtr{T,AS.Global}, i::Integer=1, args...) where {T} =
    recurse_pointer_invocation(unsafe_cached_load, p+sizeof(T)*Int(i-one(i)),
                               CachedLoadPointers, 1, args...)


text:
"""
    startree_newick(n, l)

String for the Newick parenthetical description of the star tree
with `n` tips, and all branch lengths equal to `l`.
"""
function startree_newick(n::Integer, ell=1.0)
    return "(" * join(["t$i:$ell" for i=1:n], ",") * ");"
end

"""
    symmetrictree_newick(n::Int, ell::Real, i=1)

String for the Newick parenthetical description of a symmetric tree with
2^n tips, numbered from i to i+2^n-1. All branch lengths are set equal to `ell`.
The tree can be created later by reading the string with [`readTopology`](@ref).
"""
function symmetrictree_newick(n::Int, ell::Real, i=1::Int)
    tree = "A$(i-1+2^n):$(ell)"
    if n==0 return("("*"A$(i-1+2^n):0"*");") end
    for k in 1:(n-1)
        tree = "(" * tree * "," * tree * "):$(ell)"
    end
    tree = "(" * tree * "," * tree * ");"
    # Rename tips
    for k in (2^n-1):-1:1
        tree = replace(tree, "A$(i-1+k+1):$(ell)"=>"A$(i-1+k):$(ell)", count=k)
    end
    return(tree)
end

"""
    symmetricnet_newick(n::Int, i::Int, j::Int, γ::Real, l::Real)

Newick string for a network with a symmetric major tree with 2^n tips,
numbered from 1 to 2^n. All the branch lengths of the major tree are set to `l`.
One hybrid branch, going from level i to level j is added, cutting in half
each initial edge in the tree. The new edge has length `l` and inheritance `γ`.
"""
function symmetricnet_newick(n::Int, i::Int, j::Int, gamma::Real, ell::Real)
    (n < i || i < j || j < 1) && error("must have n >= i >= j > 0")
    # Underlying tree
    tree = symmetrictree_newick(n, ell)
    ## start hyb
    op = "(A1:$(ell)"
    clo = "A$(2^(i-1)):$(ell)):$(ell)"
    for k in 3:i
        op = "("*op
        clo = clo*"):$(ell)"
    end
    clobis = clo[1:(length(clo)-length("$(ell)"))]*"$(0.5*ell)"
    tree = replace(tree,  op => "(#H:$(ell)::$(gamma),"*op)
    tree = replace(tree, clo => clobis*"):$(0.5*ell)")
    ## end hyb
    op = "A$(2^(i-1)+1):$(ell)"
    clo = "A$(2^(i-1) + 2^(j-1)):$(ell)"
    for k in 2:j
        op = "("*op
        clo = clo*"):$(ell)"
    end
    clobis = clo[1:(length(clo)-length("$(ell)"))]*"$(0.5*ell)"
    tree = replace(tree, op  => "("*op)
    tree = replace(tree, clo => clobis*")#H:$(0.5*ell)::$(1-gamma)")
    return(tree)
end

"""
    symmetricnet_newick(n::Int, h::Int, γ::Real)

Create a string for a symmetric network with 2^n tips, numbered from 1 to 2^n,
with a symmetric major tree, whose branch lengths are all equal.
2^(n-h) hybrids are added from level h to h-1 "symmetrically".
The network is time-consistent and ultrametric, with a total height of 1.
"""
function symmetricnet_newick(n::Int, h::Int, gamma::Real, i=1::Int)
    (n < h || h < 2) && error("must have n >= h > 1.")
    # length of branch
    ell = 1.0/n
    # Element net
    net = symmetricnet_newick(h, h, h-1, gamma, ell)
    # Iterate
    if n==h return(net) end
    net = chop(net)
    for k in 1:(n-h)
        net = "(" * net * ":$(ell)," * net * ":$(ell))"
    end
    net = net * ";"
    # Rename hybrids
    net = replace(net, "H" => "H$(2^(n-h))")
    for k in (2^(n-h)-1):-1:1
        net = replace(net, "H$(k+1)" => "H$(k)", count=2*k)
    end
    # Rename tips
    for k in (2^h):-1:1
        net = replace(net, "A$(k):" => "A$(i-1+2^n):")
    end
    for k in (2^n-1):-1:1
        net = replace(net, "A$(i-1+k+1):" => "A$(i-1+k):", count=k)
    end
    return(net)
end


text:
struct MixSkewT{T<:AbstractVector{<:Real}, A<:Dirichlet,
                B<:Union{Normal, <:OrderedNormalMeanPrior},
                C<:Real, D<:LogNormal, E<:Normal, F<:Gamma} <: MCMC.Model
  y::T
  K::Int
  eta::A
  mu::B
  a_omega::C
  nu::D  # If want GMM, set `nu = LogNormal(m, 0.0)`, where `m > 10000`.
  psi::E  # If want no skew, set `psi = Normal(0, 0)`.
  tau::F
end

has_ordered_mu(m::MixSkewT) = m.mu isa OrderedNormalMeanPrior

function MixSkewT(y, K; eta=Dirichlet(K, 1/K), mu=Normal(), a_omega=2.5,
                  nu=LogNormal(2, .5), psi=Normal(0, 3), tau=Gamma(0.5, 1),
                  skew::Bool=true, tdist::Bool=true)
  skew || (psi = Normal(0, 0))
  tdist || (nu = LogNormal(20000, 0))
  return MixSkewT(y, K, eta, mu, a_omega, nu, psi, tau)
end

uses_skew(m::MixSkewT) = (m.psi !== Normal(0, 0))
uses_tdist(m::MixSkewT) = !(m.nu.μ > 10000 && m.nu.σ == 0)

function MCMC.make_init_state(m::MixSkewT)
  N = length(m.y)
  eta = rand(m.eta)
  lambda = rand(Categorical(eta), N)
  tau = rand(m.tau)
  mu = if has_ordered_mu(m)
    # cumsum(rand(m.mu))
    cumsum(mean.(m.mu.priors))
  else
    rand(m.mu, m.K)
  end
  omega = rand(InverseGamma(m.a_omega, tau), m.K)
  nu = rand(m.nu, m.K)
  psi = rand(m.psi, m.K)

  if uses_tdist(m)
    v = rand.(Gamma.(nu[lambda] / 2, 2 ./ nu[lambda]))
  else
    v = one.(m.y)
  end

  if uses_skew(m)
    zeta = rand.(truncated.(Normal.(0, 1 ./ sqrt.(v)), 0, Inf))
  else
    zeta = zero.(m.y)
  end

  sigma = @.scalefromaltskewt(sqrt(omega), psi)
  phi = @.skewfromaltskewt(sqrt(omega), psi)

  return (eta=eta, lambda=lambda, tau=tau, mu=mu, omega=omega, nu=nu, psi=psi,
          v=v, zeta=zeta, sigma=sigma, phi=phi)
end

"""
    make_sampler(m::MixSkewT; init, skew::Bool=true, tdist::Bool=true)

Note that when `skew=false` and `tdist=false`, then we have a Gaussian mixture
model.
"""
function make_sampler(m::MixSkewT; init=nothing)
  init === nothing && (init = MCMC.make_init_state(m))

  if uses_tdist(m)
    cond_nu = RWM(:nu, logprob_nu, mvarwm(init.nu), bijector=Bijectors.Log{1}())
    _update_v = update_v
  else
    cond_nu = Conditional(:nu, (m, s) -> fill(20000.0, m.K))
    _update_v = (m, s) -> one.(m.y)
  end

  if uses_skew(m)
    _update_zeta = update_zeta
    _update_psi = update_psi
  else
    _update_zeta = (m, s) -> zero.(m.y)
    _update_psi = (m, s) -> zeros(m.K)
  end

  if has_ordered_mu(m)
    _update_mu = (model, s) -> let
      s2 = s.omega[s.lambda] ./ s.v
      g = model.y - s.psi[s.lambda] .* s.zeta
      update_mu(model.mu, s.mu, g, s2, s.lambda)
    end
  else
    _update_mu = update_mu
  end

  return Gibbs(m,
               Conditional(:lambda, update_lambda),
               Conditional(:v, _update_v),
               Conditional(:zeta, _update_zeta),
               Conditional(:eta, update_eta),
               Conditional(:mu, _update_mu),
               Conditional(:tau, update_tau),
               Conditional(:omega, update_omega),
               Conditional(:psi, _update_psi),
               cond_nu,
               Conditional((:sigma, :phi), update_sigma_phi))
end

function print_model_info(m::MixSkewT)
  println("Model info:")
  for key in fieldnames(MixSkewT)
    if key == :y
      println("N: ", length(m.y))
    else
      println("$(key): ", getfield(m, key))
    end
  end
end

function make_callback(model::MCMC.Model; nburn, nsamps, thin)
  return function callback(chain, state, sample, i, metrics, iterator)
    flush(stdout)
    if i == 1
      metrics[:loglike] = Float64[]
    elseif i > nburn && mod(i, thin) == 0
      ll = loglike(model, state)
      append!(metrics[:loglike], ll)
      if iterator isa MCMC.ProgressBar
        MCMC.ProgressBars.set_postfix(iterator, loglike=round(ll, digits=3))
      end
    end
  end
end

"""
Finding a good starting value is important when using OrderedNormalMeanPrior.
"""
function find_good_seed(model::MCMC.Model; seeds=1:10, nsamps=1, nburn=200,
                        thin=1, kwargs...)
  # Preallocate results in Dict{seed => mean(loglike)}.
  results = Dict((seed, 0.0) for seed in seeds)
  callback = make_callback(model, nsamps=nsamps, nburn=nburn, thin=thin)

  for seed in MCMC.ProgressBar(seeds)
    Random.seed!(seed)
    spl = make_sampler(model)
    chain, metrics = mcmc(spl, nsamps, nburn=nburn, thin=thin,
                          callback=callback, progress=false, kwargs...)
    results[seed] = mean(metrics[:loglike])
  end

  best_seed = argmax(results)
  return best_seed, results[best_seed]
end


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("Foundations of machine learning")

# ## Learning with a single neuron
#
# In this notebook, we'll build a neuron that classifies an image as an apple or as a banana using multiple features from the image. We'll **train** our neuron using data from many images that have already been correctly categorised; our neuron will thereby **learn** what parameters to use, by minimizing a loss function using gradient descent.
#
# We'll do this with almost the simplest possible neuron, namely one that takes just two inputs:

include(datapath("scripts/draw_neural_net.jl"))
number_inputs, number_neurons = 2, 1
draw_network([number_inputs, number_neurons])

# ### Roadmap
#
# To use a neuron with two inputs, we will use just two of the data (numbers) for each image, say columns 3 and 4, the mean amount of red and the mean amount of green; each data point will then be a 2-dimensional vector, and the data points lie on a two-dimensional plane. We will have many input data, labelled with an index $i$. We will denote the $i$th data point as  $\mathbf{x}^{(i)}$.
#
# The goal is that our neuron will take a single point on the two-dimensional plane as its input, and should return a single output that **classifies** it as an apple ($0$) or a banana ($1$).
#
# To do so, it must "**learn**" the correct values of its parameters $\mathbf{w}$ and $b$. For this learning to take place, we'll need **labels** for each of the input data, which identify them as either an apple (0) or as a banana (1).
#
# These labels will, in turn, allow us to create a loss function, which will allow our algorithm to learn to determine if a given choice of parameters does a good or a poor job of classifying our fruit images.
#
# **So what do we have left to do?**
#
# The above might sound complex, but luckily we can break it down into a series of actionable steps:
#
# 1. Clean our input data (amounts of red and green) to get it into a usable format;
# 2. Create a sequence of labels that we can use to identify correct and incorrect classifications;
# 3. Define a loss function that contains parameters;
# 4. Implement an algorithm to pick parameters for our neuron by minimizing the loss function with respect to the parameters;
# 5. Use all of the above to train our neuron how to classify images as apples or bananas!

#-

# #### Note:
#
# Note that *in general we cannot expect that a single neuron will be adequate for classification.*
#
# If a single neuron struggles to classify our images, we may need to use a more complicated neural network structure (which corresponds to using a more complicated function).

#-

# To do this, we need to work with and *clean* some real data. Let's get started!

#-

# ## Loading in some data

#-

# Let's load in some real data! We'll use data that we have prepared from photos of apples and bananas; it turns out to be stored on disk in data files as "tab-separated values". We can read this data in with the `CSV.jl` package, as follows.

#nb ;head $(datapath("data/Apple_Golden_1.dat"))
#jl run(`head $(datapath("data/Apple_Golden_1.dat"))`)

#-

#nb ;head $(datapath("data/bananas.dat"))
#jl run(`head $(datapath("data/bananas.dat"))`)

#-

using CSV, DataFrames

apples = DataFrame(CSV.File(datapath("data/Apple_Golden_1.dat"), delim='\t', allowmissing=:none, normalizenames=true))
bananas = DataFrame(CSV.File(datapath("data/bananas.dat"), delim='\t', allowmissing=:none, normalizenames=true));

# CSV.jl defaults to creating a `DataFrame` from DataFrames.jl to store the information:

apples

#-

typeof(apples)

#-

names(apples)

# Each of the two data sets is stored in a `DataFrame` (from the `DataFrames.jl` package).

#-

# ## Cleaning the data

#-

# Usually it will be necessary to "clean" the data in some way, i.e. pre-process it, before it can be used for whichever task you are interested in.

#-

# Our next *meta*-exercise will be to collect all the data from columns 3 and 4 into a *single* Julia vector `x` (of which each entry is itself a vector), and the labels into a single vector `y`. Let's do this in a series of steps!

#-

# #### Exercise 1
#
# First, let's practice grabbing a single column of a `DataFrame`. To grab a column, you can index into the `DataFrame` with the name of the column you want, passed as a symbol. In Julia, symbols are names preceded by a `:`. For example, we could grab the "height" column from `apples` by indexing into `apples` with the symbol, `:height`:
#
# ```julia
# apples[:height]
# ```
#
# Index into `apples` to get the "red" column. What is the type of the object returned? How many entries does it have?
#
# A) Array, 5

# C) Array, 64

# E) Array, 492

# F) Float64, 492

#-

# #### Solution
#
# To find the answer, grab the `red` column via

apples[:red]

# #### Exercise 2
#
# We can grab an individual entry of a `DataFrame` by specifying the row index of the entry and the column symbol. For example, to access the height of the 4th image of an apple, we would execute
#
# ```julia
# apples[4, :height]
# ```
#
# How much red is there in the 63rd image of a banana?

#-

# #### Solution
#
# We determine the amount of red in the 63rd banana via

bananas[63, :red]

# #### Exercise 3
#
# We want to reorganize data from the 3rd and 4th columns of `apples` and `bananas` to put that data in a single array. Let's start by organizing the data from the 3rd and 4th columns of `apples` into a single array, `x_apples`. Create `x_apples` such that there is a single element in `x_apples` for each image in `apples`. The $i^\mathrm{th}$ element in `x_apples` should be a `Vector`, i.e. a 1D `Array`, with two elements - the amount of red and the amount of blue in the $i^\mathrm{th}$ image from `apples`.
#
# Similarly create the `Array` `x_bananas` using data from `bananas`.

#-

# #### Solution
#
# Create `x_apples` and `x_bananas` via

x_apples = [[row.red, row.green] for row in eachrow(apples)]
x_bananas = [[row.red, row.green] for row in eachrow(bananas)]

# **Tests**

@assert isapprox(x_apples[404][1],  0.6774070822358786)
@assert isapprox(x_apples[404][2], 0.5832823897939495)
@assert isapprox(x_bananas[404][1], 0.5384643155305847)
@assert isapprox(x_bananas[404][2], 0.4738283806042557)

# #### Exercise 4
#
# Next we want to combine the elements of `x_apples` and `x_bananas` into a single array, `xs`. `xs` should contain, first, all the elements of `x_apples`, and then all the elements of `x_bananas`. Use the `vcat` function to create `xs`.

#-

# #### Solution
#
# The following will produce the desired `xs`:

xs = vcat(x_apples, x_bananas)

# **Tests**

isapprox(xs[404][1], 0.6774070822358786)
isapprox(xs[404][2], 0.5832823897939495)
isapprox(xs[808][1], 0.5771131704577359)
isapprox(xs[808][2], 0.4821497596511228)

# #### Exercise 5
#
# If you've gotten this far, our data is in the format we want for learning. Now we need labels! We want to store a label (either a `0` or a `1` for every apple or banana image in our data set in an array called `ys`. Recall that "0" refers to an apple and "1" refers to a banana.
#
# Create an array `ys` where the $i^\mathrm{th}$ element of `ys` is a `0` if the $i^\mathrm{th}$ element of `xs` is an apple, and where the $i^\mathrm{th}$ element of `ys` is a `1` if the $i^\mathrm{th}$ element of `xs` is a banana.

#-

# #### Solution
#
# You can generate `ys` via
# ```julia
# ys = [0 for i in 1:length(x_apples)]
# ys = vcat(ys, [1 for i in 1:length(x_bananas)])
# ```
#
# **Tests**
#
# ```julia
# ys[1] == 0
# ys[492] == 0
# ys[493] == 1
# ys[982] == 1
# ```

ys = vcat(fill(0, size(x_apples)), fill(1, size(x_bananas)))

# #### Exercise 6
#
# Add data points for all apple and banana images in our data sets to a plot using `scatter`. Plot data points for apples in one color and use a different color for banana data points.
#
# Hint: You may want to use the `first` and `last` functions.

#-

# #### Solution
#
# The following code will produce the desired plot:
#
# ```julia
# scatter(first.(x_apples), last.(x_apples), m=:cross, label="apples")
# scatter!(first.(x_bananas), last.(x_bananas), label="bananas")
# ```

using Plots
scatter(first.(x_apples), last.(x_apples), m=:cross, label="apples")
scatter!(first.(x_bananas), last.(x_bananas), label="bananas")

# ## "Learning" by hand

#-

# Intuitively, looking at the plot of the data, we see that it should be "easy" to find a function that separates the data into bananas on one side and apples on the other: we just need to draw a straight line that divides the two clouds of data. We can do this "by hand", as follows.

#-

# In the code below, the neuron will learn a function of the form $\sigma(\mathbf{w} \cdot \mathbf{x} + b)$. Since $\sigma$ looks like a smooth version of a step function, we can think of $\sigma$ classifying based on whether the value of its output argument is less than `0.5` or greater than `0.5`.
#
# **Game**: Use the interactive visualization to find suitable values of $\mathbf{w}$ and $b$ such that the hyperplane $\sigma(w_1 x_1 + w_2 x_2 + b) = 0.5$ divides the data. This is the same as the hyperplane for which $w_1 x_1 + w_2 x_2 + b = 0$ ! (Note that there are many such values!)

#-

# We can solve for $x_2$ via
#
# $$x_2 = -(w_1 x_1 + b) / w_2,$$
#
# and use this to draw the corresponding hyperplane.

#-

w1 = 3.0 # Try manipulating w1 to have values between -2 and 6
w2 = 3.0 # Try manipulating w2 to have values between -2 and 6
b  = -2.5 # Try manipulating b  to have values between -5 and 0
scatter(first.(x_apples), last.(x_apples), m=:cross, label="apples", xlim=(.3, .9), ylim=(.2,.8))
scatter!(first.(x_bananas), last.(x_bananas), label="bananas")
plot!(x -> -(w1*x + b) / w2)
## w1 = 2.22
## w2 = 1.05
## b  = -2.0

# ## How can the neuron *learn* to classify the data?

#-

# We are now ready for our first experience of **machine learning**: we will let the neuron learn automatically by processing data and tuning model parameters accordingly (the process we call "learning")!

#-

# For given values of the parameters $w_1$, $w_2$ and $b$, the function $f_{\mathbf{w}, b}$ maps a vector of length $2$ to a number between $0$ and $1$ (due to the definition of $\sigma$). Now we want to have a neuron *learn* suitable values of these parameters.
#
# We want to discover (learn!) the parameters such that $f$ models the relationship between the data we explored above about the fruit images, and outputs a classification of the fruit: $0$ if it corresponds to an apple, and $1$ if it is a banana.
#
# So the neuron's input will be a vector of 2 pieces of information about an image; let's call the data about the $i$th image $\mathbf{x}^{(i)}$.
# We also are given the label that says which type of fruit it is: $0$ for an apple, and $1$ for a banana; let's call this *desired* output number $y^{(i)}$. When we feed in the $i$th data, $\mathbf{x}^{(i)}$, we want the neuron to give an output that is *as close as possible* to the desired output $y^{(i)}$; i.e. it should **minimize** the mean-squared distance
#
# $$L_i = [f_{\mathbf{w}, b}(\mathbf{x}^{(i)}) - y^{(i)} ]^2.$$

#-

# However, now we see a key difference from what we did previously: the neuron should vary its parameters in such a way that it manages to minimize this distance for *all* of the input data, simultaneously!
#
# How can we express this mathematically? We once again define a loss function, $L(\mathbf{w}, b)$, which tells us "how wrong" we are when the parameters take on the given values, and then **minimize** this loss function with respect to all of its parameters.
#
# One way to take account of all the data at once is to use the "mean-squared error" loss function, which is the mean (squared) over all the differences between the output of the network, $f_{\mathbf{w}, b}(\mathbf{x}^{(i)})$ on the $i$th data, and the desired output $y^{(i)}$:
#
# $$L_\mathrm{total}(\mathbf{w}, b) = \frac{1}{N} \sum_i L_i = \frac{1}{N} \sum_i [f_{\mathbf{w}, b}(\mathbf{x}^{(i)}) - y^{(i)} ]^2,$$
#
# where $N$ is the total number of data in the training set.
#
# Why do we choose this particular loss function? Because the minimum possible value of this loss function is $0$ (since it is a sum of squares), and this is reached only when the neural network perfectly predicts the output. If we can find a way to minimize this loss function, we will get as close as possible to this perfect prediction. (In general, though, we won't be able to get an exact prediction.)

#-

# ## Minimizing the loss function: *stochastic* gradient descent

#-

# We already know how to minimise loss functions on a computer: we just calculate the gradient, and do gradient descent! But here we hit a problem: the function $L_\mathrm{total}$ usually has a *lot* of terms, and so calculating the gradient of that function will be very time-consuming.
#
# Instead, we will use a variant of gradient descent, called *stochastic* gradient descent. Here, the idea is that we will not use the complete loss function; instead, at each step we will choose a random data point, number $i$, and do a step of gradient descent for the partial loss function $L_i$ *corresponding to only that data point*.

#-

# **Exercise 7:**
#
# Write functions for the partial loss function `L(w, b, x, y)`.
#
# To do this, recall
#
# $$
# \mathbf{x} = \begin{pmatrix} x_1 \\ x_2 \end{pmatrix};
# \qquad
# \mathbf{w} = \begin{pmatrix} w_1 \\ w_2 \end{pmatrix};
# \qquad
# f_{\mathbf{w}, b}(\mathbf{x}) = \sigma(\mathbf{w} \cdot \mathbf{x} + b),$$
#
# and declare `f(x, w, b)` as in notebook 8.

#-

# #### Solution

## Setup:
σ(x) = 1 / (1 + exp(-x))
f(x, w, b) = σ(w'x + b)
## The loss function:
L(w, b, x, y) = (f(x, w, b) - y)^2

# **Tests**

@assert isapprox(L([-1.5, 1.5], .5, [.3, .4], 1), 0.11764182271544733)
@assert isapprox(L([-1, 1], 0, [-.4, -.3], 0), 0.27560314728604807)

# #### Exercise 8
#
# Write a function for the gradient of `L`, i.e. `∇L(w, b, x, y)`, with respect to the parameters $(w_1, w_2, b)$, using finite differences. $∇L$ will be a vector with one component per parameter:
#
# $$∇L = \left( \frac{\partial L}{\partial w_1}, \frac{\partial L}{\partial w_2}, \frac{\partial L}{\partial b} \right).$$

#-

# #### Solution

## A separate function for each partial:
partial_w1(w, b, x, y, h = 0.001) = (L(w .+ [h, 0], b, x, y) - L(w, b, x, y))/h
partial_w2(w, b, x, y, h = 0.001) = (L(w .+ [0, h], b, x, y) - L(w, b, x, y))/h
partial_b(w, b, x, y, h = 0.001) = (L(w, b + h, x, y) - L(w, b, x, y))/h

## The combined gradient:
∇L(w, b, x, y) = [partial_w1(w, b, x, y), partial_w2(w, b, x, y), partial_b(w, b, x, y)]

# **Tests**

@assert isapprox(∇L([-1.5, 1.5], .5, [.3, .4], 1)[1], -0.04636839045066288)
@assert isapprox(∇L([-1.5, 1.5], .5, [.3, .4], 1)[2], -0.061821518678797927)
@assert isapprox(∇L([-1.5, 1.5], .5, [.3, .4], 1)[3], -0.15450877075851022)
@assert isapprox(∇L([-1, 1], 0, [-.4, -.3], 0)[1], -0.10472487950158316)
@assert isapprox(∇L([-1, 1], 0, [-.4, -.3], 0)[2], -0.07854532958612159)
@assert isapprox(∇L([-1, 1], 0, [-.4, -.3], 0)[3], 0.2618900851394379)

# #### Exercise 9
#
# Implement stochastic gradient descent in the function `stochastic_gradient_descent(L, w, b, xs, ys, N = 1000)`. Use this function to minimize the function $L_\mathrm{total}$.
#
# The algorithm: For each of `N` steps, randomly select an index $i$ into the vector `xs` storing your image data. Calculate the gradient of the cost function, $L_i$, for this image and update each of the parameters, $p_j$, of $L_i$ according to
#
# $$p_j = p_j - 0.01 * ∇L_j$$
#
# (Here, $j$ signifies the $j^{th}$ parameter of $L$ and similarly the $j^{th}$ component of $∇L$.)
#
# `stochastic_gradient_descent` should return the updated values for vector $\mathbf{w}$ and scalar $b$.
#
# Optional: Keep track of the value of $L_\mathrm{total}$ over time if you want to visualize the learning process.

#-

# #### Solution

function stochastic_gradient_descent(∇L, w, b, xs, ys, N=1000)

    η = 0.01

    for step in 1:N

        i = rand(1:length(xs))  # choose a data point

        x = xs[i]
        y = ys[i]

        grad = ∇L(w, b, x, y)
        w = w .- η .* grad[1:2]
        b = b - η * grad[3]
    end

    return w, b

end

# **Tests**

w, b = ([0.736019, 0.518619], 0.03408673204690893)
w, b = stochastic_gradient_descent(∇L, w, b, xs, ys, 1e6)
@assert -18.8 < w[1] < -18.4 # roughly -18.6
@assert -19.4 < w[2] < -19.0 # roughly -19.2
@assert 21 < b < 23 # roughly 22

# #### Exercise 10
#
# Use the values of `w` and `b` from the last exercise to see how `f` is classifying a couple of the images in out dataset.
#
# In particular, calculate `f` using the 1st and the 900th image in `xs`. For which image is the output of `f` closer to the value of its label?
#
# A) The output of `f` for the 1st image in `xs` is closer to its label
#
# B) The output of `f` for the 900th image in `xs` is closer to its label.

#-

# #### Solution
#
# A)
#
# ```julia
# f(xs[900], w, b) # yields ~0.9565428318826446
# ys[900] - f(xs[900], w, b) # yields ~.041 (4% difference)
# f(xs[1], w, b) # yields ~0.027378537417953305
# ys[1] - f(xs[1], w, b) # yields ~ -0.030 (3% difference)
# ```

@show f(xs[1], w, b) # yields ~0.027378537417953305
@show ys[1] - f(xs[1], w, b) # yields ~ -0.030 (3% difference)

@show f(xs[900], w, b) # yields ~0.9565428318826446
@show ys[900] - f(xs[900], w, b) # yields ~.041 (4% difference);

# We see that *with sufficient training*, the single neuron is approximately able to learn the function for most of the data, but not all:

#-

# #### Exercise 11
#
# Use the `maximum` function to determine the maximum squared distance of the prediction from the true value. (For each image, this formula is $y_i - f_{w, b}(x_i)$.)

#-

# #### Solution

@show maximum((ys .- f.(xs, [w], [b])).^2)

# #### Exercise 12
#
# Use `w` and `b` from stochastic gradient descent to draw the function that the network has learned, as before, as the hyperplane $w_1 x + w_2 y + b = 0$. Overlay this with the data.
#
#
# Does this hyperplane correctly separate the data? (I.e., is the data for all apples on one side of the line, and is the data for all bananas on the other side of the line?)
#
# A) Yes
# B) No

#-

# #### Solution
#
# A) Yes

using Plots; gr()

plot(x ->(-w[1]*x - b) / w[2])
scatter!(first.(x_apples), last.(x_apples), m=:cross)
scatter!(first.(x_bananas), last.(x_bananas))
xlims!(0.45, 0.75)
ylims!(0.35, 0.75)

# Since the neuron has two input variables and one output variable, we can also use a "heatmap" to visualize the value of $f_{\mathbf{w}, b}(x_1, x_2)$ as a function of the input data pairs $(x_1, x_2)$:

heatmap(0:0.01:1, 0:0.01:1, (x,y) -> f([x,y], w, b))

scatter!(first.(x_apples), last.(x_apples), m=:cross)
scatter!(first.(x_bananas), last.(x_bananas))

xlims!(0.45, 0.75)
ylims!(0.35, 0.75)

# We see that the single neuron has indeed *learned* to separate the data using a hyperplane.

max_stars_repo_path: Courses/Foundations of machine learning/1100.Intro-to-neurons-quiz.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
#=
demo_HDP_MultinomialDirichlet.jl

A demo for Hierarchical Dirichlet Process mixture models with MultinomialDirichlet
Bayesian components. This demo uses auxiliary variable method for for inference.

28/07/2015
Adham Beyki, odinay@gmail.com
=#

using BNP
srand(123)

true_gg = 1.0
true_aa = 0.5
n_groups = 100
n_group_j   = 100 * ones(Int, n_groups)
KK_truncation = 25

beta_tilde = rand(Distributions.Beta(1, true_gg), KK_truncation)
beta_tilde[KK_truncation] = 1.0
my_beta = BNP.stick_breaking(beta_tilde)

pi_tilde = zeros(Float64, n_groups, KK_truncation)
my_pi = zeros(Float64, n_groups, KK_truncation)
for jj = 1:n_groups
    for kk = 1:KK_truncation-1
        gamma1 = true_aa * my_beta[kk]
        gamma2 = true_aa * (1 - sum(my_beta[1:kk]))
        pi_tilde[jj, kk] = rand(Distributions.Beta(gamma1, gamma2))
    end
    pi_tilde[jj, KK_truncation] = 1.0
    my_pi[jj, :] = BNP.stick_breaking(pi_tilde[jj, :][:])
end

true_zz = Array(Vector{Int}, n_groups)
for jj = 1:n_groups
    true_zz[jj] = zeros(Int, n_group_j[jj])
    for ii = 1:n_group_j[jj]
        kk = sample(my_pi[jj, :][:])
        true_zz[jj][ii] = kk
    end
end

zz_flat = Int[]
for jj = 1:n_groups
	append!(zz_flat, true_zz[jj])
end

zz_unq = sort(unique(zz_flat))
true_KK = length(zz_unq)
for jj = 1:n_groups
	for kk = 1:true_KK
		idx = find(x -> x==zz_unq[kk], true_zz[jj])
		true_zz[jj][idx] = kk
	end
end

n_groups   = 50
n_group_j  = 100 * ones(Int, n_groups)
vocab_size = 25

# constructing the topics
true_topics = BNP.gen_bars(true_KK, vocab_size, 0)

# constructing the observations and labels
xx = Array(Vector{Int}, n_groups)
true_nn = zeros(Int, n_groups, true_KK)

for jj = 1:n_groups
    xx[jj] = zeros(Int, n_group_j[jj])
    for ii = 1:n_group_j[jj]
        kk = true_zz[jj][ii]
        xx[jj][ii] = sample(true_topics[kk, :][:])
        true_nn[jj, kk] += 1
    end
end

## ------- inference -------- ##
# constructing the Bayesian component
dd = vocab_size
aa = 0.5*dd
q0 = MultinomialDirichlet(dd, aa)

# constructing the  model
dhdp_KK_init = 1
KK_truncation = 15
dhdp_gg = 1.0
dhdp_g1 = 0.1
dhdp_g2 = 0.1
dhdp_aa = 0.5
dhdp_a1 = 0.1
dhdp_a2 = 0.1
dhdp_aw = 2.0
dhdp_bw = 3.0
dhdp = dHDP(q0, dhdp_KK_init, dhdp_gg, dhdp_g1, dhdp_g2, dhdp_aa, dhdp_a1, dhdp_a2, dhdp_aw, dhdp_bw)

# sampling
zz = Array(Vector{Int}, n_groups)
for jj = 1:n_groups
    zz[jj] = ones(Int, n_group_j[jj])
    zz[jj] = rand(1:dhdp.KK, n_group_j[jj])
end

n_burnins   = 500
n_lags      = 1
n_samples   = 500
sample_hyperparam = false
n_internals = 10
store_every = 1000
results_path = ""
filename    = "demo_dHDP_MultinomialDirichlet"

KK_list, KK_dict, w_tilde = BNP.truncated_gibbs_sampler(dhdp, KK_truncation, xx, zz,
    n_burnins, n_lags, n_samples, sample_hyperparam, n_internals, store_every, results_path, filename)

KK_hist = hist(KK_list, 0.5:maximum(KK_list)+0.5)[2]
candidate_KK = indmax(KK_hist)
posterior_components, nn = posterior(dhdp, xx, KK_dict[candidate_KK], candidate_KK)


text:
using Oceananigans.Grids: AbstractGrid

const ω̂₁ = 5/18
const ω̂ₙ = 5/18
const ε₂ = 1e-20

# Here in the future we can easily add UpwindBiasedFifthOrder
const BoundPreservingScheme = PositiveWENO

@inline function div_Uc(i, j, k, grid, advection::BoundPreservingScheme, U, c)

    div_x = bounded_tracer_flux_divergence_x(i, j, k, grid, advection, U.u, c)
    div_y = bounded_tracer_flux_divergence_y(i, j, k, grid, advection, U.v, c)
    div_z = bounded_tracer_flux_divergence_z(i, j, k, grid, advection, U.w, c)

    return 1/Vᶜᶜᶜ(i, j, k, grid) * (div_x + div_y + div_z)
end

# Support for Flat directions
@inline bounded_tracer_flux_divergence_x(i, j, k, ::AbstractGrid{FT, Flat, TY, TZ}, args...) where {FT, TY, TZ} = zero(FT)
@inline bounded_tracer_flux_divergence_y(i, j, k, ::AbstractGrid{FT, TX, Flat, TZ}, args...) where {FT, TX, TZ} = zero(FT)
@inline bounded_tracer_flux_divergence_z(i, j, k, ::AbstractGrid{FT, TX, TY, Flat}, args...) where {FT, TX, TY} = zero(FT)

@inline function bounded_tracer_flux_divergence_x(i, j, k, grid, advection::BoundPreservingScheme, u, c)

    lower_limit = advection.bounds[1]
    upper_limit = advection.bounds[2]

    cᵢⱼ  = c[i, j, k]

    c₊ᴸ =  _left_biased_interpolate_xᶠᵃᵃ(i+1, j, k, grid, advection, c)
    c₊ᴿ = _right_biased_interpolate_xᶠᵃᵃ(i+1, j, k, grid, advection, c)
    c₋ᴸ =  _left_biased_interpolate_xᶠᵃᵃ(i,   j, k, grid, advection, c)
    c₋ᴿ = _right_biased_interpolate_xᶠᵃᵃ(i,   j, k, grid, advection, c)

    p̃   =  (cᵢⱼ - ω̂₁ * c₋ᴿ - ω̂ₙ * c₊ᴸ) / (1 - 2ω̂₁)
    M   = max(p̃, c₊ᴸ, c₋ᴿ)
    m   = min(p̃, c₊ᴸ, c₋ᴿ)
    θ   = min(abs((upper_limit - cᵢⱼ)/(M - cᵢⱼ + ε₂)), abs((lower_limit - cᵢⱼ)/(m - cᵢⱼ + ε₂)), one(grid))

    c₊ᴸ = θ * (c₊ᴸ - cᵢⱼ) + cᵢⱼ
    c₋ᴿ = θ * (c₋ᴿ - cᵢⱼ) + cᵢⱼ

    return Axᶠᶜᶜ(i+1, j, k, grid) * upwind_biased_product(u[i+1, j, k], c₊ᴸ, c₊ᴿ) -
           Axᶠᶜᶜ(i,   j, k, grid) * upwind_biased_product(u[i,   j, k], c₋ᴸ, c₋ᴿ)
end

@inline function bounded_tracer_flux_divergence_y(i, j, k, grid, advection::BoundPreservingScheme, v, c)

    lower_limit = advection.bounds[1]
    upper_limit = advection.bounds[2]

    cᵢⱼ  = c[i, j, k]

    c₊ᴸ =  _left_biased_interpolate_yᵃᶠᵃ(i, j+1, k, grid, advection, c)
    c₊ᴿ = _right_biased_interpolate_yᵃᶠᵃ(i, j+1, k, grid, advection, c)
    c₋ᴸ =  _left_biased_interpolate_yᵃᶠᵃ(i, j,   k, grid, advection, c)
    c₋ᴿ = _right_biased_interpolate_yᵃᶠᵃ(i, j,   k, grid, advection, c)

    p̃   =  (cᵢⱼ - ω̂₁ * c₋ᴿ - ω̂ₙ * c₊ᴸ) / (1 - 2ω̂₁)
    M   = max(p̃, c₊ᴸ, c₋ᴿ)
    m   = min(p̃, c₊ᴸ, c₋ᴿ)
    θ   = min(abs((upper_limit - cᵢⱼ)/(M - cᵢⱼ + ε₂)), abs((lower_limit - cᵢⱼ)/(m - cᵢⱼ + ε₂)), one(grid))

    c₊ᴸ = θ * (c₊ᴸ - cᵢⱼ) + cᵢⱼ
    c₋ᴿ = θ * (c₋ᴿ - cᵢⱼ) + cᵢⱼ

    return Ayᶜᶠᶜ(i, j+1, k, grid) * upwind_biased_product(v[i, j+1, k], c₊ᴸ, c₊ᴿ) -
           Ayᶜᶠᶜ(i, j,   k, grid) * upwind_biased_product(v[i, j,   k], c₋ᴸ, c₋ᴿ)
end

@inline function bounded_tracer_flux_divergence_z(i, j, k, grid, advection::BoundPreservingScheme, w, c)

    lower_limit = advection.bounds[1]
    upper_limit = advection.bounds[2]

    cᵢⱼ  = c[i, j, k]

    c₊ᴸ =  _left_biased_interpolate_zᵃᵃᶠ(i, j, k+1, grid, advection, c)
    c₊ᴿ = _right_biased_interpolate_zᵃᵃᶠ(i, j, k+1, grid, advection, c)
    c₋ᴸ =  _left_biased_interpolate_zᵃᵃᶠ(i, j, k,   grid, advection, c)
    c₋ᴿ = _right_biased_interpolate_zᵃᵃᶠ(i, j, k,   grid, advection, c)

    p̃   =  (cᵢⱼ - ω̂₁ * c₋ᴿ - ω̂ₙ * c₊ᴸ) / (1 - 2ω̂₁)
    M   = max(p̃, c₊ᴸ, c₋ᴿ)
    m   = min(p̃, c₊ᴸ, c₋ᴿ)
    θ   = min(abs((upper_limit - cᵢⱼ)/(M - cᵢⱼ + ε₂)), abs((lower_limit - cᵢⱼ)/(m - cᵢⱼ + ε₂)), one(grid))

    c₊ᴸ = θ * (c₊ᴸ - cᵢⱼ) + cᵢⱼ
    c₋ᴿ = θ * (c₋ᴿ - cᵢⱼ) + cᵢⱼ

    return Azᶜᶜᶠ(i, j, k+1, grid) * upwind_biased_product(w[i, j, k+1], c₊ᴸ, c₊ᴿ) -
           Azᶜᶜᶠ(i, j, k,   grid) * upwind_biased_product(w[i, j, k],   c₋ᴸ, c₋ᴿ)
end


text:
#options
#=
mutable struct Options
    method
    ulim
    llim
    vlink
    x_tol
    f_tol
    g_tol
    x_abstol::Real
    x_reltol::Real
    f_abstol::Real
    f_reltol::Real
    g_abstol::Real
    g_reltol::Real
    outer_x_tol
    outer_f_tol
    outer_g_tol
    outer_x_abstol::Real
    outer_x_reltol::Real
    outer_f_abstol::Real
    outer_f_reltol::Real
    outer_g_abstol::Real
    outer_g_reltol::Real
    f_calls_limit::Int
    g_calls_limit::Int
    h_calls_limit::Int
    allow_f_increases::Bool
    allow_outer_f_increases::Bool
    successive_f_tol::Int
    iterations::Int
    outer_iterations::Int
    store_trace::Bool
    trace_simplex::Bool
    show_trace::Bool
    extended_trace::Bool
    show_every::Int
    callback
    time_limit
end

function Options(;
method = Optim.Newton(),
ulim = nothing,
llim = nothing,
vlink = nothing,
x_tol = nothing,
f_tol = nothing,
g_tol = nothing,
x_abstol::Real = 0.0,
x_reltol::Real = 0.0,
f_abstol::Real = 0.0,
f_reltol::Real = 0.0,
g_abstol::Real = 1e-8,
g_reltol::Real = 1e-8,
outer_x_tol = 0.0,
outer_f_tol = 0.0,
outer_g_tol = nothing,
outer_x_abstol::Real = 0.0,
outer_x_reltol::Real = 0.0,
outer_f_abstol::Real = 0.0,
outer_f_reltol::Real = 0.0,
outer_g_abstol::Real = 1e-8,
outer_g_reltol::Real = 1e-8,
f_calls_limit::Int = 0,
g_calls_limit::Int = 0,
h_calls_limit::Int = 0,
allow_f_increases::Bool = true,
allow_outer_f_increases::Bool = true,
successive_f_tol::Int = 1,
iterations::Int = 1_000,
outer_iterations::Int = 1000,
store_trace::Bool = false,
trace_simplex::Bool = false,
show_trace::Bool = false,
extended_trace::Bool = false,
show_every::Int = 1,
callback = nothing,
time_limit = NaN)
    Options(
        method,
        ulim,
        llim,
        vlink,
        x_tol,
        f_tol,
        g_tol,
        x_abstol,
        x_reltol,
        f_abstol,
        f_reltol,
        g_abstol,
        g_reltol,
        outer_x_tol,
        outer_f_tol,
        outer_g_tol,
        outer_x_abstol,
        outer_x_reltol,
        outer_f_abstol,
        outer_f_reltol,
        outer_g_abstol,
        outer_g_reltol,
        f_calls_limit,
        g_calls_limit,
        h_calls_limit,
        allow_f_increases,
        allow_outer_f_increases,
        successive_f_tol,
        iterations,
        outer_iterations,
        store_trace,
        trace_simplex,
        show_trace,
        extended_trace,
        show_every,
        callback,
        time_limit)
end
=#


text:
using StatsBase
using Base.Test

X = randn(3, 8)

Z1 = X .- mean(X, 1)
Z2 = X .- mean(X, 2)

w1 = rand(3)
w2 = rand(8)

wv1 = weights(w1)
wv2 = weights(w2)

Z1w = X .- mean(X, wv1, 1)
Z2w = X .- mean(X, wv2, 2)

## reference results

S1 = Z1'Z1
S2 = Z2 * Z2'

Sz1 = X'X
Sz2 = X * X'

S1w = Z1w' * diagm(w1) * Z1w
S2w = Z2w * diagm(w2) * Z2w'

Sz1w = X' * diagm(w1) * X
Sz2w = X * diagm(w2) * X'

## scattermat

if VERSION < v"0.5.0-dev+679"
    @test_approx_eq scattermat(X) S1
    @test_approx_eq scattermat(X; vardim=2) S2

    @test_approx_eq scattermat(X; mean=0) Sz1
    @test_approx_eq scattermat(X; mean=0, vardim=2) Sz2

    @test_approx_eq scattermat(X; mean=mean(X,1)) S1
    @test_approx_eq scattermat(X; mean=mean(X,2), vardim=2) S2

    @test_approx_eq scattermat(X; mean=zeros(1,8)) Sz1
    @test_approx_eq scattermat(X; mean=zeros(3), vardim=2) Sz2

    ## weighted scatter mat

    @test_approx_eq scattermat(X, wv1) S1w
    @test_approx_eq scattermat(X, wv2; vardim=2) S2w

    @test_approx_eq scattermat(X, wv1; mean=0) Sz1w
    @test_approx_eq scattermat(X, wv2; mean=0, vardim=2) Sz2w

    @test_approx_eq scattermat(X, wv1; mean=mean(X, wv1, 1)) S1w
    @test_approx_eq scattermat(X, wv2; mean=mean(X, wv2, 2), vardim=2) S2w

    @test_approx_eq scattermat(X, wv1; mean=zeros(1,8)) Sz1w
    @test_approx_eq scattermat(X, wv2; mean=zeros(3), vardim=2) Sz2w
else
    @test_approx_eq scattermat(X) S1
    @test_approx_eq scattermat(X, 2) S2

    @test_approx_eq StatsBase.scattermatm(X, 0) Sz1
    @test_approx_eq StatsBase.scattermatm(X, 0, 2) Sz2

    @test_approx_eq StatsBase.scattermatm(X, mean(X,1)) S1
    @test_approx_eq StatsBase.scattermatm(X, mean(X,2), 2) S2

    @test_approx_eq StatsBase.scattermatm(X, zeros(1,8)) Sz1
    @test_approx_eq StatsBase.scattermatm(X, zeros(3), 2) Sz2

    ## weighted scatter mat

    @test_approx_eq scattermat(X, wv1) S1w
    @test_approx_eq scattermat(X, wv2, 2) S2w

    @test_approx_eq StatsBase.scattermatm(X, 0, wv1) Sz1w
    @test_approx_eq StatsBase.scattermatm(X, 0, wv2, 2) Sz2w

    @test_approx_eq StatsBase.scattermatm(X, mean(X, wv1, 1), wv1) S1w
    @test_approx_eq StatsBase.scattermatm(X, mean(X, wv2, 2), wv2, 2) S2w

    @test_approx_eq StatsBase.scattermatm(X, zeros(1,8), wv1) Sz1w
    @test_approx_eq StatsBase.scattermatm(X, zeros(3), wv2, 2) Sz2w
end

# weighted covariance

if VERSION < v"0.5.0-dev+679"
    @test_approx_eq cov(X, wv1) S1w ./ sum(wv1)
    @test_approx_eq cov(X, wv2; vardim=2) S2w ./ sum(wv2)

    @test_approx_eq cov(X, wv1; mean=0) Sz1w ./ sum(wv1)
    @test_approx_eq cov(X, wv2; mean=0, vardim=2) Sz2w ./ sum(wv2)

    @test_approx_eq cov(X, wv1; mean=mean(X, wv1, 1)) S1w ./ sum(wv1)
    @test_approx_eq cov(X, wv2; mean=mean(X, wv2, 2), vardim=2) S2w ./ sum(wv2)

    @test_approx_eq cov(X, wv1; mean=zeros(1,8)) Sz1w ./ sum(wv1)
    @test_approx_eq cov(X, wv2; mean=zeros(3), vardim=2) Sz2w ./ sum(wv2)
else
    @test_approx_eq cov(X, wv1) S1w ./ sum(wv1)
    @test_approx_eq cov(X, wv2, 2) S2w ./ sum(wv2)

    @test_approx_eq Base.covm(X, 0, wv1) Sz1w ./ sum(wv1)
    @test_approx_eq Base.covm(X, 0, wv2, 2) Sz2w ./ sum(wv2)

    @test_approx_eq Base.covm(X, mean(X, wv1, 1), wv1) S1w ./ sum(wv1)
    @test_approx_eq Base.covm(X, mean(X, wv2, 2), wv2, 2) S2w ./ sum(wv2)

    @test_approx_eq Base.covm(X, zeros(1,8), wv1) Sz1w ./ sum(wv1)
    @test_approx_eq Base.covm(X, zeros(3), wv2, 2) Sz2w ./ sum(wv2)
end

# mean_and_cov
if VERSION < v"0.5.0-dev+679"
    (m, C) = mean_and_cov(X; vardim=1)
    @test m == mean(X, 1)
    @test C == cov(X; vardim=1)

    (m, C) = mean_and_cov(X; vardim=2)
    @test m == mean(X, 2)
    @test C == cov(X; vardim=2)

    (m, C) = mean_and_cov(X, wv1; vardim=1)
    @test m == mean(X, wv1, 1)
    @test C == cov(X, wv1; vardim=1)

    (m, C) = mean_and_cov(X, wv2; vardim=2)
    @test m == mean(X, wv2, 2)
    @test C == cov(X, wv2; vardim=2)
else
    (m, C) = mean_and_cov(X, 1)
    @test m == mean(X, 1)
    @test C == cov(X, 1)

    (m, C) = mean_and_cov(X, 2)
    @test m == mean(X, 2)
    @test C == cov(X, 2)

    (m, C) = mean_and_cov(X, wv1, 1)
    @test m == mean(X, wv1, 1)
    @test C == cov(X, wv1, 1)

    (m, C) = mean_and_cov(X, wv2, 2)
    @test m == mean(X, wv2, 2)
    @test C == cov(X, wv2, 2)
end


text:
@everywhere begin
    import AARMBEMC.AARMBEMCvdW.vdWenergy
end

vdWenergy(ARGS)


text:
# Automatically generated test file
function discontinuous1()

g = [
0.000000000000000,
0.828539361054709,
1.231479419283098,
1.542809041582063,
1.793807989999907,
1.993807989999906,
2.142809041582063,
2.231479419283098,
2.228539361054709,
1.800000000000000,
-6.000000000000000,
2.000000000000000,
10.000000000000000,
18.000000000000000,
26.000000000000000,
-6.000000000000000,
2.000000000000000,
10.000000000000000,
18.000000000000000,
26.000000000000000]

ζ_vec = [873.808, 250.088, 71.5766, 20.4856, 5.86308, 1.67804, 0.480265, 0.137454, 0.0393402, 0.0112594]

I_sols = [
[1, 20],
[1, 18, 20],
[1, 15, 16, 20],
[1, 12, 15, 16, 20],
[1, 10, 11, 15, 16, 20],
[1, 5, 10, 11, 15, 16, 20],
[1, 3, 8, 10, 11, 15, 16, 20],
Int[],
[1, 2, 4, 7, 9, 10, 11, 15, 16, 20]]

f_sols = [
1102.237767513050130,
873.808100414683622,
544.329717092552983,
57.308828172005178,
1.311071990397977,
0.282505163300357,
0.069492548637299,
0.018264420672040,
0.004190223377918]

return g, ζ_vec, I_sols, f_sols

end


text:
module JuliaZH

import Base.Docs: DocStr

const keywords = Dict{Symbol,DocStr}()

"""
源列表
"""
const mirrors = Dict{String,String}()

mirrors["BFSU"] = "https://mirrors.bfsu.edu.cn/julia/static"
mirrors["default"] = "https://pkg.julialang.org"

"""
    set_mirror(mirror_name="BFSU")

设置镜像源为 `mirror_name`。目前可以选择的有：

- `"BFSU"`: TUNA支持的境内源
- `"default"`: Julia官方默认源
"""
function set_mirror(mirror_name = "BFSU")
    ENV["JULIA_PKG_SERVER"] = mirrors[mirror_name]
    return
end

const regex_PKG_SERVER = r"^\s*[^#]*\s*(ENV\[\"JULIA_PKG_SERVER\"\]\s*=\s*)\"([\w\.:\/]*)\""

"""
    generate_startup(mirror_name::String="BFSU")
自动产生将镜像设置到 `mirror_name` 的 `startup.jl` 文件。默认值是 `BFSU
` 的源。
"""
function generate_startup(mirror_name::String = "BFSU")
    default_dot_julia = first(DEPOT_PATH)
    config_path = joinpath(default_dot_julia, "config")
    if !ispath(config_path)
        mkpath(config_path)
    end

    startup_path = joinpath(config_path, "startup.jl")
    if ispath(startup_path)
        startup_lines = readlines(startup_path)
    else
        startup_lines = String[]
    end

    new_upstream = mirrors[mirror_name]
    new_line = "ENV[\"JULIA_PKG_SERVER\"] = \"$(new_upstream)\""

    pkg_matches = map(x->match(regex_PKG_SERVER, x), startup_lines)
    pkg_indices = findall(x->!isnothing(x), pkg_matches)
    if isempty(pkg_indices)
        @info "添加 PkgServer" 服务器地址=new_upstream 配置文件=config_path
        append!(startup_lines, ["", "# 以下这一行由 JuliaCN 自动生成", new_line, ""])
    else
        # only modify the last match
        idx = last(pkg_indices)
        old_upstream = pkg_matches[idx].captures[2]

        is_upstream_unchanged = occursin(new_upstream, old_upstream) || occursin(old_upstream, new_upstream)
        if !is_upstream_unchanged
            @info "更新 PkgServer" 新服务器地址=new_upstream 原服务器地址=old_upstream 配置文件=config_path
            startup_lines[idx] = new_line
        end
    end

    write(startup_path, join(startup_lines, "\n"))
    return
end

function __init__()
    # set to Chinese env by default
    ENV["REPL_LOCALE"] = "zh_CN"
    set_mirror("BFSU")
    if ccall(:jl_generating_output, Cint, ()) == 0
        include(joinpath(@__DIR__, "i18n_patch.jl"))
        include(joinpath(@__DIR__, "basedocs.jl"))
    end
end

end # module


text:
# This file is a part of NeighborhoodApproximationIndex.jl

using InvertedFiles: getcachepositions

"""
    search(idx::KnrIndex, q, res::KnnResult; ksearch=idx.opt.ksearch, ordering=idx.ordering, pools=getpools(D))

Searches nearest neighbors of `q` inside the `index` under the distance function `dist`.
"""
function search(idx::KnrIndex, q, res::KnnResult; pools=getpools(idx), ksearch=idx.opt.ksearch)
    enc = getencodeknnresult(ksearch, pools)
    search(idx.centers, q, enc)
    ifpools = getpools(idx.invfile)
    Q = prepare_posting_lists_for_querying(idx.invfile, enc, ifpools)
    P = getcachepositions(length(Q), ifpools)
    search_(idx, q, enc, Q, P, res, idx.ordering)
end

function search_(idx::KnrIndex, q, _, Q, P_, res::KnnResult, ::DistanceOrdering)
    dist = idx.dist
    cost = umerge(Q, P_) do L, P, _
        @inbounds objID = L[1].I[P[1]]
        @inbounds push!(res, objID, evaluate(dist, q, idx[objID]))
    end

    (res=res, cost=cost)
end

function search_(idx::KnrIndex, q, enc, Q, P_, res::KnnResult, ordering::DistanceOnTopKOrdering)
    enc = reuse!(enc, ordering.top)
    search(idx.invfile, Q, P_, 1) do objID, d
        @inbounds push!(enc, objID, d)
    end

    for objID in idview(enc)
        @inbounds push!(res, objID, evaluate(idx.dist, q, idx[objID]))
    end

    (res=res, cost=length(enc))
end

function search_(idx::KnrIndex, q, _, Q, P_, res::KnnResult, ::InternalDistanceOrdering)
    cost = search(idx.invfile, Q, P_, 1) do objID, d
        @inbounds push!(res, objID, d)
    end

    (res=res, cost=cost)
end

max_stars_repo_name: sadit/NeighborhoodApproximationIndex

max_issues_repo_name: sadit/NeighborhoodApproximationIndex

max_forks_repo_name: sadit/NeighborhoodApproximationIndex


text:
@system Death(Stage) begin
    dieable => true ~ flag
    dead => begin
        #FIXME implement Count trait equivalent
        #self.pheno.plant.count.dropped_leaves >= self.pheno.leaves_initiated
        false
    end ~ flag
    dying(dieable & !dead) ~ flag

    # #FIXME postprocess similar to @produce?
    # def finish(self):
    #     #FIXME record event?
    #     GDD_sum = self.pheno.gdd_recorder.rate
    #     T_grow = self.pheno.growing_temperature
    #     print(f"* Death: time = {self.time}, GDDsum = {GDD_sum}, Growing season T = {T_grow}")
end


text:
# periods.jl

function Base.:+(x::QDate, y::Dates.CompoundPeriod)
    for p in y.periods
        if isa(p, Dates.Week)
            x += Dates.Day(value(p) * 7)
        else
            x += p
        end
    end
    return x
end

function Base.:-(x::QDate, y::Dates.CompoundPeriod)
    for p in y.periods
        if isa(p, Dates.Week)
            x -= Dates.Day(value(p) * 7)
        else
            x -= p
        end
    end
    return x
end


text:
# Parameters for the test script (set in runtests.jl)
#ENV["JULIA_DEBUG"] = "HierarchicalTemporalMemory"
plot_enabled= false

using HierarchicalTemporalMemory, BenchmarkTools, CSV, Printf, Lazy, Test
using Plots; gr()
import Random: seed!, bitrand

using StatsBase: median
seed!(0)

include("utils/utils.jl")

plot_mase(data,pred,pred_timesteps)= begin
  windowLength= 240
  errormetric= zeros(length(data)-windowLength)
  for w=1:length(data)-windowLength
    errormetric[w]= mase(data[w:w+windowLength-1], pred[w:w+windowLength-1], pred_timesteps)
  end
  crit_t1= 1:500; crit_t2= 3000:3200
  plot(plot(errormetric, label="10-day MASE"),
       plot([data pred], label=["timeseries","prediction"], legend=:none),
       plot(plot(crit_t1, [data[crit_t1] pred[crit_t1]], legend= :none),
            plot(crit_t2, [data[crit_t2] pred[crit_t2]], legend= :none),
            layout= (1,2)),
       layout= (3,1))|> display
  @info @sprintf("Min 10-day MASE: %.2f\n",minimum(errormetric))
end
display_evaluation(t,sp,sp_activity,spDims)= @info("t=$t")
process_data!(tN,data,encParams,sp,tm,decoder)=
  for t in 1:tN
    z,a,power_bucket= _process_sp(t,tN,data,encParams,sp,display_evaluation)
    A,Π,B= _process_tm(t,tN, tm,a)
    prediction= predict!(decoder,Π,power_bucket)
    likelyPred= reverse_simpleArithmetic(prediction,"highmean",encParams.power_p)
    history_enc[:,t]= z;    history_SP[:,t]= a
    history_TMout[:,t]= A;  history_TMpred[:,t]= Π
    history_decodedPred[:,t]= prediction; history_likelyPred[t]= likelyPred

    global avg_burst= ((t-1)*avg_burst+count(B)/length(B))/t
  end

prediction_timesteps=1
inputDims= ((15,6,3).*25,)
spDims= (1600,)
k= 8
@info "creating Spatial Pooler"
sp= SpatialPooler(SPParams(
      szᵢₙ= map(sum,inputDims), szₛₚ=spDims,
      γ=600,
      s=0.03,
      prob_synapse=0.85,
      θ_stimulus_activate=5,
      p⁺_01= 0.20,
      p⁻_01= 0.12,
      β=3,
      Tboost=350,
      enable_local_inhibit=false,
      enable_boosting=true))
@info "creating Temporal Memory"
tm= TemporalMemory(TMParams(
      Nc=prod(spDims),
      k=k,
      θ_stimulus_activate=14,
      θ_stimulus_learn=12,
      synapseSampleSize=35,
      p⁺_01=0.24,
      p⁻_01=0.08,
      LTD_p⁻_01= 0.012
     ))

Ncol= prod(spDims); Ncell= Ncol*k
# Define input data
data,tN= read_gympower("test_data/gym_power_benchmark-extended.csv")
encParams= initenc_powerDay(data.power_hourly_kw, data.hour, data.is_weekend,
                 encoder_size=inputDims[1], w=(34,35,35))
decoder= SDRClassifier(Ncell,encParams.power_p.buckets,
                  α=0.09, buffer_length=prediction_timesteps)
# Histories
history_enc= falses(map(sum,inputDims)|>prod,tN)
history_SP= falses(Ncol,tN)
history_TMout=  falses(Ncell,tN)
history_TMpred= falses(Ncell,tN)
history_decodedPred= zeros(encParams.power_p.buckets,tN)
history_likelyPred= zeros(tN)
avg_burst= 0

process_data!(tN,data,encParams,sp,tm,decoder)

errormetric= mase(data.power_hourly_kw[400:end], history_likelyPred[400:end],prediction_timesteps)
@info @sprintf("Prediction MASE: %.3f\n", errormetric)

avg_TMout_sparsity= mapslices(x->count(x)./length(x),history_TMout,dims=1)'|>median
plot_enabled && plot_mase(data.power_hourly_kw, history_likelyPred, prediction_timesteps)
@info @sprintf("avg_TMout_sparsity: %.3f%%\n", 100*avg_TMout_sparsity)
@info @sprintf("avg_burst: %.3f%%\n", 100*avg_burst)

# Sanity check: If this is not true, there's something very wrong with the temporal memory
@test errormetric < 1.15


text:
Pkg.add("Images")
Pkg.add("DataFrames")
using Images
using DataFrames

#typeData could be either "train" or "test.
#labelsInfo should contain the IDs of each image to be read
#The images in the trainResized and testResized data files
#are 20x20 pixels, so imageSize is set to 400.
#path should be set to the location of the data files.

function read_data(typeData, labelsInfo, imageSize, path)
 #Intialize x matrix
 x = zeros(size(labelsInfo, 1), imageSize)

 for (index, idImage) in enumerate(labelsInfo["ID"])
  #Read image file
  nameFile = "$(path)/$(typeData)Resized/$(idImage).Bmp"
  img = imread(nameFile)

  #Convert img to float values
  temp = float32sc(img)

  #Convert color images to gray images
  #by taking the average of the color scales.
  if ndims(temp) == 3
   temp = mean(temp.data, 1)
  end

  #Transform image matrix to a vector and store
  #it in data matrix
  x[index, :] = reshape(temp, 1, imageSize)
 end
 return x
end

imageSize = 400 # 20 x 20 pixel

#Set location of data files, folders
path = ...

#Read information about training data , IDs.
labelsInfoTrain = readtable("$(path)/trainLabels.csv")

#Read training matrix
xTrain = read_data("train", labelsInfoTrain, imageSize, path)

#Read information about test data ( IDs ).
labelsInfoTest = readtable("$(path)/sampleSubmission.csv")

#Read test matrix
xTest = read_data("test", labelsInfoTest, imageSize, path)

#Get only first character of string (convert from string to character).
#Apply the function to each element of the column "Class"
yTrain = map(x -> x[1], labelsInfoTrain["Class"])

#Convert from character to integer
yTrain = int(yTrain)

xTrain = xTrain'
xTest = xTest'

addprocs(3)

@everywhere function euclidean_distance(a, b)
 distance = 0.0
 for index in 1:size(a, 1)
  distance += (a[index]-b[index]) * (a[index]-b[index])
 end
 return distance
end

@everywhere function get_k_nearest_neighbors(x, i, k)
 nRows, nCols = size(x)
 imageI = Array(Float32, nRows)
 for index in 1:nRows
  imageI[index] = x[index, i]
 end
 imageJ = Array(Float32, nRows)
 distances = Array(Float32, nCols)
 for j in 1:nCols
  for index in 1:nRows
   imageJ[index] = x[index, j]
  end
  distances[j] = euclidean_distance(imageI, imageJ)
 end
 sortedNeighbors = sortperm(distances)
 kNearestNeighbors = sortedNeighbors[2:k+1]
 return kNearestNeighbors
end

@everywhere function assign_label(x, y, k, i)
 kNearestNeighbors = get_k_nearest_neighbors(x, i, k)
 counts = Dict{Int, Int}()
 highestCount = 0
 mostPopularLabel = 0
 for n in kNearestNeighbors
  labelOfN = y[n]
  if !haskey(counts, labelOfN)
   counts[labelOfN] = 0
  end
  counts[labelOfN] += 1
  if counts[labelOfN] > highestCount
   highestCount = counts[labelOfN]
   mostPopularLabel = labelOfN
  end
 end
 return mostPopularLabel
end

tic()
k = 1
sumValues = @parallel (+) for i in 1:size(xTrain, 2)
 assign_label(xTrain, yTrain, k, i) == yTrain[i, 1]
end
loofCvAccuracy = sumValues / size(xTrain, 2)
println("The LOOF-CV accuracy of 1NN is $(loofCvAccuracy)")
toc()

#Similar to function assign_label.
#Only changes are commented
@everywhere function assign_label_each_k(x, y, maxK, i)
 kNearestNeighbors = get_k_nearest_neighbors(x, i, maxK)

 #The next array will keep the labels for each value of k
 labelsK = zeros(Int, 1, maxK)

 counts = Dict{Int, Int}()
 highestCount = 0
 mostPopularLabel = 0

 #We need to keep track of the current value of k
 for (k, n) in enumerate(kNearestNeighbors)
  labelOfN = y[n]
  if !haskey(counts, labelOfN)
   counts[labelOfN] = 0
  end
  counts[labelOfN] += 1
  if counts[labelOfN] > highestCount
   highestCount = counts[labelOfN]
   mostPopularLabel = labelOfN
  end
  #Save current most popular label
  labelsK[k] = mostPopularLabel
 end
 #Return vector of labels for each k
 return labelsK
end

tic()
maxK = 20 #Any value can be chosen
yPredictionsK = @parallel (vcat) for i in 1:size(xTrain, 2)
 assign_label_each_k(xTrain, yTrain, maxK, i)
end
for k in 1:maxK
 accuracyK = mean(yTrain .== yPredictionsK[:, k])
 println("The LOOF-CV accuracy of $(k)-NN is $(accuracyK)")
end
toc()

@everywhere function get_k_nearest_neighbors(xTrain, imageI, k)
 nRows, nCols = size(xTrain)
 imageJ = Array(Float32, nRows)
 distances = Array(Float32, nCols)
 for j in 1:nCols
  for index in 1:nRows
   imageJ[index] = xTrain[index, j]
  end
  distances[j] = euclidean_distance(imageI, imageJ)
 end
 sortedNeighbors = sortperm(distances)
 kNearestNeighbors = sortedNeighbors[1:k]
 return kNearestNeighbors
end

@everywhere function assign_label(xTrain, yTrain, k, imageI)
 kNearestNeighbors = get_k_nearest_neighbors(xTrain, imageI, k)
 counts = Dict{Int, Int}()
 highestCount = 0
 mostPopularLabel = 0
 for n in kNearestNeighbors
  labelOfN = yTrain[n]
  if !haskey(counts, labelOfN)
   counts[labelOfN] = 0
  end
  counts[labelOfN] += 1 #add one to the count
  if counts[labelOfN] > highestCount
   highestCount = counts[labelOfN]
   mostPopularLabel = labelOfN
  end
 end
 return mostPopularLabel
end

println("Running kNN on test data")
tic()
k = 3 # The CV accuracy shows this value to be the best.
yPredictions = @parallel (vcat) for i in 1:size(xTest, 2)
 nRows = size(xTrain, 1)
 imageI = Array(Float32, nRows)
 for index in 1:nRows
  imageI[index] = xTest[index, i]
 end
 assign_label(xTrain, yTrain, k, imageI)
end
toc()

#Convert integer predictions to character
labelsInfoTest["Class"] = char(yPredictions)

#Save predictions
writetable("$(path)/juliaKNNSubmission.csv", labelsInfoTest, separator=',', header=true)

println("Submission file saved in $(path)/juliaKNNSubmission.csv")


text:
# ----------- #
# Blink setup #
# ----------- #

type ElectronDisplay <: AbstractPlotlyDisplay
    divid::Base.Random.UUID
    w::Nullable{Any}
    js_loaded::Bool
end

typealias ElectronPlot SyncPlot{ElectronDisplay}

ElectronDisplay(divid::Base.Random.UUID) = ElectronDisplay(divid, Nullable(), false)
ElectronDisplay(p::Plot) = ElectronDisplay(p.divid)
ElectronPlot(p::Plot) = ElectronPlot(p, ElectronDisplay(p.divid))

fork(jp::ElectronPlot) = ElectronPlot(fork(jp.plot), ElectronDisplay())

isactive(ed::ElectronDisplay) = isnull(ed.w) ? false : Blink.active(get(ed.w))

Base.close(ed::ElectronDisplay) = isactive(ed) && close(get(ed.w))

function get_window(p::ElectronPlot; kwargs...)
    w, h = size(p.plot)
    get_window(p.view; width=w, height=h, kwargs...)
end

function get_window(ed::ElectronDisplay; kwargs...)
    if !isnull(ed.w) && active(get(ed.w))
        w = get(ed.w)
    else
        w = get_window(Dict(kwargs))
        ed.w = w
        ed.js_loaded = false  # can't have js if we made a new window
    end
    w
end

get_window(opts::Dict) = Juno.isactive() ? Juno.Atom.blinkplot() : Window(opts)

js_loaded(ed::ElectronDisplay) = ed.js_loaded

function loadjs(ed::ElectronDisplay)
    if !ed.js_loaded
        Blink.loadjs!(get_window(ed), _mathjax_cdn_path)
        Blink.load!(get_window(ed), _js_path)
        ed.js_loaded = true
    end
end

Base.display(p::ElectronPlot; show = true, resize::Bool=_autoresize[1]) =
  display_blink(p; show = show, resize = resize)

function display_blink(p::ElectronPlot; show=true, resize::Bool=_autoresize[1])
    w = get_window(p, show=show)
    loadjs(p.view)
    lowered = JSON.lower(p.plot)
    if !resize
        @js w begin
            trydiv = document.getElementById($(string(p.plot.divid)))
            if trydiv == nothing
                @var gd = document.createElement("div")
                gd.id = $(string(p.plot.divid))
                document.body.appendChild(gd)
            else
                @var gd = trydiv
            end
            @var _ = Plotly.newPlot(gd, $(json(lowered[:data])),
            $(json(lowered[:layout])),
            d("showLink"=> false))
            _.then(()->Promise.resolve())
        end
    else
        # remove the width and height fields on the layout or
        # the autoresize won't work
        pop!(p.plot.layout.fields, :width, nothing)
        pop!(p.plot.layout.fields, :height, nothing)
        magic = """
        <script>
        (function() {
            var d3 = Plotly.d3
            var WIDTH_IN_PERCENT_OF_PARENT = 100
            var HEIGHT_IN_PERCENT_OF_PARENT = 100;
            var gd3 = d3.select('body')
            .append('div').attr("id", "$(p.plot.divid)")
            .style({
                width: WIDTH_IN_PERCENT_OF_PARENT + '%',
                'margin-left': (100 - WIDTH_IN_PERCENT_OF_PARENT) / 2 + '%',
                height: HEIGHT_IN_PERCENT_OF_PARENT + 'vh',
                'margin-top': (100 - HEIGHT_IN_PERCENT_OF_PARENT) / 2 + 'vh'
            });
            var gd = gd3.node();
            var data = $(json(p.plot.data));
            var layouts = $(json(p.plot.layout));
            Plotly.newPlot(gd, data, layouts);
            window.onresize = function() {
            Plotly.Plots.resize(gd);
            };
        })();
        </script>
        """
        Blink.body!(w, magic)
    end
    p.plot
end

## API Methods for ElectronDisplay
function _img_data(p::ElectronPlot, format::String; show::Bool=false)
    _formats = ["png", "jpeg", "webp", "svg"]
    if !(format in _formats)
        error("Unsupported format $format, must be one of $_formats")
    end

    opened_here = !isactive(p.view)
    opened_here && display(p; show=show, resize=false)

    out = @js p.view begin
        ev = Plotly.Snapshot.toImage(this, d("format"=>$format))
        @new Promise(resolve -> ev.once("success", resolve))
    end
    opened_here && close(p)
    out
end

function svg_data(p::ElectronPlot, format="png")
    opened_here = !isactive(p.view)
    opened_here && display(p; show=false, resize=false)

    out = @js p.view Plotly.Snapshot.toSVG(this, $format)

    opened_here && close(p)

    return out
end

Base.close(p::ElectronPlot) = close(p.view)

function Blink.js(p::ElectronDisplay, code::JSString; callback=true)
    if !isactive(p)
        return
    end

    Blink.js(get_window(p),
             :(Blink.evalwith(document.getElementById($(string(p.divid))),
                              $(Blink.jsstring(code)))),
             callback=callback)
end

Blink.js(p::ElectronPlot, code::JSString; callback=true) =
    Blink.js(p.view, code; callback=callback)

# Methods from javascript API (docstrings found in api.jl)
relayout!(p::ElectronDisplay, update::Associative=Dict(); kwargs...) =
    @js_ p Plotly.relayout(this, $(merge(update, prep_kwargs(kwargs))))

restyle!(p::ElectronDisplay, ind::Int, update::Associative=Dict(); kwargs...) =
    @js_ p Plotly.restyle(this, $(merge(update, prep_kwargs(kwargs))), $(ind-1))

restyle!(p::ElectronDisplay, inds::AbstractVector{Int}, update::Associative=Dict(); kwargs...) =
    @js_ p Plotly.restyle(this, $(merge(update, prep_kwargs(kwargs))), $(inds-1))

restyle!(p::ElectronDisplay, update=Dict(); kwargs...) =
    @js_ p Plotly.restyle(this, $(merge(update, prep_kwargs(kwargs))))

addtraces!(p::ElectronDisplay, traces::AbstractTrace...) =
    @js_ p Plotly.addTraces(this, $traces)

addtraces!(p::ElectronDisplay, where::Int, traces::AbstractTrace...) =
    @js_ p Plotly.addTraces(this, $traces, $(where-1))

deletetraces!(p::ElectronDisplay, traces::Int...) =
    @js_ p Plotly.deleteTraces(this, $(collect(traces)-1))

movetraces!(p::ElectronDisplay, to_end::Int...) =
    @js_ p Plotly.moveTraces(this, $(collect(to_end)-1))

movetraces!(p::ElectronDisplay, src::AbstractVector{Int}, dest::AbstractVector{Int}) =
    @js_ p Plotly.moveTraces(this, $(src-1), $(dest-1))

redraw!(p::ElectronDisplay) =
    @js_ p Plotly.redraw(this)

# unexported (by plotly.js) api methods
extendtraces!(ed::ElectronDisplay, update::Associative=Dict(),
              indices::Vector{Int}=[1], maxpoints=-1;) =
    @js_ ed Plotly.extendTraces(this, $(prep_kwargs(update)), $(indices-1), $maxpoints)

prependtraces!(ed::ElectronDisplay, update::Associative=Dict(),
               indices::Vector{Int}=[1], maxpoints=-1;) =
    @js_ ed Plotly.prependTraces(this, $(prep_kwargs(update)), $(indices-1), $maxpoints)


text:
import Oracle
using Dates

@assert VERSION >= v"0.7-"

@assert isfile(joinpath(@__DIR__, "..", "test","credentials.jl")) """
Before running tests, create a file `test/credentials.jl` with the content:

username = "your-username"
password = "your-password"
connect_string = "your-connect-string"
auth_mode = Oracle.ORA_MODE_AUTH_DEFAULT # or Oracle.ORA_MODE_AUTH_SYSDBA if user is SYSDBA
"""
include(joinpath(@__DIR__, "..", "test","credentials.jl"))

conn = Oracle.Connection(username, password, connect_string, auth_mode=auth_mode)

Oracle.execute(conn, "CREATE TABLE TB_BENCH_EXECUTE_MANY ( ID NUMBER(15,0) NULL, FLT NUMBER(15,4) NULL )")

let
    # JIT
    stmt = Oracle.Stmt(conn, "INSERT INTO TB_BENCH_EXECUTE_MANY ( ID, FLT ) VALUES ( :1, :2 )")
    stmt[1] = 1
    stmt[2] = 0.5
    Oracle.execute(stmt)
    Oracle.rollback(conn)
    Oracle.close(stmt)
end

let
    # JIT
    Oracle.execute(conn, "INSERT INTO TB_BENCH_EXECUTE_MANY ( ID, FLT ) VALUES ( :1, :2 )", [ [1, 2], [0.5, 1.5] ])
    Oracle.rollback(conn)
end

NUM_ROWS = 1_000

column_1 = [ i for i in 1:NUM_ROWS ]
column_2 = .5 * column_1

println("one row at a time")
@time let
    stmt = Oracle.Stmt(conn, "INSERT INTO TB_BENCH_EXECUTE_MANY ( ID, FLT ) VALUES ( :1, :2 )")
    for i in 1:NUM_ROWS
        stmt[1] = column_1[i]
        stmt[2] = column_2[i]
        Oracle.execute(stmt)
    end
    Oracle.rollback(conn)
    Oracle.close(stmt)
end

println("execute many")
@time let
    Oracle.execute(conn, "INSERT INTO TB_BENCH_EXECUTE_MANY ( ID, FLT ) VALUES ( :1, :2 )", [ column_1, column_2 ])
    Oracle.rollback(conn)
end

Oracle.execute(conn, "DROP TABLE TB_BENCH_EXECUTE_MANY")