text:
export GameTree, GameNode, GameAction
export train, selfplay, transition

mutable struct GameAction
    id::Int
    P::Float32
    nodeid::Int
    N::Int
    Q::Float32
end

function GameAction(id::Int, P::Float32)
    GameAction(id, P, 0, 0, Float32(0))
end

function update!(act::GameAction, value::Float32)
    act.Q = (act.N * act.Q + value) / (act.N + 1)
    act.N += 1
end

mutable struct GameNode
    state
    p
    v
    r::Float32
    actions::Vector{GameAction}
end

isfinal(node::GameNode) = isempty(node.actions)

function nextaction_search(node::GameNode)
    @assert !isempty(node.actions)
    length(node.actions) == 1 && return node.actions[1]
    sqrtN = sqrt(sum(a -> a.N, node.actions))
    cpuct = 1
    maxact, maxU = node.actions[1], typemin(Float32)

    for act in node.actions
        U = act.Q + cpuct * act.P * sqrtN / (act.N+1)
        if U > maxU
            maxU = U
            maxact = act
        end
    end
    maxact
end

function nextaction_play(node::GameNode, tau::Int)
    @assert !isempty(node.actions)
    length(node.actions) == 1 && return node.actions[1]

    if tau == 0
        maxact = node.actions[1]
        for act in node.actions
            if act.N > maxact.N
                maxact = act
            end
        end
        maxact
    elseif tau == 1
        N = sum(a -> a.N, node.actions)
        r = N * rand()
        total = 0
        for act in node.actions
            total += act.N
            total > r && return act
        end
        node.actions[end]
    else
        throw("Invalid tau: $tau")
    end
end

function getpi(node::GameNode)
    N = sum(a -> a.N, node.actions)
    @assert N > 0
    pi = zeros(node.p.data)
    for a in node.actions
        pi[a.id] = a.N / N
    end
    pi
end

mutable struct GameTree{T}
    state2id::Dict{T,Int}
    nodes::Vector{GameNode}
    model
end

function GameTree(initstate::T, model) where T
    tree = GameTree(Dict{T,Int}(), GameNode[], model)
    get!(tree, initstate)
    tree
end

function Base.get!(tree::GameTree{T}, state::T) where T
    id = get!(tree.state2id, state, length(tree.state2id)+1)
    if id > length(tree.nodes)
        p, v = tree.model(state)
        actids = legalactions(state)
        if isempty(actids) # final state
            r = reward(state)
            actions = Int[]
        else
            actdict = Dict(x => x for x in actids)
            for i = 1:length(p.data)
                haskey(actdict,i) || (p.data[i] = -1024)
            end

            P = softmax(p.data[actids])
            r = typemin(Float32)
            actions = Array{GameAction}(length(actids))
            for i = 1:length(actids)
                actions[i] = GameAction(actids[i], P[i])
            end
        end
        node = GameNode(state, p, v, r, actions)
        push!(tree.nodes, node)
        @assert length(tree.nodes) == length(tree.state2id)
    end
    id
end

function train(tree::GameTree, niters::Int, nplays::Int, nsearches::Int)
    opt = SGD()
    for iter = 1:niters
        println("iter: $iter")
        samples = GameNode[]
        prog = Progress(nplays)
        for i = 1:nplays
            nodes = selfplay(tree, nsearches)
            append!(samples, nodes)
            ProgressMeter.next!(prog)
        end
        opt.rate = 0.01 / length(samples) * sqrt(length(samples))

        p = concat(2, map(s -> s.p, samples))
        pi = cat(2, map(s -> getpi(s), samples)...)
        v = concat(2, map(s -> s.v, samples))
        r = map(s -> s.r, samples)
        r = reshape(r, 1, length(r))
        l = softmax_crossentropy(Var(pi), p)
        l += mse(v, Var(r))

        loss = sum(l.data) / length(samples)
        println("Loss:\t$loss")

        params = gradient!(l)
        foreach(opt, params)
    end
end

function selfplay(tree::GameTree, nsearches::Int)
    initstate = tree.nodes[1].state
    tree = GameTree(initstate, tree.model)
    node = tree.nodes[1]
    nodes = GameNode[]

    count = 1
    while !isfinal(node)
        push!(nodes, node)
        for _ = 1:nsearches
            search!(tree, node)
        end
        tau = count <= 3 ? 1 : 0
        # tau = 0
        act = nextaction_play(node, tau)
        node = tree.nodes[act.nodeid]
        count += 1
    end

    finalnode = node
    r = -finalnode.r
    for i = length(nodes):-1:1
        nodes[i].r = r
        r = -r
    end
    # push!(nodes, finalnode)
    nodes
end

function search!(tree::GameTree, rootnode::GameNode)
    node = rootnode
    actions = []
    while !isfinal(node)
        act = nextaction_search(node)
        push!(actions, act)
        if act.nodeid == 0
            state = transition(node.state, act.id)
            nodeid = get!(tree, state)
            act.nodeid = nodeid
            node = tree.nodes[nodeid]
            break
        else
            node = tree.nodes[act.nodeid]
        end
    end

    #value = isfinal(node) ? node.r : node.v.data[1]
    #for act in actions
    #    update!(act, value)
    #end
    v = isfinal(node) ? node.r : node.v.data[1]
    for i = length(actions):-1:1
        update!(actions[i], v)
        v = -v
    end
end


text:
######################################################################
# mcts.jl
#
# like the first one, but graph-based
######################################################################

using POMDPs
using MCTS

using FEBOL

# State type
struct TS
    x::Float64      # x-position of UAV
    y::Float64      # y-position of UAV
    h::Float64      # heading of UAV

    b::ParticleCollection{TargetTuple}

    r::Float64
    kids::Matrix{Int}
end
TS(x, b) = TS(x[1], x[2], x[3], b, 0.0, zeros(Int,4,2))
TS(x, b, r) = TS(x[1], x[2], x[3], b, r, zeros(Int,4,2))
TS(x, b, r, arr) = TS(x[1], x[2], x[3], b, r, arr)

const TA = NTuple{3,Float64}

struct Graph
    belief_nodes::Vector{TS}
end

mutable struct BMDP{S,TC,TB} <: MDP{TS, Int}
    sensor::S
    lambda::Float64

    action_list::Vector{TA}
    discount_factor::Float64    # does this do anything (not really)
    cache::TC
    pf::TB
    discrete_b::Matrix{Float64}
end
POMDPs.discount(::BMDP) = 1.0
function BMDP(lambda::Float64, tc, f)
    action_list = make_action_list(5, 8, [-20, 0, 20])
    discount_factor = 0.99

    return BMDP(FOV(), lambda, action_list, discount_factor, tc, f, zeros(50,50))
end

POMDPs.actions(bmdp::BMDP) = 1:length(bmdp.action_list)

function POMDPs.generate_sr(bmdp::BMDP, s::TS, a::Int, rng::MersenneTwister)
    na = length(bmdp.action_list)

    # TODO: handle heading and bounds correctly
    # update UAV state with action
    p = new_pose( (s.x,s.y,s.h), bmdp.action_list[a] )

    # sample a state from the belief
    #random_particle = sample_state(s.b)
    random_particle = rand(Base.GLOBAL_RNG, s.b)

    # generate a possible observation
    o = observe(random_particle, bmdp.sensor, p)

    # now check if this observation has an associated belief
    if s.kids[a,o+1] == 0
        # create new belief
        bp = update_b(bmdp.pf, s.b, p, o)

        # compute reward
        penalty = bmdp.lambda * FEBOL.nmac_penalty(bp, p, 15., bmdp.pf.n)
        r = penalty + cheap_entropy(bp, bmdp.discrete_b, 200.0)
        #r = penalty + cheap_entropy(bp, 200.0, 50)

        sp = TS(p, bp, -r, zeros(Int,na,2))

        # put sp into cache and update kids list
        push!(bmdp.cache, sp)
        s.kids[a, o+1] = length(bmdp.cache)
    else
        sp = bmdp.cache[s.kids[a, o+1]]
    end

    return sp, sp.r
end

# ok, let's create a FEBOL policy
struct MCTSPolicy{BM,P} <: FEBOL.Policy
    action_list::Vector{TA}
    bmdp::BM
    policy::P
    np::Int     # number of particles in reduced pf
end
function MCTSPolicy(f::PF;
                    lambda::Float64=0.0,
                    n_iterations::Int=50,
                    depth::Int=5
                   )
    action_list = make_action_list(5, 8, [-20, 0, 20])
    na = length(action_list)
    s = TS( (0.0,0.0,0.0), f._b, 0.0, zeros(Int, na, 2))
    bmdp = BMDP(lambda, [s], f)

    gr = MothRollout(action_list)
    solver = MCTSSolver(n_iterations=n_iterations,
                        depth=depth,
                        exploration_constant=5.0)
                        #estimate_value=RolloutEstimator(gr)
                       #)
    policy = solve(solver, bmdp)

    return MCTSPolicy(bmdp.action_list, bmdp, policy, f.n)
end
function FEBOL.action(m::SearchDomain, x::Vehicle, o, f::AbstractFilter, p::MCTSPolicy)
    npf = sample_states(f, p.np)

    # create cache and first node
    na = length(p.action_list)
    s = TS( (x.x,x.y,x.heading), npf, 0.0, zeros(Int, na, 2))
    p.bmdp.cache = [s]

    # determine the action
    a_idx = POMDPs.action(p.policy, s)

    return p.action_list[a_idx]
end

# just step to it
struct MothRollout <: POMDPs.Policy
    action_list::Vector{TA}
end
function POMDPs.action(p::MothRollout, s::TS)
    x = (s.x, s.y, s.h)
    dmin = Inf
    amin = 1
    tx = centroid(s.b)
    ai = 1
    for a in p.action_list
        xp = new_pose(x, a)
        d = FEBOL.get_distance2(x, tx)
        if dmin < d
            dmin = d
            amin = ai
        end
        ai += 1
    end

    return amin
end

#end #module

ext: jl


text:
"""
# Summary ScalePath
    mutable struct ScalePath <: Scaler
# Fields
    scale::VarOrNil
Applies scalar multiplication over a N-dimensional input.
This scalar called `scale` is a learnable parameter.
"""
mutable struct ScalePath <: Scaler
    scale::VarOrNil
    function ScalePath(scalar::AbstractFloat; ndims::Int, type::Type=Array{Float32})
        @assert ndims >= 1 "ndims >= 1 shall be met, but got ndims=$ndims"
        shape = ntuple(i->1, ndims)
        scale = Variable{type}(Zeros(type, shape) .+ eltype(type)(scalar), true, true, true);
        new(scale)
    end
    function ScalePath()
        new(nothing)
    end
end

function clone(this::ScalePath; type::Type=Array{Float32})
    cloned = ScalePath()
    cloned.scale = clone(this.scale, type=type)
    return cloned
end

function Base.show(io::IO, m::ScalePath)
    print(io, "ScalePath(scale=$(m.scale.value[1]); type=$(typeof(m.scale.value)))")
end

function paramsof(m::ScalePath)
    params = Vector{Variable}(undef,1)
    params[1] = m.scale
    return params
end

function xparamsof(m::ScalePath)
    xparams = Vector{XVariable}(undef,1)
    xparams[1] = ('w', m.scale)
    return xparams
end

function nparamsof(m::ScalePath)
    return 1
end

function bytesof(m::ScalePath, unit::String="MB")
    return blocksize(sizeof(m.scale), uppercase(unit))
end

function forward(m::ScalePath, x::Variable{T}) where T
    k = m.scale
    y = Variable{T}(ᵛ(x) .* ᵛ(k), x.backprop)

    if x.backprop
        function ScalePathBackward()
            if need2computeδ!(x) δ(x) .+=     δ(y) .* ᵛ(k)  end
            if need2computeδ!(k) δ(k) .+= sum(δ(y) .* ᵛ(x)) end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, ScalePathBackward)
    end
    return y
end

function predict(m::ScalePath, x::AbstractArray)
    k = ᵛ(m.scale)
    return ᵛ(x) .* k
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
## --- node update function (integrator) -----------
function upd(s,intr,slf,n)
    sss = similar(s)
    inp = intr*s
    for (i=1:n)
        sss[i] = (inp[i]>0)?1:(inp[i]<0)?0:(slf[i]==0)?0:s[i]
    end
    return sss
end

fname = "/tmp/jtb/sample_nws/randomnkp_one/11.29.15/1.dat"

edges = readdlm(fname,' ',Int64);
n = maximum(edges[:,1:2])
interc=zeros(Int64,n,n)
for i=1:length(edges[:,1])
	interc[edges[i,1],edges[i,2]]=edges[i,3]
end

selfact = diag(interc)                  # yeast specific
interc -= diagm(selfact)                # interaction matrix

