text:
module MutableStrings

using Compat

import Base.convert, Base.promote_rule, Base.endof, Base.getindex, Base.setindex!, Base.sizeof, Base.search, Base.rsearch, Base.string, Base.print, Base.write, Base.map!, Base.reverse!, Base.matchall, Base.next, Base.length, Base.isvalid

export MutableASCIIString, MutableUTF8String, MutableString, uppercase!, lowercase!, replace!, map!, ucfirst!, lcfirst!, reverse!, matchall, search, rsearch, getindex, setindex!
export convert, promote_rule, next, length, endof, sizeof, isvalid, string

immutable MutableASCIIString <: DirectIndexString
    data::Array{Uint8,1}
end

immutable MutableUTF8String <: String
    data::Array{Uint8,1}
end

typealias MutableString Union(MutableASCIIString, MutableUTF8String)

sizeof{T<:MutableString}(s::T) = sizeof(s.data)

matchall(re::Regex, str::MutableString, overlap::Bool=false) = matchall(re, convert(UTF8String, str), overlap)

replace!(s::MutableString, pat, r) = replace!(s, pat, r, 0)

# won't be required if MutableASCIIString is made part of ByteString
search(str::MutableString, re::Regex) = search(str, re, 1)
function search(str::MutableString, re::Regex, idx::Integer)
    if idx > nextind(str,endof(str))
        throw(BoundsError())
    end
    opts = re.options & Base.PCRE.EXECUTE_MASK
    Base.compile(re)
    m, n = exec(re.regex, re.extra, str, idx-1, opts, true)
    isempty(m) ? (0:-1) : ((m[1]+1):prevind(str,m[2]+1))
end

# won't be required if MutableASCIIString is made part of ByteString
exec(regex::Ptr{Void}, extra::Ptr{Void}, str::MutableString, offset::Integer, options::Integer, cap::Bool) = exec(regex, extra, str, 0, offset, sizeof(str), options, cap)

function exec(regex::Ptr{Void}, extra::Ptr{Void}, str::MutableString, shift::Integer, offset::Integer, len::Integer, options::Integer, cap::Bool)
    if offset < 0 || len < offset || len+shift > sizeof(str)
        error(BoundsError)
    end
    ncap = Base.PCRE.info(regex, extra, Base.PCRE.INFO_CAPTURECOUNT, Int32)
    ovec = Array(Int32, 3(ncap+1))
    n = ccall((:pcre_exec, :libpcre), Int32,
              (Ptr{Void}, Ptr{Void}, Ptr{Uint8}, Int32, Int32, Int32, Ptr{Int32}, Int32),
              regex, extra, pointer(str.data,shift+1), len, offset, options, ovec, length(ovec))
    if n < -1
        error("error $n")
    end
    cap ? ((n > -1 ? ovec[1:2(ncap+1)] : Array(Int32,0)), ncap) : n > -1
end

print{T <: MutableString}(io::IO, s::T) = (write(io, s);nothing)
write{T <: MutableString}(io::IO, s::T) = write(io, s.data)

convert(::Type{UTF8String},         s::MutableASCIIString) = UTF8String(s.data)
convert(::Type{ASCIIString},        s::MutableASCIIString) = ASCIIString(s.data)
convert(::Type{MutableASCIIString}, s::ASCIIString)        = MutableASCIIString(copy(s.data))
convert(::Type{MutableASCIIString}, s::UTF8String)         = convert(MutableASCIIString, copy(s.data))
convert(::Type{MutableASCIIString}, a::Array{Uint8,1})     = is_valid_ascii(a) ? MutableASCIIString(a) : error("invalid ASCII sequence")
convert(::Type{MutableASCIIString}, s::MutableASCIIString) = s
convert(::Type{MutableASCIIString}, s::AbstractString)     = convert(MutableASCIIString, bytestring(s))

convert(::Type{UTF8String},         s::MutableUTF8String)  = UTF8String(s.data)
convert(::Type{MutableUTF8String},  s::ByteString)         = MutableUTF8String(copy(s.data))
convert(::Type{MutableUTF8String},  a::Array{Uint8,1})     = is_valid_utf8(a) ? MutableUTF8String(a) : error("invalid UTF8 sequence")
convert(::Type{MutableUTF8String},  s::MutableUTF8String)  = s
convert(::Type{MutableUTF8String},  s::AbstractString)     = convert(MutableUTF8String, bytestring(s))

promote_rule(::Type{UTF8String},  ::Type{MutableUTF8String})  = UTF8String
promote_rule(::Type{UTF8String},  ::Type{MutableASCIIString}) = UTF8String
promote_rule(::Type{ASCIIString}, ::Type{MutableASCIIString}) = ASCIIString

include("mutable_ascii.jl")
include("mutable_utf8.jl")

end


text:
using Replay

@assert !Base.get_have_color() "please run julia with --color=no e.g. `julia --color=no examples/disable_color/app.jl`"

instructions = [
    """
    function greet(msg::String="Hello!")
        println(msg)
    end
    """,
    "greet()",
    """greet("Hello, World!")""",
]

replay(instructions, use_ghostwriter=false, cmd="--color=no")
replay(instructions, use_ghostwriter=true, cmd="--color=no")

ext: jl


text:
using Test, gmsh

@testset "test_t1.jl" begin
    include("../examples/t1.jl")
    @test filesize("t1.msh") > 0
end
@testset "test_t2.jl" begin
    include("../examples/t2.jl")
    @test filesize("t2.msh") > 0
end
@testset "test_t3.jl" begin
    include("../examples/t3.jl")
    @test filesize("t3.msh") > 0
end
@testset "test_t4.jl" begin
    include("../examples/t4.jl")
    @test filesize("t4.msh") > 0
end
@testset "test_t5.jl" begin
    include("../examples/t5.jl")
    @test filesize("t5.msh") > 0
end
#@testset "test_t10.jl" begin
#    include("../examples/t10.jl")
#    @test filesize("t10.msh") > 0
#end
#@testset "test_t16.jl" begin
#    include("../examples/t16.jl")
#    @test filesize("t16.msh") > 0
#end


text:
# Autogenerated wrapper script for Osi_jll for x86_64-w64-mingw32-cxx03
export libOsi, libOsiCommonTests

