text:
#=

 Unit tests from libieeep1788 for reverse interval operations
 (Original author: Marco Nehmeier)
 converted into portable ITL format by Oliver Heimlich.

 Copyright 2013-2015 Marco Nehmeier (nehmeier@informatik.uni-wuerzburg.de)
 Copyright 2015-2017 Oliver Heimlich (oheim@posteo.de)

 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at

     http://www.apache.org/licenses/LICENSE-2.0

 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.

=#
#Language imports

#Test library imports
using Test

#Arithmetic library imports
using IntervalArithmetic

#Preamble
setprecision(53)
setprecision(Interval, Float64)
setrounding(Interval, :tight)
# Set full format, and show decorations
@format full

@testset "minimal_sqr_rev_test" begin

end

@testset "minimal_sqr_rev_bin_test" begin

end

@testset "minimal_sqr_rev_dec_test" begin

end

@testset "minimal_sqr_rev_dec_bin_test" begin

end

@testset "minimal_abs_rev_test" begin

end

@testset "minimal_abs_rev_bin_test" begin

end

@testset "minimal_abs_rev_dec_test" begin

end

@testset "minimal_abs_rev_dec_bin_test" begin

end

@testset "minimal_pown_rev_test" begin

end

@testset "minimal_pown_rev_bin_test" begin

end

@testset "minimal_pown_rev_dec_test" begin

end

@testset "minimal_pown_rev_dec_bin_test" begin

end

@testset "minimal_sin_rev_test" begin

end

@testset "minimal_sin_rev_bin_test" begin

end

@testset "minimal_sin_rev_dec_test" begin

end

@testset "minimal_sin_rev_dec_bin_test" begin

end

@testset "minimal_cos_rev_test" begin

end

@testset "minimal_cos_rev_bin_test" begin

end

@testset "minimal_cos_rev_dec_test" begin

end

@testset "minimal_cos_rev_dec_bin_test" begin

end

@testset "minimal_tan_rev_test" begin

end

@testset "minimal_tan_rev_bin_test" begin

end

@testset "minimal_tan_rev_dec_test" begin

end

@testset "minimal_tan_rev_dec_bin_test" begin

end

@testset "minimal_cosh_rev_test" begin

end

@testset "minimal_cosh_rev_bin_test" begin

end

@testset "minimal_cosh_rev_dec_test" begin

end

@testset "minimal_cosh_rev_dec_bin_test" begin

end

@testset "minimal_mul_rev_test" begin

end

@testset "minimal_mul_rev_ten_test" begin

end

@testset "minimal_mul_rev_dec_test" begin

end

@testset "minimal_mul_rev_dec_ten_test" begin

end
# FactCheck.exitstatus()


text:
# 2bit Record
# ===========

mutable struct Record
    filled::Bool
    dnasize::UInt32
    blockcount::UInt32
    blockstarts::Vector{UInt32}
    blocksizes::Vector{UInt32}
    maskedblockcount::UInt32
    maskedblockstarts::Vector{UInt32}
    maskedblocksizes::Vector{UInt32}
    reserved::UInt32
    packeddna::Vector{UInt8}
end

"""
    TwoBit.Record()

Create an unfilled 2bit record.
"""
function Record()
    return Record(
        false,
        # dnasize-blocksizes
        0, 0, UInt32[], UInt32[],
        # maskedblockcount-reserved
        0, UInt32[], UInt32[], 0,
        # packeddna
        UInt8[])
end

function initialize!(record::Record)
    record.filled = false
    record.dnasize = 0
    record.blockcount = 0
    empty!(record.blockstarts)
    empty!(record.blocksizes)
    record.maskedblockcount = 0
    empty!(record.maskedblockstarts)
    empty!(record.maskedblocksizes)
    record.reserved = 0
    return record
end

function isfilled(record::Record)
    return record.filled
end

function Base.:(==)(record1::Record, record2::Record)
    if isfilled(record1) == isfilled(record2) == true
        return (
            record1.dnasize           == record2.dnasize           &&
            record1.blockcount        == record2.blockcount        &&
            record1.blockstarts       == record2.blockstarts       &&
            record1.blocksizes        == record2.blocksizes        &&
            record1.maskedblockcount  == record2.maskedblockcount  &&
            record1.maskedblockstarts == record2.maskedblockstarts &&
            record1.maskedblocksizes  == record2.maskedblocksizes  &&
            record1.reserved          == record2.reserved          &&
            memcmp(pointer(record1.packeddna), pointer(record2.packeddna), cld(record1.dnasize, 4)) == 0)
    else
        return isfilled(record1) == isfilled(record2) == false
    end
end

function Base.show(io::IO, record::Record)
    print(io, summary(record), ':')
    if isfilled(record)
        println(io)
        println(io, "         length: ", record.dnasize)
        println(io, "       sequence: ", truncate(sequence(record), 40))
        println(io, "       N blocks: ", record.blockcount)
          print(io, "  masked blocks: ", record.maskedblockcount)
    else
        print(io, " <not filled>")
    end
end

function truncate(seq, len)
    if length(seq) > len
        return "$(seq[1:len-1])…"
    else
        return string(seq)
    end
end

# Accessors
# ---------

function hassequence(record::Record)
    return isfilled(record)
end

"""
    sequence([::Type{S},] record::Record)::S

Get the sequence of `record` as `S`.

`S` can be either `BioSequences.ReferenceSequence` or `BioSequences.DNASequence`.
If `S` is omitted, the default type is `BioSequences.ReferenceSequence`.
"""
function sequence(record::Record)
    return sequence(BioSequences.ReferenceSequence, record)
end

function sequence(::Type{BioSequences.ReferenceSequence}, record::Record)
    checkfilled(record)
    data = decode_sequence(record.packeddna, record.dnasize, 2, twobit2refseq_table)
    nmask = falses(record.dnasize)
    for i in 1:record.blockcount
        nmask[record.blockstarts[i] .+ (1:record.blocksizes[i])] .= true
    end
    return BioSequences.ReferenceSequence(data, BioSequences.NMask(nmask), 1:record.dnasize)
end

function sequence(::Type{BioSequences.DNASequence}, record::Record)
    checkfilled(record)
    data = decode_sequence(record.packeddna, record.dnasize, 4, twobit2dnaseq_table)
    seq = BioSequences.DNASequence(data, 1:record.dnasize, false)
    for i in 1:record.blockcount
        for j in record.blockstarts[i] .+ (1:record.blocksizes[i])
            seq[j] = BioSequences.DNA_N
        end
    end
    return seq
end

function BioCore.sequence(record::Record)
    return sequence(record)
end

function BioCore.hassequence(record::Record)
    return hassequence(record)
end

"""
    maskedblocks(record::Record)::Vector{UnitRange{Int}}

Get the masked blocks.
"""
function maskedblocks(record::Record)
    checkfilled(record)
    blocks = Vector{UnitRange{Int}}(undef, record.maskedblockcount)
    for i in 1:lastindex(blocks)
        blocks[i] = record.maskedblockstarts[i] .+ (1:record.maskedblocksizes[i])
    end
    return blocks
end

function decode_sequence(packeddna, seqlen, nbits, table)
    @assert nbits ∈ (2, 4)
    data = zeros(UInt64, cld(seqlen, div(64, nbits)))
    stop = BioSequences.BitIndex(seqlen, nbits)
    i = BioSequences.BitIndex(1, nbits)
    j = 1
    while i ≤ stop
        data[BioSequences.index(i)] |= table[Int(packeddna[j])+1] << BioSequences.offset(i)
        i += 4 * nbits
        j += 1
    end
    return data
end

const twobit2refseq_table = let
    # T: 00, C: 01, A: 10, G: 11
    f(x) = x == 0b00 ? UInt64(3) :
           x == 0b01 ? UInt64(1) :
           x == 0b10 ? UInt64(0) :
           x == 0b11 ? UInt64(2) : error()
    tcag = 0b00:0b11
    tbl = UInt64[]
    for x in tcag, y in tcag, z in tcag, w in tcag
        push!(tbl, f(x) | f(y) << 2 | f(z) << 4 | f(w) << 6)
    end
    tbl
end

const twobit2dnaseq_table = let
    # T: 00, C: 01, A: 10, G: 11
    f(x) = x == 0b00 ? convert(UInt64, BioSequences.DNA_T) :
           x == 0b01 ? convert(UInt64, BioSequences.DNA_C) :
           x == 0b10 ? convert(UInt64, BioSequences.DNA_A) :
           x == 0b11 ? convert(UInt64, BioSequences.DNA_G) : error()
    tcag = 0b00:0b11
    tbl = UInt64[]
    for x in tcag, y in tcag, z in tcag, w in tcag
        push!(tbl, f(x) | f(y) << 4 | f(z) << 8 | f(w) << 12)
    end
    tbl
end

function checkfilled(record)
    if !isfilled(record)
        throw(ArgumentError("unfilled 2bit record"))
    end
end

function memcmp(p1::Ptr, p2::Ptr, n::Integer)
    return ccall(:memcmp, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), p1, p2, n)
end

"""
`WriteRecord{S,T}` is a type holding a named sequence of type `S`, along with
a mask.

Mostly to make writing 2bit sequences painless with the depreceation of
SeqRecord.
I can't see this existing much beyond the addition of masked sequences.
Creates a sequence record suitable for writing to 2bit, in a similar way that
FASTA and FASTQ records are created before being written to file. i.e. by calling
a `Record` method on a name, some sequence, and some masks.
"""
mutable struct WriteRecord{S<:BioSequences.Sequence}
    name::String
    seq::S
    masks::Union{Vector{UnitRange{Int}}, Nothing}
end

"""
Record(name::AbstractString, seq::BioSequences.Sequence, masks::Union{Vector{UnitRange{Int}}, Nothing} = nothing)

Prepare a record for writing to a 2bit formatted file.

Needs a `name`, a `sequence`, and (optionally) `masks`: a vector of
ranges that delineate masked regions of sequence.
"""
function Record(name::AbstractString,
                seq::BioSequences.Sequence,
                masks::Union{Vector{UnitRange{Int}}, Nothing} = nothing)
    return WriteRecord(string(name), seq, masks)
end

ext: jl


text:
@testset "Tropic Coordinates" begin
    # example from arXiv:1709.02647, p.13
    pd1 = Dict(0 => diagram(1.0=>3.0, 3.0=>4.0))
    pd2 = Dict(0 => diagram(2.0=>4.0))
    pds = map(ComputationalHomology.flatten, [pd1, pd2])
    d,n,m,mv = ComputationalHomology.tropicstats(pds)
    @test m == 3
    @test tropic(pds[1]) == [2, 3, 4, 6, 7]
    xs = tropic(pds)
    @test xs[:,1] == [2, 3, 4, 6, 7]
    @test xs[:,2] == [2, 2, 4, 4, 4]

    # with half-open interval
    img = reshape([0,0,2,2,2,3,0,0,0,2,0,0,1,3,0,0,0,0,2,1,2,0,0,0,0,0,0,4,1,0,0,0,0,0,1,1,1,1,0,0,0,0,2,0,0,2,1,0,0,0,2,0,0,0,2,0,0,0,1,2,2,3,0,0],8,8)
    flt = filtration(BitmapComplex(maximum(img) .- img))
    pdn = diagram(flt)
    @test tropic(pdn) == [4, 6, 4, 8, 8]

    # img = 255 .- readdlm("misc/du8.csv", ',', UInt8)
    # x = BitmapComplex(img) |> filtration |> diagram |> tropic
end


text:
#
## Array methods
#

function call(::Type{Sidef_Types_Array_Array}, v::Any...)
    Sidef_Types_Array_Array(Any[v...])
end

function getindex(a::Sidef_Types_Array_Array, i::Sidef_Types_Number_Number)

    a = a.value
    if (length(a) <= i.value)
        return NIL
    end

    a[round(Int, i.value+1)]
end

function setindex!(a::Sidef_Types_Array_Array, v::Any, i::Sidef_Types_Number_Number)
    #setindex!(a.value, v, i.value)

    i = i.value
    a = a.value
    len = length(a)

    if (len <= i)
        for _ in len:i+1
            push!(a, NIL)
        end
    end

    a[round(Int, i+1)] = v
end

function is_empty(a::Sidef_Types_Array_Array)
    (length(a.value) == 0) ? TRUE : FALSE
end

function len(a::Sidef_Types_Array_Array)
    Sidef_Types_Number_Number(length(a.value))
end

function ft(a::Sidef_Types_Array_Array, from::Sidef_Types_Number_Number, to::Sidef_Types_Number_Number)
    Sidef_Types_Array_Array((a.value)[from.value+1:to.value+1])
end

function ft(a::Sidef_Types_Array_Array, from::Sidef_Types_Number_Number)
    Sidef_Types_Array_Array((a.value)[from.value+1:length(a.value)])
end

function each(a::Sidef_Types_Array_Array, block::Sidef_Types_Block_Block)
    for item in a.value
        block(item)
    end
end

function map(a::Sidef_Types_Array_Array, block::Sidef_Types_Block_Block)
    arr = Any[]
    for item in a.value
        push!(arr, block(item))
    end
    Sidef_Types_Array_Array(arr)
end

function grep(a::Sidef_Types_Array_Array, block::Sidef_Types_Block_Block)
    arr = Any[]
    for item in a.value
        convert(Bool, block(item)) && push!(arr, item)
    end
    Sidef_Types_Array_Array(arr)
end

function join(a::Sidef_Types_Array_Array)
    str = ""
    for item in a.value
        str *= string(item)
    end
    Sidef_Types_String_String(str)
end

function join(a::Sidef_Types_Array_Array, sep::Sidef_Types_String_String)
    str = ""
    a = a.value
    sep = sep.value
    len = length(a)
    for i in 1:len-1
        str *= string(a[i]) * sep
    end
    if (len > 0)
        str *= string(a[end])
    end
    Sidef_Types_String_String(str)
end

function string(a::Sidef_Types_Array_Array)
    str = "["
    a = a.value
    len = length(a)
    for i in 1:len-1
        str *= string(a[i]) * ", "
    end
    if (len > 0)
        str *= string(a[end])
    end
    str * "]"
end

max_stars_repo_name: trizen/Sidef-Deparse-Julia

max_issues_repo_name: trizen/Sidef-Deparse-Julia

max_forks_repo_name: trizen/Sidef-Deparse-Julia