fname = "/tmp/jtb/sample_nws/randomnkp_one/11.29.15/1.att"
a = readdlm(fname,' ',Any);
fps = int(a[:,1:n])
basin = a[:,n+1]
cyclen = int(a[:,n+2])
coop = a[:,n+3]
cycle = a[:,n+6]
s = fps[1,1:n]'
s = upd(s,interc,selfact,n)
println(s')

max_stars_repo_name: akabakcioglu/network-dynamics

max_issues_repo_name: akabakcioglu/network-dynamics

max_forks_repo_name: akabakcioglu/network-dynamics


text:
module InfoZIP

export open_zip, create_zip

have_infozip() = haskey(ENV, "HAVE_INFOZIP") || try
    occursin(r"^Copyright.* Info-ZIP", read(`zip -h`, String)) &&
    occursin(r"^UnZip.*by Info-ZIP", read(`unzip -h`, String))
catch ex
    if ex isa Base.IOError ex.code == Base.UV_ENOENT
        return false
    end
    rethrow(ex)
end

if have_infozip()
    include("info_zip.jl")
else
    @warn "InfoZIP falling back to ZipFile.jl backend!"
    include("zip_file.jl")
end

end # module


text:
"""
$(SIGNATURES)

Log jacobian abs determinant via automatic differentiation. For testing.
"""
AD_logjac(t::VectorTransform, x, vec_y) =
    logabsdet(ForwardDiff.jacobian(x -> vec_y(transform(t, x)), x))[1]

AD_logjac(t::ScalarTransform, x) =
    log(abs(ForwardDiff.derivative(x -> transform(t, x), x)))

"""
$(SIGNATURES)

Test transformation `t` with random values, `N` times.

`is_valid_y` checks the result of the transformation.

`vec_y` converts the result to a vector, for checking the log Jacobian with
automatic differentiation.

`test_inverse` determines whether the inverse is tested.
"""
function test_transformation(t::AbstractTransform, is_valid_y;
                             vec_y = identity, N = 1000, test_inverse = true)
    for _ in 1:N
        x = t isa ScalarTransform ? randn() : randn(dimension(t))
        if t isa ScalarTransform
            @test random_arg(t) isa Float64
        else
            y = random_arg(t)
            @test y isa Vector{Float64} && length(y) == dimension(t)
        end
        x isa ScalarTransform && @test dimension(x) == 1
        y = @inferred transform(t, x)
        @test is_valid_y(y)
        @test transform(t, x) == y
        y2, lj = @inferred transform_and_logjac(t, x)
        @test y2 == y
        if t isa ScalarTransform
            @test lj ≈ AD_logjac(t, x)
        else
            @test lj ≈ AD_logjac(t, x, vec_y)
        end
        if test_inverse
            x2 = inverse(t, y)
            @test x ≈ x2
            ι = inverse(t)
            @test x ≈ ι(y)
        end
    end
end

"""
$(SIGNATURES)

Elements (strictly) above the diagonal as a vector.
"""
function vec_above_diagonal(U::UpperTriangular{T}) where T
    n = size(U, 1)
    index = 1
    x = fill(T(NaN), unit_triangular_dimension(n))
    for col in 1:n
        for row in 1:(col-1)
            x[index] = U[row, col]
            index += 1
        end
    end
    x
end

"""
$(SIGNATURES)

Check if the argument `U` makes a valid correlation matrix `U'*U`.
"""
function is_valid_corr_cholesky(U::UpperTriangular)
    Ω = U'*U
    all(isapprox.(diag(Ω), 1)) && all(@. abs(Ω) ≤ (1+√eps()))
end


text:
using Test

@testset "ImageShow" begin
include("writemime.jl")
end


text:
module SpecialMatrices

import Base:  getindex, size, *

struct JordanBlock{T<:Number} <: AbstractMatrix{T}
  blocks::Vector{Int}
  diag::T
end

size(A::JordanBlock) = Tuple([1; 1] * sum(A.blocks))

function getindex(A::JordanBlock{T}, i::Int, j::Int) where T <: Number
  if i == j
    A.diag
  elseif i > j || j > i + 1
    zero(T)
  else
    if i in cumsum(A.blocks)
      zero(T)
    else
      one(T)
    end
  end
end

function *(J::JordanBlock, A::AbstractMatrix)
  n = size(J, 2)
  na, ma = size(A)
  n == na || error("dimension mismatch")
  M = A * J.diag
  for i = 1:n-1
    if J[i,i+1] != 0
      M[i,:] += A[i+1,:]
    end
  end
  M
end

end #module

max_stars_repo_name: KlausC/LinAlgBigFloat

max_issues_repo_name: KlausC/LinAlgBigFloat

max_forks_repo_name: KlausC/LinAlgBigFloat


text:
function calculate_transcription_rates(state::Array{Float64,1}, data_dictionary::Dict{String,Any}, vargs...)

    # initialize -
    transcription_kinetics_array = Float64[]

    # get stuff from the data_dictionary/biophysical_dictionary -
    species_symbol_type_array = data_dictionary["species_symbol_type_array"]
    gene_coding_length_array = data_dictionary["gene_coding_length_array"]
    time_constant_modifier_array = data_dictionary["time_constant_modifier_array"]
    species_concentration_array = data_dictionary["species_concentration_array"]
    biophysical_dictionary = data_dictionary["biophysical_dictionary"]
    characteristic_length = biophysical_dictionary["characteristic_transcript_length"]
    transcription_elongation_rate = biophysical_dictionary["transcription_elongation_rate"]
    characteristic_initiation_time = biophysical_dictionary["characteristic_initiation_time_transcription"]
    KX = biophysical_dictionary["transcription_saturation_constant"]
    RNAPII_concentration = biophysical_dictionary["RNAPII_concentration"]
    species_symbol_type_array = biophysical_dictionary["species_symbol_type_array"]

    # which indices are mRNA (the production of transcription)?
    idx_gene = findall(x->x==:GENE,species_symbol_type_array)

    # compute the kinetic limit -
    for (gene_index, gene_length) in enumerate(gene_coding_length_array)

        # how much gene do we have?
        # for host_type == :cell_free this is in concentration units -
        G = species_concentration_array[idx_gene[gene_index]]

        # what is the correction index?
        correction_index = idx_gene[gene_index]
        time_scale_parameter = time_constant_modifier_array[correction_index]

        # compute kE -
        kE = transcription_elongation_rate*(1/gene_length)

        # compute kI -
        kI = (1/characteristic_initiation_time)

        # compute the tau factor -
        tau_factor = (kE/kI)*time_scale_parameter

        # Compute the rate -
        sat_term = (G/(KX*tau_factor+(1+tau_factor)*G))
        value = kE*RNAPII_concentration*sat_term*(3600) # nmol/gDW-hr

        # push -
        push!(transcription_kinetics_array, value)
    end

    # return -
    return transcription_kinetics_array
end

max_stars_repo_name: varnerlab/CFMG

max_issues_repo_name: varnerlab/CFMG

max_forks_repo_name: varnerlab/CFMG


text:
type SmolyakInterp
    grid::SmolyakGrid
    f_on_grid::Vector{Float64}
    coefs::Vector{Float64}
end

SmolyakInterp(sg::SmolyakGrid, fgrid::AbstractVector) =
    SmolyakInterp(sg, fgrid, find_coefs(sg, fgrid))

# Use pre-computed LU fact of B to get coefficients
find_coefs(sg::SmolyakGrid, f_on_grid::AbstractVector) = sg.B_fact \ f_on_grid

# Convenience method to update coefficients.
update_coefs!(si::SmolyakInterp, f_on_grid) =
    (si.coefs = find_coefs(si.grid, f_on_grid))

# Interpolation
function interpolate!(out::AbstractVector, new_B::AbstractMatrix,
                      si::SmolyakInterp, pts::AbstractMatrix)
    d = size(pts, 2)

    d != si.grid.d && error("pts should have $d columns")

    cube_pts = dom2cube(pts, si.grid.lb, si.grid.ub)
    build_B!(new_B, d, si.grid.mu, cube_pts, si.grid.pinds)
    A_mul_B!(out, new_B, si.coefs)
end

function interpolate!(out::AbstractVector, si::SmolyakInterp, pts::AbstractMatrix)
    new_B = ones(size(pts, 1), size(si.grid.pinds, 1))
    interpolate!(out, new_B, si, pts)
end

interpolate(si::SmolyakInterp, pts::AbstractMatrix) =
    interpolate!(Array(Float64, size(pts, 1)), si, pts)

Base.call(si::SmolyakInterp, pts::AbstractMatrix) = interpolate(si, pts)

max_issues_repo_name: EconForge/Smolyak

max_forks_repo_name: EconForge/Smolyak


text:
function e008(N)
    return max_substring(long_number,N)
end

function solution008()
    N = 13
    answer = e008(N)
    println("The answer to problem 008 is $answer")
end

function max_substring(aString, substr_length)
    max::Int128 = 0
    for i = 1:length(aString)-substr_length+1
        substr = aString[i:i+substr_length-1]
        number::Int128 = 1
        for num in substr
            number *= parse(Int128,num)
        end
        if number > max
            max = number
        end
    end
    return max
end

const long_number = replace("73167176531330624919225119674426574742355349194934
96983520312774506326239578318016984801869478851843
85861560789112949495459501737958331952853208805511
12540698747158523863050715693290963295227443043557
66896648950445244523161731856403098711121722383113
62229893423380308135336276614282806444486645238749
30358907296290491560440772390713810515859307960866
70172427121883998797908792274921901699720888093776
65727333001053367881220235421809751254540594752243
52584907711670556013604839586446706324415722155397
53697817977846174064955149290862569321978468622482
83972241375657056057490261407972968652414535100474
82166370484403199890008895243450658541227588666881
16427171479924442928230863465674813919123162824586
17866458359124566529476545682848912883142607690042
24219022671055626321111109370544217506941658960408
07198403850962455444362981230987879927244284909188
84580156166097919133875499200524063689912560717606
05886116467109405077541002256983155200055935729725
71636269561882670428252483600823257530420752963450", "\n" => "")

max_stars_repo_name: mithodin/ttpe

max_issues_repo_name: mithodin/ttpe

max_forks_repo_name: mithodin/ttpe


text:
using Neurthino
using HDF5

@test Neurthino.cpphases(0) == 0
@test Neurthino.cpphases(1) == 0
@test Neurthino.cpphases(2) == 0
@test Neurthino.cpphases(3) == 1
@test Neurthino.cpphases(4) == 3

osc = Neurthino.OscillationParameters(3)
mixingangle!(osc, 1=>2, 0.59)
Neurthino.mixingangle!(osc, 1=>3, 0.15)
Neurthino.mixingangle!(osc, 2=>3, 0.84)
Neurthino.masssquareddiff!(osc, 2=>3, -2.523e-3)
Neurthino.masssquareddiff!(osc, 1=>2,-7.39e-5)
U = Neurthino.PMNSMatrix(osc)
H = Neurthino.Hamiltonian(osc)

@test isvalid(osc)

test_values = Neurthino.oscprob(U, H, 1, 1.6e4)[Energy=1, Baseline=1]
@test test_values[1,1] ≈ 0.142 atol=0.01
@test test_values[1,2] ≈ 0.420 atol=0.01
@test test_values[1,3] ≈ 0.436 atol=0.01
@test test_values[2,2] ≈ 0.257 atol=0.01
@test test_values[2,3] ≈ 0.323 atol=0.01
@test test_values[3,3] ≈ 0.240 atol=0.01

Neurthino.masssquareddiff!(osc, 3=>2, 2.523e-3)
Neurthino.masssquareddiff!(osc, 2=>1, 7.39e-5)

test_values = Neurthino.oscprob(osc, 1, 1.6e4)[Energy=1, Baseline=1]
@test test_values[1,1] ≈ 0.142 atol=0.01
@test test_values[1,2] ≈ 0.420 atol=0.01
@test test_values[1,3] ≈ 0.436 atol=0.01
@test test_values[2,2] ≈ 0.257 atol=0.01
@test test_values[2,3] ≈ 0.323 atol=0.01
@test test_values[3,3] ≈ 0.240 atol=0.01

Neurthino.cpphase!(osc, 1=>3, 234 * π / 180)

test_values = Neurthino.oscprob(osc, 1, 1.6e4; anti=true)[Energy=1, Baseline=1]
@test test_values[1,1] ≈ 0.142 atol=0.01
@test test_values[1,2] ≈ 0.479 atol=0.01
@test test_values[1,3] ≈ 0.378 atol=0.01
@test test_values[2,2] ≈ 0.289 atol=0.01
@test test_values[2,3] ≈ 0.462 atol=0.01
@test test_values[3,3] ≈ 0.159 atol=0.01

@test_logs (:warn, "Mass squared difference fields (partially) overdetermined!") Neurthino.masssquareddiff!(osc, 3=>1, 1)

for i in 3:100
    osc_params_dims = Neurthino.OscillationParameters(i)
    @test_throws ErrorException("Mass squared difference with equal index cannot be modified.") setΔm²!(osc_params_dims, 1=>1, 1)
    for j in 1:i-1
        @test_throws ErrorException("Mass squared differences not fully determined!") Hamiltonian(osc_params_dims)
        Neurthino.masssquareddiff!(osc_params_dims, j=>j+1, 1.0)
    @test_logs (:warn, "Mass squared difference fields (partially) overdetermined!") Neurthino.masssquareddiff!(osc, 3=>1, 1)
    end
end

osc = OscillationParameters(2)
mixingangle!(osc, 1=>2, 0.59)
masssquareddiff!(osc, 1=>2,-7.39e-5)
test_values = Neurthino.oscprob(osc, 1, 10000)[Energy=1, Baseline=1]
@test test_values[1,1] ≈ 0.4454 atol=0.01

osc = OscillationParameters(4)
mixingangle!(osc, 1=>2, 0.59)
mixingangle!(osc, 1=>3, 0.15)
mixingangle!(osc, 2=>3, 0.84)
mixingangle!(osc, 1=>4, 0.2)
mixingangle!(osc, 2=>4, 0.3)
mixingangle!(osc, 3=>4, 0.4)
masssquareddiff!(osc, 2=>3, -2.523e-3)
masssquareddiff!(osc, 1=>2,-7.39e-5)
masssquareddiff!(osc, 3=>4, -1)

h5open("data/refdata.h5", "r") do file
    # Nu-Fit v5.0 Values
    osc_nh = OscillationParameters(3);
    mixingangle!(osc_nh, 1=>2, 5.836e-1);
    mixingangle!(osc_nh, 1=>3, 1.496e-1);
    mixingangle!(osc_nh, 2=>3, 8.587e-1);
    cpphase!(osc_nh, 1=>3, 197 * π / 180);
    masssquareddiff!(osc_nh, 2=>3, -2.517e-3);
    masssquareddiff!(osc_nh, 1=>2, -7.42e-5);

    vac_baselines = collect(1:1.:1000)
    vac_energy = 0.01

    U_nh = PMNSMatrix(osc_nh)
    H_nh = Hamiltonian(osc_nh)

    data_vac_nh = Pνν(U_nh, H_nh, vac_energy, vac_baselines)[Energy=1];
    refdata = read(file, "vacuum/prob_nh")
    @test data_vac_nh ≈ refdata atol=0.01
end


text:
using MultipleViewGeometry, Test, LinearAlgebra

𝐊₁ = Matrix{Float64}(I, 3, 3)
𝐑₁ = Matrix{Float64}(I, 3, 3)
𝐭₁ = [1.0, 1.0, 1.0]
𝐊₂ = Matrix{Float64}(I, 3, 3)
𝐑₂ = Matrix{Float64}(I, 3, 3)
𝐭₂ = [2.0, 2.0, 2.0]

@test construct(FundamentalMatrix(),𝐊₁,𝐑₁,𝐭₁,𝐊₂,𝐑₂,𝐭₂) == vec2antisym([-1,-1,-1])


text:
using SparseArrays
using Test
using Languages
using TextAnalysis
using WordTokenizers
using Serialization

println("Running tests:")

include("coom.jl")
include("tokenizer.jl")
include("ngramizer.jl")
include("document.jl")
include("metadata.jl")
include("corpus.jl")
include("preprocessing.jl")
include("dtm.jl")
include("stemmer.jl")
include("tf_idf.jl")
include("lda.jl")
include("lsa.jl")
include("summarizer.jl")
include("bayes.jl")
include("taggingschemes.jl")
include("evaluation_metrics.jl")
include("LM.jl")


text:
@with_kw mutable struct Feedback
    t0::Float64=time()                      # start time
    t1::Float64=time()                      # time for duration estimation
    nans_detected::Bool=false               # did NaNs occur in the simulation?
    progress_txt::Union{IOStream,Nothing}   # txt is a Nothing in case of no output
    output::Bool                            # output to netCDF?
    i::Int=0                                # time step increment
    nt::Int                                 # number of time steps
    nout::Int                               # number of time steps with output
    run_id::Int                             # run identification number
    runpath::String                         # output path plus run????/
end

"""Returns a human readable string representing seconds in terms of days, hours, minutes or seconds."""
function readable_secs(secs::Real)
    days = Int(floor(secs/3600/24))
    hours = Int(floor((secs/3600) % 24))
    minutes = Int(floor((secs/60) % 60))
    seconds = Int(floor(secs%3600%60))
    secs1f = @sprintf "%.1fs" secs%3600%60
    secs2f = @sprintf "%.2fs" secs%3600%60

    if days > 0
        return "$(days)d, $(hours)h"
    elseif hours > 0
        return "$(hours)h, $(minutes)min"
    elseif minutes > 0
        return "$(minutes)min, $(seconds)s"
    elseif seconds > 10
        return secs1f
    else
        return secs2f
    end
end

"""Estimates the total time the model integration will take."""
function duration_estimate(feedback::Feedback,S::ModelSetup)

    @unpack t1,i,nt,output,progress_txt = feedback
    @unpack nadvstep = S.grid

    time_per_step = (time()-t1) / (i-nadvstep)
    time_total = Int(round(time_per_step*nt))
    time_to_go = Int(round(time_per_step*(nt-i)))

    s1 = "Model integration will take approximately "*readable_secs(time_total)*","
    s2 = "and is hopefully done on "*Dates.format(now() + Dates.Second(time_to_go),Dates.RFC1123Format)

    print("\r\u1b[K")
    println(s1)     # print inline
    println(s2)
    if output == 1  # print in txt
        write(progress_txt,"\n"*s1*"\n")
        write(progress_txt,s2*"\n")
        flush(progress_txt)
    end
end

"""Returns a boolean whether the prognostic variables contains a NaN."""
function nan_detection!(Prog::PrognosticVars,feedback::Feedback)

    #TODO include a check for Posits
    #TODO include check for sst

    @unpack u,v,η,sst = Prog

    n_nan = sum(isnan.(u)) + sum(isnan.(v)) + sum(isnan.(η)) + sum(isnan.(sst))
    if n_nan > 0
        feedback.nans_detected = true
    end
end

"""Initialises the progress txt file."""
function feedback_init(S::ModelSetup)
    @unpack output,init_starti = S.parameters
    @unpack nt,nout = S.grid

    if output

        @unpack Ndays,initial_cond,init_run_id,T,bc,α = S.parameters
        @unpack nx,ny,Δ = S.grid

        run_id,runpath = get_run_id_path(S)

        txt = open(joinpath(runpath,"progress.txt"),"w")
        s = "Starting ShallowWaters run $run_id on "*Dates.format(now(),Dates.RFC1123Format)
        println(s)
        write(txt,s*"\n")
        write(txt,"Model will integrate $(Ndays)days at a resolution of $(nx)x$(ny) with Δ=$(Δ/1e3)km\n")
        write(txt,"Initial conditions are ")
        if initial_cond == "rest"
            write(txt,"rest.\n")
        else
            if init_starti > 0
                write(txt,"time step $init_starti of run $init_run_id.\n")
            else
                write(txt,"last time step of run $init_run_id.\n")
            end
        end
        write(txt,"Boundary conditions are $bc with lbc=$α.\n")
        write(txt,"Number format is "*string(T)*".\n")
        write(txt,"\nAll data will be stored in $runpath\n")

        # Parameter.txt
        ptxt = open(joinpath(runpath,"parameter.txt"),"w")
        print(ptxt,S.parameters)
        close(ptxt)
    else
        println("Starting ShallowWaters on "*Dates.format(now(),Dates.RFC1123Format)*" without output.")
        txt = nothing
        run_id = -1
        runpath = ""
    end

    return Feedback(progress_txt=txt,output=output,nt=nt,nout=nout,run_id=run_id,runpath=runpath)
end

"""Feedback function that calls duration estimate, nan_detection and progress."""
function feedback!(Prog::PrognosticVars,feedback::Feedback,S::ModelSetup)

    @unpack nadvstep = S.grid
    @unpack i = feedback

    if feedback.output
        if i == nadvstep # measure time after tracer advection executed once
            feedback.t1 = time()
        elseif i == 2*nadvstep
            # after the tracer advection executed twice or at least 50 steps
            duration_estimate(feedback,S)
        end
    end

    @unpack i, nans_detected, nout, output = feedback

    if !nans_detected
        if i % nout == 0    # only check for nans when output is produced
            nan_detection!(Prog,feedback)
            if feedback.nans_detected
                println(" NaNs detected at time step $i")
                if output
                    @unpack progress_txt = feedback
                    write(progress_txt," NaNs detected at time step $i")
                    flush(progress_txt)
                end
            end
        end
    end

    if feedback.i > 100      # show percentage only after duration is estimated
        progress!(feedback)
    end
end

"""Finalises the progress txt file."""
function feedback_end!(feedback::Feedback)
    @unpack output,t0,progress_txt = feedback

    s = " Integration done in "*readable_secs(time()-t0)*"."
    println(s)
    if output
        write(progress_txt,"\n"*s[2:end]*"\n")  # close txt file with last output
        flush(progress_txt)
    end
end

"""Converts time step into percent for feedback."""
function progress!(feedback::Feedback)

    @unpack i,nt,progress_txt,output = feedback

    if ((i+1)/nt*100 % 1) < (i/nt*100 % 1)  # update every 1 percent steps.
        percent = Int(round((i+1)/nt*100))
        print("\r\u1b[K")
        print("$percent%")
        if output && (percent % 5 == 0) # write out only every 5 percent step.
            write(progress_txt,"\n$percent%")
            flush(progress_txt)
        end
    end
end


text:
function bulk_modulus(l, m)
    return l + 2m / 3
end

function lame_lambda(k, m)
    return k - 2m / 3
end

function displacement(alpha, x)
    u1 = alpha * x[2] * sin(pi * x[1])
    u2 = alpha * (x[1]^3 + cos(pi * x[2]))
    return [u1, u2]
end

function body_force(lambda, mu, alpha, x)
    b1 = alpha * (lambda + 2mu) * pi^2 * x[2] * sin(pi * x[1])
    b2 =
        -alpha * (6mu * x[1] + (lambda + mu) * pi * cos(pi * x[1])) +
        alpha * (lambda + 2mu) * pi^2 * cos(pi * x[2])
    return [b1, b2]
end

function stress_field(lambda, mu, alpha, x)
    s11 =
        (lambda + 2mu) * alpha * pi * x[2] * cos(pi * x[1]) -
        lambda * alpha * pi * sin(pi * x[2])
    s22 =
        -(lambda + 2mu) * alpha * pi * sin(pi * x[2]) +
        lambda * alpha * pi * x[2] * cos(pi * x[1])
    s12 = alpha * mu * (3x[1]^2 + sin(pi * x[1]))
    return [s11, s22, s12]
end

function onboundary(x, L, W)
    return x[2] ≈ 0.0 || x[1] ≈ L || x[2] ≈ W || x[1] ≈ 0.0
end

function assemble_linear_system(
    mesh,
    basis,
    cellquads,
    facequads,
    interfacequads,
    stiffness,
    penalty,
    eta,
    displacementscale,
)

    L, W = CutCellDG.widths(mesh)
    lambda1, mu1 = CutCellDG.lame_coefficients(stiffness, +1)

    sysmatrix = CutCellDG.SystemMatrix()
    sysrhs = CutCellDG.SystemRHS()

    CutCellDG.assemble_displacement_bilinear_forms!(
        sysmatrix,
        basis,
        cellquads,
        stiffness,
        mesh,
    )
    CutCellDG.assemble_interelement_condition!(
        sysmatrix,
        basis,
        facequads,
        stiffness,
        mesh,
        penalty,
        eta,
    )
    CutCellDG.assemble_coherent_interface_condition!(
        sysmatrix,
        basis,
        interfacequads,
        stiffness,
        mesh,
        penalty,
        eta,
    )
    CutCellDG.assemble_penalty_displacement_bc!(
        sysmatrix,
        sysrhs,
        x -> displacement(displacementscale, x),
        basis,
        facequads,
        stiffness,
        mesh,
        x -> onboundary(x, L, W),
        penalty,
    )
    CutCellDG.assemble_body_force!(
        sysrhs,
        x -> body_force(lambda1, mu1, displacementscale, x),
        basis,
        cellquads,
        mesh,
    )

    matrix = CutCellDG.sparse_displacement_operator(sysmatrix, mesh)
    rhs = CutCellDG.displacement_rhs_vector(sysrhs, mesh)

    return matrix, rhs
end

function nodal_displacement(
    distancefunc,
    stiffness,
    nelmts,
    basis,
    numqp,
    penaltyfactor,
    eta,
)
    L, W = 1.0, 1.0
    lambda1, mu1 = CutCellDG.lame_coefficients(stiffness, +1)
    lambda2, mu2 = CutCellDG.lame_coefficients(stiffness, -1)

    dx = 1.0 / nelmts
    penalty = penaltyfactor / dx * 0.5 * (lambda1 + mu1 + lambda2 + mu2)
    displacementscale = 0.1

    mesh = CutCellDG.DGMesh([0.0, 0.0], [L, W], [nelmts, nelmts], basis)
    levelset = InterpolatingPolynomial(1, basis)
    levelsetcoeffs = CutCellDG.levelset_coefficients(distancefunc, mesh)

    cutmesh = CutCellDG.CutMesh(mesh, levelset, levelsetcoeffs)
    cellquads =
        CutCellDG.CellQuadratures(cutmesh, levelset, levelsetcoeffs, numqp)
    interfacequads =
        CutCellDG.InterfaceQuadratures(cutmesh, levelset, levelsetcoeffs, numqp)
    facequads =
        CutCellDG.FaceQuadratures(cutmesh, levelset, levelsetcoeffs, numqp)

    mergedwithcell, hasmergedcells = CutCellDG.merge_tiny_cells_in_mesh!(
        cutmesh,
        cellquads,
        facequads,
        interfacequads,
        tinyratio=0.2
    )
    mergedmesh = CutCellDG.MergedMesh(cutmesh, mergedwithcell)

    matrix, rhs = assemble_linear_system(
        mergedmesh,
        basis,
        cellquads,
        facequads,
        interfacequads,
        stiffness,
        penalty,
        eta,
        displacementscale,
    )

    nodaldisplacement = matrix \ rhs

    return nodaldisplacement, mergedmesh, cellquads, facequads, interfacequads
end


text:
# This file assumes that example_setup.jl has been run already!!!

# nll and derivatives in the format Ipopt wants. This is also a simple template
# for exploiting redundant computations in ForwardDiff.
#
# This file expects example_setup.jl to have been run first, as well as the
# "vecc" object to have been created already.
#
# Note that there is a little bit of gaming about using SVectors, which really
# reduces the number of allocations when you use AD, thus helping with the
# multithreading. Or, at least, I think it should. Admittedly, it doesn't seem
# to make a huge difference. But that just means I have a problem to take care
# of somewhere else.

using ForwardDiff, StaticArrays, SIMDDualNumbers
const RESULTS = Dict{UInt64, DiffResults.MutableDiffResult}()

__nll(p) = nll(vecc, SVector{2,eltype(p)}(p))
function addkey!(p)
  newres = DiffResults.HessianResult(p)
  ForwardDiff.hessian!(newres, __nll, p)
  RESULTS[hash(p)] = newres
end

function _nll(p)
  haskey(RESULTS, hash(p)) && DiffResults.value(RESULTS[hash(p)])
  __nll(p)
end

function grad!(p, store)
  if haskey(RESULTS, hash(p))
    DiffResults.gradient(RESULTS[hash(p)])
    store .= collect(DiffResults.gradient(RESULTS[hash(p)]))
  else
    addkey!(p)
    store .= collect(DiffResults.gradient(RESULTS[hash(p)]))
  end
end

function hess(p)
  haskey(RESULTS, hash(p)) && DiffResults.hessian(RESULTS[hash(p)])
  addkey!(p)
  collect(DiffResults.hessian(RESULTS[hash(p)]))
end


text:
using AbstractTrees
import AbstractTrees: children, printnode
import Base: start, done, next

immutable Directory
    path::String
end

immutable File
    path::String
end

children(f::File) = ()

immutable DirectoryListing
    d::Directory
    names::Vector{String}
end

children(d::Directory) = DirectoryListing(d,readdir(d.path))

start(l::DirectoryListing) = start(l.names)
function next(l::DirectoryListing, i)
    (v,i) = next(l.names,i)
    path = joinpath(l.d.path,v)
    (isdir(path) ? Directory(path) : File(path), i)
end
done(l::DirectoryListing, i) = done(l.names, i)

# Pretty printing
printnode(io::IO, d::Directory) = Base.print_with_color(:blue, io, basename(d.path))
printnode(io::IO, f::File) = print(io, basename(f.path))


text:
@testset "Trie" begin

    t=Trie{Int}()

    t["amy"]=56
    t["ann"]=15
    t["emma"]=30
    t["rob"]=27
    t["roger"]=52

    @test haskey(t, "roger")
    @test get(t,"rob",nothing) == 27
    @test sort(keys(t)) == ["amy", "ann", "emma", "rob", "roger"]
    @test t["rob"] == 27
    @test sort(keys_with_prefix(t,"ro")) == ["rob", "roger"]

    # constructors
    ks = ["amy", "ann", "emma", "rob", "roger"]
    vs = [56, 15, 30, 27, 52]
    kvs = collect(zip(ks, vs))
    @test isa(Trie(ks, vs), Trie{Int})
    @test isa(Trie(kvs), Trie{Int})
    @test isa(Trie(Dict(kvs)), Trie{Int})
    @test isa(Trie(ks), Trie{Nothing})

    # path iterator
    t0 = t
    t1 = t0.children['r']
    t2 = t1.children['o']
    t3 = t2.children['b']
    @test collect(path(t, "b")) == [t0]
    @test collect(path(t, "rob")) == [t0, t1, t2, t3]
    @test collect(path(t, "robb")) == [t0, t1, t2, t3]
    @test collect(path(t, "ro")) == [t0, t1, t2]
    @test collect(path(t, "roa")) == [t0, t1, t2]

end # @testset Trie


text:
module HistogramThresholding

using Images
using LinearAlgebra

abstract type HistogramDistribution end
struct Unimodal <: HistogramDistribution end
struct Bimodal <: HistogramDistribution end

abstract type ThresholdAlgorithm end
struct Otsu <: ThresholdAlgorithm end
struct UniThreshold <: ThresholdAlgorithm end
struct Moments <: ThresholdAlgorithm end

include("otsu.jl")
include("unimodal.jl")
include("moments.jl")

export
	# main functions
  find_threshold,
	Otsu,
	UniThreshold,
	Moments

end # module


text:
"""
    get_path_to_example_config_files()::String

Returns the path to example config files provided by the package.
"""
function get_path_to_example_config_files()::String
    return joinpath(@__DIR__, "../examples/example_config_files")
end

const SSD_examples = Dict{Symbol,String}()
SSD_examples[:Coax] = (
    joinpath(@__DIR__, "../examples/example_detector_config_files/public_Coax_config.json")
)
SSD_examples[:InvertedCoax] = (
    joinpath(@__DIR__, "../examples/example_detector_config_files/public_ivc_config.json")
)
SSD_examples[:InvertedCoaxInCryostat] = (
    joinpath(@__DIR__, "../examples/example_detector_config_files/public_ivc_cryostat_config.json")
)
SSD_examples[:BEGe] = (
    joinpath(@__DIR__, "../examples/example_detector_config_files/public_SegBEGe_config.json")
)
SSD_examples[:CGD] = (
    joinpath(@__DIR__, "../examples/example_detector_config_files/public_CGD_config.json")
)
SSD_examples[:Spherical] = (
    joinpath(@__DIR__, "../examples/example_detector_config_files/public_spherical_detector_config.json")
)
SSD_examples[:SigGen] = (
    joinpath(@__DIR__, "../examples/example_detector_config_files/public_ppc_config_SigGen.config")
)
SSD_examples[:InfiniteParallelPlateCapacitor] = (
    joinpath(@__DIR__, "../examples/example_detector_config_files/infinite_parallel_plate_capacitor.json")
)
SSD_examples[:InfiniteCoaxialCapacitor] = (
    joinpath(@__DIR__, "../examples/example_detector_config_files/infinite_coaxial_capacitor.json")
)
SSD_examples[:InfiniteCoaxialCapacitorCartesianCoords] = (
    joinpath(@__DIR__, "../examples/example_detector_config_files/infinite_coaxial_capacitor_cartesian_coords.json")
)
SSD_examples[:Hexagon] = (
    joinpath(@__DIR__, "../examples/example_detector_config_files/minimum_hexagon_config.json")
)


text:
push!(LOAD_PATH,"../src/")
using ImageSmooth, TestImages, MosaicViews
using Documenter

DocMeta.setdocmeta!(ImageSmooth, :DocTestSetup, :(using ImageSmooth); recursive=true)

makedocs(;
    modules=[ImageSmooth],
    authors="Johnny Chen <johnnychen94@hotmail.com>, JKay Chen <jkaychen0327@hotmail.com>",
    repo="https://github.com/johnnychen94/ImageSmooth.jl/blob/{commit}{path}#{line}",
    sitename="ImageSmooth.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://johnnychen94.github.io/ImageSmooth.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Using ImageSmooth" => "usage.md",
        "Algorithms" => "algorithms.md",
        "Package References" => "reference.md",
    ],
)

deploydocs(;
    repo="github.com/johnnychen94/ImageSmooth.jl",
)


text:
################################################################################

# set hyperparameters
@with_kw mutable struct NNParams
  η::Float64                   = 1e-3             # learning rate
  epochs::Int                  = 25               # number of epochs
  batchsize::Int               = 1000             # batch size for training
  throttle::Int                = 5                # throttle timeout
  device::Function             = gpu              # set as gpu, if gpu available
  σ::Function                  = leakyrelu        # learning function
  λ::Int                       = 64               # hidden layer on autoencoder
end;

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

# hidden Markov model parameters
hmmParams = HMMParams(
  penalty = 200,
  distance = euclDist,
  minimumFrequency = 20,
  verbosity = false,
)

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

max_stars_repo_name: DanielRivasMD/MindReader

max_issues_repo_name: DanielRivasMD/MindReader

max_forks_repo_name: DanielRivasMD/MindReader


text:
using Ramnet.Encoders

using ..Approximators: ContinousActionPolicy
using ..Buffers: MultiStepDynamicBuffer, add!

mutable struct ContinousFunctionalPG{OS,AS,T<:Real,O<:StaticArray{Tuple{OS},T,1},A<:StaticArray{Tuple{AS},T,1},E<:AbstractEncoder{T}} <: AbstractAgent{O,A}
    γ::Float64
    policy::ContinousActionPolicy{OS,AS,T,O,A,E}
    buffer::MultiStepDynamicBuffer{O,A}
    observation::Union{Nothing,O}
    done::Bool
    action::Union{Nothing,A}
    min_obs::Union{Nothing,O}
    max_obs::Union{Nothing,O}
    # rng::MersenneTwister
end

function ContinousFunctionalPG(::Type{O}, ::Type{A}, n, start_learning_rate, end_learning_rate, learning_rate_decay, epochs, discount, encoder::E, start_cov::Float64, end_cov::Float64, cov_decay::Int; seed::Union{Nothing,UInt} = nothing) where {OS,AS,T<:Real,O<:StaticArray{Tuple{OS},T,1},A<:StaticArray{Tuple{AS},T,1},E<:AbstractEncoder{T}}
    # !isnothing(seed) && seed < 0 && throw(DomainError(seed, "Seed must be non-negative"))
    # rng = isnothing(seed) ? MersenneTwister() : MersenneTwister(seed)

    n < 1 && throw(DomainError(n, "Tuple size must be at least 1"))

    policy = ContinousActionPolicy(O, A, n, encoder; start_cov, end_cov, cov_decay, start_learning_rate, end_learning_rate, learning_rate_decay, epochs, partitioner = :uniform_random, seed)
    buffer = MultiStepDynamicBuffer{O,A}()

    ContinousFunctionalPG{OS,AS,T,O,A,E}(discount, policy, buffer, nothing, false, nothing, nothing, nothing)
end

function ContinousFunctionalPG(env, encoder::E; tuple_size, start_learning_rate, end_learning_rate, learning_rate_decay, epochs, discount, start_cov, end_cov, cov_decay, seed = nothing) where {T<:Real,E<:AbstractEncoder{T},C<:AbstractMatrix{T}}
    # ContinousFunctionalPG{action_length(env),T,observation_type(env),action_type(env),E,C}(n, observation_length(env), η, γ, encoder, cov; seed)
    ContinousFunctionalPG(observation_type(env), action_type(env), tuple_size, start_learning_rate, end_learning_rate, learning_rate_decay, epochs, discount, encoder, start_cov, end_cov, cov_decay; seed)
end

function observe!(agent::ContinousFunctionalPG{OS,AS,T,O,A,E}, observation::O) where {OS,AS,T,O,A,E}
    agent.observation = observation
    agent.done = false

    # Finding the edges of the domain
    if isnothing(agent.min_obs)
        agent.min_obs = observation
    else
        agent.min_obs = min.(agent.min_obs, observation)
    end

    if isnothing(agent.max_obs)
        agent.max_obs = observation
    else
        agent.max_obs = max.(agent.max_obs, observation)
    end

    nothing
end

# TODO: All observe! methods are the same for all agents. Generalize.
# TODO: This method does not need to take in the action
function observe!(agent::ContinousFunctionalPG{OS,AS,T,O,A,E}, action::A, reward::Float64, observation::O, done::Bool) where {OS,AS,T,O,A,E}
    # REMEMBER: you are forcing reward clipping here!!
    add!(agent.buffer, agent.observation, agent.action, clamp(reward, -1.0, 1.0))

    agent.observation = observation
    agent.done = done

    # Finding the edges of the domain
    agent.min_obs = min.(agent.min_obs)
    agent.max_obs = max.(agent.max_obs)

    nothing
end

function select_action(agent::ContinousFunctionalPG{OS,AS,T,O,A,E}, observation::O; deterministic::Bool = false) where {OS,AS,T,O,A,E}
    return select_action(agent.policy, observation; deterministic)
end

function select_action!(agent::ContinousFunctionalPG{OS,AS,T,O,A,E}, observation::O) where {OS,AS,T,O,A,E}
    agent.action = select_action(agent, observation)

    return agent.action
end

function update!(agent::ContinousFunctionalPG)
    if agent.done
        G = 0.0
        for transition in Iterators.reverse(agent.buffer)
            G = transition.reward + agent.γ * G

            update!(agent.policy, transition.observation, transition.action, G)
        end

        reset!(agent.buffer)
    end

    nothing
end

function Agents.reset!(agent::ContinousFunctionalPG; seed::Union{Nothing,UInt} = nothing)
    reset!(agent.policy; seed)
    reset!(agent.buffer)

    agent.observation = nothing
    agent.action = nothing
    agent.done = false

    nothing
end


text:
function create_historical_scenarios(container:: FXForwardMarketDataContainer; valuation_date:: Union{Date, Nothing}=nothing, time_horizon:: Integer = 1, lookback:: Integer = 250)

    scenario_dates_raw = collect(keys(first(container.fx_forward_curves).second.curves))
    if valuation_date === nothing
        idx = length(scenario_dates_raw)
    else
        idx = findlast(scenario_dates_raw .<= valuation_date)
    end
    scenario_dates = @view scenario_dates_raw[1:idx]

    fx_scenarios = Dict(key => _calculate_returns(val, time_horizon, lookback, scenario_dates, idx)
        for (key, val) in pairs(container.fx_forward_curves))

    disc_scenarios = Dict(key => _calculate_returns(val, time_horizon, lookback, scenario_dates, idx)
        for (key, val) in pairs(container.discount_curves))

    return FXForwardMarketDataContainer(fx_scenarios, disc_scenarios)
end

function _calculate_returns(ts:: T, time_horizon:: Integer, lookback:: Integer, ts_dates, idx:: Int) where {T <: AbstractTimeSeries}

    ts_values_raw = collect(values(ts.curves))
    ts_values = @view ts_values_raw[1:idx]

    base_scenario = last(ts_values)
    values_in_lookback = @view ts_values[end-lookback-time_horizon:end-1] # excluding base scenario

    shifts_raw = values_in_lookback ./ lag(values_in_lookback, time_horizon)
    shifts = [x for x in shifts_raw if !ismissing(x)] # the lag function adds `missing` values at the beginning of the time series, they are removed here.

    scenarios = shifts .* base_scenario
    push!(scenarios, base_scenario) # last scenario is base scenario

    # the scenarios are from the end of the time period
    scenario_dates =  ts_dates[end-length(scenarios)+1:end]

    return T(SortedDict(scenario_dates .=> scenarios))

end

max_issues_repo_name: lungben/MarketRiskPrototype

max_forks_repo_name: lungben/MarketRiskPrototype


text:
@recipe function f(iterate::PararealIterate; label="")
    P = length(iterate)
    for n = 1:P
        @series begin
            label := n != P ? "" : label
            iterate[n]
        end
    end
    # solves a bug in Plots with OrdinaryDiffEq
    # primary := false
    # ()
end

@recipe function f(solution::PararealSolution)
    @↓ saveiterates = solution
    @series begin
        if saveiterates
            @↓ alliterates = solution
            return alliterates[end]
        else
            @↓ lastiterate = solution
            return lastiterate
        end
    end
end

@userplot CONVERGENCE
@recipe function f(h::CONVERGENCE; label="")
    if length(h.args) == 1 && h.args[1] isa PararealSolution
        solution = h.args[1]
    else
        error("convergence must be given PararealSolution. Got $(typeof(h.args))")
    end
    framestyle --> :box
    gridalpha --> 0.2
    linewidth --> 1.5
    markershape --> :circle
    markersize --> 3
    minorgrid --> 0.1
    minorgridstyle --> :dash
    seriestype --> :path
    tick_direction --> :out
    @series begin
        label := label
        solution.errors
    end
end


text:
#!/usr/bin/env julia

using SnaggedNecklace

using ArgParse

s = ArgParseSettings()
s.autofix_names = true
@add_arg_table! s begin
    "--mass"
        metavar = "M"
        help = "mass (g/mol)"
        arg_type = Float64
        required = true
    "--epsilon"
        metavar = "E"
        help = "LJ epsilon (kJ/mol)"
        arg_type = Float64
        required = true
    "--sigma"
        metavar = "S"
        help = "LJ sigma (nm)"
        arg_type = Float64
        required = true
    "--l-max"
        metavar = "L"
        help = "basis truncation"
        arg_type = Int
        required = true
    "--dr"
        metavar = "D"
        help = "grid size"
        arg_type = Float64
        required = true
    "--grid-max"
        metavar = "G"
        help = "grid truncation"
        arg_type = Float64
        required = true
    "--temp"
        metavar = "T"
        help = "temperature"
        arg_type = Float64
    "--temp-P"
        metavar = "TP"
        help = "temperature multiplied by P"
        arg_type = Float64
    "-P"
        metavar = "P"
        help = "number of beads"
        arg_type = Int
        required = true
    "--rho-file"
        metavar = "PATH"
        help = "path to density output file"
        arg_type = String
end
c = parse_args(ARGS, s, as_symbols=true)

const mass = c[:mass]
const epsilon = c[:epsilon]
const sigma = c[:sigma]
const l_max = c[:l_max]
const dr = c[:dr]
const grid_max = c[:grid_max]

if isnothing(c[:temp]) == isnothing(c[:temp_P])
    error("One of --temp or --temp-P must be specified.")
elseif !isnothing(c[:temp])
    const temperature = c[:temp]
else
    const temperature = c[:temp_P] / c[:P]
end

const P = c[:P]

const rho_file = c[:rho_file]

pot(r) = 4.0 * epsilon * ((sigma / r)^12 - (sigma / r)^6)

sys = NmmSystem(pot, mass, temperature)
np = NmmParameters(l_max, dr, grid_max)
grid, rho_diag = nmm(sys, np, P)

if !isnothing(rho_file)
    open(rho_file, "w") do f
        for i in eachindex(grid, rho_diag)
            println(f, "$(grid[i]) $(rho_diag[i])")
        end
    end
end


text:
using MultivariatePolynomials,
      DynamicPolynomials,
      JuMP,
      PolyJuMP,
      SumOfSquares,
      MathOptInterfaceMosek,
      SemialgebraicSets
using Reachability
using MathematicalSystems: PolynomialContinuousSystem, InitialValueProblem
const ∂ = differentiate

export XFZ18

"""
    XFZ18 <: ContinuousPost

Implementation of the reachability algorithm for the class of polynomial ODEs
with uncertain initial states (see [1], abbreviated `XFZ18`).
This method consists of reducing the Hamilton-Jacobi-Bellman equation to a
hierarchy of semidefinite programs that are solved using an SDP solver.

### Fields

- `options` -- an `Options` structure that holds the algorithm-specific options

### Algorithm

We refer to [1] for technical details.

- [1] *Xue, B., Fränzle, M., & Zhan, N. (2018, April). Under-Approximating Reach
Sets for Polynomial Continuous Systems. In [Proceedings of the 21st International
Conference on Hybrid Systems: Computation and Control (part of CPS Week)
(pp. 51-60). ACM.](https://dl.acm.org/citation.cfm?id=3178133)*
"""
struct XFZ18 <: ContinuousPost
    options::Options

    function XFZ18(𝑂::Options)

        # merge!(𝑂, :relaxation_degree)

        return new(𝑂)
    end
end

# build the SOS problem
function build_sos(𝒮, opt)
    T = opt[:T]

    # scale dynamics
    T = opt[:T]
    f = T * 𝒮.s.p

    # define polynomial symbolic variables
    vars = variables(𝒮.s) # TODO: variables(𝒮)
    @polyvar t

    k = opt[:relaxation_degree]

    # monomial vector up to order k
    # 0 <= sum_i alpha_i <= k, if alpha_i is the exponent of x_i
    X = monomials(vars, 0:k)
    XT = monomials([vars; t], 0:k)

    # create a SOS JuMP model to solve with Mosek
    solver = opt[:solver]
    model = SOSModel(with_optimizer(solver))

    # add unknown Φ to the model
    @variable(model, Φ, Poly(XT))

    # jacobian
    ∂t = α -> ∂(α, t)
    ∂xf = α -> ∂(α, x₁) * f[1] + ∂(α, x₂) * f[2]
    LΦ = ∂t(Φ) + ∂xf(Φ)

    # Φ(x, t) at time 0
    Φ₀ = subs(Φ, t => 0.)

    # scalar variable
    @variable(model, ϵ)

    dom1 = @set t*(T-t) >= 0 && g >= 0
    dom2 = @set g >= 0
    @constraint(model, ϵ >= 0.)
    @constraint(model, LΦ ∈ SOSCone(), domain = dom1)
    @constraint(model, ϵ - LΦ ∈ SOSCone(), domain = dom1)
    @constraint(model, Φ₀ - V₀ ∈ SOSCone(), domain = dom2)
    @constraint(model, ϵ + V₀ - Φ₀ ∈ SOSCone(), domain = dom2)

    @objective(model, Min, ϵ)
    return model
end

# solve model, check feasibility and return polynomials
function solve_sos(model; verbose=true)
    optimize!(model)

    if verbose
        println("Relaxation order : k = $k")
        println("JuMP.termination_status(model) = ", JuMP.termination_status(model))
        println("JuMP.primal_status(model) = ", JuMP.primal_status(model))
        println("JuMP.dual_status(model) = ", JuMP.dual_status(model))
        println("JuMP.objective_bound(model) = ", JuMP.objective_bound(model))
        println("JuMP.objective_value(model) = ", JuMP.objective_value(model))
    end

    # TODO: error if it fails?
end

function extract_approximations(model, 𝑂)

    # time horizon TODO : check consistency w/rescaling
    T = 𝑂[:T]

    # Recovering the solution:
    ϵopt = JuMP.objective_value(model)

    # Punder <= 0  TODO: @set Punder <= 0 ?
    Punder = subs(JuMP.value(model[:Φ]), t => T)

    # Pover <= 0   TODO: @set Pover <= 0 ?
    Pover = subs(JuMP.value(model[:Φ]), t => T) - ϵopt * (T+1)

    return (ϵopt, Punder, Pover)
end

function post(𝒫::XFZ18, 𝒮::AbstractSystem, 𝑂::Options)

    # dynamics
    @assert 𝒮.s isa PolynomialContinuousSystem

    # construct sum-of-squares problem
    model = build_sos(𝒫, 𝒮, 𝑂)

    # solve the sum-of-squares optimization
    solve_sos(model, 𝑂)

    # extract under and over approximations
    (ϵopt, Punder, Pover) = extract_approximations(model)

    # returns the polynomial under and overapproximations of the reach set
    # for any t ∈ [0, T]
    return (Punder, Pover)
end


text:
# adaptors and type aliases for working with array wrappers

using LinearAlgebra

permutation(::PermutedDimsArray{T,N,perm}) where {T,N,perm} = perm

export WrappedArray

adapt_structure(to, A::SubArray) =
      SubArray(adapt(to, Base.parent(A)), adapt(to, parentindices(A)))
adapt_structure(to, A::PermutedDimsArray) =
      PermutedDimsArray(adapt(to, Base.parent(A)), permutation(A))
adapt_structure(to, A::Base.ReshapedArray) =
      Base.reshape(adapt(to, Base.parent(A)), size(A))
@static if isdefined(Base, :NonReshapedReinterpretArray)
    adapt_structure(to, A::Base.NonReshapedReinterpretArray) =
          Base.reinterpret(Base.eltype(A), adapt(to, Base.parent(A)))
    adapt_structure(to, A::Base.ReshapedReinterpretArray) =
          Base.reinterpret(reshape, Base.eltype(A), adapt(to, Base.parent(A)))
else
    adapt_structure(to, A::Base.ReinterpretArray) =
          Base.reinterpret(Base.eltype(A), adapt(to, Base.parent(A)))
end
@eval function adapt_structure(to, A::Base.LogicalIndex{T}) where T
      # prevent re-calculating the count of booleans during LogicalIndex construction
      mask = adapt(to, A.mask)
      $(Expr(:new, :(Base.LogicalIndex{T, typeof(mask)}), :mask, :(A.sum)))
end

adapt_structure(to, A::LinearAlgebra.Adjoint) =
      LinearAlgebra.adjoint(adapt(to, Base.parent(A)))
adapt_structure(to, A::LinearAlgebra.Transpose) =
      LinearAlgebra.transpose(adapt(to, Base.parent(A)))
adapt_structure(to, A::LinearAlgebra.LowerTriangular) =
      LinearAlgebra.LowerTriangular(adapt(to, Base.parent(A)))
adapt_structure(to, A::LinearAlgebra.UnitLowerTriangular) =
      LinearAlgebra.UnitLowerTriangular(adapt(to, Base.parent(A)))
adapt_structure(to, A::LinearAlgebra.UpperTriangular) =
      LinearAlgebra.UpperTriangular(adapt(to, Base.parent(A)))
adapt_structure(to, A::LinearAlgebra.UnitUpperTriangular) =
      LinearAlgebra.UnitUpperTriangular(adapt(to, Base.parent(A)))
adapt_structure(to, A::LinearAlgebra.Diagonal) =
      LinearAlgebra.Diagonal(adapt(to, Base.parent(A)))
adapt_structure(to, A::LinearAlgebra.Tridiagonal) =
      LinearAlgebra.Tridiagonal(adapt(to, A.dl), adapt(to, A.d), adapt(to, A.du))
adapt_structure(to, A::LinearAlgebra.Symmetric) =
      LinearAlgebra.Symmetric(adapt(to, Base.parent(A)))

# we generally don't support multiple layers of wrappers, but some occur often
# and are supported by Base aliases like StridedArray.

WrappedReinterpretArray{T,N,Src} =
      Base.ReinterpretArray{T,N,<:Any,<:Union{Src,SubArray{<:Any,<:Any,Src}}}

WrappedReshapedArray{T,N,Src} =
      Base.ReshapedArray{T,N,<:Union{Src,
                                     SubArray{<:Any,<:Any,Src},
                                     WrappedReinterpretArray{<:Any,<:Any,Src}}}

WrappedSubArray{T,N,Src} =
      SubArray{T,N,<:Union{Src,
                           WrappedReshapedArray{<:Any,<:Any,Src},
                           WrappedReinterpretArray{<:Any,<:Any,Src}}}

"""
    WrappedArray{T,N,Src,Dst}

Union-type that encodes all array wrappers known by Adapt.jl. Typevars `T` and `N` encode
the type and dimensionality of the resulting container.

Two additional typevars are used to encode the parent array type: `Src` when the wrapper
uses the parent array as a source, but changes its properties (e.g.
`SubArray{T,1,Array{T,2}` changes `N`), and `Dst` when those properties are copied and thus
are identical to the destination wrapper's properties (e.g. `Transpose{T,Array{T,N}}` has
the same dimensionality as the inner array). When creating an alias for this type, e.g.
`WrappedSomeArray{T,N} = WrappedArray{T,N,...}` the `Dst` typevar should typically be set to
`SomeArray{T,N}` while `Src` should be more lenient, e.g., `SomeArray`.

Only use this type for dispatch purposes. To convert instances of an array wrapper, use
[`adapt`](@ref).
"""
WrappedArray{T,N,Src,Dst} = Union{
      #SubArray{T,N,<:Src},
      Base.LogicalIndex{T,<:Src},
      PermutedDimsArray{T,N,<:Any,<:Any,<:Src},
      #Base.ReshapedArray{T,N,<:Src},
      #Base.ReinterpretArray{T,N,<:Any,<:Src},

      LinearAlgebra.Adjoint{T,<:Dst},
      LinearAlgebra.Transpose{T,<:Dst},
      LinearAlgebra.LowerTriangular{T,<:Dst},
      LinearAlgebra.UnitLowerTriangular{T,<:Dst},
      LinearAlgebra.UpperTriangular{T,<:Dst},
      LinearAlgebra.UnitUpperTriangular{T,<:Dst},
      LinearAlgebra.Diagonal{T,<:Dst},
      LinearAlgebra.Tridiagonal{T,<:Dst},
      LinearAlgebra.Symmetric{T,<:Dst},

      WrappedReinterpretArray{T,N,<:Src},
      WrappedReshapedArray{T,N,<:Src},
      WrappedSubArray{T,N,<:Src},
}

# XXX: this Union is a hack:
# - only works with one level of wrapping
# - duplication of Src and Dst typevars (without it, we get `WrappedGPUArray{T,N,AT{T,N}}`
#   not matching `SubArray{T,1,AT{T,2}}`, and leaving out `{T,N}` makes it impossible to
#   match e.g. `Diagonal{T,AT}` and get `N` out of that). alternatively, computed types
#   would make it possible to do `SubArray{T,N,<:AT.name.wrapper}` or `Diagonal{T,AT{T,N}}`.
#
# ideally, Base would have, e.g., `Transpose <: WrappedArray`, and we could use
# `Union{SomeArray, WrappedArray{<:Any, <:SomeArray}}` for dispatch.
# https://github.com/JuliaLang/julia/pull/31563

# accessors for extracting information about the wrapper type
ndims(::Type{<:Base.LogicalIndex}) = 1
ndims(::Type{<:LinearAlgebra.Adjoint}) = 2
ndims(::Type{<:LinearAlgebra.Transpose}) = 2
ndims(::Type{<:LinearAlgebra.LowerTriangular}) = 2
ndims(::Type{<:LinearAlgebra.UnitLowerTriangular}) = 2
ndims(::Type{<:LinearAlgebra.UpperTriangular}) = 2
ndims(::Type{<:LinearAlgebra.UnitUpperTriangular}) = 2
ndims(::Type{<:LinearAlgebra.Diagonal}) = 2
ndims(::Type{<:LinearAlgebra.Tridiagonal}) = 2
ndims(::Type{<:WrappedArray{<:Any,N}}) where {N} = N

eltype(::Type{<:WrappedArray{T}}) where {T} = T  # every wrapper has a T typevar

for T in [:(Base.LogicalIndex{<:Any,<:Src}),
          :(PermutedDimsArray{<:Any,<:Any,<:Any,<:Any,<:Src}),
          :(WrappedReinterpretArray{<:Any,<:Any,<:Src}),
          :(WrappedReshapedArray{<:Any,<:Any,<:Src}),
          :(WrappedSubArray{<:Any,<:Any,<:Src})]
    @eval begin
        parent(::Type{<:$T}) where {Src} = Src.name.wrapper
    end
end
parent(::Type{<:WrappedArray{<:Any,<:Any,<:Any,Dst}}) where {Dst} = Dst.name.wrapper


text:
using EnhancedGJK
using BenchmarkTools
using MeshIO, FileIO  # v0.3.2
using CoordinateTransformations, StaticArrays
import GeometryTypes: HyperRectangle, HyperSphere, Vec, Point, HomogenousMesh
using MeshCat, Colors

# Dimensions: (0.264 x 0.14 x 0.0885 m)
c1 = load(string(pwd(),"/test/meshes/r_foot_chull.obj"))
c2 = load(string(pwd(),"/test/meshes/r_foot_chull.obj"))
# trans = Translation(SVector{3,Float64}(rand(3,1)))
trans = Translation(SVector{3,Float64}(0,0,0.1))
poseA = trans
trans = Translation(SVector{3,Float64}(0,0,0))
poseB = trans

# Test
result = gjk(c1, c2, poseA, poseB)
ref_result = gjk_original(c1, c2, poseA, poseB)

println(result.termination, "\n", result.iterations)
println(ref_result.closest_point_in_body, "\n", ref_result.iterations)

println(separation_distance(result))
println(separation_distance(ref_result))

@btime gjk(c1, c2, poseA, poseB)
@btime gjk_original(c1, c2, poseA, poseB)

# Auto-differentiation through gjk
q = [0.,0.,0.1,0.,0.,0.]
∇f = ForwardDiff.gradient(gjk_diff, q)
println(∇f)

# Visualization
vis = Visualizer()
open(vis)

setobject!(vis["mesh/1"], c1)
setobject!(vis["mesh/2"], c2)
settransform!(vis["mesh/1"], poseA)
settransform!(vis["mesh/2"], poseB)

# Add arrows for the closest points

# Red is original
dist_vec1 = ArrowVisualizer(vis["dist/a"])
setobject!(dist_vec1, MeshLambertMaterial(color=RGBA{Float32}(1, 0, 0, 0.5)))
settransform!(dist_vec1, Point(poseA(ref_result.closest_point_in_body.a)),
    Point(poseB(ref_result.closest_point_in_body.b)))  # poseB transformation

# Green is SV
dist_vec2 = ArrowVisualizer(vis["dist/b"])
setobject!(dist_vec2, MeshLambertMaterial(color=RGBA{Float32}(0, 1, 0, 0.5)))
settransform!(dist_vec2, Point(poseA(result.closest_point_in_body.a)),
    Point(poseB(result.closest_point_in_body.b)))  # poseB transformation


text:
#TDS2002B scope driver

mutable struct TDS2002B
    handle
	initialized
    instr_dict
end

TDS2002B(handle) = TDS2002B(
                handle,
				false,
                Dict(
                    "CH1"=>"CH1",
                    "CH2"=>"CH2",
                    "Time"=>"HORizontal:MAIn:SCAle",
                    "Source"=>"TRIGger:MAIn:EDGE:SOUrce",
                    "Level"=>"TRIGger:MAIn:LEVel",
                    "Mode"=>"TRIGger:MAIn:MODe",
					"AUTO"=>"AUTO",
					"NORMAL"=>"NORMAL",
                    "Meas_ch"=>"SOUrce?",
                    "Meas_type"=>"TYPe?",
                    "Meas_value"=>"VALue?",
                    "MEAS1"=>1,
                    "MEAS2"=>2,
                    "MEAS3"=>3,
                    "MEAS4"=>4,
                    "MEAS5"=>5,
                    )
                    )

#Misc
function connect!(obj::TDS2002B, handle)
    obj.handle = handle
    obj.connected = 1
    return nothing
    end
function disconnect!(obj::TDS2002B)
    obj.handle = 0
    obj.connected = 0
    return nothing
    end

reset_instr(obj::TDS2002B) = write(obj.handle, "*RST")
clear_instr(obj::TDS2002B) = write(obj.handle, "*CLS")
get_idn(obj::TDS2002B) = query(obj.handle, "*IDN?")

#Get OSC Settings
function get_vertical_scale(obj::TDS2002B, ch)
    ch = obj.instr_dict[ch]
    cmd = "$ch:SCAle?"
    value = parse(Float64, strip(query(obj.handle, cmd)))
end

function get_horizontal_scale(obj::TDS2002B, time)
    time = obj.instr_dict[time]
    cmd = "$time?"
    value = parse(Float64, strip(query(obj.handle, cmd)))
end

function get_ch_position(obj::TDS2002B, ch)
    cmd = "$ch:POSition?"
    value = parse(Float64, strip(query(obj.handle, cmd)))
end

function set_ch_position(obj::TDS2002B, ch, lev)
    cmd = "$ch:POSition $lev"
    value = write(obj.handle, cmd)
end

# Channel coupling
function set_ch_coupling(obj::TDS2002B, ch, cpl)
	write(obj.handle, "$ch:COUPling $cpl")
end
function get_ch_coupling(obj::TDS2002B, ch)
	strip(query(obj.handle, "$ch:COUPling?"))
end

function get_trig_data(obj::TDS2002B, trig)
    trig = obj.instr_dict[trig]
    cmd = "$trig?"
    value = strip(query(obj.handle, cmd))
end

function get_meas_data(obj::TDS2002B, Meas_Nr1)
    nr = obj.instr_dict[Meas_Nr1]
    #cmd = "*OPC?"
    #value = strip(query(obj.handle, cmd))
    cmd = "MEASUrement:MEAS$nr:VALUE?"
    value = parse(Float64, strip(query(obj.handle, cmd)))
end

# set OSC settings low access functions
function set_vertical_scale(obj::TDS2002B, ch, value)
    ch = obj.instr_dict[ch]
    cmd = "$ch:SCAle $value"
    write(obj.handle, cmd)
end

function set_horizontal_scale(obj::TDS2002B, value)
    cmd = "HORizontal:MAIn:SCAle $value"
    write(obj.handle, cmd)
end

function set_trig_ch(obj::TDS2002B, ch)
    ch = obj.instr_dict[ch]
    cmd = "TRIGger:MAIn:EDGE:SOUrce $ch"
    write(obj.handle, cmd)
end

function set_trig_level(obj::TDS2002B, level)
    cmd = "TRIGger:MAIn:LEVel $level"
    write(obj.handle, cmd)
end

function set_trig_mode(obj::TDS2002B, mode)
	mode = obj.instr_dict[mode]
    cmd = "TRIGger:MAIn:MODe $mode"            # AUTO or NORMal
    write(obj.handle, cmd)
end

function set_meas_ch(obj::TDS2002B, Meas_Nr1, ch)
    nr = obj.instr_dict[Meas_Nr1]
    ch = obj.instr_dict[ch]
    cmd = "MEASUrement:MEAS$nr:SOUrce $ch"
    write(obj.handle, cmd)
end

function set_meas_type(obj::TDS2002B, Meas_Nr1, Type)
    nr = obj.instr_dict[Meas_Nr1]
    cmd = "MEASUrement:MEAS$nr:TYPe $Type"
    write(obj.handle, cmd)
end

function set_meas(obj::TDS2002B, Meas_Nr1, ch, Type)
    set_meas_ch(obj,Meas_Nr1, ch)
    set_meas_type(obj, Meas_Nr1, Type)
end

# select channel for acq and config
function conf_acq_ch(obj::TDS2002B, ch)
	ch = obj.instr_dict[ch]
    write(obj.handle, "SELect:$ch ON")
	write(obj.handle, "DATa:SOUrce $ch;DATa:ENCdg SRIBINARY;DATa:WIDth 1;DATa:STARt 1;DATa:STOP 2500")
	write(obj.handle, "CURVE?")
end

function auto_set(obj::TDS2002B)
	write(obj.handle, "AUTOSet EXECute")
	@info "Scope auto-set in progress"
	sleep(1)
end

function Trigger_Aquistion(obj::TDS2002B, ch)
	ch = obj.instr_dict[ch]

	y_buffer = zeros(UInt8, 2506)
	viRead!(obj.handle, y_buffer)

	gain = parse(Float64, query(obj.handle, "WFMPre:YMUlt?"))
	#offs = parse(Float64, query(obj.handle, "WFMPre:YOFf?")) # THIS IS A PIECE OF SHIT COMMAND!!!!
	offs = parse(Float64, query(obj.handle, "$ch:POSition?"))
	dt = parse(Float64, query(obj.handle, "WFMPre:XINcr?"))

	y = Float64.(reinterpret(Int8, y_buffer[7:end])) .* gain .- offs
	t = collect(0:dt:dt*(length(y)-1))

	return t, y
end

max_stars_repo_name: iuliancioarca/TIVM

max_issues_repo_name: iuliancioarca/TIVM

max_forks_repo_name: iuliancioarca/TIVM


text:
using Test, StaticRanges, Dates
using StaticRanges: can_set_first, can_set_last, can_set_step, has_step, can_set_length,
    stephi, steplo, refhi, reflo, eqmax, eqmin, ltmax, ltmin, gtmax, gtmin, group_max,
    group_min, min_of_group_max, max_of_group_min, ordmin, ordmax, next_type, prev_type,
    Unordered, set_ref!, set_offset!, set_lendiv!

# Uniqueness methods
using StaticRanges: all_unique, Uniqueness, AllUnique, NotUnique, UnkownUnique

using StaticRanges: ArrayInterface.ismutable

using Base: OneTo, step_hp
using Base.Order

include("continuity_tests.jl")
include("uniqueness_tests.jl")
include("order_tests.jl")
include("unitrange_tests.jl")
include("twiceprecision.jl")
include("mutate.jl")
include("find.jl")
include("range_interface.jl")
include("promotion.jl")
include("broadcast.jl")
include("onetorange.jl")
include("unitrange_tests.jl")
include("steprange_tests.jl")
include("intersect.jl")
include("reverse.jl")

include("steprangelen_test.jl")

include("linrange_test.jl")
include("size_tests.jl")
include("step_tests.jl")
include("length_tests.jl")

for frange in (mrange, srange)
    @testset "$frange" begin
        @test_throws ArgumentError frange(1)
        @test_throws ArgumentError frange(1, step=1)
        @test_throws ArgumentError frange(nothing)
        @test_throws ArgumentError frange(nothing, length=1)
        @test_throws ArgumentError frange(1, step=1, stop=1, length=1)
        if frange == mrange
            # cannot infer a static parameter from construction
            @testset "colon" begin
                @inferred(frange(10, step=1, stop=0))
                @inferred(frange(1, step=.2, stop=2))
                @inferred(frange(1., step=.2, stop=2.))
                @inferred(frange(2, step=-.2, stop=1))
                @inferred(frange(1, 0))
                @inferred(frange(0.0, -0.5))
            end
        end
        @testset "indexing" begin
            L32 = frange(Int32(1), stop=Int32(4), length=4)
            L64 = frange(Int64(1), stop=Int64(4), length=4)
            @test @inferred(L32[1]) === 1.0 && @inferred(L64[1]) === 1.0
            @test L32[2] == 2 && L64[2] == 2
            @test L32[3] == 3 && L64[3] == 3
            @test L32[4] == 4 && L64[4] == 4
            @test frange(1.0, stop=2.0, length=2)[1] === 1.0
            @test frange(1.0f0, stop=2.0f0, length=2)[1] === 1.0f0
            @test frange(Float16(1.0), stop=Float16(2.0), length=2)[1] === Float16(1.0)

            let r = frange(5, step=-1, stop=1)
                @test r[1]==5
                @test r[2]==4
                @test r[3]==3
                @test r[4]==2
                @test r[5]==1
            end
            @test @inferred(frange(0.1, step=0.1, stop=0.3)[2]) === 0.2
            @test @inferred(frange(0.1f0, step=0.1f0, stop=0.3f0)[2]) === 0.2f0

            @test frange(1, 5)[1:4] == 1:4
            @test frange(1.0, 5)[1:4] == 1.0:4
            @test frange(2, 6)[1:4] == 2:5
            @test frange(1, 6)[2:5] == 2:5
            @test frange(1, 6)[2:2:5] == 2:2:4
            @test frange(1, step=2, stop=13)[2:6] == 3:2:11
            @test frange(1, step=2, stop=13)[2:3:7] == 3:6:13

            @test isempty(frange(1, 4)[5:4])
            #@test_throws BoundsError frange(1:10)[8:-1:-2]
        end
       @testset "keys/values" begin
            keytype_is_correct(r) = keytype(r) == eltype(keys(r))
            valtype_is_correct(r) = valtype(r) == eltype(values(r))
            @test keytype_is_correct(frange(1, 3))
            @test keytype_is_correct(frange(1, step=.3, stop=4))
            @test keytype_is_correct(frange(.11, step=.1, stop=.3))
            @test keytype_is_correct(frange(Int8(1), Int8(5)))
            @test keytype_is_correct(frange(Int16(1), Int8(5)))
            @test keytype_is_correct(frange(Int16(1), step=Int8(3), stop=Int8(5)))
            @test keytype_is_correct(frange(Int8(1), step=Int16(3), stop=Int8(5)))
            @test keytype_is_correct(frange(Int8(1), step=Int8(3), stop=Int16(5)))
            @test keytype_is_correct(frange(Int64(1), Int64(5)))
            @test keytype_is_correct(frange(Int64(1), Int64(5)))
            @test keytype_is_correct(frange(Int128(1), Int128(5)))
            @test valtype_is_correct(1:3)
            @test valtype_is_correct(1:.3:4)
            @test valtype_is_correct(.1:.1:.3)
            @test valtype_is_correct(frange(Int8(1), Int8(5)))
            @test valtype_is_correct(frange(Int16(1), Int8(5)))
            @test valtype_is_correct(frange(Int16(1), step=Int8(3), stop=Int8(5)))
            @test valtype_is_correct(frange(Int8(1), step=Int16(3), stop=Int8(5)))
            @test valtype_is_correct(frange(Int8(1), step=Int8(3), stop=Int16(5)))
            @test valtype_is_correct(frange(Int64(1), Int64(5)))
            @test valtype_is_correct(frange(Int64(1), Int64(5)))
            @test valtype_is_correct(frange(Int128(1), Int128(5)))

            if frange isa typeof(mrange)
                @test keytype_is_correct(OneToMRange(4))
                @test keytype_is_correct(OneToMRange(Int32(4)))
                @test valtype_is_correct(OneToMRange(4))
                @test valtype_is_correct(OneToMRange(Int32(4)))
            else
                @test keytype_is_correct(OneToSRange(4))
                @test keytype_is_correct(OneToSRange(Int32(4)))
                @test valtype_is_correct(OneToSRange(4))
                @test valtype_is_correct(OneToSRange(Int32(4)))
            end
        end

        @testset "findall(::Base.Fix2{typeof(in)}, ::Array)" begin
            @test findall(in(3:20), [5.2, 3.3]) == findall(in(Vector(3:20)), [5.2, 3.3])

            let span = frange(5, 20),
                r = frange(-7, step=3, stop=42)
                @test findall(in(span), r) == 5:10

                r = frange(15, step=-2, stop=-38)
                @test findall(in(span), r) == 6:-1:1
            end
        end

        if VERSION > v"1.2"
            @testset "issubset" begin
                @test issubset(frange(1, 3), 1:typemax(Int)) #32461
                @test issubset(frange(1, 3), 1:3)
                @test issubset(frange(1, 3), 1:4)
                @test issubset(frange(1, 3), 0:3)
                @test issubset(frange(1, 3), 0:4)
                @test !issubset(frange(1, 5), 2:5)
                @test !issubset(frange(1, 5), 1:4)
                @test !issubset(frange(1, 5), 2:4)
                @test issubset(frange(1, step=3, stop=10), 1:10)
                @test !issubset(frange(1, 10), 1:3:10)

                @test issubset(1:3, frange(1, typemax(Int))) #32461
                @test issubset(1:3, frange(1, 3))
                @test issubset(1:3, frange(1, 4))
                @test issubset(1:3, frange(0, 3))
                @test issubset(1:3, frange(0, 4))
                @test !issubset(1:5, frange(2, 5))
                @test !issubset(1:5, frange(1, 4))
                @test !issubset(1:5, frange(2, 4))
                @test issubset(1:3:10, frange(1, 10))
                @test !issubset(1:10, frange(1, step=3, stop=10))
                if frange isa typeof(mrange)
                    @test issubset(OneToMRange(5), OneToMRange(10))
                    @test !issubset(OneToMRange(10), OneToMRange(5))
                    @test issubset(OneToMRange(5), OneToMRange(10))
                    @test !issubset(OneToMRange(10), OneToMRange(5))
                else
                    @test issubset(OneToSRange(5), OneToSRange(10))
                    @test !issubset(OneToSRange(10), OneToSRange(5))
                    @test issubset(OneToSRange(5), OneToSRange(10))
                    @test !issubset(OneToSRange(10), OneToSRange(5))
                end

            end
        end

        @testset "in" begin
            @test 0 in frange(UInt(0), step=100, stop=typemax(UInt))
            @test last(frange(UInt(0), step=100, stop=typemax(UInt))) in frange(UInt(0), step=100, stop=typemax(UInt))
            @test -9223372036854775790 in frange(-9223372036854775790, step=100, stop=9223372036854775710)
            @test -9223372036854775690 in frange(-9223372036854775790, step=100, stop=9223372036854775710)
            @test -90 in frange(-9223372036854775790, step=100, stop=9223372036854775710)
            @test 10 in frange(-9223372036854775790, step=100, stop=9223372036854775710)
            @test 110 in frange(-9223372036854775790, step=100, stop=9223372036854775710)
            @test 9223372036854775610 in frange(-9223372036854775790, step=100, stop=9223372036854775710)
            @test 9223372036854775710 in frange(-9223372036854775790, step=100, stop=9223372036854775710)

            @test !(3.5 in frange(1, 5))
            @test (3 in frange(1, 5))
            @test (3 in frange(5, step=-1, stop=1))
            #@test (3 in 3+0*(1:5))
            #@test !(4 in 3+0*(1:5))

            let r = frange(0.0, step=0.01, stop=1.0)
                @test (r[30] in r)
            end
            let r = frange((-4*Int64(maxintfloat(Int === Int32 ? Float32 : Float64))), 5)
                @test (3 in r)
                @test (3.0 in r)
            end

            @test !(1 in frange(1, 0))
            @test !(1.0 in frange(1.0, 0.0))
        end

        @testset "in() works across types, including non-numeric types (#21728)" begin
            @test 1//1 in frange(1, 3)
            @test 1//1 in frange(1.0, 3.0)
            @test !(5//1 in frange(1, 3))
            @test !(5//1 in frange(1.0, 3.0))
            @test Complex(1, 0) in frange(1, 3)
            @test Complex(1, 0) in frange(1.0, 3.0)
            @test Complex(1.0, 0.0) in frange(1, 3)
            @test Complex(1.0, 0.0) in frange(1.0, 3.0)
            @test !(Complex(1, 1) in frange(1, 3))
            @test !(Complex(1, 1) in frange(1.0, 3.0))
            @test !(Complex(1.0, 1.0) in frange(1, 3))
            @test !(Complex(1.0, 1.0) in frange(1.0, 3.0))
            @test !(π in frange(1, 3))
            @test !(π in frange(1.0, 3.0))
            @test !("a" in frange(1, 3))
            @test !("a" in frange(1.0, 3.0))
            @test !(1 in frange(Date(2017, 01, 01), step=Dates.Day(1), stop=Date(2017, 01, 05)))
            @test !(Complex(1, 0) in frange(Date(2017, 01, 01), step=Dates.Day(1), stop=Date(2017, 01, 05)))
            @test !(π in frange(Date(2017, 01, 01), step=Dates.Day(1), stop=Date(2017, 01, 05)))
            @test !("a" in frange(Date(2017, 01, 01), step=Dates.Day(1), stop=Date(2017, 01, 05)))
        end
        @testset "sums of ranges" begin
            @test sum(frange(1, 100)) == 5050
            @test sum(frange(0, 100)) == 5050
            @test sum(frange(-100, 100)) == 0
            @test sum(frange(0, step=2, stop=100)) == 2550
        end
        @testset "Tricky sums of StepRangeLen #8272" begin
            @test sum(frange(10000., step=-0.0001, stop=0)) == 5.00000005e11
            @test sum(frange(0, step=0.001, stop=1)) == 500.5
            @test sum(frange(0, step=0.000001, stop=1)) == 500000.5
            @test sum(frange(0, step=0.1, stop=10)) == 505.
        end
        @testset "operations between ranges and arrays" begin
            @test all(([frange(1, 5);] + (frange(5, step=-1, stop=1))) .== 6)
            @test all(((frange(5, step=-1, stop=1)) + [frange(1, 5);]) .== 6)
            @test all(([frange(1, 5);] - (frange(1, 5))) .== 0)
            @test all((frange(1, 5) - [frange(1, 5);]) .== 0)
        end

        @testset "broadcasted operations with scalars" begin
            @test broadcast(-, frange(1, 3)) == -1:-1:-3
            @test broadcast(-, frange(1, 3), 2) == -1:1
            @test broadcast(-, frange(1, 3), 0.25) == 1-0.25:3-0.25
            @test broadcast(+, frange(1, 3)) == 1:3
            @test broadcast(+, frange(1, 3), 2) == 3:5
            @test broadcast(+, frange(1, 3), 0.25) == 1+0.25:3+0.25
            @test broadcast(+, frange(1, step=2, stop=6), 1) == 2:2:6
            @test broadcast(+, frange(1, step=2, stop=6), 0.3) == 1+0.3:2:5+0.3
            @test broadcast(-, frange(1, step=2, stop=6), 1) == 0:2:4
            @test broadcast(-, frange(1, step=2, stop=6), 0.3) == 1-0.3:2:5-0.3
            @test broadcast(-, 2, frange(1, 3)) == 1:-1:-1
        end

        @testset "loops involving typemin/typemax" begin
            n = 0
            s = 0
            # loops ending at typemax(Int)
            for i = frange((typemax(Int)-1), typemax(Int))
                s += 1
                @test s <= 2
            end
            @test s == 2

            s = 0
            for i = frange((typemax(Int)-2), (typemax(Int)-1))
                s += 1
                @test s <= 2
            end
            @test s == 2

            s = 0
            for i = frange(typemin(Int), (typemin(Int)+1))
                s += 1
                @test s <= 2
            end
            @test s == 2

            # loops covering the full range of integers
            s = 0
            for i = frange(typemin(UInt8), typemax(UInt8))
                s += 1
            end
            @test s == 256

            s = 0
            for i = frange(typemin(UInt), typemax(UInt))
                i == 10 && break
                s += 1
            end
            @test s == 10

            s = 0
            for i = frange(typemin(UInt8), step=one(UInt8), stop=typemax(UInt8))
                s += 1
            end
            @test s == 256

            s = 0
            for i = frange(typemin(UInt), step=1, stop=typemax(UInt))
                i == 10 && break
                s += 1
            end
            @test s == 10

            # loops past typemax(Int)
            n = 0
            s = Int128(0)
            for i = frange(typemax(UInt64)-2, typemax(UInt64))
                n += 1
                s += i
            end
            @test n == 3
            @test s == 3*Int128(typemax(UInt64)) - 3

            # loops over empty ranges
            s = 0
            for i = 0xff:0x00
                s += 1
            end
            @test s == 0

            s = 0
            for i = frange(Int128(typemax(Int128)), Int128(typemin(Int128)))
                s += 1
            end
            @test s == 0
        end

        function range_fuzztests(::Type{T}, niter, nrange) where {T}
            for i = 1:niter, n in nrange
                strt, Δ = randn(T), randn(T)
                Δ == 0 && continue
                stop = strt + (n-1)*Δ
                # `n` is not necessarily unique s.t. `strt + (n-1)*Δ == stop`
                # so test that `length(strt:Δ:stop)` satisfies this identity
                # and is the closest value to `(stop-strt)/Δ` to do so
                lo = hi = n
                while strt + (lo-1)*Δ == stop; lo -= 1; end
                while strt + (hi-1)*Δ == stop; hi += 1; end
                m = clamp(round(Int, (stop-strt)/Δ) + 1, lo+1, hi-1)
                r = strt:Δ:stop
                @test m == length(r)
                @test strt == first(r)
                @test Δ == step(r)
                @test_skip stop == last(r)
                l = range(strt, stop=stop, length=n)
                @test n == length(l)
                @test strt == first(l)
                @test stop  == last(l)
            end
        end

        @testset "range fuzztests for $T" for T = (Float32, Float64,)
            range_fuzztests(T, 2^15, frange(1, 5))
        end

        @testset "range with very large endpoints for type $T" for T = (Float32, Float64)
            largeint = Int(min(maxintfloat(T), typemax(Int)))
            a = floatmax()
            for i = 1:5
                @test [frange(a, stop=a, length=1);] == [a]
                @test [frange(-a, stop=-a, length=1);] == [-a]
                b = floatmax()
                for j = 1:5
                    @test [frange(-a, stop=b, length=0);] == []
                    @test [frange(-a, stop=b, length=2);] == [-a,b]
                    @test [frange(-a, stop=b, length=3);] == [-a,(b-a)/2,b]
                    @test [frange(a, stop=-b, length=0);] == []
                    @test [frange(a, stop=-b, length=2);] == [a,-b]
                    @test [frange(a, stop=-b, length=3);] == [a,(a-b)/2,-b]
                    for c = largeint-3:largeint
                        s = range(-a, stop=b, length=c)
                        @test first(s) == -a
                        @test last(s) == b
                        @test length(s) == c
                        s = range(a, stop=-b, length=c)
                        @test first(s) == a
                        @test last(s) == -b
                        @test length(s) == c
                    end
                    b = prevfloat(b)
                end
                a = prevfloat(a)
            end
        end

        @testset "ranges with very small endpoints for type $T" for T = (Float32, Float64)
            z = zero(T)
            u = eps(z)
            @test first(frange(u, stop=u, length=0)) == u
            @test last(frange(u, stop=u, length=0)) == u
            @test first(frange(-u, stop=u, length=0)) == -u
            @test last(frange(-u, stop=u, length=0)) == u
            @test [frange(-u, stop=u, length=0);] == []
            @test [frange(-u, stop=-u, length=1);] == [-u]
            @test [frange(-u, stop=u, length=2);] == [-u,u]
            @test [frange(-u, stop=u, length=3);] == [-u,0,u]
            @test first(frange(-u, stop=-u, length=0)) == -u
            @test last(frange(-u, stop=-u, length=0)) == -u
            @test first(frange(u, stop=-u, length=0)) == u
            @test last(frange(u, stop=-u, length=0)) == -u
            @test [frange(u, stop=-u, length=0);] == []
            @test [frange(u, stop=u, length=1);] == [u]
            @test [frange(u, stop=-u, length=2);] == [u,-u]
            @test [frange(u, stop=-u, length=3);] == [u,0,-u]
            v = frange(-u, stop=u, length=12)
            @test length(v) == 12
            @test [-3u:u:3u;] == [frange(-3u, stop=3u, length=7);] == [-3:3;].*u
            @test [3u:-u:-3u;] == [frange(3u, stop=-3u, length=7);] == [3:-1:-3;].*u
        end

        @testset "tricky floating-point ranges" begin
            for (start, step, stop, len) in ((1, 1, 3, 3),
                                             (0, 1, 3, 4),
                                             (3, -1, -1, 5),
                                             (1, -1, -3, 5),
                                             (0, 1, 10, 11),
                                             (0, 7, 21, 4),
                                             (0, 11, 33, 4),
                                             (1, 11, 34, 4),
                                             (0, 13, 39, 4),
                                             (1, 13, 40, 4),
                                             (11, 11, 33, 3),
                                             (3, 1, 11, 9),
                                             (0, 10, 55, 0),
                                             (0, -1, 5, 0),
                                             (0, 10, 5, 0),
                                             (0, 1, 5, 0),
                                             (0, -10, 5, 0),
                                             (0, -10, 0, 1),
                                             (0, -1, 1, 0),
                                             (0, 1, -1, 0),
                                             (0, -1, -10, 11))
                r = frange(start/10, step=step/10, stop=stop/10)
                a = Vector(frange(start, step=step, stop=stop))./10
                ra = Vector(r)

                @test r == a
                @test isequal(r, a)

                @test r == ra
                @test isequal(r, ra)

                @test hash(r) == hash(a)
                @test hash(r) == hash(ra)

                if len > 0
                    l = frange(start/10, stop=stop/10, length=len)
                    la = Vector(l)

                    @test a == l
                    @test r == l
                    @test isequal(a, l)
                    @test isequal(r, l)

                    @test l == la
                    @test isequal(l, la)

                    @test hash(l) == hash(a)
                    @test hash(l) == hash(la)
                end
            end

            @test 1.0:1/49:27.0 == range(1.0, stop=27.0, length=1275) == [49:1323;]./49
            @test isequal(1.0:1/49:27.0, range(1.0, stop=27.0, length=1275))
            @test isequal(frange(1.0, step=1/49, stop=27.0), Vector(49:1323)./49)
            @test hash(1.0:1/49:27.0) == hash(frange(1.0, stop=27.0, length=1275)) == hash(Vector(frange(49, 1323))./49)

            @test [frange(prevfloat(0.1), step=0.1, stop=0.3);] == [prevfloat(0.1), 0.2, 0.3]
            @test [frange(nextfloat(0.1), step=0.1, stop=0.3);] == [nextfloat(0.1), 0.2]
            @test [frange(prevfloat(0.0), step=0.1, stop=0.3);] == [prevfloat(0.0), 0.1, 0.2]
            @test [frange(nextfloat(0.0), step=0.1, stop=0.3);] == [nextfloat(0.0), 0.1, 0.2]
            @test [frange(0.1, step=0.1, stop=prevfloat(0.3));] == [0.1, 0.2]
            @test [frange(0.1, step=0.1, stop=nextfloat(0.3));] == [0.1, 0.2, nextfloat(0.3)]
            @test [frange(0.0, step=0.1, stop=prevfloat(0.3));] == [0.0, 0.1, 0.2]
            @test [frange(0.0, step=0.1, stop=nextfloat(0.3));] == [0.0, 0.1, 0.2, nextfloat(0.3)]
            @test [frange(0.1, step=prevfloat(0.1), stop=0.3);] == [0.1, 0.2, 0.3]
            @test [frange(0.1, step=nextfloat(0.1), stop=0.3);] == [0.1, 0.2]
            @test [frange(0.0, step=prevfloat(0.1), stop=0.3);] == [0.0, prevfloat(0.1), prevfloat(0.2), 0.3]
            @test [frange(0.0, step=nextfloat(0.1), stop=0.3);] == [0.0, nextfloat(0.1), nextfloat(0.2)]
        end

        @testset "overflowing sums (see #5798)" begin
            if Sys.WORD_SIZE == 64
                @test sum(frange(Int128(1), 10^18)) == div(10^18 * (Int128(10^18)+1), 2)
                @test sum(frange(Int128(1), 10^18-1)) == div(10^18 * (Int128(10^18)-1), 2)
            else
                @test sum(frange(Int64(1), 10^9)) == div(10^9 * (Int64(10^9)+1), 2)
                @test sum(frange(Int64(1), 10^9-1)) == div(10^9 * (Int64(10^9)-1), 2)
            end
        end

        @testset "issue #20373 (unliftable ranges with exact end points)" begin
            @test [3*0.05:0.05:0.2;]    == [frange(3*0.05, stop=0.2, length=2);]   == [3*0.05,0.2]
            @test [0.2:-0.05:3*0.05;]   == [frange(0.2, stop=3*0.05, length=2);]   == [0.2,3*0.05]
            @test [-3*0.05:-0.05:-0.2;] == [frange(-3*0.05, stop=-0.2, length=2);] == [-3*0.05,-0.2]
            @test [-0.2:0.05:-3*0.05;]  == [frange(-0.2, stop=-3*0.05, length=2);] == [-0.2,-3*0.05]
        end

        @testset "issue #6973" begin
            r1 = frange(1.0, step=0.1, stop=2.0)
            r2 = frange(1.0f0, step=0.2f0, stop=3.0f0)
            r3 = frange(1, step=2, stop=21)
            @test r1 + r1 == 2*r1
            @test r1 + r2 == 2.0:0.3:5.0
            @test (r1 + r2) - r2 == r1
            @test r1 + r3 == convert(StepRangeLen{Float64}, r3) + r1
            @test r3 + r3 == 2 * r3
        end

       @testset "Inexact errors on 32 bit architectures. #22613" begin
            @test first(frange(log(0.2), stop=log(10.0), length=10)) == log(0.2)
            @test last(frange(log(0.2), stop=log(10.0), length=10)) == log(10.0)
            # not used internally for StaticRanges
            #@test length(Base.floatrange(-3e9, 1.0, 1, 1.0)) == 1
        end

        # issue #7426
        @test [frange(typemax(Int), step=1, stop=typemax(Int));] == [typemax(Int)]

        @testset "issue #7387" begin
            for r in (frange(0, 1), frange(0.0, 1.0))
                local r
                @test [r .+ im;] == [r;] .+ im
                @test [r .- im;] == [r;] .- im
                @test [r * im;] == [r;] * im
                @test [r / im;] == [r;] / im
            end
        end

        # near-equal ranges
        @test frange(0.0, step=0.1, stop=1.0) != 0.0f0:0.1f0:1.0f0

        #issue #7484
        let r7484 = frange(0.1, step=0.1, stop=1)
            @test [reverse(r7484);] == reverse([r7484;])
        end

        # issue #2959
        @test frange(1.0, 1.5) == 1.0:1.0:1.5 == 1.0:1.0

        @testset "comparing UnitRanges and OneTo" begin
            @test frange(1, step=2, stop=10) == 1:2:10 != 1:3:10 != 1:3:13 != frange(2, step=3, stop=13) == 2:3:11 != frange(2, 11)
            @test frange(1, step=1, stop=10) == 1:10 == 1:10 == OneToMRange(10) == OneToSRange(10)
            @test 1:10 != frange(2, 10) != 2:11 != Base.OneTo(11)
            @test OneToMRange(10) != OneToSRange(11) != frange(1, 10)
        end

        @testset "issue #7114" begin
            let r = frange(-0.004532318104333742, step=1.2597349521122731e-5, stop=0.008065031416788989)
                @test length(r[1:end-1]) == length(r) - 1
                @test isa(r[1:2:end],AbstractRange) && length(r[1:2:end]) == div(length(r)+1, 2)
                @test r[3:5][2] ≈ r[4]
                @test r[5:-2:1][2] ≈ r[3]
                @test_throws BoundsError r[0:10]
                @test_throws BoundsError r[1:10000]
            end

            let r = frange(1/3, stop=5/7, length=6)
                @test length(r) == 6
                @test r[1] == 1/3
                @test abs(r[end] - 5/7) <= eps(5/7)
            end

            let r = frange(0.25, stop=0.25, length=1)
                @test length(r) == 1
                @test_throws ArgumentError frange(0.25, stop=0.5, length=1)
            end
        end

        # Preservation of high precision upon addition
        let r = range(-0.1, step=0.1, stop=0.3) + broadcast(+, -0.3:0.1:0.1, 1e-12)
            @test r[3] == 1e-12
        end

        @testset "range with 1 or 0 elements (whose step length is NaN)" begin
            @test issorted(frange(1, stop=1, length=0))
            @test issorted(frange(1, stop=1, length=1))
        end
        @testset "indexing range with empty range (#4309)" begin
            @test frange(3, 6)[5:4] == 7:6
            @test_throws BoundsError frange(3, 6)[5:5]
            @test_throws BoundsError frange(3, 6)[5]
            @test frange(0, step=2, stop=10)[7:6] == 12:2:10
            @test_throws BoundsError (0:2:10)[7:7]
        end
        # indexing with negative ranges (#8351)
        for a=AbstractRange[3:6, frange(0, step=2, stop=10)], b=AbstractRange[frange(0, 1), frange(2, step=-1, stop=0)]
            @test_throws BoundsError a[b]
        end
        @testset "sort/sort!/partialsort" begin
            @test sort(frange(1, 2)) == UnitRange(1,2)
            @test sort!(frange(1, 2)) == UnitRange(1,2)
            @test sort(frange(1, 10), rev=true) == 10:-1:1
            @test sort(frange(-3, 3), by=abs) == [0,-1,1,-2,2,-3,3]
            @test partialsort(frange(1, 10), 4) == 4
        end

        if frange == mrange
            # this works with `srange` but takes way too long to run.
            function loop_range_values(::Type{T}) where T
                for a = -5:25,
                    s = [-5:-1; 1:25; ],
                    d = 1:25,
                    n = -1:15

                    denom = convert(T, d)
                    strt = convert(T, a)/denom
                    Δ     = convert(T, s)/denom
                    stop  = convert(T, (a + (n - 1) * s)) / denom
                    vals  = T[a:s:(a + (n - 1) * s); ] ./ denom
                    r = strt:Δ:stop
                    @test [r;] == vals
                    @test [frange(strt, stop=stop, length=length(r));] == vals
                    n = length(r)
                    @test [r[1:n];] == [r;]
                    @test [r[2:n];] == [r;][2:end]
                    @test [r[1:3:n];] == [r;][1:3:n]
                    @test [r[2:2:n];] == [r;][2:2:n]
                    @test [r[n:-1:2];] == [r;][n:-1:2]
                    @test [r[n:-2:1];] == [r;][n:-2:1]
                end
            end

            @testset "issue #7420 for type $T" for T = (Float32, Float64,) # BigFloat),
                loop_range_values(T)
            end
        end

        # comparing and hashing ranges
        @testset "comparing and hashing ranges" begin
            Rs = AbstractRange[]
            for r in (frange(1, 1),
                      frange(1, step=1, stop=1),
                      frange(1, 2),
                      frange(1, step=1, stop=2),
                      map(Int32, frange(1, step=3, stop=17)),
                      map(Int64, frange(1, step=3, stop=17)),
                      frange(1, 0),
                      frange(1, step=-1, stop=0),
                      frange(17, 0, step=-3),
                      frange(0.0, step=0.1, stop=1.0),
                      map(Float32, frange(0.0, step=0.1, stop=1.0)),
                      map(Float32, ifelse(frange == mrange, LinMRange(0.0, 1.0, 11), LinSRange(0.0, 1.0, 11))),
                      frange(1.0, step=eps(), stop=1.0) .+ 10eps(),
                      frange(9007199254740990., step=1.0, stop=9007199254740994),
                      frange(0, stop=1, length=20),
                      map(Float32, frange(0, stop=1, length=20)))
                local r
                ar = Vector(r)
                @test r == ar
                @test isequal(r,ar)
                @test hash(r) == hash(ar)
                for s in Rs
                    as = Vector(s)
                    @test isequal(r,s) == (hash(r)==hash(s))
                    @test (r==s) == (ar==as)
                end
            end
        end
    end
end

@testset "AbstractStepRangeLen" begin
    @test StaticRanges.floatmrange(1.0, 1.0, 10, 2.0) == Base.floatrange(1.0, 1.0, 10, 2.0)
    @test StaticRanges.floatsrange(1.0, 1.0, 10, 2.0) == Base.floatrange(1.0, 1.0, 10, 2.0)
end
#=
#@test 1.0:(.3-.1)/.1 == 1.0:2.0
=#