text:
export Locations, CtrlLocations, AbstractLocations, LocationError, merge_locations

const LocationStorageTypes = Union{Int,NTuple{N,Int} where N,UnitRange{Int},Vector{Int}}

abstract type AbstractLocations end

"""
    merge_locations(locations...)

Construct a new `Locations` by merging two or more existing locations.
"""
merge_locations(x::AbstractLocations, y::AbstractLocations, locations::AbstractLocations...) =
    merge_locations(merge_locations(x, y), locations...)

"""
    Locations <: AbstractLocations

Type to annotate locations in quantum circuit.

    Locations(x)

Create a `Locations` object from a raw location statement. Valid storage types are:

- `Int`: single position
- `NTuple{N, Int}`: a list of locations
- `UnitRange{Int}`: contiguous locations

Other types will be converted to the storage type via `Tuple`.
"""
struct Locations{T<:LocationStorageTypes} <: AbstractLocations
    storage::T

    Locations(x::T) where {T<:LocationStorageTypes} = new{T}(x)
end

# skip it if x is a location
Locations(x::Locations) = x
Locations(xs...) = Locations(xs)
Locations(x::NTuple{N,T}) where {N,T} = throw(LocationError("expect Int, got $T"))

Base.@propagate_inbounds Base.getindex(l::Locations, idx...) = getindex(l.storage, idx...)
Base.length(l::Locations) = length(l.storage)
Base.iterate(l::Locations) = iterate(l.storage)
Base.iterate(l::Locations, st) = iterate(l.storage, st)
Base.eltype(::Type{T}) where {T<:Locations} = Int
Base.eltype(x::Locations) = Int
Base.show(io::IO, x::Locations) = print(io, "Locations(", x.storage, ")")
Base.Tuple(x::Locations) = (x.storage...,)

struct LocationError <: Exception
    msg::String
end

function merge_locations(l1::Locations, l2::Locations)
    Locations((l1.storage..., l2.storage...))
end

# location mapping
# TODO: preserve sign when indexing
# TODO: provide a @inlocation macro via Expr(:meta, :inlocation, true) so when we compile to Julia functions
#       we can use unsafe_mapping directly
@inline unsafe_mapping(parent::Locations{Int}, sub::Locations{Int}) = parent
@inline unsafe_mapping(parent::Locations{Int}, sub::Locations{NTuple{N,Int}}) where {N} = parent
@inline unsafe_mapping(parent::Locations{Int}, sub::Locations{UnitRange{Int}}) = parent
@inline unsafe_mapping(parent::Locations{NTuple{N,Int}}, sub::Locations{Int}) where {N} =
    Locations(@inbounds parent[sub.storage])
@inline unsafe_mapping(parent::Locations{NTuple{N,Int}}, sub::Locations{NTuple{M,Int}}) where {N,M} =
    Locations(map(x -> @inbounds(parent[x]), sub.storage))
@inline unsafe_mapping(parent::Locations{NTuple{N,Int}}, sub::Locations{UnitRange{Int}}) where {N} =
    Locations(@inbounds parent[sub.storage])
@inline unsafe_mapping(parent::Locations{UnitRange{Int}}, sub::Locations{Int}) =
    Locations(@inbounds parent[sub.storage])
@inline unsafe_mapping(parent::Locations{UnitRange{Int}}, sub::Locations{NTuple{N,Int}}) where {N} =
    Locations(map(x -> @inbounds(parent[x]), sub.storage))
@inline unsafe_mapping(parent::Locations{UnitRange{Int}}, sub::Locations{UnitRange{Int}}) =
    Locations(@inbounds parent[sub.storage])

map_error(parent, sub) = throw(LocationError("got $sub in parent space $parent"))

@inline function map_check(parent::Locations{Int}, sub::Locations{Int})
    sub.storage == 1 || map_error(parent, sub)
end

@inline function map_check(parent::Locations{Int}, sub::Locations{Tuple{Int}})
    sub.storage[1] == 1 || map_error(parent, sub)
end

@inline function map_check(parent::Locations{Int}, sub::Locations{NTuple{N,Int}}) where {N}
    map_error(parent, sub)
end

@inline function map_check(parent::Locations{Int}, sub::Locations{UnitRange{Int}})
    (length(sub) == 1) && (sub.storage.start == 1) || map_error(parent, sub)
end

@inline function map_check(parent::Locations{NTuple{N,Int}}, sub::Locations{Int}) where {N}
    1 <= sub.storage <= N || map_error(parent, sub)
end

@inline function map_check(
    parent::Locations{NTuple{N,Int}},
    sub::Locations{NTuple{M,Int}},
) where {N,M}
    all(x -> (1 <= x <= N), sub.storage) || map_error(parent, sub)
end

@inline function map_check(parent::Locations{NTuple{N,Int}}, sub::Locations{UnitRange{Int}}) where {N}
    (1 <= sub.storage.start) && (sub.storage.stop <= N) || map_error(parent, sub)
end

@inline function map_check(parent::Locations{UnitRange{Int}}, sub::Locations{Int})
    1 <= sub.storage <= length(parent) || map_error(parent, sub)
end

@inline function map_check(parent::Locations{UnitRange{Int}}, sub::Locations{NTuple{N,Int}}) where {N}
    all(x -> (1 <= x <= length(parent)), sub.storage) || map_error(parent, sub)
end

@inline function map_check(parent::Locations{UnitRange{Int}}, sub::Locations{UnitRange{Int}})
    (1 <= sub.storage.start) && (sub.storage.stop <= length(parent)) || map_error(parent, sub)
end

# comparing
function Base.:(==)(l1::Locations, l2::Locations)
    length(l1) == length(l2) || return false
    return all(l1.storage .== l2.storage)
end

# CtrlLocations
struct CtrlLocations{T<:LocationStorageTypes} <: AbstractLocations
    storage::Locations{T}
    configs::BitVector
end

# skip itself
CtrlLocations(x::CtrlLocations) = x
CtrlLocations(x::Locations) = CtrlLocations(x, trues(length(x)))
CtrlLocations(x::LocationStorageTypes, cfg::Tuple) =
    CtrlLocations(Locations(x), BitVector(map(Bool, cfg)))
CtrlLocations(xs...) = CtrlLocations(Locations(xs...))

Base.length(l::CtrlLocations) = length(l.storage)

function Base.show(io::IO, x::CtrlLocations)
    print(io, "CtrlLocations(")
    if all(x.configs)
        print(io, x.storage.storage)
    else
        join(io, map((l, c) -> c ? string(l) : "!" * string(l), x.storage.storage, x.configs), ", ")
    end
    print(io, ")")
end

function merge_locations(l1::CtrlLocations, l2::CtrlLocations)
    CtrlLocations(merge_locations(l1.storage, l2.storage), vcat(l1.configs, l2.configs))
end

# NOTE: CtrlLocations can not be mapped by Locations
@inline unsafe_mapping(parent::Locations, sub::CtrlLocations) =
    CtrlLocations(unsafe_mapping(parent, sub.storage), sub.configs)
@inline map_check(parent::Locations, sub::CtrlLocations) = map_check(parent, sub.storage)

Base.:(==)(l1::CtrlLocations, l2::CtrlLocations) =
    (l1.storage == l2.storage) && (l1.configs == l2.configs)

# parent location has to be Locations since CtrlLocations can't be parent
@inline function Base.getindex(parent::Locations, sub::AbstractLocations)
    map_check(parent, sub)
    return unsafe_mapping(parent, sub)
end


text:
## Day 19: Beacon Scanner ##
############################

# Wow, this is an interesting problem. Definitely not just "take known algorithm".

# Calculating the "signature" of a scanner can help me match scanners.
# This is the set of all differences of the beacons seen by the scanner.
# That is translation invariant, but still not rotation invariant.

# Rotation invariant would be the "square distance signature", so I better start
# with that.

function p(data::Vector{String})::Vector{Tuple{Int,Vector}}
    scanners = []
    current_scanner_id = nothing
    current_scanner_seen = []
    for line in data
        m = match(r"--- scanner ([0-9]+) ---", line)
        if m !== nothing
            if current_scanner_id !== nothing
                push!(scanners, (parse(Int, current_scanner_id), current_scanner_seen))
            end
            current_scanner_id = m[1]
            current_scanner_seen = []
        else
            m = match(r"(-?[0-9]+),(-?[0-9]+),(-?[0-9]+)", line)
            if m !== nothing
                push!(current_scanner_seen, (parse(Int, m[1]), parse(Int, m[2]), parse(Int, m[3])))
            end
        end
    end
    push!(scanners, (parse(Int, current_scanner_id), current_scanner_seen))
    scanners
end

test = p(readlines(open("test-19")))
input = p(readlines(open("input-19")))

function sig(sensor_data)::Set{Tuple{Int,Int,Int}}
    sig = Set()
    # To make sure there is no issue with the direction the deltas are pointing
    # to, I include them in both directions.
    for b1 in sensor_data, b2 in sensor_data
        if b1 == b2
            continue
        end
        candidate = (b1[1] - b2[1], b1[2] - b2[2], b1[3] - b2[3])
        push!(sig, candidate)
    end
    sig
end

rot_x_r1((x, y, z)) = (x, y, z)
rot_x_r2((x, y, z)) = (x, -z, y)
rot_x_r3((x, y, z)) = (x, -y, -z)
rot_x_r4((x, y, z)) = (x, z, -y)
rot_y_r1((x, y, z)) = (z, x, y)
rot_y_r2((x, y, z)) = (-y, x, z)
rot_y_r3((x, y, z)) = (-z, x, -y)
rot_y_r4((x, y, z)) = (y, x, -z)
rot_z_r1((x, y, z)) = (y, z, x)
rot_z_r2((x, y, z)) = (-z, y, x)
rot_z_r3((x, y, z)) = (-y, -z, x)
rot_z_r4((x, y, z)) = (z, -y, x)

function rotate_all(sig::Set)::Vector
    # Rotations can be specified by showing where (1,0,0) maps to + 1 of 4 rotations
    # Since we already put both signs in the signature, we can restrict x to
    # point in only 3 directions.
    result = []
    push!(result, Set(rot_x_r1.(sig)))
    push!(result, Set(rot_x_r2.(sig)))
    push!(result, Set(rot_x_r3.(sig)))
    push!(result, Set(rot_x_r4.(sig)))
    push!(result, Set(rot_y_r1.(sig)))
    push!(result, Set(rot_y_r2.(sig)))
    push!(result, Set(rot_y_r3.(sig)))
    push!(result, Set(rot_y_r4.(sig)))
    push!(result, Set(rot_z_r1.(sig)))
    push!(result, Set(rot_z_r2.(sig)))
    push!(result, Set(rot_z_r3.(sig)))
    push!(result, Set(rot_z_r4.(sig)))
end

function look_for_matches(data)
    # Look for possible matches
    for i = 1:length(data)
        sig1_base = sig(data[i][2])
        # Iterate all 12 orientations
        for sig1 in rotate_all(sig1_base)
            for j = (i+1):length(data)
                sig2 = sig(data[j][2])
                @show i, j, length(intersect(sig1, sig2))
            end
        end
    end
end

ext: jl


text:
#####
##### Memory Buffer:
##### Datastructure to collect self-play experience
#####

"""
    TrainingSample{State}

Type of a training sample. A sample features the following fields:
- `s::State` is the state
- `π::Vector{Float64}` is the recorded MCTS policy for this position
- `z::Float64` is the discounted reward cumulated from state `s`
- `t::Float64` is the (average) number of moves remaining before the end of the game
- `n::Int` is the number of times the state `s` was recorded

As revealed by the last field `n`, several samples that correspond to the
same state can be merged, in which case the `π`, `z` and `t`
fields are averaged together.
"""
struct TrainingSample{State}
  s :: State
  π :: Vector{Float64}
  z :: Float64
  t :: Float64
  n :: Int
end

sample_state_type(::Type{<:TrainingSample{S}}) where S = S

"""
    MemoryBuffer{State}

A circular buffer to hold memory samples.

# Constructor

    MemoryBuffer{State}(size, experience=[])

"""
mutable struct MemoryBuffer{State}
  buf :: CircularBuffer{TrainingSample{State}}
  cur_batch_size :: Int
  function MemoryBuffer{S}(size, experience=[]) where S
    buf = CircularBuffer{TrainingSample{S}}(size)
    append!(buf, experience)
    new{S}(buf, 0)
  end
end

"""
    get_experience(::MemoryBuffer{S}) where S :: Vector{TrainingSample{S}}

Return all samples in the memory buffer.
"""
get_experience(mem::MemoryBuffer) = mem.buf[:]

last_batch(mem::MemoryBuffer) = mem.buf[end-cur_batch_size(mem)+1:end]

cur_batch_size(mem::MemoryBuffer) = min(mem.cur_batch_size, length(mem))

new_batch!(mem::MemoryBuffer) = (mem.cur_batch_size = 0)

function Base.empty!(mem::MemoryBuffer)
  empty!(mem.buf)
  mem.cur_batch_size = 0
end

Base.length(mem::MemoryBuffer) = length(mem.buf)

