text:
########################################################
# Bilevel Parameter Learning via Nonsmooth Trust Region
########################################################

using SparseArrays

using AlgTools
using AlgTools.Util
using AlgTools.LinOps
import AlgTools.Iterate

using VariationalImaging.GradientOps
using VariationalImaging.Util

using LinearAlgebra

export bilevel_learn

#############
# Data types
#############

ImageSize = Tuple{Integer,Integer}
Image = Array{Float64,2}
Primal = Image
Dual = Array{Float64,3}
Parameter = Union{Real,AbstractArray}
Dataset = Tuple{Array{Float64,3},Array{Float64,3}}

#########################
# Iterate initialisation
#########################

function init_rest(x,learning_function::Function,Δ,ds)
    x̄ = copy(x)
    u,fx,gx = learning_function(x,ds)
    ū = copy(u)
    fx̄ = copy(fx)
    gx̄ = copy(gx)
    B = ZeroOp{typeof(x)}()
    #B = IdOp{typeof(x)}()
    return x, x̄, u, ū, fx, gx, fx̄, gx̄, Δ, B
end

############
# Auxiliary Functions
############

function cauchy_point_box(x::Real,Δ,g,B)
    Δmax = 10.0
    γ = min(1,Δmax/norm₂(g))
    t = 0
    gᵗBg = VariationalImaging.Util._dot(g,B(g))
    if gᵗBg ≤ 0 # Negative curvature detected
        t = (Δ/10.0)*γ
    else
        t = min(norm₂(g)^2/gᵗBg,(Δ/10.0)*γ)
    end
    d = -t*g
    x_ = x + d
    if x_ <= 0
       x_ = eps()
    end
    return x_-x
end

function cauchy_point_box(x::AbstractArray,Δ,g,B)
    Δmax = 10.0
    γ = min(1,Δmax/norm₂(g))
    t = 0
    gᵗBg = VariationalImaging.Util._dot(g,B(g))
    if gᵗBg ≤ 0 # Negative curvature detected
        t = (Δ/10.0)*γ
    else
        t = min(norm₂(g)^2/gᵗBg,(Δ/10.0)*γ)
    end
    d = -t*g
    x_ = x + d
    Px = clamp!(x_,eps(),Inf)
    return Px-x
end

############
# Algorithm
############

function bilevel_learn(ds :: Dataset,
    learning_function::Function;
    xinit :: Parameter,
    iterate = Iterate.simple_iterate,
    params::NamedTuple)

    ################################
    # Extract and set up parameters
    ################################

    η₁, η₂ = params.η₁, params.η₂
    β₁, β₂ =  params.β₁, params.β₂
    Δ₀ = params.Δ₀

    ######################
    # Initialise iterates
    ######################

    x, x̄, u, ū, fx, gx, fx̄, gx̄, Δ, B = init_rest(xinit,learning_function,Δ₀,ds)

    ####################
    # Run the algorithm
    ####################

    v = iterate(params) do verbose :: Function

        #println("x=$(norm(x,1)/length(x)), Δ=$Δ, gx=$gx")

        p = cauchy_point_box(x,Δ,gx,B) # solve tr subproblem

        x̄ = x + p  # test new point

        ū,fx̄,gx̄ = learning_function(x̄,ds)
        ρ = (-VariationalImaging.Util._dot(p,gx) - 0.5*VariationalImaging.Util._dot(p,B(p)))/(fx-fx̄) # pred/ared
        println("ρ=$ρ")

        if ρ < η₁               # radius update
            Δ = β₁*Δ
        elseif ρ > η₂
            Δ = β₂*Δ
        else
            Δ = β₁*Δ
        end

        if ρ > η₁
            x = x̄
            u = ū
            fx = fx̄
            gx = gx̄
        end

        ################################
        # Give function value if needed
        ################################
        v = verbose() do
            x, u[:,:,1], fx, norm₂(gx), Δ # just show the first image on the dataset
            #fx,u[:,:,1]
        end

        v
    end

    return x, u, v
end

max_stars_repo_name: dvillacis/BPLDenoising

max_issues_repo_name: dvillacis/BPLDenoising

max_forks_repo_name: dvillacis/BPLDenoising


text:
function time_series_intraday_extended(symbol::String, interval::String="60min", slice::String="year1month1"; client = GLOBAL[], parser = "default")
    @argcheck in(interval, ["1min", "5min", "15min", "30min", "60min"])
    sliceMatch = match(r"year(?<year>\d+)month(?<month>\d+)", slice)
    @argcheck !Compat.isnothing(sliceMatch)
    @argcheck parse(Int, sliceMatch["year"]) > 0
    @argcheck parse(Int, sliceMatch["year"]) < 3
    @argcheck parse(Int, sliceMatch["month"]) > 0
    @argcheck parse(Int, sliceMatch["month"]) < 13
    params = Dict(
        "function"=>"TIME_SERIES_INTRADAY_EXTENDED",
        "symbol"=>symbol,
        "interval"=>interval,
        "slice"=>slice,
        "apikey"=>key(client)
    )
    uri = _build_uri(client.scheme, client.host, "query", params)
    data = retry(_get_request, delays=Base.ExponentialBackOff(n=3, first_delay=5, max_delay=1000))(uri)
    p = _parser(parser, "csv")
    return p(data)
end

function time_series_intraday(symbol::String, interval::String="1min"; client = GLOBAL[], outputsize::String="compact", datatype::Union{String, Nothing}=nothing, parser = "default")
    @argcheck in(interval, ["1min", "5min", "15min", "30min", "60min"])
    @argcheck in(outputsize, ["compact", "full"])
    @argcheck in(datatype, ["json", "csv", nothing])
    params = Dict(
        "function"=>"TIME_SERIES_INTRADAY",
        "symbol"=>symbol,
        "interval"=>interval,
        "outputsize"=>outputsize,
        "datatype"=>datatype,
        "apikey"=>key(client)
    )
    uri = _build_uri(client.scheme, client.host, "query", params)
    data = retry(_get_request, delays=Base.ExponentialBackOff(n=3, first_delay=5, max_delay=1000))(uri)
    p = _parser(parser, datatype)
    return p(data)
end

for func in (:daily, :daily_adjusted, :weekly, :weekly_adjusted, :monthly, :monthly_adjusted)
    x = "time_series_$(func)"
    fname = Symbol(x)
    @eval begin
        function ($fname)(symbol::String; client = GLOBAL[], outputsize::String="compact", datatype::Union{String, Nothing}=nothing, parser = "default")
            @argcheck in(outputsize, ["compact", "full"])
            @argcheck in(datatype, ["json", "csv", nothing])
            params = Dict(
                "function"=>uppercase($x),
                "symbol"=>symbol,
                "outputsize"=>outputsize,
                "datatype"=>isnothing(datatype) ? "csv" : datatype,
                "apikey"=>key(client)
            )
            uri = _build_uri(client.scheme, client.host, "query", params)
            data = retry(_get_request, delays=Base.ExponentialBackOff(n=3, first_delay=5, max_delay=1000))(uri)
            p = _parser(parser, datatype)
            return p(data)
        end

        export $fname
    end
end

ext: jl


text:
# test without resolution
A = ones(Int, 8,8)
edge_count = count_edges(A)
@test length(edge_count) == 0

A = [[ones(Int, 4,4) 2*ones(Int, 4,4)]; [3*ones(Int, 4,4) 4*ones(Int, 4,4)]]
edge_count = count_edges(A)
@test edge_count[(1,2)] == 4
@test edge_count[(1,3)] == 4
@test edge_count[(2,4)] == 4
@test edge_count[(3,4)] == 4

A = ones(Int, 8,8)
A[2:7,2:7] = 2*ones(Int, 6,6)
A[3:6,3:6] = 3*ones(Int, 4,4)
A[4:5,4:5] = 4*ones(Int, 2,2)
edge_count = count_edges(A)
@test edge_count[(1,2)] == 24
@test edge_count[(2,3)] == 16
@test edge_count[(3,4)] == 8

A = ones(Int, 0,0)
edge_count = count_edges(A)
@test length(edge_count) == 0

A = ones(Int, 4,4,4)
A[2:3,2:3,2:3] = 2*ones(Int, 2,2,2)
edge_count = count_edges(A)
@test edge_count[(1,2)] == 24

A = ones(Int, 4,4,4,4)
A[2:3,2:3,2:3,2:3] = 2*ones(Int, 2,2,2,2)
edge_count = count_edges(A)
@test edge_count[(1,2)] == 64

# test with resolution
A = [[ones(Int, 4,4) 2*ones(Int, 4,4)]; [3*ones(Int, 4,4) 4*ones(Int, 4,4)]]
edge_count = count_edges(A, [1,2])
@test edge_count[(1,2)] == 4*2
@test edge_count[(1,3)] == 4*1
@test edge_count[(2,4)] == 4*1
@test edge_count[(3,4)] == 4*2

A = ones(Int, 4,4,4)
A[2:3,2:3,2:3] = 2*ones(Int, 2,2,2)
edge_count = count_edges(A, [1,2,2])
@test edge_count[(1,2)] == 40

A = ones(Int, 4,4,4)
A[2:3,2:3,2:3] = 2*ones(Int, 2,2,2)
@test_throws AssertionError count_edges(A, [2])


text:
type MixPC <: PairCop
    cops::Array{PairCop, 1}
    wgts::FloatVec
end


text:
function scale_range(vmin, vmax, n=1, threshold=100)
    dv = abs(vmax - vmin)  # > 0 as nonsingular is called before.
    meanv = (vmax + vmin) / 2
    offset = if abs(meanv) / dv < threshold
        0.0
    else
        copysign(10 ^ (log10(abs(meanv)) ÷ 1), meanv)
    end
    scale = 10 ^ (log10(dv / n) ÷ 1)
    scale, offset
end

function _staircase(steps)
    n = length(steps)
    result = Vector{Float64}(undef, 2n)
    for i in 1:(n-1)
        @inbounds result[i] = 0.1 * steps[i]
    end
    for i in 1:n
        @inbounds result[i+(n-1)] = steps[i]
    end
    result[end] = 10 * steps[2]
    return result
    # [0.1 .* steps[1:end-1]; steps; 10 .* steps[2]]
end

struct EdgeInteger
    step::Float64
    offset::Float64

    function EdgeInteger(step, offset)
        if step <= 0
            error("Step must be positive")
        end
        new(step, abs(offset))
    end
end

function closeto(e::EdgeInteger, ms, edge)
    tol = if e.offset > 0
        digits = log10(e.offset / e.step)
        tol = max(1e-10, 10 ^ (digits - 12))
        min(0.4999, tol)
    else
        1e-10
    end
    abs(ms - edge) < tol
end

function le(e::EdgeInteger, x)
    # 'Return the largest n: n*step <= x.'
    d, m = divrem(x, e.step)
    if closeto(e, m / e.step, 1)
        d + 1
    else
        d
    end
end

function ge(e::EdgeInteger, x)
    # 'Return the smallest n: n*step >= x.'
    d, m = divrem(x, e.step)
    if closeto(e, m / e.step, 0)
        d
    else
        d + 1
    end
end

"""
A cheaper function that tries to come up with usable tick locations for a given value range
"""
function locateticks(vmin, vmax, n_ideal::Int, _integer::Bool = false, _min_n_ticks::Int = 2)

    _steps = (1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0)
    _extended_steps = _staircase(_steps)

    scale, offset = scale_range(vmin, vmax, n_ideal)

    _vmin = vmin - offset
    _vmax = vmax - offset

    raw_step = (_vmax - _vmin) / n_ideal

    steps = _extended_steps .* scale

    if _integer
        # For steps > 1, keep only integer values.
        filter!(steps) do i
            (i < 1) || (abs(i - round(i)) < 0.001)
        end
    end

    #istep = np.nonzero(steps >= raw_step)[0][0]
    istep = findfirst(1:length(steps)) do i
        @inbounds return steps[i] >= raw_step
    end
    ticks = 1.0:0.1:0.0
    for istep in istep:-1:1
        step = steps[istep]

        if _integer && (floor(_vmax) - ceil(_vmin) >= _min_n_ticks - 1)
            step = max(1, step)
        end
        best_vmin = (_vmin ÷ step) * step

        # Find tick locations spanning the vmin-vmax range, taking into
        # account degradation of precision when there is a large offset.
        # The edge ticks beyond vmin and/or vmax are needed for the
        # "round_numbers" autolimit mode.
        edge = EdgeInteger(step, offset)
        low = le(edge, _vmin - best_vmin)
        high = ge(edge, _vmax - best_vmin)
        ticks = (low:high) .* step .+ best_vmin
        # Count only the ticks that will be displayed.
        # nticks = sum((ticks .<= _vmax) .& (ticks .>= _vmin))

        # manual sum because broadcasting was slow
        nticks = 0
        for t in ticks
            if _vmin <= t <= _vmax
                nticks += 1
            end
        end

        if nticks >= _min_n_ticks
            break
        end
    end

    ticks = ticks .+ offset #(first(ticks) + offset):step(ticks):(last(ticks) + offset)
    vals = filter(x -> vmin <= x <= vmax, ticks)

    # for some reason, the values coming out of this method sometimes are ever so slightly off from
    # their intended decimal representation, causing printed values like
    # 0.08000000001 instead of 0.08
    # so here we round off the numbers to the required number of digits after the decimal point
    exponent = floor(Int, minimum(log10.(abs.(diff(vals)))))
    round.(vals, digits = max(0, -exponent+1))
end

function locateticks(vmin, vmax, width_px, ideal_tick_distance::Float32, _integer::Bool = false, _min_n_ticks::Int = 2)
    # how many ticks would ideally fit?
    n_ideal = round(Int, width_px / ideal_tick_distance) + 1
    locateticks(vmin, vmax, n_ideal, _integer, _min_n_ticks)
end


text:
module NiTE

using ..OpenNI2
using Cxx

import ..OpenNI2: libNiTE2
Libdl.dlopen(libNiTE2, Libdl.RTLD_GLOBAL)

import Base: read, read!

nite_header_path = replace(dirname(libNiTE2), "\/Redist", "\/Include")
@assert isdir(nite_header_path)
addHeaderDir(nite_header_path, kind=C_System)

cxx"""
#include <NiTE.h>
"""

### NiteEnums ###