text:
# ------------------------------------------------------------------------------------------
# # Distributed Arrays
# ------------------------------------------------------------------------------------------

using JuliaRunClient

initializeCluster(12)

@everywhere using DistributedArrays

a=dzeros(12,12)

a=dfill(5.0, 12, 12)

b=drand(1000,1000)

procs(b)

@time sum(b)

convert(Array{Float64}, a)

?DArray

@everywhere function par(I)
    # create our local patch
    # I is a tuple of intervals, each interval is
    # regarded as a 1D array with integer entries
    # size(I[1], 1) gives the number of entries in I[1]
    # size(I[2], 1) gives the number of entries in I[2]
    d=(size(I[1], 1), size(I[2], 1))
    m = fill(myid(), d)
    return m
end

m = DArray(par, (800, 800))

m.indexes

rank(m)

mm = @spawnat 2 rank(localpart(m))
fetch(mm)

# ------------------------------------------------------------------------------------------
# #### Credits
#   * http://www.csd.uwo.ca/~moreno/cs2101a_moreno/Parallel_computing_with_Julia.pdf
#   *
# ------------------------------------------------------------------------------------------

releaseCluster()

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
#!/usr/bin/env julia
using CSVFiles, DataFrames, LinearAlgebra
using Clustering # kmeans
meta_data = DataFrame(load(File(format"TSV", "data/meta_data.tsv.gz")));
categorical!(meta_data, :dataset);
data_mat = DataFrame(load(File(format"TSV", "data/pcs.tsv.gz")));

vars_use = "dataset"