"""
    push_game!(mem::MemoryBuffer, trace::Trace, gamma)

Collect samples out of a game trace and add them to the memory buffer.

Here, `gamma` is the reward discount factor.
"""
function push_game!(mem::MemoryBuffer, trace, gamma)
  n = length(trace)
  wr = 0.
  for i in reverse(1:n)
    wr = gamma * wr + trace.rewards[i]
    s = trace.states[i]
    π = trace.policies[i]
    wp = GI.white_playing(GameType(trace), s)
    z = wp ? wr : -wr
    t = float(n - i + 1)
    push!(mem.buf, TrainingSample(s, π, z, t, 1))
  end
  mem.cur_batch_size += n
end

function merge_samples(samples)
  s = samples[1].s
  π = mean(e.π for e in samples)
  z = mean(e.z for e in samples)
  n = sum(e.n for e in samples)
  t = mean(e.t for e in samples)
  return eltype(samples)(s, π, z, t, n)
end

# Merge samples that correspond to identical states
function merge_by_state(samples)
  Sample = eltype(samples)
  State = sample_state_type(Sample)
  dict = Dict{State, Vector{Sample}}()
  sizehint!(dict, length(samples))
  for s in samples
    if haskey(dict, s.s)
      push!(dict[s.s], s)
    else
      dict[s.s] = [s]
    end
  end
  return [merge_samples(ss) for ss in values(dict)]
end

function apply_symmetry(Game, sample, (symstate, aperm))
  mask = GI.actions_mask(Game(sample.s))
  symmask = GI.actions_mask(Game(symstate))
  π = zeros(eltype(sample.π), length(mask))
  π[mask] = sample.π
  π = π[aperm]
  @assert iszero(π[.~symmask])
  π = π[symmask]
  return typeof(sample)(
    symstate, π, sample.z, sample.t, sample.n)
end

function augment_with_symmetries(Game, samples)
  symsamples = [apply_symmetry(Game, s, sym)
    for s in samples for sym in GI.symmetries(Game, s.s)]
  return [samples ; symsamples]
end


text:
export ruleSPMomentConstraintOutG

"""
Evaluate the expectation ∫q(x; η)g(x)dx for a given eta
"""
function expectedValue(eta::Float64, m_bw::Float64, V_bw::Float64, g::Function)
    Z_tilde = gaussianQuadrature(x->exp(eta*g(x)),              m=m_bw, v=V_bw) # Obtain normalizing constant for constrained q(x)
    G_tilde = gaussianQuadrature(x->exp(eta*g(x))*g(x)/Z_tilde, m=m_bw, v=V_bw) # Evaluate the expectation of g(x) w.r.t. the constrained q(x)

    return G_tilde
end

"""
Obtain the first and second moment of the constrained q(x; η) for a given η
"""
function constrainedMoments(eta::Float64, m_bw::Float64, V_bw::Float64, g::Function)
    Z_tilde    = gaussianQuadrature(x->exp(eta*g(x)),             m=m_bw, v=V_bw) # Obtain normalizing constant for constrained q(x)
    m_tilde    = gaussianQuadrature(x->exp(eta*g(x))*x/Z_tilde,   m=m_bw, v=V_bw) # Compute mean
    V_tilde_nc = gaussianQuadrature(x->exp(eta*g(x))*x^2/Z_tilde, m=m_bw, v=V_bw) # Compute non-central variance
    V_tilde    = V_tilde_nc - m_tilde^2 # Compute central variance

    return (m_tilde, V_tilde)
end

#-------------
# Update rules
#-------------

function ruleSPMomentConstraintOutG(msg_out::Message{<:Gaussian, Univariate}, g::Function, G::Float64, eta_init::Float64)
    (m_bw, V_bw) = unsafeMeanCov(msg_out.dist)
    (xi_bw, W_bw) = unsafeWeightedMeanPrecision(msg_out.dist)

    eta_star = find_zero(eta->expectedValue(eta, m_bw, V_bw, g)-G, eta_init) # Find optimal eta
    (m_tilde, V_tilde) = constrainedMoments(eta_star, m_bw, V_bw, g) # Compute moments of constrained q(x) with optimal eta

    # Convert approximated marginal to canonical form
    W_tilde = inv(V_tilde)
    xi_tilde = W_tilde*m_tilde

    return Message(Univariate, GaussianWeightedMeanPrecision, xi=xi_tilde-xi_bw, w=W_tilde-W_bw) # Compute forward message
end

#---------------------------
# Custom inbounds collectors
#---------------------------

function collectSumProductNodeInbounds(node::MomentConstraint, entry::ScheduleEntry)
    inbounds = Any[]

    interface_to_schedule_entry = current_inference_algorithm.interface_to_schedule_entry
    for node_interface in node.interfaces
        inbound_interface = ultimatePartner(node_interface)
        if isa(inbound_interface.node, Clamp)
            # Hard-code outbound message of constant node in schedule
            push!(inbounds, assembleClamp!(inbound_interface.node, Message))
        else
            # Collect message from previous result
            push!(inbounds, interface_to_schedule_entry[inbound_interface])
        end
    end

    # Push function and value to calling signature; these need to be defined in the scope of the user
    push!(inbounds, Dict{Symbol, Any}(:g => node.g,
                                      :keyword => false))
    push!(inbounds, Dict{Symbol, Any}(:G => node.G,
                                      :keyword => false))
    push!(inbounds, Dict{Symbol, Any}(:eta_init => node.eta_init,
                                      :keyword => false))

    return inbounds
end


text:
using MLDatasets
using BSON
using BSON: @load
using Plots
using Statistics