for name in [
    # JointType
    :JOINT_HEAD,
    :JOINT_NECK,
    :JOINT_LEFT_SHOULDER,
    :JOINT_RIGHT_SHOULDER,
    :JOINT_LEFT_ELBOW,
    :JOINT_RIGHT_ELBOW,
    :JOINT_LEFT_HAND,
    :JOINT_RIGHT_HAND,
    :JOINT_TORSO,
    :JOINT_LEFT_HIP,
    :JOINT_RIGHT_HIP,
    :JOINT_LEFT_KNEE,
    :JOINT_RIGHT_KNEE,
    :JOINT_LEFT_FOOT,
    :JOINT_RIGHT_FOOT,

    # SkeletonState
    :SKELETON_NONE,
    :SKELETON_CALIBRATING,
    :SKELETON_TRACKED,
    :SKELETON_CALIBRATION_ERROR_NOT_IN_POSE,
    :SKELETON_CALIBRATION_ERROR_HANDS,
    :SKELETON_CALIBRATION_ERROR_HEAD,
    :SKELETON_CALIBRATION_ERROR_LEGS,
    :SKELETON_CALIBRATION_ERROR_TORSO,

    # Status
    :STATUS_OK,
    :STATUS_ERROR,
    :STATUS_BAD_USER_ID,
    :STATUS_OUT_OF_FLOW,

    # PostType
    :POSE_PSI,
    :POSE_CROSSED_HANDS,

    # GestureType
    :GESTURE_WAVE,
    :GESTURE_CLICK,
    :GESTURE_HAND_RAISE,
    ]
    cppname = string("nite::", name)
    ex = Expr(:macrocall, Symbol("@icxx_str"), string(cppname, ";"))
    @eval begin
        global const $name = $ex
        @assert isa($name, Cxx.CppEnum)
    end
end

@inline function checkStatus(rc)
    if rc != STATUS_OK
        error("[NiTE error (code: $(rc.val))]:\n$(OpenNI2.getExtendedError())")
    end
    rc
end

### nite::NiTE::xxx ###

function initialize()
    rc = icxx"nite::NiTE::initialize();"
    checkStatus(rc)
end
shutdown() = icxx"nite::NiTE::shutdown();"
getVersion() = icxx"nite::NiTE::getVersion();"

### UserTracker ###

type UserTracker{T}
    handle::T
end
function UserTracker()
    handle = icxx"nite::UserTracker();"
    UserTracker(handle)
end

function create(tracker::UserTracker, device::OpenNI2.DevicePtr)
    icxx"$(tracker.handle).create($(device.handle).get());"
end

destroy(tracker::UserTracker) = icxx"$(tracker.handle).destroy();"

function setSkeletonSmoothingFactor(tracker::UserTracker, v)
    icxx"$(tracker.handle).setSkeletonSmoothingFactor($v);"
end

function convertJointCoordinatesToDepth(tracker::UserTracker, x, y, z)
    outX = Cfloat[1]
    outY = Cfloat[1]
    icxx"$(tracker.handle).convertJointCoordinatesToDepth(
            $x,$y,$z,$(pointer(outX)),$(pointer(outY)));"
    outX[1],outY[1]
end

for f in [
    :startSkeletonTracking,
    :stopSkeletonTracking,
    ]
    ex = Expr(:macrocall, Symbol("@icxx_str"), "\$(tracker.handle).$f(\$id);")
    @eval $f(tracker::UserTracker, id) = $ex
end

### UserMap ###

type UserMap{T}
    handle::T
end

typealias UserId Int16 # short int

for f in [
    :getWidth,
    :getHeight,
    :getStride,
    ]
    ex = Expr(:macrocall, Symbol("@icxx_str"), "\$(umap.handle).$f();")
    @eval $f(umap::UserMap) = $ex
end

function getPixels(umap::UserMap)
    rawpixels = icxx"$(umap.handle).getPixels();"::Ptr{UserId}
    unsafe_wrap(Array, rawpixels, getWidth(umap) * getHeight(umap))
end

### UserTrackerFrameRef ###

type UserTrackerFrameRef{T}
    handle::T
end
function UserTrackerFrameRef()
    handle = icxx"nite::UserTrackerFrameRef();"
    UserTrackerFrameRef(handle)
end

for f in [
    :isValid,
    :release,
    :getFloorConfidence,
    :getFloor,
    :getDepthFrame,
    :getTimestamp,
    :getFrameIndex,
    ]
    ex = Expr(:macrocall, Symbol("@icxx_str"), "\$(frame.handle).$f();")
    @eval $f(frame::UserTrackerFrameRef) = $ex
end

function getUserMap(frame::UserTrackerFrameRef)
    UserMap(icxx"$(frame.handle).getUserMap();")
end

function readFrame(tracker::UserTracker, frameRef::UserTrackerFrameRef)
    rc = icxx"$(tracker.handle).readFrame(&$(frameRef.handle));"
    checkStatus(rc)
end
function read!(tracker::UserTracker, frameRef::UserTrackerFrameRef)
    readFrame(tracker, frameRef)
end
function read(tracker::UserTracker)
    frame = UserTrackerFrameRef()
    read!(tracker, frame)
end

typealias UserData Union{cxxt"nite::UserData", cxxt"nite::UserData&"}

for f in [
    :getId,
    :getBoundingBox,
    :getCenterOfMass,
    :isNew,
    :isVisible,
    :isLost,
    :getSkeleton,
    ]
    ex = Expr(:macrocall, Symbol("@icxx_str"), "\$user.$f();")
    @eval $f(user::UserData) = $ex
end

# Wrap nite::Array<nite::UserData> into a Julia type to give array-like access
type UserDataArray{T}
    array::T
end
Base.start(ar::UserDataArray) = 0
Base.next(ar::UserDataArray,i) = (ar[i], i+1)
Base.done(ar::UserDataArray,i) = i >= length(ar)
Base.getindex(ar::UserDataArray,i) = icxx"($(ar.array))[$i];"
Base.length(ar::UserDataArray) = icxx"$(ar.array).getSize();"

function getUsers(frame::UserTrackerFrameRef)
    users = icxx"$(frame.handle).getUsers();"
    UserDataArray(users)
end

### Skeleton ###

typealias Skeleton Union{cxxt"nite::Skeleton", cxxt"nite::Skeleton&"}

function getJoint(s::Skeleton, typ)
    icxx"$s.getJoint($typ);"
end
function getJoint(s::Skeleton, typ::Integer)
    icxx"$s.getJoint(static_cast<nite::JointType>($typ));"
end

getState(s::Skeleton) = icxx"$s.getState();"

### SkeletonJoint ###

typealias SkeletonJoint Union{cxxt"nite::SkeletonJoint", cxxt"nite::SkeletonJoint&"}

for f in [
    :getType,
    :getPosition,
    :getPositionConfidence,
    :getOrientation,
    :getOrientationConfidence,
    ]
    ex = Expr(:macrocall, Symbol("@icxx_str"), "\$s.$f();")
    @eval $f(s::SkeletonJoint) = $ex
end

### Point3f ###

typealias Point3f Union{cxxt"nite::Point3f", cxxt"nite::Point3f&"}

function Base.show(io::IO, p::Point3f)
    x = icxx"$p.x;"
    y = icxx"$p.y;"
    z = icxx"$p.z;"
    # println(io, string(typeof(p)));
    println(io, "nite::Point3f")
    print(io, "(x,y,z): ");
    print(io, (x,y,z))
end

### Quaternion ###

typealias Quaternion Union{cxxt"nite::Quaternion", cxxt"nite::Quaternion&"}

function Base.show(io::IO, p::Quaternion)
    x = icxx"$p.x;"
    y = icxx"$p.y;"
    z = icxx"$p.z;"
    w = icxx"$p.w;"
    # println(io, string(typeof(p)));
    println(io, "nite::Quaternion")
    print(io, "(x,y,z,w): ");
    print(io, (x,y,z,w))
end

end # module


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2020
#
# This is a Julia package for discrete event simulation
#
using DiscreteEvents

println("... basic tests: printing  ...")
resetClock!(𝐶)
println(𝐶)
println(𝐶, stdout, 1, 2, 3)
DiscreteEvents.prettyClock(true)
str = "Clock 1: state=:idle, t=0.0, Δt=0.01, prc:0\n  scheduled ev:0, cev:0, sampl:0\n"
@test repr(𝐶) == str
# DiscreteEvents.prettyClock(false)
# str = "Clock{Array{DiscreteEvents.ClockChannel,1}}(1, DiscreteEvents.ClockChannel[], DiscreteEvents.Idle(), 0.0, , 0.01, DiscreteEvents.Schedule(DataStructures.PriorityQueue{DiscreteEvents.DiscreteEvent,Float64,Base.Order.ForwardOrdering}(), DiscreteEvents.DiscreteCond[], DiscreteEvents.Sample[]), Dict{Any,Prc}(), Channel[], 0.0, 0.0, 0.0, 0, 0)"
# @test repr(𝐶) == str
DiscreteEvents.prettyClock(true)
@test tau() == 0

println("... basic tests: only events  ...")
ex1 = :(1+1)
ex2 = :(1+2)
e() = 123
f(a) = a+3
g(a) = a+4
h(a, b; c = 1, d = 2) = a + b + c + d
i(; a = 1, b = 2) = a + b
j(x) = x == :unknown

@test DiscreteEvents._evaluate(fun(e)) == 123
@test DiscreteEvents._evaluate(fun(f, 1)) == 4
@test DiscreteEvents._evaluate(fun(h, 1, 2, c=3, d=4)) == 10

a = 11; b = 12; c = 13; d = 14;
sf1 = fun(h, a, b, c=c, d=d)
sf2 = fun(h, :a, :b, c=:c, d=:d)
sf3 = fun(h, a, b, c=c, d=d)
sf4 = fun(h, :a, :b, c=:c, d=:d)
sf5 = fun(i, a=a, b=b)
sf6 = fun(event!, fun(f, a), 1, cid=2, spawn=true)
@test sf1() == 50
@test sf2() == 50
@test sf3() == 50
@test sf4() == 50
@test sf5() == 23
sf6()
ev = DiscreteEvents._nextevent(𝐶)
@test ev.ex() == 14
a = 21; b = 22; c = 23; d = 24;
@test sf1() == 50
@test sf2() == 90
@test fun(h, :a, 2, c=:c, d=4)() == 50
@test fun(j, :unknown)()
@test fun(<=, fun(tau), 1)()

@test fun(i, a=10, b=20)() == 30

# one expression
ev = DiscreteEvents.DiscreteEvent(:(1+1), 10.0, nothing, 1)
@test DiscreteEvents._evaluate(ev.ex) == 2
@test ev.t == 10

# two expressions
ev = DiscreteEvents.DiscreteEvent((:(1+1), :(1+2)), 15.0, nothing, 1)
@test DiscreteEvents._evaluate(ev.ex) == (2, 3)
@test ev.t == 15

# one fun
ev = DiscreteEvents.DiscreteEvent(fun(f, 1), 10.0, nothing, 1)
@test DiscreteEvents._evaluate(ev.ex) == 4

# two funs
ev = DiscreteEvents.DiscreteEvent((fun(f, 1), fun(g, 1)), 10.0, nothing, 1)
@test DiscreteEvents._evaluate(ev.ex) == (4, 5)

# expressions and funs mixed in a tuple
ev = DiscreteEvents.DiscreteEvent((:(1+1), fun(g,2), :(1+2), fun(f, 1)), 10.0, nothing, 1)
@test sum([DiscreteEvents._evaluate(ex) for ex in ev.ex if ex isa Function]) == 10
@test sum([eval(ex) for ex in ev.ex if isa(ex, Expr)]) == 5
@test DiscreteEvents._evaluate(ev.ex) == (2, 6, 3, 4)

@test DiscreteEvents._scale(0) == 1
@test DiscreteEvents._scale(pi*1e7) == 1e7

clk = Clock()  # set up clock without sampling
@test_warn "undefined transition" DiscreteEvents.step!(clk, clk.state, DiscreteEvents.Resume())
DiscreteEvents.init!(clk)
@test clk.state == DiscreteEvents.Idle()
@test tau(clk) == 0
clk = Clock(t0=100)
@test tau(clk) == 100
# @test_warn "nothing to _evaluate" incr!(clk)


text:
module OpenTelemetryInstrumentationGenie

using OpenTelemetryAPI
using Genie
using URIs
using TOML

const ACTION_BEFORE_INSTRUMENT = :ACTION_BEFORE_INSTRUMENT

