text:
module Traversals

# This module does not export much. The primary reasons are that exporting would increase the
# likelihood of name conflicts, and because the functions here should be used only by developers
# to create new graph algorithms that rely on breadth-first or depth-first search traversals.

using LightGraphs
using SimpleTraits
"""
    abstract type TraversalAlgorithm

`TraversalAlgorithm` is the abstract type used to specify breadth-first traversal (`BFS`) or
depth-first traversal (`DFS`) for the various traversal functions.
"""
abstract type TraversalAlgorithm end

"""
    abstract type AbstractTraversalState

`AbstractTraversalState` is the abstract type used to hold mutable states
for various traversal algorithms (see [`traverse_graph!`](@ref)).

When creating concrete types, you should override the following functions where relevant. These functions
are listed in order of occurrence in the traversal:
- [`initfn!(<:AbstractTraversalState, u::Integer)`](@ref): runs prior to traversal, used to set initial state.
- [`previsitfn!(<:AbstractTraversalState, u::Integer)`](@ref): runs prior to neighborhood discovery for vertex `u`.
- [`visitfn!(<:AbstractTraversalState, u::Integer, v::Integer)`](@ref): runs for each neighbor `v` (newly-discovered or not) of vertex `u`.
- [`newvisitfn!(<:AbstractTraversalState, u::Integer, v::Integer)`](@ref): runs when a new neighbor `v` of vertex `u` is discovered.
- [`postvisitfn!(<:AbstractTraversalState, u::Integer)`](@ref): runs after neighborhood discovery for vertex `u`.
- [`postlevelfn!(<:AbstractTraversalState)`](@ref): runs after each traversal level.

Each of these functions should return a boolean. If the return value of the function is `false`, the traversal will return the state
immediately. Otherwise, the traversal will continue.

For better performance, use the `@inline` directive and make your functions branch-free.
"""
abstract type AbstractTraversalState end

"""
    initfn!(state, u)

Modify [`AbstractTraversalState`](@ref) `state` on initialization of traversal
with source vertices, and return `true` if successful; `false` otherwise.
`initfn!` will be called once for each vertex passed to [`traverse_graph!`](@ref).
"""
@inline initfn!(::AbstractTraversalState, u) = true

"""
    previsitfn!(state, u)

Modify [`AbstractTraversalState`](@ref) `state` before examining neighbors of vertex `u`,
and return `true` if successful; `false` otherwise.
"""
@inline previsitfn!(::AbstractTraversalState, u) = true

"""
    newvisitfn!(state, u, v)

Modify [`AbstractTraversalState`](@ref) `state` when the first edge between `u` and `v` is encountered,
and return `true` if successful; `false` otherwise.
"""
@inline newvisitfn!(::AbstractTraversalState, u, v) = true

"""
    visitfn!(state, u, v)

Modify [`AbstractTraversalState`](@ref) `state` when the edge between `u` and `v` is encountered,
and return `true` if successful; `false` otherwise. Note: `visitfn!` may be called multiple times
per edge, depending on the traversal algorithm, for a function that operates on the first occurrence
only, use [`newvisitfn!`](@ref).
"""
@inline visitfn!(::AbstractTraversalState, u, v) = true

"""
    postvisitfn!(state, u)

Modify [`AbstractTraversalState`](@ref) `state` after having examined all neighbors of vertex `u`,
and return `true` if successful; `false` otherwise.
"""
@inline postvisitfn!(::AbstractTraversalState, u) = true
"""
    postlevelfn!(state)

Modify [`AbstractTraversalState`](@ref) `state` before moving to the next vertex in the traversal algorithm,
and return `true` if successful; `false` otherwise.
"""
@inline postlevelfn!(::AbstractTraversalState) = true

##############
# functions common to both BFS and DFS
##############
"""
    traverse_graph!(g, s, alg, state, neighborfn=outneighbors)
    traverse_graph!(g, ss, alg, state, neighborfn=outneighbors)

Traverse a graph `g` from source vertex `s` / vertices `ss` keeping track of `state`. Return `true` if
traversal finished normally; `false` if one of the visit functions returned `false` (see )

"""
traverse_graph!(g::AbstractGraph, s::Integer, alg, state, neighborfn) = traverse_graph!(g, [s], alg, state, neighborfn)
traverse_graph!(g::AbstractGraph, s::Integer, alg, state) = traverse_graph!(g, [s], alg, state)

struct VisitState{T<:Integer} <: AbstractTraversalState
    visited::Vector{T}
end

@inline function initfn!(s::VisitState, u)
    push!(s.visited, u)
    return true
end
@inline function newvisitfn!(s::VisitState, u, v)
    push!(s.visited, v)
    return true
end

"""
    visited_vertices(g, s, alg)
    visited_vertices(g, ss, alg)

Return a vector representing the vertices of `g` visited in order by [`TraversalAlgorithm`](@ref) `alg`
starting at vertex `s` (vertices `ss`).
"""
function visited_vertices(
    g::AbstractGraph{U},
    ss::AbstractVector,
    alg::TraversalAlgorithm
    ) where U<:Integer

    v = Vector{U}()
    sizehint!(v, nv(g))  # actually just the largest connected component, but we'll use this.
    state = VisitState(v)
    traverse_graph!(g, ss, alg, state)

    return state.visited
end

visited_vertices(g::AbstractGraph, s::Integer, alg::TraversalAlgorithm) = visited_vertices(g, [s], alg)

mutable struct ParentState{T<:Integer} <: AbstractTraversalState
    parents::Vector{T}
end

@inline function newvisitfn!(s::ParentState, u, v)
    s.parents[v] = u
    return true
end

"""
    parents(g, s, alg, neighborfn=outneighbors)

Return a vector of parent vertices indexed by vertex using [`TraversalAlgorithm`](@ref) `alg` starting with
vertex `s`. If `neighborfn` is specified, use the corresponding neighbor-generation function ([`inneighbors`](@ref)r
and [`outneighbors`](@ref) are valid values).

### Performance
This implementation is designed to perform well on large graphs. There are
implementations which are marginally faster in practice for smaller graphs,
but the performance improvements using this implementation on large graphs
can be significant.
"""
function parents(g::AbstractGraph{T}, s::Integer, alg::TraversalAlgorithm, neighborfn=outneighbors) where T
    parents = zeros(T, nv(g))
    state = ParentState(parents)

    traverse_graph!(g, s, alg, state, neighborfn)
    return state.parents
end

include("bfs.jl")
include("dfs.jl")

"""
    tree(p)

Return a directed acyclic graph based on a [`parents`](@ref) vector `p`.
"""
function tree(p::AbstractVector{T}) where T <: Integer
    n = T(length(p))
    t = DiGraph{T}(n)
    for (v, u) in enumerate(p)
        if u > zero(T) && u != v
            add_edge!(t, u, v)
        end
    end
    return t
end

"""
    tree(g, s, alg, neighborfn=outneighbors)

Return a directed acyclic graph based on traversal of the graph `g` starting with source vertex `s`
using algorithm `alg` with neighbor function `neighborfn`.
"""
function tree(g::AbstractGraph, s::Integer, alg::TraversalAlgorithm, neighborfn::Function=outneighbors)
    p = parents(g, s, alg, neighborfn)
    return tree(p)
end

export visited_vertices, parents, distances, topological_sort, tree

end  # module


text:
using NBInclude

@testset "notebook" begin
    @nbinclude(joinpath(@__DIR__, "..", "demo.ipynb"))
end

ext: jl


text:
push!(LOAD_PATH,"../src/")
using SnpParser

Documentermakedocs(
         sitename = "SnpParser.jl",
         modules  = [SnpParser],
         pages=[
                "Home" => "index.md"
               ])

deploydocs(;
    repo="github.com/Sedictious/SnpParser.jl",
)


text:
using ProtoStructs
using Test

@testset "ProtoStructs.jl" begin
    include("test_ProtoStruct.jl")
end


text:
"""
```
sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) + sum(pconvac[c] for c in bus_convs) - pd - gs*1^2
sum(q[a] for a in bus_arcs) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == sum(qg[g] for g in bus_gens) + sum(qconvac[c] for c in bus_convs) - qd + bs*1^2
```
"""
function constraint_power_balance_ac(pm::_PM.AbstractDCPModel, n::Int,  i::Int, bus_arcs, bus_arcs_dc, bus_gens, bus_convs_ac, bus_loads, bus_shunts, pd, qd, gs, bs)
    p = _PM.var(pm, n, :p)
    pg = _PM.var(pm, n, :pg)
    pconv_ac = _PM.var(pm, n, :pconv_ac)
    pconv_grid_ac = _PM.var(pm, n, :pconv_tf_fr)
    v = 1

    JuMP.@constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(pconv_grid_ac[c] for c in bus_convs_ac)  == sum(pg[g] for g in bus_gens) - sum(pd[d] for d in bus_loads) - sum(gs[s] for s in bus_shunts)*v^2)
end

"""
Creates Ohms constraints for DC branches

```
p[f_idx] + p[t_idx] == 0)
```
"""
function constraint_ohms_dc_branch(pm::_PM.AbstractDCPModel, n::Int,  f_bus, t_bus, f_idx, t_idx, r, p)
    p_dc_fr = _PM.var(pm, n, :p_dcgrid, f_idx)
    p_dc_to = _PM.var(pm, n, :p_dcgrid, t_idx)

    JuMP.@constraint(pm.model, p_dc_fr + p_dc_to == 0)
end
"`vdc[i] == vdcm`"
function constraint_dc_voltage_magnitude_setpoint(pm::_PM.AbstractDCPModel, n::Int,  i, vdcm)
    # not used
end

function variable_dcgrid_voltage_magnitude(pm::_PM.AbstractDCPModel; kwrags...)
    # not used nw::Int=pm.cnw
end

function constraint_dc_branch_current(pm::_PM.AbstractDCPModel, n::Int,  f_bus, f_idx, ccm_max, p)
# do nothing
end

#################### TNEP constraints #################
"""
```
sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) + sum(pconvac[c] for c in bus_convs) - pd - gs*1^2
sum(q[a] for a in bus_arcs) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == sum(qg[g] for g in bus_gens) + sum(qconvac[c] for c in bus_convs) - qd + bs*1^2
```
"""
function constraint_power_balance_ac_dcne(pm::_PM.AbstractDCPModel, n::Int, i::Int, bus_arcs, bus_arcs_dc, bus_gens, bus_convs_ac, bus_convs_ac_ne, bus_loads, bus_shunts, pd, qd, gs, bs)
    p = _PM.var(pm, n, :p)
    pg = _PM.var(pm, n, :pg)
    pconv_grid_ac_ne = _PM.var(pm, n, :pconv_tf_fr_ne)
    pconv_grid_ac = _PM.var(pm, n, :pconv_tf_fr)
    pconv_ac = _PM.var(pm, n, :pconv_ac)
    pconv_ac_ne = _PM.var(pm, n, :pconv_ac_ne)
    v = 1
    JuMP.@constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(pconv_grid_ac[c] for c in bus_convs_ac) + sum(pconv_grid_ac_ne[c] for c in bus_convs_ac_ne)  == sum(pg[g] for g in bus_gens) - sum(pd[d] for d in bus_loads) - sum(gs[s] for s in bus_shunts)*v^2)
end