function run_harmony(data_mat, meta_data, vars_use, sigma,
    tau, block_size, max_iter_harmony, max_iter_cluster, epsilon_cluster,
    epsilon_harmony)
    println("Sigma is $sigma")
    N = size(meta_data, 1)
    n_clust = minimum([N ÷ 30, 100])
    theta = 2.
    lamb = 0.1

    sigma = fill(sigma, N)

    labels = meta_data[!, Symbol(vars_use)]
    categories = size(levels(labels),1)

    phi = fill(0., (categories, N))

    for (index, value) in enumerate(levels(labels))
        phi[index,:] = (value .== labels)
    end

    N_b = sum(phi, dims=2)
    Pr_b = N_b / N

    B_vec = [categories]

    println("Bveec is $B_vec")

    theta = fill(theta, size(B_vec,1))

    lamb = fill(lamb, size(B_vec,1))

    lamb_mat = Diagonal(pushfirst!(lamb, 0))

    padding = fill(1., N)

    phi_moe = vcat(padding', phi)

    harmony(data_mat, phi, phi_moe, Pr_b, sigma, theta, n_clust, lamb_mat)

end

function harmony(data_mat, phi, phi_moe, Pr_b, sigma, theta, n_clust, lamb_mat)
    max_iter_cluster = 200
    epsilon_cluster = 1e-5
    epsilon_harmony = 1e-4
    block_size = 0.05

    # normalize data_mat
    data_mat_cos = data_mat ./ (maximum.(eachcol(data_mat)))'

    # calculate norm
    # why not use LinearAlgebra.norm?
    data_mat_norm = mapslices(norm, convert(Matrix, data_mat_cos), dims=1)

    data_mat_cos = data_mat_cos ./ data_mat_norm

    M = convert(Matrix, data_mat_cos)

    d, n = size(data_mat)
    B = size(phi, 1)

    window_size = 3

    # initialize cluster and then harmonize
    println("Initializing cluster with n_clust: $n_clust")
    init_cluster(M, n_clust, sigma, Pr_b, phi)
    harmonize()
end

function init_cluster(M, num_clust, sigma, Pr_b, phi)
    # nb: kmeans works with $d x $n data matrix - hence the transpose.
    km = kmeans(M', num_clust, maxiter=10)

    @assert nclusters(km) == num_clust

    Y = km.centers'

    println("Normalizing...")
    # normalize
    Y_norm = mapslices(norm, convert(Matrix, Y), dims=1)
    Y = Y ./ Y_norm
    println("Begin cluster assignment")
    # assign cluster probabilities
    dist_mat = 2 * (1 .- (Y * M'))

    R = -dist_mat
    sigma = reshape(sigma, (length(sigma), 1));
    R = R ./ sigma';
    max = maximum.(eachcol(R));
    R = R .- reshape(max, (length(max), 1))';
    R = exp.(R)
    summy = sum.(eachcol(R))
    R = R ./ reshape(summy, (length(summy), 1))';
    println("Assignment done!")
    # (3) batch diversity statistics
    E = sum.(eachrow(R)) .* Pr_b'
    O = R * phi'

    println("E shape ", size(E))
    println("O shape ", size(O))

    compute_objective()

end

function harmonize(iter_harmony=10)
    for i in 1:10
        println("This is iteration $i out of $iter_harmony")

        # step 1: clustering
        cluster()

        # step 2: regress out covariates
        Z_cos, Z_corr, W, Phi_Rk = moe_correct_ridge(Z_orig, Z_cos, Z_corr,
        R, W, K, Phi_Rk, Phi_moe, lamb)

        # step 3: check for convergence
        #if check_convergence()
        #    println("Harmony converged after $i iteration(s)")
        #    break
    end
end

function cluster()
end

function moe_correct_ridge()
end

function check_convergence()
end

function compute_objective()
    kmeans_error = 5
end

function safe_entropy(x)
    f(j) = log(j)
    y = x .* map.(f, x)
    if any(isnan.(y))
        y[isnan.(y)] .= 0.
    end
    return y
end

run_harmony(data_mat, meta_data, vars_use, 0.1, 0, 0.05, 10, 200, 1e-5, 1e-4)


text:
module MyExample

# Write your package code here.
print("Hello world")
include("extra_file.jl")

export my_f
end


text:
#Adding Julia Scripts to Jupyter on Windows

Pkg.add("IJulia")

max_stars_repo_name: PacktPublishing/Learning-Jupyter

max_issues_repo_name: PacktPublishing/Learning-Jupyter

max_forks_repo_name: PacktPublishing/Learning-Jupyter


text:
function construct_billiard(polysp::AbstractArray{T, 2}, sp::AbstractArray{T, 2},
    ep::AbstractArray{T, 2}, clr::Vector{Symbol}) where {T<:AbstractFloat}
    polyep = circshift(polysp, (0, -1))
    polyn = [0.0 -1.0; 1.0 0.0]*(polyep - polysp)

    world = []
    for i = 1:size(polysp, 2)
        wall = FiniteWall(polysp[:,i], polyep[:,i], polyn[:,i])
        push!(world, wall)
    end

    n = [0.0 -1.0; 1.0 0.0]*(ep - sp)

    for i = 1:length(clr)
        gate = VirtualGate(sp[:,i], ep[:,i], n[:,i], clr[i])
        push!(world, gate)
        gate = VirtualGate(sp[:,i], ep[:,i], -n[:,i], clr[i])
        push!(world, gate)
    end
    bt = Billiard(world)

    return bt
end

function construct_billiard(polys_sp::AbstractArray{Array{T, 2}, 1}, sp::AbstractArray{T, 2},
    ep::AbstractArray{T, 2}, clr::Vector{Symbol}) where {T<:AbstractFloat}
    polys_ep = []
    for i = 1:length(polys_sp)
        push!(polys_ep, circshift(polys_sp[i], (0, -1)))
    end
    polysp = hcat(polys_sp...)
    polyep = hcat(polys_ep...)
    polyn = [0.0 -1.0; 1.0 0.0]*(polyep - polysp)

    world = []
    for i = 1:size(polysp, 2)
        wall = FiniteWall(polysp[:,i], polyep[:,i], polyn[:,i])
        push!(world, wall)
    end

    n = [0.0 -1.0; 1.0 0.0]*(ep - sp)

    for i = 1:length(clr)
        gate = VirtualGate(sp[:,i], ep[:,i], n[:,i], clr[i])
        push!(world, gate)
        gate = VirtualGate(sp[:,i], ep[:,i], -n[:,i], clr[i])
        push!(world, gate)
    end
    bt = Billiard(world)

    return bt
end

function construct_controlmap(num_gates::T) where {T<:Signed}
    u = Dict()
    for i = 1:num_gates
        push!(u, clr[i]=>0)
    end
    return u
end

max_stars_repo_name: gbarsih/Representations-in-Robotics

max_issues_repo_name: gbarsih/Representations-in-Robotics

max_forks_repo_name: gbarsih/Representations-in-Robotics


text:
function rw(a)
    n = length(a)
    a[1] = rand() - 0.5
    for i = 2:n
        a[i] = a[i-1] + rand() - 0.5
    end
    a
end

function fib(n::Int64)
  if (n <= 2)
     1
  else
     fib(n-1) + fib(n-2)
  end
end

max_stars_repo_name: oeclint/BIDMach

max_issues_repo_name: oeclint/BIDMach

max_forks_repo_name: oeclint/BIDMach


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""QueuingConfiguration holds the configuration parameters for queuing

    IoK8sApiFlowcontrolV1alpha1QueuingConfiguration(;
        handSize=nothing,
        queueLengthLimit=nothing,
        queues=nothing,
    )

    - handSize::Int32 : `handSize` is a small positive number that configures the shuffle sharding of requests into queues.  When enqueuing a request at this priority level the request's flow identifier (a string pair) is hashed and the hash value is used to shuffle the list of queues and deal a hand of the size specified here.  The request is put into one of the shortest queues in that hand. `handSize` must be no larger than `queues`, and should be significantly smaller (so that a few heavy flows do not saturate most of the queues).  See the user-facing documentation for more extensive guidance on setting this field.  This field has a default value of 8.
    - queueLengthLimit::Int32 : `queueLengthLimit` is the maximum number of requests allowed to be waiting in a given queue of this priority level at a time; excess requests are rejected.  This value must be positive.  If not specified, it will be defaulted to 50.
    - queues::Int32 : `queues` is the number of queues for this priority level. The queues exist independently at each apiserver. The value must be positive.  Setting it to 1 effectively precludes shufflesharding and thus makes the distinguisher method of associated flow schemas irrelevant.  This field has a default value of 64.
"""
mutable struct IoK8sApiFlowcontrolV1alpha1QueuingConfiguration <: SwaggerModel
    handSize::Any # spec type: Union{ Nothing, Int32 } # spec name: handSize
    queueLengthLimit::Any # spec type: Union{ Nothing, Int32 } # spec name: queueLengthLimit
    queues::Any # spec type: Union{ Nothing, Int32 } # spec name: queues

    function IoK8sApiFlowcontrolV1alpha1QueuingConfiguration(;handSize=nothing, queueLengthLimit=nothing, queues=nothing)
        o = new()
        validate_property(IoK8sApiFlowcontrolV1alpha1QueuingConfiguration, Symbol("handSize"), handSize)
        setfield!(o, Symbol("handSize"), handSize)
        validate_property(IoK8sApiFlowcontrolV1alpha1QueuingConfiguration, Symbol("queueLengthLimit"), queueLengthLimit)
        setfield!(o, Symbol("queueLengthLimit"), queueLengthLimit)
        validate_property(IoK8sApiFlowcontrolV1alpha1QueuingConfiguration, Symbol("queues"), queues)
        setfield!(o, Symbol("queues"), queues)
        o
    end
end # type IoK8sApiFlowcontrolV1alpha1QueuingConfiguration

const _property_map_IoK8sApiFlowcontrolV1alpha1QueuingConfiguration = Dict{Symbol,Symbol}(Symbol("handSize")=>Symbol("handSize"), Symbol("queueLengthLimit")=>Symbol("queueLengthLimit"), Symbol("queues")=>Symbol("queues"))
const _property_types_IoK8sApiFlowcontrolV1alpha1QueuingConfiguration = Dict{Symbol,String}(Symbol("handSize")=>"Int32", Symbol("queueLengthLimit")=>"Int32", Symbol("queues")=>"Int32")
Base.propertynames(::Type{ IoK8sApiFlowcontrolV1alpha1QueuingConfiguration }) = collect(keys(_property_map_IoK8sApiFlowcontrolV1alpha1QueuingConfiguration))
Swagger.property_type(::Type{ IoK8sApiFlowcontrolV1alpha1QueuingConfiguration }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiFlowcontrolV1alpha1QueuingConfiguration[name]))}
Swagger.field_name(::Type{ IoK8sApiFlowcontrolV1alpha1QueuingConfiguration }, property_name::Symbol) =  _property_map_IoK8sApiFlowcontrolV1alpha1QueuingConfiguration[property_name]

function check_required(o::IoK8sApiFlowcontrolV1alpha1QueuingConfiguration)
    true
end

function validate_property(::Type{ IoK8sApiFlowcontrolV1alpha1QueuingConfiguration }, name::Symbol, val)
end


text:
##################################################################
# Define a convenience class for working with the sparse
# MCMC sample data.
#
# Samples are stored as vectors of "changepoints":
# (idx,val) pairs where the sequence changes to
# value `val` at index `idx`, and stays constant otherwise.
##################################################################

mutable struct ChangepointVec
    changepoints::Vector
    len::Int64
    default_v::Number
end

ChangepointVec(changepoints::Vector, len::Int64) = ChangepointVec(changepoints, len, 0)

import Base: map, reduce, getindex, length, sum

length(cpv::ChangepointVec) = cpv.len

function getindex(cpv::ChangepointVec, r::UnitRange)

    # take care of this corner case
    if length(cpv.changepoints) == 0
        return ChangepointVec([], length(r), cpv.default_v)
    end

    # Index the range into the changepoints
    start_searched = searchsorted(cpv.changepoints, r.start; by=x-> x[1])
    start_idx, start_2 = start_searched.start, start_searched.stop
    stop_searched = searchsorted(cpv.changepoints, r.stop; by=x-> x[1])
    stop_1, stop_idx = stop_searched.start, stop_searched.stop

    new_changepoints = cpv.changepoints[start_idx:stop_idx]

    # Adjust the start of the vector
    if start_idx > start_2 && start_2 != 0
        pushfirst!(new_changepoints, (r.start, cpv.changepoints[start_2][2]))
    end

    # Shift their indices so that we start from 1
    new_changepoints = map(x->(x[1] - r.start + 1, x[2]), new_changepoints)

    return ChangepointVec(new_changepoints, length(r), cpv.default_v)
end

"""
We assume t in [1, length(cpv)].
"""
function getindex(cpv::ChangepointVec, t::Int)

    if length(cpv.changepoints) == 0
        return cpv.default_v
    end

    idx = searchsortedfirst(cpv.changepoints, t; by=x->x[1])

    if idx > length(cpv.changepoints)
        return cpv.changepoints[end][2]
    elseif cpv.changepoints[idx][1] == t
        return cpv.changepoints[idx][2]
    else
        if idx == 1
            return cpv.default_v
        else
            return cpv.changepoints[idx-1][2]
        end
    end
end

function map(f::Function, cpv::ChangepointVec)
    new_changepoints = map(x->(x[1],f(x[2])), cpv.changepoints)
    return ChangepointVec(new_changepoints, cpv.len, cpv.default_v)
end

function sum(cpv::ChangepointVec)

    if length(cpv.changepoints) == 0
        return length(cpv) * cpv.default_v
    end

    s = zero(cpv.default_v)
    if cpv.changepoints[1][1] > 1
        cur_t_v = (1, cpv.default_v)
    else
        cur_t_v = cpv.changepoints[1]
    end

    for cur_changepoint=1:length(cpv.changepoints)-1
        next_t_v = cpv.changepoints[cur_changepoint+1]
        s += cur_t_v[2] * (next_t_v[1] - cur_t_v[1])
        cur_t_v = next_t_v
    end
    s += cur_t_v[2] * (length(cpv) - cur_t_v[1] + 1)

    return s
end

function mult_ind(binary_expr, value)
    if Bool(binary_expr)
        return value
    else
        return 0
    end
end

function count_nonzero(cpv::ChangepointVec)

    if length(cpv.changepoints) == 0
        return mult_ind(cpv.default_v, length(cpv))
    end

    s = zero(cpv.default_v)
    if cpv.changepoints[1][1] > 1
        cur_t_v = (1, cpv.default_v)
    else
        cur_t_v = cpv.changepoints[1]
    end

    for cur_changepoint=1:length(cpv.changepoints)-1
        next_t_v = cpv.changepoints[cur_changepoint+1]
        s += mult_ind(cur_t_v[2], next_t_v[1] - cur_t_v[1])
        cur_t_v = next_t_v
    end

    s += mult_ind(cur_t_v[2], length(cpv) - cur_t_v[1] + 1)

    return s
end

function mean(cpv::ChangepointVec; is_binary::Bool=false)
    if is_binary
        return count_nonzero(cpv) / length(cpv)
    else
        return sum(cpv) / length(cpv)
    end
end

function binop(f::Function, cpva::ChangepointVec, cpvb::ChangepointVec)

    new_default = f(cpva.default_v, cpvb.default_v)
    new_len = length(cpva)
    ncp_a = length(cpva.changepoints)
    ncp_b = length(cpvb.changepoints)
    new_changepoints = Vector{Tuple}()

    # Handle this corner case
    if new_len == 0
        return ChangepointVec([], new_len, new_default)
    end

    cur_t = 1
    cur_cp_a = 0
    cur_cp_b = 0
    cur_v_a = cpva.default_v
    cur_v_b = cpvb.default_v

    while cur_t < new_len

        if cur_cp_a < ncp_a
            #next_t_a, next_v_a = cpva.changepoints[cur_cp_a+1]
            next_cpa = cpva.changepoints[cur_cp_a+1]
            next_t_a = next_cpa[1]
            next_v_a = next_cpa[2]
        else
            next_t_a = new_len
            next_v_a = cur_v_a
        end
        if cur_cp_b < ncp_b
            next_t_b, next_v_b = cpvb.changepoints[cur_cp_b+1]
        else
            next_t_b = new_len
            next_v_b = cur_v_b
        end

        if next_t_a <= next_t_b
            cur_t = next_t_a
            cur_cp_a += 1
            cur_v_a = next_v_a
        end
        if next_t_b <= next_t_a
            cur_t = next_t_b
            cur_cp_b += 1
            cur_v_b = next_v_b
        end

        push!(new_changepoints, (cur_t, f(cur_v_a, cur_v_b)))
    end

    return ChangepointVec(new_changepoints, new_len, new_default)
end

function to_vec(cpv::ChangepointVec)
    vec = Vector{Float64}(undef, cpv.len)
    if cpv.changepoints[1][1] > 1
        vec[1:cpv.changepoints[1][1]] .= cpv.default_v
    end
    for (i, cp) in enumerate(cpv.changepoints[1:end-1])
        vec[cp[1]:cpv.changepoints[i+1][1]] .= cp[2]
    end
    last = cpv.changepoints[end]
    vec[last[1]:end] .= last[2]
    return vec
end

max_stars_repo_name: gitter-lab/ssps

max_issues_repo_name: gitter-lab/ssps

max_forks_repo_name: gitter-lab/ssps


text:
using IbnrClaims
using Test

@testset "IbnrClaims.jl" begin
    si = [1996 463 157 158 232 115; 2762 289 456 273 407 0; 2898 163 432 283 0 0; 2930 315 217 0 0 0; 2474 395 0 0 0 0; 2684 0 0 0 0 0]
    sa = [1996 2459 2616 2774 3006 3121; 2762 3051 3507 3780 4187 0; 2898 3061 3493 3776 0 0; 2930 3245 3462 0 0 0; 2474 2869 0 0 0 0; 2684 0 0 0 0 0]
    pr = [3095,5640,5416,5069,4657,5352]
    lr = [0.85,0.85,0.85,0.85,0.85,0.85]

    ti = Incremental(si)
    ta = Accumulated(sa)

    @test ibnrcl(ta) ≈ 3885.51761789107
    @test ibnrbf(ta, pr, lr) ≈ 4128.66896034168

    @test ti.claims == tdeacum(ta).claims
    @test ta.claims == taccum(ti).claims
    @test ti.claims == tdeacum(taccum(ti)).claims
    @test ta.claims == taccum(tdeacum(ta)).claims
    @test incurred(ta) == [3121, 4187, 3776, 3462, 2869, 2684]

    @test tfactors(ti) ≈ [1.1244257274, 1.1068043331, 1.0742512479, 1.09749771132, 1.03825681969, 1.0]
    @test tfactors(ta) ≈ [1.1244257274, 1.1068043331, 1.0742512479, 1.09749771132, 1.03825681969, 1.0]

end


text:
"""
    st_dim(ga::AbstractGeoArray; mid::Vector{Int} = [1, 1])
    st_dim(b::bbox, cellsize::T) where {T <: Real}
    st_dim(ga::AbstractGeoArray, dim::Symbol; mid::Vector{Int} = [1, 1])

Get dimensions of x and y
"""
function st_dim(ga::AbstractGeoArray; mid::Vector{Int} = [1, 1])
    if length(mid) == 1; mid = [mid, mid]; end

    cellsize_x = ga.f.linear[1]
    cellsize_y = abs(ga.f.linear[4])
    cellsize_y2 = ga.f.linear[4]

    delta = [cellsize_x, cellsize_y]/2 .* mid

    rbbox = st_bbox(ga)
    x = rbbox.xmin + delta[1]:cellsize_x:rbbox.xmax
    y = rbbox.ymin + delta[2]:cellsize_y:rbbox.ymax
    if cellsize_y2 < 0; y = reverse(y); end

    x, y
end

function st_dim(b::bbox, cellsize::T) where {T <: Real}
    lon = b.xmin + cellsize/2 : cellsize : b.xmax
    lat = reverse(b.ymin + cellsize/2 : cellsize : b.ymax)
    lon, lat # return
end

st_dim_x(ga::AbstractGeoArray; mid::Vector{Int} = [1, 1]) = st_dim(ga; mid)[1]
st_dim_y(ga::AbstractGeoArray; mid::Vector{Int} = [1, 1]) = st_dim(ga; mid)[2]

function st_dim(ga::AbstractGeoArray, dim::Symbol; mid::Vector{Int} = [1, 1])
    if dim==:x
        ci = st_dim_x(ga; mid = mid)
    elseif dim==:y
        ci = st_dim_y(ga; mid = mid)
    else
        error("Use :x or :y as second argument")
    end
    return ci
end

export meshgrid, st_dim


text:
using EffectiveWaves
using LinearAlgebra
using Plots
pyplot()

medium = Medium(ρ=1.0, c=1.0)
species = [Specie(ρ=0.1,r=1.0, c=0.01, volfrac=0.3)]

ωs = [0.5]
ω = ωs[1]
k = ω/medium.c
tol = 1e-8

kp = 1.0+0.4im
kps = collect(0.2:0.2:1.0) .+ collect(0.3:0.2:1.1) .* im
hankel_order = ho = 2

wavesystemPlane = wavematrix3DPlane(ωs[1], medium, species; tol = 1e-5, hankel_order=hankel_order)
detP(kp) = det(wavesystemPlane(kp))

# dispersion = dispersion_function(ω, medium, species; tol = low_tol, dim=dim, symmetry = :plane, hankel_order = ho)

# dim = 3
# k_effs = wavenumbers_path(ω, medium, species;
#     num_wavenumbers=7, dim=dim, symmetry = :plane,
#     hankel_order = ho, tol=tol)

# detsplanes = det.(wavesystemPlane.(kps))
# MPs = wavesystemPlane.(kps)

wavesystem3D = wavematrix3D(ωs[1], medium, species; tol = 1e-5, hankel_order=hankel_order);
detR(kp) = sqrt(det(wavesystem3D(kp)))

wavesystem3D = wavematrix3D(ωs[1], medium, species; tol = 1e-5, hankel_order=1);
detR(kp) = sqrt(det(wavesystem3D(kp)))

# wavesystem3D_2 = wavematrix3D_allocate(ωs[1], medium, species; tol = 1e-4, hankel_order=hankel_order);
# detR2(kp) = det(wavesystem3D_2(kp))

# @time detR.(kps);
#   1.996288 seconds (469.46 k allocations: 52.257 MiB, 1.06% gc time)
#
# @time detR.(kps);
#   2.103966 seconds (327.72 k allocations: 45.012 MiB, 0.39% gc time)
#
# @time detR.(kps);
# 0.084732 seconds (32.86 k allocations: 3.510 MiB)

# Effective medium
    dim = 3; T= Float64
    β = medium.ρ*medium.c^2 # medium bulk modulus
    φ = sum(volume_fraction.(species; dim = dim))
    β_eff = one(T)/((1-φ)/β + sum(volume_fraction(s; dim = dim)/(s.ρ*s.c^2) for s in species))
    ρ_frac = sum(volume_fraction(s; dim = dim)*(medium.ρ - s.ρ)/(medium.ρ + 2*s.ρ) for s in species)

    ρ_eff = medium.ρ*(1 - ρ_frac)/(1 + T(dim - 1) * ρ_frac)
    eff_medium = Medium(ρ=ρ_eff, c=sqrt(β_eff/ρ_eff))
    k_low1 = ω/eff_medium.c

    ρ_eff = medium.ρ*(1 - ρ_frac)/(1 + (ho/(2ho+3))*ρ_frac)
    eff_medium = Medium(ρ=ρ_eff, c=sqrt(β/ρ_eff))
    k_low2 = ω/eff_medium.c

    ρ_eff = medium.ρ*(1 - ρ_frac)/(1 + ((ho-1)/(2ho+1))*ρ_frac)
    eff_medium = Medium(ρ=ρ_eff, c=sqrt(β_eff/ρ_eff))
    k_low3 = ω/eff_medium.c

kp0s = wavenumbers(ω, medium, species; dim = dim, symmetry = :plane, hankel_order = ho,
    num_wavenumbers=20, tol = tol, mesh_points = 14, mesh_size = 3.0)

kp_vecs = reduce_kvecs([[real(kp),imag(kp)] for kp in kp0s], 1e-5)
kps = [kp[1]+im*kp[2] for kp in kp_vecs]
kps = kps[1:30]

# inds = findall(abs.(detP.(kps)) .< 1e-6)
# kps = kp0s[inds]

# For low frequency
    # k_lows = [k_low1,k_low2,k_low3]
    # k1 = minimum(real.(k_lows[[1,3]]))
    # k2 = maximum(real.(k_lows))
    # x = LinRange(0.9*k1,1.1*k2,110)
    # y = k2 .* LinRange(-0.01,0.2,90)

# Given some kps
    # dk1 = real(kps[1]) / 4
    # dk2 = imag(kps[1]) / 4
    k1 = maximum(abs.(real.(kps)))
    k2 = imag(kps[end])
    x = LinRange(-k1,k1,6600)
    y = LinRange(-0.01,k2,160)

X = repeat(x',length(y),1)
Y = repeat(y,1,length(x))

ZPs = map( (x,y) -> abs(detP(x+y*im)),X,Y)
# ZRs = map( (x,y) -> abs(detR(x+y*im)),X,Y)

minR = 1e-18
minP = 0.8

# ZRs2 = [ (abs(ZRs[ind]) < minR ? ZRs[ind] : NaN ) for ind in LinearIndices(ZRs)];
ZPs2 = [ (abs(ZPs[ind]) < minP ? ZPs[ind] : NaN ) for ind in LinearIndices(ZPs)];

# h1 = heatmap(x,y,ZPs, xlab = "Re k*", ylab = "Im k*", title="Planewave det(I + G), ka~1", clims = (0.,1.0))

# h2 = heatmap(x,y,ZRs2, xlab = "Re k*", ylab = "Im k*", title="General det(I + G), ka~$(ω*species[1].r)", clims = (0.,minR))
h1 = heatmap(x,y,ZPs2, xlab = "Re k*", ylab = "Im k*", title="Planewave det(I + G), ka~$(ω*species[1].r)", clims = (0.,minP))
scatter!(kps)

# Z = map( (x,y) -> abs(((x^2-1)^2*(y*im - 1)*(y*im - 2))),X,Y)
# Z = [ (abs(Z[ind]) < 1.0 ? Z[ind] : NaN ) for ind in LinearIndices(Z)]
# heatmap(x,y,Z, xlab = "Re k*", ylab = "Im k*", clims = (0.,1.0))

abs.(detsplanes)
abs.(dets3D)

h1 = heatmap(x,y,ZPs, xlab = "Re k*", ylab = "Im k*", title="Planewave det(I + G), ka~$(ω*species[1].r)", clims = (0.,1.0))
scatter!(kps)


text:
# Perform actomyosin network simulations
# Alex Tam, 12/10/2020

# Load packages
using Revise # Modify code without restarting Julia
using Parameters # Tools for storing parameters in data structures
using LinearAlgebra # Matrix and vector operations
using StatsBase # Sampling without replacement
using Plots # Plotting library
using Plots.PlotMeasures # Enable re-sizing plot margins
using LaTeXStrings # Display LaTeX output
using ForwardDiff # Automatic differentiation tools
using Optim # Optimisation routines
using LineSearches # Line searches for optimisation
using Printf # C-style printing macros
using Distributions # Probability distributions
using DelimitedFiles # Read and write delimited data

# Include code from external files
include("model_parameters.jl")
include("actomyosin_network.jl")
include("State.jl")
include("initial_condition.jl")
include("Actin_Filament.jl")
include("periodic.jl")
include("Cross_Link.jl")
include("Myosin_Motor.jl")
include("intersection_search.jl")
include("network_statistics.jl")
include("network_force.jl")
include("draw_network.jl")
include("turnover.jl")
include("optimise_network.jl")
include("energy.jl")

function run_simulations()
    par_ref = 0.04;
    nTrials = 1;
    for i = 1:length(par_ref)
        par = par_ref[i]; # Select parameter
        bulk_stress = Vector{Float64}(); # Pre-allocate time-averaged bulk stress
        curvature = Vector{Float64}(); # Pre-allocate curvature
        index = Vector{Float64}(); # Pre-allocate integrated tension
        filament_speed = Vector{Float64}(); # Pre-allocate actin filament node speed
        motor_speed = Vector{Float64}(); # Pre-allocate motor head speed
        angle_roc = Vector{Float64}(); # Pre-allocate rate of change of angles
        motor_pos = Vector{Float64}(); # Pre-allocate time-averaged motor position
        motor_angle = Vector{Float64}(); # Pre-allocate time-averaged motor angle
        for j = 1:nTrials
            # Specify parameters
            parN = Numerical_Parameters(); # Initialise struct of numerical parameters
            parA = Actin_Properties(k_off = par); # Initialise struct of actin filament properties
            parM = Myosin_Properties(); # Initialise struct of myosin motor properties

            # Run simulations
            @time state, af, mm, xl, Bulk_Stress, Curvature, Index, Filament_Speed, Motor_Speed, Angle_ROC, Motor_Pos, Motor_Angle, Bins, Hist = actomyosin_network(parN, parA, parM, i, j);
            push!(bulk_stress, Bulk_Stress); # Store time-averaged net stress
            push!(curvature, Curvature); # Store time-averaged integrated curvature
            push!(index, Index); # Store time-averaged two-filament index
            push!(filament_speed, Filament_Speed); # Store time-averaged filament node speed
            push!(motor_speed, Motor_Speed); # Store time-averaged motor head speed
            push!(angle_roc, Angle_ROC); # Store time-averaged rate of change of angle
            push!(motor_pos, Motor_Pos); # Store time-averaged motor position
            push!(motor_angle, Motor_Angle); # Store time-averaged motor angle
            writedlm("bulk_stress-$par.txt", bulk_stress); # Write time-averaged bulk stress per trial to file
            writedlm("curvature-$par.txt", curvature); # Write time-averaged curvature per trial to file
            writedlm("index-$par.txt", index); # Write time-averaged two-filament index per trial to file
            writedlm("filament_speed-$par.txt", filament_speed); # Write time-averaged actin filament node speed per trial to file
            writedlm("motor_speed-$par.txt", motor_speed); # Write time-averaged motor head speed per trial to file
            writedlm("angle_roc-$par.txt", angle_roc); # Write time-averaged angle rate of change to file
            writedlm("motor_pos-$par.txt", motor_pos); # Write time-averaged motor position to file
            writedlm("motor_angle-$par.txt", motor_angle); # Write time-averaged motor angle to file
            writedlm("pd_end_bins-$par-trial-$j.txt", Bins); # Write paried-distances bins to file
            writedlm("pd_end_counts-$par-trial-$j.txt", Hist); # Write paired-distances counts to file
        end
    end
end

run_simulations();


text:
function EarthStandardAtmosphereRefraction()
    return EarthStandardAtmosphereRefraction(())
end

function EarthStandardAtmosphereRefraction(arg0::jdouble, arg1::jdouble)
    return EarthStandardAtmosphereRefraction((jdouble, jdouble), arg0, arg1)
end

function get_pressure(obj::EarthStandardAtmosphereRefraction)
    return jcall(obj, "getPressure", jdouble, ())
end

function get_refraction(obj::EarthStandardAtmosphereRefraction, arg0::jdouble)
    return jcall(obj, "getRefraction", jdouble, (jdouble,), arg0)
end

function get_temperature(obj::EarthStandardAtmosphereRefraction)
    return jcall(obj, "getTemperature", jdouble, ())
end

function set_pressure(obj::EarthStandardAtmosphereRefraction, arg0::jdouble)
    return jcall(obj, "setPressure", void, (jdouble,), arg0)
end

function set_temperature(obj::EarthStandardAtmosphereRefraction, arg0::jdouble)
    return jcall(obj, "setTemperature", void, (jdouble,), arg0)
end


text:
function load_result(res_dir)
    tree = Core.ResDirTree(res_dir);
    return Core.load_data(tree);
end

max_stars_repo_name: josePereiro/MyResultsManager

max_issues_repo_name: josePereiro/MyResultsManager

max_forks_repo_name: josePereiro/MyResultsManager


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Libepoxy_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Libepoxy")
JLLWrappers.@generate_main_file("Libepoxy", UUID("42c93a91-0102-5b3f-8f9d-e41de60ac950"))
end  # module Libepoxy_jll


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

using Test, Profile, Serialization

Profile.clear()
Profile.init()

let iobuf = IOBuffer()
    for fmt in (:tree, :flat)
        Test.@test_logs (:warn, r"^There were no samples collected\.") Profile.print(iobuf, format=fmt, C=true)
        Test.@test_logs (:warn, r"^There were no samples collected\.") Profile.print(iobuf, [0x0000000000000001], Dict(0x0000000000000001 => [Base.StackTraces.UNKNOWN]), format=fmt, C=false)
    end
end

@noinline function busywait(t, n_tries)
    iter = 0
    init_data = Profile.len_data()
    while iter < n_tries && Profile.len_data() == init_data
        iter += 1
        tend = time() + t
        while time() < tend end
    end
end

@profile busywait(1, 20)

let r = Profile.retrieve()
    mktemp() do path, io
        serialize(io, r)
        close(io)
        open(path) do io
            @test isa(deserialize(io), Tuple{Vector{UInt},Dict{UInt64,Vector{Base.StackTraces.StackFrame}}})
        end
    end
end

let iobuf = IOBuffer()
    Profile.print(iobuf, format=:tree, C=true)
    str = String(take!(iobuf))
    @test !isempty(str)
    truncate(iobuf, 0)
    Profile.print(iobuf, format=:tree, maxdepth=2)
    str = String(take!(iobuf))
    @test !isempty(str)
    truncate(iobuf, 0)
    Profile.print(iobuf, format=:flat, C=true)
    str = String(take!(iobuf))
    @test !isempty(str)
    truncate(iobuf, 0)
    Profile.print(iobuf)
    @test !isempty(String(take!(iobuf)))
    truncate(iobuf, 0)
    Profile.print(iobuf, format=:flat, sortedby=:count)
    @test !isempty(String(take!(iobuf)))
    Profile.clear()
    @test isempty(Profile.fetch())
    @test Profile.callers("\\") !== nothing
    @test Profile.callers(\) !== nothing
    # linerange with no filename provided should fail
    @test_throws ArgumentError Profile.callers(\; linerange=10:50)
end

# issue #13229
module I13229
using Test, Profile
global z = 0
@timed @profile for i = 1:5
    function f(x)
        return x + i
    end
    global z = f(i)
end
@test z == 10
end

@testset "setting sample count and delay in init" begin
    n_, delay_ = Profile.init()
    @test n_ == 1_000_000
    def_delay = Sys.iswindows() ? 0.01 : 0.001
    @test delay_ == def_delay
    Profile.init(n=1_000_001, delay=0.0005)
    n_, delay_ = Profile.init()
    @test n_ == 1_000_001
    @test delay_ == 0.0005
    Profile.init(n=1_000_000, delay=def_delay)
end

max_stars_repo_name: smallnamespace/julia

max_issues_repo_name: smallnamespace/julia

max_forks_repo_name: smallnamespace/julia


text:
using IntervalArithmetic, IntervalRootFinding

struct Quad{T}
    a::Interval{T}
    b::Interval{T}
    c::Interval{T}
    x::Vector{Interval{T}}
end

@testset "Quadratic Interval Equations" begin
    rts = Quad{Float64}[]
    push!(rts, Quad(interval(1, 5), interval(2, 4), interval(0, 2), [interval(-4, -2), interval(-0, -0)]))
    push!(rts, Quad(interval(-1, 1), interval(-2, 2), interval(-1, 1), [interval(-∞, -1), interval(-1, -1), interval(-1, ∞)]))
    push!(rts, Quad(interval(1, 2), interval(1, 3), interval(2, 6), [interval(-2, -1)]))
    push!(rts, Quad(interval(1, 1), interval(2, 2), interval(1, 1), [interval(-1, -1), interval(-1, -1)]))

    for i in 1:length(rts)
            @test quadratic_roots(rts[i].a, rts[i].b, rts[i].c) == rts[i].x
    end
end


text:
using Base.Test
using Polynomials
using TrajOpt

@testset "HermiteBasis" begin

    """ Check basis function values and derivatives """
    function check_endpoints(basis)
        all_pass = true
        for i = 1:length(basis.left)

            # Check i-1st derivative at left endpoint
            all_pass &= polyder( basis.left[i],  i-1 )(0.0) == 1.0
            all_pass &= polyder( basis.right[i], i-1 )(0.0) == 0.0

            # Check i-1st derivative at right endpoint
            all_pass &= polyder( basis.left[i],  i-1 )(1.0) == 0.0
            all_pass &= polyder( basis.right[i], i-1 )(1.0) == 1.0

        end
        return all_pass
    end

    H1 = HermiteBasis(degree = 1)
    H3 = HermiteBasis(degree = 3)
    H5 = HermiteBasis(degree = 5)
    @test_throws ArgumentError HermiteBasis(degree = 2)

    # Verify properties of the basis functions
    @test check_endpoints(H1)
    @test check_endpoints(H3)
    @test check_endpoints(H5)

    # Verify basis function are stored in a vector
    @test size(H3.left)  == (2,)
    @test size(H5.left)  == (3,)
    @test size(H5.right) == (3,)

end

@testset "QuadratureRule" begin

    """ Check quadrature rule accuracy """
    function check_quadrature(qr::QuadratureRule)

        npt = length(qr.locations)
        test_function = Poly(ones(2*npt))

        # Compute approximate integral between [-1.0, 1.0] via quadrature
        samples = polyval(test_function, qr.locations)
        approx = dot(qr.weights, samples)

        # Compute integral between [-1.0, 1.0] via symbolic integration
        integral = polyint(test_function)
        exact = integral(1.0) - integral(-1.0)

        return exact ≈ approx

    end

    Q1 = QuadratureRule(num_points = 1)
    Q2 = QuadratureRule(num_points = 2)
    Q3 = QuadratureRule(num_points = 3)
    @test_throws ArgumentError QuadratureRule(num_points = 4)

    # Verify properties of the basis functions
    @test check_quadrature(Q1)
    @test check_quadrature(Q2)
    @test check_quadrature(Q3)

    # Verify basis function are stored in a vector
    @test size(Q2.locations)  == (2,)
    @test size(Q3.locations)  == (3,)
    @test size(Q3.weights)    == (3,)

end

@testset "Test Cases" begin

    """
    # Test Case 1: Constant Acceleration in 1D w/ Full Initial Conditions

    This is a simple test to verify that the stiffness matrix, interpolation
    matrices and boundary conditions are implemented correctly. In the case of a
    constant linear acceleration, the load vector is independent of the solution and
    the trajectory may be directly computed by inverting the K matrix onto the load
    vector (no linearization or iteration required). Since the exact solution is a
    quadratic polynomial (x = x0 + v0*t + 0.5*a*t^2), the discretized solution
    should match the exact solution for a HermiteBasis with degree >= 3.

    """
    function test_case_1()

        const x0 = 3.0  # Initial displacement [m]
        const v0 = 2.5  # Initial velocity     [m/s]
        const a  = 1.0  # Acceleration         [m/s^2]
        const Δt = 4.0  # Integration Period   [s]

        # Define time discretization and solution method
        time_points = linspace(0.0, Δt, 9)
        qr = QuadratureRule(num_points = 2)
        hb = HermiteBasis(degree = 3)

        # Setup stiffness matrix
        K = assemble_stiffness_matrix(hb, time_points)
        apply_bc!(K, BCTypes.dirichlet, 1)
        apply_bc!(K, BCTypes.neumann,   1)

        # Setup load vector
        (HL,HR) = compute_interpolation_matrices(hb, qr, time_points)
        #F = assemble_load_vector(HL, HR, qr.weights, a)
        #apply_bc!(F, BCTypes.dirichlet, 1, x0)
        #apply_bc!(F, BCTypes.neumann,   1, v0)

        # Solve
        #x = K \ F

        return true #x[end][pos] == x0 + v0*Δt + 0.5*a*Δt^2

    end
    @test test_case_1()

end


text:
# Duffing_orbit_diagram_vs_period_zooom_v1.jl
# ===========================================

using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, OrdinaryDiffEq
using StaticArrays
using Random
using Statistics

function Duffing_ODE(du,u,p,t)
    d, a, omega = p
    du[1] =  u[2]
    du[2] =  - u[1] - u[1]*u[1]*u[1] - d*u[2] + a*sin(omega*t)
end

function Poincare_orbit(d,a,omega,y_init,ntrans,nattr,pcol)

    p = [d a omega]
    period = 2*pi / omega

 # transient to attractor
 # ----------------------
    tspan = (0., ntrans*period)
    prob = ODEProblem(Duffing_ODE, y_init, tspan, p)
    sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

  # attractor  with root of Poincare map
  # ---------
    tspan = (0., nattr*period)
    y_init = sol.u[end]
    prob = ODEProblem(Duffing_ODE, y_init, tspan, p)
    sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8, saveat = period/2)

    u1vec = [u[1] for u in sol.u]
    npts = length(u1vec)
    for j = 2:2:npts
        u1vec[j] = -u1vec[j]
    end

  # plot
  # ----
    plt.scatter(period*ones(npts),u1vec,s=1, c = pcol)

    return y_init

end # of function

# parameters and initial values
# -----------------------------
d = 0.1
avec = [19.7 19.75]
n_T = 1500

ntrans = 800   # no. of transient periods
nattr = 64     # no. of points on the attractor
y_init = [1 0]

# open plot
# ---------
fig = figure()

# subplot loop
# ------------
@time for iplt = 1:2

    if iplt  == 1
        a = avec[1]
        T_min = 4.1 # 1
        T_max = 5.8 # 32
        T_vec = LinRange(T_min, T_max, n_T)
    else
        a = avec[2]
        T_min = 5.06 # 1
        T_max = 5.42 # 32
        T_vec = LinRange(T_max, T_min, n_T)
    end

    subplot(1,2,iplt)

    xlabel(L"T = 2\pi / \omega")
    ylabel(L"u_1")
    pcol = COLORS[1]

#  control parameter loop (period)
#  ----------------------
   for i_T = 1:length(T_vec)
       period = T_vec[i_T]
       omega = 2*pi / period
       y_init =  Poincare_orbit(d,a,omega,y_init,ntrans,nattr,pcol)
   end # of period loop

#   tight_layout()
#   ax = gca()
Tbmin = 5.1
Tbmax = 5.3
u1bmin = -0.94
u1bmax = -0.7

#  plot zoom box
if iplt == 1
#    xbox = [Tbmin, Tbmax, Tbmax , Tbmin, Tbmin]
#    ybox = [u1bmin, u1bmin, u1bmax, u1bmax, u1bmin]
#    plt.plot(xbox,ybox,linewidth=1.0,color="k")
   ax = gca()
   ax.set_xlim([T_min,T_max])
   ax.set_ylim([-1.7,1.4])
end
if iplt == 2
    ax = gca()
    ax.set_xlim([Tbmin,Tbmax])
    ax.set_ylim([u1bmin,u1bmax])
end

end # of subplot loop

# %%

add_identifiers!(fig)
fig.tight_layout(pad=0.3)

wsave(plotsdir("9", "duffing_orbitdiagram_zoom"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
# Copyright 2017 Raytheon BBN Technologies
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the license at
#
#		http://www.apache.org/licenses/LICENSE-2.0
using Test

function in_range(x, minx, maxx)
  return (x >= minx) && (x <= maxx)
end

function check_fit_result(true_values::Dict{String, Float64}, result::Qlab.FitResult)
    @test Set(keys(true_values)) == Set(keys(result.fit_params))
    @test Set(keys(true_values)) == Set(keys(result.errors))

    for k in keys(true_values)
        println("$(k): $(result.fit_params[k]) => $(true_values[k]) ± $(result.errors[k])")
        @test in_range(result.fit_params[k], true_values[k]-10*result.errors[k], true_values[k]+10*result.errors[k])
    end
end

#y = a exp(-t/T₁) + b
function test_T1_fit()
    #generate some fake data
    vals = Dict("a" => 0.87, "T" => 20., "b" => 0.1)
    t = collect(range(0., stop=100., length=201))
    y = vals["a"] * exp.(-t/vals["T"]) .+ vals["b"] .+ 0.01 * randn((201,))
    result = Qlab.fit_t1(t, y)
    check_fit_result(vals, result)
end

#y = ax + b.
function test_line_fit()
    vals = Dict("a" => -0.4, "b" => 2.1)
    x = range(-4, stop=4, length=201)
    y = vals["a"] * x .+ vals["b"] .+ 0.05 * randn((201,))
    result = Qlab.fit_line(x,y)
    check_fit_result(vals, result)
end

#a*exp(-t ./ T).*cos(2πf .* t + ϕ) + b
function test_ramsey_fit()
    p = Dict("a"=>0.8, "T"=>12.3, "f"=>0.55, "ϕ"=>π/12, "b"=>-0.15)
    t = range(0, stop=25, length=201)
    y = p["a"]*exp.(-t ./ p["T"]).*cos.(2π*p["f"] .* t .+ p["ϕ"]) .+ p["b"] + 0.01 * randn((201,))
    result = Qlab.fit_ramsey(t,y)
    check_fit_result(p, result)
end

test_T1_fit()
test_line_fit()
test_ramsey_fit()


text:
using MLJLinearModels
using BenchmarkTools, Random, LinearAlgebra
DO_COMPARISONS = false; include("../testutils.jl")

n, p = 50_000, 500
((X, y, θ), (X1, y1, θ1)) = generate_binary(n, p; seed=525)

λ = 5.0
lr = LogisticRegression(λ)

# XXX need to check whether this is doing the right thing or not... might
# be worth comparing with your own implementation using a simple backtracking
# line search or your own implementation of Hanger-Zhang

# newtoncg = NewtonCG()
# θ_ncg = fit(lr, X, y1, solver=newtoncg)

lbfgs = LBFGS()

# by far the fastest
@btime θ_lbfgs = fit($lr, $X, $y1, solver=$lbfgs)

# <old>       1.421 s (2866 allocations: 142.33 MiB)
# with cache  1.357 s (1912 allocations:  30.15 MiB)
# May 12, 2020:: 1.762 s (10696 allocations: 112.38 MiB)
#             :: 1.563 s (10193 allocations: 29.68 MiB)


text:
@testset "Variational Student-T Processes" begin
    N = 50

    x = rand(N)
    X = collect(reshape(x, :, 1))
    y = rand(N)
    ν = 3.0
    k = SqExponentialKernel()
    l = LaplaceLikelihood()
    vi = AnalyticVI()
    m = VStP(x, y, k, l, vi, ν)
    @test_throws ErrorException VStP(x, y, k, l, vi, 0.5)
    @test_throws ErrorException VStP(x, y, k, l, QuadratureVI(), 0.5)
    @test repr(m) ==
          "Variational Student-T Process with a Laplace likelihood (β=1.0) infered by Analytic Variational Inference "
    train!(m, 20)
end


text:
# Testing algos generating an invalid input.
function funcfunc(x::Array{Float64,1})
    # The first coordinate convergence to 4.0 by 1 unit per iterate.
    output = Array{Float64,1}(x)
    if abs(x[1] - 4.0) <= 1.0
        output[1] = 4.0
    elseif x[1] > 4.0
        output[1] = x[1] - 1.0
    else
        output[1] = x[1] + 1.0
    end
    # The second does aitken convergence to 2.3
    output[2] = x[2] + (2.3-x[2])/2.0
    return output
end
Inputs = [19.0,10.0]
fp = fixed_point(funcfunc, Inputs; Algorithm = Aitken)
fp.FailedEvaluation_.Error_ == InputInfsDetected
fp.FailedEvaluation_.Input_[1] == -Inf
!isinf(fp.FailedEvaluation_.Input_[2])
# Now fixing with replace element
fp = fixed_point(funcfunc, Inputs; Algorithm = Aitken, ReplaceInvalids = ReplaceElements)
fp.TerminationCondition_ == ReachedConvergenceThreshold
# Now fixing with replace element
fp2 = fixed_point(funcfunc, Inputs; Algorithm = Aitken, ReplaceInvalids = ReplaceVector)
fp2.TerminationCondition_ == ReachedConvergenceThreshold
fp.Iterations_ != fp2.Iterations_


text:
using BinaryBuilder

name = "CGNS"
version = v"4.3.0"

sources = [
    ArchiveSource("https://github.com/CGNS/CGNS/archive/refs/tags/v$(version).tar.gz",
                  "7709eb7d99731dea0dd1eff183f109eaef8d9556624e3fbc34dc5177afc0a032"),
]

script = raw"""
cd ${WORKSPACE}/srcdir/CGNS*
mkdir build && cd build
H5LIB=""
if [[ "${target}" == *-mingw* ]]; then
    H5LIB="-DHDF5_hdf5_LIBRARY_RELEASE=$(ls ${WORKSPACE}/destdir/bin/libhdf5-*.${dlext})"
fi
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
      -DCMAKE_BUILD_TYPE=Release \
      -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}" \
      ${H5LIB} ..
make -j${nproc}
make install
"""

platforms = [
    Platform("x86_64", "linux"),
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("x86_64", "macos"),
    Platform("x86_64", "windows"),
    Platform("i686", "windows"),
    Platform("aarch64", "macos"),
]

products = [
    LibraryProduct("libcgns", :libcgns),
    ExecutableProduct("cgnscheck", :cgnscheck),
    ExecutableProduct("cgnscompress", :cgnscompress),
    ExecutableProduct("cgnsconvert", :cgnsconvert),
    ExecutableProduct("cgnsdiff", :cgnsdiff),
    ExecutableProduct("cgnslist", :cgnslist),
    ExecutableProduct("cgnsnames", :cgnsnames),
]

dependencies = [
    Dependency("HDF5_jll"),
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
@testset "broadcast" begin
    i = HashIndices([1,2,3,4,5])

    @test issetequal(pairs((i .+ 1)::HashDictionary), [1=>2, 2=>3, 3=>4, 4=>5, 5=>6])
    @test issetequal(pairs(Base.Broadcast.broadcasted(+, i, 1)::BroadcastedDictionary), [1=>2, 2=>3, 3=>4, 4=>5, 5=>6])

    @test issetequal(pairs((i .+ i)::HashDictionary), [1=>2, 2=>4, 3=>6, 4=>8, 5=>10])
    @test issetequal(pairs(Base.Broadcast.broadcasted(+, i, i)::BroadcastedDictionary), [1=>2, 2=>4, 3=>6, 4=>8, 5=>10])

    @test issetequal(pairs((i .+ copy(i))::HashDictionary), [1=>2, 2=>4, 3=>6, 4=>8, 5=>10])
    @test issetequal(pairs(Base.Broadcast.broadcasted(+, i, copy(i))::BroadcastedDictionary), [1=>2, 2=>4, 3=>6, 4=>8, 5=>10])

    d = i .+ 1

    @test issetequal(pairs((d .+ 1)::HashDictionary), [1=>3, 2=>4, 3=>5, 4=>6, 5=>7])
    @test issetequal(pairs(Base.Broadcast.broadcasted(+, d, 1)::BroadcastedDictionary), [1=>3, 2=>4, 3=>5, 4=>6, 5=>7])

    @test issetequal(pairs((d .+ d)::HashDictionary), [1=>4, 2=>6, 3=>8, 4=>10, 5=>12])
    @test issetequal(pairs(Base.Broadcast.broadcasted(+, d, d)::BroadcastedDictionary), [1=>4, 2=>6, 3=>8, 4=>10, 5=>12])

    @test issetequal(pairs((d .+ copy(d))::HashDictionary), [1=>4, 2=>6, 3=>8, 4=>10, 5=>12])
    @test issetequal(pairs(Base.Broadcast.broadcasted(+, d, copy(d))::BroadcastedDictionary), [1=>4, 2=>6, 3=>8, 4=>10, 5=>12])

    d2 = similar(d)
    d2 .= d .+ d
    @test issetequal(pairs(d2), [1=>4, 2=>6, 3=>8, 4=>10, 5=>12])
end


text:
import Base.show
import Base.convert
import Base.getindex
#using Compat

abstract Rule

# by default no children
get_children(rule::Rule) = []

type Grammar
  rules::Dict{Symbol, Rule}
end

# empty rule is also accepted and never consumes
type EmptyRule <: Rule
  name
  action

  function EmptyRule()
    return new("")
  end

  function EmptyRule(name::AbstractString)
    return new(name)
  end
end

type ParserData
  # map of parsers to use
  parsers

  ParserData(parsers) = new(parsers)
end

function parseDefinition(name::AbstractString, sym::Symbol, pdata::ParserData)
  fn = get(pdata.parsers, sym, nothing)

  if fn !== nothing
    return fn(name, pdata, nothing)
  end

  # if not found, just return the symbol
  return sym
end

# general case, just return value
expand_names(value) = value

# if it's a symbol, check that it matches, and if so, convert it
function expand_names(sym::Symbol)
  m = match(r"_(\d+)", string(sym))
  if m !== nothing
    i = parse(Int, m.captures[1])

    return i == 0 ? :(value) : :(children[$i])
  end
  return sym
end

# if it's an expression, recursively go through tree and
# transform all symbols that match '_i'
function expand_names(expr::Expr)
  new_args = [expand_names(arg) for arg in expr.args]
  return Expr(expr.head, new_args...)
end

# FIXME: There is a weird mismatch in the parsers and parseDefinition..
# should they really be different? if not, parseDefinition might need
# to be changed to allow for arrays to be passed in
function parseDefinition(name::AbstractString, expr::Expr, pdata::ParserData)
  rule = EmptyRule()

  # if it's a macro (e.g. r"regex", then we want to expand it first)
  if expr.head === :macrocall
    # FIXME: using an evil eval
    rule = parseDefinition(name, eval(expr), pdata)
  elseif expr.head === :curly
    rule = parseDefinition(name, expr.args[1], pdata)
    rule.action = expand_names(expr.args[2])
  else
    parser = get(pdata.parsers, expr.args[1], nothing)

    if parser !== nothing
      rule = parser(name, pdata, expr.args[2:end])
    end
  end

  return rule
end

function collect_rules(rule::Rule, lst::Array)
  push!(lst, rule)

  for child in get_children(rule)
    collect_rules(child, lst)
  end

  return lst
end

function parseGrammar(grammar_name::Symbol, expr::Expr, pdata::ParserData)
  code = Any[]
  push!(code, :(rules = Dict()))

  all_rules = Rule[]

  for definition in expr.args[2:2:end]
    name = string(definition.args[1])
    ref_by_name = Expr(:ref, :rules, name)

    rule = parseDefinition(name, definition.args[2], pdata)
    push!(code, :(rules[$name] = $rule))
    all_rules = collect_rules(rule, all_rules)
  end

  for rule in all_rules
    if typeof(rule.action) != Function
      dot = Expr(:(.), rule, QuoteNode(:action))
      fn = Expr(:(->),
        Expr(:tuple, :rule, :value, :first, :last, :children),
        Expr(:block, rule.action))
      push!(code, Expr(:escape, Expr(:(=), dot, fn)))
    end
  end

  push!(code, :($(esc(grammar_name)) = Grammar(rules)))

  return Expr(:block, code...)
end


text:
module RiceBPH

using Agents: run!

Model = include("Model.jl")
Replication = include("Replication.jl")
PostProcess = include("PostProcess.jl")

function run_model(params, laststep=2880; seed=rand(1:2000))
    model = Model.init_model(; seed=seed, params...)
    return run!(
        model,
        Model.agent_step!,
        model.model_step!,
        laststep;
        adata=Model.adata,
        mdata=Model.mdata,
    )
end

function replication(params, n; seed_offset=0)
    return Replication.replication(
        Model.init_model,
        Model.agent_step!,
        Model.model_step!,
        2880,
        n;
        seed_offset=seed_offset,
        adata=Model.adata,
        mdata=Model.mdata,
        post_process=Model.post_process,
        params...,
    )
end

if !isempty(PROGRAM_FILE) && realpath(PROGRAM_FILE) == @__FILE__
    Dashboard = include("Dashboard.jl")
    try
        Dashboard.start()
    catch e
        showerror(stdout, e, catch_backtrace())
    end
end

end # module

max_stars_repo_name: ndgnuh/RiceBPH

max_issues_repo_name: ndgnuh/RiceBPH

max_forks_repo_name: ndgnuh/RiceBPH


text:
DMS_FSTRINGS = FormatExpr.([
    "{1:d} {2:d} {3}",
    "{1:d}:{2:d}:{3}",
    "{1:d}d{2:d}m{3}s",
    "{1:d}°{2:d}'{3}\"",
    "{1:d}°{2:d}′{3}″",
    "{1:d} {2:d}:{3}"
])

HMS_FSTRINGS = FormatExpr.([
    "{1:d} {2:d} {3}",
    "{1:d}:{2:d}:{3}",
    "{1:d}h{2:d}m{3}s",
    "{1:d}h{2:d}'{3}\"",
    "{1:d}h{2:d}′{3}″",
    "{1:d} {2:d}:{3}"
])

@testset "parsing" for fstring in DMS_FSTRINGS
    @test all(randdms(rng, 1000)) do dms
        d, m, s = dms
        str = format(fstring, d, m, s)
        res = parse_dms(str)
        @show dms
        @show str
        @show res
        return all(dms .≈ res)
    end

    @test all(randdms(rng, 1000)) do dms
        d, m, s = dms
        str = format(fstring, d, m, s)
        t1 = dms2rad(str) == parse_dms(str) |> dms2rad
        t2 = dms2deg(str) == parse_dms(str) |> dms2deg
        t3 = dms2ha(str) == parse_dms(str) |> dms2ha
        return t1 && t2 && t3
    end
end

@testset "no-ops" begin
    val = randn(rng)
    @test dms2rad(val) === val
    @test dms2deg(val) === val
    @test dms2ha(val) === val
    @test hms2rad(val) === val
    @test hms2deg(val) === val
    @test hms2ha(val) === val
end

@testset "parsing" for fstring in HMS_FSTRINGS
    @test all(randhms(rng, 1000)) do hms
        h, m, s = hms
        str = format(fstring, h, m, s)
        res = parse_hms(str)
        return all(hms .≈ res)
    end

    @test all(randhms(rng, 1000)) do hms
        h, m, s = hms
        str = format(fstring, h, m, s)
        t1 = hms2rad(str) == parse_hms(str) |> hms2rad
        t2 = hms2deg(str) == parse_hms(str) |> hms2deg
        t3 = hms2ha(str) == parse_hms(str) |> hms2ha
        return t1 && t2 && t3
    end
end
# test-macro-throws
# https://discourse.julialang.org/t/exceptions-in-macros-in-julia-0-7-1-0/14145/2
macro tmt(typ,expr)
    quote
       @test_throws $typ begin
          try
             $expr
          catch e
             rethrow(e.error)
          end
       end
    end
 end

@testset "macros" begin
    str = "12:37:34.2344"
    dms = parse_dms(str)
    @test all(dms .== [12, 37, 34.2344])
    @test dms"12:37:34.2344" == dms"12:37:34.2344"rad ≈ dms2rad(dms)
    @test dms"12:37:34.2344"deg ≈ dms2deg(dms)
    @test dms"12:37:34.2344"ha ≈ dms2ha(dms)
    @tmt ErrorException @dms_str("12:37:34.2344", "invalid")

    hms = parse_hms(str)
    @test all(hms .== [12, 37, 34.2344])
    @test hms"12:37:34.2344" == hms"12:37:34.2344"rad ≈ hms2rad(hms)
    @test hms"12:37:34.2344"deg ≈ hms2deg(hms)
    @test hms"12:37:34.2344"ha ≈ hms2ha(hms)
    @tmt ErrorException @hms_str("12:37:34.2344", "invalid")
end

@testset "negatives" begin
    @test parse_dms("-0:0:1") == (-0.0, 0.0, 1.0)
    @test parse_hms("-0:0:1") == (-0.0, 0.0, 1.0)
    @test parse_dms("- 0:0:1") == (-0.0, 0.0, 1.0)
    @test parse_hms("- 0:0:1") == (-0.0, 0.0, 1.0)
    @test parse_dms("+0:0:1") == (0.0, 0.0, 1.0)
    @test parse_hms("+0:0:1") == (0.0, 0.0, 1.0)
end

@testset "directions" begin
    @test parse_dms("1:0:0N") == (1.0, 0.0, 0.0)
    @test parse_dms("1:0:0S") == (-1.0, 0.0, 0.0)
    @test parse_dms("-1:0:0S") == (1.0, 0.0, 0.0)
    @test parse_dms("1:0:0E") == (1.0, 0.0, 0.0)
    @test parse_dms("1:0:0W") == (-1.0, 0.0, 0.0)
    @test parse_dms("-1:0:0W") == (1.0, 0.0, 0.0)

    @test parse_hms("1:0:0E") == (1.0, 0.0, 0.0)
    @test parse_hms("1:0:0W") == (-1.0, 0.0, 0.0)
    @test parse_hms("-1:0:0W") == (1.0, 0.0, 0.0)
    @test parse_dms("1:0:0N") == (1.0, 0.0, 0.0)
    @test parse_hms("1:0:0S") == (-1.0, 0.0, 0.0)
    @test parse_dms("-1:0:0S") == (1.0, 0.0, 0.0)
end

@testset "astropy examples" begin
    @test dms2deg("10.2345d") ≈ 10.2345
    @test dms2deg("1:2:30.43") ≈ 1.04178611
    @test hms2ha("1 2 0") ≈ 1.03333333
    @test dms2deg("1°2′3″") ≈ dms2deg("1°2′3″N") ≈ 1.03416667
    @test dms2deg("1d2m3.4s") ≈ -dms2deg("1d2m3.4sS") ≈ 1.03427778
    @test hms2ha("-1h2m3s") ≈ -hms2ha("-1h2m3sW") ≈ -1.03416667
    @test dms2deg((-1, 2, 3)) ≈ -1.03416667
end

@testset "partial representations" begin
    @test parse_dms("10.234") == (10.234, 0.0, 0.0)
    @test parse_dms("10.234 0.1") == (10.234, 0.1, 0.0)
    @test parse_dms("10.234 0.1 0.3") == (10.234, 0.1, 0.3)
    @test parse_hms("10.234") == (10.234, 0.0, 0.0)
    @test parse_hms("10.234 0.1") == (10.234, 0.1, 0.0)
    @test parse_hms("10.234 0.1 0.3") == (10.234, 0.1, 0.3)
end


text:
#!/Applications/Julia-0.6.1.app/Contents/Resources/julia/bin/julia
# Nearest Neighbor classifier, neural network based, and employing Hamming distances
# This classifier works with any mix of categorical or continuous variables, any number of classes,
# and can handle missing values.

__precompile__()

module HamNNNNclassifier

export main, train, test, readdatafile, printdatafiledescription, hamclass, generatetraintesttable, printattributeranking, anneal, glass

using NullableArrays, DataFrames, CategoricalArrays, CSV, Query, StatsBase, ArgParse

include("HamNNNNclassifierUtilities.jl")
include("HamNNNNclassifierio.jl")
include("HamNNNNclassifierReports.jl")

global version = v"0.1.1"
# v0.1.1 discontinues using the DataTables package, which has been deprecated. Uses DataFrames instead.
# also, need to use NullableArrays specifically, in order to have dropnull.

# Terminology:
# instances: cases; records; examples
#     instances can be grouped into sets, eg training set, test set, validation set
#     corresponds to a table row in a tabular data base
# Dataset
#     consists of a number of instances (cases, or examples)
# features: variables; fields; attributes
#     corresponds to a column in a tabular data base. The feature name is the column header
#     one of the features will be the class feature or class variable (also called the target variable, or
# 	class attribute)
#     features can be combined in various ways, to create new features
# feature extraction
#     the process of finding those features which provide the best classification performance
# classes
#     the set of values that the class feature can take
# parameters
#     eg, the k in k-nearest-neighbors
#     often, the work involved in applying machine learning to a problem is to find appropriate or optimal values for the parameters used by a given ML methodology
# Hamming distance (may also be called "overlap metric"
# bins
#     the number of bins or slices to be applied for a given continuous feature
# binning
#     the process of converting a continuous feature into a discrete feature
# discrete (as applied to a feature or variable)
#     nominal or ordinal data
# continuous (as applied to a feature or variable)
#     real-valued
#     ordinal data with a range greater than a certain parameter may also be treated as continous data

# NEWER ORANGE FORMAT:
# Prefixed attributes contain a one- or two-lettered prefix, followed by '#'
# and the attribute name. The first letter of the prefix can be:
# 	'm' for meta-attributes
# 	'i' to ignore the attribute
# 	'c' to define the class attribute
#
# the second letter of the prefix can be:
# 	'D' for discrete
# 	'C' for continuous
# 	'S' for string
# 	'B' for basket
# if there are no prefixes for an attribute (ie just the attribute name)
# then the attribute will be treated as discrete, unless the actual values
# are numbers, in which case it will be treated as continuous.
#
# OLDER ORANGE FORMAT:
# the information about variable type, etc is contained in two rows, following the row with
# attribute names:
# in the second line:
# 	'd' or 'discrete' or a list of values: denotes a discrete attribute
# 	'c' or 'continuous': denotes a continuous attribute
# 	'string' denotes a string variable, which we ignore
# 	'basket': these are continuous-valued meta attributes; ignore
# 	it may also contain a string of values separated by spaces. Use these
# 	as the values for a discrete attribute.
# the third line contains optional flags:
# 	'i' or 'ignore'
# 	'c' or 'class': there can only be one class attribute. If none is found,
# 	 use the last attribute as the class attribute.
# 	'm' or 'meta': meta attribute, eg weighting information; ignore
# 	'-dc' followed by a value: indicates how a don't care is represented.

# # Design Principles
# To work with the MLBase Julia package, it will be helpful to have functions which can be used
# as parameters for MLBase functions (see the MLBase documentation), eg for cross_validate,
# estfun –
# The estimation function, which takes a vector of training indices as input and returns a learned model, as:
# model = estfun(train_inds)
# evalfun –
# The evaluation function, which takes a model and a vector of testing indices as input and returns a score that indicates the goodness of the model, as
# score = evalfun(model, test_inds)

function hamclass()
	println("hamclass...")
end

function train()
	println("training...")
end

function test()
	println("testing...")
end

function parse_commandline()
	s = ArgParseSettings()

	@add_arg_table s begin
		"--opt1"
			help = "an option with an argument"
		"--opt2", "-o"
			help = "another option with an argument"
			arg_type = Int
			default = 0
		"--flag1"
			help = "an option without argument, i.e. a flag"
			action = :store_true
		"arg1"
			help = "a positional argument"
			required = true
	end

	return parse_args(s)
end

# function main()
# 	parsed_args = parse_commandline()
# 	println("Parsed args:")
# 	for (arg,val) in parsed_args
# 		println("  $arg  =>  $val")
# 	end
# end

paramDict = Dict

dataFile="/Users/henryolders/hammingnn/explore/test.tab"

function main()
	printdatafiledescription(dataFile)
	printattributeranking(dataFile)
	# train_x, train_codes, train_y = generatetraintesttable(dataFile)
	# @show train_x train_codes train_y
	train()
	test()
end

function anneal()
	printdatafiledescription("/Users/henryolders/hammingnn/explore/anneal.tab")
	printattributeranking("/Users/henryolders/hammingnn/explore/anneal.tab")
end

function glass()
	printdatafiledescription("/Users/henryolders/hammingnn/explore/glass.tab")
	printattributeranking("/Users/henryolders/hammingnn/explore/glass.tab")
end

# uncomment the line below, to run this module from the command line.
main()

end # module


text:
using StatisticalRethinking, CmdStan
#gr(size=(600,600));

ProjDir = rel_path("..", "scripts", "04")
cd(ProjDir)

howell1 = CSV.read(rel_path("..", "data", "Howell1.csv"), delim=';')
df = convert(DataFrame, howell1);

df2 = filter(row -> row[:age] >= 18, df);
first(df2, 5)

weightsmodel = "
data {
 int < lower = 1 > N; // Sample size
 vector[N] height; // Predictor
 vector[N] weight; // Outcome
}

parameters {
 real alpha; // Intercept
 real beta; // Slope (regression coefficients)
 real < lower = 0 > sigma; // Error SD
}

model {
 height ~ normal(alpha + weight * beta , sigma);
}

generated quantities {
}
";

stanmodel = Stanmodel(name="weights", monitors = ["alpha", "beta", "sigma"],model=weightsmodel,
  output_format=:mcmcchains);

heightsdata = Dict("N" => length(df2[:height]), "height" => df2[:height],
  "weight" => df2[:weight]);

rc, chn, cnames = stan(stanmodel, heightsdata, ProjDir, diagnostics=false,
  summary=false, CmdStanDir=CMDSTAN_HOME)

chn.value[1:5,:,1]

p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 4)
nvals = [10, 50, 150, 352];

for i in 1:length(nvals)
  N = nvals[i]
  heightsdataN = [
    Dict("N" => N, "height" => df2[1:N, :height], "weight" => df2[1:N, :weight])
  ]
  rc, chnN, cnames = stan(stanmodel, heightsdataN, ProjDir, diagnostics=false,
    summary=false, CmdStanDir=CMDSTAN_HOME)

  xi = 30.0:0.1:65.0
  rws, vars, chns = size(chnN)
  alpha_vals = convert(Vector{Float64}, reshape(chnN.value[:, 1, :], (rws*chns)))
  beta_vals = convert(Vector{Float64}, reshape(chnN.value[:, 2, :], (rws*chns)))

  p[i] = scatter(df2[1:N, :weight], df2[1:N, :height], leg=false,
    color=:darkblue, xlab="weight")
  for j in 1:N
    yi = alpha_vals[j] .+ beta_vals[j]*xi
    plot!(p[i], xi, yi, title="N = $N", color=:lightgrey)
  end
  scatter!(p[i], df2[1:N, :weight], df2[1:N, :height], leg=false,
    color=:darkblue, xlab="weight")
end
plot(p..., layout=(2, 2))

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
## Julia implementation of a ConditionalRBM from Graham Taylor's PhD Thesis

## Links:
##     Thesis - http://www.cs.nyu.edu/~gwtaylor/thesis/Taylor_Graham_W_200911_PhD_thesis.pdf
##     FCRBM - http://www.cs.toronto.edu/~fritz/absps/fcrbm_icml.pdf

# Input data layout for 2 batches and 3 steps of history
#
# NOTE: in the code below the history is referred to as `cond`
# since technically you can condition on any features not just
# previous time steps.
#
#      batch 1               batch 2
# |--------------------|--------------------|-----
# |  current visible   |  current visible   | ...
# |--------------------|--------------------|-----
# |  history step 1    |  history step 1    | ...
# |  history step 2    |  history step 2    | ...
# |  history step 3    |  history step 3    | ...
# |--------------------|--------------------|-----

import StatsBase: predict

@runonce mutable struct ConditionalRBM{T,V,H} <: AbstractRBM{T,V,H}
    W::Matrix{T}  # standard weights
    A::Matrix{T}  # autoregressive params (vis to vis)
    B::Matrix{T}  # hidden params(vis to hid)
    vbias::Vector{T}
    hbias::Vector{T}
    dyn_vbias::Array{T}
    dyn_hbias::Array{T}
end

function ConditionalRBM(T::Type, V::Type, H::Type,
                        n_vis::Int, n_hid::Int, n_cond::Int; sigma=0.01)
    ConditionalRBM{T,V,H}(
        map(T, rand(Normal(0, sigma), n_hid, n_vis)),
        map(T, rand(Normal(0, sigma), n_vis, n_cond)),
        map(T, rand(Normal(0, sigma), n_hid, n_cond)),
        zeros(T, n_vis),
        zeros(T, n_hid),
        zeros(T, n_vis),
        zeros(T, n_hid),
    )
end

function ConditionalRBM(V::Type, H::Type,
                        n_vis::Int, n_hid::Int, n_cond::Int; sigma=0.01)
    ConditionalRBM(Float64, V, H, n_vis, n_hid, n_cond;
                   sigma=sigma)
end

function ConditionalRBM(T::Type, V::Type, H::Type, n_vis::Int, n_hid::Int;
                        steps=5, sigma=0.01)
    ConditionalRBM(T, V, H, n_vis, n_hid, (n_vis * steps);
                   sigma=sigma)
end

function ConditionalRBM(V::Type, H::Type, n_vis::Int, n_hid::Int;
                        steps=5, sigma=0.01)
    ConditionalRBM(Float64, V, H, n_vis, n_hid, (n_vis * steps);
                   sigma=sigma)
end

function Base.show(io::IO, crbm::ConditionalRBM{T,V,H}) where {T,V,H}
    n_vis = size(crbm.vbias, 1)
    n_hid = size(crbm.hbias, 1)
    n_cond = size(crbm.A, 2)
    print(io, "ConditionalRBM{$T,$V,$H}($n_vis, $n_hid, $n_cond)")
end

function split_vis(crbm::ConditionalRBM, vis::Mat)
    vis_size = length(crbm.vbias)

    curr = vis[1:vis_size, :]
    cond = vis[(vis_size + 1):(vis_size + size(crbm.A, 2)), :]

    return curr, cond
end

function dynamic_biases!(crbm::ConditionalRBM, cond::Mat)
    crbm.dyn_vbias = crbm.A * cond .+ crbm.vbias
    crbm.dyn_hbias = crbm.B * cond .+ crbm.hbias
end

function hid_means(crbm::ConditionalRBM, vis::Mat)
    p = crbm.W * vis .+ crbm.dyn_hbias
    return logistic(p)
end

function vis_means(crbm::ConditionalRBM, hid::Mat)
    p = crbm.W' * hid .+ crbm.dyn_vbias
    return logistic(p)
end

function gradient_classic(crbm::ConditionalRBM{T}, X::Mat{T},
                       ctx::Dict) where T
    vis, cond = split_vis(crbm, X)
    sampler = @get_or_create(ctx, :sampler, persistent_contdiv)
    v_pos, h_pos, v_neg, h_neg = sampler(crbm, vis, ctx)
    n_obs = size(vis, 2)
    # updating weight matrix W
    dW = @get_array(ctx, :dW_buf, size(crbm.W), similar(crbm.W))
    # same as: dW = ((h_pos * v_pos') - (h_neg * v_neg')) / n_obs
    gemm!('N', 'T', T(1 / n_obs), h_neg, v_neg, T(0.0), dW)
    gemm!('N', 'T', T(1 / n_obs), h_pos, v_pos, T(-1.0), dW)

    # updating vis to vis matrix A
    dA = @get_array(ctx, :dA_buf, size(crbm.A), similar(crbm.A))
    # same as: dW = (h_pos * v_pos') - (h_neg * v_neg')
    gemm!('N', 'T', T(1 / n_obs), v_neg, cond, T(0.0), dA)
    gemm!('N', 'T', T(1 / n_obs), v_pos, cond, T(-1.0), dA)

    # updating hid to hid matrix A
    dB = @get_array(ctx, :dB_buf, size(crbm.B), similar(crbm.B))
    gemm!('N', 'T', T(1 / n_obs), h_neg, cond, T(0.0), dB)
    gemm!('N', 'T', T(1 / n_obs), h_pos, cond, T(-1.0), dB)

    # gradient for vbias and hbias
    db = dropdims(sum(v_pos, dims = 2) - sum(v_neg, dims = 2), dims = 2) ./ n_obs
    dc = dropdims(sum(h_pos, dims = 2) - sum(h_neg, dims = 2), dims = 2) ./ n_obs
    return dW, dA, dB, db, dc
end

function grad_apply_learning_rate!(crbm::ConditionalRBM{T},
                                   X::Mat{T},
                                   dtheta::Tuple, ctx::Dict) where T
    dW, dA, dB, db, dc = dtheta
    lr = T(@get(ctx, :lr, 0.1))
    # same as: dW *= lr
    scal!(length(dW), lr, dW, 1)
    scal!(length(dA), lr, dA, 1)
    scal!(length(dB), lr, dB, 1)
    scal!(length(db), lr, db, 1)
    scal!(length(dc), lr, dc, 1)
end

function grad_apply_momentum!(crbm::ConditionalRBM{T}, X::Mat{T},
                              dtheta::Tuple, ctx::Dict) where T
    dW, dA, dB, db, dc = dtheta
    momentum = @get(ctx, :momentum, 0.9)
    dW_prev = @get_array(ctx, :dW_prev, size(dW), zeros(T, size(dW)))
    # same as: dW += momentum * dW_prev
    axpy!(momentum, dW_prev, dW)
end

function grad_apply_weight_decay!(rbm::ConditionalRBM,
                                     X::Mat,
                                     dtheta::Tuple, ctx::Dict)
    # The decay penalty should drive all weights toward
    # zero by some small amount on each update.
    dW, dA, dB, db, dc = dtheta
    decay_kind = @get_or_return(ctx, :weight_decay_kind, nothing)
    decay_rate = @get(ctx, :weight_decay_rate,
                      throw(ArgumentError("If using :weight_decay_kind, weight_decay_rate should also be specified")))
    is_l2 = @get(ctx, :l2, false)
    if decay_kind == :l2
        # same as: dW -= decay_rate * W
        axpy!(-decay_rate, rbm.W, dW)
        axpy!(-decay_rate, rbm.A, dA)
        axpy!(-decay_rate, rbm.B, dB)
    elseif decay_kind == :l1
        # same as: dW -= decay_rate * sign(W)
        axpy!(-decay_rate, sign(rbm.W), dW)
        axpy!(-decay_rate, sign(rbm.A), dA)
        axpy!(-decay_rate, sign(rbm.B), dB)
    end
end

function grad_apply_sparsity!(rbm::ConditionalRBM{T}, X::Mat,
                              dtheta::Tuple, ctx::Dict) where T
    # The sparsity constraint should only drive the weights
    # down when the mean activation of hidden units is higher
    # than the expected (hence why it isn't squared or the abs())
    dW, dA, dB, db, dc = dtheta
    cost = @get_or_return(ctx, :sparsity_cost, nothing)
    vis, cond = split_vis(rbm, X)
    target = @get(ctx, :sparsity_target,
                  throw(ArgumentError("If :sparsity_cost is used, :sparsity_target should also be defined")))
    curr_sparsity = mean(hid_means(rbm, vis))
    penalty = T(cost * (curr_sparsity - target))
    add!(dW, -penalty)
    add!(dA, -penalty)
    add!(dB, -penalty)
    add!(db, -penalty)
    add!(dc, -penalty)
end

function update_weights!(crbm::ConditionalRBM, dtheta::Tuple, ctx::Dict)
    dW, dA, dB, db, dc = dtheta
    axpy!(1.0, dW, crbm.W)
    axpy!(1.0, dA, crbm.A)
    axpy!(1.0, dB, crbm.B)
    crbm.vbias += db
    crbm.hbias += dc
    # save previous dW
    dW_prev = @get_array(ctx, :dW_prev, size(dW), similar(dW))
    copyto!(dW_prev, dW)
end

function update_classic!(crbm::ConditionalRBM, X::Mat,
                            dtheta::Tuple, ctx::Dict)
    # apply gradient updaters. note, that updaters all have
    # the same signature and are thus composable
    grad_apply_learning_rate!(crbm, X, dtheta, ctx)
    grad_apply_momentum!(crbm, X, dtheta, ctx)
    grad_apply_weight_decay!(crbm, X, dtheta, ctx)
    grad_apply_sparsity!(crbm, X, dtheta, ctx)
    # add gradient to the weight matrix
    update_weights!(crbm, dtheta, ctx)
end

function free_energy(crbm::ConditionalRBM, vis::Mat)
    vb = sum(vis .* crbm.dyn_vbias, dims = 1)
    Wx_b_log = sum(log.(1 .+ exp.(crbm.W * vis .+ crbm.dyn_hbias)), dims = 1)
    result = - vb - Wx_b_log
    tofinite!(result)

    return result
end

function fit_batch!(crbm::ConditionalRBM, X::Mat, ctx = Dict())
    grad = @get_or_create(ctx, :gradient, gradient_classic)
    upd = @get_or_create(ctx, :update, update_classic!)
    curr, cond = split_vis(crbm, X)
    dynamic_biases!(crbm, cond)
    dtheta = grad(crbm, X, ctx)
    upd(crbm, X, dtheta, ctx)
    return crbm
end

function fit(crbm::ConditionalRBM{T}, X::Mat, opts::Dict{Any,Any}) where T
    @assert minimum(X) >= 0 && maximum(X) <= 1
    ctx = copy(opts)
    n_examples = size(X, 2)
    batch_size = @get(ctx, :batch_size, 100)
    batch_idxs = split_evenly(n_examples, batch_size)
    if @get(ctx, :randomize, false)
        batch_idxs = sample(batch_idxs, length(batch_idxs); replace=false)
    end
    n_epochs = @get(ctx, :n_epochs, 10)
    scorer = @get_or_create(ctx, :scorer, pseudo_likelihood)
    reporter = @get_or_create(ctx, :reporter, TextReporter())
    for epoch=1:n_epochs
        epoch_time = @elapsed begin
            for (batch_start, batch_end) in batch_idxs
                # BLAS.gemm! can't handle sparse matrices, so cheaper
                # to make it dense here
                batch = Matrix(X[:, batch_start:batch_end])
                batch = ensure_type(T, batch)
                fit_batch!(crbm, batch, ctx)
            end
        end
        curr, cond = split_vis(crbm, X)
        dynamic_biases!(crbm, cond)

        # We convert to a dense matrix, to avoid changing the the n_obs if
        # X is a sparse matrix
        score = scorer(crbm, Matrix(curr))
        report(reporter, crbm, epoch, epoch_time, score)
    end

    return crbm
end

fit(crbm::ConditionalRBM{T}, X::Mat; opts...) where {T} = fit(crbm, X, Dict{Any,Any}(opts))

function transform(crbm::ConditionalRBM{T}, X::Mat) where T
    curr, cond = split_vis(crbm, ensure_type(T, X))
    dynamic_biases!(crbm, cond)
    return hid_means(crbm, curr)
end

function generate(crbm::ConditionalRBM{T}, X::Mat; n_gibbs=1) where T
    curr, cond = split_vis(crbm, ensure_type(T, X))
    dynamic_biases!(crbm, cond)
    return gibbs(crbm, curr; n_times=n_gibbs)[3]
end

generate(crbm::ConditionalRBM{T}, vis::Vec; n_gibbs=1) where {T} =
    generate(crbm, reshape(ensure_type(T, vis), length(vis), 1); n_gibbs=n_gibbs)

function predict(crbm::ConditionalRBM{T}, cond::Mat; n_gibbs=1) where T
    cond = ensure_type(T, cond)
    @assert size(cond, 1) == size(crbm.A, 2)

    curr = view(cond, 1:length(crbm.vbias), :)
    vis = vcat(curr, cond)

    return generate(crbm, vis; n_gibbs=n_gibbs)
end

predict(crbm::ConditionalRBM{T}, cond::Vec; n_gibbs=1) where {T} =
    predict(crbm, reshape(ensure_type(T, cond), length(cond), 1); n_gibbs=n_gibbs)

predict(crbm::ConditionalRBM{T}, vis::Vec, cond::Vec{T}; n_gibbs=1) where {T} =
    predict(crbm, reshape(vis, length(vis), 1),
            reshape(ensure_type(T, cond), length(cond), 1); n_gibbs=n_gibbs)


text:
module SparsePaulis

# TODO - learn how to allow user to choose between Set and IntSet.

export SparsePauli

"""
For multiplying Paulis by complex numbers, we have a dict that
basically gives us the complex log:

    const PHASES = Dict{Complex, UInt8}(1 => 0, im => 1, -1 => 2, -im => 3)
"""
const PHASES = Dict{Complex, UInt8}(1 => 0, im => 1, -1 => 2, -im => 3)

"""
Parametric type representing sparse Pauli operators on an unknown
number of qubits. Uses two sets internally, `x_set` and `z_set`, as
well as an integer between 0 and 3 which represents the phase.

Written out, the entire n-qubit Pauli operator would look like this:

i^{`phase_int`} * prod_{j in `x_set`} X_j  * prod_{k in `z_set`} Z_j

Recommended reading: Dehaene and De Moor, PRA 2004
"""
type SparsePauli{T}
    x_set::Set{T}
    z_set::Set{T}
    phase_int::UInt8
    function SparsePauli{T}(x_set::Set{T}, z_set::Set{T}, phase_int::Integer) where T
        new(x_set, z_set, UInt8(mod(phase_int, 4)))
    end #function
end # type

function SparsePauli{T}(x_set::Set{T}, z_set::Set{T}) where T
    SparsePauli{T}(x_set, z_set, 0x00)
end #function

function SparsePauli{T}() where T
    SparsePauli{T}(Set{T}([]), Set{T}([]))
end #function

SparsePauli() = SparsePauli{Any}()

"""
Spits out an X/Y/Z depending on what sets elem is a member of.

Note: I've decided that elem has to be the same type as the
referenced SparsePauli, for performance.
"""
function char(elem::T, pauli::SparsePauli{T}) where T
    if in(elem, pauli.x_set)
        return in(elem, pauli.z_set) ? 'Y' : 'X'
    else
        return in(elem, pauli.z_set) ? 'Z' : 'I'
    end #if
end #function

"""
Copies a SparsePauli by returning a new SparsePauli whose x_set, z_set
and phase_int are all `copy`d.
"""
function copy(pauli::SparsePauli{T}) where T
    SparsePauli{T}(copy(pauli.x_set), copy(pauli.z_set), copy(phase_int))
end #function

end # module


text:
module TestAttitude

using Test
using LinearAlgebra
using Flight.Attitude
using Flight.Quaternions

export test_attitude

logrange(x1, x2, n) = (10^y for y in range(log10(x1), log10(x2), length=n))

function lograngesym(x1, x2, n)
    log_semi = collect(logrange(x1, x2, n))
    log_sym = [-reverse(log_semi); log_semi]
    return (i for i in log_sym)
end

function test_attitude()
    @testset verbose = true "Attitude" begin
        @testset verbose = true "RQuat" begin test_RQuat() end
        @testset verbose = true "RMatrix" begin test_RMatrix() end
        @testset verbose = true "RAxAng" begin test_RAxAng() end
        @testset verbose = true "REuler" begin test_REuler() end
    end
end

function test_RQuat()

        v_ab = [0.7, -1, 0.1, 1]
        v_bc = [1.7, 2, -1, 1]
        q_ab = RQuat(v_ab)
        q_bc = RQuat(v_bc)
        q_ab_neg = RQuat(-v_ab)
        x_c = [2, 3, -1.1]

    @testset "Constructors" begin
        #RQuat accepts any AbstractVector of length 4 and normalizes it, no need to
        #pass a UnitQuat instance
        @test q_ab._u ≈ RQuat(UnitQuat(v_ab))._u
        @test_throws DimensionMismatch RQuat([1,2,3])
        @test RQuat()._u[:] == [1,0,0,0]
    end

    @testset "Normalization" begin
        r_bad = RQuat(UnitQuat([1,2,3,4], normalization = false))
        @test abs(norm(r_bad._u)-1) > 1e-2
        @test norm(normalize(r_bad)) ≈ 1
        # normalize!(r_bad)
        # @test norm(r_bad) ≈ 1
    end

    @testset "Operators" begin
        #equality
        @test q_ab != q_bc
        @test q_ab == q_ab_neg
        @test q_ab == q_ab

        #isapprox
        @test q_ab ≈ q_ab_neg
        @test q_ab ≈ q_ab
        @test !(q_ab ≈ q_bc)

        #transformation, inversion & composition
        @test q_ab * (q_bc * x_c) ≈ (q_ab ∘ q_bc) * x_c
        @test x_c ≈ q_bc' * (q_bc * x_c)

        #callable notation
        @test q_ab(q_bc(x_c)) == q_ab * (q_bc * x_c)
        @test x_c ≈ q_bc'(q_bc(x_c))

        #time derivative
        ω_ab_b = [10, -4, 2]
        @test dt(q_ab, ω_ab_b) == 0.5 * (q_ab._u * FreeQuat(imag = ω_ab_b))

    end

end

function test_RAxAng()

        u_ab_bad = [5, -2, 1]
        u_ab = normalize(u_ab_bad)
        μ_ab = 1.214
        u_μ_ab = (u_ab, μ_ab)

        r_ab = RAxAng(u_ab, μ_ab)
        r_bc = RAxAng([-2,-3, -1], 0.72)
        x_c = [-1, 2, 3]

    @testset "Constructors" begin

        #inner constructor
        @test RAxAng(u_ab_bad, μ_ab, normalization = false).axis == u_ab_bad
        @test RAxAng(u_ab_bad, μ_ab).axis ≈ normalize(u_ab)

        #tuple constructor
        @test r_ab.axis == RAxAng(u_μ_ab).axis && r_ab.angle == RAxAng(u_μ_ab).angle

        #zero argument constructor
        @test RAxAng().axis == [1, 0, 0] && RAxAng().angle == 0

    end

    #these actually build upon quaternion conversion, but that's an
    #implementation detail we should not care about, so we can test these before
    @testset "Operators" begin

        #equality must account for different axis and angle sign combinations, and
        #2kπ increments
        @test RAxAng(u_ab, μ_ab) == RAxAng(u_ab, μ_ab)
        @test RAxAng(u_ab, μ_ab) ≈ RAxAng(-u_ab, -μ_ab)
        @test RAxAng(-u_ab, μ_ab) ≈ RAxAng(u_ab, -μ_ab)
        @test RAxAng(u_ab, μ_ab) ≈ RAxAng(u_ab, μ_ab + 4π)
        @test RAxAng(u_ab, 0) ≈ RAxAng([1, 2, 5], 0)
        @test RAxAng(u_ab, μ_ab) != RAxAng(u_ab, -μ_ab)
        @test !(RAxAng(u_ab, μ_ab) ≈ RAxAng(u_ab, -μ_ab))

        #approximate equality
        @test RAxAng(u_ab, μ_ab) ≈ RAxAng(u_ab, μ_ab + 1e-8)
        @test !(RAxAng(u_ab, μ_ab) ≈ RAxAng(u_ab, -μ_ab))

        #transformation, composition & inversion
        @test r_ab * (r_bc * x_c) ≈ (r_ab ∘ r_bc) * x_c
        @test x_c ≈ r_bc' * (r_bc * x_c)

    end

    @testset "Conversions" begin

        #conversions to and from RQuat (REuler tested in the REuler testset)

        #null rotation
        @test RQuat(RAxAng([1, 2, -1], 0)) ≈ RQuat([1,0,0,0])
        @test RQuat(RAxAng([3, 1, 1], 2π)) ≈ RQuat([1,0,0,0])

        #complete [-2π, 2π] angle range focusing on  small angles
        u = normalize([7, -5, 2])
        μ_array = lograngesym(1e-12, 2π, 10)
        r_array = collect(RAxAng(u, μ) for μ in μ_array)
        r_array_test = collect(r |> RQuat |> RAxAng for r in r_array)
        @test all(r_array .≈ r_array_test)

    end

    @testset "Mixed Operators" begin

        q_ab = RQuat(r_ab); q_bc = RQuat(r_bc)

        #approximate equality
        @test r_ab ≈ q_ab #some loss of precision should be expected
        @test !(r_ab ≈ q_bc) #some loss of precision should be expected

        #transformation, composition & inversion
        @test q_ab(r_bc(x_c)) ≈ r_ab(q_bc(x_c))
        @test (q_ab ∘ r_bc)(x_c) ≈ (r_ab ∘ q_bc)(x_c)
        @test r_ab' ≈ q_ab'
    end
end

function test_REuler()

    @testset "Constructors" begin

        ψ = 1.459; θ = -1.122; φ = 0.454
        r1 = REuler(ψ, θ, φ) #individual arguments
        r2 = REuler((ψ, θ, φ)) #tuple
        r3 = REuler(θ = θ) #tuple
        r4 = REuler() #zero argument
        @test r1.ψ == r2.ψ && r1.θ == r2.θ && r1.φ == r2.φ
        @test r3.θ == θ && r3.ψ == r3.φ == 0
        @test r4.ψ == r4.θ == r4.φ == 0

        @test_throws AssertionError REuler(θ = 1.7)
        @test_throws AssertionError REuler(ψ = 1.65π)
        @test_throws AssertionError REuler(φ = 1.65π)

    end

    @testset "Operators" begin

        r_ab = REuler(-1.234, √2, √3)
        r_bc = REuler(1.34, -√2, √2)
        x_c = [-1, 2, 3]

        #equality
        @test r_ab == r_ab
        @test r_ab != r_bc

        #approximate equality
        @test r_bc ≈ REuler(1.34000001, -√2, √2)
        @test REuler(ψ = π) ≈ REuler(ψ = -π)
        @test REuler(φ = π) ≈ REuler(φ = -π)
        @test !(r_ab ≈ r_bc)

        #transformation, composition & inversion
        @test r_ab * (r_bc * x_c) ≈ (r_ab ∘ r_bc) * x_c
        @test x_c ≈ r_bc' * (r_bc * x_c)

    end

    @testset "Conversions" begin

        #conversions to and from RQuat and RAxAng

        #null rotation
        @test RQuat(REuler()) ≈ RQuat([1,0,0,0])

        #cover Euler angle ranges including interval bounds but avoiding gimbal lock
        ψ_range = range(-π, π, length = 10)
        θ_range = range(-(π/2 - 0.001), π/2 - 0.001, length = 10)
        φ_range = range(-π, π, length = 10)
        r_array = vec([REuler(i) for i in Iterators.product(ψ_range, θ_range, φ_range)])

        r_array_test = [r |> RQuat |> REuler |> RAxAng |> REuler for r in r_array]
        # r_failed = r_array[r_array .≈ r_array_test]
        # r_test_failed = r_array_test[r_array .≈ r_array_test]
        @test all(r_array .≈ r_array_test)

    end

    @testset "Mixed Operators" begin

        #equality
        r_ab = REuler(0.1, -1, 0.6); r_bc = REuler(-1.8, 0.2, -1)
        q_ab = RQuat(r_ab); q_bc = RQuat(r_bc)
        x_c = [2,3,10]

        #approximate equality
        @test r_ab ≈ q_ab #some loss of precision should be expected
        @test !(r_ab ≈ q_bc) #some loss of precision should be expected

        #transformation, composition & inversion
        @test q_ab * (r_bc * x_c) ≈ r_ab * (q_bc * x_c)
        @test (q_ab ∘ r_bc) * x_c ≈ (r_ab ∘ q_bc) * x_c
        @test r_ab' ≈ q_ab'
    end

end

function test_RMatrix()

    A = [0.590969 0.288519 0.902137;
        0.657604 0.927693 0.218535;
        0.370707  0.481724  0.10268]

    A_orth = qr(A).Q
    r_ab = RMatrix(A_orth)
    r_bc = RMatrix(sin.(A_orth)) #generate a different matrix and let the constructor normalize
    x_c = [-1, 2, 3]

    @testset "Constructors" begin

        A_dist = A_orth + 1e-6*sin.(A)
        r_orth = RMatrix(A_orth, normalization = false)
        r_dist = RMatrix(A_dist, normalization = false)
        @test det(r_orth) ≈ 1
        @test !(det(r_dist) ≈ 1)
        @test det(RMatrix(A_dist)) ≈ 1 #automatic normalization on construction

        r_norm = normalize(r_dist)
        @test det(r_norm) ≈ 1
        # r_norm_copy = RMatrix(A_dist, normalization = false)
        # normalize!(r_norm_copy)
        # @test det(r_norm_copy) ≈ 1

        #error after renormalization
        r_err = r_orth ∘ r_norm'
        @test RAxAng(r_err).angle < 1e-7

    end

    @testset "Operators" begin

        #equality
        @test r_ab == r_ab
        @test r_ab != r_bc

        #approximate equality
        @test r_ab ≈ RMatrix(r_ab._mat .+ 1e-10, normalization = false)
        @test !(r_ab ≈ r_bc)

        #transformation, composition & inversion
        @test r_ab * (r_bc * x_c) ≈ (r_ab ∘ r_bc) * x_c
        @test r_ab * (r_bc * x_c) ≈ (r_ab * r_bc) * x_c
        @test x_c ≈ r_bc' * (r_bc * x_c)

        #time derivative
        ω_ab_b = [10, -4, 2]
        @test dt(r_ab, ω_ab_b) == r_ab._mat * Attitude.skew(ω_ab_b)
    end

    @testset "Conversions" begin

        #cover Euler angle ranges including interval bounds but avoiding gimbal lock
        ψ_range = range(-π, π, length = 10)
        θ_range = range(-(π/2 - 0.001), π/2 - 0.001, length = 10)
        φ_range = range(-π, π, length = 10)
        r_array = vec([REuler(i) for i in Iterators.product(ψ_range, θ_range, φ_range)])

        #this should exercise all cases of RMatrix -> RQuat conversion
        r_array_test = [r |> RMatrix |> RQuat |> RMatrix |> REuler |> RMatrix |> RAxAng |> RMatrix for r in r_array]
        @test all(r_array .≈ r_array_test)

    end

    @testset "Mixed Operators" begin

        #equality
        q_ab = RQuat(r_ab)
        q_bc = RQuat(r_bc)

        #approximate equality
        @test r_ab ≈ q_ab #some loss of precision should be expected
        @test !(r_ab ≈ q_bc) #some loss of precision should be expected

        #transformation, composition & inversion
        @test q_ab * (r_bc * x_c) ≈ r_ab * (q_bc * x_c)
        @test (q_ab ∘ r_bc) * x_c ≈ (r_ab ∘ q_bc) * x_c
        @test r_ab' ≈ q_ab'
    end

end

end #module