const PKG_VERSION =
    VersionNumber(TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"])
const INSTRUMENTATION_INFO = InstrumentationInfo(
    name = string(@__MODULE__),
    version = PKG_VERSION,
    schema_url = "https://oolong.dev/OpenTelemetry.jl/dev/OpenTelemetryInstrumentationGenie/",
)

const GENIE_TRACER = Ref{Tracer}()
const GENIE_METRICS = Ref{Counter{UInt}}()

function instrument(req, resp, params)
    if haskey(params, Genie.PARAMS_ROUTE_KEY)
        r = params[Genie.PARAMS_ROUTE_KEY]
        a = r.action
        r.action =
            () -> with_context(extract(req.headers, TraceContextTextMapPropagator())) do
                with_span(
                    "Genie $(r.path)",
                    GENIE_TRACER[];
                    kind = SPAN_KIND_SERVER,
                    attributes = Dict{String,TAttrVal}(
                        "http.method" => r.method,
                        "http.target" => req.target,
                        "http.request_content_length" => length(req.body),
                        "http.route" => r.path,
                    ),
                ) do
                    res = a() |> Genie.Router.to_response
                    s = current_span()
                    s["http.status_code"] = string(res.status)
                    s["http.response_content_length"] = length(res.body)
                    GENIE_METRICS[](; route = r.path, status = Int(res.status))
                    r.action = a  # recover action
                    res
                end
            end
    end
    req, resp, params
end

"""
## Metrics

| Meter Name | Instrument Name | Instrument Type | Unit | Dimensions        | Description         |
|:---------- |:--------------- |:--------------- |:---- |:----------------- |:------------------- |
| `http`     | `request`       | `Counter{UInt}` |      | `route`, `status` | Number of requests. |

## Spans

We try to follow [semantic conventions for HTTP
spans](https://github.com/open-telemetry/opentelemetry-specification/blob/main/specification/trace/semantic_conventions/http.md)
here.

The following attributes are added on span creation:

  - `http.method`
  - `http.target`
  - `http.request_content_length`
  - `http_route`

The following attributes are added when response is received:

  - `http.status_code`
  - `http.response_content_length`
"""
function init(;
    tracer_provider = global_tracer_provider(),
    meter_provider = global_meter_provider(),
)
    if instrument in Genie.Router.content_negotiation_hooks
        # already initialized
    else
        push!(Genie.Router.content_negotiation_hooks, instrument)
    end

    GENIE_METRICS[] = Counter{UInt}(
        "Genie",
        Meter(
            "request";
            provider = meter_provider,
            instrumentation_info = INSTRUMENTATION_INFO,
        );
        unit = "",
        description = "Number of requests received.",
    )
    GENIE_TRACER[] =
        Tracer(; provider = tracer_provider, instrumentation_info = INSTRUMENTATION_INFO)
end

function __init__()
    init()
end

end # module


text:
using Oceananigans
using Oceananigans.Units
using OceanTurbulenceParameterEstimation
using Distributions
using DataDeps
using GLMakie

using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities:
    CATKEVerticalDiffusivity,
    SurfaceTKEFlux,
    MixingLength

mixing_length = MixingLength(Cᴬu   = 0.0,
                             Cᴬc   = 0.0,
                             Cᴬe   = 0.0,
                             CᴷRiᶜ = 2.0,
                             Cᵟu   = 0.5,
                             Cᵟc   = 0.5,
                             Cᵟe   = 0.5,
                             Cᴷuʳ  = 0.0,
                             Cᴷcʳ  = 0.0,
                             Cᴷeʳ  = 0.0)

catke = CATKEVerticalDiffusivity(; mixing_length)

function perturbation_prior(θ★, ϵ=0.5)
    L = (1 - ϵ) * θ★
    U = (1 + ϵ) * θ★
    return ScaledLogitNormal(bounds=(L, U))
end

priors = NamedTuple(name => perturbation_prior(θ_constant_Ri[name]) for name in keys(θ_variable_Ri))
free_parameters = FreeParameters(priors)

case_path(case) = @datadep_str("six_day_suite_4m/$(case)_instantaneous_statistics.jld2")

cases = ["free_convection",
         "weak_wind_strong_cooling",
         "strong_wind_weak_cooling",
         "strong_wind",
         "strong_wind_no_rotation"]

field_names = (:u, :v, :T, :e)
regrid_size = nothing
observations = [SyntheticObservations(case_path(case); field_names, regrid_size)
                for case in cases]

observation = first(observations)
α = observation.metadata.parameters.thermal_expansion_coefficient
equation_of_state = LinearEquationOfState(; α)

Nensemble = 3
buoyancy = SeawaterBuoyancy(; equation_of_state, constant_salinity=true)

simulation = ensemble_column_model_simulation(observations;
                                              Nensemble,
                                              buoyancy,
                                              architecture = CPU(),
                                              tracers = (:T, :e),
                                              closure = catke)

simulation.Δt = 1.0
simulation.stop_time = 2days

progress(sim) = @info "Iter: $(iteration(sim)), time: $(prettytime(sim))"
simulation.callbacks[:progress] = Callback(progress, IterationInterval(100))

model = simulation.model
simulation.output_writers[:fields] =
    JLD2OutputWriter(model, merge(model.velocities, model.tracers, model.diffusivity_fields),
                     prefix = "catke_test_simulation",
                     schedule = TimeInterval(10minutes),
                     force = true)

Qᵘ = simulation.model.velocities.u.boundary_conditions.top.condition
Qᵀ = simulation.model.tracers.T.boundary_conditions.top.condition
dTdz = simulation.model.tracers.T.boundary_conditions.bottom.condition

for (case, obs) in enumerate(observations)
    @show case cases[case]
    @show obs.metadata.parameters.momentum_flux
    @show obs.metadata.parameters.temperature_flux
    @show obs.metadata.parameters.buoyancy_flux
    @show f = obs.metadata.parameters.coriolis_parameter

    view(Qᵘ, :, case) .= obs.metadata.parameters.momentum_flux
    view(Qᵀ, :, case) .= obs.metadata.parameters.temperature_flux
    view(dTdz, :, case) .= obs.metadata.parameters.dθdz_deep
    view(simulation.model.coriolis, :, case) .= Ref(FPlane(f=f))
end

calibration = InverseProblem(observations, simulation, free_parameters)

forward_run!(calibration, [θ_variable_Ri, θ_constant_Ri]; suppress=false)

output_path = simulation.output_writers[:fields].filepath
u = FieldTimeSeries(output_path, "u")
v = FieldTimeSeries(output_path, "v")
T = FieldTimeSeries(output_path, "T")
e = FieldTimeSeries(output_path, "e")

z = znodes(u)

fig = Figure(resolution=(1200, 800))

ax_T = [Axis(fig[2, j], ylabel = "z (m)") for j = 1:5]
ax_u = [Axis(fig[3, j], ylabel = "z (m)") for j = 1:5]
ax_e = [Axis(fig[4, j], ylabel = "z (m)") for j = 1:5]

[Label(fig[1, j], replace(cases[j], "_" => " ")) for j=1:length(cases)]

Label(fig[2, 6], "Temperature", tellheight=false)
Label(fig[3, 6], "Velocity", tellheight=false)
Label(fig[4, 6], "Turbulent \n Kinetic \n Energy \n", tellheight=false)

Nt = length(u.times)
slider = Slider(fig[0, :], range=1:Nt, horizontal=true, startvalue=1)
n = slider.value

for j = 1:length(observations)
    obs = observations[j]
    u_truth = @lift interior(obs.field_time_serieses.u[$n])[1, 1, :]
    v_truth = @lift interior(obs.field_time_serieses.v[$n])[1, 1, :]
    T_truth = @lift interior(obs.field_time_serieses.T[$n])[1, 1, :]
    e_truth = @lift interior(obs.field_time_serieses.e[$n])[1, 1, :]

    lines!(ax_T[j], T_truth, z, color=(:gray23, 0.4), linewidth=6,   label="LES")
    lines!(ax_u[j], u_truth, z, color=(:gray23, 0.4), linewidth=6,   label="u, LES")
    lines!(ax_u[j], v_truth, z, color=(:gray23, 0.4), linewidth=2, label="v, LES")
    lines!(ax_e[j], e_truth, z, color=(:gray23, 0.4), linewidth=6,   label="LES")

    i = 1
    #for i = 1:Nensemble
        u★ = @lift interior(u[$n])[i, j, :]
        v★ = @lift interior(v[$n])[i, j, :]
        T★ = @lift interior(T[$n])[i, j, :]
        e★ = @lift interior(e[$n])[i, j, :]

        lines!(ax_T[j], T★, z, color=:purple1, linewidth=2,   label="$i CATKE")
        lines!(ax_u[j], u★, z, color=:purple1, linewidth=2,   label="u, $i CATKE")
        lines!(ax_u[j], v★, z, color=:purple1, linewidth=1, label="v, $i CATKE")
        lines!(ax_e[j], e★, z, color=:purple1, linewidth=2,   label="$i CATKE")
   #end
end

for j = 1:length(observations)

    if j == 1
        hidespines!(ax_T[j], :t, :r, :b)
        hidespines!(ax_u[j], :t, :r, :b)
        hidespines!(ax_e[j], :t, :r, :b)
        hidexdecorations!(ax_T[j])
        hidexdecorations!(ax_u[j])
        hidexdecorations!(ax_e[j])
        ax_T[j].ygridvisible[] = false
        ax_u[j].ygridvisible[] = false
        ax_e[j].ygridvisible[] = false
    else
        hidedecorations!(ax_T[j])
        hidedecorations!(ax_u[j])
        hidedecorations!(ax_e[j])
        hidespines!(ax_T[j])
        hidespines!(ax_u[j])
        hidespines!(ax_e[j])
    end

    xlims!(ax_T[j], 19.6, 20.05)
    xlims!(ax_u[j], -0.6, 0.6)
    xlims!(ax_e[j], -0.01, 0.01)
    #axislegend(ax_T[j], position=:rb)
    #axislegend(ax_u[j], position=:rb)
    #axislegend(ax_e[j], position=:rb)
end

display(fig)

max_stars_repo_name: glwagner/LocalOceanClosureCalibration


text:
immutable SVM <: RegERM
    X::Matrix               # n x m matrix of n m-dimensional training examples
    y::Vector               # 1 x n vector with training classes
    n::Int                  # number of training examples
    m::Int                  # number of features
    kernel::Symbol          # kernel function
    regression_type::Symbol # ordinal, binomial, multinomial
    params::Hyperparameters # hyperparameters (e.g. λ)
end

function SVM(X::Matrix, y::Vector; kernel::Symbol=:linear, λ::Float64=0.1)
    check_arguments(X, y, :binomial)
    SVM(X, y, size(X)..., kernel, :binomial, RegularizationParameters(λ))
end

function SVM(X::Matrix, y::Vector, params::Hyperparameters; kernel::Symbol=:linear)
    check_hyperparameters(params)
    check_arguments(X, y, :binomial)
    SVM(X, y, size(X)..., kernel, :binomial, params)
end

methodname(::SVM) = "Support Vector Machine"
loss(::SVM) = HingeLoss()
regularizer(SVM::SVM, w::Vector) = L2reg(w, SVM.params.λ)


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

include("libgit2.jl")
include("online.jl")

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
module Archibald
    using Distributions, StatsBase, LinearAlgebra
    using DataFrames, DataFramesMeta

    export MarkovChain, MarkovChainOverSample
    export MetropolisHastings!, MetropolisHastings
    export tidybayes

    abstract type Markovian end
    """
    We create a MarkovChain object with a vector of parameters, a data vector, and
    parametric type target_pdf which will be some sort of function. Using parametric
    types is much more efficient than creating a field with abstract type Function
    here.
    """
    struct MarkovChain{F<:Function} <: Markovian
        parameters::Vector
        data::Vector
        target_pdf::F
    end

    struct MarkovChainOverSample{F<:Function} <: Markovian
        single_sample_parameters::Vector
        S::Int64
        oversample_parameters::Matrix
        data::Vector
        target_pdf::F
    end
    include("samplers.jl")
    include("utils.jl")
end

max_stars_repo_name: EdJeeOnGitHub/Archibald

max_issues_repo_name: EdJeeOnGitHub/Archibald

max_forks_repo_name: EdJeeOnGitHub/Archibald


text:
module PlasmaParameters

using Unitful, PhysicalConstants.CODATA2018

export electron_gyrofrequency, ion_gyrofrequency,
    electron_plasma_frequency, ion_plasma_frequency,
    electron_trapping_rate, ion_trapping_rate,
    electron_collision_rate, ion_collision_rate,
    electron_debroglie_length, classical_minimum_approach_dist,
    electron_gryoradius, ion_gyroradius,
    debye_length, magnetic_debye_length,
    electron_thermal_velocity, ion_thermal_velocity,
    ion_sound_velocity, alfven_velocity

include("fundamental_plasma_params.jl")

end


text:
abstract type AbstractBlock{T<:Real} <: AbstractMatrix{T} end

const ZIGZAG = @SMatrix [
     1  2  6  7 15 16 28 29
     3  5  8 14 17 27 30 43
     4  9 13 18 26 31 42 44
    10 12 19 25 32 41 45 54
    11 20 24 33 40 46 53 55
    21 23 34 39 47 52 56 61
    22 35 38 48 51 57 60 62
    36 37 49 50 58 59 63 64
]

const ZIGZAG_INDICES = @SVector [findfirst(isequal(i), ZIGZAG) for i = 1:64]

Base.CartesianIndices(b::AbstractBlock) = ZIGZAG_INDICES

Base.getindex(b::AbstractBlock, y::Integer, x::Integer) = getindex(b, ZIGZAG[y, x])
Base.size(b::AbstractBlock) = (8, 8)
Base.length(b::AbstractBlock) = 64

abstract type Block{T<:Real} <: AbstractBlock{T} end

struct DCTBlock{T<:Real} <: Block{T}
    coeffs::SVector{64,T}
end

Base.getindex(b::DCTBlock, i::Integer) = b.coeffs[i]


text:
""" `pp = QuadraticProgram(A, b, C, d, z=zeros(..), P=I; kwargs...)`
    Type for problem
    `min_x 1/2xᵀPx + zᵀx, s.t Ax=b, Cx≤d`
    where `P` is positive definite.

    Keyword arguments:
    `semidefinite=true` Indicates that the dual might be semidefinite, which enables iterative refinement
    `ϵ = sqrt(eps(T))`  Relaxation used in iterative refinement (H+ϵI)
    `smartstart=true`   Enables a smart guess on initial active set
    `scaling=true`      Scales the constraints to have unit row norm

    Stores: `A,b,C,d,z,P,PF,sol,boxQP` and some temporary variables.
    `sol` is the solution after calling `solve!(qp)`,
    `PF` is a factorization of `P`,
    `boxQP::BoxConstrainedQP` representes the dual problem to be solved, with quadratic cost
    `G=[A*P⁻¹*A' A*P⁻¹*C'; C*P⁻¹*A' C*P⁻¹*C']`,
    linear cost `c = [A*z + b; C*z+d]`,
    and constraints `0≤xᵢ ∀i>size(A,1)`.
"""
struct QuadraticProgram{T, GT<:AbstractCholeskySpecial{T}, VT<:AbstractVector{T}, MT<:AbstractMatrix{T}, PT, PFT}
    A::MT
    C::MT
    b::VT
    d::VT
    z::VT
    P::PT
    PF::PFT
    tmp3::VT
    sol::VT
    scaling::Bool
    scaleA::VT
    scaleC::VT
    boxQP::BoxConstrainedQP{T,GT,VT}
end

function QuadraticProgram(A::MT, b::VT, C::MT, d::VT, z::VT=fill(zero(T), size(A,2)), P=I; semidefinite=true, ϵ = sqrt(eps(T)), smartstart=true, scaling=true) where {T, VT<:AbstractVector{T}, MT<:AbstractMatrix{T}}
    m = size(A,1)
    n = size(C,1)
    # Build matrix a bit more efficient
    #GQ = [A*P⁻¹*A' A*P⁻¹*C'; C*P⁻¹*A' C*P⁻¹*C']
    A = copy(A)
    C = copy(C)
    b = copy(b)
    d = copy(d)
    if scaling
        scaleA = vec(sqrt.(sum(abs2, A, dims=2)))
        scaleC = vec(sqrt.(sum(abs2, C, dims=2)))
        A .= A ./ scaleA
        C .= C ./ scaleC
        b .= b ./ scaleA
        d .= d ./ scaleC
    else
        # Empty vectors, make sure not to use them
        scaleA = VT(undef,0)
        scaleC = VT(undef,0)
    end

    GQ = similar(A, m+n, m+n)
    GQ11 = view(GQ, 1:m,1:m)
    GQ12 = view(GQ, 1:m, (m+1):(m+n))
    GQ22 = view(GQ, (m+1):(m+n), (m+1):(m+n))

    # Dispatch on type of P to get efficient solve of
    # PF = factorize(P), PiAt = PF\A', PiCt = PF\C'
    PF, PiAt, PiCt = factorizesolve(P,A,C)

    mul!(GQ11,  A, PiAt)
    mul!(GQ12, A, PiCt)

    mul!(GQ22, C, PiCt)
    GQ[(m+1):(m+n), 1:m] .= GQ12'

    # Dual linear cost
    dualc = similar(b, m+n)     # Set in update!
    duald = fill(zero(T), n)       # Always zero

    # Dual problem is simple Semidefinite Quadratic, with 0≤xᵢ ∀i>m
    boxQP = BoxConstrainedQP(Hermitian(GQ), dualc, duald; semidefinite=semidefinite, ϵ = ϵ)

    sol = similar(z)
    tmp3 = similar(z)
    QP = QuadraticProgram{T,typeof(boxQP.G), VT, MT, typeof(P), typeof(PF)}(
        A, C, b, d, z,
        P, PF,
        tmp3, sol, scaling, scaleA, scaleC, boxQP)

    # TODO smartstart already at first factorization
    # Set dual linear cost
    update!(QP)
    if smartstart
        run_smartstart(boxQP)
    end
    return QP
end

# Get status from boxQP
function Base.getproperty(QP::QuadraticProgram, s::Symbol)
    if s == :status
        return getproperty(QP.boxQP, s)
    end

    return getfield(QP, s)
end

# General safe implementation, Sparse LDLT doesn't support PF\A'
function factorizesolve(P, A, C)
    PF = factorize(P)
    PiAt = PF\Matrix(A')
    PiCt = PF\Matrix(C')
    return PF, PiAt, PiCt
end

function factorizesolve(P::UniformScaling{Bool}, A, C)
    if P.λ != true
        error("P not positive definite")
    end
    return P, A', C'
end

function factorizesolve(P::UniformScaling, A, C)
    # Send back Adjoint instread?
    PF = P
    PiAt = Matrix(A')
    PiCt = Matrix(C')
    PiAt ./= P.λ
    PiCt ./= P.λ
    return P, PiAt, PiCt
end

function factorizesolve(P::Matrix, A, C)
    PF = factorize(P)
    # Should work for dense case
    PiAt = PF\A'
    PiCt = PF\C'
    return PF, PiAt, PiCt
end

function solve!(QP::QuadraticProgram)
    reset!(QP.status)
    # Solve dual problem
    xk = solve!(QP.boxQP)

    m = size(QP.A,1)
    n = size(QP.C,1)
    mul!(QP.sol, QP.C', view(xk, (m+1):(m+n)))

    mul!(QP.tmp3, QP.A', view(xk, 1:m))
    # sol = -P\(z + A'μλ[1:m] + C'μλ[(m+1):(m+n)])
    QP.sol .= .- QP.z .- QP.tmp3 .- QP.sol
    # TODO Efficiency
    sol = QP.PF\QP.sol
    QP.sol .= sol

    #1/2xᵀPx - zᵀx
    # TODO Efficiency
    tmp = QP.P*sol
    return sol, dot(tmp,sol)/2 + dot(sol, QP.z)
end

function update!(QP::QuadraticProgram; b=nothing, d=nothing, z=QP.z)
    # TODO ?

    # If updating b or z, do scaling before saving
    if  b !== nothing
        if QP.scaling
            QP.b .= b ./ QP.scaleA
        else
            QP.b .= b
        end
    end
    if d !== nothing
        if QP.scaling
            QP.d .= d ./ QP.scaleC
        else
            QP.d .= d
        end
    end

    # We need to update c
    m = size(QP.A,1)
    n = size(QP.C,1)
    c1 = view(QP.boxQP.c, 1:m)
    c2 = view(QP.boxQP.c, (m+1):(m+n))

    # TODO efficiency
    Piz = QP.PF\z
    mul!(c1, QP.A, Piz)
    c1 .= c1 .+ QP.b # c1 = A*z + b
    mul!(c2, QP.C, Piz)
    c2 .= c2 .+ QP.d # c2 = C*z + d

    QP.z .= z

    return
end

max_stars_repo_name: mfalt/QPDAS

max_issues_repo_name: mfalt/QPDAS

max_forks_repo_name: mfalt/QPDAS


text:
function star_attack(board::UInt64, ui::UInt64)
    free_squares, edge_squares = orthogonal_attack(board, ui)
    free_squares, edge_squares  = cross_attack(free_squares, edge_squares,
        board, ui)

    return free_squares, edge_squares
end

function star_attack(free_squares::Array{UInt64,1},
    edge_squares::Array{UInt64,1}, occ::UInt64, ui::UInt64)
    free_squares, edge_squares = orthogonal_attack(free_squares, edge_squares,
        occ, ui)
    free_squares, edge_squares = cross_attack(free_squares, edge_squares, occ,
        ui)

    return free_squares, edge_squares
end


text:
@testset "KW" begin
    @test @kw(a=2) === FastKeywords.KW{:a, Int}(2)
    @test (@kw a=2) === FastKeywords.KW{:a, Int}(2)

    @test FastKeywords.name(@kw(a=2)) === :a
    @test get(@kw(a=2)) === 2
end


text:
module Noop

export @stopwatch

macro stopwatch(_, code)
    esc(code)
end

macro stopwatch(code)
    esc(code)
end

end


text:
using OrdinaryDiffEq
prob = ODEProblem((u,p,t) -> -u,BigFloat(1.0),(0.0,1.0))
solve(prob,Tsit5())
solve(prob,KenCarp4())


text:
using FileIO

# ***** Delimiters *****
@enum Delimiter space=1 comma=2
DELIMITER = Dict{Delimiter, Char}(
    space::Delimiter => ' ',
    comma::Delimiter => ',')


text:
function get_configuration(obj::ItrfVersionProvider, arg0::JString, arg1::jint)
    return jcall(obj, "getConfiguration", ITRFVersionLoader_ITRFVersionConfiguration, (JString, jint), arg0, arg1)
end


text:
using Flux: TrackedVector, TrackedMatrix
using Flux.Tracker: TrackedVecOrMat, TrackedReal
using Dates, Formatting
import Base: convert

macro noopwhen(condition, expression)
    quote
        if !($condition)
            $expression
        end
    end |> esc
end

function countmap(x::Vector{T}; d::Int=-1) where T <: Signed
    # approx. 10x faster than StatsBase for small arrays of ints.
    # Note that it always enumerates each possible integer up to max.
    d = d > 0 ? d : maximum(x)
    out = zeros(Int64, d)
    for cx in x
        out[cx] += 1
    end
    return out
end

function countmap_labelled(x::Vector{T}) where T <: Signed
    # as countmap, but not array of contiguous labels: a labelled
    # vector is also returned, and any key of zero count is removed.
    sort!(x)
    lbls = unique(x)
    d = length(lbls)
    out = zeros(Int64, d)
    px = x[1]
    k = 1
    for cx in x
        if cx != px
            px = cx
            k += 1
        end
        out[k] += 1
    end
    return out, lbls
end

invert_index(x::Vector{T}) where T <: Signed = sortperm(x)

# useful for pipes
dropdim1(x::Union{Number, TrackedReal}) = x
dropdim1(x::Union{VecOrMat, TrackedVecOrMat}) = dropdims(x, dims=1)

dropdim2(x::Union{Matrix, TrackedMatrix}) = dropdims(x, dims=2)
dropdim2(x::Union{Vector, TrackedVector}) = x
dropdim2(x::Union{Number, TrackedReal}) = x

# converting NamedTuple -> Dict fails if multiple types in the NamedTuple.
# The *collect* in the second zip arg here is the work-around.
convert(::Type{Dict}, nt::NamedTuple) = Dict(k => v for (k,v) in zip(keys(nt), collect(values(nt))))

# issue #29560 mkborregaard solution to repelem/rep like behaviour for vectors.
# Will likely be obselete when PR accepted.
#
# NOTICE THAT THE FIRST FUNCTION WILL RETURN A 2D MATRIX EVEN IF MATRIX GIVEN.
# (Hence my vector function below.)
function Base.repeat(a::AbstractMatrix, m::AbstractVector{<:Integer}, n::Integer=1)
    o, p = size(a,1), size(a,2)
    length(m) == o || throw(ArgumentError("m must be an AbstractVector{<:Integer} of same dimension as a"))
    b = similar(a, sum(m), p*n)
    for j=1:n
        d = (j-1)*p+1
        R = d:d+p-1
        c = 0
        for i ∈ eachindex(m)
            r = m[i]
            v = a[i,:]
            while r > 0
                b[c+=1, R] = v
                r -= 1
            end
        end
    end
    return b
end

function Base.repeat(a::AbstractVector, m::AbstractVector{<:Integer})
    o, p = size(a,1), size(a,2)
    length(m) == o || throw(ArgumentError("m must be an AbstractVector{<:Integer} of same dimension as a"))
    b = similar(a, sum(m))
    c = 0
    for i ∈ eachindex(m)
        r = m[i]
        v = a[i]
        while r > 0
            b[c+=1] = v
            r -= 1
        end
    end
    return b
end

function construct_unique_filename(filestem; path="./", date_fmt="yyyy_mm_dd", ext="")
    fnm_base = filestem * Dates.format(now(), date_fmt)
    fnm_uid = 0
    _construct_fnm(uid) = joinpath(path, fnm_base * "_" * format("{:03d}", uid) * ext)
    while isfile(_construct_fnm(fnm_uid))
        fnm_uid += 1
    end
    return _construct_fnm(fnm_uid)
end

# BSON.bson("tmptest.bson",
#     par=Float32[0.1 2.345 6.78 9.0; 3.14 15.92 65.35 89.7; 2.71 8.2 8.18 2.8459; 4.66 9.2 0.16 0.91])
# bson_raw = BSON.parse("tmptest.bson")[:par]   # ok... not raw exactly

function parse_serialized_array(bson_raw)
    type_string = bson_raw[:type][:name][end]
    T = eval(Meta.parse(type_string))

    @assert isconcretetype(T) "The strategy here assumes the array is a concrete type."

    io = IOBuffer()
    write(io, bson_raw[:data])
    seekstart(io)

    flat = T[]
    while true
        try
            element = read(io, T)
            push!(flat, element)
        catch e
            if e isa EOFError
                break
            end
            throw(e)
        end
    end
    reshape(flat, bson_raw[:size]...)
end


text:
using Revise; using RHEOS; using BenchmarkTools; using FunctionWrappers:FunctionWrapper;

RheoFloat = Float64

struct TT
   G::FunctionWrapper{RheoFloat, Tuple{RheoFloat, Vector{RheoFloat}}}
   J::FunctionWrapper{RheoFloat, Tuple{RheoFloat, Vector{RheoFloat}}}

   Ga::FunctionWrapper{Vector{RheoFloat}, Tuple{Vector{RheoFloat}, Vector{RheoFloat}}}
   Ja::FunctionWrapper{Vector{RheoFloat}, Tuple{Vector{RheoFloat}, Vector{RheoFloat}}}
end

function test1()
   # test base case
   p=(:a1, :a2)
   unpack_expr = Meta.parse(string(join(string.(p), ","), ", = a"))
   G = :(a1*t + a2)
   J = :(a2*t + a1)
   @eval return(TT(
      ((t, a) -> ($unpack_expr; $G)) |> FunctionWrapper{RheoFloat, Tuple{RheoFloat, Vector{RheoFloat}}},
      ((t, a) -> ($unpack_expr; $J)) |> FunctionWrapper{RheoFloat, Tuple{RheoFloat, Vector{RheoFloat}}},

      ((ta, a) -> ($unpack_expr; [$G for t in ta])) |> FunctionWrapper{Vector{RheoFloat}, Tuple{Vector{RheoFloat}, Vector{RheoFloat}}},
      ((ta, a) -> ($unpack_expr; [$J for t in ta])) |> FunctionWrapper{Vector{RheoFloat}, Tuple{Vector{RheoFloat}, Vector{RheoFloat}}},
      ))
end

function test2()
   # check if quote ... end affects speed
   p=(:a1, :a2)
   unpack_expr = Meta.parse(string(join(string.(p), ","), ", = a"))
   G = quote a1*t + a2 end
   J = quote a2*t + a1 end
   @eval return(TT(
      ((t, a) -> ($unpack_expr; $G)) |> FunctionWrapper{RheoFloat, Tuple{RheoFloat, Vector{RheoFloat}}},
      ((t, a) -> ($unpack_expr; $J)) |> FunctionWrapper{RheoFloat, Tuple{RheoFloat, Vector{RheoFloat}}},

      ((ta, a) -> ($unpack_expr; [$G for t in ta])) |> FunctionWrapper{Vector{RheoFloat}, Tuple{Vector{RheoFloat}, Vector{RheoFloat}}},
      ((ta, a) -> ($unpack_expr; [$J for t in ta])) |> FunctionWrapper{Vector{RheoFloat}, Tuple{Vector{RheoFloat}, Vector{RheoFloat}}},
      ))
end

function test3()
   # check if begin ... end affects speed
   p=(:a1, :a2)
   unpack_expr = Meta.parse(string(join(string.(p), ","), ", = a"))
   G = :(a1*t + a2)
   J = :(a2*t + a1)
   @eval return(TT(
      ((t, a) -> begin $unpack_expr; $G end) |> FunctionWrapper{RheoFloat, Tuple{RheoFloat, Vector{RheoFloat}}},
      ((t, a) -> begin $unpack_expr; $J end) |> FunctionWrapper{RheoFloat, Tuple{RheoFloat, Vector{RheoFloat}}},

      ((ta, a) -> begin $unpack_expr; [$G for t in ta] end) |> FunctionWrapper{Vector{RheoFloat}, Tuple{Vector{RheoFloat}, Vector{RheoFloat}}},
      ((ta, a) -> begin $unpack_expr; [$J for t in ta] end) |> FunctionWrapper{Vector{RheoFloat}, Tuple{Vector{RheoFloat}, Vector{RheoFloat}}},
      ))
end

function test4()
   # check if . broadcast has any affect on speed
   p=(:a1, :a2)
   unpack_expr = Meta.parse(string(join(string.(p), ","), ", = a"))
   G = :(a1*t + a2)
   J = :(a2*t + a1)
   @eval return(TT(
      ((t, a) -> ($unpack_expr; $G)) |> FunctionWrapper{RheoFloat, Tuple{RheoFloat, Vector{RheoFloat}}},
      ((t, a) -> ($unpack_expr; $J)) |> FunctionWrapper{RheoFloat, Tuple{RheoFloat, Vector{RheoFloat}}},

      ((ta, a) -> ($unpack_expr; [$G for t in ta])) |> FunctionWrapper{Vector{RheoFloat}, Tuple{Vector{RheoFloat}, Vector{RheoFloat}}},
      ((ta, a) -> ($unpack_expr; [$J for t in ta])) |> FunctionWrapper{Vector{RheoFloat}, Tuple{Vector{RheoFloat}, Vector{RheoFloat}}},
      ))
end

params = [1.0, 1.0];
testvector = Vector(0.0:0.001:100.0)

println("Testing broadcast over single-entry function wrapper\n=======================================================")
t1 = test1()
println("Array argument function wrapper")
@btime t1.Ga($testvector, $params)
# 76.687 μs // 781 KiB
println("Broadcast over single-entry function wrapper")
paramsTuple = ([1.0, 1.0],);
@btime t1.G.($testvector, $paramsTuple)
# 798 μs // 781 KiB
println("Test complete.\n\n\n")
# Array argument function wrapper always ~1 order of magnitude faster.

println("Testing quote ... end vs :( ... ) syntax\n=======================================================")
t2 = test2()
println("quote ... end syntax")
@btime t2.Ga($testvector, $params)
# 77 μs // 781 KiB
println(":( ) syntax")
paramsTuple = ([1.0, 1.0],);
@btime t1.Ga($testvector, $params)
# 77 μs // 781 KiB
println("Test complete.\n\n\n")
# Alternating very slight difference, negligible. Stick with quote ... end for now.

println("Testing begin ... end vs ( ... ) syntax\n=======================================================")
t3 = test3()
println("begin ... end syntax")
@btime t3.Ga($testvector, $params)
# 76 μs // 781 KiB
println("( ) syntax")
paramsTuple = ([1.0, 1.0],);
@btime t1.Ga($testvector, $params)
# 76 μs // 781 KiB
println("Test complete.\n\n\n")
# Alternating very slight difference, negligible. Stick with begin ... end for now.

println("Testing . broadcast vs array comprehension syntax\n=======================================================")
t4 = test4()
println("dot broadcast")
@btime t4.Ga($testvector, $params)
#  μs //  KiB
println("array comprehension")
paramsTuple = ([1.0, 1.0],);
@btime t1.Ga($testvector, $params)
#  μs //  KiB
println("Test complete.\n\n\n")
#

max_stars_repo_name: JuliaRheology/RheoBenchAndTest

max_issues_repo_name: JuliaRheology/RheoBenchAndTest

max_forks_repo_name: JuliaRheology/RheoBenchAndTest


text:
module TestTimestepArrays

using Mimi
using Test

import Mimi:
    FixedTimestep, VariableTimestep, TimestepVector, TimestepMatrix, TimestepArray, next_timestep, hasvalue,
    isuniform, first_period, last_period, first_and_step

function reset_time_val(arr, vals::Array{T, 1}) where T
    arr[:] = vals
end

function reset_time_val(arr, vals::Array{T, 2}) where T
    arr[:,:] = vals
end

function reset_time_val(arr, vals::Array{T, 3}) where T
    arr[:,:,:] = vals
end

## quick check of isuniform
@test isuniform([]) == false
@test isuniform([1]) == true
@test isuniform([1,2,3]) == true
@test isuniform([1,2,3,5]) == false
@test isuniform(1) == true

@test first_and_step([2,3,4]) == (2,1)
@test first_and_step([1:2:10...]) == (1,2)

# shorthand for stuff used a lot
idx1 = TimestepIndex(1)
idx2 = TimestepIndex(2)
idx3 = TimestepIndex(3)
idx4 = TimestepIndex(4)

#------------------------------------------------------------------------------
# 1. Test TimestepVector - Fixed Timestep
#------------------------------------------------------------------------------

#1a.  test constructor, lastindex, and length (with both
# matching years and mismatched years)

x = TimestepVector{FixedTimestep{2000, 1}, Int}([9, 10, 11, 12])
@test length(x) == 4
@test lastindex(x) == 4

time_dim_val = [9, 10, 11, 12]
temp_dim_val = [100, 101, 102, 103]

#1b.  test hasvalue, getindex, and setindex! (with both matching years and
# mismatched years)

# Using a colon in the time dimension
@test x[:] == time_dim_val
x[:] = temp_dim_val
@test x[:] == temp_dim_val
reset_time_val(x, time_dim_val)

# TimestepValue and TimestepIndex Indexing
@test x[idx1] == time_dim_val[1]
@test x[idx1 + 1] == time_dim_val[2]
@test x[idx4] == time_dim_val[4]
@test_throws ErrorException x[TimestepIndex(5)]

x[idx1] = temp_dim_val[1]
@test x[idx1] == temp_dim_val[1]
reset_time_val(x, time_dim_val)

@test x[TimestepValue(2000)] == time_dim_val[1]
@test x[TimestepValue(2000; offset = 1)] == time_dim_val[2]
@test x[TimestepValue(2000) + 1] == time_dim_val[2]
@test_throws ErrorException x[TimestepValue(2005)]
@test_throws ErrorException x[TimestepValue(2004)+1]

x[TimestepValue(2000)] = temp_dim_val[1]
@test x[TimestepValue(2000)] == temp_dim_val[1]
reset_time_val(x, time_dim_val)

# AbstractTimestep Indexing
t = FixedTimestep{2001, 1, 3000}(1)

@test hasvalue(x, t)
@test !hasvalue(x, FixedTimestep{2000, 1, 2012}(10))
@test x[t] == 10

t2 = next_timestep(t)
@test x[t2] == time_dim_val[3]
x[t2] = temp_dim_val[3]
@test x[t2] == temp_dim_val[3]
reset_time_val(x, time_dim_val)

t3 = FixedTimestep{2000, 1, 2003}(1)
@test x[t3] == time_dim_val[1]
x[t3] = temp_dim_val[1]
@test x[t3] == temp_dim_val[1]
reset_time_val(x, time_dim_val)

# Deprecated int indexing should still run
@test x[3] == time_dim_val[3]
x[TimestepIndex(3)] = temp_dim_val[3]
@test x[TimestepIndex(3)] == temp_dim_val[3]
reset_time_val(x, time_dim_val)

#------------------------------------------------------------------------------
# 2. Test TimestepVector - Variable Timestep
#------------------------------------------------------------------------------

years = (2000, 2005, 2015, 2025)
x = TimestepVector{VariableTimestep{years}, Int}([9, 10, 11, 12])

time_dim_val = [9, 10, 11, 12]
temp_dim_val = [100, 101, 102, 103]

#2a.  test hasvalue, getindex, and setindex! (with both matching years and
# mismatched years)

# Using a colon in the time dimension
@test x[:] == time_dim_val
x[:] = temp_dim_val
@test x[:] == temp_dim_val
reset_time_val(x, time_dim_val)

# TimestepValue and TimestepIndex Indexing
@test x[idx1] == time_dim_val[1]
@test x[idx1 + 1] == time_dim_val[2]
@test x[idx4] == time_dim_val[4]
@test_throws ErrorException x[TimestepIndex(5)]

x[idx1] = temp_dim_val[1]
@test x[idx1] == temp_dim_val[1]
reset_time_val(x, time_dim_val)

@test x[TimestepValue(2000)] == time_dim_val[1]
@test x[TimestepValue(2000; offset = 1)] ==  time_dim_val[2]
@test x[TimestepValue(2015)] ==  time_dim_val[3]
@test x[TimestepValue(2015; offset = 1)] ==  time_dim_val[4]
@test x[TimestepValue(2000) + 1] ==  time_dim_val[2]
@test_throws ErrorException x[TimestepValue(2014)]
@test_throws ErrorException x[TimestepValue(2025)+1]

x[TimestepValue(2015)] = temp_dim_val[3]
@test x[TimestepValue(2015)] == temp_dim_val[3]
reset_time_val(x, time_dim_val)

# AbstractTimestep Indexing
y2 = Tuple([2005:5:2010; 2015:10:3000])
t = VariableTimestep{y2}()

@test hasvalue(x, t)
@test !hasvalue(x, VariableTimestep{years}(time_dim_val[2]))
@test x[t] == time_dim_val[2]

t2 =  next_timestep(t)
@test x[t2] == time_dim_val[3]
x[t2] = temp_dim_val[3]
@test x[t2] == temp_dim_val[3]
reset_time_val(x, time_dim_val)

t3 = VariableTimestep{years}()
@test x[t3] == time_dim_val[1]
x[t3] = temp_dim_val[1]
@test x[t3] == temp_dim_val[1]
reset_time_val(x, time_dim_val)

@test x[TimestepIndex(3)] == time_dim_val[3]
x[TimestepIndex(3)] = temp_dim_val[3]
@test x[TimestepIndex(3)] == temp_dim_val[3]
reset_time_val(x, time_dim_val)

#------------------------------------------------------------------------------
# 3. Test TimestepMatrix - Fixed Timestep
#------------------------------------------------------------------------------

for ti = 1:2

    #3a.  test constructor (with both matching years
    # and mismatched years)

    y = TimestepMatrix{FixedTimestep{2000, 1}, Int, ti}(collect(reshape(1:8, 4, 2)))
    z = TimestepMatrix{FixedTimestep{2000, 2}, Int, ti}(collect(reshape(1:8, 4, 2)))

    time_dim_val = collect(reshape(1:8, 4, 2))
    temp_dim_val = collect(reshape(100:107, 4, 2))

    #3b.  test hasvalue, getindex, and setindex! (with both matching years and
    # mismatched years)

    # Using a colon in the time dimension
    y[:,:] = temp_dim_val
    @test y[:,:] == temp_dim_val
    y[:,:] = time_dim_val # reset
    if ti == 1
        @test y[:,1] == time_dim_val[:,1]
        y[:, 1] = temp_dim_val[:,1]
        @test y[:, 1] == temp_dim_val[:,1]
    else
        @test y[1,:] == time_dim_val[1,:]
        y[1,:] = temp_dim_val[1,:]
        @test y[1,:] == temp_dim_val[1,:]
    end
    reset_time_val(y, time_dim_val)

    # TimestepValue and TimestepIndex Indexing
    if ti == 1
        @test y[idx1, 1] == time_dim_val[1,1]
        @test y[idx1, 2] == time_dim_val[1,2]
        @test y[idx1 + 1, 1] == time_dim_val[2,1]
        @test y[idx4, 2] == time_dim_val[4,2]
        @test_throws ErrorException y[TimestepIndex(5), 2]

        y[idx1, 1] = temp_dim_val[1]
        @test y[idx1, 1] == temp_dim_val[1]
        reset_time_val(y, time_dim_val)

        @test y[TimestepValue(2000), 1] == time_dim_val[1, 1]
        @test y[TimestepValue(2000), 2] == time_dim_val[1, 2]
        @test y[TimestepValue(2001), :] == time_dim_val[2, :]
        @test y[TimestepValue(2000; offset = 1), 1] == time_dim_val[2,1]
        @test y[TimestepValue(2000) + 1, 1] == time_dim_val[2,1]
        @test_throws ErrorException y[TimestepValue(2005), 1]
        @test_throws ErrorException y[TimestepValue(2004)+1, 1]

        y[TimestepValue(2000), 1] = temp_dim_val[1]
        @test y[TimestepValue(2000), 1] == temp_dim_val[1]
        reset_time_val(y, time_dim_val)
    else
        @test y[1, idx1] == time_dim_val[1,1]
        @test y[2, idx1] == time_dim_val[2, 1]
        @test y[1, idx1 + 1] == time_dim_val[1, 2]
        @test y[2, TimestepIndex(2)] == time_dim_val[2,2]
        @test_throws ErrorException y[2, TimestepIndex(3)]

        y[1, idx1] = temp_dim_val[1]
        @test y[1, idx1] == temp_dim_val[1]
        reset_time_val(y, time_dim_val)

        @test y[1, TimestepValue(2000)] == time_dim_val[1,1]
        @test y[2, TimestepValue(2000)] == time_dim_val[2,1]
        @test y[:, TimestepValue(2001)] == time_dim_val[:,2]
        @test y[1, TimestepValue(2000; offset = 1)] == time_dim_val[1,2]
        @test y[1, TimestepValue(2000) + 1] == time_dim_val[1,2]
        @test_throws ErrorException y[1, TimestepValue(2003)]
        @test_throws ErrorException y[1, TimestepValue(2002)+1]

        y[1, TimestepValue(2000)] = temp_dim_val[1]
        @test y[1, TimestepValue(2000)] == temp_dim_val[1]
        reset_time_val(y, time_dim_val)
    end

    # AbstractTimestep Indexing
    t = FixedTimestep{2001, 1, 3000}(1)
    @test hasvalue(y, t, 1)
    @test !hasvalue(y, FixedTimestep{2000, 1, 3000}(10), 1)

    if ti == 1
        @test y[t,1] == time_dim_val[2,1]
        @test y[t,2] == time_dim_val[2,2]

        t2 = next_timestep(t)
        @test y[t2,1] == time_dim_val[3,1]
        @test y[t2,2] == time_dim_val[3,2]
        y[t2, 1] = temp_dim_val[3,1]
        @test y[t2, 1] == temp_dim_val[3,1]
        reset_time_val(y, time_dim_val)

        t3 = FixedTimestep{2000, 1, 2005}(1)
        @test y[t3, 1] == time_dim_val[1,1]
        @test y[t3, 2] == time_dim_val[1,2]
        y[t3, 1] = temp_dim_val[1,1]
        @test y[t3,1] == temp_dim_val[1,1]
        reset_time_val(y, time_dim_val)

        #3c.  interval wider than 1 using z from above
        t = FixedTimestep{1980, 2, 3000}(11)

        @test z[t,1] == time_dim_val[1,1]
        @test z[t,2] == time_dim_val[1,2]

        t2 = next_timestep(t)
        @test z[t2,1] == time_dim_val[2,1]
        @test z[t2,2] == time_dim_val[2,2]

    else
        @test y[1, t] == time_dim_val[1,2]
        @test y[2, t] == time_dim_val[2,2]

        t2 = FixedTimestep{2000, 1, 2005}(1)
        @test y[1, t2] == time_dim_val[1,1]
        @test y[2, t2] == time_dim_val[2,1]
        y[1, t2] = temp_dim_val[1,1]
        @test y[1, t2] == temp_dim_val[1,1]
        reset_time_val(y, time_dim_val)

        #3c.  interval wider than 1 using z from above
        t = FixedTimestep{1980, 2, 3000}(11)

        @test z[1, t] == time_dim_val[1,1]
        @test z[2, t] == time_dim_val[2,1]

        t2 = next_timestep(t)
        @test z[1, t2] == time_dim_val[1,2]
        @test z[2, t2] == time_dim_val[2,2]
    end

    @test y[idx1, 2] == time_dim_val[1,2]
    @test z[idx1, 2] == time_dim_val[1,2]
    y[idx1, 2] = temp_dim_val[1]
    z[idx1, 2] = temp_dim_val[1]
    @test y[idx1, 2] == z[idx1, 2] == temp_dim_val[1]

    reset_time_val(y, time_dim_val)
    reset_time_val(z, time_dim_val)
end

#------------------------------------------------------------------------------
# 4. Test TimestepMatrix - Variable Timestep
#------------------------------------------------------------------------------

for ti = 1:2

    #4a.  test constructor (with both matching years
    # and mismatched years)

    if ti == 1
        years = Tuple([2000:5:2005; 2015:10:2025])
    else
        years = (2000, 2005)
    end
    y = TimestepMatrix{VariableTimestep{years}, Int, ti}(collect(reshape(1:8, 4, 2)))

    time_dim_val = collect(reshape(1:8, 4, 2))
    temp_dim_val = collect(reshape(100:107, 4, 2))

    #4b.  test hasvalue, getindex, setindex!, and lastindex (with both matching years and
    # mismatched years)

    # Using a colon in the time dimension
    y[:,:] = temp_dim_val
    @test y[:,:] == temp_dim_val
    y[:,:] = time_dim_val # reset
    if ti == 1
        @test y[:,1] == time_dim_val[:,1]
        y[:, 1] = temp_dim_val[:,1]
        @test y[:, 1] == temp_dim_val[:,1]
    else
        @test y[1,:] == time_dim_val[1,:]
        y[1,:] = temp_dim_val[1,:]
        @test y[1,:] == temp_dim_val[1,:]
    end
    reset_time_val(y, time_dim_val)

    # TimestepValue and TimestepIndex Indexing
    if ti == 1
        @test y[idx1, 1] == time_dim_val[1,1]
        @test y[idx1, 2] == time_dim_val[1,2]
        @test y[idx1 + 1, 1] == time_dim_val[2,1]
        @test y[idx4, 2] == time_dim_val[4,2]
        @test_throws ErrorException y[TimestepIndex(5),2]

        y[idx1, 1] = temp_dim_val[1,1]
        @test y[idx1, 1] == temp_dim_val[1,1]
        reset_time_val(y, time_dim_val)

        @test y[TimestepValue(2000), 1] == time_dim_val[1,1]
        @test y[TimestepValue(2000), 2] == time_dim_val[1,2]
        @test y[TimestepValue(2000; offset = 1), 1] == time_dim_val[2,1]
        @test y[TimestepValue(2000) + 1, 1] == time_dim_val[2,1]
        @test y[TimestepValue(2015), 1] == time_dim_val[3,1]
        @test y[TimestepValue(2015) + 1, 2] == time_dim_val[4,2]
        @test_throws ErrorException y[TimestepValue(2006), 1]
        @test_throws ErrorException y[TimestepValue(2025)+1, 1]

        y[TimestepValue(2015), 1] = temp_dim_val[3,1]
        @test y[TimestepValue(2015), 1] == temp_dim_val[3,1]
        reset_time_val(y, time_dim_val)
    else
        @test y[1, idx1] == time_dim_val[1,1]
        @test y[2, idx1] == time_dim_val[2,1]
        @test y[1, idx1 + 1] == time_dim_val[1,2]
        @test y[2, TimestepIndex(2)] == time_dim_val[2,2]
        @test_throws ErrorException y[2, TimestepIndex(3)]

        y[1, idx1] = temp_dim_val[1,1]
        @test y[1, idx1] == temp_dim_val[1,1]
        reset_time_val(y, time_dim_val)

        @test y[1, TimestepValue(2000)] == time_dim_val[1,1]
        @test y[2, TimestepValue(2000)] == time_dim_val[2,1]
        @test y[1, TimestepValue(2000; offset = 1)] == time_dim_val[1,2]
        @test y[1, TimestepValue(2000) + 1] == time_dim_val[1,2]
        @test y[1, TimestepValue(2005)] == time_dim_val[1,2]
        @test_throws ErrorException y[1, TimestepValue(2006)]
        @test_throws ErrorException y[1, TimestepValue(2005)+1]

        y[1, TimestepValue(2005)] = temp_dim_val[1,2]
        @test y[1, TimestepValue(2005)] == temp_dim_val[1,2]
        reset_time_val(y, time_dim_val)
    end

    # AbstractTimestep Indexing
    t = VariableTimestep{Tuple([2005:5:2010; 2015:10:3000])}()
    if ti == 1
        @test hasvalue(y, t, time_dim_val[1])
        @test !hasvalue(y, VariableTimestep{years}(10))
        @test y[t,1] == time_dim_val[2,1]
        @test y[t,2] == time_dim_val[2,2]

        t2 = next_timestep(t)
        @test y[t2,1] == time_dim_val[3,1]
        @test y[t2,2] == time_dim_val[3,2]
        y[t2, 1] = temp_dim_val[3,1]
        @test y[t2, 1] == temp_dim_val[3,1]
        reset_time_val(y, time_dim_val)

        t3 = VariableTimestep{years}()
        @test y[t3, 1] == time_dim_val[1,1]
        @test y[t3, 2] == time_dim_val[1,2]
        y[t3, 1] = temp_dim_val[1,1]
        @test y[t3,1] == temp_dim_val[1,1]
        reset_time_val(y, time_dim_val)

    else
        @test hasvalue(y, t, time_dim_val[1])
        @test !hasvalue(y, VariableTimestep{years}(10))
        @test y[1,t] == time_dim_val[1,2]
        @test y[2,t] == time_dim_val[2,2]

        t3 = VariableTimestep{years}()
        @test y[1, t3] == time_dim_val[1,1]
        @test y[2, t3] == time_dim_val[2,1]
        y[1,t3] = temp_dim_val[1,1]
        @test y[1,t3] == temp_dim_val[1,1]
        reset_time_val(y, time_dim_val)
    end

    @test y[idx1,2] == time_dim_val[1,2]
    y[idx1,2] = temp_dim_val[1,2]
    @test y[idx1,2] == temp_dim_val[1,2]
    reset_time_val(y, time_dim_val)
end

#------------------------------------------------------------------------------
# 5. Test TimestepArray methods (3 dimensional)
#------------------------------------------------------------------------------

for ti = 1:2

    years = Tuple([2000:5:2005; 2015:10:2025])
    arr_fixed = TimestepArray{FixedTimestep{2000, 5}, Int, 3, ti}(collect(reshape(1:64, 4, 4, 4)))
    arr_variable = TimestepArray{VariableTimestep{years}, Int, 3, ti}(collect(reshape(1:64, 4, 4, 4)))

    time_dim_val = collect(reshape(1:64, 4, 4, 4))
    temp_dim_val = collect(reshape(100:163, 4, 4, 4))

    # Using a colon in the time dimension
    if ti == 1
        @test arr_fixed[:,1,1] == arr_variable[:,1,1] == time_dim_val[:,1,1]
        @test arr_fixed[:,2,3] == arr_variable[:,2,3] == time_dim_val[:,2,3]
        arr_fixed[:,1,1] = temp_dim_val[:,1,1]
        arr_variable[:,1,1] = temp_dim_val[:,1,1]
        @test arr_fixed[:,1,1] == arr_variable[:,1,1] == temp_dim_val[:,1,1]
        arr_fixed[:,:,2] = temp_dim_val[:,:,2]
        arr_variable[:,:,2] = temp_dim_val[:,:,2]
        @test arr_fixed[:,:,2] == arr_variable[:,:,2] == temp_dim_val[:,:,2]
        arr_fixed[:,:,:] = temp_dim_val
        arr_variable[:,:,:] = temp_dim_val
        @test arr_fixed[:,:,:] == arr_variable[:,:,:] == temp_dim_val[:,:,:]
    else
        @test arr_fixed[1,:,1] == arr_variable[1,:,1] == time_dim_val[1,:,1]
        @test arr_fixed[2,:,3] == arr_variable[2,:,3] == time_dim_val[2,:,3]
        arr_fixed[1,:,1] = temp_dim_val[1,:,1]
        arr_variable[1,:,1] = temp_dim_val[1,:,1]
        @test arr_fixed[1,:,1] == arr_variable[1,:,1] == temp_dim_val[1,:,1]
        arr_fixed[:,:,2] = temp_dim_val[:,:,2]
        arr_variable[:,:,2] = temp_dim_val[:,:,2]
        @test arr_fixed[:,:,2] == arr_variable[:,:,2] == temp_dim_val[:,:,2]
        arr_fixed[:,:,:] = temp_dim_val
        arr_variable[:,:,:] = temp_dim_val
        @test arr_fixed[:,:,:] == arr_variable[:,:,:] == temp_dim_val[:,:,:]
    end
    reset_time_val(arr_fixed, time_dim_val)
    reset_time_val(arr_variable, time_dim_val)

    @test_throws ErrorException arr_fixed[TimestepValue(2000)]
    @test_throws ErrorException arr_variable[TimestepValue(2000)]

    # Indexing with single TimestepIndex
    if ti == 1
        @test arr_fixed[idx1, 1, 1] == arr_variable[idx1, 1, 1] == time_dim_val[1,1,1]
        @test arr_fixed[TimestepIndex(3), 3, 3] == arr_variable[TimestepIndex(3), 3, 3] == time_dim_val[3,3,3]

        arr_fixed[idx1, 1, 1] = temp_dim_val[1,1,1]
        arr_variable[idx1, 1, 1] = temp_dim_val[1,1,1]
        @test arr_fixed[idx1, 1, 1] == arr_variable[idx1, 1, 1] == temp_dim_val[1,1,1]
        reset_time_val(arr_fixed, time_dim_val)
        reset_time_val(arr_variable, time_dim_val)

        @test_throws ErrorException arr_fixed[idx1]
        @test_throws ErrorException arr_variable[idx1]

        # Indexing with Array{TimestepIndex, N}
        @test arr_fixed[TimestepIndex.([1,3]), 1, 1] == time_dim_val[[1,3], 1, 1]
        @test arr_variable[TimestepIndex.([2,4]), 1, 1] == time_dim_val[[2,4], 1, 1]

        # Indexing with Array{TimestepIndex, N} created by Colon syntax
        @test arr_fixed[idx1:TimestepIndex(3), 1, 1] == time_dim_val[[1:3...], 1, 1]
        @test arr_fixed[idx1:2:TimestepIndex(3), 1, 1] == time_dim_val[[1:2:3...], 1, 1]

        # Indexing with single TimestepValue
        @test arr_fixed[TimestepValue(2000), 1, 1] == arr_variable[TimestepValue(2000), 1, 1] == time_dim_val[1,1,1]
        @test arr_fixed[TimestepValue(2010), 3, 3] == arr_variable[TimestepValue(2015), 3, 3] == time_dim_val[3,3,3]

        arr_fixed[TimestepValue(2000), 1, 1] = time_dim_val[1,1,1]
        arr_variable[TimestepValue(2000), 1, 1] = time_dim_val[1,1,1]
        @test arr_fixed[TimestepValue(2000), 1, 1] == arr_variable[TimestepValue(2000), 1, 1] == time_dim_val[1,1,1]
        reset_time_val(arr_fixed, time_dim_val)
        reset_time_val(arr_variable, time_dim_val)

        # Indexing with Array{TimestepValue, N}
        @test arr_fixed[TimestepValue.([2000, 2010]), 1, 1] == time_dim_val[[1,3],1,1]
        @test arr_variable[TimestepValue.([2000, 2005, 2025]), 1, 1] == time_dim_val[[1,2,4],1,1]

        arr_fixed[TimestepValue.([2000, 2010]), 1, 1] = temp_dim_val[[1,3],1,1]
        arr_variable[TimestepValue.([2000, 2005, 2025]), 1, 1] = temp_dim_val[[1,2,4],1,1]
        @test arr_fixed[TimestepValue.([2000, 2010]), 1, 1] == temp_dim_val[[1,3],1,1]
        @test arr_variable[TimestepValue.([2000, 2005, 2025]), 1, 1] == temp_dim_val[[1,2,4],1,1]

        reset_time_val(arr_fixed, time_dim_val)
        reset_time_val(arr_variable, time_dim_val)

        @test arr_fixed[idx1,2,3] == time_dim_val[1,2,3]
        @test arr_variable[idx1,2,3] == time_dim_val[1,2,3]
        arr_fixed[idx1,2,3] = temp_dim_val[1,2,3]
        arr_variable[idx1,2,3] = temp_dim_val[1,2,3]
        @test arr_fixed[idx1,2,3] == arr_variable[idx1,2,3] == temp_dim_val[1,2,3]

    else

        @test arr_fixed[1, idx1, 1] == arr_variable[1, idx1, 1] == time_dim_val[1,1,1]
        @test arr_fixed[3, TimestepIndex(3), 3] == arr_variable[3, TimestepIndex(3), 3] == time_dim_val[3,3,3]

        arr_fixed[1, idx1, 1] = temp_dim_val[1,1,1]
        arr_variable[1, idx1, 1] = temp_dim_val[1,1,1]
        @test arr_fixed[1, idx1, 1] == arr_variable[1, idx1, 1] == temp_dim_val[1,1,1]
        reset_time_val(arr_fixed, time_dim_val)
        reset_time_val(arr_variable, time_dim_val)

        # Indexing with Array{TimestepIndex, N}
        @test arr_fixed[1, TimestepIndex.([1,3]), 1] == time_dim_val[1, [1,3], 1]
        @test arr_variable[1, TimestepIndex.([2,4]), 1] == time_dim_val[1, [2,4], 1]

        # Indexing with Array{TimestepIndex, N} created by Colon syntax
        @test arr_fixed[1, idx1:TimestepIndex(3), 1] == time_dim_val[1, [1:3...], 1]
        @test arr_fixed[1, idx1:2:TimestepIndex(3), 1] == time_dim_val[1, [1:2:3...], 1]

        # Indexing with single TimestepValue
        @test arr_fixed[1, TimestepValue(2000), 1] == arr_variable[1, TimestepValue(2000), 1] == time_dim_val[1,1,1]
        @test arr_fixed[3, TimestepValue(2010), 3] == arr_variable[3, TimestepValue(2015), 3] == time_dim_val[3,3,3]

        arr_fixed[1, TimestepValue(2000), 1] = temp_dim_val[1,1,1]
        arr_variable[1, TimestepValue(2000), 1] = temp_dim_val[1,1,1]
        @test arr_fixed[1, TimestepValue(2000), 1] == arr_variable[1, TimestepValue(2000), 1] == temp_dim_val[1,1,1]
        reset_time_val(arr_fixed, time_dim_val)
        reset_time_val(arr_variable, time_dim_val)

        # Indexing with Array{TimestepValue, N}
        @test arr_fixed[1, TimestepValue.([2000, 2010]), 1] == time_dim_val[1, [1,3],1]
        @test arr_variable[1, TimestepValue.([2000, 2005, 2025]), 1] == time_dim_val[1,[1,2,4],1]
        arr_fixed[1, TimestepValue.([2000, 2010]), 1] = temp_dim_val[1, [1,3],1]
        arr_variable[1, TimestepValue.([2000, 2005, 2025]), 1] = temp_dim_val[1,[1,2,4],1]
        @test arr_fixed[1, TimestepValue.([2000, 2010]), 1] == temp_dim_val[1, [1,3],1]
        @test arr_variable[1, TimestepValue.([2000, 2005, 2025]), 1] == temp_dim_val[1,[1,2,4],1]

        reset_time_val(arr_fixed, time_dim_val)
        reset_time_val(arr_variable, time_dim_val)

        @test arr_fixed[1,idx2,3] == time_dim_val[1,2,3]
        @test arr_variable[1,idx2,3] == time_dim_val[1,2,3]
        arr_fixed[1,idx2,3] = temp_dim_val[1,2,3]
        arr_variable[1,idx2,3] = temp_dim_val[1,2,3]
        @test arr_fixed[1,idx2,3] == arr_variable[1,idx2,3] == temp_dim_val[1,2,3]
    end

    reset_time_val(arr_fixed, time_dim_val)
    reset_time_val(arr_variable, time_dim_val)
end

# other methods
time_dim_val = collect(reshape(1:64, 4, 4, 4))

x_years = Tuple(2000:5:2015) #fixed
y_years = Tuple([2000:5:2005; 2015:10:2025]) #variable

x_vec = TimestepVector{FixedTimestep{2000, 5}, Int}(time_dim_val[:,1,1])
x_mat = TimestepMatrix{FixedTimestep{2000, 5}, Int, 1}(time_dim_val[:,:,1])
y_vec = TimestepVector{VariableTimestep{y_years}, Int}(time_dim_val[:,2,2])
y_mat = TimestepMatrix{VariableTimestep{y_years}, Int, 1}(time_dim_val[:,:,2])

@test first_period(x_vec) == first_period(x_mat) == x_years[1]
@test first_period(y_vec) == first_period(y_mat) == y_years[1]
@test last_period(x_vec) == last_period(x_mat) == x_years[end]
@test last_period(y_vec) == last_period(y_mat) == y_years[end]

@test size(x_vec) == size(y_vec) == (4,)
@test size(x_mat) == size(y_mat) == (4,4)

@test ndims(x_vec) == ndims(y_vec) == 1
@test ndims(x_mat) == ndims(y_mat) == 2

@test eltype(x_vec) == eltype(y_vec) == eltype(y_vec) == eltype(y_mat) == eltype(time_dim_val)

#------------------------------------------------------------------------------
# 6. Test that getindex for TimestepArrays doesn't allow access to `missing`
#       values during `run` that haven't been computed yet.
#------------------------------------------------------------------------------

@defcomp foo begin
    par1 = Parameter(index=[time])
    var1 = Variable(index=[time])
    function run_timestep(p, v, d, t)
        if is_last(t)
            v.var1[t] = 0
        else
            v.var1[t] = p.par1[t+1]   # This is where the error will be thrown, if connected to an internal variable that has not yet been computed.
        end
    end
end

@defcomp bar begin
    par2 = Parameter(index=[time])
    var2 = Variable(index=[time])
    function run_timestep(p, v, d, t)
        if is_last(t)
            v.var2[t] = 0
        else
            v.var2[t] = p.par2[t+1]   # This is where the error will be thrown, if connected to an internal variable that has not yet been computed.
        end
    end
end

years = 2000:2010

m = Model()
set_dimension!(m, :time, years)
add_comp!(m, foo, :first)
add_comp!(m, bar, :second)
connect_param!(m, :second => :par2, :first => :var1)
set_param!(m, :first, :par1, 1:length(years))

@test_throws MissingException run(m)

# Check broadcast assignment to underlying array
x = Mimi.TimestepVector{Mimi.FixedTimestep{2005,10}, Float64}(zeros(10))
x[:] .= 10
@test all(x.data .== 10)

#------------------------------------------------------------------------------
# 7. Test TimestepArrays with time not as the first dimension
#------------------------------------------------------------------------------

@defcomp gdp begin
    growth = Parameter(index=[regions, foo, time, 2])   # test that time is not first but not last
    gdp = Variable(index=[regions, foo, time, 2])
    gdp0 = Parameter(index=[regions, foo, 2])

    pgrowth = Parameter(index=[regions, 3, time])       # test time as last
    pop = Variable(index=[regions, 3, time])

    mat = Parameter(index=[regions, time])              # test time as last for a matrix
    mat2 = Variable(index=[regions, time])

    function run_timestep(p, v, d, ts)
        if is_first(ts)
            v.gdp[:, :, ts, :] = (1 .+ p.growth[:, :, ts, :]) .* p.gdp0
            v.pop[:, :, ts] = zeros(2, 3)
        else
            v.gdp[:, :, ts, :] = (1 .+ p.growth[:, :, ts, :]) .* v.gdp[:, :, ts-1, :]
            v.pop[:, :, ts] = v.pop[:, :, ts-1] .+ p.pgrowth[:, :, ts]
        end
        v.mat2[:, ts] = p.mat[:, ts]
    end
end

time_index = 2000:2100
regions = ["OECD","non-OECD"]
nsteps=length(time_index)

m = Model()
set_dimension!(m, :time, time_index)
set_dimension!(m, :regions, regions)
set_dimension!(m, :foo, 3)
add_comp!(m, gdp)
set_param!(m, :gdp, :gdp0, [3; 7] .* ones(length(regions), 3, 2))
set_param!(m, :gdp, :growth, [0.02; 0.03] .* ones(length(regions), 3, nsteps, 2))
set_leftover_params!(m, Dict{String, Any}([
    "pgrowth" => ones(length(regions), 3, nsteps),
    "mat" => rand(length(regions), nsteps)
]))
run(m)

@test size(m[:gdp, :gdp]) == (length(regions), 3, length(time_index), 2)

@test all(!ismissing, m[:gdp, :gdp])
@test all(!ismissing, m[:gdp, :pop])
@test all(!ismissing, m[:gdp, :mat2])

end #module


text:
abstract type TagBase end

struct PayloadGroup
    range::NTuple{2,Int32}
    ratio::Float64
end

PayloadGroup(range::Vector{<:Integer}, ratio::T) where {T} =
    PayloadGroup((range[1], range[2]), get_ratio(ratio))

struct DataStreamConfig
    name::String
    payload_groups::Vector{PayloadGroup}
end

struct DataSetConfig
    name::String
    basename::String
    datagen_backend_tag::TagBase
    io_backend_tag::TagBase
    nsteps::Int32
    compute_seconds::Float64
    streams::Vector{DataStreamConfig}
end

abstract type DataObject end

mutable struct DataVector <: DataObject
    vec::Vector{Float64}
end

DataVector() = DataVector(Float64[])

mutable struct DataSet
    cfg::DataSetConfig
    curr_step::Int32
    timestamp::Float64
    streams::Vector{DataObject}
end

DataSet(cfg::DataSetConfig) = DataSet(cfg, 1, time(), Float64[])


text:
using Test
using Flux
using Serialization
using ConvOptDL: solve_qp_batch_with_back, solve_qp_batch

# https://optimization.mccormick.northwestern.edu/index.php/Quadratic_programming#Numerical_example
Qs = reshape([[3.0 1.0]; [1.0 1.0]], 2, 2, 1)
ps = reshape([1.0, 6.0], 2, 1)
Gs = reshape([[-2.0 -3.0]; [-1.0 0.0]; [0.0 -1.0]], 3, 2, 1)
hs = reshape([-4.0, 0.0, 0.0], 3, 1)
As = reshape([1.0 1.0], 1, 2, 1)
bs = reshape([1.5], 1, 1)

@testset "QP Adjoint Test 1" begin
    _, back = solve_qp_batch_with_back(Qs, ps, Gs, hs, As, bs)
    ΔX = ones(2, 1)
    Δλs = zeros(3, 1)
    Δνs = zeros(1, 1)
    Δ = (ΔX, Δλs, Δνs)
    ΔQs, Δps, ΔGs, Δhs, ΔAs, Δbs, _ = back(Δ)
    ref_grads = Serialization.deserialize("ref_grads.jls")
    @test ΔQs ≈ ref_grads["dQs"] atol=1e-16
    @test Δps ≈ ref_grads["dps"] atol=1e-16
    @test ΔGs ≈ ref_grads["dGs"] atol=1e-14
    @test Δhs ≈ ref_grads["dhs"] atol=1e-16
    @test ΔAs ≈ ref_grads["dAs"] atol=1e-5
    @test Δbs ≈ ref_grads["dbs"] atol=1e-16
end

function dummy(Qs, ps, Gs, hs, As, bs)
    X, _, _ = solve_qp_batch(Qs, ps, Gs, hs, As, bs)
    sum(X)
end

@testset "QP Adjoint Test 2" begin
    ΔQs, Δps, ΔGs, Δhs, ΔAs, Δbs = gradient(dummy , Qs, ps, Gs, hs, As, bs)
    ref_grads = Serialization.deserialize("ref_grads.jls")
    @test ΔQs ≈ ref_grads["dQs"] atol=1e-16
    @test Δps ≈ ref_grads["dps"] atol=1e-16
    @test ΔGs ≈ ref_grads["dGs"] atol=1e-14
    @test Δhs ≈ ref_grads["dhs"] atol=1e-16
    @test ΔAs ≈ ref_grads["dAs"] atol=1e-5
    @test Δbs ≈ ref_grads["dbs"] atol=1e-16
end


text:
function crossover(parent1::NEATInd, parent2::NEATInd, cfg::Dict)
    # Get dicts {innovation_number: connection}
    p1_genes=Dict()
    p2_genes=Dict()

    # Add genes or parent 1
    for c in parent1.connections
        p1_genes[c.innovation]= c
    end

    # Add genes or parent 2
    for c in parent2.connections
        p2_genes[c.innovation]= c
    end

    # Get IDs of all neurons used by the individual
    n_in = cfg["n_in"]
    n_out = cfg["n_out"]
    neuron_numbers = Array(1:n_in + n_out) # In / out neurons

    # Child connections
    child_connections::Array{Connection}(undef, 0)
    for i in 1:cfg["innovation_max"]
        new_con = nothing
        # Both parents have the gene: pick random
        if i in keys(p1_genes) && i in keys(p2_genes)
            if rand()>0/5
                new_con =p1_genes[i]
            else
                new_con =p2_genes[i]
            end

        # Only parent1 has the gene: pick if its fitness is higher
        elseif i in keys(p1_genes) and parent1.fitness >= parent2.fitness
            new_con =p1_genes[i]

        # Only parent2 has the gene: pick if its fitness is higher
        elseif i in keys(p2_genes) and parent2.fitness >= parent1.fitness
            new_con =p2_genes[i]
        end

        # Add the connection to the child, add the neuron IDs
        if new_con != nothing
            push!(child_connections, new_con)
            if !(new_con.in_node in neuron_numbers)
                push!(neuron_numbers, new_con.in_node)
            end
            # Add out_neuron if not already in
            if !(new_con.out_node in neuron_numbers)
                push!(neuron_numbers, new_con.out_node)
            end
        end
    end

    # Child neurons
    neurons = Array{Neuron}(undef, 0)
    for i in 1:length(neuron_numbers)
        push!(neurons, Neuron(0.0, 0.0, false))
    end
    # Child fitness
    fitness = -Inf .* ones(cfg["d_fitness"])
    NEATInd(n_in, n_out, 0, child_connections, neurons, fitness)

end


text:
for F in (:signbit, :sign)
    @eval begin
        $F(x::T) where {T<:MarkableFloat64} = $F(Float64(x))
        $F(x::T) where {T<:MarkableFloat32} = $F(Float32(x))
        $F(x::T) where {T<:MarkableFloat16} = $F(Float16(x))
    end
end

for F in (:abs, :float)
    @eval begin
        $F(x::T) where {T<:MarkableFloat64} = T($F(Float64(x)))
        $F(x::T) where {T<:MarkableFloat32} = T($F(Float32(x)))
        $F(x::T) where {T<:MarkableFloat16} = T($F(Float16(x)))
    end
end

for F in (:flipsign, :copysign)
    @eval begin
        $F(x::T, y::T) where {T<:MarkableFloat64} = T($F(Float64(x), Float64(y)))
        $F(x::T, y::T) where {T<:MarkableFloat32} = T($F(Float32(x), Float32(y)))
        $F(x::T, y::T) where {T<:MarkableFloat16} = T($F(Float16(x), Float16(y)))

        $F(x::T1, y::T2) where {T1<:MarkableFloat64, T2<:MarkableFloat64} =
            promote_type(T1,T2)($F(Float64(x), Float64(y)))
        $F(x::T1, y::T2) where {T1<:MarkableFloat32, T2<:MarkableFloat32} =
            promote_type(T1,T2)($F(Float32(x), Float32(y)))
        $F(x::T1, y::T2) where {T1<:MarkableFloat16, T2<:MarkableFloat16} =
            promote_type(T1,T2)($F(Float16(x), Float16(y)))

        $F(x::T, y::R) where {T<:MarkableFloat64, R<:AbstractFloat} = T($F(Float64(x), Float64(y)))
        $F(x::T, y::R) where {T<:MarkableFloat32, R<:AbstractFloat} = T($F(Float32(x), Float32(y)))
        $F(x::T, y::R) where {T<:MarkableFloat16, R<:AbstractFloat} = T($F(Float16(x), Float16(y)))

        $F(x::R, y::T) where {T<:MarkableFloat64, R<:AbstractFloat} = T($F(Float64(x), Float64(y)))
        $F(x::R, y::T) where {T<:MarkableFloat32, R<:AbstractFloat} = T($F(Float32(x), Float32(y)))
        $F(x::R, y::T) where {T<:MarkableFloat16, R<:AbstractFloat} = T($F(Float16(x), Float16(y)))
    end
end


text:
typealias AVecI AbstractVector{Int}
typealias AMatI AbstractMatrix{Int}
typealias VecI Vector{Int}
typealias MatI Matrix{Int}
typealias IntIterable AVec{Int}
typealias FloatIterable AVecF

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

immutable TransposeView{T} <: AbstractMatrix{T}
  mat::Matrix{T}
end
Base.getindex(tv::TransposeView, i::Integer, j::Integer) = tv.mat[j,i]
Base.setindex!{T}(tv::TransposeView{T}, val::T, i::Integer, j::Integer) = (tv.mat[j,i] = val)
Base.length(tv::TransposeView) = length(tv.mat)
Base.size(tv::TransposeView) = (ncols(tv.mat), nrows(tv.mat))

# import Base: *, ctranspose
# function *{T<:Real}(tv::TransposeView{T}, v::AVec{T})
#   res = similar(v, T, nrows(tv))
#   for i in 1:nrows(tv)
#     res[i] = dot(col(tv.mat, i), v)
#   end
#   res
# end

# ctranspose(tv::TransposeView) = tv.mat

####################################################

donothing(x...) = return
nop(x) = x
returntrue(x) = true

unzip(x) = [p[1] for p in x], [p[2] for p in x]

sizes(x) = map(size, x)

mapf(functions, obj) = [f(obj) for f in functions]

row(M::AMat, i::Integer) = rowvec_view(M, i)
col(M::AMat, i::Integer) = view(M, :, i)
column(M::AMat, i::Integer) = view(M, :, i)

rows(M::AMat, rng::AVecI) = view(M, rng, :)
cols(M::AMat, rng::AVecI) = view(M, :, rng)
columns(M::AMat, rng::AVecI) = view(M, :, rng)

row!{T}(M::AMat{T}, i::Integer, v::AVec{T}) = (M[i,:] = v; nothing)
col!{T}(M::AMat{T}, i::Integer, v::AVec{T}) = (M[:,i] = v; nothing)
column!{T}(M::AMat{T}, i::Integer, v::AVec{T}) = (M[:,i] = v; nothing)

function if0then1(x::Real)
  if isnan(x) || x == 0
    one(typeof(x))
  else
    x
  end
end

####################################################

getPctOfInt(pct::Float64, T::Int) = round(Int, max(0., min(1., pct)) * T)

function splitRange(pct::Float64, T::Int)
  lastin = getPctOfInt(pct,T)
  1:lastin, lastin+1:T
end

function splitMatrixRows(mat::Matrix, pct::Float64)
  rng1, rng2 = splitRange(pct, nrows(mat))
  rows(mat,rng1), rows(mat,rng2)
end

function apply(f::Function, A::AbstractArray)
  for x in A
    f(x)
  end
end

function foreach(A::AbstractArray, f::Function, fs::Function...)
  for x in A
    f(x)
  end
  for g in fs
    foreach(A, g)
  end
  A
end

function getLinspace(n, h)
  h = n > 1 ? h : 0
  linspace(-h, h, n)
end


text:
const JNI_VERSION_1_1 = convert(Cint, 0x00010001)
const JNI_VERSION_1_2 = convert(Cint, 0x00010002)
const JNI_VERSION_1_4 = convert(Cint, 0x00010004)
const JNI_VERSION_1_6 = convert(Cint, 0x00010006)
const JNI_VERSION_1_8 = convert(Cint, 0x00010008)

const JNI_TRUE = convert(Cchar, 1)
const JNI_FALSE = convert(Cchar, 0)

# Return Values
const JNI_OK           = convert(Cint, 0)               #/* success */
const JNI_ERR          = convert(Cint, -1)              #/* unknown error */
const JNI_EDETACHED    = convert(Cint, -2)              #/* thread detached from the VM */
const JNI_EVERSION     = convert(Cint, -3)              #/* JNI version error */
const JNI_ENOMEM       = convert(Cint, -4)              #/* not enough memory */
const JNI_EEXIST       = convert(Cint, -5)              #/* VM already created */
const JNI_EINVAL       = convert(Cint, -6)              #/* invalid arguments */

function javahome_winreg()
    try
        keypath = "SOFTWARE\\JavaSoft\\Java Runtime Environment"
        value = querykey(WinReg.HKEY_LOCAL_MACHINE, keypath, "CurrentVersion")
        keypath *= "\\"*value
        return querykey(WinReg.HKEY_LOCAL_MACHINE, keypath, "JavaHome")
    catch
        try
            keypath = "SOFTWARE\\JavaSoft\\Java Development Kit"
            value = querykey(WinReg.HKEY_LOCAL_MACHINE, keypath, "CurrentVersion")
            keypath *= "\\"*value
            return querykey(WinReg.HKEY_LOCAL_MACHINE, keypath, "JavaHome")
        catch
            error("Cannot find an installation of Java in the Windows Registry.  Please install Java from https://java.com/en/download/manual.jsp and subsequently reload JavaCall.")
        end
    end
end

@static Sys.isunix() ? (global const libname = "libjvm") : (global const libname = "jvm")

function findjvm()
    javahomes = Any[]
    libpaths = Any[]

    if haskey(ENV,"JAVA_HOME")
        push!(javahomes,ENV["JAVA_HOME"])
    else
        @static Sys.iswindows() ? ENV["JAVA_HOME"] = javahome_winreg() : nothing
        @static Sys.iswindows() ? push!(javahomes, ENV["JAVA_HOME"]) : nothing
    end

    if isfile("/usr/libexec/java_home")
        push!(javahomes,chomp(read(`/usr/libexec/java_home`, String)))
    end

    if isdir("/usr/lib/jvm/default-java/")
        push!(javahomes, "/usr/lib/jvm/default-java/")
    end

    push!(libpaths, pwd())
    for n in javahomes
        @static if Sys.iswindows()
            push!(libpaths, joinpath(n, "bin", "server"))
            push!(libpaths, joinpath(n, "jre", "bin", "server"))
            push!(libpaths, joinpath(n, "bin", "client"))
        end
        @static if Sys.islinux()
            if Sys.WORD_SIZE==64
                push!(libpaths, joinpath(n, "jre", "lib", "amd64", "server"))
                push!(libpaths, joinpath(n, "lib", "amd64", "server"))
			elseif Sys.WORD_SIZE==32
                push!(libpaths, joinpath(n, "jre", "lib", "i386", "server"))

                push!(libpaths, joinpath(n, "lib", "i386", "server"))
             end
         end
        push!(libpaths, joinpath(n, "jre", "lib", "server"))
        push!(libpaths, joinpath(n, "lib", "server"))
    end

    ext = @static Sys.iswindows() ? "dll" : (@static Sys.isapple() ? "dylib" : "so")
    ext = "."*ext

    try
        for n in libpaths
            libpath = joinpath(n,libname*ext);
            if isfile(libpath)
                if Sys.iswindows()
                    bindir = dirname(dirname(libpath))
                    m = filter(x -> occursin(r"msvcr(?:.*).dll",x), readdir(bindir))
                    Libdl.dlopen(joinpath(bindir,m[1]))
                end
                global libjvm = Libdl.dlopen(libpath)
                println("Loaded $libpath")
                return
            end
        end
    catch
    end

    errorMsg =
        [
         "Cannot find java library $libname$ext\n",
         "Search Path:"
         ];
    for path in libpaths
        push!(errorMsg,"\n   $path")
    end
    throw(JavaCallError(reduce(*,errorMsg)))
end

struct JavaVMOption
    optionString::Ptr{UInt8}
    extraInfo::Ptr{Nothing}
end

struct JavaVMInitArgs
    version::Cint
    nOptions::Cint
    options::Ptr{JavaVMOption}
    ignoreUnrecognized::Cchar
end

@static Sys.isunix() ? (const sep = ":") : nothing
@static Sys.iswindows() ? (const sep = ";") : nothing
cp = OrderedSet{String}()
opts = OrderedSet{String}()

function addClassPath(s::String)
    if isloaded()
        @warn("JVM already initialised. This call has no effect")
        return
    end
    if s==""; return; end
    if endswith(s, "/*") && isdir(s[1:end-2])
        for x in s[1:end-1] .* readdir(s[1:end-2])
            endswith(x, ".jar") && push!(cp, x)
        end
        return
    end
    push!(cp, s)
    return
end

addOpts(s::String) = isloaded() ? @warn("JVM already initialised. This call has no effect") : push!(opts, s)

function init()
    if isempty(cp)
        init(opts)
    else
        ccp = collect(cp)
        init(vcat(collect(opts), reduce((x,y)->string(x,sep,y),"-Djava.class.path=$(ccp[1])",ccp[2:end])))
    end
end

isloaded() = isdefined(JavaCall, :jnifunc) && isdefined(JavaCall, :penv) && penv != C_NULL

assertloaded() = isloaded() ? nothing : throw(JavaCallError("JVM not initialised. Please run init()"))
assertnotloaded() = isloaded() ? throw(JavaCallError("JVM already initialised")) : nothing

# Pointer to pointer to pointer to pointer alert! Hurrah for unsafe load
function init(opts)
    assertnotloaded()
    opt = [JavaVMOption(pointer(x), C_NULL) for x in opts]
    ppjvm = Array{Ptr{JavaVM}}(undef, 1)
    ppenv = Array{Ptr{JNIEnv}}(undef, 1)
    vm_args = JavaVMInitArgs(JNI_VERSION_1_6, convert(Cint, length(opts)),
                             convert(Ptr{JavaVMOption}, pointer(opt)), JNI_TRUE)
    res = ccall(create, Cint, (Ptr{Ptr{JavaVM}}, Ptr{Ptr{JNIEnv}}, Ptr{JavaVMInitArgs}), ppjvm, ppenv,
                Ref(vm_args))
    res < 0 && throw(JavaCallError("Unable to initialise Java VM: $(res)"))
    global penv = ppenv[1]
    global pjvm = ppjvm[1]
    jnienv = unsafe_load(penv)
    jvm = unsafe_load(pjvm)
    global jvmfunc = unsafe_load(jvm.JNIInvokeInterface_)
    global jnifunc = unsafe_load(jnienv.JNINativeInterface_) #The JNI Function table
    return
end

function destroy()
    if (!isdefined(JavaCall, :penv) || penv == C_NULL)
        throw(JavaCallError("Called destroy without initialising Java VM"))
    end
    res = ccall(jvmfunc.DestroyJavaVM, Cint, (Ptr{Nothing},), pjvm)
    res < 0 && throw(JavaCallError("Unable to destroy Java VM"))
    global penv=C_NULL; global pjvm=C_NULL;
end


text:
using Dierckx
using ASE
using JuLIP
using BenchmarkTools
using Main.Threads
using DataDeps
using Profile
using ProfileView
using Interpolations
using LinearAlgebra
using Test

register(DataDep(
    "JuLIP_testpots",
    "A few EAM potentials for testing",
    "https://www.dropbox.com/s/leub1c9ft1mm9fg/JuLIP_data.zip?dl=1",
    post_fetch_method = file -> run(`unzip $file`)
    ))
test_pots = joinpath(datadep"JuLIP_testpots", "JuLIP_data") * "/"
Ni = test_pots * "Ni.eam.fs"

eam0 = EAM(Ni)
eam1 = EAM(Ni; s = 1e-3)
eam = eam1
display(eam0.ϕ[1].f.spl)
display(eam1.ϕ[1].f.spl)

@show Main.Threads.nthreads()
@show JuLIP.nthreads()

##

spl_dx = eam0.ϕ[1].f.spl
knots = spl_dx.t[5:end-4]
@assert norm(diff(diff(knots))) < 1e-12
h = knots[2]-knots[1]
knots = (knots[1]-h):h:(knots[end]+3*h/2)
vals = spl_dx.(knots)
spl_i = interpolate(vals, BSpline(Cubic(Line(OnGrid()))))
spl_i = scale(spl_i, knots)

##

r = 2.5 + rand() - 0.5
spl_dx(r)
spl_i(r)

@info("Dierckx: eval, grad")
@btime ($spl_dx)($r)
@btime Dierckx.derivative($spl_dx, $r)

@info("Interpolations, eval, grad")
@btime ($spl_i)($r)
@btime Interpolations.gradient($spl_i, $r)[1]

##

function test(r)
    val_dx = spl_dx(r)
    val_i = spl_i(r)
    err = minimum([ abs(val_dx - val_i),
                    abs(val_dx - val_i) / (abs(val_i) + eps()) ]
                  )
    JuLIP.Testing.print_tf(@test err < 1e-12)
end

@info("test for correctness")
@info("  [1] random evaluation")
for _ = 1:100
    r = knots[1] + rand() * (knots[end]-knots[1])
    test(r)
end

@info("test for correctness")
@info("  [2] close to cutoff")
for _ = 1:100
    r = knots[end] - rand() * 1e-4
    test(r)
end

@info("test for correctness")
@info("  [3] even closer to cutoff")
for _ = 1:100
    r = knots[end] - rand() * 1e-12
    test(r)
end