"""
```
sum(p[a] for a in bus_arcs) sum(p_ne[a] for a in bus_arcs_ne) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) + sum(pconvac[c] for c in bus_convs) - pd - gs*1^2
```
"""
function constraint_power_balance_acne_dcne(pm::_PM.AbstractDCPModel, n::Int, i::Int, bus_arcs, bus_arcs_ne, bus_arcs_dc, bus_gens, bus_convs_ac, bus_convs_ac_ne, bus_loads, bus_shunts, pd, qd, gs, bs)
    p = _PM.var(pm, n, :p)
    pg = _PM.var(pm, n, :pg)
    pconv_grid_ac_ne = _PM.var(pm, n, :pconv_tf_fr_ne)
    pconv_grid_ac = _PM.var(pm, n, :pconv_tf_fr)
    pconv_ac = _PM.var(pm, n, :pconv_ac)
    pconv_ac_ne = _PM.var(pm, n, :pconv_ac_ne)
    p_ne = _PM.var(pm,n, :p_ne)
    v = 1
    JuMP.@constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_ne[a] for a in bus_arcs_ne) + sum(pconv_grid_ac[c] for c in bus_convs_ac) + sum(pconv_grid_ac_ne[c] for c in bus_convs_ac_ne)  == sum(pg[g] for g in bus_gens) - sum(pd[d] for d in bus_loads) - sum(gs[s] for s in bus_shunts)*v^2)
end

"""
Creates Ohms constraints for candidate DC branches

```
p[f_idx] + p[t_idx] == 0)
```
"""
function constraint_ohms_dc_branch_ne(pm::_PM.AbstractDCPModel, n::Int, f_bus, t_bus, f_idx, t_idx, r, p)
    p_dc_fr_ne = _PM.var(pm, n, :p_dcgrid_ne, f_idx)
    p_dc_to_ne = _PM.var(pm, n, :p_dcgrid_ne, t_idx)
    JuMP.@constraint(pm.model, p_dc_fr_ne + p_dc_to_ne == 0)
end

function variable_dcgrid_voltage_magnitude_ne(pm::_PM.AbstractDCPModel; nw::Int=pm.cnw, bounded::Bool = true, report::Bool=true)
    # not used
end

function add_dc_bus_voltage_setpoint_ne(sol, pm::_PM.AbstractDCPModel)
    _PM.add_setpoint!(sol, pm, "busdc_ne", "vm", :vdcm_ne, status_name="Vdc", inactive_status_value = 4)
    for (i, bus) in sol["busdc_ne"]
        sol["busdc_ne"]["$i"]["vm"] = 1.0
    end
end


text:
using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra
using IterativeSolvers

# The user can create his own orthogonalization function to use in IAR
function doubleGS_function!(VV, vv, h)
    h[:]=VV'*vv; vv[:]=vv-VV*h; g=VV'*vv; vv[:]=vv-VV*g;
    h[] = h[]+g[]; β=opnorm(vv); vv[:]=vv/β; return β
end
# Then it is needed to create a type to access to this function
abstract type DoubleGS <: IterativeSolvers.OrthogonalizationMethod end
# And then introduce a function dispatch for this new type in order to use
# the defined orthogonalization function
import IterativeSolvers.orthogonalize_and_normalize!
function orthogonalize_and_normalize!(V,v,h,::Type{DoubleGS})
    doubleGS_function!(V, v, h) end