using CoinUtils_jll
using OpenBLAS32_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Osi")
JLLWrappers.@declare_library_product(libOsi, "libOsi-1.dll")
JLLWrappers.@declare_library_product(libOsiCommonTests, "libOsiCommonTests-1.dll")
function __init__()
    JLLWrappers.@generate_init_header(CoinUtils_jll, OpenBLAS32_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libOsi,
        "bin\\libOsi-1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libOsiCommonTests,
        "bin\\libOsiCommonTests-1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#**************************************************************************************
# Config.jl
# =============== part of the GeoEfficiency.jl package.
#
# this file contain the configurable parameters of the program.
#
#**************************************************************************************

#------------------ ---- set_global_consts -----------------------------

# source_file                                                   # default value

##Physics_Model.jl##
##Input_Consolejl##

##Input_Batch.jl##
const dataFolder    = string(".", @__MODULE__)                 # string(".", @__MODULE__) = ".GeoEfficiency"
const dataDir       = joinpath(homedir(), dataFolder)          # joinpath(homedir(), dataFolder)

##Calculations.jl##
const integrate     = begin using QuadGK; QuadGK.quadgk; end    # begin using QuadGK; QuadGK.quadgk; end
const relativeError = 1.0E-4                                    # 1.0E-4
const absoluteError = eps(1.0)                                  # eps(1.0)

##Output_Interface.jl##
const resultsFolder = "results"                                 # "results"

"""
set the default value for the global variable `_max_batch`
"""
const max_display   =  20                                       # 20


text:
using Pidoh
using Test
using SafeTestsets

@safetestset "BitString Tests" begin
    include("bitstringtests.jl")
end
@safetestset "Instance Tests" begin
    include("instancetests.jl")
end
@safetestset "Simulated Annealing Tests" begin
    include("satests.jl")
end

@safetestset "Correctness of Algorithms Tests" begin
    include("algorithmstests.jl")
end

@testset "Pidoh.jl" begin
    # Write your tests here.
end

max_stars_repo_name: DTUComputeTONIA/SDRLS


text:
######
# Data-Driven Sets for JuMPeR
######

module DDUS

# The JuMPeR oracle interface
using JuMP
using Compat
import JuMPeR
import JuMPeR: AbstractOracle, registerConstraint, setup, generateCut, generateReform

include("helpers.jl")
export suppFcn
include("FBOracle.jl")
include("UIOracle.jl")
include("UCSOracle.jl")
include("LCXOracle.jl")
include("UMOracle.jl")
#include("UDYOracle.jl")  #Needs SDP support in JuMPeR ...  add later

# Implements the JuMPeR oracle interface for all oracles
include("common.jl")

end # module


text:
# push!(LOAD_PATH, abspath(@__DIR__,".."))

using LinearAlgebra,Statistics,StatsBase,Distributions
using Plots,NamedColors ; theme(:dark) #; plotlyjs();
using SparseArrays
using SpikingRNNs; const global S = SpikingRNNs
using Test
using BenchmarkTools
using ProgressMeter
import FileIO: save # save raster png

using InvertedIndices

using Random; Random.seed!(0)

##

using Plots ; theme(:dark)
w = [ 1. , 1.2, 5, 10, 11.] .+ 1.0
wtemp = 100.0

probs = @. exp(-w/wtemp)
probs ./= sum(probs)

bar(w,probs;leg=false,color=:white)

##

##
dt = 0.1E-3
myτe = 20E-3 # seconds
myτi = 20E-3 # seconds
τrefr= 1E-3 # refractoriness
vth_e = -20.0   # mV
vthexp = -52.0 # actual threshold for spike-generation
vth_i = vthexp
eifslope = 2.0
Cap = 300.0 #capacitance mF
v_rest_e = -60.0
v_rest_i = -60.0
v_rev_e = 0.0
v_rev_i = -75.0
v_leak_e = v_rest_e
v_leak_i = v_rest_i
v_reset_e = v_rest_e
v_reset_i = v_rest_i

# synaptic kernel
taueplus = 6E-3 #e synapse decay time
taueminus = 1E-3 #e synapse rise time

tauiplus = 2E-3 #i synapse decay time
tauiminus = 0.5E-3 #i synapse rise time

const Ne = 200
nt_e = let sker = S.SKExpDiff(taueplus,taueminus)
  sgen = S.SpikeGenEIF(vthexp,eifslope)
  S.NTLIFConductance(sker,sgen,myτe,Cap,
   vth_e,v_reset_e,v_rest_e,τrefr,v_rev_e)
end
ps_e = S.PSLIFConductance(nt_e,Ne)

# step like function for half the population

function genfun(t::Real,i::Integer)
  min = 100.0
  max = 200.0
  period = 1.0
  if i>100
    return min
  else
    trem = rem(t,period)
    return trem < 0.5period ? min : max
  end
end
function genfun_upper(t::Real,i::Integer)
  min = 100.0
  max = 300.0
  if i>100
    return min
  else
    return max
  end
end

nt_in = let
  sker = S.SKExpDiff(taueplus,taueminus)
  sgen = S.SGPoissonFExact(genfun,genfun_upper)
  S.NTInputConductance(sgen,sker,v_rev_e)
end
in_weight = 10.0
ps_in = S.PSInputPoissonConductanceExact(nt_in,in_weight,Ne)

##
pop_e = S.Population(ps_e,(S.FakeConnection(),ps_in))
ntw = S.RecurrentNetwork(dt,pop_e)

Ttot = 4.0
# record spiketimes and internal potential
krec = 1
n_e_rec = Ne
t_wup = 0.0
rec_state_e = S.RecStateNow(ps_e,krec,dt,Ttot;idx_save=collect(1:n_e_rec),t_warmup=t_wup)
rec_spikes_e = S.RecSpikes(ps_e,500.0,Ttot;idx_save=collect(1:n_e_rec),t_warmup=t_wup)

## Run

times = (0:ntw.dt:Ttot)
nt = length(times)
# clean up
S.reset!.([rec_state_e,rec_spikes_e])
S.reset!.([ps_e,ps_in])
# initial conditions
ps_e.state_now .= v_reset_e

@time begin
  @showprogress 5.0 "network simulation " for (k,t) in enumerate(times)
    rec_state_e(t,k,ntw)
    rec_spikes_e(t,k,ntw)
    S.dynamics_step!(t,ntw)
  end
end

S.add_fake_spikes!(1.0vth_e,rec_state_e,rec_spikes_e)
myrast = S.raster_png(0.01,rec_spikes_e;Nneurons=Ne,Ttot=Ttot)
save("/tmp/rast.png",myrast)


text:
using kNNRecommenders
using Test

function test()
    println("Loading MovieLens 1M data...")
    mldata = MovieLens()
    load_data(mldata)
    split_data(mldata, random_state=220)

    transform_to_rating_matrix(mldata, 5)

    println("Performing ItemkNN")
    knn = kNNRecommender(260, 0.0, :tfidf)
    knn.rating_matrix = mldata.rating_matrix_train
    compute_similarity_matrix(knn)

    recalls = evaluate_recall(knn, 20, mldata.rating_matrix_test)
    mean_racall = sum(recalls) / length(recalls)
    @show mean_racall

    return mean_racall ≈ 0.2774901060050168
end

@test test()


text:
#Run an all-pair shortest path using dijkstra, minimizing time and not costs
function shortestPaths!(pb::TaxiProblem)
    pb.sp = shortestPaths(pb.network, pb.roadTime, pb.roadCost)
end

function shortestPaths(n::Network, roadTime::SparseMatrixCSC{Float64, Int},
                                   roadCost::SparseMatrixCSC{Float64, Int})

  nLocs  = length( vertices(n))
  pathTime = Array(Float64, (nLocs,nLocs))
  pathCost = Array(Float64, (nLocs,nLocs))
  previous = Array(Int, (nLocs,nLocs))

  for i in 1:nLocs
    parents, dists, costs = custom_dijkstra(n,i, roadTime, roadCost)
    previous[i,:] = parents
    pathTime[i,:] = dists
    pathCost[i,:] = costs
  end
  return ShortPaths(pathTime, pathCost, previous)
end

#Compute the table of the next locations on the shortest paths
#next[i, j] = location after i when going to j
function nextLoc(n::Network, sp::ShortPaths, roadTime::SparseMatrixCSC{Float64, Int})
  nLocs = size(sp.previous,1)
  next = Array(Int, (nLocs,nLocs))
  for i in 1:nLocs, j in 1:nLocs
    if i == j
      next[i,i] = i
    else
      minTime = Inf
      mink = 0
      for n in out_neighbors(n,i)
        if roadTime(i,n) + sp.traveltime[n,j] < minTime
          mink = n
          minTime = roadTime(i,n) + sp.traveltime[n,j]
        end
      end
      next[i,j] = mink
    end
  end
  return next
end

# TEMPORARY FIX
# Dijkstra algorithm

function custom_dijkstra(
    g::SimpleGraph,
    src::Int,
    edge_dists::AbstractArray{Float64, 2},
    edge_costs::AbstractArray{Float64, 2})

    # find number of vertices
    nvg = nv(g)
    # initialize return types
    dists = fill(typemax(Float64), nvg)
    costs = fill(typemax(Float64), nvg)
    parents = zeros(Int, nvg)
    visited = falses(nvg)
    # Create mutable binary heap and associated hashmap
    h = DijkstraEntry{Float64}[]
    sizehint(h, nvg)
    H = mutable_binary_minheap(h)
    hmap = zeros(Int, nvg)
    dists[src] = 0.0
    costs[src] = 0.0
    # Add source node to heap
    ref = push!(H, DijkstraEntry{Float64}(src, dists[src], costs[src]))
    hmap[src] = ref
    # As long as all edges have not been explored
    while !isempty(H)
        # Retrieve closest element to source
        hentry = pop!(H)
        u = hentry.vertex
        # Look at all its neighbors and update relevant distances if necessary
        for v in out_neighbors(g,u)
            if dists[u] == typemax(Float64)
                alt = typemax(Float64)
                alt2 = typemax(Float64)
            else
                alt = dists[u] + edge_dists[u,v]
                alt2 = costs[u] + edge_costs[u,v]
            end
            if !visited[v]  # If a vertex has never been visited, push it to the heap
                dists[v] = alt
                costs[v] = alt2
                parents[v] = u
                visited[v] = true
                ref = push!(H, DijkstraEntry{Float64}(v, alt, alt2))
                hmap[v] = ref
            else            # If a vertex has been visited, decrease its key if distance estimate decreases
                if alt < dists[v]
                    dists[v] = alt
                    costs[v] = alt2
                    parents[v] = u
                    update!(H, hmap[v], DijkstraEntry{Float64}(v, alt, alt2))
                end
            end
        end
    end

    dists[src] = 0
    costs[src] = 0.0
    parents[src] = 0

    return parents, dists, costs
end

max_stars_repo_name: sebmart/NYC-timing

max_issues_repo_name: sebmart/NYC-timing

max_forks_repo_name: sebmart/NYC-timing


text:
function VestaTheme()
    Theme(
        ambient   = Vec3f0(0.1),
        diffuse   = Vec3f0(0.9),
        specular  = Vec3f0(0.8),
        shininess = 32f0
    )
end

function setup_axis!(scene, obj)
    if obj isa Node
        throw(MethodError(setup_axis!, (scene, obj)))
    else
        return setup_axis!(scene, const_lift(obj))
    end
end
function setup_axis!(ax::Axis, bs::Node{BS}) where {BS <: AbstractBandstructure}
    ax.xlabel = "Momentum"
    ax.ylabel = "Energy"
    on(bs) do bs
        k_point_labels = label.(path(bs))
        Es = energies(bs)
        ticks = cumsum(vcat(1, [length(Es[i][1]) for i in eachindex(Es)]))
        ax.xticks[] = (ticks, k_point_labels)
    end
    # Trigger update of ticks
    bs[] = bs[]
    ax.ygridvisible = false

    ax
end
function setup_axis!(ls::LScene, bz::Node{BZ}) where {BZ <: AbstractBrillouinZone}
    setup_axis!(ls.scene, bz)

end
function setup_axis!(s::Scene, bz::Node{BZ}) where {BZ <: AbstractBrillouinZone}
    s.show_axis = false
    s.update_limits = false
    on(bz) do bz
        min, max = extrema(corners(bz))
        center = mean(corners(bz))
        width = 1.2maximum(max .- min)
        s.limits[] = FRect(Point((center .- 0.5width)...), Vec{length(center)}(width))
    end
    bz[] = bz[]
end


text:
"""
This file checks whether QuadGK.quadgk() can integrate numerically to arbitrary precision using BigFloat types. To check, I look at the integral:

    ∫_0^1 cosh(xt) dt

which has a known closed analytic form:

    sinh(x)/x,  x!=0

which we can use for comparison. x has to be parsed as a BigFloat type using
BigFloat("\$x") (remove space in actual code) NOT BigFloat(x) which would just convert a Float64 type to a BigFloat (and so would only be accurate to 64 bits). To check that passing BigFloat("\$x") into sinh(x)/x gives an accurate arbitrary precision answer (and not plagued by rounding errors), I checked with wolframalpha answers which does produce arb. prec. answers and they matched and so passing BigFloat("\$x") into sinh(x)/x gives a good check of whether quadgk above works to arb. prec. too.

I confirmed that quadgk can give arb prec answers (I just using it incorrectly previously to give Float64 accurate results), so now using the correct method, I convert hyperbolic_integral_two() from check_hyperbolic_integral.jl:

    β/2 ∫_0^1 (1 - cosh(Ωβ(1-x)/2)) cosh(xβ/2) / (a^2 - β^2 x^2/4 - bcosh(vβx/2))^{3/2} dx

into arb prec version and see what it gives. Spoilers: it's super slow, further confirmation that it's now arb prec.
"""

"""
This function just evaluates the integral:

    ∫_0^1 cosh(xt) dt

using quadgk. To make it arbitrary precision, the limits have to be parsed properly into a BigFloat type using BigFloat("#") and the absolute tolerance of the integration has to be set to the machine accuracy of your required precision. This is obtained using eps(x::BigFloat) (e.g. eps(Float64) ~ 2.22e-16) and set with atol = eps(x::BigFloat). The precision of the BigFloats can then be changed using setprecision(BigFloat, precision::Int) externally.
"""
function cosh_integral(x)
    err = eps(x)
    integrand(t) = cosh(x * t)
    integral = quadgk(t -> integrand(t), BigFloat("0"), BigFloat("1"), atol = err)
end

"""
Just the exact closed analytic answer to the cosh integral used for comparison. x has to be parsed as a BigFloat("#") before passing into the function to evaluate it at higher precision properly.
"""
function exact_form(x)
    sinh(x) / x
end

"""
Check outputs of the two above functions.
"""

setprecision(BigFloat, 128) # set the precision of the BigFloat types.
x_range = BigFloat("5")

# for bits = 128 and x_range = BigFloat("5") if get:
# 0.005659 seconds (13.49 k allocations: 768.068 KiB)
@time exact = [exact_form(x) for x in x_range]
# exact = fill(14.84064211555775179540189439921291311993)
@show(exact)
# 0.374177 seconds (1.00 M allocations: 46.879 MiB, 3.29% gc time)
@time integral = [cosh_integral(x) for x in x_range]
# fill((14.84064211555775179540189439921291311988, 2.350988701644575015937473074444491355637e-38))
@show(integral)
# Can see that the integral uses large memory allocations. Last two digits differing is common for BigFloat precisions (natural rounding error from float arithmatic) but we can see that the answers match to 128 bit precision.

"""
Using the arbitrary precision quadgk method used above, apply this to the hyperbolic integral I have been trying to solve i.e. hyperbolic_integral_one() in check_hyperbolic_integrals.jl. I actually use hyperbolic_integral_two() here as the change of variables to the limits [0,1] makes it easier for quadgk to evaluate at higher βs.

This allocates a lot of memory to evaluating the integral. So I did use @code_typewarn to check that there are no type-stabilities, but it did not find any so I think this large memory allocation is just a result of small atol set by the machine precision. For 128 bits this is atol ~ 1.93e-34.
"""
function hyperbolic_integral_two(Ω, β, α, v, w)

    err = eps(Ω)

    # Initialise constants.
    R = (v^2 - w^2) / (w^2 * v)
    a = sqrt(β^2 / 4 + R * β * coth(β * v / 2))
    b = R * β / sinh(β * v / 2)

    coefficient = 2 * α * β^(3 // 2) * v^3 / (3 * sqrt(π) * sinh(β / 2) * w^3)

    integrand(x) = (1 - cosh(Ω * β * (1 - x) / 2)) * cosh(x * β / 2) / (a^2 - β^2 * x^2 / 4 - b * cosh(v * β * x / 2))^(3 // 2)

    integral = quadgk(x -> integrand(x), BigFloat("0.0"), BigFloat("1.0"), atol = err)

    return coefficient * β * integral[1] / 2
end

"""
Check output and timing of the above function.
"""

setprecision(BigFloat, 128)
Ω_range = [BigFloat("12.01")]
β = BigFloat("4.0")
α = BigFloat("7.0")
v = BigFloat("5.8")
w = BigFloat("1.6")
# For Ω_range = [BigFloat("12.01")], β = BigFloat("4.0"), α = BigFloat("7.0"), v = BigFloat("5.8"), w = BigFloat("1.6") I get:
# 209.267907 seconds (1.02 G allocations: 37.016 GiB, 5.82% gc time)
@time hyp = hyperbolic_integral_two.(Ω_range, β, α, v, w)
# hyp = BigFloat[-7.214150399279151977978320029665100466133e+09]
@show(hyp)
# This takes a LONG TIME to evaluate, with a massive 1.02 G allocations to allocate 37.016 GiB memory!


text:
"""
# module SimpleTAGMLTokenizer

- Julia version:
- Author: bramb
- Date: 2019-02-12

# Examples

```jldoctest
julia>using SimpleTAGMLTokenizer

julia>tokens = tokenize("[r>tagml example<r]")

```
"""
module SimpleTAGMLTokenizer

export tokenize,
    Token,
    TokenType,OPENTAG,CLOSETAG,TEXT,ERROR

@enum(TokenType,
    OPENTAG,
    CLOSETAG,
    TEXT,
    ERROR
)

struct Token
    content::String
    type::TokenType
end

function tokenize(tagml::String)
    tokens = Token[]
    token_buf = IOBuffer()
    token_type = ERROR
    for char in tagml
        if (char == '[' || char == '<')
            _push_token!(tokens, token_buf, token_type)
            token_type = (char == '[') ? OPENTAG : CLOSETAG
            print(token_buf,char)
        elseif (char == '>' || char == ']')
            print(token_buf,char)
            _push_token!(tokens, token_buf, token_type)
            token_type = TEXT
        else
            print(token_buf,char)
        end
    end
    token_content = String(take!(token_buf))
    if (token_content != "")
        token = Token(token_content,ERROR)
        push!(tokens,token)
    end
    return tokens
end

function _push_token!(tokens, token_buf, token_type)
    token_content = String(take!(token_buf))
    if (token_content != "")
        token = Token(token_content,token_type)
        push!(tokens,token)
    end
    token_buf = IOBuffer()
end

end

max_stars_repo_name: brambg/HyperGraphTools

max_issues_repo_name: rhdekker/HyperGraphTools

max_forks_repo_name: rhdekker/HyperGraphTools


text:
using JuliaWebAPIPlugins
using Base.Test
using Logging

mytest(v1::AbstractString, v2::AbstractString) = mytest(parse(Int,v1), parse(Int,v2))
mytest(v1::Int, v2::Int) = v1 * v2

function mytestexit()
    info("exiting mytest")
end

function test_apiclntsrvr()
    configure(Dict(:debug_level => Logging.DEBUG))
    srvr = APIServer(mytest, mytestexit, "mytestserver")
    @async run(srvr)

    result = srvrcall("mytestserver", "mytest", 10, 20)
    @test result == 200
end

!isempty(ARGS) && (ARGS[1] == "--run") && test_apiclntsrvr()

max_stars_repo_name: tanmaykm/JuliaWebAPIPlugins

max_issues_repo_name: tanmaykm/JuliaWebAPIPlugins

max_forks_repo_name: tanmaykm/JuliaWebAPIPlugins


text:
module Gaius

import LinearAlgebra
import LoopVectorization
import StructArrays
import UnsafeArrays
import VectorizationBase

import VectorizationBase: stridedpointer

using LinearAlgebra: Adjoint, Transpose
using LoopVectorization: @avx
using StructArrays: StructArray
using UnsafeArrays: @uviews, UnsafeArray
using VectorizationBase: AbstractStridedPointer, gesp, vload, vstore!

export t_blocked_mul
export t_blocked_mul!

include("global_constants.jl")
include("types.jl")

include("block_operations.jl")
include("check_compatible_sizes.jl")
include("choose_block_size.jl")
include("kernels.jl")
include("matmul.jl")
include("public_mul.jl")

include("init.jl") # `Gaius.__init__()` is defined in this file

end # module Gaius


text:
#=
  Original Author: Anastasia Yendiki

  Copyright © 2022 The General Hospital Corporation (Boston, MA) "MGH"

  Terms and conditions for use, reproduction, distribution and contribution
  are found in the 'FreeSurfer Software License Agreement' contained
  in the file 'LICENSE' found in the FreeSurfer distribution, and here:

  https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense

  Reporting: freesurfer@nmr.mgh.harvard.edu
=#

using LinearAlgebra, Statistics, DelimitedFiles

export RUMBASD, rumba_rec, rumba_peaks!, rumba_write

"Container for outputs of a RUMBA-SD fit"
struct RUMBASD
  fodf::MRI
  fgm::MRI
  fcsf::MRI
  peak::Vector{MRI}
  gfa::MRI
  var::MRI
  snr_mean::Float32
  snr_std::Float32
end

"""
    RUMBAwork{T}

Pre-allocated workspace for RUMBA-SD reconstruction computations

- `T::DataType`         : Data type for computations (default: `Float32`)
- `nmask::Int`          : Number of voxels in brain mask
- `ndir::Int`           : Number of unique diffusion-encoding directions
- `ncomp::Int`          : Number of components in fODF signal decomposition
- `nvox::NTuple{3,Int}` : Dimensions of image volume (or nothing if TV not used)
- `isort::Vector{Vector{Int}}`   : Fiber ODF peak vertex indices (sorted)
- `fodf_peak::Vector{Vector{T}}` : Fiber ODF amplitudes at local peaks
- `dodf_mat::Matrix{T}` : Work matrices in diffusion ODF space
- `dodf_sig_mat::Matrix{T}`
- `Iratio::Matrix{T}`
- `fodf_mat::Matrix{T}` : Work matrices in fiber ODF space
- `rl_mat::Matrix{T}`
- `rl_mat2::Matrix{T}`
- `tv_mat::Matrix{T}`
- `tv_vol::Vector{Array{T,3}}` : Work volumes in image space
- `Gx_vol::Vector{Array{T,3}}`
- `Gy_vol::Vector{Array{T,3}}`
- `Gz_vol::Vector{Array{T,3}}`
- `Div_vol::Vector{Array{T,3}}`
- `λ::Array{T,3}`
- `σ2_vec::Matrix{T}`   : Noise variance and SNR by voxel
- `snr_vec::Matrix{T}`
"""
struct RUMBAwork{T}
  nmask::Int
  ndir::Int
  ncomp::Int
  nvox::NTuple{3,Int}
  isort::Vector{Vector{Int}}
  fodf_peak::Vector{Vector{T}}
  dodf_mat::Matrix{T}
  dodf_sig_mat::Matrix{T}
  Iratio::Matrix{T}
  fodf_mat::Matrix{T}
  rl_mat::Matrix{T}
  rl_mat2::Matrix{T}
  tv_mat::Matrix{T}
  tv_vol::Vector{Array{T,3}}
  Gx_vol::Vector{Array{T,3}}
  Gy_vol::Vector{Array{T,3}}
  Gz_vol::Vector{Array{T,3}}
  Div_vol::Vector{Array{T,3}}
  λ::Array{T,3}
  σ2_vec::Matrix{T}
  snr_vec::Matrix{T}

  function RUMBAwork(nmask::Int, ndir::Int, ncomp::Int,
                     nvox::Union{NTuple{3,Int}, Nothing}=nothing,
                     T::DataType=Float32)

    # Work vectors for fiber ODF peak search
    isort        = [Vector{Int}(undef, ncomp-2) for tid in 1:Threads.nthreads()]
    fodf_peak    = [Vector{T}(undef, ncomp-2)   for tid in 1:Threads.nthreads()]

    # Work matrices in diffusion ODF space
    dodf_mat     = Matrix{T}(undef, ndir, nmask)
    dodf_sig_mat = Matrix{T}(undef, ndir, nmask)
    Iratio       = Matrix{T}(undef, ndir, nmask)

    # Work matrices in fiber ODF space
    fodf_mat     = Matrix{T}(undef, ncomp, nmask)
    rl_mat       = Matrix{T}(undef, ncomp, nmask)
    rl_mat2      = Matrix{T}(undef, ncomp, nmask)
    tv_mat       = Matrix{T}(undef, ncomp, nmask)

    # Work volumes in image space
    isnothing(nvox) && (nvox = (0,0,0))
    tv_vol       = [Array{T,3}(undef, nvox) for tid in 1:Threads.nthreads()]
    Gx_vol       = [Array{T,3}(undef, nvox) for tid in 1:Threads.nthreads()]
    Gy_vol       = [Array{T,3}(undef, nvox) for tid in 1:Threads.nthreads()]
    Gz_vol       = [Array{T,3}(undef, nvox) for tid in 1:Threads.nthreads()]
    Div_vol      = [Array{T,3}(undef, nvox) for tid in 1:Threads.nthreads()]
    λ            = Array{T,3}(undef, nvox)

    # Noise variance and SNR by voxel
    σ2_vec       = Matrix{T}(undef, 1, nmask)
    snr_vec      = Matrix{T}(undef, 1, nmask)

    new{T}(
      nmask::Int,
      ndir::Int,
      ncomp::Int,
      nvox::NTuple{3,Int},
      isort::Vector{Vector{Int}},
      fodf_peak::Vector{Vector{T}},
      dodf_mat::Matrix{T},
      dodf_sig_mat::Matrix{T},
      Iratio::Matrix{T},
      fodf_mat::Matrix{T},
      rl_mat::Matrix{T},
      rl_mat2::Matrix{T},
      tv_mat::Matrix{T},
      tv_vol::Vector{Array{T,3}},
      Gx_vol::Vector{Array{T,3}},
      Gy_vol::Vector{Array{T,3}},
      Gz_vol::Vector{Array{T,3}},
      Div_vol::Vector{Array{T,3}},
      λ::Array{T,3},
      σ2_vec::Matrix{T},
      snr_vec::Matrix{T}
    )
  end
end

"""
    tensor_model(φ::Number, θ::Number, λ::Vector, b::Vector, g::Matrix, s0::Number)

    Compute the expected diffusion-weighted signal in a voxel, assuming that
    diffusion can be modeled by a tensor with orientation angles `φ`, `θ` and
    eigenvalues `λ` and that the signal was acquired with b-values `b` and
    gradient vectors `g`, and that the non-diffusion-weighted signal is `s0`.
"""
function tensor_model(φ::Number, θ::Number, λ::Vector, b::Vector, g::Matrix, s0::Number)

  if length(λ) != 3
    error("Length of diffusivity vector " * string(λ) * " must be 3")
  end

  R = ang2rot(φ, θ)
  D = R * Diagonal(λ) * R'

  S = s0 .* exp.(-b .* diag(g*D*g'))

  return S
end

"""
    besseli_ratio(nu::Integer, z::Float32)

Evaluate the ratio of the modified Bessel functions of the first kind,
of orders `nu` and `nu`-1, at the points in `z`.

Instead of direct evaluation, i.e., besseli(`nu`,`z`) / besseli(`nu`-1,`z`),
compute with Perron's continued fraction equation, which is an order of
magnitude faster.

Perron's continued fraction is substantially superior than Gauss' continued
fraction when z >> nu, and only moderately inferior otherwise
(Walter Gautschi and Josef Slavik, Math. Comp., 32(143):865-875, 1978).
"""
function besseli_ratio(nu::Integer, z::T) where T <: AbstractFloat

  return z / ((2*nu + z) -
               ((2*nu+1)*z /
               (2*z + (2*nu+1) -
               ((2*nu+3)*z /
               ((2*nu+2) + 2*z - ((2*nu+5) * z / ((2*nu+3) + 2*z)))))))
end

"""
    Gradient operator
"""
function sd_grad!(Gx_vol::Array{T, 3}, Gy_vol::Array{T, 3}, Gz_vol::Array{T, 3}, fODF_vol::Array{T, 3}) where T <: AbstractFloat

  @views Gx_vol .= fODF_vol[[2:end; end], :, :] .- fODF_vol
  @views Gy_vol .= fODF_vol[:, [2:end; end], :] .- fODF_vol
  @views Gz_vol .= fODF_vol[:, :, [2:end; end]] .- fODF_vol
end

"""
    Divergence operator
"""
function sd_div!(Div_vol::Array{T, 3}, Gx_vol::Array{T, 3}, Gy_vol::Array{T, 3}, Gz_vol::Array{T, 3}) where T <: AbstractFloat

  @views Div_vol[2:end-1,:,:]  .=  Gx_vol[2:end-1,:,:] .- Gx_vol[1:end-2, :, :]
  @views Div_vol[1,:,:]        .=  Gx_vol[1,:,:]	# boundaries
  @views Div_vol[end,:,:]      .= -Gx_vol[end-1,:,:]

  @views Div_vol[:,2:end-1,:] .+=  Gy_vol[:,2:end-1,:] .- Gy_vol[:, 1:end-2, :]
  @views Div_vol[:,1,:]       .+=  Gy_vol[:,1,:]	# boundaries
  @views Div_vol[:,end,:]     .+= -Gy_vol[:,end-1,:]

  @views Div_vol[:,:,2:end-1] .+=  Gz_vol[:,:,2:end-1] .- Gz_vol[:, :, 1:end-2]
  @views Div_vol[:,:,1]       .+=  Gz_vol[:,:,1]	# boundaries
  @views Div_vol[:,:,end]     .+= -Gz_vol[:,:,end-1]
end

"""
    Total variation (TV) regularization term

Read the fODF amplitudes for a single vertex from a 3D volume in the workspace
and compute the TV regularization term in place, overwriting the volume
"""
function rumba_tv!(W::RUMBAwork{T}) where T <: AbstractFloat

  tid = Threads.threadid()

  # Compute spatial gradients
  sd_grad!(W.Gx_vol[tid], W.Gy_vol[tid], W.Gz_vol[tid], W.tv_vol[tid])

  # Normalize spatial gradients
  W.tv_vol[tid] .= sqrt.(W.Gx_vol[tid].^2 .+ W.Gy_vol[tid].^2 .+
                                             W.Gz_vol[tid].^2 .+ eps(T))
  W.Gx_vol[tid] ./= W.tv_vol[tid]
  W.Gy_vol[tid] ./= W.tv_vol[tid]
  W.Gz_vol[tid] ./= W.tv_vol[tid]

  # Compute divergence
  sd_div!(W.Div_vol[tid], W.Gx_vol[tid], W.Gy_vol[tid], W.Gz_vol[tid])

  # Compute TV term (abs/eps to ensure values > 0)
  W.tv_vol[tid] .= 1 ./ (abs.(1 .- W.λ .* W.Div_vol[tid]) .+ eps(T))
end

"""
Initialize RUMBA-SD estimates
"""
function rumba_sd_initialize!(W::RUMBAwork{T}, fodf_0::Vector{T}, kernel::Matrix{T}, signal_mat::AbstractArray{T}, λ::T) where T <: AbstractFloat

  fill!(W.fodf_mat, 1)
  W.fodf_mat .*= fodf_0

  fill!(W.dodf_mat, 1)
  W.dodf_mat .*= (kernel * fodf_0)

  fill!(W.λ, λ)
  fill!(W.σ2_vec, λ)

  W.dodf_sig_mat .= (signal_mat .* W.dodf_mat) ./ W.σ2_vec

  fill!(W.tv_mat, 1)
end

"""
    fodf:     Volume fractions of anisotropic compartments
    fgm:      Volume fraction of GM isotropic component
    fcsf:     Volume fraction of CSF isotropic component
    var:      Noise variance
    snr_mean: Estimated mean snr
    snr_std:  Estimated SNR standard deviation
"""
function rumba_sd_iterate!(W::RUMBAwork{T}, signal_mat::AbstractArray{T}, kernel::Array, ind_mask::Vector{Int}, iter::Int, n_order::Integer, coil_combine::String, ipat_factor::Integer, use_tv::Bool) where T <: AbstractFloat

  fzero = T(0)
  ε = eps(T)

  ndir, ncomp = size(kernel)

  # -------------------- R-L deconvolution term ------------------------ #
  # Ratio of modified Bessel functions of order n_order and n_order-1
  W.Iratio .= besseli_ratio.(n_order, W.dodf_sig_mat)

  mul!(W.rl_mat,  kernel', signal_mat .* W.Iratio)
  mul!(W.rl_mat2, kernel', W.dodf_mat)
  W.rl_mat ./= (W.rl_mat2 .+= ε)

  # -------------------- TV regularization term ------------------------ #
  @time if use_tv
    Threads.@threads for icomp in 1:ncomp
      tid = Threads.threadid()

      # Embed fODF amplitudes in brain mask
      fill!(W.tv_vol[tid], 0)
      W.tv_vol[tid][ind_mask] = view(W.fodf_mat, icomp, :)

      # Compute TV term in place
      rumba_tv!(W)

      # Extract TV term from brain mask
      W.tv_mat[icomp, :] = view(W.tv_vol[tid], ind_mask)
    end
  end

  # ------------------------- Update estimate -------------------------- #
  # Enforce positivity
  W.fodf_mat .= max.(W.fodf_mat .* W.rl_mat .* W.tv_mat, fzero)

  if iter <= 100		&& false
    # Energy preservation at each step that included the bias on the s0 image
    # Only used to stabilize recovery in early iterations
    cte = sqrt(1 + 2*n_order*mean(W.σ2_vec))
    cte = sqrt(1 + n_order*mean(W.σ2_vec))
    cte = 1
    W.fodf_mat .= cte .* W.fodf_mat ./ (sum(W.fodf_mat, dims=1) .+ ε)
  end

  mul!(W.dodf_mat, kernel, W.fodf_mat)
  W.dodf_sig_mat .= (signal_mat .* W.dodf_mat) ./ W.σ2_vec

  # --------------------- Noise variance estimate ---------------------- #
  W.Iratio .= (signal_mat.^2 + W.dodf_mat.^2) ./ 2 .-
              (W.σ2_vec .* W.dodf_sig_mat) .* W.Iratio
  W.σ2_vec .= sum(W.Iratio , dims=1) ./ (n_order * ndir)

  # Assume that estimate of σ is in interval [1/snr_min, 1/snr_max],
  # where snr_min = 8 and snr_max = 80
  clamp!(W.σ2_vec, T((1/80)^2), T((1/8)^2))

  W.snr_vec .= 1 ./ sqrt.(W.σ2_vec)

  # ----------------- Update regularization parameter λ ----------------- #
  if use_tv
    if ipat_factor == 1
      # Penalize all voxels equally, assuming equal variance in all voxels
      # For low levels of noise, enforce a minimum level of regularization
      fill!(W.λ, max(mean(W.σ2_vec), T((1/30)^2)))

    elseif ipat_factor > 1
      # Adaptive spatial regularization, assuming spatially inhomogeneous
      # variance, e.g., tissue dependent or due to parallel imaging
      # (in the future, λ could be low-pass filtered for robust estimation)
      fill!(W.λ, 0)
      W.λ[ind_mask] = W.σ2_vec
    end
  end
end

"""
Find fODF peaks, given an fODF amplitude vector and an amplitude threshold.

Return the number of peaks found.
"""
function rumba_peaks!(W::RUMBAwork{T}, fodf::AbstractVector{T}, f_iso::T, idx_neig::Vector{Vector{Int}}, thr::T) where T <: AbstractFloat

  tid = Threads.threadid()

  # Use higher threshold in voxels with high f_iso
  #thr_xyz = thr + f_iso
  thr_xyz = thr / (1 - f_iso)

  # Find local maxima of ODF within a neighborhood around each vertex
  W.fodf_peak[tid] .= fodf

  thr_abs = thr_xyz * maximum(fodf)

  # Compare the fODF amplitude at each vertex to the amplitudes of its neighbors
  for ivert in 1:length(fodf)
    if fodf[ivert] < thr_abs ||
       fodf[ivert] .<= maximum(@view fodf[idx_neig[ivert]])
      W.fodf_peak[tid][ivert] = 0
    end
  end

  # Sort peaks by decreasing amplitude
  sortperm!(W.isort[tid], W.fodf_peak[tid], rev=true)

  return length(findall(W.fodf_peak[tid] .> 0))
end

"""
    rumba_rec(dwi::MRI{T}, mask::MRI, odf_dirs::ODF=sphere_724, niter::Integer=600, λ_para::T=T(1.7*10^-3), λ_perp::T=T(0.2*10^-3), λ_csf::T=T(3.0*10^-3), λ_gm::T=T(0.8*10^-4), ncoils::Integer=1, coil_combine::String="SMF-SENSE", ipat_factor::Integer=1, use_tv::Bool=true) where T <: AbstractFloat

Perform robust and unbiased model-based spherical deconvolution (RUMBA-SD)
reconstruction of DWIs, and return a `RUMBASD` structure.

If you use this method, please cite:
Erick J. Canales-Rodríguez, et al. (2015). Spherical deconvolution of multichannel diffusion MRI data with non-Gaussian noise models and spatial regularization. PLoS ONE, 10(10), e0138910. https://doi.org/10.1371/journal.pone.0138910

# Arguments
- `dwi::MRI{T}`: A series of DWIs, stored in an `MRI` structure with valid `.bvec`
  and `.bval` fields
- `mask::MRI`: A brain mask volume, stored in an `MRI` structure
- `odf_dirs::ODF=sphere_724`: The vertices and faces of the ODF tessellation,
  stored in an `ODF` structure
- `λ_para::T=T(1.7*10^-3)`: Axial diffusivity in white-matter
  voxels with a single fiber population
- `λ_perp::T=T(0.2*10^-3)`: Radial diffusivity in white-matter
  voxels with a single fiber population
- `λ_csf::T=T(3.0*10^-3)`: Mean diffusivity in CSF voxels
- `λ_gm::T=T(0.8*10^-4)`: Mean diffusivity in gray-matter voxels
- `ncoils::Integer=1`: Number of receive coil elements, if the DWIs were
  collected with parallel imaging, or 1 otherwise
- `coil_combine::String="SMF-SENSE"`: Method that was used to combine signals
  from receive coil elements (possible values: "SMF-SENSE" or "SoS-GRAPPA";
  has no effect if `ncoils` equals 1)
- `ipat_factor::Integer=1`: Acceleration factor, if the DWIs were collected
  with parallel imaging, or 1 otherwise
- `use_tv::Bool=true`: If true, include a total-variation regularization term
  in the FOD reconstruction

# Output
In the `RUMBASD` structure:
- `.fodf`: Volume fractions of anisotropic compartments (one per vertex)
- `.fgm`: Volume fraction of GM isotropic compartment
- `.fcsf`: Volume fraction of CSF isotropic compartment
- `.peak`: Orientation vectors of the 5 peak ODF amplitudes
- `.gfa`: Generalized fractional anisotropy
- `.var`: Estimated noise variance
- `.snr_mean`: Estimated mean SNR
- `.snr_std`: Estimated SNR standard deviation

"""
function rumba_rec(dwi::MRI{Array{T,4}}, mask::MRI, odf_dirs::ODF=sphere_724, niter::Integer=600, λ_para::T=T(1.7*10^-3), λ_perp::T=T(0.2*10^-3), λ_csf::T=T(3.0*10^-3), λ_gm::T=T(0.8*10^-4), ncoils::Integer=1, coil_combine::String="SMF-SENSE", ipat_factor::Integer=1, use_tv::Bool=true) where T <: AbstractFloat

  if isempty(dwi.bval)
    error("Missing b-value table from input DWI structure")
  end

  if isempty(dwi.bvec)
    error("Missing gradient table from input DWI structure")
  end

  n_order = 1

  if coil_combine == "SoS-GRAPPA"
    n_order = ncoils
  elseif coil_combine != "SMF-SENSE"
    error("Unknown coil combine mode " * coil_combine)
  end

  if ipat_factor < 1
    error("iPAT factor must be a positive integer")
  end

  nx, ny, nz = size(mask.vol)
  nxyz = nx * ny * nz

  ind_mask = findall(vec(mask.vol) .> 0)
  nmask = length(ind_mask)

  ib0 = (dwi.bval .== minimum(dwi.bval))
  ndir = length(findall(.!ib0)) + 1

  # Extract diffusion signal from brain mask
  # Rearrange signals so that average low-b volume is first, followed by DWIs
  signal_mat = Array{T, 2}(undef, ndir, nmask)

  @views signal_mat[1, :] .=
         mean(max.(dwi.vol[:, :, :, ib0], 0), dims=4)[ind_mask]
  @views signal_mat[2:end, :] .=
         reshape(max.(dwi.vol[:, :, :, .!ib0], 0), nxyz, ndir-1)[ind_mask, :]'

  # Divide DWIs by average low-b volume
  @views signal_mat[2:end, :] ./= signal_mat[1, :]'
  @views signal_mat[isnan.(signal_mat)] .= 0

  # signal=1 if b=0
  @views signal_mat[1, :] .= T.(signal_mat[1, :] .> 0)
  @views signal_mat[signal_mat .> 1] .= 1

  # Normalize gradient vectors
  g = vcat([0 0 0], dwi.bvec[.!ib0, :] ./
                    sqrt.(sum(dwi.bvec[.!ib0, :].^2, dims=2)))
  b = vec(vcat([0], dwi.bval[.!ib0, :]))

  # Vertices and faces on the sphere where the fODF will be sampled
  nvert = size(odf_dirs.vertices, 1)
  nvert = div(nvert,2)			# Will work on half sphere
  half_vertices  = odf_dirs.vertices[1:nvert, :]

  # Find all vertices within an angular neigborhood
  if odf_dirs == sphere_724 || odf_dirs == sphere_642
    ang_neig = 12.5
  elseif odf_dirs == sphere_362
    ang_neig = 16
  end

  cos_ang = half_vertices * half_vertices'
  cos_ang[cos_ang .>  1] .=  1
  cos_ang[cos_ang .< -1] .= -1

  ang = acosd.(cos_ang)
  ang = min.(ang, 180 .- ang)

  isneig = (ang .< ang_neig)
  isneig[diagind(isneig)] .= 0
  idx_neig = findall.(eachslice(isneig, dims=1))

  #
  # Generate reconstruction kernel based on multi-tensor model
  #
  Kernel = Matrix{T}(undef, length(b), nvert+2)

  s0 = 1		# Assuming unit (normalized) signal

  # One radially symmetric tensor compartment for each fODF vertex
  (φ, θ) = cart2sph(odf_dirs.vertices[nvert+1:end,1],
                    odf_dirs.vertices[nvert+1:end,2],
                    odf_dirs.vertices[nvert+1:end,3])
  θ .= -θ

  for ivert in 1:nvert
    Kernel[:, ivert] .= tensor_model(φ[ivert], θ[ivert],
                                     [λ_para, λ_perp, λ_perp], b, g, s0)
  end

  # One isotropic tensor compartment for CSF
  Kernel[:, nvert+1] .= tensor_model(0, 0, [λ_csf, λ_csf, λ_csf], b, g, s0)

  # One isotropic tensor compartment for GM
  Kernel[:, nvert+2] .= tensor_model(0, 0, [λ_gm, λ_gm, λ_gm], b, g, s0)

  # Total number of compartments
  ncomp = nvert+2

  #
  # Initialize fODFs (uniform distribution on the half sphere)
  #
  fodf_0 = ones(T, nvert+2)
  	fodf_0 .= fodf_0 ./ (2*nvert+2)
  fodf_0 .= fodf_0 ./ sum(fodf_0)

  # Allocate workspace
  W = RUMBAwork(nmask, ndir, ncomp, use_tv ? (nx,ny,nz) : (0,0,0))

  # Initialize estimates
  σ0 = T(1/15)
  λ0 = σ0^2
  snr_mean = snr_std = T(0)

  rumba_sd_initialize!(W, fodf_0, Kernel, signal_mat, λ0)

  #
  # Reconstruct fODFs
  #
  @time for iter in 1:niter
    println("Iteration " * string(iter) * " of " * string(niter))

    @time rumba_sd_iterate!(W, signal_mat, Kernel, ind_mask, iter,
                      n_order, coil_combine, ipat_factor, use_tv)

    snr_mean = mean(W.snr_vec)
    snr_std  =  std(W.snr_vec; mean=snr_mean)

    println("Estimated mean SNR (s0/σ) = " * string(snr_mean) *
            " (+-) " * string(snr_std))

    println("Number of coils = " * string(n_order))

    println("Mean sum(fODF) = " * string(mean(sum(W.fodf_mat, dims=1))))
  end

  # Energy preservation
  W.fodf_mat ./= (sum(W.fodf_mat, dims=1) .+ eps(T))

  #
  # Embed estimates in brain mask
  #
  fodf = MRI(mask, nvert, T)
  fcsf = MRI(mask, 1, T)
  fgm  = MRI(mask, 1, T)
  var  = MRI(mask, 1, T)
  gfa  = MRI(mask, 1, T)

  # Volume fractions of anisotropic WM compartments
  for icomp in 1:ncomp-2
    @views fodf.vol[ind_mask .+ (icomp-1)*nxyz] .= W.fodf_mat[icomp, :]
  end

  # Volume fraction of isotropic CSF compartment
  fcsf.vol[ind_mask] = W.fodf_mat[end-1, :]

  # Volume fraction of isotropic GM compartment
  fgm.vol[ind_mask] = W.fodf_mat[end, :]

  f_iso = fgm.vol + fcsf.vol

  # Add isotropic components to ODF and normalize to sum=1
  fodf.vol .= fodf.vol .+ f_iso

  fodf.vol .= fodf.vol ./ sum(fodf.vol, dims=4)
  fodf.vol[isnan.(fodf.vol)] .= 0;

  # Noise variance
  var.vol[ind_mask] = W.σ2_vec

  #
  # Compute GFA
  #
  gfa.vol .= std(fodf.vol, dims=4) ./ sqrt.(mean(fodf.vol.^2, dims=4))
  gfa.vol[isnan.(gfa.vol)] .= 0

  #
  # Extract ODF peaks
  #
  npeak = 5			# Max number of peaks per voxel
  fthresh = T(.1)		# Min volume fraction to retain a peak

  peak = Vector{MRI}(undef, npeak);

  for ipeak in 1:npeak
    peak[ipeak] = MRI(mask, 3, T)
  end

  Threads.@threads for iz in 1:nz
    for iy in 1:ny
      for ix in 1:nx
        (mask.vol[ix, iy, iz] == 0) && continue

        tid = Threads.threadid()

        nvalid = rumba_peaks!(W, fodf.vol[ix, iy, iz, :], f_iso[ix, iy, iz],
                              idx_neig, fthresh)

        n = min(nvalid, npeak)

        fnorm = (1 - f_iso[ix, iy, iz]) /
                sum(fodf.vol[ix, iy, iz, W.isort[tid][1:n]])

        for ipeak in 1:n
          @views peak[ipeak].vol[ix, iy, iz, :] .=
                   half_vertices[W.isort[tid][ipeak], :] .*
                   (fodf.vol[ix, iy, iz, W.isort[tid][ipeak]] * fnorm)
        end
      end
    end
  end

  return RUMBASD(fodf, fgm, fcsf, peak, gfa, var, snr_mean, snr_std)
end

"""
    rumba_write(rumba::RUMBASD, basename::String)

Write the volumes from a `RUMBASD` structure that was created by `rumba_rec()`
to files whose names start with the specified base name.
"""
function rumba_write(rumba::RUMBASD, basename::String)

  for var in fieldnames(RUMBASD)
    vartype = fieldtype(RUMBASD, var)

    if vartype == MRI
      fname = basename * "_" * string(var) * ".nii.gz"
      mri_write(getfield(rumba, var), fname)
    elseif vartype == Vector{MRI}
      for ivol in 1:length(getfield(rumba, var))
        fname = basename * "_" * string(var) * string(ivol) * ".nii.gz"
        mri_write(getfield(rumba, var)[ivol], fname)
      end
    else
      fname = basename * "_" * string(var) * ".txt"
      writedlm(fname, getfield(rumba, var), ' ')
    end
  end
end

max_stars_repo_name: freesurfer/jupysurf

max_issues_repo_name: freesurfer/jupysurf

max_forks_repo_name: freesurfer/jupysurf


text:
import GitLab
using Base.Test

myauth = GitLab.authenticate(ENV["GITLAB_AUTH"]) # don't hardcode your access tokens!
println("Authentication successful")
options = Dict("private_token" => myauth.token)

myrepo = GitLab.repo_by_name("TestProject1"; headers=options)

merge_requests, page_data = GitLab.pull_requests(myrepo; params=options)
## @show merge_requests

merge_request = GitLab.pull_request(myrepo, get(merge_requests[1].id); params=options)
## @show merge_request
@test get(merge_request.title) == "test"

println("Merge Request Tests Done !!!")


text:
"""
    compute_E_M_Ising(σ::Array{Int8, 2})

Computes the Ising model internal energy and magnetization for the input spin configuration `σ`.

Arguments
* `σ`     : A 2D array representing a square lattice in the Ising model
Returns
* `(E, M)`: The system's energy and magnetization
"""
function compute_E_M_Ising(σ::Array{Int8,2})
    E::Int = 0
    M::Int = 0
    L = size(σ, 1)
    for j = 1:L, i = 1:L
        E -= σ[i, j] * (σ[plus(i, L), j] + σ[i, plus(j, L)])
        M += σ[i, j]
    end

    return (E, M)
end

"""
    compute_E_Potts(σ::Array{Int8, 2})

Computes the Potts model internal energy for the input spin configuration `σ`.

Arguments
* `σ` : A 2D array representing a square lattice in the Potts model
Returns
* `E` : The system's energy
"""
function compute_E_Potts(σ::Array{Int8,2})
    E::Int = 0
    L = size(σ, 1)
    for j = 1:L, i = 1:L
        E -= Int(σ[i, j] == σ[plus(i, L), j]) + Int(σ[i, j] == σ[i, plus(j, L)])
    end

    return E
end

"""
    compute_E_M_XY(σ::Array{NTuple{2, Float64}, 2})

Computes the XY model internal energy and magnetization for the input spin configuration `σ`.

Arguments
* `σ`          : A 2D array representing a square lattice in the Ising model
Returns
* `(E, Mx, My)`: The system's energy, x-magnetization, and y-magnetization
"""
function compute_E_M_XY(σ::Array{NTuple{2,Float64},2})
    E = 0.0
    Mx = 0.0
    My = 0.0
    L = size(σ, 1)
    for j = 1:L, i = 1:L
        E -= dot(σ[i, j], σ[plus(i, L), j]) + dot(σ[i, j], σ[i, plus(j, L)])
        Mx += σ[i, j][1]
        My += σ[i, j][2]
    end

    return (E, Mx, My)
end

"""
    update_observables!(model::Ising)

Pushes the current energy and magnetization to their respective arrays.

Arguments
* `model`: Ising type
"""
function update_observables!(model::Ising)
    E, M = compute_E_M_Ising(model.σ)
    push!(model.observables.E, E)
    push!(model.observables.M, M)
end

"""
    update_observables!(model::Potts)

Pushes the current energy and magnetization to their respective arrays.

Arguments
* `model`: Potts type
"""
function update_observables!(model::Potts)
    E = compute_E_Potts(model.σ)
    push!(model.observables.E, E)
end

"""
    update_observables!(model::XY)

Pushes the current energy and magnetization to their respective arrays.

Arguments
* `model`: XY type
"""
function update_observables!(model::XY)
    E, Mx, My = compute_E_M_XY(model.σ)
    push!(model.observables.E, E)
    push!(model.observables.Mx, Mx)
    push!(model.observables.My, My)
end

"""
    compute_observables_statistics!(model::AbstractMagnetismModel)

Computes the average energy per spin, average magnetization per spin,
heat capacity per spin, magnetic susceptibility per spin, Binder parameter,
and the integrated autocorrelation times (computed using the binning method).
"""
function compute_observables_statistics!(model::AbstractMagnetismModel)
    V = model.params.L^2
    model.observables.statistics[:e] = mean(model.observables.E) / V
    model.observables.statistics[:e²] = mean(model.observables.E .^ 2) / V^2
    model.observables.statistics[:c] = model.params.β^2 * var(model.observables.E) / V

    if isa(model, Ising)
        model.observables.statistics[:m] = mean(abs.(model.observables.M)) / V
        model.observables.statistics[:m²] = mean(model.observables.M .^ 2) / V^2
        model.observables.statistics[:χ] =
            model.params.β * var(abs.(model.observables.M)) / V
        model.observables.statistics[:U] =
            1 - mean(model.observables.M .^ 4) / 3 / mean(model.observables.M .^ 2)^2
    elseif isa(model, XY)
        model.observables.statistics[:mx] = mean(model.observables.Mx) / V
        model.observables.statistics[:mx²] = mean(model.observables.Mx .^ 2) / V^2
        model.observables.statistics[:my] = mean(model.observables.My) / V
        model.observables.statistics[:my²] = mean(model.observables.My .^ 2) / V^2
        model.observables.statistics[:m²] =
            mean(model.observables.Mx .^ 2 .+ model.observables.My .^ 2) / V^2
        model.observables.statistics[:χ] = model.observables.statistics[:m²] * V
    end

    model.observables.statistics[:τ_int] =
        compute_autocorrelation_time(model, collect(1:400))
end


text:
export dcdp_tmpcntrfact, vwEUP, prDrill_infill!, serial_counterfact_all!

struct dcdp_tmpcntrfact{T<:Real}
    rin::Array{T,3}
    rex::Array{T,3}
    ubUfull::Array{T,3}
    q::Array{T,3}
end

function dcdp_tmpcntrfact(p::dcdp_primitives{FF,T}) where {FF,T}
    nz = _nz(p)
    nψ = _nψ(p)
    nd = _nd(p)
    rin = Array{T,3}(nz,nψ,nd)
    rex = similar(rin)
    ubUfull = similar(rin)
    q = similar(rin)
    return dcdp_tmpcntrfact{T}(rin,rex,ubUfull,q)
end

function vwEUP(sEU::AbstractArray{T,N}, sP::AbstractArray{T,N2}, typidx::Integer...) where {T,N,N2}
    ntyps = length(typidx)

    dimz = N - ntyps - 2

    szP = size(sP)
    zdims = szP[1:dimz]
    nψ, dmaxp1, nS = szP[dimz+(1:3)]
    typdims = szP[dimz+4:end]
    (zdims..., nψ, nS, typdims...) == size(sEU) || throw(DimensionMismatch())

    colons(n) = ntuple((x)-> Colon(), n)
    nz = prod(zdims)

    @views EU = reshape(sEU[colons(dimz+2)..., typidx...], nz, nψ, nS)
    @views P  = reshape(sP[ colons(dimz+3)..., typidx...], nz, nψ, dmaxp1, nS)

    return EU, P
end

function prDrill_infill!(
    EV::AbstractArray3, uin::AbstractArray4, uex::AbstractArray3, ubVfull::AbstractArray3,
    EU::AbstractArray3, rin::AbstractArray3, rex::AbstractArray3, ubUfull::AbstractArray3,
    P::AbstractArray4, q::AbstractArray3,

    lse::AbstractMatrix, tmp::AbstractMatrix, IminusTEVp::AbstractMatrix,

    wp::AbstractUnitProblem, Πz::AbstractMatrix,

    Πψtmp::AbstractMatrix, ψspace::AbstractVector, σ::Real, β::Real
    )

    nz,nψ,nS = size(EV)
    dmaxp1 = _dmax(wp)+1

    # ------------------------ size checks ----------------------------------

    (nz,nψ,nS) == size(EU)                           || throw(DimensionMismatch())
    (nz,nψ,dmaxp1,nS) == size(P)                     || throw(DimensionMismatch())
    (nz,nψ,dmaxp1,2) == size(uin)                    || throw(DimensionMismatch())
    (nz,nz) == size(IminusTEVp) == size(Πz)          || throw(DimensionMismatch())
    (nz,nψ) == size(tmp) == size(lse)                || throw(DimensionMismatch())
    (nz,nψ,dmaxp1) == size(rin) == size(rex) == size(q) || throw(DimensionMismatch())
    (nz,nψ,dmaxp1) == size(ubVfull) == size(ubUfull) || throw(DimensionMismatch())
    nψ == size(Πψtmp,1) == size(Πψtmp,2) == length(ψspace) || throw(DimensionMismatch())

    # ------------------------ compute things ----------------------------------

    # initialize P
    P .= 0.0

    # terminal conditions
    EV[:,:,end] .= 0.0
    EU[:,:,end] .= 0.0

    # ----------- Infill drilling ----------------
    for i in ind_inf(wp)
        idxd, idxs, horzn, s = wp_info(wp, i)

        @views EV0  = EV[:,:,i]
        @views EU0  = EU[:,:,i]
        @views ubV  = ubVfull[:,:,idxd]
        @views ubU  = ubUfull[:,:,idxd]
        @views ubV .= uin[:,:,idxd,1+s.d1] .+ β .* EV[:,:,idxs]
        @views ubU .= rin[:,:,idxd] .+ β .* EU[:,:,idxs]

        # Make probabilities
        softmax3!(ubV, lse, tmp)
        P[:,:,idxd,i] .= ubV

        # Form expected revenue = (I⊗Πz)∑_d q_d .* ( u .+ β .* EV)
        sumprod!(lse, ubV, ubU)
        A_mul_B_md!(EU0, Πz, lse, 1)

        # If infinite horizon, handle this
        if horzn == :Infinite
            for j in 1:nψ
                @views update_IminusTVp!(IminusTEVp, Πz, β, ubV[:,j,1])
                fact = lufact(IminusTEVp)
                ldiv!(fact, EU0[:,j])
            end
        end
    end

    # ----------- Learning integration ----------------
    d2plus  = 2:_dmax(wp)+1  # TODO: would be good to make this more general.
    lrn2inf = inf_fm_lrn(wp)
    exp2lrn = exploratory_learning(wp)

    # integrate over U
    @views ubU1 = ubUfull[:,:,d2plus]
    @views EU1  = EU[ :,:,lrn2inf]
    _βΠψ!(Πψtmp, ψspace, σ, β)
    A_mul_B_md!(ubU1, Πψtmp, EU1, 2)
    @views ubU1 .+= rex[:,:,d2plus]

    # ubVtilde = u[:,:,2:dmaxp1] + β * Πψ ⊗ I * EV[:,:,2:dmaxp1]
    ubVfull[:,:,d2plus] .= uex[:,:,d2plus] .+ β .* EV[:,:,exp2lrn]

    # ----------- Exploratory drilling ----------------
    for i in explore_state_inds(wp)
        idxd, idxs, horzn, s = wp_info(wp, i)
        ip = action0(wp,i)

        @views EU0  = EU[:,:,i]

        @views ubV  = ubVfull[:,:,idxd]
        @views ubU  = ubUfull[:,:,idxd]
        @views qvw  = q[:,:,idxd]

        @views ubV[:,:,1] .= 0.0 + β .* EV[:,:,ip]
        @views ubU[:,:,1] .= 0.0 + β .* EU[:,:,ip]

        softmax3!(qvw, lse, tmp, ubV)
        P[:,:,idxd,i] .= qvw
        sumprod!(lse, qvw, ubU)
        A_mul_B_md!(EU0, Πz, lse, 1)
    end

end

function prDrill_infill!(evs::dcdp_Emax, tmpv::dcdp_tmpvars, prim::dcdp_primitives, EU::AbstractArray3, P::AbstractArray4, tmpc::dcdp_tmpcntrfact, σ::Real)
    return prDrill_infill!(
        evs.EV, tmpv.uin, tmpv.uex, tmpv.ubVfull,
        EU,     tmpc.rin, tmpc.rex, tmpc.ubUfull,
        P, tmpc.q,
        tmpv.lse, tmpv.tmp, tmpv.IminusTEVp,
        prim.wp, prim.Πz,
        tmpv.Πψtmp, prim.ψspace, σ, prim.β
    )
end

function prDrill_infill!(evs::dcdp_Emax, t::dcdp_tmpvars, p::dcdp_primitives{FF}, EU::AbstractArray3, P::AbstractArray4, tmpc::dcdp_tmpcntrfact, θ::AbstractVector, itype::Tuple) where {FF}
    θt, σ = _θt(θ, p), _σv(θ)
    pdct = makepdct(p, Val{:u})
    fillflows!(t, p, θt, σ, itype...)
    fillflowrevs!(FF, flowrev, tmpc.rin, tmpc.rex, θt, σ, pdct, itype...)
    prDrill_infill!(evs, t, p, EU, P, tmpc, σ)
end

function serial_counterfact_all!(sev::SharedEV, tmpv::dcdp_tmpvars, prim::dcdp_primitives, sEU::AbstractArray, sP::AbstractArray, tmpc::dcdp_tmpcntrfact, θ::AbstractVector)
    for Idx in CartesianIndices( length.(sev.itypes) )
        evs, typs = dcdp_Emax(sev, Idx.I...)
        EU,P = vwEUP(sEU, sP, Idx.I...)
        itype = getindex.(sev.itypes, Idx.I)
        prDrill_infill!(evs, tmpv, prim, EU, P, tmpc, θ, itype)
    end
end


text:
#=
    There are nine general cases for processing a*b
                   b <= 0    b <0< b    b >= 0
       a <=  0    LteLte      LteZer    LteGte
       a <0< a    ZerLte      ZerZer    ZerGte
       a >=  0    GteLte      GteZer    GteGte
=#

function (*){S<:Sculpt,W<:Sculpt,Q<:Qualia,C<:Clay}(a::Flex{S,Q,C}, b::Flex{W,Q,C})
    z = zero(C)
    if     a.hi <= z
       if     b.hi <= z
           mulLteLte(a,b)
       elseif b.lo >= z
           mulLteGte(a,b)
       else   # b straddles 0
           mulLteZer(a,b)
       end
    elseif a.lo >= z
       if     b.hi <= z
           mulGteLte(a,b)
       elseif b.lo >= z
           mulGteGte(a,b)
       else   # b straddles 0
           mulGteZer(a,b)
       end
    else   # a straddles 0
       if     b.hi <= z
           mulZerLte(a,b)
       elseif b.lo >= z
           mulZerGte(a,b)
       else   # b straddles 0
           mulZerZer(a,b)
       end
    end
end

for (fn,loa,lob,hia,hib) in [ (:mulLteLte, :(a.hi), :(b.hi), :(a.lo), :(b.lo)),
                              (:mulLteGte, :(a.lo), :(b.hi), :(a.hi), :(b.lo)),
                              (:mulLteZer, :(a.lo), :(b.hi), :(a.lo), :(b.lo)),
                              (:mulGteLte, :(a.hi), :(b.lo), :(a.lo), :(b.hi)),
                              (:mulGteGte, :(a.lo), :(b.lo), :(a.hi), :(b.hi)),
                              (:mulGteZer, :(a.hi), :(b.lo), :(a.hi), :(b.hi)),
                              (:mulZerLte, :(a.hi), :(b.lo), :(a.lo), :(b.lo)),
                              (:mulZerGte, :(a.lo), :(b.hi), :(a.hi), :(b.hi)),
                            ]
  @eval begin
    function ($fn){S<:Sculpt,W<:Sculpt,Q<:Qualia,C<:Clay}(a::Flex{S,Q,C}, b::Flex{W,Q,C})
        aLoIsOpen, aHiIsOpen = boundries(S)
        bLoIsOpen, bHiIsOpen = boundries(W)
        sculpting = boundries( (aLoIsOpen|bLoIsOpen), (aHiIsOpen|bHiIsOpen) )

        lo = (*)(($loa), ($lob), RoundDown)
        hi = (*)(($hia), ($hib), RoundUp)

        Flex{sculpting,Q,C}(lo, hi)
    end
  end
end

(*){S<:Sculpt,Q<:Qualia,C<:Clay}(a::Flex{S,Q,C}, b::C) = (*)(a, Flex{S,Q,C}(b))
(*){S<:Sculpt,Q<:Qualia,C<:Clay}(a::C, b::Flex{S,Q,C}) = (*)(Flex{S,Q,C}(a), b)

(*){S<:Sculpt,Q<:Qualia,C<:Clay}(a::Bool, b::Flex{S,Q,C}) = a ? b : Flex{S,Q,C}(0.0) # quash ambig notice
(*){S<:Sculpt,Q<:Qualia,C<:Clay,T<:Real}(a::Flex{S,Q,C}, b::T) = (*)(a, Flex{S,Q,C}(convert(C,b)))
(*){S<:Sculpt,Q<:Qualia,C<:Clay,T<:Real}(a::T, b::Flex{S,Q,C}) = (*)(Flex{S,Q,C}(convert(C,a)), b)

(*){S<:Sculpt,T<:Sculpt,Q<:Qualia,R<:Qualia,C<:Clay,D<:Clay}(a::Flex{S,Q,C}, b::Flex{T,R,D}) = (*)(promote(a,b)...)


text:
using RayTraceGR

using DifferentialEquations
using LinearAlgebra
using StaticArrays
using Test

const BigRat = Rational{BigInt}

@testset "Minkowski metric" begin
    T = BigRat
    metric = minkowski

    x = SVector{D,T}(0, 0, 0, 0)
    g = metric(x)

    detg = det(g)
    gu = inv(g)
    detgu = det(gu)

    @test detg * detgu == 1
    @test g * gu == I

    g1,dg = dmetric(metric, x)
    @test g1 == g
    @test all(==(0), dg)

    Γ = christoffel(metric, x)
    @test all(==(0), Γ)
end

@testset "Kerr-Schild metric" for i in 1:7
    T = Float32
    tol = eps(T)^(T(3)/4)
    metric = kerr_schild

    ix = i & 1
    iy = i & 2
    iz = i & 4
    x = SVector{D,T}(0, 2ix, 2iy, 2iz)

    g = metric(x)
    @test !any(isnan, g)

    detg = det(g)
    gu = inv(g)
    detgu = det(gu)

    @test abs(detg * detgu - 1) <= tol
    @test maximum(abs.(g * gu - I)) <= tol

    g1,dg = dmetric(metric, x)
    @test maximum(abs.(g - metric(x))) <= tol

    Γ = christoffel(metric, x)
    @test !any(isnan, Γ)
end

# @testset "rays" begin
#     T = Float32
#     tol = eps(T)^(T(3)/4)
#     metric = minkowski
#     x = SVector{D,T}(0, 0, 0, 0)
#     u = SVector{D,T}(-1, 1, 0, 0)
#     p = Pixel{T}(x, u, zeros(SVector{3,T}))
#     objs = Object{T}[]
#     condition(s, λ, integrator) = min_distance(objs, s)
#     affect!(integrator) = terminate!(integrator)
#     cb = ContinuousCallback(condition, affect!)
#     p = trace_ray(metric, objs, cb, p)
#     # @test maximum(abs.(p.rgb - [10, 0, 0])) <= tol
#     @test maximum(abs.(p.rgb - [1, 0, 0])) <= tol
# end


text:
@testset "bisect" begin
    @testset "basic" begin
        @system SBisect(Controller) begin
            x(x) => x - 1 ~ bisect(lower=0, upper=2)
        end
        s = instance(SBisect)
        @test s.x' == 1
    end

    @testset "unit" begin
        @system SBisectUnit(Controller) begin
            x(x) => x - u"1m" ~ bisect(lower=u"0m", upper=u"2m", u"m")
        end
        s = instance(SBisectUnit)
        @test s.x' == u"1m"
    end

    @testset "eval unit" begin
        @system SBisectEvalUnit(Controller) begin
            f(x) => (x/1u"s" - 1u"m/s") ~ track(u"m/s")
            x(f) ~ bisect(lower=0, upper=2, u"m", evalunit=u"m/s")
        end
        s = instance(SBisectEvalUnit)
        @test s.x' == 1u"m"
        @test Cropbox.unittype(s.x) == u"m"
        @test Cropbox.evalunit(s.x) == u"m/s"
    end

    @testset "equal" begin
        @system SBisectEqual(Controller) begin
            x1(x1) => (x1 - 1) ~ bisect(lower=0, upper=2)
            x2(x2) => (x2 - 1 ⩵ 0) ~ bisect(lower=0, upper=2)
        end
        s = instance(SBisectEqual)
        @test s.x1' == s.x2'
    end
end


text:
# Sample from arbitrary arrays

function sample{T <: Real}(a::AbstractArray, probs::Vector{T})
    i = rand(Categorical(probs))
    return a[i]
end

function sample(a::AbstractArray)
    n = length(a)
    probs = ones(n) ./ n
    return sample(a, probs)
end


text:
# Define all enumeration of UHD

@cenum uhd_error::UInt32 begin
    UHD_ERROR_NONE = 0
    UHD_ERROR_INVALID_DEVICE = 1
    UHD_ERROR_INDEX = 10
    UHD_ERROR_KEY = 11
    UHD_ERROR_NOT_IMPLEMENTED = 20
    UHD_ERROR_USB = 21
    UHD_ERROR_IO = 30
    UHD_ERROR_OS = 31
    UHD_ERROR_ASSERTION = 40
    UHD_ERROR_LOOKUP = 41
    UHD_ERROR_TYPE = 42
    UHD_ERROR_VALUE = 43
    UHD_ERROR_RUNTIME = 44
    UHD_ERROR_ENVIRONMENT = 45
    UHD_ERROR_SYSTEM = 46
    UHD_ERROR_EXCEPT = 47
    UHD_ERROR_BOOSTEXCEPT = 60
    UHD_ERROR_STDEXCEPT = 70
    UHD_ERROR_UNKNOWN = 100
end

@cenum uhd_rx_metadata_error_code_t::UInt32 begin
    UHD_RX_METADATA_ERROR_CODE_NONE = 0
    UHD_RX_METADATA_ERROR_CODE_TIMEOUT = 1
    UHD_RX_METADATA_ERROR_CODE_LATE_COMMAND = 2
    UHD_RX_METADATA_ERROR_CODE_BROKEN_CHAIN = 4
    UHD_RX_METADATA_ERROR_CODE_OVERFLOW = 8
    UHD_RX_METADATA_ERROR_CODE_ALIGNMENT = 12
    UHD_RX_METADATA_ERROR_CODE_BAD_PACKET = 15
end

@cenum uhd_async_metadata_event_code_t::UInt32 begin
    UHD_ASYNC_METADATA_EVENT_CODE_BURST_ACK = 1
    UHD_ASYNC_METADATA_EVENT_CODE_UNDERFLOW = 2
    UHD_ASYNC_METADATA_EVENT_CODE_SEQ_ERROR = 4
    UHD_ASYNC_METADATA_EVENT_CODE_TIME_ERROR = 8
    UHD_ASYNC_METADATA_EVENT_CODE_UNDERFLOW_IN_PACKET = 16
    UHD_ASYNC_METADATA_EVENT_CODE_SEQ_ERROR_IN_BURST = 32
    UHD_ASYNC_METADATA_EVENT_CODE_USER_PAYLOAD = 64
end

@cenum uhd_sensor_value_data_type_t::UInt32 begin
    UHD_SENSOR_VALUE_BOOLEAN = 98
    UHD_SENSOR_VALUE_INTEGER = 105
    UHD_SENSOR_VALUE_REALNUM = 114
    UHD_SENSOR_VALUE_STRING = 115
end

@cenum uhd_tune_request_policy_t::UInt32 begin
    UHD_TUNE_REQUEST_POLICY_NONE = 78
    UHD_TUNE_REQUEST_POLICY_AUTO = 65
    UHD_TUNE_REQUEST_POLICY_MANUAL = 77
end

@cenum uhd_stream_mode_t::UInt32 begin
    UHD_STREAM_MODE_START_CONTINUOUS = 97
    UHD_STREAM_MODE_STOP_CONTINUOUS = 111
    UHD_STREAM_MODE_NUM_SAMPS_AND_DONE = 100
    UHD_STREAM_MODE_NUM_SAMPS_AND_MORE = 109
end

@cenum uhd_log_severity_level_t::UInt32 begin
    UHD_LOG_LEVEL_TRACE = 0
    UHD_LOG_LEVEL_DEBUG = 1
    UHD_LOG_LEVEL_INFO = 2
    UHD_LOG_LEVEL_WARNING = 3
    UHD_LOG_LEVEL_ERROR = 4
    UHD_LOG_LEVEL_FATAL = 5
end

max_stars_repo_name: RGerzaguet/UHDBindings

max_issues_repo_name: RGerzaguet/UHDBindings


text:
function EasyStream.fit(model::MClassifier, verbosity::Int, X, Y)
    return fit(model, verbosity, X, Y)
end

function EasyStream.predict(classifier::MClassifier, fitresult, instance)
    return predict(classifier, fitresult, instance)
end

function EasyStream.updatePredict(classifier::MClassifier, fitresult, instance)
    return update_predict(classifier, fitresult, instance)
end


text:
@test isempty(lintstr("""
function f(x)
    local y::Float64 = 1
    x + y
end
"""))

s = """
function f(c::Char)
    x = Int8(c)
    x = Int16(x)
    x = Int32(x)
    x = Int64(x)
    x = Int(x)
    x = Rational(x)
    x = float(x)
    x = Complex(x)
    @lintpragma("Info type x")
    return x
end
"""
msgs = lintstr(s)
@test_broken msgs[1].code == :W545
@test_broken contains(msgs[1].message, "previously used variable has apparent type")
@test msgs[end].code == :I271
@test_broken contains(msgs[end].message, "typeof(x) == Complex")

s = """
function f()
    x = rand()
    @lintpragma("Info type x")
    return x
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(x) == Float64")

s = """
function f()
    x = rand(3)
    y = rand(Bool, 5, 5)
    @lintpragma("Info type x")
    @lintpragma("Info type y")
    return x, y
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(x) == Array{Float64,1}")
@test msgs[2].code == :I271
@test contains(msgs[2].message, "typeof(y) == Array{Bool,2}")

s = """
function f()
    d = Dict{Symbol,Int}(:a=>1, :b=>2)
    x = d[:a]
    x = 1.0
    return x
end
"""
msgs = lintstr(s)
@test_broken msgs[1].code == :W545
@test_broken contains(msgs[1].message, "previously used variable has apparent type Int64, but " *
    "now assigned Float64")

s = """
function f(arr::Array{Any,1})
    x = arr[1]::Int64
    @lintpragma("Info type x")
    return x
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(x) == Int")

s = """
g(x) = x

function f()
    x = g
    @lintpragma("Info type x")
    return x
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(x) == Function")

s = """
module MyModule
end

function f()
    x = MyModule
    @lintpragma("Info type x")
    return x
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(x) == Module")

s = """
function f()
    z = Complex(1.0, 1.0)
    z.re = 2.0
    z
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E525
@test msgs[1].variable == "z"
@test contains(msgs[1].message, "is of an immutable type Complex")

#= TODO: the warning here should be on a = Array{Int32, n}, not the E521
s = """
n = 32
a = Array{Int32, n} # bug is here
for i in 1:n
    a[i] = i
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E521
@test msgs[1].variable == "a"
@test contains(msgs[1].message, "apparent type Type")
=#

include("E522.jl")

s = """
function f()
    d = Dict(:a=>1, :b=>"")
    @lintpragma("Info type d")
    return d
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(d) == Dict")

s = """
function f()
    d = Dict{Symbol,Any}(:a=>1, :b=>"")
    @lintpragma("Info type d")
    return d
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(d) == Dict")

s = """
function f(n)
    a = Array(Float64, (1,2,3))
    @lintpragma("Info type a")
    c = Array(Float64, 1,2,3)
    @lintpragma("Info type c")
    d = zeros(Float64, (1,2))
    @lintpragma("Info type d")
    e1 = zeros(2)
    @lintpragma("Info type e1")
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(a) == Array{Float64,3}")
@test contains(msgs[2].message, "typeof(c) == Array{Float64,3}")
@test contains(msgs[3].message, "typeof(d) == Array{Float64,2}")
@test contains(msgs[4].message, "typeof(e1) == Array{Float64,1}")

s = """
function f()
    a = Array(Float64, (1,2,3))
    s = size(a)
    @lintpragma("Info type s")
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(s) == Tuple{Int64,Int64,Int64}")

s = """
function f()
    a = Complex{Float64}[]
    @lintpragma("Info type a")
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(a) == Array{Complex{Float64},1}")

s = """
    Complex(0.0,0.0) == 0
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
    Complex(1.0,0.0) > 0
"""
msgs = lintstr(s)
@test msgs[1].code == :W542
@test contains(msgs[1].message, "comparing apparently incompatible type")

s = """
s = Union(Int,Double)
"""
msgs = lintstr(s)
@test msgs[1].code == :W441

s = """
a = 1.0
a /= 2
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f(x::ANY...)
   x[1]
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function typed_hcat(A::AbstractVecOrMat...)
   A[1]
end
"""
msgs = lintstr(s)
@test isempty(msgs)

@testset "E539" begin
    # assigning error to a variable
    @test messageset(lintstr("""
    x = throw(ArgumentError("error!"))
    """)) == Set([:E539])

    @test messageset(lintstr("""
    x = 1 + "x"
    """)) == Set([:E422, :E539])

    @test isempty(lintstr("""
    x = 1 + 1 == 2 ? "OK" : error("problem")
    """))

    @test messageset(lintstr("""
    x, y = error()
    """)) == Set([:E539])

    @test messageset(lintstr("""
    κ = sqrt("x")
    """)) == Set([:E539])
end


text:
module TestEC2Utils
using AWS.EC2
using AWS
using Base.Test

function runtests(env, config)
    desc_avail_zones(env)

    ami = get(config, "ami-generic", "")
    insttype = get(config, "insttype", "")
    keyname = get(config, "keyname", "")
    ownertag = get(config, "ownertag", "")
    nametag = get(config, "nametag", "")

    for v in (ami, insttype, keyname, ownertag, nametag)
        if isempty(v)
            warn("skipping further EC2 utils test, not all parameters configured")
            return
        end
    end

    info("launching instances...")
    instances = launch_n_ec2(env, 2, ami, insttype, keyname, ownertag, nametag)
    info("launched instances: ", instances)

    info("getting all running instances...")
    running_instances = get_running_instances_by_owner(env, ownertag)
    info("all running instances: ", running_instances)
    @test length(running_instances) >= length(instances)

    hostnames = get_hostnames(env, instances)
    info("hostnames of launched instances: ", hostnames)

    info("terminating instances...")
    terminate_instances_by_owner(env, ownertag)
end

function check_resp(resp, resp_obj_type=Any)
    info("response: ", resp.obj)
    @test 200 <= resp.http_code <= 206
    @test isa(resp, EC2Response)
    @test isa(resp.obj, resp_obj_type)
    @test !isa(resp.obj, EC2Error)
end

function desc_avail_zones(env)
    resp = DescribeAvailabilityZones(env)
    check_resp(resp, DescribeAvailabilityZonesResult)
    resp.obj
end

function terminate_instances(env, instances)
    resp = TerminateInstances(env; instanceId=instances)
    check_resp(resp, TerminateInstancesResult)
    instances
end

function terminate_instances_by_owner(env, owner::String)
    instances = get_instances_by_owner(env, owner)
    println("$owner has the following [$instances] instances")
    terminate_instances(env, instances)
end

get_running_instances_by_owner(env, owner::String) = get_instances_by_owner(env, owner; running_only=true)

function get_instances_by_owner(env, owner::String; running_only::Bool=false)
    instances = String[]
    resp = DescribeInstances(env; filter=[AWS.EC2.Filter(; name="tag:Owner", value=[owner])])
    check_resp(resp, DescribeInstancesResult)
    reservs = resp.obj.reservationSet
    for reserv in reservs
        instancesSet = reserv.instancesSet
        for i in instancesSet
            if !running_only || (running_only && (i.instanceState.code == 16))
                push!(instances, i.instanceId)
            end
        end
    end
    instances
end

function check_running(env, chk_instances)
    new_chk_instances = String[]
    resp = DescribeInstanceStatus(env; instanceId=chk_instances, includeAllInstances=true)
    check_resp(resp, DescribeInstanceStatusResult)
    statuses = resp.obj.instanceStatusSet
    for status in statuses
        println("Status of $(status.instanceId) is $(status.instanceState.code):$(status.instanceState.name)")
        if status.instanceState.code != 16
            push!(new_chk_instances, status.instanceId)
        end
    end
    new_chk_instances
end

function get_hostnames(env, instances)
    names = Vector{Tuple{String,String}}()
    resp = DescribeInstances(env; instanceId=instances)
    check_resp(resp, DescribeInstancesResult)
    reservs = resp.obj.reservationSet
    for reserv in reservs
        instancesSet = reserv.instancesSet
        for i in instancesSet
            push!(names, (i.instanceId, i.dnsName))
        end
    end
    names
end

function launch_n_ec2(env, n::Int, ami::String, insttype::String, keyname::String, ownertag::String, nametag::String)
    # use "c1.medium" or "m1.small"
    resp = RunInstances(env; imageId=ami, instanceType=insttype, minCount=n, maxCount=n, keyName=keyname)
    check_resp(resp, RunInstancesResult)
    instances = String[]
    for inst in resp.obj.instancesSet
        push!(instances, inst.instanceId)
    end
    println("Launched instances : $instances")

    # Tag the instances..also test ec2_basic API..
    resp = CreateTags(env; resourceId=instances, tag=[Tag(;key="Name",value=nametag), Tag(;key="Owner",value=ownertag)])
    check_resp(resp, CreateTagsResult)
    @test resp.obj._return
    println("Tagged instances : $instances")

    # Wait for the instances to come to a running state....
    start = time()
    chk_instances = check_running(env, instances)
    while (((time() - start) < 300.0) && (length(chk_instances) > 0))
        println("All instances not yet in a running state. Waiting and trying again....")
        sleep(5.0)

        chk_instances = check_running(env, chk_instances)
    end

    if (length(chk_instances) > 0)
        println("All instances not yet in a running state. Please check after some time.")
    end
    instances
end
end # module TestEC2Utils


text:
@doc raw"""
    Spin(θ)

Spin layer, with external fields `θ`.
The energy of a layer with units ``s_i`` is given by:

```math
E = -\sum_i \theta_i s_i
```

where each spin ``s_i`` takes values ``\pm 1``.
"""
struct Spin{N, T, A <: AbstractArray{T,N}} <: AbstractLayer{N}
    θ::A
end
Spin(::Type{T}, n::Int...) where {T} = Spin(zeros(T, n...))
Spin(n::Int...) = Spin(Float64, n...)

Base.repeat(l::Spin, n::Int...) = Spin(repeat(l.θ, n...))

free_energies(layer::Spin, inputs::Union{Real,AbstractArray} = 0) = spin_free.(layer.θ .+ inputs)
mode_from_inputs(layer::Spin, inputs::Union{Real,AbstractArray} = 0) = ifelse.(layer.θ .+ inputs .> 0, Int8(1), Int8(-1))
mean_from_inputs(layer::Spin, inputs::Union{Real,AbstractArray} = 0) = tanh.(layer.θ .+ inputs)
mean_abs_from_inputs(layer::Spin, inputs::Union{Real,AbstractArray} = 0) = Ones{Int8}(size(layer))
std_from_inputs(layer::Spin, inputs::Union{Real,AbstractArray} = 0) = sqrt.(var_from_inputs(layer, inputs))

function var_from_inputs(layer::Spin, inputs::Union{Real,AbstractArray} = 0)
    μ = mean_from_inputs(layer, inputs)
    return @. (1 - μ) * (1 + μ)
end

function meanvar_from_inputs(layer::Spin, inputs::Union{Real,AbstractArray} = 0)
    μ = mean_from_inputs(layer, inputs)
    ν = @. (1 - μ) * (1 + μ)
    return μ, ν
end

function sample_from_inputs(layer::Spin, inputs::Union{Real,AbstractArray} = 0)
    θ = layer.θ .+ inputs
    u = rand!(similar(θ))
    return spin_rand.(θ, u)
end

function spin_free(θ::Real)
    abs_θ = abs(θ)
    return -abs_θ - log1pexp(-2abs_θ)
end

spin_rand(θ::Real, u::Real) = ifelse(u < logistic(2θ), Int8(1), Int8(-1))


text:
using Circuitscape

const FILE_CELLMAP_ASC = joinpath("data_files", "cellmap.asc")
const FILE_POINTS_ASC = joinpath("data_files", "points.asc")

hab = Circuitscape.HabitatRasterMap(FILE_CELLMAP_ASC, true)
@assert hab.nodata == -9999
@assert hab.hdr.ncols == 10
@assert hab.hdr.nrows == 10
@assert hab.hdr.nodata == -9999
@assert hab.cond[1,1] == 1.0
@assert hab.cond[1,2] == 0.5

fp = Circuitscape.FocalPoints(FILE_POINTS_ASC)
@assert size(fp.data) == (10,3)
@assert fp.data[3,1] == 3
@assert fp.data[10,2] == 10


text:
module NESInterface

# Import the NES Interface
include("NESFunctions.jl")
end # module


text:
function test_same_layout(layout1, layout2)
    @test all(layout1.positions .== layout2.positions)
    @test all(layout1.scales .== layout2.scales)

    for (elem1, elem2) in zip(layout1.elements, layout2.elements)
        if elem1 isa TeXElement
            @test elem1 == elem2
        else
            test_same_layout(elem1, elem2)
        end
    end
end

@testset "Layout" begin
    @testset "Decorated" begin
        expr = manual_texexpr((:decorated, 'x', 'b', 't'))
        layout = tex_layout(expr, FontFamily())
        @test length(layout.elements) == 3
        @test layout.positions[2][2] < layout.positions[3][2]
        @test layout.scales[1] == 1
        @test layout.scales[2] < 1
        @test layout.scales[3] < 1

        expr = manual_texexpr((:decorated, 'y', nothing, 't'))
        layout = tex_layout(expr, FontFamily())
        @test layout.elements[2] == Space(0)

        expr = manual_texexpr((:decorated, 'z', 'b', nothing))
        layout = tex_layout(expr, FontFamily())
        @test layout.elements[3] == Space(0)
    end

    @testset "Delimited" begin
        expr = manual_texexpr((:delimited, '(', L"\sum_a^b", ')'))
        layout = tex_layout(expr, FontFamily())

        hs = inkheight.(layout.elements) .* layout.scales
        @test hs[1] >= hs[2]
        @test hs[3] >= hs[2]
    end

    @testset "Font" begin
        expr = manual_texexpr((:char, 'u'))
        texchar = tex_layout(expr, FontFamily())
        @test isa(texchar, TeXChar)
        @test texchar.font.style_name == "10 Italic"

        expr = manual_texexpr((:font, :rm, 'u'))
        texchar = tex_layout(expr, FontFamily())
        @test isa(texchar, TeXChar)
        @test texchar.font.style_name == "10 Regular"
    end

    @testset "Group" begin
        expr = manual_texexpr((:group, 'a', 'b', 'c'))
        layout = tex_layout(expr, FontFamily())
        @test length(layout.elements) == 3
        @test length(layout.positions) == 3
        @test length(layout.scales) == 3
        @test all([pos[2] == 0 for pos in layout.positions])
        @test all(layout.scales .== 1)

        expr = manual_texexpr((:group, 'x', (:group, 'y')))
        layout = tex_layout(expr, FontFamily())
        subexpr = manual_texexpr((:group, 'y'))
        sublayout = tex_layout(subexpr, FontFamily())
        @test length(layout.elements) == 2
        @test length(layout.positions) == 2
        @test length(layout.scales) == 2
        test_same_layout(sublayout, layout.elements[2])
    end

    @testset "Space" begin
        elems = generate_tex_elements(L"\quad a")
        @test length(elems) == 1
        char, pos, size = first(elems)
        @test pos[1] == 1

        elems = generate_tex_elements(L"\qquad b")
        @test length(elems) == 1
        char, pos, size = first(elems)
        @test pos[1] == 2
    end
end

@testset "Generate elements" begin
    elems = generate_tex_elements("a + b")
    @test length(elems) == 3

    elems = generate_tex_elements(raw"{{a + b} - {c * d}}")
    @test length(elems) == 7

    # Check the following does not error
    tex = L"\lim_{α →\infty} A^j v_{(a + b)_k}^i \sqrt{2} x!= \sqrt{\frac{1+2}{4+a+x}}\int_{0}^{2π} \sin(x) dx"
    generate_tex_elements(tex)

    tex = L"Momentum $p_x$ (a.u.)"
    generate_tex_elements(tex)

    elems = generate_tex_elements(L"Time $t_0$")
    @test length(elems) == 7
end


text:
module LGCompressed

using GraphIO.CodecZlib
using LightGraphs
using LightGraphs: AbstractGraphFormat

import LightGraphs: loadgraph, loadgraphs, savegraph

export LGCompressedFormat

struct LGCompressedFormat <: AbstractGraphFormat end

function savegraph(fn::AbstractString, g::AbstractGraph, gname::AbstractString,
    format::LGCompressedFormat)
    io = open(fn, "w")
    try
        io = GzipCompressorStream(io)
        return savegraph(io, g, gname, LightGraphs.LGFormat())
    catch
        rethrow()
    finally
        close(io)
    end
end

savegraph(fn::AbstractString, g::AbstractGraph, format::LGCompressedFormat) =
    savegraph(fn, g, "graph", format)

function savegraph(fn::AbstractString, d::Dict{T,U},
    format::LGCompressedFormat) where T <: AbstractString where U <: AbstractGraph
    io = open(fn, "w")
    try
        io = GzipCompressorStream(io)
        return savegraph(io, d, LightGraphs.LGFormat())
    catch
        rethrow()
    finally
        close(io)
    end
end

# savegraph(fn::AbstractString, d::Dict; compress) = savegraph(fn, d, LGCompressedFormat())

loadgraph(fn::AbstractString, gname::AbstractString, format::LGCompressedFormat) =
    loadgraph(fn, gname, LightGraphs.LGFormat())

loadgraph(fn::AbstractString, format::LGCompressedFormat) = loadgraph(fn, "graph", LightGraphs.LGFormat())

loadgraphs(fn::AbstractString, format::LGCompressedFormat) = loadgraphs(fn, LGFormat())

end # module


text:
module CenteredAxes
export Centered

"""
	Centered

Given an array axis of length `N`, return `-N÷2:N÷2-1` for even `N` and `-N÷2:N÷2` for odd `N`.
This may be used to construct or shift an axis of an `OffsetArray`.

# Examples
```julia
julia> A = Float64.(reshape(1:9, 3, 3))
3×3 Array{Float64,2}:
 1.0  4.0  7.0
 2.0  5.0  8.0
 3.0  6.0  9.0

julia> OffsetArray(A, Centered(), Centered())
3×3 OffsetArray(::Array{Float64,2}, -1:1, -1:1) with eltype Float64 with indices -1:1×-1:1:
 1.0  4.0  7.0
 2.0  5.0  8.0
 3.0  6.0  9.0

julia> OffsetArray(A, Centered(), :)
3×3 OffsetArray(::Array{Float64,2}, -1:1, 1:3) with eltype Float64 with indices -1:1×1:3:
 1.0  4.0  7.0
 2.0  5.0  8.0
 3.0  6.0  9.0

julia> A = Float64.(reshape(1:4, 2, 2))
2×2 Array{Float64,2}:
 1.0  3.0
 2.0  4.0

julia> OffsetArray(A, Centered(), Centered())
2×2 OffsetArray(::Array{Float64,2}, -1:0, -1:0) with eltype Float64 with indices -1:0×-1:0:
 1.0  3.0
 2.0  4.0
```
"""
struct Centered end

@inline function Base.to_indices(A, ax::Tuple{Any, Vararg{Any}}, I::Tuple{Centered, Vararg{Any}})
	ax1 = first(ax)
	Lby2 = length(ax1)÷2
	axc = -Lby2:length(ax1) - Lby2 - 1
	(axc, to_indices(A, Base.tail(ax), Base.tail(I))...)
end

end


text:
export ArduinoSurveillanceUnitExternalTemp, ArduinoSurveillanceUnitExternalTempParams

Base.@kwdef struct ArduinoSurveillanceUnitExternalTempParams <: DeviceParams
  portAdress::String
  commandStart::String = "!"
  commandEnd::String = "*"

  pause_ms::Int = 30
  timeout_ms::Int = 500
  delim::String = "#"
  baudrate::Integer = 9600
  ndatabits::Integer = 8
  parity::SPParity = SP_PARITY_NONE
  nstopbits::Integer = 1
end

ArduinoSurveillanceUnitExternalTempParams(dict::Dict) = params_from_dict(ArduinoSurveillanceUnitExternalTempParams, dict)
Base.@kwdef mutable struct ArduinoSurveillanceUnitExternalTemp <: ArduinoSurveillanceUnit
  "Unique device ID for this device as defined in the configuration."
  deviceID::String
  "Parameter struct for this devices read from the configuration."
  params::ArduinoSurveillanceUnitExternalTempParams
  "Flag if the device is optional."
	optional::Bool = false
  "Flag if the device is present."
  present::Bool = false
  "Vector of dependencies for this device."
  dependencies::Dict{String,Union{Device,Missing}}

  ard::Union{SimpleArduino, Nothing} = nothing # Use composition as multiple inheritance is not supported
end

sendCommand(su::ArduinoSurveillanceUnitExternalTemp, cmdString::String) = sendCommand(su.ard, cmdString)

neededDependencies(::ArduinoSurveillanceUnitExternalTemp) = [ArduinoTemperatureSensor] # could in theory be generic temp sensor
optionalDependencies(::ArduinoSurveillanceUnitExternalTemp) = []

function init(su::ArduinoSurveillanceUnitExternalTemp)
  @info "Initializing ArduinoSurveillanceUnitExternalTemp with ID $(su.deviceID)"
  sp = SerialPort(su.params.portAdress)
  open(sp)
	set_speed(sp, su.params.baudrate)
	set_frame(sp, ndatabits=su.params.ndatabits, parity=su.params.parity, nstopbits=su.params.nstopbits)
	# set_flow_control(sp,rts=rts,cts=cts,dtr=dtr,dsr=dsr,xonxoff=xonxoff)
  sleep(2)
  flush(sp)
  write(sp, "!ConnectionEstablished*#")
  response = readuntil(sp, Vector{Char}(su.params.delim), su.params.timeout_ms);
  @info response
  if (response == "ArduinoSurveillanceV1" || response == "ArduinoSurveillanceV2"  )
    @info "Connection to ArduinoSurveillanceUnit established"
    sd = SerialDevice(sp, su.params.pause_ms, su.params.timeout_ms, su.params.delim, su.params.delim)
    su.ard = SimpleArduino(;commandStart = su.params.commandStart, commandEnd = su.params.commandEnd, delim = su.params.delim, sd = sd)
  else
    throw(ScannerConfigurationError(string("Connected to wrong Device", response)))
  end
  su.present = true
end

getTemperatureSensor(su::ArduinoSurveillanceUnitExternalTemp) = dependency(su, ArduinoTemperatureSensor)

function getTemperatures(su::ArduinoSurveillanceUnitExternalTemp; names::Bool=false)
  sensor = getTemperatureSensor(su)
  return getTemperatures(sensor, names = names)
end

function getStatus(su::ArduinoSurveillanceUnitExternalTemp)
  status = sendCommand(su,"GET:STATS")
  return status
end

function resetDAQ(su::ArduinoSurveillanceUnitExternalTemp)
  ACQ = sendCommand(su,"RESET:RP")
  checkACQ(su, ACQ)
end

hasResetDAQ(su::ArduinoSurveillanceUnitExternalTemp) = true


text:
using Test, FundStabABM, Random, StatsBase, StatsPlots, Distributions
using ProgressMeter

include("types_test.jl")
include("params_test.jl")
include("functions_test.jl")


text:
using Distributions, CmdStan, Random, Parameters, StatsFuns
ProjDir = @__DIR__
cd(ProjDir)

stream = open("Logit_Normal.stan","r")
Model = read(stream,String)
close(stream)
seed = 343

stanmodel = Stanmodel(Sample(save_warmup=false, num_warmup=1000,
  num_samples=1000, thin=1), nchains=1, name="Logit_Normal", model=Model,
  printsummary=false, output_format=:mcmcchains, random=CmdStan.Random(seed))

Nsub = 3000
N = 400
ps = logistic.(rand(Normal(-2, 0.8), Nsub))
y = rand.(Binomial.(N, ps))

stanInput = Dict("N"=>N, "Nsub"=>Nsub, "y"=>y)
@time rc, chain, cnames = stan(stanmodel, stanInput, ProjDir)

max_stars_repo_path: temp/Run_Logit_Normal_Sim_Stan .jl

max_issues_repo_path: temp/Run_Logit_Normal_Sim_Stan .jl

max_forks_repo_path: temp/Run_Logit_Normal_Sim_Stan .jl


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

JLLWrappers.@generate_main_file_header("MPFR")
JLLWrappers.@generate_main_file("MPFR", UUID("3a97d323-0669-5f0c-9066-3539efd106a3"))
end  # module MPFR_jll


text:
include("noise_processes.jl")

function polyMap(Θ,z)
  g = 0.
  for i in 1:1:length(Θ)
    g += Θ[i]*z^(i-1)
  end
  return g
end

function noisypolynomial(x₀, θ, noiseproc;n=1, seed=4)
    Random.seed!(seed)
    x = zeros(n)
    x[1] = polyMap(θ,x₀) + rand(noiseproc)
    for i in 2:n
        x[i] = polyMap(θ,x[i-1]) + rand(noiseproc)
    end
    return x
end

max_stars_repo_name: cmerkatas/GSBR

max_issues_repo_name: cmerkatas/GSBR

max_forks_repo_name: cmerkatas/GSBR


text:
using ApproxFun, BandedMatrices, BlockBandedMatrices, LinearAlgebra, BlockArrays, MultivariateOrthogonalPolynomials, Test
using SparseVolterraExamples

#############################
## Tests for Volterra integral operator
## with Clenshaw Jacobi expansion of the kernel.
####
## Includes tests for 1 to x and 1 to 1-x limits of integration.
#############################

@testset "Set #1.1 - Volterra, Clenshaw sin kern, flip = true" begin
    K(x,y)=sin(x+y);
    f=Fun(y->y, Jacobi(0,1,0..1));
    @test triVolterraFullKernelFunP01(f,K,25,true) ≈ Fun(x->-x*cos(2*x)-sin(x)+sin(2*x), Jacobi(1,0,0..1));
end

@testset "Set #1.2 - Volterra, Clenshaw sin kern, flip = false, val test" begin
    K(x,y)=sin(x+y);
    val = rand(0..1);
    f=Fun(y->y, Jacobi(0,1,0..1));
    @test triVolterraFullKernelFunP01(f,K,30,false)(val) ≈ Fun(x->((-1)*cos(1)+x*cos(1)+sin(1)-sin(x)), Jacobi(1,0,0..1))(val);
end

@testset "Set #2 - Volterra, Clenshaw sinh kern, flip = true" begin
    K(x,y)=(y-x)*sinh(x);
    f=Fun(y->y, Jacobi(0,1,0..1));
    @test triVolterraFullKernelFunP01(f,K,25,true) ≈ Fun(x->-1/6*x^3*sinh(x), Jacobi(1,0,0..1));
end

@testset "Set #3 - Volterra, Clenshaw sinh kern, flip = false" begin
    K(x,y)=(y-x)*sinh(x^2);
    f=Fun(y->y, Jacobi(0,1,0..1));
    @test triVolterraFullKernelFunP01(f,K,25,false) ≈ Fun(x->(sinh(x^2))/3-3/2*x*sinh(x^2)+2*x^2*sinh(x^2)-5/6*x^3*sinh(x^2), Jacobi(0,1,0..1));
end


text:
function test_argparse()
    install_download_args = argparse(["install", "download", "0.3"])
    @test install_download_args == Dict(
        "%COMMAND%" => "install",
        "install" => Dict(
            "%COMMAND%" => "download",
            "download" => Dict(
                "labels" => AbstractString[],
                "version" => v"0.3",
            )
        )
    )

    install_link_args = argparse(["install", "link", "/path/to/julia", "--labels", "julia-src"])
    @test install_link_args == Dict(
        "%COMMAND%" => "install",
        "install" => Dict(
            "%COMMAND%" => "link",
            "link" => Dict(
                "labels" => AbstractString["julia-src"],
                "dir" => "/path/to/julia"
            )
        )
    )

    create_args1 = argparse(["create"])
    @test create_args1 == Dict(
        "%COMMAND%" => "create",
        "create" => Dict(
            "dir" => "",
            "requirements" => "",
            "name" => "",
            "julia-version" => "",
            "req-type" => :REQUIRE,
        )
    )

    create_args2 = argparse(
        [
            "create", "/path/to/playground",
            "--name", "myplayground",
            "--julia-version", "julia-0.3",
            "--requirements", "/path/to/requirements",
            "--req-type", "DECLARE"
        ]
    )
    @test create_args2 == Dict(
        "%COMMAND%" => "create",
        "create" => Dict(
            "dir" => "/path/to/playground",
            "requirements" => "/path/to/requirements",
            "name" => "myplayground",
            "julia-version" => "julia-0.3",
            "req-type" => :DECLARE,
        )
    )

    activate_args1 = argparse(["activate"])
    @test activate_args1 == Dict(
        "%COMMAND%" => "activate",
        "activate" => Dict(
            "dir" => "",
            "name" => ""
        )
    )

    activate_args2 = argparse(["activate", "--name", "myplayground"])
    @test activate_args2 == Dict(
        "%COMMAND%" => "activate",
        "activate" => Dict(
            "dir" => "",
            "name" => "myplayground"
        )
    )

    exec_args = argparse(["exec", "ls -al", "--name", "myplayground"])
    @test exec_args == Dict(
        "%COMMAND%" => "exec",
        "exec" => Dict(
            "cmd" => "ls -al",
            "dir" => "",
            "name" => "myplayground"
        )
    )

    list_args = argparse(["list", "--show-links"])
    @test list_args == Dict(
        "%COMMAND%" => "list",
        "list" => Dict(
            "show-links" => true
        )
    )

    clean_args = argparse(["clean"])
    @test clean_args == Dict(
        "%COMMAND%" => "clean",
        "clean" => Dict()
    )

    rm_args = argparse(["rm", "myplayground"])
    @test rm_args == Dict(
        "%COMMAND%" => "rm",
        "rm" => Dict(
            "dir" => "",
            "name" => "myplayground"
        )
    )
end

test_argparse()


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCellDG
include("../useful_routines.jl")

function test_cut_simple_tension()
    x0 = [0.0, 0.0]
    widths = [4.0, 1.0]
    nelements = [2, 2]

    interfacepoint = [1.0, 0.0]
    interfaceangle = 10.0
    interfacenormal = [cosd(interfaceangle), sind(interfaceangle)]

    lambda, mu = 1.0, 2.0
    stiffness = CutCellDG.HookeStiffness(lambda, mu, lambda, mu)
    penalty = 10.0
    eta = 1
    dx = 0.1
    e11 = dx / 4.0
    e22 = -lambda / (lambda + 2mu) * e11
    dy = e22

    polyorder = 2
    numqp = required_quadrature_order(polyorder)

    elasticitybasis = LagrangeTensorProductBasis(2, polyorder)
    levelsetbasis = HermiteTensorProductBasis(2)
    # quad = tensor_product_quadrature(2,4)
    dim,nf = size(interpolation_points(levelsetbasis))
    refpoints = interpolation_points(elasticitybasis)

    cgmesh = CutCellDG.CGMesh(x0, widths, nelements, nf)
    mesh = CutCellDG.DGMesh(x0, widths, nelements, refpoints)

    levelset = CutCellDG.LevelSet(
        x -> plane_distance_function(x, interfacenormal, interfacepoint),
        cgmesh,
        levelsetbasis,
    )

    cutmesh = CutCellDG.CutMesh(mesh, levelset)

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

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

    CutCellDG.assemble_displacement_bilinear_forms!(
        sysmatrix,
        elasticitybasis,
        cellquads,
        stiffness,
        cutmesh,
    )
    CutCellDG.assemble_interelement_condition!(
        sysmatrix,
        elasticitybasis,
        facequads,
        stiffness,
        cutmesh,
        penalty,
        eta,
    )
    CutCellDG.assemble_penalty_displacement_component_bc!(
        sysmatrix,
        sysrhs,
        x -> 0.0,
        elasticitybasis,
        facequads,
        stiffness,
        cutmesh,
        x -> x[1] ≈ 0.0,
        [1.0, 0.0],
        penalty,
    )
    CutCellDG.assemble_penalty_displacement_component_bc!(
        sysmatrix,
        sysrhs,
        x -> 0.0,
        elasticitybasis,
        facequads,
        stiffness,
        cutmesh,
        x -> x[2] ≈ 0.0,
        [0.0, 1.0],
        penalty,
    )
    CutCellDG.assemble_penalty_displacement_component_bc!(
        sysmatrix,
        sysrhs,
        x -> dx,
        elasticitybasis,
        facequads,
        stiffness,
        cutmesh,
        x -> x[1] ≈ widths[1],
        [1.0, 0.0],
        penalty,
    )
    CutCellDG.assemble_coherent_interface_condition!(
        sysmatrix,
        elasticitybasis,
        interfacequads,
        stiffness,
        cutmesh,
        penalty,
        eta,
    )

    op = CutCellDG.sparse_displacement_operator(sysmatrix, cutmesh)
    rhs = CutCellDG.displacement_rhs_vector(sysrhs, cutmesh)
    #
    sol = op \ rhs
    disp = reshape(sol, 2, :)

    nodalcoordinates = CutCellDG.nodal_coordinates(cutmesh)
    nodalcoordinates = hcat(nodalcoordinates, nodalcoordinates[:, 1:18])

    testdisp = copy(nodalcoordinates)
    testdisp[1, :] .*= e11
    testdisp[2, :] .*= e22

    @test allapprox(disp, testdisp, 1e3eps())
end

test_cut_simple_tension()