function prepare_data(file_name)
    X = Iris.features()
    y_orig = Iris.labels()

    y_name = unique(y_orig)
    y = zeros(length(y_orig))
    for i in 1:length(y_name)
        y[y_orig .== y_name[i]] .= i
    end

    X = Matrix(X')
    y = Int64.(y[:])

    BSON.bson(file_name, X=X, y=y, y_name=y_name)
end

function normalize_data(X)
    col_mean = mean(X, dims=1)
    col_std = std(X, dims=1)
    return (X .- col_mean) ./ col_std
end

file_name = joinpath("src", "lecture_08", "data", "iris.bson")
prepare_data(file_name)

max_stars_repo_name: rzehumat/JuliaCourse

max_issues_repo_name: rzehumat/JuliaCourse

max_forks_repo_name: rzehumat/JuliaCourse


text:
## LICENSE
##   Copyright (c) 2018 GEOPHYSTECH LLC
##
##   Licensed under the Apache License, Version 2.0 (the "License");
##   you may not use this file except in compliance with the License.
##   You may obtain a copy of the License at
##
##       http://www.apache.org/licenses/LICENSE-2.0
##
##   Unless required by applicable law or agreed to in writing, software
##   distributed under the License is distributed on an "AS IS" BASIS,
##   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
##   See the License for the specific language governing permissions and
##   limitations under the License.
##
## initial release by Andrey Stepnov, email: a.stepnov@geophsytech.ru

"""
**ON GRID**

- On grid without Dl data:
```
gmpe_mf2013(eq::Earthquake,config::Params_mf2013,grid::Array{Point_vs30};min_val::Number=0,Dl::Number=250,Xvf::Number=0)
```
- On grid with Dl data
```
gmpe_mf2013(eq::Earthquake,config::Params_mf2013,grid::Array{Point_vs30_dl};min_val::Number=0,Xvf::Number=0)
```
where `min_val=0`, `Xvf=0` [km] by default. `Dl=250` [km] by default in case of grid pass without `Dl` data.

Output will be 1-d `Array{Point_<ground_motion_type>_out}` with points based on input grid with `ground_motion > min_val` (`pga`,`psa` is Acceleration of gravity in percent (%g) rounded to ggg.gg, `pgv` in cm/s^2)

**without grid**
```
gmpe_as2008(eq::Earthquake,config::Params_mf2013;VS30::Number=350,distance::Int64=1000,Dl::Number=250,Xvf::Number=0)
```
where `VS30=30` [m/s^2], `distance=1000` [km], `Xvf=0` [km], `Dl=250` [km] by default.

Output will be 1-d `Array{Float64}` with `1:distance` values of `pga`,`psa` (that is Acceleration of gravity in percent (%g) rounded to ggg.gg) OR pgv [cm/s].

**EXAMPLES:**
```
gmpe_mf2013(eq,config_mf2013,grid) # for PGA,PGV,PSA on GRID
gmpe_mf2013(eq,config_mf2013) # for PGA,PGV,PSA without input grid
```

**Model parameters**

Please, see `examples/morikawa-fujiwara-2013.conf`
"""
## Morikawa Fujiwara (2013) PGA,PGV,PSA modeling ON GRID, Dl = constant
function gmpe_mf2013(eq::Earthquake,config::Params_mf2013,grid::Array{Point_vs30};min_val::Number=0,Dl::Number=250,Xvf::Number=0)
  vs30_row_num = length(grid[:,1])
  # define magnitude and epicenter
  eq.moment_mag == 0 ? magnitude = eq.local_mag : magnitude = eq.moment_mag
  magnitude = min(magnitude,config.Mw0)
  epicenter = LatLon(eq.lat, eq.lon)
  # define g_global
  g_global = 9.81
  # define a*Mw'
  aMw = config.a*magnitude
  # define Gd out of loop by grid points
  Gd = config.pd*log10(max(config.Dlmin,Dl)/config.D0)
  # init output_data
  if config.ground_motion_type == "PGA"
    output_data = Array{Point_pga_out}(undef,0)
    out_type = Point_pga_out
  elseif config.ground_motion_type == "PGV"
    output_data = Array{Point_pgv_out}(undef,0)
    out_type = Point_pgv_out
  elseif config.ground_motion_type == "PSA"
    output_data = Array{Point_psa_out}(undef,0)
    out_type = Point_psa_out
  end
  # main cycle by grid points
  for i=1:vs30_row_num
    # rrup the same as X in Morikawa Fujiwara 2013 formulae
    # eq.depth the same as D in Morikawa Fujiwara 2013 formulae
    current_point = LatLon(grid[i].lat,grid[i].lon)
    r_rup = sqrt((distance(current_point,epicenter)/1000)^2 + eq.depth^2)
    # \logA where A in cm/s^2 (pga,psa) or cm/s (pgv)
    log_A = aMw + config.b*r_rup + config.c - log10(r_rup + config.d*10^(config.e*magnitude))
    # Amplification by Deep Sedimentary Layers
    log_Agd = log_A + Gd
    # vs30 amplification
    Gs = config.ps*log10(min(config.Vsmax,grid[i].vs30)/config.V0)
    log_Ags = log_Agd + Gs
    # ASID
    if config.ASID == true
      AI = config.gamma*Xvf*(eq.depth - 30)
      log_AI = log_Ags + AI
      A = 10^(log_AI) # pga and psa in cm/c^2, pgv in cm/s
    else
      A = 10^(log_Ags)
    end
    # output depend on type of motion
    if config.ground_motion_type == "PGA" || config.ground_motion_type == "PSA"
      motion = round(((A/100)/g_global * 100),digits=2) ## convert cm/c^2 to %g
    elseif config.ground_motion_type == "PGV"
      motion = A
    end
    if motion >= min_val
      output_data = push!(output_data, out_type(grid[i].lon,grid[i].lat,motion))
    end
  end
  return output_data
end
## Morikawa Fujiwara (2013) PGA,PGV,PSA modeling ON GRID, Dl provided in GRID
function gmpe_mf2013(eq::Earthquake,config::Params_mf2013,grid::Array{Point_vs30_dl};min_val::Number=0,Xvf::Number=0)
  vs30_row_num = length(grid[:,1])
  # define magnitude and epicenter
  eq.moment_mag == 0 ? magnitude = eq.local_mag : magnitude = eq.moment_mag
  magnitude = min(magnitude,config.Mw0)
  epicenter = LatLon(eq.lat, eq.lon)
  # define g_global
  g_global = 9.81
  # define a*Mw'
  aMw = config.a*magnitude
  # init output_data
  if config.ground_motion_type == "PGA"
    output_data = Array{Point_pga_out}(undef,0)
    out_type = Point_pga_out
  elseif config.ground_motion_type == "PGV"
    output_data = Array{Point_pgv_out}(undef,0)
    out_type = Point_pgv_out
  elseif config.ground_motion_type == "PSA"
    output_data = Array{Point_psa_out}(undef,0)
    out_type = Point_psa_out
  end
  # main cycle by grid points
  for i=1:vs30_row_num
    # rrup the same as X in Morikawa Fujiwara 2013 formulae
    # eq.depth the same as D in Morikawa Fujiwara 2013 formulae
    current_point = LatLon(grid[i].lat,grid[i].lon)
    r_rup = sqrt((distance(current_point,epicenter)/1000)^2 + eq.depth^2)
    # \logA where A in cm/s^2 (pga,psa) or cm/s (pgv)
    log_A = aMw + config.b*r_rup + config.c - log10(r_rup + config.d*10^(config.e*magnitude))
    # Amplification by Deep Sedimentary Layers
    Gd = config.pd*log10(max(config.Dlmin,grid[i].dl)/config.D0)
    log_Agd = log_A + Gd
    # vs30 amplification
    Gs = config.ps*log10(min(config.Vsmax,grid[i].vs30)/config.V0)
    log_Ags = log_Agd + Gs
    # ASID
    if config.ASID == true
      AI = config.gamma*Xvf*(eq.depth - 30)
      log_AI = log_Ags + AI
      A = 10^(log_AI) # pga and psa in cm/c^2, pgv in cm/s
    else
      A = 10^(log_Ags)
    end
    # output depend on type of motion
    if config.ground_motion_type == "PGA" || config.ground_motion_type == "PSA"
      motion = round(((A/100)/g_global * 100),digits=2) ## convert cm/c^2 to %g
    elseif config.ground_motion_type == "PGV"
      motion = A
    end
    if motion >= min_val
      output_data = push!(output_data, out_type(grid[i].lon,grid[i].lat,motion))
    end
  end
  return output_data
end
## Morikawa Fujiwara (2013) PGA modeling ON GRID, Dl = constant
function gmpe_mf2013(eq::Earthquake,config::Params_mf2013;VS30::Number=350,distance::Number=1000,Dl::Number=250,Xvf::Number=0)
  # define magnitude
  eq.moment_mag == 0 ? magnitude = eq.local_mag : magnitude = eq.moment_mag
  magnitude = min(magnitude,config.Mw0)
  # define g_global
  g_global = 9.81
  # define a*Mw'
  aMw = config.a*magnitude
  # define Gd out of loop
  Gd = config.pd*log10(max(config.Dlmin,Dl)/config.D0)
  # init output array
  output_data = Array{Float64}(undef,0)
  # main cycle
  for i=1:distance
    # rrup the same as X in Morikawa Fujiwara 2013 formulae
    # eq.depth the same as D in Morikawa Fujiwara 2013 formulae
    r_rup = sqrt((i)^2 + eq.depth^2)
    # \logA where A in cm/s^2 (pga,psa) or cm/s (pgv)
    log_A = aMw + config.b*r_rup + config.c - log10(r_rup + config.d*10^(config.e*magnitude))
    # Amplification by Deep Sedimentary Layers
    log_Agd = log_A + Gd
    # vs30 amplification
    Gs = config.ps*log10(min(config.Vsmax,VS30)/config.V0)
    log_Ags = log_Agd + Gs
    # ASID
    if config.ASID == true
      AI = config.gamma*Xvf*(eq.depth - 30)
      log_AI = log_Ags + AI
      A = 10^(log_AI) # pga and psa in cm/c^2, pgv in cm/s
    else
      A = 10^(log_Ags)
    end
    # output depend on type of motion
    if config.ground_motion_type == "PGA" || config.ground_motion_type == "PSA"
      motion = round(((A/100)/g_global * 100),digits=2) ## convert cm/c^2 to %g
    elseif config.ground_motion_type == "PGV"
      motion = A
    end
    output_data = push!(output_data, motion)
  end
  return output_data
end


text:
SCALE = 12
EdgesPerVertex = 16

Nmax = 2 ^ SCALE

M = EdgesPerVertex * Nmax

Nfile = 8

myFiles = 1:Nfile

for i = myFiles
    tic()
    fname = "data/" * string(i)

    srand(i)#Reset Rand Seed

    v1,v2 = KronGraph500NoPerm(SCALE,round(Int,EdgesPerVertex/Nfile))

    rowStr = join(v1,',')

    colStr = join(v2,',')
    valStr = join(Array{Int}(ones(length(v1))),',')

    fidRow = open(fname * "r.txt","w")
    fidCol = open(fname * "c.txt","w")
    fidVal = open(fname * "v.txt","w")

    write(fidRow,rowStr * ",")
    close(fidRow)

    write(fidCol,colStr * ",")
    close(fidCol)

    write(fidVal,valStr * ",")
    close(fidVal)
    pTime = toq()
    println("Time: ", pTime)
    println("Edges/sec: ", string(NumStr(rowStr) / pTime))
end


text:
include("cimgui_h.jl")
include("cimgui.jl")

function render()
	x::Array{Bool,1} = [true]
	size::ImVec2 = ImVec2(0,0)
	igBegin("window", pointer(x), 0)
	igButton("BUTTON2", ImVec2(100,100))
	igButton("BUTTON", ImVec2(0,0))
	igEnd()
end

setupAll()
setRenderCallback(render)
mainLoop()
destroyAll()


text:
# Recreates the analysis from the original data files

datapath = "dataset"
projectname = "MEG3"
experimentdirnames = ["FKBP12-mTOR"]

samplenames1 = ["E", "G"]
samplenames23 = ["E", "F"]
samplenames4 = ["C", "F"]
samplenames5 = ["A", "B"]

nreplicates = 5
ncells = 10

mc_iterations = 10000

outputdir = "output"
mkpath(outputdir)
outputdatapath = joinpath(outputdir, "resultspos.jld2")

using Distributed
currentworkers = addprocs(exeflags = "--project")
@everywhere using SMLMAssociationAnalysis_NCB
using Printf
using LocalizationMicroscopy
using Statistics
using FileIO

experimentresults = Vector{Vector{Vector{Result}}}[]
for experimentdirname ∈ experimentdirnames
    println("Starting experiment $experimentdirname.")
    experimentpath = joinpath(datapath, experimentdirname)
    replicateresults = Vector{Vector{Result}}[]
    for i ∈ 1:nreplicates
        sampleresults = Vector{Result}[]
        println("    Starting replicate $i.")
        replicatepath = joinpath(experimentpath, "Replicate $i")
        samplenames = i == 5 ? samplenames5 :
                      i == 4 ? samplenames4 :
                      i == 1 ? samplenames1 :
                               samplenames23
        for samplename ∈ samplenames
            results = Result[]
            println("        Starting sample $samplename.")
            for j ∈ 1:ncells
                jcell = j
                if i == 1 && samplename == "E"
                    jcell += 10
                end

                println("            Starting cell $jcell.")
                cellpath = joinpath(replicatepath, "$samplename $(Printf.@sprintf("%03i", jcell)).bin.txt")
                localizations = loadlocalizations(cellpath, LocalizationMicroscopy.nikonelementstext)
                ch1_name = i == 2 ? "488" : "647"
                ch2_name = i == 2 ? "647" : "488"  # the channels were reversed for replicate 2. Unsure if this caused other issues.

                ch1_startframe = 1
                ch2_startframe = 11001

                ch1_task = remotecall(getmolecules, currentworkers[1], localizations,
                    ch1_name,
                    ch1_startframe,
                    11000,
                    100,
                    10,
                    34.2,
                    500,
                    200,
                )
                ch2_task = remotecall(getmolecules, currentworkers[2], localizations,
                    ch2_name,
                    ch2_startframe,
                    11000,
                    100,
                    10,
                    34.2,
                    500,
                    200,
                )

                if i == 2
                    ch1_task, ch2_task = ch2_task, ch1_task
                end

                ch1_molecules, ch1_localizations = fetch(ch1_task)
                ch2_molecules, ch2_localizations = fetch(ch2_task)

                ch1_neighbors, ch2_neighbors, distances = exclusivenearestneighbors(ch1_molecules, ch2_molecules)

                percentileranks = montecarloaffinity(
                    ch1_molecules,
                    ch2_molecules,
                    ch1_neighbors,
                    ch2_neighbors,
                    distances,
                    800,
                    mc_iterations
                )

                positivecontrol_percentileranks, positivecontrol_distances = simulate100(ch1_molecules, ch2_molecules, ch1_neighbors, ch2_neighbors, 80, 800, mc_iterations)
                # I found that the variance in running the negative control drops smoothly as a function of neighbor counts.
                # This heuristic cuts down on computational time without sacrificing quality of the correction.
                negative_iterations = length(distances) < 250 ? 30 :
                                      length(distances) < 500 ? 20 :
                                      length(distances) < 1000 ? 10 :
                                      length(distances) < 2000 ? 5 :
                                      1
                negativecontrol = [simulate0(ch1_molecules, ch2_molecules, 800, 10000) for i in 1:negative_iterations]
                negativecontrol_percentileranks = [first(i) for i in negativecontrol]
                negativecontrol_distances = [last(i) for i in negativecontrol]

                mediandistance = length(distances) == 0 ? NaN : median(distances)

                println("                Done: $(length(distances)) neighbors from $(length(ch1_molecules)) and $(length(ch2_molecules)) molecules, $(length(ch1_localizations)) and $(length(ch2_localizations)) localizations; median distance $mediandistance")
                ch1_data = ChannelData(ch1_name, ch1_molecules, ch1_neighbors)
                ch2_data = ChannelData(ch2_name, ch2_molecules, ch2_neighbors)
                result = Result(projectname, experimentdirname, 1, cellpath, i,
                                [ch1_data, ch2_data], distances, mediandistance, percentileranks,
                                positivecontrol_distances, negativecontrol_distances,
                                positivecontrol_percentileranks, negativecontrol_percentileranks)

                push!(results, result)
            end
            push!(sampleresults, results)
        end
        push!(replicateresults, sampleresults)
    end
    push!(experimentresults, replicateresults)
end

save(outputdatapath, "experimentresults", experimentresults)

rmprocs(currentworkers)


text:
"""
    parameters(type)

Extracts the type-parameters of the `type`.

e.g. `parameters(Foo{A, B, C}) == [A, B, C]`
"""
parameters(sig::UnionAll) = parameters(sig.body)
parameters(sig::DataType) = sig.parameters
parameters(sig::Union) = Base.uniontypes(sig)


text:
# Magnetic field topology in fusion research.

using Elliptic, SpecialFunctions

"""
    getB_bottle(x, y, z, distance, a, b, I1, I2) -> Vector{Float}

Get magnetic field from a magnetic bottle.
Reference: https://en.wikipedia.org/wiki/Magnetic_mirror#Magnetic_bottles

# Arguments
- `x,y,z::Float`: location in [m].
- `distance::Float`: distance between solenoids in [m].
- `a::Float`: radius of each side coil in [m].
- `b::Float`: radius of central coil in [m].
- `I1::Float`: current in the solenoid times number of windings in side coils.
- `I2::Float`: current in the central solenoid times number of windings in the
central loop.
"""
function getB_bottle(x, y, z, distance, a, b, I1, I2)

   r = √(x^2 + y^2) # distance from z-axis

   # 1st loop
   z₁ = z + 0.5*distance
   k = √(4*r*a / (z₁^2 + (a + r)^2) )
   K, E = ellipke(k)
   Bz1 = μ₀*I1 / (2π*√(z₁^2+(a+r)^2)) * ((a^2-z₁^2-r^2)/(z₁^2+(r-a)^2)*E + K)
   Br1 = μ₀*z₁*I1/(2π*r*√(z₁^2+(a+r)^2))*((z₁^2+r^2+a^2)/(z₁^2+(r-a)^2)*E - K)
   Bx1 = Br1 * x / r
   By1 = Br1 * y / r

   # 2nd loop
   z₂ = z - 0.5*distance
   k = √(4*r*a / (z₂^2 + (a + r)^2) )
   K, E = ellipke(k)
   Bz2 = μ₀*I1 / (2π*√(z₂^2+(a+r)^2)) * ((a^2-z₂^2-r^2)/(z₂^2+(r-a)^2)*E + K)
   Br2 = μ₀*z₂*I1/(2π*r*√(z₂^2+(a+r)^2))*((z₂^2+r^2+a^2)/(z₂^2+(r-a)^2)*E - K)
   Bx2 = Br2 * x / r
   By2 = Br2 * y / r

   # central loop
   z₃ = z
   k = √(4*r*b / (z₃^2 + (b+r)^2) )
   K, E = ellipke(k)
   Bz3 = μ₀*I2 / (2π*√(z₃^2+(b+r)^2)) * ((b^2-z₃^2-r^2)/(z₃^2+(r-b)^2)*E + K)
   Br3 = μ₀*z₃*I2/(2π*r*√(z₃^2+(b+r)^2))*((z₃^2+r^2+b^2)/(z₃^2+(r-b)^2)*E - K)
   Bx3 = Br3 * x / r
   By3 = Br3 * y / r

   # total magnetic field
   if x == 0.0 && y == 0.0
      Bx = 0.0
      By = 0.0
      Bz = Bz1 + Bz2 + Bz3
   else
      Bx = Bx1 + Bx2 + Bx3
      By = By1 + By2 + By3
      Bz = Bz1 + Bz2 + Bz3
   end
   [Bx, By, Bz]
end

"""
    getB_tokamak(x, y, z, a, b, ICoils, IPlasma)

Get the magnetic field from a Tokamak topology consists of 16 coils.
Original: https://github.com/BoschSamuel/Simulation-of-a-Tokamak-Fusion-Reactor/blob/master/Simulation2.m
# Arguments
- `x,y,z::Float`: location in [m].
- `a::Float`: radius of each coil in [m].
- `b::Float`: radius of central region in [m].
- `ICoil::Float`: current in the coil times number of windings.
- `IPlasma::Float`: current of the plasma?
"""
function getB_tokamak(x, y, z, a, b, ICoils, IPlasma)

   a *= 2

   Bx, By, Bz = 0.0, 0.0, 0.0

   # magnetic field of the coils
   for i = 0:15
      θ = π/16 + i*π/8 # angle between the i-th coil and the x-axis

      if abs(sin(θ)) > 0.01
         r1_ = x/(cos(θ)-sin(θ)*tan(atan(y, x) - θ)) - a - b
         r1 = √(r1_^2 + z^2)
         z1 = ((b + a)*cos(θ) + r1_*cos(θ) - x) / sin(θ)
      else
         r1_ = x - b - a
         r1 = √(r1_^2 + z^2)
         z1 = y
      end

      k = √(4r1*a / (z1^2 + (a + r1)^2))

      K, E = ellipke(k)
      # Bz1_ is the magnetic field in the coil frame
      Bz1_ = μ₀*ICoils/(2π*√(z1^2+(a+r1)^2))*((a^2-z1^2-r1^2)/(z1^2+(r1-a)^2)*E+K)
      # Br1_ is the magnetic field in the coil frame
      Br1_ = μ₀*z1*ICoils/(2π*r1*√(z1^2+(a+r1)^2))*((z1^2+r1^2+a^2)/(z1^2+(r1-a)^2)*E-K)
      # normal coordinates
      Bx1 = -sin(θ)*Bz1_ + Br1_*r1_/r1*cos(θ)
      By1 = cos(θ)*Bz1_ + sin(θ)*Br1_*r1_/r1
      Bz1 = Br1_ * z / r1

      # add the field of a single coil to the total field
      Bx += Bx1
      By += By1
      Bz += Bz1

      if abs(Bx) < 5e-12
         Bx = 0.0
      end
      if abs(By) < 5e-12
         By = 0.0
      end
      if abs(Bz) < 5e-12
         Bz = 0.0
      end
   end

   # magnetic field of the plasma current
   σ = a/3 # parameter of the Gauss curve
   ϕ = atan(y, x)
   # distance to centre of plasma ring
   distance = √( z^2 + (x - (a + b)*cos(ϕ))^2 + (y - (a + b)*sin(ϕ))^2 )
   I2_r_plasma = IPlasma * erf(distance/(σ*√2))

   r = hypot(x, y)
   k = √(4r*(a+b)/(z^2+((a+b)+r)^2))
   K, E = ellipke(k)
   Bz_plasma = μ₀*I2_r_plasma/(2π*√(z^2+((a+b)+r)^2))*(((a+b)^2-z^2-r^2)/(z^2+(r-(a+b))^2)*E+K)
   Br_plasma = μ₀*z*I2_r_plasma/(2π*r*√(z^2+(b+r)^2))*((z^2+r^2+(a+b)^2)/(z^2+(r-(a+b))^2)*E-K)
   Bx_plasma = Br_plasma*x/r
   By_plasma = Br_plasma*y/r

   if distance > 0.0001
      Bx += Bx_plasma
      By += By_plasma
      Bz += Bz_plasma
   end
   [Bx, By, Bz]
end


text:
#!/usr/bin/env juliaplain.dash

using ColorSchemes
for cs in ColorSchemes.colorschemes
    println(first(cs))
    for color in last(cs).colors
        println(color)
    end
end

max_stars_repo_name: NightMachinary/calendarheatmap

max_issues_repo_name: NightMachinary/calendarheatmap

max_forks_repo_name: NightMachinary/calendarheatmap


text:
fname = joinpath(dirname(@__FILE__), "..", "data", "test.sfp")

s = read_sfp(fname)

@test length(s)  == 3
@test label(s) == ["Fp1", "Fpz", "Fp2"]
@test EEG.x(s) == [-27.747648, -0.085967, 27.676888]
@test EEG.y(s) == [98.803864, 103.555275, 99.133354]
@test EEG.z(s) == [34.338360, 34.357265, 34.457005]
@test typeof(s[1].coordinate) == EEG.Talairach


text:
@doc raw"""
    rafterydiag(
        x::AbstractVector{<:Real}; q=0.025, r=0.005, s=0.95, eps=0.001, range=1:length(x)
    )

Compute the Raftery and Lewis diagnostic [^Raftery1992]. This diagnostic is used to
determine the number of iterations required to estimate a specified quantile `q` within a
desired degree of accuracy.  The diagnostic is designed to determine the number of
autocorrelated samples required to estimate a specified quantile $\theta_q$, such that
$\Pr(\theta \le \theta_q) = q$, within a desired degree of accuracy. In particular, if
$\hat{\theta}_q$ is the estimand and $\Pr(\theta \le \hat{\theta}_q) = \hat{P}_q$ the
estimated cumulative probability, then accuracy is specified in terms of `r` and `s`, where
$\Pr(q - r < \hat{P}_q < q + r) = s$. Thinning may be employed in the calculation of the
diagnostic to satisfy its underlying assumptions. However, users may not want to apply the
same (or any) thinning when estimating posterior summary statistics because doing so results
in a loss of information. Accordingly, sample sizes estimated by the diagnostic tend to be
conservative (too large).

Furthermore, the argument `r` specifies the margin of error for estimated cumulative
probabilities and `s` the probability for the margin of error. `eps` specifies the tolerance
within which the probabilities of transitioning from initial to retained iterations are
within the equilibrium probabilities for the chain. This argument determines the number of
samples to discard as a burn-in sequence and is typically left at its default value.

[^Raftery1992]: A L Raftery and S Lewis. Bayesian Statistics, chapter How Many Iterations in the Gibbs Sampler? Volume 4. Oxford University Press, New York, 1992.
"""
function rafterydiag(
    x::AbstractVector{<:Real}; q=0.025, r=0.005, s=0.95, eps=0.001, range=1:length(x)
)
    nx = length(x)
    phi = sqrt(2.0) * SpecialFunctions.erfinv(s)
    nmin = ceil(Int, q * (1.0 - q) * (phi / r)^2)
    if nmin > nx
        @warn "At least $nmin samples are needed for specified q, r, and s"
        kthin = -1
        burnin = total = NaN
    else
        dichot = Int[(x .<= StatsBase.quantile(x, q))...]
        kthin = 0
        bic = 1.0
        local test , ntest
        while bic >= 0.0
            kthin += 1
            test = dichot[1:kthin:nx]
            ntest = length(test)
            temp = test[1:(ntest - 2)] + 2 * test[2:(ntest - 1)] + 4 * test[3:ntest]
            trantest = reshape(StatsBase.counts(temp, 0:7), 2, 2, 2)
            g2 = 0.0
            for i1 in 1:2, i2 in 1:2, i3 in 1:2
                tt = trantest[i1, i2, i3]
                if tt > 0
                    fitted =
                        sum(trantest[:, i2, i3]) * sum(trantest[i1, i2, :]) /
                        sum(trantest[:, i2, :])
                    g2 += 2.0 * tt * log(tt / fitted)
                end
            end
            bic = g2 - 2.0 * log(ntest - 2.0)
        end

        tranfinal = StatsBase.counts(test[1:(ntest - 1)] + 2 * test[2:ntest], 0:3)
        alpha = tranfinal[3] / (tranfinal[1] + tranfinal[3])
        beta = tranfinal[2] / (tranfinal[2] + tranfinal[4])
        kthin *= step(range)
        m = log(eps * (alpha + beta) / max(alpha, beta)) / log(abs(1.0 - alpha - beta))
        burnin = kthin * ceil(m) + first(range) - 1
        n = ((2.0 - alpha - beta) * alpha * beta * phi^2) / (r^2 * (alpha + beta)^3)
        keep = kthin * ceil(n)
        total = burnin + keep
    end
    return (
        thinning=kthin, burnin=burnin, total=total, nmin=nmin, dependencefactor=total / nmin
    )
end


text:
using Distributions
using DistributionsAD
using Zygote
using Gen
struct DistributionsBacked{T} <: Gen.Distribution{T}
    to_dist
    has_arg_grads :: Vector{Bool}
    has_out_grad :: Bool
end

Gen.random(d::DistributionsBacked, args...) = rand(d.to_dist(args...))
Gen.logpdf(d::DistributionsBacked{T}, v::T, args...) where T  = Distributions.logpdf(d.to_dist(args...), v)
Gen.logpdf_grad(d::DistributionsBacked{T}, v::T,
            args...) where T = Zygote.gradient((v, args...) -> Distributions.logpdf(d.to_dist(args...), v), v, args...)
(d::DistributionsBacked)(args...) = Gen.random(d, args...)
Gen.has_output_grad(d::DistributionsBacked) = d.has_out_grad
Gen.has_argument_grads(d::DistributionsBacked) = d.has_arg_grads

max_stars_repo_name: jamesonquinn/Gen

max_issues_repo_name: jamesonquinn/Gen

max_forks_repo_name: jamesonquinn/Gen


text:
const GPUStructArray0{T, N} =
    StructArray{T, N, <: Tuple{Vararg{GPUArrays.AnyGPUArray}}}
const GPUStructArray1{T, N} =
    StructArray{T, N, <: NamedTuple{U, <: Tuple{Vararg{Union{GPUArrays.AnyGPUArray, GPUStructArray0}}}}} where {U}
const GPUStructArray{T, N} = Union{GPUStructArray0{T, N}, GPUStructArray1{T, N}}

# These two similar functions can be removed once
# <https://github.com/JuliaArrays/StructArrays.jl/pull/94>
# is accepted.
function Base.similar(s::StructArray{T,N,C}, ::Type{T}, sz::NTuple{M,Int64}) where {T,N,M,C<:Union{Tuple,NamedTuple}}
    return StructArray{T}(map(typ -> similar(typ, sz),
                              StructArrays.components(s)))
end

function Base.similar(s::StructArray{T,N,C}, S::Type, sz::NTuple{M,Int64}) where {T,N,M,C<:Union{Tuple,NamedTuple}}
    # If not specified, we don't really know what kind of array to use for each
    # interior type, so we just pick the first one arbitrarily. If users need
    # something else, they need to be more specific.
    f1 = StructArrays.components(s)[1]
    if isstructtype(S)
        return StructArrays.buildfromschema(typ -> similar(f1, typ, sz), S)
    else
        return similar(f1, S, sz)
    end
end

# The following broadcast code is slightly modified from the code found at
# <https://github.com/JuliaArrays/StructArrays.jl/issues/150>.
const GPUStore = Tuple{Vararg{GPUArrays.BroadcastGPUArray}}
const NamedGPUStore = NamedTuple{Name,<:GPUStore} where {Name}
const StructGPUArray = StructArray{T,N,<:Union{GPUStore,NamedGPUStore}} where {T,N}
## backend for StructArray
GPUArrays.backend(A::StructGPUArray) =
    GPUArrays.backend(StructArrays.components(A))
GPUArrays.backend(t::GPUStore) = GPUArrays.backend(typeof(t))
GPUArrays.backend(nt::NamedGPUStore) =
    GPUArrays.backend(typeof(nt).parameters[2])
function GPUArrays.backend(::Type{T}) where {T<:GPUStore}
    bs = GPUArrays.backend.(tuple(T.parameters...))
    I = all(map(isequal(first(bs)), bs))
    I || throw("device error")
    GPUArrays.backend(T.parameters[1])
end

## copy from GPUArrays
@inline function Base.copyto!(dest::StructGPUArray,
                              bc::Broadcast.Broadcasted{Nothing})
    axes(dest) == axes(bc) || Broadcast.throwdm(axes(dest), axes(bc))
    isempty(dest) && return dest
    bc′ = Broadcast.preprocess(dest, bc)

    # grid-stride kernel
    function broadcast_kernel(ctx, dest, bc′, nelem)
        for i in 1:nelem
            I = GPUArrays.@cartesianidx(dest, i)
            @inbounds dest[I] = bc′[I]
        end
        return
    end
    heuristic = GPUArrays.launch_heuristic(GPUArrays.backend(dest),
                                           broadcast_kernel, dest, bc′, 1)
    config = GPUArrays.launch_configuration(GPUArrays.backend(dest),
                                            heuristic, length(dest),
                                            typemax(Int))
    GPUArrays.gpu_call(broadcast_kernel, dest, bc′, config.elements_per_thread;
                       threads=config.threads, blocks=config.blocks)

    return dest
end

function Base.similar(bc::Broadcast.Broadcasted{StructArrays.StructArrayStyle{S}},
                      ::Type{T}) where {S<:CUDA.CuArrayStyle,T}
    if isstructtype(T)
        return StructArrays.buildfromschema(typ -> similar(CuArray{typ},
                                                           axes(bc)), T)
    else
        return similar(CuArray{T}, axes(bc))
    end
end

# We follow GPUArrays approach of coping the whole array to the host when
# outputting a StructArray backed by GPU arrays.
convert_to_cpu(xs) = adapt(Array, xs)
function Base.print_array(io::IO, X::GPUStructArray{<:Any,0})
    X = convert_to_cpu(X)
    isassigned(X) ? show(io, X[]) : print(io, undef_ref_str)
end
Base.print_array(io::IO, X::GPUStructArray{<:Any,1}) =
    Base.print_matrix(io, convert_to_cpu(X))
Base.print_array(io::IO, X::GPUStructArray{<:Any,2}) where {T} =
    Base.print_matrix(io, convert_to_cpu(X))
Base.print_array(io::IO, X::GPUStructArray{<:Any,<:Any}) =
    Base.show_nd(io, convert_to_cpu(X), Base.print_matrix, true)

@kernel function fill_kernel!(A, x)
    I = @index(Global)
    @inbounds A[I] = x
end

function Base.fill!(A::GPUStructArray, x)
    event = Event(device(A))
    event = fill_kernel!(device(A), 256)(A, x, ndrange = length(A),
                                         dependencies = (event, ))
    wait(event)
end

function device(s::StructArray)
    ds = map(device, StructArrays.components(s))
    I = all(map(isequal(first(ds)), ds))
    I || throw("device error")
    return first(ds)
end

function Tullio.storage_type(S::StructArray)
    return Tullio.storage_type(StructArrays.components(S)...)
end

components(S::StructArray) = StructArrays.components(S)


text:
#=
recorder.jl:

Author: Peter Wolf (opus111@gmail.com)

Port of the FastAI Recorder to Julia

This code is inspired by FastAI, but differs from it in important ways

Recorders are Callbacks added to a Learner.  They keep a log of
statistics (losses and smooth_losses) and history (opts) during training

This design is significantly different from the original
Recorder and Metrics objects in Python

https://github.com/fastai/fastai2/blob/master/fastai2/learner.py

The documentation is copied from here

https://dev.fast.ai/learner#Recorder
=#
module Recorder
using FastAI
using FastAI: AbstractCallback,AbstractLearner
using Infiltrator
"""
    Recorder(learn::Learner; train_loss = true, train_smooth_loss = true,
                             validate_loss = true, validate_smooth_loss = true)

Container for [`Learner`](@ref) statistics (e.g. lr, loss and metrics) during training.
Statistics are indexed by name, epoch and batch.
For example to get the smoothed training loss for epoch 2, batch 3, we would call
```julia
recorder["TrainSmoothLoss", 2, 3]
```
To get the entire history of smooth training loss, one would call
```julia
recorder["TrainSmoothLoss", :, :]
```
"""
mutable struct Smoother
    alpha::Real
    val::Real
    first::Bool
end
Smoother(alpha) = Smoother(alpha, 0.0, true)

reset!(asl::Smoother) = asl.first=true

function (asl::Smoother)(value)
    if asl.first
        asl.first = false
        asl.val = value
    else
        asl.val = asl.alpha*asl.val+(1-asl.alpha)*value
    end
    return asl.val
end

"""
AbstractRecorder is shared code

All Recorders have a log value, where the log supports getindex and setindex
"""
abstract type AbstractRecorder <: AbstractCallback end
FastAI.before_fit(rec::AbstractRecorder,lrn::AbstractLearner, epoch_count, batch_size) = rec.log = fill(NaN,(epoch_count,batch_size))
Base.getindex(rec::AbstractRecorder,idx...) = rec.log[idx...]

"""
    TrainLoss

Record a log of training loss
"""
mutable struct TrainLoss <: AbstractRecorder
    log::Union{Nothing,Array{Real,2}}
end
TrainLoss()=TrainLoss(nothing)
FastAI.batch_train_loss(rec::TrainLoss,lrn::AbstractLearner, epoch, batch, loss) = rec.log[epoch,batch] = loss

"""
    ValidateLoss

Record a log of validation loss
"""
mutable struct ValidateLoss <: AbstractRecorder
    log::Union{Nothing,Array{Real,2}}
end
ValidateLoss()=ValidateLoss(nothing)
FastAI.batch_validate_loss(rec::ValidateLoss,lrn::AbstractLearner, epoch, batch, loss) = rec.log[epoch,batch] = loss

abstract type AbstractSmoothRecorder <: AbstractRecorder end

function FastAI.before_fit(rec::AbstractSmoothRecorder,lrn::AbstractLearner, epoch_count, batch_size)
     rec.log = fill(NaN,(epoch_count,batch_size))
     reset!(rec.smooth)
end

"""
    SmoothTrainLoss

Record a smoothed log of training loss
"""
mutable struct SmoothTrainLoss <: AbstractSmoothRecorder
    smooth::Smoother
    log::Union{Nothing,Array{Real,2}}
end
SmoothTrainLoss(alpha=0.98)=SmoothTrainLoss(Smoother(alpha),nothing)
FastAI.batch_train_loss(rec::SmoothTrainLoss,lrn::AbstractLearner, epoch, batch, loss) = rec.log[epoch,batch] = rec.smooth(loss)

"""
    SmoothValidateLoss

Record a smoothed log of validation loss
"""
mutable struct SmoothValidateLoss <: AbstractSmoothRecorder
    smooth::Smoother
    log::Union{Nothing,Array{Real,2}}
end
SmoothValidateLoss(alpha=0.98)=SmoothValidateLoss(Smoother(alpha),nothing)
FastAI.batch_validate_loss(rec::SmoothValidateLoss,lrn::AbstractLearner, epoch, batch, loss) = rec.log[epoch,batch] = rec.smooth(loss)

end


text:
#############################################################
#
# Direct Messages Functions
#
#############################################################

function get_direct_messages(; options=Dict{AbstractString, AbstractString}())

    r = get_oauth("https://api.twitter.com/1.1/direct_messages.json", options)

    #If 200 response, create Array{TWEETS,1}, otherwise return raw response
    return r.status == 200 ? to_TWEETS(Requests.json(r)) : r

end

function get_direct_messages_sent(; options=Dict{AbstractString, AbstractString}())

    r = get_oauth("https://api.twitter.com/1.1/direct_messages/sent.json", options)

    #If 200 response, create Array{TWEETS,1}, otherwise return raw response
    return r.status == 200 ? to_TWEETS(Requests.json(r)) : r

end

function get_direct_messages_show(id::AbstractString; options=Dict{AbstractString, AbstractString}())

    #Add required parameter(s) to options dict
    options["id"] = id

    r = get_oauth("https://api.twitter.com/1.1/direct_messages/show.json", options)

    #If 200 response, create Array{TWEETS,1}, otherwise return raw response
    return r.status == 200 ? to_TWEETS(Requests.json(r)) : r

end

function post_direct_messages_destroy(id::AbstractString; options=Dict{AbstractString, AbstractString}())

    options["id"] = id

    r = post_oauth("https://api.twitter.com/1.1/direct_messages/destroy.json", options)

    #If 200 response, create Array{TWEETS,1}, otherwise return raw response
    return r.status == 200 ? to_TWEETS(Requests.json(r)) : r

end

function post_direct_messages_send(text::AbstractString; options=Dict{AbstractString, AbstractString}())

    #Add status into options Dict
    options["text"] = text

    r = post_oauth("https://api.twitter.com/1.1/direct_messages/new.json", options)

    #If 200 response, create Array{TWEETS,1}, otherwise return raw response
    return r.status == 200 ? to_TWEETS(Requests.json(r)) : r

end


text:
using Zygote, Test
using Zygote: pullback, @adjoint

macro test_inferred(ex)
  :(let res = nothing
    @test begin
      res = @inferred $ex
      true
    end
    res
  end) |> esc
end

trace_contains(st, func, file, line) = any(st) do fr
  func in (nothing, fr.func) && endswith(String(fr.file), file) &&
    fr.line == line
end

bad(x) = x
@adjoint bad(x) = x, Δ -> error("bad")

function badly(x)
  x = x + 1
  x = bad(x)
  return x
end

y, back = pullback(badly, 2)
@test y == 3
@test_throws Exception back(1)
bt = try back(1) catch e stacktrace(catch_backtrace()) end

@test trace_contains(bt, nothing, "compiler.jl", 20)
@test trace_contains(bt, :badly, "compiler.jl", 24)

# Type inference checks

Zygote.refresh()

y, back = @test_inferred pullback(*, 2, 3)
@test_inferred(back(1))

_sincos(x) = sin(cos(x))

y, back = @test_inferred pullback(_sincos, 0.5)
@test_inferred back(1)

f(x) = 3x^2 + 2x + 1

y, back = @test_inferred pullback(f, 5)
@test y == 86
@test_inferred(back(1))

y, back = @test_inferred pullback(Core._apply, +, (1, 2, 3))
@test_inferred back(1)

# TODO fix bcast inference
# bcast(x) = x .* 5
# y, back = @test_inferred pullback(bcast, [1,2,3])
# @test_inferred back([1,1,1])

foo = let a = 4
  x -> x*a
end

@test_inferred gradient(f -> f(5), foo)

getx(x) = x.x
y, back = @test_inferred pullback(getx, (x=1,y=2.0))
@test_inferred back(1)

y, back = @test_inferred pullback(x->x[1], (5,:a))
@test_inferred back(1)

y, back = @test_inferred pullback(((a,b),) -> a, (5, 10))
@test_inferred back(1)


text:
function sciml_train(loss, θ, opt, adtype::DiffEqBase.AbstractADType = GalacticOptim.AutoZygote(), args...; kwargs...)
  optf = GalacticOptim.OptimizationFunction((x, p) -> loss(x), adtype)
  optfunc = GalacticOptim.instantiate_function(optf, θ, adtype, nothing)
  optprob = GalacticOptim.OptimizationProblem(optfunc, θ; kwargs...)
  GalacticOptim.solve(optprob, opt, args...; kwargs...)
end


text:
@testset "Canonisation and Compression" begin

D = 32
Dcut = 16

d = 2
sites = 100

T = Float64

var_tol = 1E-10
var_max_sweeps = 100

ψ = randn(MPS{T}, sites, D, d)
ϕ = randn(MPS{T}, sites, D, d)
χ = randn(MPS{T}, sites, D, d)
Φ = randn(MPS{T}, sites, D, d)

@testset "Canonisation (left)" begin
    b = canonise!(ψ, :left)
    @test is_left_normalized(ψ)
    @test dot(ψ, ψ) ≈ 1
end

@testset "Canonisation (right)" begin
    b = canonise!(ϕ, :right)
    @test is_right_normalized(ϕ)
    @test dot(ϕ, ϕ) ≈ 1
end

@testset "Copy and truncate twice" begin
    ψ̃ = copy(ψ)
    @test ψ̃ == ψ
    for (direction, predicate) ∈ ((:left, is_left_normalized), (:right, is_right_normalized))
        truncate!(ψ, direction, Dcut)
        truncate!(ψ̃, direction, Dcut)

        @test predicate(ψ)
        @test predicate(ψ̃)
        @test bond_dimension(ψ̃) == bond_dimension(ψ)
        @test all(size(A) == size(B) for (A, B) ∈ zip(ψ, ψ̃))
        @test typeof(ψ̃) == typeof(ψ)
        @test norm(ψ) ≈ norm(ψ̃) ≈ 1
        @test abs(ψ * ψ̃) ≈ abs(ψ̃ * ψ) ≈ 1
    end
end

@testset "Cauchy-Schwarz inequality (after truncation)" begin
    @test abs(dot(ϕ, ψ)) <= norm(ϕ) * norm(ψ)
end

@testset "Truncation (SVD, right)" begin
    truncate!(ψ, :right, Dcut)
    @test is_right_normalized(ψ)
    @test norm(ψ) ≈ 1
end

@testset "Truncation (SVD, left)" begin
    truncate!(ψ, :left, Dcut)
    @test is_left_normalized(ψ)
    @test norm(ψ) ≈ 1
end

@testset "<left|right>" begin
    ϵ = 1E-10
    ψ = randn(MPS{T}, sites, D, d)

    l = copy(ψ)
    r = copy(ψ)

    canonise!(l, :left)
    @test is_left_normalized(l)

    canonise!(r, :right)
    @test is_right_normalized(r)

    @test dot(l, l) ≈ 1
    @test dot(r, r) ≈ 1

    @test abs(1 - abs(dot(l, r))) < ϵ
end

@testset "Variational compression" begin
    Ψ = copy(Φ)
    canonise!(Ψ, :left)

    overlap = compress!(Φ, Dcut, var_tol, var_max_sweeps)
    #println(overlap)

    @test norm(Φ) ≈ 1
    @test is_left_normalized(Φ)
    @test is_right_normalized(Φ) == false
end

end


text:
using Openscad
using Test

@testset "2d shapes" begin
    @test (toString(Circle(10, false)) == "circle(r = 10, center = false);")
    @test (toString(Square([1, 1], true)) == "square(size = [1, 1], center = true);")
end

function test2d()
    square = Square([4, 5], true)
    square2 = Square([0, 0], true)
    tree =
    Unify([
        Circle(3, false),
        square,
        Translate([1, 1, 1], [
            Openscad.Unify([
                square,
                square2
            ])
        ]),
        End()
    ])
    tree
end

@testset "complex 2d shape" begin
    @test (toString(test2d()) == "union(){circle(r = 3, center = false);square(size = [4, 5], center = true);translate([1, 1, 1]) {union(){square(size = [4, 5], center = true);square(size = [0, 0], center = true);};};};")
end

@testset "3d shapes" begin
end


text:
make_counters(nthreads = Threads.nthreads()) =
    RecordArrays.fill((n = 0,), nthreads; align = 64)

mutable struct ThreadLocalCounter
    name::Symbol
    _module::Module
    line::LineNumberNode
    counters::typeof(make_counters(1))
end

ThreadLocalCounter(name::Symbol, m::Module, line::LineNumberNode) =
    ThreadLocalCounter(name, m, line, make_counters())

function init!(tlc::ThreadLocalCounter)
    tlc.counters = make_counters()
    return tlc
end

function inc!(tlc::ThreadLocalCounter)
    tlc.counters.n[Threads.threadid()] += 1
end
# TODO: put `@inbounds`, once `init!` is tested well?

function Base.empty!(tlc::ThreadLocalCounter)
    fill!(tlc.counters.n, 0)
    return tlc
end

const TLC_BUCKET = StaticStorages.BucketKey()

"""
    @tlc [name]

Count the time this expression is evaluated using thread-local counters.

# Examples
```jldoctest
julia> using ThreadLocalCounters

julia> hello_world() = @tlc hello_world;
```
"""
macro tlc(name::Symbol = :_default_)
    tlc = ThreadLocalCounter(name, __module__, __source__)
    StaticStorages.put!(__module__, TLC_BUCKET, tlc)

    @gensym initmodule
    initializer = quote
        module $initmodule
        const TLC = $(QuoteNode(tlc))
        __init__() = $init!(TLC)
        end
    end
    @assert initializer.head === :block
    initializer.head = :toplevel
    Base.eval(__module__, initializer)

    quote
        inc!($(esc(initmodule)).TLC)
        nothing
    end
end

Base.sum(tlc::ThreadLocalCounter) = sum(x.n for x in tlc.counters)

"""
    ThreadLocalCounters.list(; all = false)

Get a list of all thread-local counters.  The caller must ensure that no thread is accessing
the counters.

Pass `all = true` to include counters with zero counts.
"""
ThreadLocalCounters.list

function ThreadLocalCounters.list(; all = false)
    bucket = StaticStorages.getbucket(TLC_BUCKET)
    bucket === nothing && return ThreadLocalCounter[]
    cs = collect(ThreadLocalCounter, values(bucket))
    if !all
        filter!(tlc -> sum(tlc) > 0, cs)
    end
    return cs
end

"""
    ThreadLocalCounters.clear()

Reset all counters. The caller must ensure that no thread is accessing the
counters.
"""
function ThreadLocalCounters.clear()
    bucket = StaticStorages.getbucket(TLC_BUCKET)
    bucket === nothing && return nothing
    foreach(empty!, values(bucket))
end

function Base.show(io::IO, ::MIME"text/plain", tlc::ThreadLocalCounter)
    if get(io, :typeinfo, Any) !== ThreadLocalCounter
        print(io, "ThreadLocalCounter: ")
    end
    show(io, [x.n for x in tlc.counters])
    if tlc.name !== :_default_
        print(io, ' ')
        printstyled(io, tlc.name; color = :blue)
    end
    print(io, " @", tlc._module, ' ')
    print(io, tlc.line, ' ')
end


text:
module simplexpop
# These have no dependencies.
include("sort.jl")
include("hilbert.jl")

# Manipulate rasters.
include("assignweights.jl")

# Finite elements.
include("decplex.jl")
include("fibgrid.jl")
include("hexlattice.jl")

# From sort
export countedinsertionsort!

# From hilbert
export encode_hilbert_zero, decode_hilbert_zero
export encode_hilbert, decode_hilbert, hilbert_order

# From assignweights.jl
export PixelGrid, pixelgrid, DataGrid
export assignweights, intersection_area
export geo_to_transform, geo_to_xy_corner, geo_to_xy_center
export xy_bounds, pixel_containing, ij_cover_rect, block_cover_ij_rect
export pixels_of_block, corner_pixel, pixelgrid, load_single_block
export load_pixel_grid

end # module


text:
const sto3g = Any[
  # H
  [('S',
    [(3.4252509099999999, 0.15432897000000001),
     (0.62391373000000006, 0.53532813999999995),
     (0.16885539999999999, 0.44463454000000002)])],
  # He
  [('S',
    [(6.3624213899999997, 0.15432897000000001),
     (1.1589229999999999, 0.53532813999999995),
     (0.31364978999999998, 0.44463454000000002)])],
  # Li
  [('S',
    [(16.119575000000001, 0.15432897000000001),
     (2.9362007000000001, 0.53532813999999995),
     (0.79465050000000004, 0.44463454000000002)]),
   ('S',
    [(0.63628969999999996, -0.099967230000000004),
     (0.14786009999999999, 0.39951282999999999),
     (0.048088699999999998, 0.70011546999999996)]),
   ('P',
    [(0.63628969999999996, 0.15591627),
     (0.14786009999999999, 0.60768372000000004),
     (0.048088699999999998, 0.39195739000000002)])],
  # Be
  [('S',
    [(30.167871000000002, 0.15432897000000001),
     (5.4951153000000001, 0.53532813999999995),
     (1.4871927, 0.44463454000000002)]),
   ('S',
    [(1.3148331, -0.099967230000000004),
     (0.3055389, 0.39951282999999999),
     (0.099370700000000006, 0.70011546999999996)]),
   ('P',
    [(1.3148331, 0.15591627),
     (0.3055389, 0.60768372000000004),
     (0.099370700000000006, 0.39195739000000002)])],
  # B
  [('S',
    [(48.791113000000003, 0.15432897000000001),
     (8.8873622000000001, 0.53532813999999995),
     (2.4052669999999998, 0.44463454000000002)]),
   ('S',
    [(2.2369561, -0.099967230000000004),
     (0.51982050000000002, 0.39951282999999999),
     (0.16906180000000001, 0.70011546999999996)]),
   ('P',
    [(2.2369561, 0.15591627),
     (0.51982050000000002, 0.60768372000000004),
     (0.16906180000000001, 0.39195739000000002)])],
  # C
  [('S',
    [(71.616837000000004, 0.15432897000000001),
     (13.045095999999999, 0.53532813999999995),
     (3.5305122, 0.44463454000000002)]),
   ('S',
    [(2.9412493999999998, -0.099967230000000004),
     (0.68348310000000001, 0.39951282999999999),
     (0.22228990000000001, 0.70011546999999996)]),
   ('P',
    [(2.9412493999999998, 0.15591627),
     (0.68348310000000001, 0.60768372000000004),
     (0.22228990000000001, 0.39195739000000002)])],
  # N
  [('S',
    [(99.106168999999994, 0.15432897000000001),
     (18.052312000000001, 0.53532813999999995),
     (4.8856602000000002, 0.44463454000000002)]),
   ('S',
    [(3.7804559000000002, -0.099967230000000004),
     (0.87849659999999996, 0.39951282999999999),
     (0.28571439999999998, 0.70011546999999996)]),
   ('P',
    [(3.7804559000000002, 0.15591627),
     (0.87849659999999996, 0.60768372000000004),
     (0.28571439999999998, 0.39195739000000002)])],
  # O
  [('S',
    [(130.70931999999999, 0.15432897000000001),
     (23.808861, 0.53532813999999995),
     (6.4436083000000002, 0.44463454000000002)]),
   ('S',
    [(5.0331513000000001, -0.099967230000000004),
     (1.1695960999999999, 0.39951282999999999),
     (0.38038899999999998, 0.70011546999999996)]),
   ('P',
    [(5.0331513000000001, 0.15591627),
     (1.1695960999999999, 0.60768372000000004),
     (0.38038899999999998, 0.39195739000000002)])],
  # F
  [('S',
    [(166.67912999999999, 0.15432897000000001),
     (30.360811999999999, 0.53532813999999995),
     (8.2168206999999995, 0.44463454000000002)]),
   ('S',
    [(6.4648032000000004, -0.099967230000000004),
     (1.5022812000000001, 0.39951282999999999),
     (0.48858849999999998, 0.70011546999999996)]),
   ('P',
    [(6.4648032000000004, 0.15591627),
     (1.5022812000000001, 0.60768372000000004),
     (0.48858849999999998, 0.39195739000000002)])],
  # Ne
  [('S',
     [(207.01561000000001, 0.15432897000000001),
    (37.708151000000001, 0.53532813999999995),
    (10.205297, 0.44463454000000002)]),
    ('S',
     [(8.2463151000000003, -0.099967230000000004),
    (1.9162661999999999, 0.39951282999999999),
    (0.62322929999999999, 0.70011546999999996)]),
    ('P',
     [(8.2463151000000003, 0.15591627),
    (1.9162661999999999, 0.60768372000000004),
      (0.62322929999999999, 0.39195739000000002)])]
]
const basis_set_data = Dict{Any,Any}("sto3g" => sto3g);


text:
using Gtk
using Gtk.GConstants

@field GtkParent = Gtk.Window("")

@composed type GtkWindow <: AbstractGLWindow
    <: GLWindow
    GtkParent
end

function default(::Type{GtkParent}, data::Partial{GtkWindow})
    area = get!(data, Area)
    resolution = round.(Int, widths(get(data, Area)))
    parent = Gtk.Window(string(get(data, Name)), resolution..., true, true)
    Gtk.visible(parent, get(data, Visible))
    Gtk.setproperty!(parent, Symbol("is-focus"), get(data, Focused))
    parent
end

function default(::Type{NativeWindow}, data::Partial{GtkWindow})
    parent = get!(data, GtkParent)
    window = Gtk.GLArea()
    opengl_version = get!(data, OpenGLVersion)
    Gtk.gl_area_set_required_version(window, opengl_version.major, opengl_version.minor)
    GLAbstraction.new_context(window)
    push!(parent, window)
    window
end

function add!(window::GtkWindow, ::Type{Area})
    function callback_w(widget::Gtk.GtkGLArea, width::Int32, height::Int32)
        rect = window[Area]
        window[Area] = IRect(minimum(rect), width, height)
        return true
    end
    signal_connect(callback_w, window[Window], "Scroll")
    add_events(window[Window], GConstants.GdkEventMask.SCROLL)
    return
end

function add!(window::GtkWindow, ::Type{Mouse.Scroll})
    function callback(widget::Gtk.GtkGLArea, s::Gtk.GdkEventScroll)
        window[Scroll] = (s.x, s.y)
        return true
    end
    signal_connect(callback, window[NativeWindow], "scroll-event")
    add_events(window[NativeWindow], GConstants.GdkEventMask.SCROLL)
    return
end

function add!(window::GtkWindow, ::Type{Mouse.Position})
    function callback(widget::Gtk.GtkGLArea, s::Gtk.GdkEventMotion)
        window[Mouse.Position] = (s.x, s.y)
        return true
    end
    add_events(window[NativeWindow], GConstants.GdkEventMask.POINTER_MOTION)
    signal_connect(callback, window[NativeWindow], "motion-notify-event")
    return true
end

function to_mouse_button(x)
    if x == 1
        Mouse.left
    elseif x == 2
        Mouse.middle
    elseif x == 3
        Mouse.right
    else
        # TODO turn into error
        warn("Button is $x, while $(Gtk.GdkModifierType.BUTTON1)")
        Mouse.left
    end
end
function add!(window::GtkWindow, ::Type{Mouse.Buttons})
    function callback(widget::Gtk.GtkGLArea, event::Gtk.GdkEventButton)
        button = to_mouse_button(event.button)
        action = event.event_type
        set = window[Mouse.Buttons]
        if action in (GdkEventType.BUTTON_PRESS:GdkEventType.TRIPLE_BUTTON_PRESS)
            push!(set, button)
        elseif action == GdkEventType.BUTTON_RELEASE
            delete!(set, button)
        else
            warn("unknown action: $(action)")
        end
        window[Mouse.Buttons] = set # trigger setfield event!
        return true
    end
    add_events(window[NativeWindow],
        GConstants.GdkEventMask.GDK_BUTTON_PRESS_MASK |
        GConstants.GdkEventMask.GDK_BUTTON_RELEASE_MASK
    )
    signal_connect(callback, window[NativeWindow], "button_press_event")
    signal_connect(callback, window[NativeWindow], "button_release_event")
    return
end
# gint glarea_key_press_event(GtkWidget *widget, GdkEventKey *event)
# {
#   switch (event->keyval) {
#     case GDK_r:
#       g_print("Button r pressed...redrawing\n");
#       gtk_widget_draw(glarea, (GdkRectangle *)NULL);
#       break;
#     case GDK_l:
#       g_print("Button l pressed...redrawing\n");
#       gtk_widget_draw(glarea, (GdkRectangle *)NULL);
#       break;
#     case GDK_p:
#       g_print("Button p pressed...redrawing\n");
#       gtk_widget_draw(glarea, (GdkRectangle *)NULL);
#       break;
#   }
#   return (TRUE);
# }

function make_context_current(x::Gtk.GLArea)
    Gtk.make_current(x)
    GLAbstraction.make_context_current(x)
end

function gtk_render(gl_area, gdk_context)
    signal_connect(gtk_render, area, "render")
    global draw_cat, uniforms, args, vbo
    glEnable(GL_DEPTH_TEST)
    glDepthMask(GL_TRUE)
    glDepthFunc(GL_LEQUAL)
    glDisable(GL_CULL_FACE)
    GLAbstraction.enabletransparency()
    glClearColor(1,1,1,1)
    glViewport(0, 0, 200, 200)
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    args = (Vec3f0(1, 0, 0), cam[Visualize.Projection], cam[View])
    uniforms[3][1] = args # update
    draw_cat(vbo, uniforms)
    glFlush()
    return false
end


text:
function to_string(ex::SymbolicExpression)
# s expression style is easy, but not very nice
    if is_binary_operator(ex.head)
        __to_string_binary_op(ex)
    elseif is_unary_operator(ex.head)
        __to_string_unary_op(ex)
    else
        __to_string_k_ary_expression(ex)
    end
end

function to_string(s::String)
    return s
end

function __to_string_minus_op(ex::SymbolicExpression)
    if is_variable(ex.args[1])
        rhs = ex.args[1].name
    else
        rhs = to_string(ex.args[1])
    end

    do_paren = true
    if is_variable(ex.args[1]) || is_constant(ex.args[1])
        do_paren = false
    end

    left_paren  = "("
    right_paren = ")"

    if do_paren
        total_str = "-" * left_paren * rhs * right_paren
    else
        total_str = "-" * rhs
    end

    return total_str
end

function __to_string_unary_op(ex::SymbolicExpression)
    if is_minus_operator(ex.head)
        return __to_string_minus_op(ex)
    end
    arg_str   = to_string(ex.args[1])
    total_str = string(ex.head) * "(" * arg_str * ")"
    return total_str
end

function __to_string_k_ary_expression(ex)
    op = string(ex.head)

end

function __to_string_binary_op(ex::SymbolicExpression)

    L = length(ex.args)
    if L < 2
       return __to_string_unary_op(ex)
    end

    # op  = string(ex.head)
    # lhs = ex.args[1]
    # rhs = ex.args[2]

    # c op c
    # var op var
    # comp op comp -> hard case
    #
    # c op var | var op c
    #
    # c op comp

    # rule: unary minus

    # rule: const op const
    if is_constant(ex.args[1]) && is_constant(ex.args[2])
        return __to_string_const_op_const(ex)
    end

    # rule: var op var
    if is_variable(ex.args[1]) && is_variable(ex.args[2])
        return __to_string_var_op_var(ex)
    end

    # rule: c op var | var op c
    if (is_variable(ex.args[1]) || is_constant(ex.args[1])) &&
        (is_variable(ex.args[2]) || is_constant(ex.args[2]))
        return __to_string_var_op_var(ex)
    end

    # c op comp
    if (is_constant(ex.args[1]) || is_variable(ex.args[1])) && is_composite_expression(ex.args[2])
        return __to_string_const_op_comp(ex, true)
    end
    if (is_constant(ex.args[2]) || is_variable(ex.args[2])) && is_composite_expression(ex.args[1])
        return __to_string_const_op_comp(ex, false)
    end

    # rule: constant */^ var or comp both way
    if (is_product_operator(ex.head) || is_divide_operator(ex.head) || is_power_operator(ex.head)) &&
        ( is_constant(ex.args[1]) || is_constant(ex.args[2]) ||
          is_variable(ex.args[1]) || is_variable(ex.args[2]) )
        return __to_string_constant_op_composite_expression(ex)
    end

    # rule: composite + composite
    if is_plus_operator(ex.head) & is_composite_expression(ex.args[1]) & is_composite_expression(ex.args[2])
        return __to_string_comp_plus_comp(ex)
    end

    # rule: composite - composite
    if is_minus_operator(ex.head) & is_composite_expression(ex.args[1]) & is_composite_expression(ex.args[2])
        return __to_string_comp_minus_comp(ex)
    end

    # rule: composite * composite
    if is_product_operator(ex.head) & is_composite_expression(ex.args[1]) & is_composite_expression(ex.args[2])
        return __to_string_comp_prod_comp(ex)
    end

    return __to_string_fallback(ex)
end

function __to_string_fallback(ex)

    op        = string(ex.head)
    lhs       = ex.args[1]
    rhs       = ex.args[2]
    lhs_str   = to_string(lhs)
    rhs_str   = to_string(rhs)
    ws        = " "

    left_paren   = "("
    right_paren  = ")"

    return left_paren * lhs_str * right_paren * ws * op * ws * left_paren * rhs_str * right_paren
end

function __to_string_const_op_const(ex)
    op        = string(ex.head)
    ws        = " "

    if is_plus_operator(ex.head) || is_product_operator(ex.head)
        return __to_string_plus_or_prod(ex)
    end

    lhs       = ex.args[1]
    rhs       = ex.args[2]
    lhs_str   = to_string(lhs)
    rhs_str   = to_string(rhs)
    total_str = lhs_str * ws * op * ws * rhs_str

    return total_str
end

function __to_string_var_op_var(ex::SymbolicExpression)

    if is_plus_operator(ex.head) || is_product_operator(ex.head)
        return __to_string_plus_or_prod(ex)
    end

    op        = string(ex.head)
    lhs       = ex.args[1]
    rhs       = ex.args[2]
    lhs_str   = to_string(lhs)
    rhs_str   = to_string(rhs)
    ws        = is_power_operator(ex.head) ? "" : " "
    total_str = lhs_str * ws * op * ws * rhs_str

    return total_str
end

function __to_string_const_op_comp(ex, lhs_is_const::Bool)

    if is_plus_operator(ex.head) || is_product_operator(ex.head)
        return __to_string_plus_or_prod(ex)
    end

    op        = ex.head
    op_str    = string(op)
    lhs       = ex.args[1]
    rhs       = ex.args[2]
    lhs_str   = to_string(lhs)
    rhs_str   = to_string(rhs)
    ws        = " "
    left_paren   = "("
    right_paren  = ")"

    if is_plus_operator(op)
        # c + ex
        return lhs_str * ws * op_str * ws * rhs_str
    end
    if is_minus_operator(op)
        # c - ex
        # c - (a +- b)
        # paren if the next op is +-

        if lhs_is_const
            if is_plus_operator(rhs.head) || is_minus_operator(rhs.head)
                return lhs_str * ws * op_str * ws * left_paren * rhs_str * right_paren
            else
                return left_paren * lhs_str * right_paren * ws * op_str * ws * rhs_str
            end
        else
            return lhs_str * ws * op_str * ws * rhs_str
        end
    end

    rhs_is_const = !lhs_is_const

    if lhs_is_const
        if is_product_operator(op) && (is_product_operator(rhs.head) || is_divide_operator(rhs.head))
            # dont paren if the next op is also * or /
            return lhs_str * ws * op_str * ws * rhs_str
        else
            return lhs_str * ws * op_str * ws * left_paren * rhs_str * right_paren
        end
    end
    if rhs_is_const
        if is_product_operator(op) && (is_product_operator(lhs.head) || is_divide_operator(lhs.head))
            # dont paren if the next op is also * or /
            return lhs_str * ws * op_str * ws * rhs_str
        else
            return left_paren * lhs_str * right_paren * ws * op_str * ws * rhs_str
        end
    end

end

function __to_string_constant_op_composite_expression(ex::SymbolicExpression)

    if is_plus_operator(ex.head) || is_product_operator(ex.head)
        return __to_string_plus_or_prod(ex)
    end

    op  = string(ex.head)
    lhs = ex.args[1]
    rhs = ex.args[2]

    if is_constant(lhs) || is_variable(lhs)
        return __to_string_constant_op_composite_expression_constvarleft(ex)
    end

    if is_constant(rhs) || is_variable(rhs)
        return __to_string_constant_op_composite_expression_constvarright(ex)
    end
end

function __to_string_constant_op_composite_expression_constvarleft(ex)

    # a * (b +- c)
    # a / (b +- c)
    # a ^ (b +-*/ c)
    #
    # a * b */ c
    # a / b */ c

    op        = ex.head
    op_str    = string(op)
    lhs       = ex.args[1]
    rhs       = ex.args[2]
    lhs_str   = to_string(lhs)
    rhs_str   = to_string(rhs)
    ws        = " "
    left_paren   = "("
    right_paren  = ")"

    b1 = is_plus_operator(rhs.head) || is_minus_operator(rhs.head)

    do_paren = false
    if (is_product_operator(op) ||  is_divide_operator(op)) && b1
        do_paren = true
    end

    omit_ws = false
    if is_power_operator(op)
        do_paren = true
        omit_ws  = true
    end

    ws          = omit_ws ? "" : ws
    left_paren  = do_paren ? left_paren : ""
    right_paren = do_paren ? right_paren : ""

    return lhs_str * ws * op_str * ws * left_paren * rhs_str * right_paren
end

function __to_string_constant_op_composite_expression_constvarright(ex)

    # (a +- b) * c
    # (a +- b) / c
    # (a +- b) ^ c
    #
    # a * b */ c

    op        = ex.head
    op_str    = string(op)
    lhs       = ex.args[1]
    rhs       = ex.args[2]
    lhs_str   = to_string(lhs)
    rhs_str   = to_string(rhs)
    ws        = " "
    left_paren   = "("
    right_paren  = ")"

    b1 = is_plus_operator(lhs.head) || is_minus_operator(lhs.head)

    do_paren = false
    if (is_product_operator(op) ||  is_divide_operator(op)) && b1
        do_paren = true
    end

    omit_ws = false
    if is_power_operator(op)
        do_paren = true
        omit_ws  = true
    end

    ws          = omit_ws ? "" : ws
    left_paren  = do_paren ? left_paren : ""
    right_paren = do_paren ? right_paren : ""

    return left_paren * lhs_str * right_paren * ws * op_str * ws * rhs_str
end

function __to_string_plus_or_prod(ex)

    op        = string(ex.head)
    strs      = map(to_string, ex.args)
    ws        = " "
    L         = length(ex.args)
    Ltotal    = L + L-1
    total_str = Vector{String}(undef, Ltotal)

    for ii = 1:L-1
        total_str[2*ii-1] = strs[ii]
        total_str[2*ii]   = ws * op * ws
    end

    total_str[Ltotal] = strs[L]
    reduced_total_str = reduce(*, total_str)

    return reduced_total_str
end

function __to_string_comp_plus_comp(ex)
    # op        = string(ex.head)
    # lhs       = ex.args[1]
    # rhs       = ex.args[2]
    # lhs_str   = to_string(lhs)
    # rhs_str   = to_string(rhs)
    # total_str = lhs_str * ws * op * ws * rhs_str

    return __to_string_plus_or_prod(ex)
end

function __to_string_comp_minus_comp(ex)

    # a + b - (c +- d)
    # a + b - c */ d

    # if the rhs op is
    # +, - -> paren
    # *, / -> dont paren if the rhs.args are either const or var
    op        = string(ex.head)
    lhs       = ex.args[1]
    rhs       = ex.args[2]
    lhs_str   = to_string(lhs)
    rhs_str   = to_string(rhs)
    ws        = " "
    always_paren   = is_plus_operator(rhs.head) || is_minus_operator(rhs.head)
    possibly_paren = ( is_product_operator(rhs.head) ||
                       is_divide_operator(rhs.head) ) &&
                       !(is_constant(rhs.args[1]) || is_variable(rhs.args[1]))
    do_paren    = always_paren | possibly_paren
    left_paren  = do_paren ? "(" : ""
    right_paren = do_paren ? ")" : ""

    total_str = lhs_str * ws * op * ws * left_paren * rhs_str * right_paren

    return total_str
end

function __to_string_comp_prod_comp(ex)

    # special rule (not implemented yet)
    # if in one branch only * -> no paren is needed

    # (a +- b) * (c +- d)
    # a*/b * c*/d

    # if the rhs op is
    # +, - -> paren
    # *, / -> dont paren if the rhs.args are either const or var
    op        = string(ex.head)
    lhs       = ex.args[1]
    rhs       = ex.args[2]
    lhs_str   = to_string(lhs)
    rhs_str   = to_string(rhs)
    ws        = " "

    always_paren   = is_plus_operator(rhs.head) || is_minus_operator(rhs.head)
    possibly_paren = ( is_product_operator(rhs.head) ||
                       is_divide_operator(rhs.head) ) &&
                       !(is_constant(rhs.args[1]) || is_variable(rhs.args[1]))
    do_paren_right  = always_paren | possibly_paren

    always_paren   = is_plus_operator(lhs.head) || is_minus_operator(lhs.head)
    possibly_paren = ( is_product_operator(lhs.head) ||
                       is_divide_operator(lhs.head) ) &&
                       !(is_constant(lhs.args[1]) || is_variable(lhs.args[1]))
    do_paren_left  = always_paren | possibly_paren

    left_paren  = "("
    right_paren = ")"

    left_str  = do_paren_left ? (left_paren * lhs_str * right_paren) : lhs_str
    right_str = do_paren_right ? (left_paren * rhs_str * right_paren) : rhs_str
    total_str = left_str *  ws * op * ws * right_str

    return total_str
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
module Markovify

include("Tokenizer.jl")

export Model, walk, walk2, combine, makefromdict, state_with_beginning

"""
    Token{T} = Union{Symbol, T}

Tokens can be of any type. They can also include symbols `:begin` and `:end`
which are used to denote the beginning and end of a suptoken.
"""
Token{T} = Union{Symbol, T}

"""
    State{T} = Vector{Token{T}}

A state is described by a succession of tokens.
"""
State{T} = Vector{Token{T}}

"""
    TokenOccurences{T} = Dict{Token{T}, Int}

A dictionary pairing tokens (or special symbols `:begin` and `:end`) with
the number of their respective occurences.
"""
TokenOccurences{T} = Dict{Token{T}, Int}

"""
The datastructure of the Markov chain. Encodes all the different states
and the probabilities of going from one to another as a dictionary. The keys
are the states, the values are the respective
[`TokenOccurences`](@ref) dictionaries. Those are dictionaries which say
how many times was a token found *immediately* after the state.

# Fields
- `order` is the number of tokens in a [`State`](@ref)
- `nodes` is a dictionary pairing [`State`](@ref) and its respective
[`TokenOccurences`](@ref) dictionary.
"""
struct Model{T}
    order::Int
    nodes::Dict{State{T}, TokenOccurences{T}}
end

"""
    combine(chain, others)

Return a Model which is a combination of all of the models provided. All of the
arguments should have the same `order`. The nodes of all the Models are merged
using the function `merge`.
"""
function combine(chain, others...)
    return Model(chain.order, merge(chain.nodes, map(ch -> ch.nodes, others)...))
end

"""
    Model(nodes)

Return a model constructed from `nodes`. Can be used to reconstruct
a model object from its nodes, e.g. if the nodes were saved in a JSON file.
"""
function Model(nodes)
    if isempty(nodes) return nothing end
    # Pick some state and get its length; that's the order of the model
    order = length(nodes[first(keys(nodes))])
    return Model(order, nodes)
end

"""
    begseq(n)

Return the symbol `:begin` repeated `n` times. This array is then used
as a starting sequence for all suptokens.
"""
begseq(n) = fill(:begin, n)

"""
    stdweight(state, token)

A constant `1`. Used as a placeholder function in [`Model`](@ref) to represent
unbiased weight function.
"""
function stdweight(state, token)
    return 1
end

"""
    Model(suptokens::Vector{<:Vector{T}}; order=2, weight=stdweight)

Return a [`Model`](@ref) trained on an array of arrays of [`tokens`](@ref Token) (`suptokens`).
Optionally an `order` of the chain can be supplied; that is
the number of tokens in one state. A weight function of general
type `func(::State{T}, ::Token{T}) -> Int` can be supplied to be used
to bias the weights based on the state or token value.
"""
function Model(suptokens::Vector{<:Vector{T}}; order=2, weight=stdweight) where T
    nodes = Dict{State{T}, TokenOccurences{T}}()
    begin_sequence = begseq(order)
    for incomplete_tokens in suptokens
        # Tokens looks like: :begin :begin ... :begin token ... token :end
        # There will thus be a state consisting of only :begins
        # then a state with :begins and only one token etc.
        tokens = [begin_sequence; incomplete_tokens; [:end]]
        for i in 1:(length(tokens) - order)
            # Current state has a length=order
            state = tokens[i:i+order-1]
            # The token after the current state
            token = tokens[i+order]
            token_counts = get!(nodes, state, Dict{Token{T}, Int}())
            # Add a new occurence (possible with added weight)
            # of the token after the state
            token_counts[token] = get(token_counts, token, 0) + weight(state, token)
        end
    end
    return Model(order, nodes)
end

"""
    walk(model[, init_state])

Return an array of tokens obtained by a random walk through the Markov chain.
The walk starts at state `init_state` if supplied, and at state
`[:begin, :begin...]` (the length depends on
the order of the supplied `model`) otherwise. The walk ends once a special
token `:end` is reached.

See also: [`walk2`](@ref).
"""
function walk(model)
    return walker(model, begseq(model.order), [])
end

function walk(model, init_state)
    return walker(model, init_state, init_state)
end

"""
    states_with_suffix(model, init_suffix)

Return all of the states of `model` that end with `init_suffix`. If
the number of such states is 1 (or 0), the function shortens the suffix
(cuts the first token) in order to lower the requirements, and makes another try.
"""
function states_with_suffix(model, init_suffix)
    hassuffix(ar, suffix) = ar[end-length(suffix)+1:end] == suffix

    # If there is more than one state with given suffix, return the array of them
    # if not, cut the first token out of the suffix and try again
    function helper(suffix)
        states = [k for k in keys(model.nodes) if hassuffix(k, suffix)]
        # Either we got than one valid state (yay!)
        # or the suffix is already so short that we have to end
        if (length(states) > 1) || (length(suffix) <= 1)
            return states
        else
            # Shorten the suffix and try again
            return helper(suffix[2:end])
        end
    end

    return helper(init_suffix)
end

"""
    walk2(model[, init_state])

Return an array of tokens obtained by a random walk through the Markov chain.
When there is only one state following the current one (i.e. there is 100%
chance that the state will become the next one), the function shortens
the current `State` as to lower the requirements and obtain more randomness.
The `State` gets shortened until a state with at least two possible
successors is found (or until `State` is only one token long).

The walk starts at state `init_state` if supplied, and at state
`[:begin, :begin...]` (the length depends on the order of the supplied `model`)
otherwise. The walk ends once a special token `:end` is reached.

See also: [`walk`](@ref).
"""
function walk2(model)
    # First, do the normal append_token operation
    # Then choose all the states that are similar to the new state
    # but have more possible following tokens than only one
    # (or choose the new state itself if it is possibly followed by more than one token)
    # And then, chose a random state from this list
    newstate = rand ∘ (suf -> states_with_suffix(model, suf)) ∘ append_token
    return walker(model, begseq(model.order), [], newstate)
end

function walk2(model, init_state)
    # First, do the normal append_token operation
    # Then choose all the states that are similar to the new state
    # but have more possible following tokens than only one
    # (or choose the new state itself if it is possibly followed by more than one token
    # And then, chose a random state from this list
    newstate = rand ∘ (suf -> states_with_suffix(model, suf)) ∘ append_token
    return walker(model, init_state, init_state, newstate)
end

"""
    append_token(state, token)

Drop the first element in `state` and append
the `token` at the end of the `state` array.
"""
function append_token(state, token)
    return [state[2:end]; [token]]
end

"""
    walker(model, init_state, init_accum, newstate=append_token)

Return an array of tokens obtained by a random walk through the Markov chain.
The walk starts at state `init_state` and ends once a special token `:end`
is reached. A function `newstate` of general type
`func(::State{T}, ::Token{T})::State{T} where T` can be supplied
to be used to generate a new state given the old state and the following token.

This is a general function which is used by all the `walk` functions.

See also: [`walk`](@ref), [`walk2`](@ref).
"""
function walker(model, init_state, init_accum, newstate=append_token)
    function helper(state, accum)
        token = next_token(model, state)
        if token == :end
            return accum
        end
        return helper(newstate(state, token), push!(accum, token))
    end

    return helper(init_state, init_accum)
end

"""
    state_with_prefix(model, prefix; strict=false)

Attempts to return a random valid state of `model` that begins with `tokens`.
If `strict` is `false` and the `model` doesn't have any state that begins
with `tokens`, the function shortens the tokens (cuts the last token)
to lower the requirements and tries to find some valid state again.
"""
function state_with_prefix(model, prefix; strict=false)
    # The token sequence must be at most as long as the model's state
    if length(prefix) > model.order
        message =
            "The length of the prefix must be equal" *
            "to or lower than the order of the model (i.e. $(model.order))"
        throw(DomainError(tokens, message))
    end

    # If the tokens are already a valid state, just return them
    # if the token sequence is too short, just fill in :begin to make a valid state
    if haskey(model.nodes, [begseq(model.order - length(prefix)); prefix])
        return [begseq(model.order - length(prefix)); prefix]
    end

    hasprefix(ar, pref) = ar[1:length(pref)] == pref
    # Try to cut out the last element of the token sequence
    # in order to find a valid state with this given prefix
    function helper(pref, states)
        if pref == [] return nothing end
        states_with_prefix = (st for st in states if hasprefix(st, pref))
        if !isempty(states_with_prefix)
            # Return the non-empty iterator of valid states with iven prefix
            return states_with_prefix
        elseif strict
            # We didn't find any states with given prefix
            return nothing
        else
            return helper(pref[1:end-1], states)
        end
    end

    valid_states = helper(prefix, keys(model.nodes))
    if valid_states != nothing
        return rand(collect(valid_states))
    else
        return nothing
    end
end

"""
    next_token(model, state)

Return a token which will come after the current state, at random.
The probabilities of individual tokens getting choosed
are skewed by their individual values in the `TokenOccurences` dictionary
of the current `state`, that is obtained from the `model`.
"""
function next_token(model, state) where T
    # Choose a random token coming after state
    randkey(model.nodes[state])
end

"""
    randkey(dict)

Return a random key from `dict`. The probabilities of individual keys
getting chosen are skewed by their respective values.
"""
function randkey(dict)
    # Accumulate is similar to scanl in Haskell
    # It folds the array with the given function and returns
    # a list of all intermediate values
    possibility_weights = accumulate(+, collect(values(dict)))
    # Generate a random index in range 1:length(keys(dict))
    # by generating a random number and then placing it into
    # a sorted array of accomulated occurence sums
    index = indexin(possibility_weights, rand() * possibility_weights[end])
    return collect(keys(dict))[index]
end

"""
    indexin(array)

Given a sorted `array`, return the index on which `n` would be inserted in
should the insertion preserve the sorting.
"""
function indexin(array, n)
    for i in 1:length(array)
        if array[i] >= n
            return i
        end
    end
    # If we didn't return yet, n is bigger than every element of array
    return length(array) + 1
end

end


text:
#-----------------------------------------------------------

msg("DataSource: class hierachy")

@test LabeledDataSource <: DataSource
@test InMemoryLabeledDataSource <: LabeledDataSource
@test EncodedInMemoryLabeledDataSource <: InMemoryLabeledDataSource

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

msg("DataSource: abstract methods standard implementation")

immutable FaultyDataSource <: LabeledDataSource
end

source = FaultyDataSource()
@test_throws ArgumentError nobs(source)
@test_throws ArgumentError features(source)
@test_throws ArgumentError features(source, 1, 2)
@test_throws ArgumentError targets(source)
@test_throws ArgumentError targets(source, 1, 2)
@test_throws ArgumentError bias(source)

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

msg("LabeledDataSource: interface stability")

X = [1. 2. 3.;
     4. 5. 6.]
w = [1.,2.]
wn = [1. 2.;
      3. 4.]
ce1 = MultivalueClassEncoding(["V1","V2"])
ce2 = OneOfKClassEncoding(["V1","V2"])
@test_throws DimensionMismatch EncodedInMemoryLabeledDataSource(X, w, ce1)
@test_throws ArgumentError EncodedInMemoryLabeledDataSource(X, w, ce2)
@test_throws DimensionMismatch EncodedInMemoryLabeledDataSource(X, wn, ce2)
@test_throws MethodError EncodedInMemoryLabeledDataSource(X, wn, ce1)

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

msg("LabeledDataSource: constructors")

X = [1. 2. 3.;
     4. 5. 6.]
t = [1., 2., 3.]
tn = [1. 0. 0.;
      0. 0. 1.;
      0. 1. 0.;
      0. 0. 0.]

ce = MultivalueClassEncoding(["V1","V2","V3"])
ds = EncodedInMemoryLabeledDataSource(X, t, ce)
@test nobs(ds) == 3
@test nvar(ds) == 2
@test bias(ds) == 1.
@test features(ds) == X
@test features(ds, 1, 1) == [1. 4.]'
@test features(ds, 2, 2) == [2. 3.; 5. 6.]
@test groundtruth(ds) == ["V1","V2","V3"]
@test targets(ds) == t
@test targets(ds, 1, 1) == [1.]
@test targets(ds, 2, 2) == [2., 3]
@test nclasses(ds) == 3
@test labels(ds) == ["V1","V2","V3"]

ds = dataSource(X, t, ce, bias=0.)
@test typeof(ds) <: EncodedInMemoryLabeledDataSource
@test nobs(ds) == 3
@test nvar(ds) == 2
@test bias(ds) == 0.
@test features(ds) == X
@test features(ds, 1, 1) == [1. 4.]'
@test features(ds, 2, 2) == [2. 3.; 5. 6.]
@test groundtruth(ds) == ["V1","V2","V3"]
@test targets(ds) == t
@test targets(ds, 1, 1) == [1.]
@test targets(ds, 2, 2) == [2., 3]
@test nclasses(ds) == 3
@test labels(ds) == ["V1","V2","V3"]

ce = OneOfKClassEncoding(["V1","V2","V3","V4"])
ds = EncodedInMemoryLabeledDataSource(X, tn, ce, .3)
@test nobs(ds) == 3
@test nvar(ds) == 2
@test bias(ds) == .3
@test features(ds) == X
@test features(ds, 1, 1) == [1. 4.]'
@test features(ds, 2, 2) == [2. 3.; 5. 6.]
@test groundtruth(ds) == ["V1","V3","V2"]
@test targets(ds) == tn
@test targets(ds, 1, 1) == [1. 0. 0. 0.]'
@test targets(ds, 2, 2) == [0. 0.; 0. 1.; 1. 0.; 0. 0.]
@test nclasses(ds) == 4
@test labels(ds) == ["V1","V2","V3","V4"]

ds = dataSource(X, tn, ce, bias=0.)
@test typeof(ds) <: EncodedInMemoryLabeledDataSource
@test nobs(ds) == 3
@test nvar(ds) == 2
@test bias(ds) == 0.
@test features(ds) == X
@test features(ds, 1, 1) == [1. 4.]'
@test features(ds, 2, 2) == [2. 3.; 5. 6.]
@test groundtruth(ds) == ["V1","V3","V2"]
@test targets(ds) == tn
@test targets(ds, 1, 1) == [1. 0. 0. 0.]'
@test targets(ds, 2, 2) == [0. 0.; 0. 1.; 1. 0.; 0. 0.]
@test nclasses(ds) == 4
@test labels(ds) == ["V1","V2","V3","V4"]

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

msg("LabeledDataSource: splitTrainTest")

X = rand(10, 100)
t = round(rand(100) * 4) + 1
ce = MultivalueClassEncoding(["V1","V2","V3","V4","V5"])
cen = OneOfKClassEncoding(["V1","V2","V3","V4","V5"])
tn = labelencode(cen, labeldecode(ce, t))

ds = EncodedInMemoryLabeledDataSource(X, t, ce)
trainDs, testDs = splitTrainTest!(ds, p_train = .7)
@test typeof(trainDs) <: EncodedInMemoryLabeledDataSource
@test typeof(testDs) <: EncodedInMemoryLabeledDataSource
@test nobs(trainDs) + nobs(testDs) == nobs(ds)

ds = EncodedInMemoryLabeledDataSource(X, tn, cen)
trainDs, testDs = splitTrainTest!(ds, p_train = .7)
@test typeof(trainDs) <: EncodedInMemoryLabeledDataSource
@test typeof(testDs) <: EncodedInMemoryLabeledDataSource
@test nobs(trainDs) + nobs(testDs) == nobs(ds)

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

msg("LabeledDataSource: DataFrames")

#using RDatasets
#data = dataset("datasets", "mtcars")
#formula = AM ~ DRat + WT
#ds = dataSource(formula, data)
#@test bias(ds) == 1.