@testset "TIAR" begin
    n=100;
    dep=nep_gallery("dep0",n);

    @bench @testset "accuracy eigenpairs" begin
        (λ,Q)=tiar(dep,σ=1.1,γ=3,neigs=2,v=ones(n),maxit=50,tol=eps()*100,errmeasure=ResidualErrmeasure(dep));
        verify_lambdas(2, dep, λ, Q, eps()*100)
    end

    @testset "Compute as many eigenpairs as possible (neigs=Inf)" begin
        (λ,Q)=tiar(dep,σ=1.1,γ=3,neigs=Inf,v=ones(n),maxit=50,tol=eps()*100,errmeasure=ResidualErrmeasure(dep));
        verify_lambdas(7, dep, λ, Q, eps()*100)
    end

    @testset "orthogonalization" begin

    # NOW TEST DIFFERENT ORTHOGONALIZATION METHODS
    @bench @testset "DGKS" begin
        (λ,Q,Z)=tiar(dep,σ=1.1,γ=3,neigs=4,v=ones(n),maxit=50,tol=eps()*100,errmeasure=ResidualErrmeasure(dep));
        @test opnorm(Z'*Z - I) < 1e-6
     end

     @bench @testset "User provided doubleGS" begin
         (λ,Q,Z)=tiar(dep,σ=1.1,γ=3,neigs=4,v=ones(n),maxit=50,tol=eps()*100);
         @test opnorm(Z'*Z - I) < 1e-6
      end

      @bench @testset "ModifiedGramSchmidt" begin
          (λ,Q,Z)=tiar(dep,σ=1.1,γ=3,neigs=4,v=ones(n),maxit=50,tol=eps()*100);
          @test opnorm(Z'*Z - I) < 1e-6
      end

       @bench @testset "ClassicalGramSchmidt" begin
           (λ,Q,Z)=tiar(dep,σ=1.1,γ=3,neigs=4,v=ones(n),maxit=50,tol=eps()*100);
           @test opnorm(Z'*Z - I) < 1e-6
       end
    end

    # iar and tiar ara mathematically equivalent it maxit<<nep
    # verify the equivalence numerically
    @bench @testset "equivalance with IAR" begin
        (λ_tiar,Q_tiar)=tiar(dep,σ=1.1,γ=3,neigs=2,v=ones(n),maxit=50,tol=eps()*100);
        (λ_iar,Q_iar)=   iar(dep,σ=1.1,γ=3,neigs=2,v=ones(n),maxit=50,tol=eps()*100);

        # Relaxed to only test abs-values to avoid problems if iar / tiar finds
        # the other conjugate solution
        for k=1:size(λ_tiar,1);
            @test minimum(abs.(abs(λ_iar[k]) .- abs.(λ_tiar))) < 1e-6
            @test minimum(map(j-> norm(abs.(Q_tiar[:,k]/Q_tiar[1,k])-abs.(Q_iar[:,j]/Q_iar[1,k])),1:size(λ_iar,1)))<1e-6
        end

    end
    @bench @testset "Solve by projection" begin
        np=1000;

        depp=nep_gallery("dep0",np);
        nn=opnorm(compute_Mder(depp,0));
        errmeasure= (λ,v) -> norm(compute_Mlincomb(depp,λ,v))/nn;

        λ,Q = tiar(depp, σ=0, γ=3, neigs=3, v=ones(np), maxit=50,
                   tol=sqrt(eps()), check_error_every=3,
                   proj_solve=true, inner_solver_method=IARInnerSolver(),
                   errmeasure=errmeasure);

        @test errmeasure(λ[1],Q[:,1])<sqrt(eps())*10
    end

    @testset "Errors thrown" begin
        np=100;
        dep=nep_gallery("dep0",np);
        @test_throws NEPCore.NoConvergenceException (λ,Q)=tiar(dep,σ=2.0,γ=3,neigs=4,v=ones(np),maxit=5,tol=eps()*100,errmeasure=ResidualErrmeasure(dep));
    end

end


text:
"""
    horner(c,x)

Evaluate a polynomial whose coefficients are given in ascending
order in `c`, at the point `x`, using Horner's rule.
"""
function horner(c,x)
    n = length(c)
    y = c[n]
    for k in n-1:-1:1
        y = x*y + c[k]
    end
    return y
end


text:
function optdata(ex::Expr)
    allindices = getallindices(ex)
    cost = Power{:χ}(1,1)
    return Dict{Any, typeof(cost)}(i=>cost for i in allindices)
end

function optdata(optex::Expr, ex::Expr)
    if optex.head == :tuple
        isempty(optex.args) && return nothing
        args = optex.args
        if isa(args[1], Expr) && args[1].head == :call && args[1].args[1] == :(=>)
            indices = Vector{Any}(undef, length(args))
            costs = Vector{Any}(undef, length(args))
            costtype = typeof(parsecost(args[1].args[3]))
            for k = 1:length(args)
                if isa(args[k], Expr) && args[k].head == :call && args[k].args[1] == :(=>)
                    indices[k] = normalizeindex(args[k].args[2])
                    costs[k] = parsecost(args[k].args[3])
                    costtype = promote_type(costtype, typeof(costs[k]))
                else
                    error("invalid index cost specification")
                end
            end
            costs = convert(Vector{costtype}, costs)
        else
            indices = map(normalizeindex, args)
            costtype = Power{:χ,Int}
            costs = fill(Power{:χ,Int}(1,1), length(args))
        end
        return Dict{Any, costtype}(indices[k]=>costs[k] for k = 1:length(args))
    elseif optex.head == :call && optex.args[1] == :!
        allindices = unique(getallindices(ex))
        excludeind = map(normalizeindex, optex.args[2:end])
        cost = Power{:χ}(1,1)
        d = Dict{Any, typeof(cost)}(i=>cost for i in allindices)
        for i in excludeind
            d[i] = 1
        end
        return d
    else
        error("invalid index cost specification")
    end
end

# Process index cost specification for @tensoropt and friends
function parsecost(ex::Expr)
    if ex.head == :call && ex.args[1] == :*
        return *(map(parsecost, ex.args[2:end])...)
    elseif ex.head == :call && ex.args[1] == :+
        return +(map(parsecost, ex.args[2:end])...)
    elseif ex.head == :call && ex.args[1] == :-
        return -(map(parsecost, ex.args[2:end])...)
    elseif ex.head == :call && ex.args[1] == :^
        return ^(map(parsecost, ex.args[2:end])...)
    elseif ex.head == :call && ex.args[1] == :/
        return /(map(parsecost, ex.args[2:end])...)
    elseif ex.head == :call && ex.args[1] == :big
        return big(map(parsecost, ex.args[2:end])...)
    elseif ex.head == :call && ex.args[1] == :float
        return float(map(parsecost, ex.args[2:end])...)
    elseif ex.head == :call && ex.args[1] == :Int128
        return Int128(map(parsecost, ex.args[2:end])...)
    else
        error("invalid index cost specification: $ex")
    end
end
parsecost(ex::Number) = ex
parsecost(ex::Symbol) = Power{ex}(1,1)


text:
module MR1dCNN

using Base: include_package_for_output
const DIR = @__DIR__
const ARCH_PATH = DIR * "/../arch/arch.json"
using Pkg
Pkg.activate(DIR * "/..")
Pkg.status()

@info "Loading modules..."
using BSON
using CUDA
using Flux
using Flux: logitcrossentropy
using Flux.Data: DataLoader
using Flux: onehotbatch, onecold
using Logging: with_logger
using ProgressMeter: Progress, next!
using TensorBoardLogger: TBLogger, tb_overwrite
using Random, Dates, DelimitedFiles, Statistics, JSON

#load data utilities for wrangling datasets and model building
include("DataUtils.jl")
include("ModelUtilities.jl")
include("Args.jl")

export Args, train, args, test

args = Args()

function cast_param(x::A, param::T) where A<:AbstractArray{T} where T <: Real
    convert(eltype(x), param)
end

rng = MersenneTwister(args.seed)

function augment(x::A, ρ::T, loc) where A<:AbstractArray{T} where T<:Real
    if loc == gpu
        x .+ ρ * CUDA.randn(eltype(x), size(x))
    else
        x .+ ρ * randn(eltype(x), size(x))
    end
end

# training loss
function model_loss(x::A, y::B, model::Model, ρ::T, loc) where {C<:Chain, A<:AbstractArray, B<:AbstractArray, T<:Real}
    ŷs  = model(augment(x, ρ, loc))
    mean(logitcrossentropy(ŷᵢ, y) for ŷᵢ in ŷs)
end

# loss over data
function total_loss(data::DataLoader, model::Model)
    l = 0f0
    for (x, y) in data
        ŷs  = model(x)
        l += mean(logitcrossentropy(ŷᵢ, y) for ŷᵢ in ŷs)
    end
    l = l/length(data)
    return l
end

#accuracy over data
function accuracy(data::DataLoader, model::Model)
    acc = zero(Float32)
    for (x, y) in data
        ŷs = model(x)
        out = softmax(mean(ŷs))
        acc += sum(onecold(out) .== onecold(y)) * 1 / size(x,4)
    end
    acc/length(data)
end

# convert model parameters to a suitable data structure for TensorBoard logging
function fill_param_dict!(dict, m, prefix)
    if m isa Chain
        for (i, layer) in enumerate(m.layers)
            fill_param_dict!(dict, layer, prefix*"layer_"*string(i)*"/"*string(layer)*"/")
        end
    else
        for fieldname in fieldnames(typeof(m))
            val = getfield(m, fieldname)
            if val isa AbstractArray
                val = vec(val)
            end
            dict[prefix*string(fieldname)] = val
        end
    end
end

# TensorBoard callback function to log model data and loss every epoch
function TBCallback(logger, train_data, val_data, model)
  param_dict = Dict{String, Any}()
  fill_param_dict!(param_dict, model, "")
  with_logger(logger) do
    @info "model" params=param_dict log_step_increment=0
    @info "train" loss=total_loss(train_data, model) acc=accuracy(train_data, model) log_step_increment=0
    @info "validation" loss=total_loss(val_data, model) acc=accuracy(val_data, model)
  end
end

function training_function(model::Model, Xs::Flux.Data.DataLoader, params::Flux.Zygote.Params, opt::ADAM, ρ::N, loc, progress) where {N<:Real}
    for (xs, ys) in Xs
        train_loss, back = Flux.pullback(() -> model_loss(xs, ys, model, ρ, loc), params)
        grad = back(one(train_loss))
        Flux.Optimise.update!(opt, params, grad)
        next!(progress; showvalues=[(:Loss, train_loss)])
    end
end

@info "Warming up training function..."
tm = build_Model(ARCH_PATH, args.input_dims, 6)
d = DataLoader((rand(Float32, size(tm.paths[1].layers[1].layers[1].weight, 1), 128, 1, 1), onehotbatch([2], [1,2,3,4,5,6])))
training_function(tm, d, params(tm), ADAM(1e-3), Float32(.1), cpu, Progress(1))
accuracy(d, tm)

@info "Ready. Use fields in 'args' struct to change parameter settings."

# training function
function train(args::Args)
    args.seed > 0 && Random.seed!(args.seed)

    if args.cuda && has_cuda_gpu()
        loc = gpu
        @info "Training on GPU"
    else
        loc = cpu
        @info "Training on CPU"
    end

    # load data
    train_data, val_data, T = DataUtils.get_train_validation(DataUtils.data_prep(args.train_dir), readdlm(args.train_dir * "/y_train.txt", Int), args.batch_size, args.train_prop, loc, args.scale, args.shuffle)

    # initialize model
    m = build_Model(DIR * "/../arch/arch.json", args.input_dims, args.nclasses) |> loc
    best_model = build_Model(DIR * "/../arch/arch.json", args.input_dims, args.nclasses)
    #optimizer
    opt = ADAM(args.η)

    ρ = cast_param(train_data.data[1], args.ρ) |> loc

    # parameters
    ps = Flux.params(m)

    # make path for model storage and log output
    !ispath(args.save_path) && mkpath(args.save_path)

    # logging by TensorBoard.jl
    if args.tblogging
        tblogger = TBLogger(args.save_path, tb_overwrite)
    end

    #initialize tracking of accuracy and improvement
    best_acc = 0.0
    last_improvement = 0
    if loc == gpu && CUDA.functional()
        augment(randn(Float32, 2, 4) |> gpu, 0.1f0, gpu)
    end

    # training
    train_steps = 0
    @info "Starting training. Total epochs: $(args.epochs)"
    for epoch = 1:args.epochs
        @info "Epoch $(epoch)"
        progress = Progress(length(train_data))

        training_function(m, train_data, ps, opt, ρ, loc, progress)

        # calculate accuracy on validation set
        vacc = accuracy(val_data, m)
        @info "Validation set accuracy: $(vacc)"

        # log model and loss with TensorBoard
        if args.tblogging
            TBCallback(tblogger, train_data, val_data, m)
        end

        #If accuracy improves, save current model (so saved model is most performant)
        if vacc >= best_acc
            best_model = deepcopy(m)
            best_acc = vacc
            last_improvement = epoch
        end

        #If no improvement in args.lr_patience epochs, reduce learning rate
        if epoch - last_improvement >= args.lr_patience && opt.eta > 1e-6
            opt.eta /= args.γ
            @warn(" -> No improvement in $(args.lr_patience) epochs, reducing learning rate to $(opt.eta)!")
            # reset last_improvement to provide enough time to improve after reducing LR
            last_improvement = epoch
        end

        # Early stopping - Stop if no improvement in args.convergence epochs
        if epoch - last_improvement >= args.convergence
            @warn(" -> No improvement in $(args.convergence) epochs. Approximately converged.")
            break
        end
    end
    if args.save_model
        model_path = abspath(joinpath(args.save_path, "model.bson"))
        let model = best_model |> cpu, data = train_data |> cpu, args = args, transform = T
            BSON.@save model_path model args data transform
            @info "Best model saved: $(model_path)"
        end
    end
    if loc == gpu
        CUDA.reclaim()
    end
    return nothing
end

function test(model_path::String, scale::Bool=true, shuffle::Bool=true)
    saved_model = BSON.load(model_path)
    model = saved_model[:model]
    T = saved_model[:transform]

    test_data = DataUtils.get_test_set(DataUtils.data_prep(DIR*"/../data/test"), readdlm(DIR*"/../data/test/y_test.txt"), T, cpu, scale, shuffle)
    accuracy(test_data, model)
end

end #module


text:
using StringDistances, Test

@test evaluate(Levenshtein(), "", "") == 0
@test evaluate(Levenshtein(), "abc", "") == 3
@test evaluate(Levenshtein(), "", "abc") == 3
@test evaluate(Levenshtein(), "bc", "abc") == 1
@test evaluate(Levenshtein(), "kitten", "sitting") == 3
@test evaluate(Levenshtein(), "saturday", "sunday") == 3

@test evaluate(Levenshtein(), "hi, my name is", "my name is") == 4
@test evaluate(Levenshtein(), "alborgów", "amoniak") == 6

@test evaluate(DamerauLevenshtein(), "", "") == 0
@test evaluate(DamerauLevenshtein(), "abc", "") == 3
@test evaluate(DamerauLevenshtein(), "bc", "abc") == 1
@test evaluate(DamerauLevenshtein(), "fuor", "four") == 1


text:
# Taken from: https://github.com/tkf/Transducers.jl/tree/master/benchmark
if lowercase(get(ENV, "CI", "false")) == "true"
    @info "Executing in CI. Instantiating benchmark environment..."
    using Pkg
    Pkg.activate(@__DIR__)
    Pkg.instantiate()
end

using BenchmarkTools
const SUITE = BenchmarkGroup()
for file in sort([file for file in readdir(@__DIR__) if occursin(r"bench[_0-9]+(.*).jl", file)])
    m = match(r"bench[_0-9]+(.*).jl", file)

    SUITE[m[1]] = include(file)
end


text:
# test using cross-validation and a search algorithm to optimise the value of a hyper-parameter
using RatPack
using Random
using Distributions
using PlotlyJS

# generate large simulated dataset, and run base cross-validation on it, and optimise hyper-parameters
include("test_cv_2.jl")

##############################################3
# optimise theta
r = [10.0, 500.0]
ε = 1.0e-1
score(t) = score_ratings( srule, test_competitions,
                         UpdateIterate( rule=UpdateElo(;r0 = r0, K = 10.0, θ=t ), batch_size=batch_size ),
                         real_ratings
                         )
(x_star1, y_star1, error_estimate1, n_evaluations1) = linesearch( "Golden", x->score(x); range=r, ε=ε)

data1 = scatter(;x=θ_test, y=score_s, mode="lines+markers", name="misfit θ")
data2 = scatter(;x=[θ, θ], y=[s2, s2], mode="markers", marker=attr(symbol="x", size=12), name="Input value")
data3 = scatter(;x=[0.0, 0.0], y=[s3, s3], mode="markers",  marker=attr(symbol="diamond-dot", size=12), name="Deterministic")
data4 = scatter(;x=[x_star1, x_star1], y=[y_star1, y_star1], mode="markers", marker=attr(symbol="diamond-dot", size=12), name="Maximum")

data = [data1, data2, data3, data4]
layout = Layout(;
                xaxis=attr(title="θ" ),
                yaxis=attr(title="(normalised) Brier score" ),
                )
p1 = PlotlyJS.plot( data, layout )
redraw!(p1)
# interestingly the optimal choice of theta is just slightly different from the real value
# but that is just probably the result of using a random simulation

#######################################
# optimise K
r = [0.0, 40.0]
ε = 1.0e-1
score(k) = score_ratings( srule, test_competitions,
                         UpdateIterate( rule=UpdateElo(;r0 = r0, K = k, θ=θ ), batch_size=batch_size ),
                         update_ratings(UpdateIterate( rule=UpdateElo(;r0 = r0, K = k, θ=θ ), batch_size=batch_size ), initial_ratings, input_competitions)
                         )
(x_star2, y_star2, error_estimate2, n_evaluations2) = linesearch( "Golden", x->score(x); range=r, ε=ε)

data1 = scatter(;x=Ks, y=score_k, mode="lines+markers", name="score as a function of K")
data2 = scatter(;x=[0, maximum(Ks)], y=[s2, s2], mode="lines", name="upper bound")
data3 = scatter(;x=[x_star2, x_star2], y=[y_star2, y_star2], mode="markers", marker=attr(symbol="diamond-dot", size=12), name="maximum")

data = [data1, data2, data3]
layout = Layout(;
                xaxis=attr(title="K" ),
                yaxis=attr(title="(normalised) Brier score" ),
                )
p2 = PlotlyJS.plot( data, layout )
redraw!(p2)


text:
# interfacing with LinearAlgebra standard library

#
# BLAS 1
#

LinearAlgebra.rmul!(x::StridedCuArray{<:CublasFloat}, k::Number) =
  scal!(length(x), k, x)

# Work around ambiguity with GPUArrays wrapper
LinearAlgebra.rmul!(x::DenseCuArray{<:CublasFloat}, k::Real) =
  invoke(rmul!, Tuple{typeof(x), Number}, x, k)

function LinearAlgebra.BLAS.dot(DX::StridedCuArray{T}, DY::StridedCuArray{T}) where T<:Union{Float32,Float64}
    n = length(DX)
    n==length(DY) || throw(DimensionMismatch("dot product arguments have lengths $(length(DX)) and $(length(DY))"))
    dot(n, DX, DY)
end

function LinearAlgebra.BLAS.dotc(DX::StridedCuArray{T}, DY::StridedCuArray{T}) where T<:Union{ComplexF32,ComplexF64}
    n = length(DX)
    n==length(DY) || throw(DimensionMismatch("dot product arguments have lengths $(length(DX)) and $(length(DY))"))
    dotc(n, DX, DY)
end

function LinearAlgebra.BLAS.dot(DX::DenseCuArray{T}, DY::DenseCuArray{T}) where T<:Union{ComplexF32,ComplexF64}
    BLAS.dotc(DX, DY)
end

function LinearAlgebra.BLAS.dotu(DX::StridedCuArray{T}, DY::StridedCuArray{T}) where T<:Union{ComplexF32,ComplexF64}
    n = length(DX)
    n==length(DY) || throw(DimensionMismatch("dot product arguments have lengths $(length(DX)) and $(length(DY))"))
    dotu(n, DX, DY)
end

LinearAlgebra.norm(x::DenseCuArray{<:CublasFloat}) = nrm2(x)
LinearAlgebra.BLAS.asum(x::StridedCuArray{<:CublasFloat}) = asum(length(x), x)

function LinearAlgebra.axpy!(alpha::Number, x::StridedCuArray{T}, y::StridedCuArray{T}) where T<:CublasFloat
    length(x)==length(y) || throw(DimensionMismatch("axpy arguments have lengths $(length(x)) and $(length(y))"))
    axpy!(length(x), alpha, x, y)
end

function LinearAlgebra.axpby!(alpha::Number, x::StridedCuArray{T}, beta::Number, y::StridedCuArray{T}) where T<:CublasFloat
    length(x)==length(y) || throw(DimensionMismatch("axpby arguments have lengths $(length(x)) and $(length(y))"))
    axpby!(length(x), alpha, x, beta, y)
end

function LinearAlgebra.rotate!(x::StridedCuArray{T}, y::StridedCuArray{T}, c, s) where T<:CublasFloat
    nx = length(x)
    ny = length(y)
    nx==ny || throw(DimensionMismatch("rotate arguments have lengths $nx and $ny"))
    rot!(nx, x, y, c, s)
end

function LinearAlgebra.reflect!(x::StridedCuArray{T}, y::StridedCuArray{T}, c, s) where T<:CublasFloat
    nx = length(x)
    ny = length(y)
    nx==ny || throw(DimensionMismatch("reflect arguments have lengths $nx and $ny"))
    rot!(nx, x, y, c, s)
    scal!(ny, -real(one(T)), y)
    x, y
end

#
# BLAS 2
#

# GEMV

function gemv_dispatch!(Y::CuVector, A, B, alpha::Number=true, beta::Number=false)
    mA, nA = size(A)

    if nA != length(B)
        throw(DimensionMismatch("second dimension of A, $nA, does not match length of B, $(length(B))"))
    end

    if mA != length(Y)
        throw(DimensionMismatch("first dimension of A, $mA, does not match length of Y, $(length(Y))"))
    end

    if mA == 0
        return Y
    end

    if nA == 0
        return rmul!(Y, 0)
    end

    tA, dA = if A isa Transpose
        'T', parent(A)
    elseif A isa Adjoint
        'C', parent(A)
    else
        'N', A
    end

    T = eltype(Y)
    if T <: CublasFloat && A isa StridedCuArray{T} && B isa StridedCuArray{T}
        gemv!(tA, alpha, dA, B, beta, Y)
    else
        gemm_dispatch!(Y, A, B, alpha, beta)
    end
end

for NT in (Number, Real)
    # NOTE: alpha/beta also ::Real to avoid ambiguities with certain Base methods
    @eval begin
        LinearAlgebra.mul!(Y::CuVector, A::StridedCuMatrix, B::StridedCuVector, a::$NT, b::$NT) =
            gemv_dispatch!(Y, A, B, a, b)
        LinearAlgebra.mul!(Y::CuVector, A::Transpose{<:Any, <:StridedCuVecOrMat}, B::StridedCuVector, a::$NT, b::$NT) =
            gemv_dispatch!(Y, A, B, a, b)
        LinearAlgebra.mul!(Y::CuVector, A::Adjoint{<:Any, <:StridedCuVecOrMat}, B::StridedCuVector, a::$NT, b::$NT) =
            gemv_dispatch!(Y, A, B, a, b)
    end
end

# triangular

## direct multiplication/division
for (t, uploc, isunitc) in ((:LowerTriangular, 'L', 'N'),
                            (:UnitLowerTriangular, 'L', 'U'),
                            (:UpperTriangular, 'U', 'N'),
                            (:UnitUpperTriangular, 'U', 'U'))
    @eval begin
        # Multiplication
        LinearAlgebra.lmul!(A::$t{T,<:DenseCuMatrix},
                            b::StridedCuVector{T}) where {T<:CublasFloat} =
            CUBLAS.trmv!($uploc, 'N', $isunitc, parent(A), b)

        # Left division
        LinearAlgebra.ldiv!(A::$t{T,<:DenseCuMatrix},
                            B::StridedCuVector{T}) where {T<:CublasFloat} =
            CUBLAS.trsv!($uploc, 'N', $isunitc, parent(A), B)
    end
end

## adjoint/transpose multiplication ('uploc' reversed)
for (t, uploc, isunitc) in ((:LowerTriangular, 'U', 'N'),
                            (:UnitLowerTriangular, 'U', 'U'),
                            (:UpperTriangular, 'L', 'N'),
                            (:UnitUpperTriangular, 'L', 'U'))
    @eval begin
        # Multiplication
        LinearAlgebra.lmul!(A::$t{<:Any,<:Transpose{T,<:DenseCuMatrix}},
                            b::DenseCuVector{T}) where {T<:CublasFloat} =
            CUBLAS.trmv!($uploc, 'T', $isunitc, parent(parent(A)), b)
        LinearAlgebra.lmul!(A::$t{<:Any,<:Adjoint{T,<:DenseCuMatrix}},
                            b::DenseCuVector{T}) where {T<:CublasReal} =
            CUBLAS.trmv!($uploc, 'T', $isunitc, parent(parent(A)), b)
        LinearAlgebra.lmul!(A::$t{<:Any,<:Adjoint{T,<:DenseCuMatrix}},
                            b::DenseCuVector{T}) where {T<:CublasComplex} =
            CUBLAS.trmv!($uploc, 'C', $isunitc, parent(parent(A)), b)

        # Left division
        LinearAlgebra.ldiv!(A::$t{<:Any,<:Transpose{T,<:DenseCuMatrix}},
                            B::StridedCuVector{T}) where {T<:CublasFloat} =
            CUBLAS.trsv!($uploc, 'T', $isunitc, parent(parent(A)), B)
        LinearAlgebra.ldiv!(A::$t{<:Any,<:Adjoint{T,<:DenseCuMatrix}},
                            B::StridedCuVector{T}) where {T<:CublasReal} =
            CUBLAS.trsv!($uploc, 'T', $isunitc, parent(parent(A)), B)
        LinearAlgebra.ldiv!(A::$t{<:Any,<:Adjoint{T,<:DenseCuMatrix}},
                            B::StridedCuVector{T}) where {T<:CublasComplex} =
            CUBLAS.trsv!($uploc, 'C', $isunitc, parent(parent(A)), B)
    end
end

#
# BLAS 3
#

# GEMM

function gemm_dispatch!(C::CuVecOrMat, A, B, alpha::Number=true, beta::Number=false)
    if ndims(A) > 2
        throw(ArgumentError("A has more than 2 dimensions"))
    elseif ndims(B) > 2
        throw(ArgumentError("B has more than 2 dimensions"))
    end
    mA, nA = size(A,1), size(A,2)
    mB, nB = size(B,1), size(B,2)

    if nA != mB
        throw(DimensionMismatch("A has dimensions ($mA,$nA) but B has dimensions ($mB,$nB)"))
    end

    if C === A || B === C
        throw(ArgumentError("output matrix must not be aliased with input matrix"))
    end

    if mA == 0 || nA == 0 || nB == 0
        if size(C) != (mA, nB)
            throw(DimensionMismatch("C has dimensions $(size(C)), should have ($mA,$nB)"))
        end
        return LinearAlgebra.rmul!(C, 0)
    end

    tA, dA = if A isa Transpose
        'T', parent(A)
    elseif A isa Adjoint
        'C', parent(A)
    else
        'N', A
    end

    tB, dB = if B isa Transpose
        'T', parent(B)
    elseif B isa Adjoint
        'C', parent(B)
    else
        'N', B
    end

    T = eltype(C)
    if dA isa DenseCuArray && dB isa DenseCuArray &&
       gemmExComputeType(eltype(A), eltype(B), eltype(C), mA, nA, nB) !== nothing
        gemmEx!(tA, tB, alpha, dA, dB, beta, C)
    elseif T <: CublasFloat && dA isa DenseCuArray{T} && dB isa DenseCuArray{T}
        gemm!(tA, tB, alpha, dA, dB, beta, C)
    else
        GPUArrays.generic_matmatmul!(C, A, B, alpha, beta)
    end
end

for NT in (Number, Real)
    # NOTE: alpha/beta also ::Real to avoid ambiguities with certain Base methods
    @eval begin
        LinearAlgebra.mul!(C::CuMatrix, A::StridedCuVecOrMat, B::StridedCuVecOrMat, a::$NT, b::$NT) =
            gemm_dispatch!(C, A, B, a, b)

        LinearAlgebra.mul!(C::CuMatrix, A::Transpose{<:Any, <:StridedCuVecOrMat}, B::StridedCuMatrix, a::$NT, b::$NT) =
            gemm_dispatch!(C, A, B, a, b)
        LinearAlgebra.mul!(C::CuMatrix, A::StridedCuMatrix, B::Transpose{<:Any, <:StridedCuVecOrMat}, a::$NT, b::$NT) =
            gemm_dispatch!(C, A, B, a, b)
        LinearAlgebra.mul!(C::CuMatrix, A::Transpose{<:Any, <:StridedCuVecOrMat}, B::Transpose{<:Any, <:StridedCuVecOrMat}, a::$NT, b::$NT) =
            gemm_dispatch!(C, A, B, a, b)

        LinearAlgebra.mul!(C::CuMatrix, A::Adjoint{<:Any, <:StridedCuVecOrMat}, B::StridedCuMatrix, a::$NT, b::$NT) =
            gemm_dispatch!(C, A, B, a, b)
        LinearAlgebra.mul!(C::CuMatrix, A::StridedCuMatrix, B::Adjoint{<:Any, <:StridedCuVecOrMat}, a::$NT, b::$NT) =
            gemm_dispatch!(C, A, B, a, b)
        LinearAlgebra.mul!(C::CuMatrix, A::Adjoint{<:Any, <:StridedCuVecOrMat}, B::Adjoint{<:Any, <:StridedCuVecOrMat}, a::$NT, b::$NT) =
            gemm_dispatch!(C, A, B, a, b)

        LinearAlgebra.mul!(C::CuMatrix, A::Transpose{<:Any, <:StridedCuVecOrMat}, B::Adjoint{<:Any, <:StridedCuVecOrMat}, a::$NT, b::$NT) =
            gemm_dispatch!(C, A, B, a, b)
        LinearAlgebra.mul!(C::CuMatrix, A::Adjoint{<:Any, <:StridedCuVecOrMat}, B::Transpose{<:Any, <:StridedCuVecOrMat}, a::$NT, b::$NT) =
            gemm_dispatch!(C, A, B, a, b)
    end
end

# triangular

## direct multiplication/division
for (t, uploc, isunitc) in ((:LowerTriangular, 'L', 'N'),
                            (:UnitLowerTriangular, 'L', 'U'),
                            (:UpperTriangular, 'U', 'N'),
                            (:UnitUpperTriangular, 'U', 'U'))
    @eval begin
        # Multiplication
        LinearAlgebra.lmul!(A::$t{T,<:DenseCuMatrix},
                            B::DenseCuMatrix{T}) where {T<:CublasFloat} =
            CUBLAS.trmm!('L', $uploc, 'N', $isunitc, one(T), parent(A), B, B)
        LinearAlgebra.rmul!(A::DenseCuMatrix{T},
                            B::$t{T,<:DenseCuMatrix}) where {T<:CublasFloat} =
            CUBLAS.trmm!('R', $uploc, 'N', $isunitc, one(T), parent(B), A, A)

        # optimization: Base.mul! uses lmul!/rmul! with a copy (because of BLAS)
        LinearAlgebra.mul!(X::DenseCuMatrix{T}, A::$t{T,<:DenseCuMatrix},
                           B::DenseCuMatrix{T}) where {T<:CublasFloat} =
            CUBLAS.trmm!('L', $uploc, 'N', $isunitc, one(T), parent(A), B, X)
        LinearAlgebra.mul!(X::DenseCuMatrix{T}, A::DenseCuMatrix{T},
                           B::$t{T,<:DenseCuMatrix}) where {T<:CublasFloat} =
            CUBLAS.trmm!('R', $uploc, 'N', $isunitc, one(T), parent(B), A, X)

        # Left division
        LinearAlgebra.ldiv!(A::$t{T,<:DenseCuMatrix},
                            B::DenseCuMatrix{T}) where {T<:CublasFloat} =
            CUBLAS.trsm!('L', $uploc, 'N', $isunitc, one(T), parent(A), B)

        # Right division
        LinearAlgebra.rdiv!(A::DenseCuMatrix{T},
                            B::$t{T,<:DenseCuMatrix}) where {T<:CublasFloat} =
            CUBLAS.trsm!('R', $uploc, 'N', $isunitc, one(T), parent(B), A)

        # Matrix inverse
        function LinearAlgebra.inv(x::$t{T, <:CuMatrix{T}}) where T<:CublasFloat
            out = CuArray{T}(I(size(x,1)))
            $t(LinearAlgebra.ldiv!(x, out))
        end
    end
end

## adjoint/transpose multiplication ('uploc' reversed)
for (t, uploc, isunitc) in ((:LowerTriangular, 'U', 'N'),
                            (:UnitLowerTriangular, 'U', 'U'),
                            (:UpperTriangular, 'L', 'N'),
                            (:UnitUpperTriangular, 'L', 'U'))
    @eval begin
        # Multiplication
        LinearAlgebra.lmul!(A::$t{<:Any,<:Transpose{T,<:DenseCuMatrix}},
                            B::DenseCuMatrix{T}) where {T<:CublasFloat} =
            CUBLAS.trmm!('L', $uploc, 'T', $isunitc, one(T), parent(parent(A)), B, B)
        LinearAlgebra.lmul!(A::$t{<:Any,<:Adjoint{T,<:DenseCuMatrix}},
                            B::DenseCuMatrix{T}) where {T<:CublasComplex} =
            CUBLAS.trmm!('L', $uploc, 'C', $isunitc, one(T), parent(parent(A)), B, B)
        LinearAlgebra.lmul!(A::$t{<:Any,<:Adjoint{T,<:DenseCuMatrix}},
                            B::DenseCuMatrix{T}) where {T<:CublasReal} =
            CUBLAS.trmm!('L', $uploc, 'T', $isunitc, one(T), parent(parent(A)), B, B)

        LinearAlgebra.rmul!(A::DenseCuMatrix{T},
                            B::$t{<:Any,<:Transpose{T,<:DenseCuMatrix}}) where {T<:CublasFloat} =
            CUBLAS.trmm!('R', $uploc, 'T', $isunitc, one(T), parent(parent(B)), A, A)
        LinearAlgebra.rmul!(A::DenseCuMatrix{T},
                            B::$t{<:Any,<:Adjoint{T,<:DenseCuMatrix}}) where {T<:CublasComplex} =
            CUBLAS.trmm!('R', $uploc, 'C', $isunitc, one(T), parent(parent(B)), A, A)
        LinearAlgebra.rmul!(A::DenseCuMatrix{T},
                            B::$t{<:Any,<:Adjoint{T,<:DenseCuMatrix}}) where {T<:CublasReal} =
            CUBLAS.trmm!('R', $uploc, 'T', $isunitc, one(T), parent(parent(B)), A, A)

        # optimization: Base.mul! uses lmul!/rmul! with a copy (because of BLAS)
        LinearAlgebra.mul!(X::DenseCuMatrix{T}, A::$t{<:Any,<:Transpose{T,<:DenseCuMatrix}},
                           B::DenseCuMatrix{T}) where {T<:CublasFloat} =
            CUBLAS.trmm!('L', $uploc, 'T', $isunitc, one(T), parent(parent(A)), B, X)
        LinearAlgebra.mul!(X::DenseCuMatrix{T}, A::$t{<:Any,<:Adjoint{T,<:DenseCuMatrix}},
                           B::DenseCuMatrix{T}) where {T<:CublasComplex} =
            CUBLAS.trmm!('L', $uploc, 'C', $isunitc, one(T), parent(parent(A)), B, X)
        LinearAlgebra.mul!(X::DenseCuMatrix{T}, A::$t{<:Any,<:Adjoint{T,<:DenseCuMatrix}},
                           B::DenseCuMatrix{T}) where {T<:CublasReal} =
            CUBLAS.trmm!('L', $uploc, 'T', $isunitc, one(T), parent(parent(A)), B, X)
        LinearAlgebra.mul!(X::DenseCuMatrix{T}, A::DenseCuMatrix{T},
                           B::$t{<:Any,<:Transpose{T,<:DenseCuMatrix}}) where {T<:CublasFloat} =
            CUBLAS.trmm!('R', $uploc, 'T', $isunitc, one(T), parent(parent(B)), A, X)
        LinearAlgebra.mul!(X::DenseCuMatrix{T}, A::DenseCuMatrix{T},
                           B::$t{<:Any,<:Adjoint{T,<:DenseCuMatrix}}) where {T<:CublasComplex} =
            CUBLAS.trmm!('R', $uploc, 'C', $isunitc, one(T), parent(parent(B)), A, X)
        LinearAlgebra.mul!(X::DenseCuMatrix{T}, A::DenseCuMatrix{T},
                           B::$t{<:Any,<:Adjoint{T,<:DenseCuMatrix}}) where {T<:CublasReal} =
            CUBLAS.trmm!('R', $uploc, 'T', $isunitc, one(T), parent(parent(B)), A, X)

        # Left division
        LinearAlgebra.ldiv!(A::$t{<:Any,<:Transpose{T,<:DenseCuMatrix}},
                            B::DenseCuMatrix{T}) where {T<:CublasFloat} =
            CUBLAS.trsm!('L', $uploc, 'T', $isunitc, one(T), parent(parent(A)), B)
        LinearAlgebra.ldiv!(A::$t{<:Any,<:Adjoint{T,<:DenseCuMatrix}},
                            B::DenseCuMatrix{T}) where {T<:CublasReal} =
            CUBLAS.trsm!('L', $uploc, 'T', $isunitc, one(T), parent(parent(A)), B)
        LinearAlgebra.ldiv!(A::$t{<:Any,<:Adjoint{T,<:DenseCuMatrix}},
                            B::DenseCuMatrix{T}) where {T<:CublasComplex} =
            CUBLAS.trsm!('L', $uploc, 'C', $isunitc, one(T), parent(parent(A)), B)

        # Right division
        LinearAlgebra.rdiv!(A::DenseCuMatrix{T},
                            B::$t{<:Any,<:Transpose{T,<:DenseCuMatrix}}) where {T<:CublasFloat} =
            CUBLAS.trsm!('R', $uploc, 'T', $isunitc, one(T), parent(parent(B)), A)
        LinearAlgebra.rdiv!(A::DenseCuMatrix{T},
                            B::$t{<:Any,<:Adjoint{T,<:DenseCuMatrix}}) where {T<:CublasReal} =
            CUBLAS.trsm!('R', $uploc, 'T', $isunitc, one(T), parent(parent(B)), A)
        LinearAlgebra.rdiv!(A::DenseCuMatrix{T},
                            B::$t{<:Any,<:Adjoint{T,<:DenseCuMatrix}}) where {T<:CublasComplex} =
            CUBLAS.trsm!('R', $uploc, 'C', $isunitc, one(T), parent(parent(B)), A)
    end
end


text:
using DataDrivenDiffEq
using ModelingToolkit
using OrdinaryDiffEq
using LinearAlgebra
using Plots
gr()

# Create a
function pendulum(u, p, t)
    x = u[2]
    y = -9.81sin(u[1]) - 0.1u[2]^3 -0.2*cos(u[1])
    return [x;y]
end

u0 = [0.99π; -1.0]
tspan = (0.0, 20.0)
prob = ODEProblem(pendulum, u0, tspan)
sol = solve(prob, Tsit5(), saveat = 0.3)

plot(sol)

# Create the differential data
DX = similar(sol[:,:])
for (i, xi) in enumerate(eachcol(sol[:,:]))
    DX[:,i] = pendulum(xi, [], 0.0)
end

# Create a basis
@variables u[1:2]
@parameters w[1:2]
# Lots of polynomials
polys = Operation[1]
for i ∈ 1:5
    push!(polys, u.^i...)
    for j ∈ 1:i-1
        push!(polys, u[1]^i*u[2]^j)
    end
end

# And some other stuff
h = [w[1]*cos(u[1]); w[2]*sin(u[1]); u[1]*u[2]; u[1]*sin(u[2]); u[2]*cos(u[2]); polys...]

basis = Basis(h, u, parameters = w)
println(basis)

# Get the reduced basis via the sparse regression
# Thresholded Sequential Least Squares, works fine for more data
# than assumptions, converges fast but fails sometimes with too much noise
opt = STRRidge(1e-2)

res = SparseIdentificationResult[]
λs = exp10.(-5:0.1:-1)

function bla(X, DX, Ψ::Basis, thresholds::AbstractArray ; pf::ParetoFront = ParetoFront(), p::AbstractArray = [], t::AbstractVector = [], maxiter::Int64 = 10, opt::T = Optimise.STRRidge(),denoise::Bool = false, normalize::Bool = true, convergence_error = eps()) where {T <: DataDrivenDiffEq.Optimise.AbstractOptimiser, S <: AbstractSortingMethod}
    @assert size(X)[end] == size(DX)[end]
    nx, nm = size(X)
    ny, nm = size(DX)

    θ = Ψ(X, p, t)

    scales = ones(eltype(X), length(Ψ))

    ξ = zeros(eltype(X), length(Ψ))
    Ξ = zeros(eltype(X), length(thresholds), length(Ψ))
    Ξ_opt = zeros(eltype(X), length(Ψ), ny)

    iters = zeros(Int64, length(thresholds))

    denoise ? optimal_shrinkage!(θ') : nothing
    normalize ? DataDrivenDiffEq.normalize_theta!(scales, θ) : nothing

    _iter = Inf
    _thresh = Inf

    @inbounds for i in 1:ny
        for (j, threshold) in enumerate(thresholds)
            set_threshold!(opt, threshold)
            iters[j] = sparse_regression!(ξ, θ, DX[i, :]', maxiter, opt, false, false, convergence_error)
            normalize ? DataDrivenDiffEq.rescale_xi!(ξ, scales) : nothing
            pc = ParetoCandidate([norm(ξ, 0); norm(DX[i, :] .- θ'*ξ)], ξ)
            add_candidate!(pf, pc)
        end
        idx = sort!(pf)
        Ξ_opt[:, i] .= parameter(best(pf))
        iters[idx] <  _iter ? _iter = iters[idx] : nothing
        thresholds[idx] < _thresh ? _thresh = thresholds[idx] : nothing
        empty!(pf)
    end

    ## Create the evaluation
    #@inbounds for i in 1:ny
    #    pf = ParetoFront(x[:, :, i])
    #    sort!(pf)
    #    idx = best(pf)[1]
    #    iters[idx] <  _iter ? _iter = iters[idx] : nothing
    #    thresholds[idx] < _thresh ? _thresh = thresholds[idx] : nothing
    #    Ξ_opt[:, i] = Ξ[idx, i, :]
    #end

    set_threshold!(opt, _thresh)
    #return Ξ_opt
    return SparseIdentificationResult(Ξ_opt, Ψ, _iter, opt, _iter < maxiter, DX, X, p = p)
end

λs = exp10.(-5:0.1:-1)
# Use SR3 with high relaxation (allows the solution to diverge from LTSQ) and high iterations
opt = STRRidge()
wf = WeigthedExponentialSum()
pf = ParetoFront(sorting = wf)
opt = SR3(1e-3, 5.0)
Ψ = bla(sol[:,1:10], DX[:, 1:10], basis, λs, pf = pf, p = [1.0; 1.0], maxiter = 15000, opt = opt)
print_equations(Ψ)

using BenchmarkTools

Juno.@profiler Ψ = bla(sol[:,:], DX[:, :], basis, λs, pf = pf, p = [1.0; 1.0], maxiter = 10, opt = opt)
@btime Ψ = SInDy($sol[:,:], $DX[:, :], $basis, $λs, p = [1.0; 1.0], maxiter = 100, opt = $opt)

println(Ψ)
print_equations(Ψ)


text:
using Test

@testset "(32)(32) no-signaling regression" begin

using BellScenarip: LocalPolytope, ConvexPolytope

# Description:
#   Parses inequalities published at http://www.itp.tuwien.ac.at/~svozil/publ/3-3.html
#   from .txt file. Inequalities are represented in the no-signaling subspace for
#   3 inputs and 2 outputs.
#
# Input: None.
#
# Output: Array, contains 16-dim inequality vectors representing a facet inequality
#
function svozil_pitowsky_inequalities()
    filepath = "./test/regression/files/32-32_svozil_pitowsky_facets.txt"

    (inequalities) = open(filepath) do file
        lines = readlines(file)

        ids = Dict(
            "a1" => 2,
            "a2" => 3,
            "a3" => 4,
            "b1" => 5,
            "b2" => 6,
            "b3" => 7,
            "a1b1" => 8,
            "a1b2" => 9,
            "a1b3" => 10,
            "a2b1" => 11,
            "a2b2" => 12,
            "a2b3" => 13,
            "a3b1" => 14,
            "a3b2" => 15,
            "a3b3" => 16,
        )

        inequalities = []

        for id in 7:2:size(lines)[1]
            line = lines[id]

            lhs_match = match(r"^#\d+:\s(\d)\s\?", line)
            rhs_matches = collect( eachmatch(r"([+-])\s?(\d)?(a\db\d|a\d|b\d|)", line) )

            if size(rhs_matches)[1] > 0
                ineq = zeros(Int64, (1,16))
                ineq[1] = parse(Int64, "-" * lhs_match.captures[1])

                for match in rhs_matches
                    index = ids[match[3]]
                    sign = match[1]
                    multiplier = match.captures[2] === nothing ? "1" : match.captures[2]

                    ineq[index] = parse(Int64, sign * multiplier )
                end

                push!(inequalities, ineq)
            end
        end

        (inequalities)
    end

    inequalities
end

@testset "verifying polytope" begin
    dir = "./test/regression/files/"

    vertices = LocalPolytope.vertices((3,2),(3,2))

    @test size(vertices[1]) == (15,1)

    constraints = ConvexPolytope.facet_constraints(vertices, dir=dir)
    @test constraints[2] == []
    @test constraints[3] == []

    inequalities = constraints[1]
    @test size(inequalities) == (684,)
    @test size(inequalities[1]) == (1,16)
    @test issetequal(inequalities, svozil_pitowsky_inequalities())
end

end


text:
using Test
using GeoParams, LinearAlgebra, DelimitedFiles
@testset "ZirconAge.jl" begin

# declare constant variables
s2y 			= 365.0*24.0*3600.0 							# second to year
ZirconData  	= ZirconAgeData();		# default data

# Parameters for Zircon statistical analysis
n_analyses		= 300											# number of synthetic zircon analyses

# read input file (this step should be skipped, when using tracers)
filename 		= 	"Data/Tt15_st_Bl_rad5_0.0126.txt"				# path to input file
Tt_paths		=	readdlm(filename,',')
n_paths			= 	size(Tt_paths,2)-1							# -1 ro remove first column that corresponds to the time
time_years 		= 	Tt_paths[:,1]./s2y							# time fron seconds to years

# Here a matrix is constructed that only contains the temperature information:
Tt_paths_Temp 	= 	Tt_paths[:,2:end]

# first: test case in which we provide the Tt-path as matrix
prob, ages_eruptible, number_zircons, T_av_time, T_sd_time = compute_zircons_Ttpath(time_years, Tt_paths_Temp, ZirconData=ZirconData)
time_Ma, PDF_zircons, time_Ma_average, PDF_zircon_average  = zircon_age_PDF(ages_eruptible, number_zircons, bandwidth=1e5, n_analyses=n_analyses)

# add tests to check that results are consistent
@test sum(number_zircons[:,200]) == 40479.0
@test sum(number_zircons)==5.411985e6
@test  prob[100] ≈ 5.5432526143365145e-6

# A second way to generate the input is having Vector{Vector} for both time and Tt-path.
time_years_vecs = Vector{Vector{Float64}}(undef,size(Tt_paths_Temp,2));
Tt_paths_Temp_vecs		= Vector{Vector{Float64}}(undef,size(Tt_paths_Temp,2));
for i=1:size(Tt_paths_Temp,2)
	if i<400
		time_years_vecs[i] 		= time_years
		Tt_paths_Temp_vecs[i] 	= Tt_paths_Temp[:,i]
	else
		# slightly adjust the size of the vectors
		time_years_vecs[i] 		= time_years[3:end]
		Tt_paths_Temp_vecs[i] 	= Tt_paths_Temp[3:end,i]
	end

end

# the calling routine is the same, but we get one additional output vector:
time_years1, prob1, ages_eruptible1, number_zircons1, T_av_time1, T_sd_time1  = compute_zircons_Ttpath(time_years_vecs, Tt_paths_Temp_vecs)

# add tests to check that results are consistent
@test sum(number_zircons1[:,200]) == 40479.0
@test sum(number_zircons1)==5.411985e6
@test  prob1[100] ≈ 5.5432526143365145e-6

# Do the same but with a single routine that also returns the PDF's
# Note that given the randomness, you'll always get different results
time_Ma, PDF_zircons, time_Ma_average, PDF_zircon_average, time_years,
	prob2, ages_eruptible, number_zircons2, T_av_time, T_sd_time = compute_zircon_age_PDF(time_years_vecs, Tt_paths_Temp_vecs)

@test sum(number_zircons2[:,200]) == 40479.0
@test sum(number_zircons2)==5.411985e6
@test  prob2[100] ≈ 5.5432526143365145e-6

#=
    # Plot Zircon age probability distribution
	# these are the ploting routines using Makie, which is currently not a dependency of GeoParams (but may become one @ some stage)

	f = Figure()
	Axis(f[1, 1], xlabel = "Age [Myr]", ylabel = "Kernel density [ ]", title = "Zircon age probability distribution")
	for i in 1:length(PDF_zircons)
		CairoMakie.lines!(time_Ma[i]/1e6, PDF_zircons[i], color="gray66",linewidth=0.25)
	end
	CairoMakie.lines!(time_Ma_average/1e6, PDF_zircon_average, color="grey0",linewidth=2.)
	CairoMakie.xlims!(-1e5/1e6,1.5e6/1e6)
	save("Zircon_probability_plot_800kyrs.png",f)

	# plot evolution of the average temperature since the onset of magma injection
	f = Figure()
	Axis(f[1, 1], ylabel = "Temperature [°C]", xlabel = "Time [Myr]", title = "Evolution of the average temperature of the magmatic system")
	CairoMakie.lines!(time_years./1e6, T_av_time, color="grey0",linewidth=0.5)
	f
	save("Average_T.png",f)
=#

end


text:
module Flexpart

using Dates
using NCDatasets
using RecipesBase
using DataStructures
using CSV
using YAML
# using Debugger
# using PyPlot

const FlexpartPath = String

mutable struct FlexpartDir
    path::FlexpartPath
    pathnames::OrderedDict{Symbol, String}
    # FlexpartDir(path::String) = is_fp_dir(path) && new(path)
    function FlexpartDir(path::String)
        pns = try
            OrderedDict(name |> Symbol => p for (name, p) in zip(PATHNAMES_VALUES, pathnames(joinpath(path, PATHNAMES))))
        catch e
            if isa(e, SystemError)
                DEFAULT_PATHNAMES
            else
                throw(e)
            end
        end
        new(path, pns)
    end
end
Base.show(io::IO, fpdir::FlexpartDir) = print(io, fpdir.path)
Base.getindex(fpdir::FlexpartDir, name::Symbol) = fpdir.pathnames[name]
function Base.setindex!(fpdir::FlexpartDir, value::String, name::Symbol)
    fpdir.pathnames[name] = value
end

const OPTIONS_DIR = "options"
const OUTPUT_DIR = "output"
const INPUT_DIR = "input"
const AVAILABLE = "AVAILABLE"
const PATHNAMES = "pathnames"
const PATHNAMES_VALUES = ("options", "output", "input", "available")

DEFAULT_PATHNAMES = OrderedDict(
    k |> Symbol => v for (k, v) in zip(PATHNAMES_VALUES, [
        OPTIONS_DIR, OUTPUT_DIR, INPUT_DIR, AVAILABLE
    ])
)

NEEDED_FILES = [OPTIONS_DIR, OUTPUT_DIR, AVAILABLE, PATHNAMES]

const DEFAULT_FP_DIR = joinpath(@__DIR__, "flexpart_dir_template")

function is_fp_dir(path::FlexpartPath)
    files = readdir(path)
    for needed in NEEDED_FILES
        needed in files || error("$path is not a flexpart directory")
    end
    true
end

function create(name::String; force=false)
    spl = splitpath(name)
    path = length(spl) == 1 ? joinpath(pwd(), name) : name
    !force && ispath(path) && error("$path already exists. force = true is required to remove existing dir")
    cp(DEFAULT_FP_DIR, path, force=force)
    FlexpartDir(path)
end

function pathnames(fpdir::FlexpartDir)
    readlines(joinpath(fpdir.path, PATHNAMES))
end

function pathnames(path::String)
    readlines(path)
end

getdir(fpdir::FlexpartDir, type::Symbol) = isempty(fpdir.pathnames) ? joinpath(fpdir.path, DEFAULT_PATH[type]) |> abspath : joinpath(fpdir.path, fpdir[type]) |> abspath

function write(fpdir::FlexpartDir)
    open(joinpath(fpdir.path, PATHNAMES), "w") do f
        for (k, v) in fpdir.pathnames
            Base.write(f, v*"\n")
        end
    end
end

include("readgrib.jl")
include("utils.jl")
include("FpInput.jl")
include("FpIO.jl")
include("FpPlots.jl")
include("Flexextract.jl")

export
    FlexControl,
    FlexpartDir,
    FlexpartOptions,
    FlexpartOutput,
    FlexextractDir,
    Available,
    FeSource,
    MarsRequest,
    prepare,
    set!,
    set_area!,
    set_steps!,
    ncf_files,
    retrieve,
    # Releases,
    # Releases_ctrl,
    # Release,
    # Outgrid,
    # OutgridN,
    # Command,
    select,
    selected,
    update_available,
    write_options, area2outgrid, format,
    attrib,
    variables2d,
    deltamesh,
    areamesh,
    relloc,
    start_dt,
    end_dt,
    namelist2dict,
    write
end


text:
# Day 13 of Advent of Code 2021
#
# Repeatedly "fold" a matrix of 1/0 values along horizontal
# or vertical lines, and (part 1) count up the number of points
# remaining after one fold, and (part 2) visualize the matrix
# after all folds to yield a dot matrix alphabetic code.
#
# AK, 13/12/2021

infile = "test_input.txt"
infile = "input.txt"

# Read data into a sparse matrix, i.e., a dict of (x,y) pairs
points = Dict()
folds = []
for l in readlines(infile)
    if length(l) == 0
        continue
    end
    if startswith(l, "fold")
        push!(folds, l)
    else
        p = map(n -> parse(Int64, n), split(l, ","))
        points[p] = 1
    end
end

# Process each fold (for part 1, only the first)
part1 = 0
for f in folds

    # Turn "fold along x=3" into "x" and 3
    f = split(f," ")[3]
    dir, n = split(f, "=")
    n = parse(Int64, n)

    # Process each point below/right of fold
    pp = keys(points)
    if dir == "y"
        for p in pp
            if p[2] > n
                newP = [p[1], 2 * n - p[2]]
                points[newP] = 1
                delete!(points, p) # Delete old point
            end
        end
    else   # fold along X
        for p in pp
            if p[1] > n
                newP = [2 * n - p[1], p[2]]
                points[newP] = 1
                delete!(points, p) # Delete old point
            end
        end
    end

    # For part 1, report number of points after first iteration
    if part1 == 0
        global part1 = sum(values(points))
    end

end

# Show part 1, total number of points remaining after first fold
println("Part 1: ", part1)

# Output for part 2: visualize matrix to get 8 capital letters
output = []  # Each row a list of characters
for p in keys(points)
    x, y = map(n -> n + 1, p)
    while length(output) < y
        push!(output, [])
    end
    while length(output[y]) < x
        push!(output[y], " ")
    end
    output[y][x] = "X"
end

# Show the matrix
println("\nPart 2:\n")
for l in output
    println(join(l))
end


text:
using StaticArrays
using Printf
using Images
using OrderedCollections: OrderedDict
using GLMakie
using SLAM

include("kitty.jl")

"""
Visualizer for the SLAM.
Can be used to play in "live" mode as the SLAM is working
or to replay from ReplaySaver.
# Usage:
To use it in "live" preview mode, pass it as a keyword argument,
when creating SlamManager:

```julia
visualizer = Visualizer(;resolution=(900, 600), image_resolution=(370, 1226))
SlamManager(params, camera; visualizer)
```
"""
mutable struct Visualizer
    figure::Figure

    top_grid::GridLayout
    bottom_grid::GridLayout

    map_axis::Axis3
    image_axis::Makie.Axis

    positions::Observable{Vector{Point3f0}}
    position_ids::Dict{Int64, Int64}
    positions_queue::OrderedDict{Int64, Point3f0}

    image::Observable{Matrix{Gray{Float64}}}

    positions_lock::ReentrantLock
end

"""
    Visualizer(;resolution, image_resolution)

# Arguments:
- `resolution`: Initial resolution of the figure in `(width, height)` format.
- `image_resolution`: Initial resolution of the figure in `(height, width)` format.
"""
function Visualizer(;resolution, image_resolution)
    figure = Figure(;resolution)
    top_grid = figure[1, 1:2] = GridLayout()
    bottom_grid = figure[2, 1:2] = GridLayout()

    map_axis = Axis3(
        top_grid[1, 1]; aspect=:data, azimuth=-π / 2, elevation=1)
    image_axis = Makie.Axis(
        top_grid[1, 2]; aspect=DataAspect(),
        leftspinevisible=false, rightspinevisible=false,
        bottomspinevisible=false, topspinevisible=false)

    top_grid[1, 1, Top()] = Label(figure, "Map")
    top_grid[1, 2, Top()] = Label(figure, "Current Image")

    camera_dir_element = MarkerElement(;color=:red, marker="-")
    mappoint_element = MarkerElement(;color=:gray, marker="•")

    Legend(
        bottom_grid[1, 1],
        [mappoint_element, camera_dir_element],
        ["Mappoint", "Camera position"];
        orientation=:horizontal, tellheight=true)
    trim!(figure.layout)

    positions = Observable(Point3f0[])
    position_ids = Dict{Int64, Int64}()
    positions_queue = OrderedDict{Int64, Point3f0}()

    lines!(map_axis, positions; color=:red, quality=1, linewidth=2)

    image = Observable(zeros(Gray{Float64}, image_resolution))
    image!(image_axis, image)

    Visualizer(
        figure, top_grid, bottom_grid,
        map_axis, image_axis,
        positions, position_ids, positions_queue,
        image,
        ReentrantLock())
end

Base.display(v::Visualizer) = Base.display(v.figure)

"""
    set_image!(v::Visualizer, image)

Update image in the visualizer.
It should have the same dimensions as the original.
"""
function set_image!(v::Visualizer, image)
    v.image[] = copy!(v.image[], image)
end

"""
    set_position!(v::Visualizer, position)

Add new position to the camera positions. This updates the plot immediately.
"""
function set_position!(v::Visualizer, position)
    v.positions[] = push!(v.positions[], position)
    autolimits!(v.map_axis)
end

"""
    set_frame_wc!(v::Visualizer, frame_id, wc)

Add new frame `wc` to the visualizer queue.
This is used when other threads are updating the visualizer.
The queue is processed in `process_frame_wc!` method.
"""
function set_frame_wc!(v::Visualizer, frame_id, wc)
    lock(v.positions_lock) do
        base_position = SVector{4, Float64}(0, 0, 0, 1)
        position = (wc * base_position)[[1, 3, 2]]
        v.positions_queue[frame_id] = position
    end
end

"""
    process_frame_wc!(v::Visualizer)

Process pose queue. This updates the plot immediately.
"""
function process_frame_wc!(v::Visualizer)
    lock(v.positions_lock) do
        is_done = isempty(v.positions_queue)
        do_update = !is_done

        while !is_done
            frame_id, position = popfirst!(v.positions_queue)
            is_done = isempty(v.positions_queue)

            pid = get(v.position_ids, frame_id, -1)
            if pid == -1
                v.positions[] = push!(v.positions[], position)
                v.position_ids[frame_id] = length(v.positions[])
            else
                v.positions[][pid] = position
            end
        end

        if do_update
            autolimits!(v.map_axis)
            sleep(1 / 60)
        end
    end
end

function replay_kitty(kitty_dir, save_dir, sequence, n_frames)
    dataset = KittyDataset(kitty_dir, sequence; stereo=false)
    println(dataset)

    save_dir = joinpath(save_dir, "kitty-$sequence")
    isdir(save_dir) || mkdir(save_dir)
    @info "Save directory: $save_dir"

    saver = ReplaySaver()
    SLAM.load!(saver, save_dir)
    @show length(saver.positions), n_frames
    @assert length(saver.positions) == n_frames

    resolution = (900, 600)
    image_resolution = (1241, 376)
    # image_resolution = (1226, 370)
    visualizer = Visualizer(;resolution, image_resolution)
    display(visualizer)

    t1 = time()
    for i in 2:n_frames
        left_frame, right_frame = dataset[i]
        left_frame = Gray{Float64}.(left_frame)

        position = saver.positions[i - 1]
        set_image!(visualizer, rotr90(left_frame))
        set_position!(visualizer, position)

        sleep(1 / 60)
    end
    t2 = time()
    @info "Visualization took: $(t2 - t1) seconds."

    visualizer
end


text:
# Author: Soma S Dhavala
# Given a <concept,concept,score>, and <score type>, populate
# conceptSimilarity table  in learnder_db

using JSON
using Requests
using JuMP

if Pkg.installed("CQL") ==  nothing
  Pkg.clone("https://github.com/somasdhavala/CQL.jl.git")
end
using CQL
#import CQL: connect, disconnect, command, query, asyncCommand

server = "localhost";
con = CQL.connect(server);

#concept1 text,
#concept2 text,
#relation_type text,
#sim double
str_prefix = "insert into learner_db.conceptsimilaritymatrix (concept1,concept2,relation_type, sim) values "
n = 2
relation_type = "mock_rel"
for i in 1:n
  concept_i = "test_concept_$i"
  for j in 1:n
    concept_j = "test_concept_$j"
    sim = rand()
    req = str_prefix * "('$concept_i','$concept_j','$relation_type',$sim);"
    # insert the record to cassandra table
    println(req)
    query = CQL.query(con,req)
  end
end
#
query = CQL.disconnect(con)
# generate concept-ids

max_stars_repo_name: harshavardhanc/sunbird-analytics

max_issues_repo_name: harshavardhanc/sunbird-analytics

max_forks_repo_name: harshavardhanc/sunbird-analytics


text:
# query.jl

### Days of week
const 先勝,友引,先負,仏滅,大安,赤口 = 1,2,3,4,5,6
const qdaysofweek = Dict(1=>"先勝",2=>"友引",3=>"先負",4=>"仏滅",5=>"大安",6=>"赤口")

dayname(dt::Integer) = qdaysofweek[dt]
dayabbr(dt::Integer) = qdaysofweek[dt]
dayname(qdt::QDate) = qdaysofweek[dayofweek(qdt)]
dayabbr(qdt::QDate) = qdaysofweek[dayofweek(qdt)]

# Days of week from 先勝 = 1 to 赤口 = 6
function dayofweek(qdt::QDate)
    qdinfo = QREF.qref(qdt)
    mod1(qdinfo.m + qdinfo.md - 1, 6)
end

# define is先勝/is友引/is先負/is仏滅/is大安/is赤口
for (dow, nm) in qdaysofweek
    @eval ($(Symbol("is$nm")))(qdt::QDate) = dayofweek(qdt) == $dow
end

### Months
const 睦月,如月,弥生,卯月,皐月,水無月 = 1,2,3,4,5,6
const 文月,葉月,長月,神無月,霜月,師走 = 7,8,9,10,11,12
const qmonths = Dict(1=>"睦月",2=>"如月",3=>"弥生",4=>"卯月",5=>"皐月",6=>"水無月",
                     7=>"文月",8=>"葉月",9=>"長月",10=>"神無月",11=>"霜月",12=>"師走")
const leap_prefix = "閏"
monthname(dt::Integer, leap::Bool=false) = leap ? (leap_prefix * qmonths[dt]) : qmonths[dt]
@inline monthabbr(dt::Integer, leap::Bool=false) = monthname(dt, leap)
function monthname(qdt::QDate)
    qdinfo = QREF.qref(qdt)
    monthname(qdinfo.m, qdinfo.leap)
end
function monthabbr(qdt::QDate)
    qdinfo = QREF.qref(qdt)
    monthname(qdinfo.m, qdinfo.leap)
end

function daysinmonth(y::Integer, m::Integer, leap::Bool=false)
    qdinfo = QREF.rqref(y, m, leap)
    QREF.daysinmonth(qdinfo)
end
function daysinmonth(qdt::QDate)
    qdinfo = QREF.qref(qdt)
    QREF.daysinmonth(qdinfo)
end

function isleapyear(qdt::QDate)
    qdinfo = QREF.qref(qdt)
    return QREF.daysinyear(qdinfo) > 360
end

function daysinyear(qdt::QDate)
    qdinfo = QREF.qref(qdt)
    QREF.daysinyear(qdinfo)
end

function monthsinyear(qdt::QDate)
    daysinyear(qdt) ÷ 29
end

dayofyear(y::Integer, m::Integer, d::Integer) = dayofyear(y, m, false, d)
function dayofyear(y::Integer, m::Integer, l::Bool, d::Integer=1)
    QREF.dayofyear(QREF.rqref(y, m, l, d))
end
function dayofyear(qdt::QDate)
    QREF.dayofyear(QREF.qref(qdt))
end


text:
my_tests = ["examples/mcmc_mvtdist.jl", "examples/rand_mvtdist.jl", "examples/rand_tdist.jl"]

println("Running tests:")

for my_test in my_tests
    println(" * $(my_test)")
    include(my_test)
end


text:
using Test, Pkg

ENV["DATADEPS_ALWAYS_ACCEPT"] = true
ENV["JULIA_LOG_LEVEL"] = "WARN"

function include_test(_module)
    println("Starting tests for $_module")
    t = @elapsed include(joinpath(_module, "runtests.jl"))
    println("Completed tests for $_module, $(round(Int, t)) seconds elapsed")
    return nothing
end

@testset "ClimateMachine" begin
    all_tests = isempty(ARGS) || "all" in ARGS ? true : false

    function has_submodule(sm)
        any(ARGS) do a
            a == sm && return true
            first(split(a, '/')) == sm && return true
            return false
        end
    end

    for submodule in [
        "InputOutput",
        "Utilities",
        "Common",
        "Arrays",
        "Atmos",
        "Land",
        "Numerics",
        "Diagnostics",
        "Ocean",
        "Driver",
    ]
        if all_tests || has_submodule(submodule) || "ClimateMachine" in ARGS
            include_test(submodule)
        end
    end
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
import Zygote

module TheisLike
	#push!(LOAD_PATH,"/Users/dharp/source/DPFEHM.jl/src")
	import DifferentiableBackwardEuler
	import DPFEHM

	n = 51
	ns = [n, n]
	steadyhead = 0e0
	sidelength = 200
	thickness = 1.0
	coords, neighbors, areasoverlengths, volumes = DPFEHM.regulargrid2d([-sidelength, -sidelength], [sidelength, sidelength], ns, thickness)
	dirichletnodes = Int[]
	dirichleths = zeros(size(coords, 2))
	for i = 1:size(coords, 2)
		if abs(coords[1, i]) == sidelength || abs(coords[2, i]) == sidelength
			push!(dirichletnodes, i)
			dirichleths[i] = steadyhead
		end
	end

	function getQs(Qs::Vector, rs::Vector)
		return sum(map(i->getQs(Qs[i], rs[i]), 1:length(Qs)))
	end

	function getQs(Q::Number, r::Number)#this splits the Q across all the nodes that are as close as possible to r (should be at least 4)
		dists = map(i->(sqrt(sum(coords[:, i] .^ 2)) - r) ^ 2, 1:size(coords, 2))
		mindist = minimum(dists)
		goodnodes = (dists .≈ mindist)
		Qs = Q * goodnodes / sum(goodnodes)
		return Qs
	end

	function solve_numerical(Qs, T, S, t, rs)
		Ks = fill(T, length(neighbors))#note the thickness is 1
		specificstorage = fill(S, size(coords, 2))#note the thickness is 1
		Qs = getQs(Qs, rs)
		h0 = fill(steadyhead, size(coords, 2))
		p = [Qs; Ks; specificstorage]
		h_gw = DifferentiableBackwardEuler.steps(h0, f_gw, f_gw_u, f_gw_p, f_gw_t, p, 0.0, t; abstol=1e-2, reltol=1e-2)
		goodnode = div(size(coords, 2), 2) + 1
		@assert coords[:, goodnode] == [0, 0]#make sure we are looking at the right "good node"
		return h_gw[goodnode, end] - steadyhead
	end

	function unpack(p)
		@assert length(p) == length(neighbors) + 2 * size(coords, 2)
		Qs, Ks, specificstorage = p[1:size(coords, 2)], p[size(coords, 2) + 1:size(coords, 2) + length(neighbors)], p[size(coords, 2) + length(neighbors) + 1:size(coords, 2) + length(neighbors) + size(coords, 2)]
		return Qs, Ks, specificstorage
	end
	function f_gw(u, p, t)
		Qs, Ks, specificstorage = unpack(p)
		return DPFEHM.groundwater_residuals(u, Ks, neighbors, areasoverlengths, dirichletnodes, dirichleths, Qs, specificstorage, volumes)
	end
	function f_gw_u(u, p, t)
		Qs, Ks, specificstorage = unpack(p)
		retval = DPFEHM.groundwater_h(u, Ks, neighbors, areasoverlengths, dirichletnodes, dirichleths, Qs, specificstorage, volumes)
		return retval
	end
	function f_gw_p(u, p, t)
		Qs, Ks, specificstorage = unpack(p)
		J1 = DPFEHM.groundwater_Qs(u, Ks, neighbors, areasoverlengths, dirichletnodes, dirichleths, Qs, specificstorage, volumes)
		J2 = DPFEHM.groundwater_Ks(u, Ks, neighbors, areasoverlengths, dirichletnodes, dirichleths, Qs, specificstorage, volumes)
		J3 = DPFEHM.groundwater_specificstorage(u, Ks, neighbors, areasoverlengths, dirichletnodes, dirichleths, Qs, specificstorage, volumes)
		return hcat(J1, J2, J3)
	end
	f_gw_t(u, p, t) = zeros(length(u))
end

#Qs = [0.031688, -0.031688]
#T = 0.1
#S = 0.01
#t = 60 * 60 * 24.0 * 30
#rs = [sqrt(2) * 125, sqrt(2) * 25]
#@time @show TheisLike.solve_numerical(Qs, T, S, t, rs)
#@time gradient = Zygote.gradient(Qs->TheisLike.solve_numerical(Qs, T, S, t, rs), Qs)[1]
#@show gradient


text:
struct FlippedMesh{U,D1,T} <: AbstractMesh{U,D1,T}
    mesh::AbstractMesh{U,D1,T}
end

cells(m::FlippedMesh) = (fliporientation(c) for c in cells(m.mesh))
vertices(m::FlippedMesh) = m.mesh.vertices

Base.:-(m::AbstractMesh) = FlippedMesh(m)
Base.:-(m::FlippedMesh) = m.mesh


text:
using NanoclusterModeler
using DifferentialEquations
using ODEInterfaceDiffEq
using Queryverse
using DataFrames
using JLD2, FileIO

# Set parameters for model
imax = 68 # Maximum size cluster considered in model
prefactor = 2.08 # Scaling factor in Ns
Ns = NanoclusterModeler.Ns(imax, prefactor) # Number of surface sites on the cluster
Mplus = 0.05 # Concentration in mM
L     = 6.00 # Concentration in mM
c0 = NanoclusterModeler.c0(imax, Ns, Mplus, L) # Initialize 2D matrix containing intial concentrations of M+ and L
kpij = NanoclusterModeler.kplusij(imax, Ns) # Define 2D matrix containing fraction of free sites on surface of NCs
kmij = NanoclusterModeler.kminusij(imax, Ns) # Define 2D matrix containing fraction of occupied sites on surface of NCs
kp  = 1.0e+3
kb  = 1.0e+5
kub = 1.0e-7
kn  = 1.0e+1
kg1 = 1.0e+4
kd1 = 1.0e-9
kg2 = 1.0e+4
kd2 = 1.0e-9
ka  = 1.0e-3
ke  = 1.0e+3
kc  = 1.0e+3
odes = NanoclusterModeler.odes!

# Set file prefix for output
file = "kp1e3_kb1e5_kn1e1_kg1e4_ka1e-3_ke1e3_kc1e3_M+005_Ns208"

# Time is seconds
tspan = (0.0, 1000000.0)

# Parameters being passed to ode function
# The order in which these are stored is important
p   = (imax, Ns, kpij, kmij,
    kp, kb, kub, kn, kg1, kd1, kg2, kd2, ka, ke, kc)

# Define ode problem and solve
prob = ODEProblem(odes, c0, tspan, p)
sol  = solve(prob, radau(), reltol=1e-10, abstol=1e-10)

# Post processing of data
rows, columns, tlast = size(sol)
t = sol.t

# If solution is negative, it is out of the tolerance. Set to zero.
for t = 1:tlast
    for j = 1:columns
        for i = 1:rows
            if sol[i,j,t] < 0.0
                sol[i,j,t] = 0.0
            end
        end
    end
end

Mplus = sol[1,1,:]
M     = sol[3,2,:]
L     = sol[2,3,:]
ML    = sol[3,3,:]
C_20  = sol[4,2,:]
C_21  = sol[4,3,:]
C_22  = sol[4,4,:]
C_23  = sol[4,5,:]
C_2   = hcat(C_20, C_21, C_22, C_23)

C_tot = zeros(Float64, tlast)
C_bar = zeros(Float64, imax, tlast)
for t = 1:tlast
    C_tot[t]   = sum(sol[3:rows-1, 2:columns-1, t])
    C_bar[:,t] = sum(sol[3:rows-1, 2:columns-1, t], dims=2)
end

file_time = string(file, "_time.csv")
df_t = DataFrame([sol.t, Mplus, M, L, ML, C_20, C_21, C_22, C_23, C_tot])
rename!(df_t, [:t, :Mplus, :M, :L, :ML, :C_20, :C_21, :C_22, :C_23, :C_tot])
save(file_time, df_t)

file_Cbar = string(file, "_cbar.csv")
C_bar_T = transpose(C_bar)
df_Cbar = DataFrame(C_bar_T)
save(file_Cbar, df_Cbar)

P_j = zeros(Float64, imax, Ns[imax]+1, tlast)
for t = 1:tlast
    for j = 2:columns
        for i = 4:rows
            if sol[i,j,t] >= 1e-10
                P_j[i-2,j-1,t] = sol[i,j,t]/C_bar[i-2,t]
            end
        end
    end
end

file_P2D_t = string(file, "_prob_j_t.jld2")
@save file_P2D_t P_j

file_sol = string(file, "_full_sol.jld2")
@save file_sol sol


text:
using Base.Test
include("Palindrome_detection.jl")

# Simple test
@test palindrome("abcdcba")
@test !palindrome("abd")

# Test sets
@testset "palindromes" begin
    @test palindrome("aaaaa")
    @test palindrome("abcba")
    @test palindrome("1")
    @test palindrome("12321")
end

@testset "non-palindromes" begin
    @test !palindrome("abc")
    @test !palindrome("a11")
    @test !palindrome("012")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module TestSyntaxHighlighter

using Test

using OhMyREPL
import  OhMyREPL.Passes.SyntaxHighlighter.SYNTAX_HIGHLIGHTER_SETTINGS

using Tokenize

OhMyREPL.colorscheme!("Monokai16")

b = IOBuffer()
str = "function :foobar, foobar # foobar"
OhMyREPL.test_pass(b, SYNTAX_HIGHLIGHTER_SETTINGS, str)

println("Original string: ", str)
println("Highlighted string: ", String(take!(b)))

println()

str = "(function :foobar, foobar @foobar())# foobar"
OhMyREPL.test_passes(b, OhMyREPL.PASS_HANDLER, str, 3)

println("Original string: ", str)
println("Highlighted string: ", String(take!(b)))
println()

end


text:
__precompile__()
module InfluxDB

export InfluxServer, create_db, query
import Base: write
using HTTP, JSON, DataFrames

include("query.jl")
include("write.jl")
include("server.jl")
include("database.jl")

end # module InfluxDB


text:
using FaultTolerantControl
using Transducers
using Plots

function main()
    N = 4
    d = 2
    θ = [[0, 0], [0.3, 0.1], [0.5, 0.2], [0.7, 0.3], [1, 1]]
    θ_x = θ |> Map(_θ -> _θ[1]) |> collect
    θ_y = θ |> Map(_θ -> _θ[2]) |> collect
    t0 = 0.0
    tf = 1.0
    curve = Bezier(θ, t0, tf)
    curve_params = t0:0.01:tf
    points = curve_params |> Map(curve) |> collect
    points_x = points |> Map(point -> point[1]) |> collect
    points_y = points |> Map(point -> point[2]) |> collect
    fig = plot(points_x, points_y; label="Bezier curve")
    plot!(fig, θ_x, θ_y; st=:scatter, label="control points")
end