text:
# Following work by @jiahao, we compute character widths using a combination of
#   * advance widths from GNU Unifont (advance width 512 = 1 en)
#   * UAX 11: East Asian Width
#   * a few exceptions as needed
# Adapted from http://nbviewer.ipython.org/gist/jiahao/07e8b08bf6d8671e9734
#
# Requires Julia (obviously) and FontForge.

#############################################################################
CharWidths = Dict{Int,Int}()

#############################################################################
# Use ../libutf8proc for category codes, rather than the one in Julia,
# to minimize bootstrapping complexity when a new version of Unicode comes out.
catcode(c) = ccall((:utf8proc_category,"../libutf8proc"), Cint, (Int32,), c)

# use Base.UTF8proc module to get category codes constants, since
# we won't change these in utf8proc.
import Base.UTF8proc

#############################################################################
# Use a default width of 1 for all character categories that are
# letter/symbol/number-like.  This can be overriden by Unifont or UAX 11
# below, but provides a useful nonzero fallback for new codepoints when
# a new Unicode version has been released but Unifont hasn't been updated yet.

zerowidth = Set{Int}() # categories that may contain zero-width chars
push!(zerowidth, UTF8proc.UTF8PROC_CATEGORY_CN)
push!(zerowidth, UTF8proc.UTF8PROC_CATEGORY_MN)
push!(zerowidth, UTF8proc.UTF8PROC_CATEGORY_MC)
push!(zerowidth, UTF8proc.UTF8PROC_CATEGORY_ME)
push!(zerowidth, UTF8proc.UTF8PROC_CATEGORY_SK)
push!(zerowidth, UTF8proc.UTF8PROC_CATEGORY_ZS)
push!(zerowidth, UTF8proc.UTF8PROC_CATEGORY_ZL)
push!(zerowidth, UTF8proc.UTF8PROC_CATEGORY_ZP)
push!(zerowidth, UTF8proc.UTF8PROC_CATEGORY_CC)
push!(zerowidth, UTF8proc.UTF8PROC_CATEGORY_CF)
push!(zerowidth, UTF8proc.UTF8PROC_CATEGORY_CS)
push!(zerowidth, UTF8proc.UTF8PROC_CATEGORY_CO)
for c in 0x0000:0x110000
    if catcode(c) ∉ zerowidth
        CharWidths[c] = 1
    end
end

#############################################################################
# Widths from GNU Unifont

#Read sfdfile for character widths
function parsesfd(filename::AbstractString, CharWidths::Dict{Int,Int}=Dict{Int,Int}())
    state=:seekchar
    lineno = 0
    codepoint = width = nothing
    for line in readlines(open(filename))
        lineno += 1
        if state==:seekchar         #StartChar: nonmarkingreturn
            if contains(line, "StartChar: ")
                codepoint = nothing
                width = nothing
                state = :readdata
            end
        elseif state==:readdata #Encoding: 65538 -1 2, Width: 1024
            contains(line, "Encoding:") && (codepoint = parse(Int, split(line)[3]))
            contains(line, "Width:") && (width = parse(Int, split(line)[2]))
            if codepoint!=nothing && width!=nothing && codepoint >= 0
                w=div(width, 512) # 512 units to the en
                if w > 0
                    # only add nonzero widths, since (1) the default is zero
                    # and (2) this circumvents some apparent bugs in Unifont
                    # (https://savannah.gnu.org/bugs/index.php?45395)
                    CharWidths[codepoint] = w
                end
                state = :seekchar
            end
        end
    end
    CharWidths
end
CharWidths=parsesfd("unifont.sfd", CharWidths)
CharWidths=parsesfd("unifont_upper.sfd", CharWidths)

#############################################################################
# Widths from UAX #11: East Asian Width
#   .. these take precedence over the Unifont width for all codepoints
#      listed explicitly as wide/full/narrow/half-width

for line in readlines(open("EastAsianWidth.txt"))
    #Strip comments
    (isempty(line) || line[1] == '#') && continue
    precomment = split(line, '#')[1]
    #Parse code point range and width code
    tokens = split(precomment, ';')
    length(tokens) >= 2 || continue
    charrange = tokens[1]
    width = strip(tokens[2])
    #Parse code point range into Julia UnitRange
    rangetokens = split(charrange, "..")
    charstart = parse(UInt32, "0x"*rangetokens[1])
    charend = parse(UInt32, "0x"*rangetokens[length(rangetokens)>1 ? 2 : 1])

    #Assign widths
    for c in charstart:charend
        if width=="W" || width=="F" # wide or full
            CharWidths[c]=2
        elseif width=="Na"|| width=="H" # narrow or half
            CharWidths[c]=1
        end
    end
end

#############################################################################
# A few exceptions to the above cases, found by manual comparison
# to other wcwidth functions and similar checks.

for c in keys(CharWidths)
    cat = catcode(c)

    # make sure format control character (category Cf) have width 0,
    # except for the Arabic characters 0x06xx (see unicode std 6.2, sec. 8.2)
    if cat==UTF8proc.UTF8PROC_CATEGORY_CF && c ∉ [0x0601,0x0602,0x0603,0x06dd]
        CharWidths[c]=0
    end

    # Unifont has nonzero width for a number of non-spacing combining
    # characters, e.g. (in 7.0.06): f84,17b4,17b5,180b,180d,2d7f, and
    # the variation selectors
    if cat==UTF8proc.UTF8PROC_CATEGORY_MN
        CharWidths[c]=0
    end

    # We also assign width of zero to unassigned and private-use
    # codepoints (Unifont includes ConScript Unicode Registry PUA fonts,
    # but since these are nonstandard it seems questionable to recognize them).
    if cat==UTF8proc.UTF8PROC_CATEGORY_CO || cat==UTF8proc.UTF8PROC_CATEGORY_CN
        CharWidths[c]=0
    end

    # for some reason, Unifont has width-2 glyphs for ASCII control chars
    if cat==UTF8proc.UTF8PROC_CATEGORY_CC
        CharWidths[c]=0
    end
end

#By definition, should have zero width (on the same line)
#0x002028 '
' category: Zl name: LINE SEPARATOR/
#0x002029 '
' category: Zp name: PARAGRAPH SEPARATOR/
CharWidths[0x2028]=0
CharWidths[0x2029]=0

#By definition, should be narrow = width of 1 en space
#0x00202f ' ' category: Zs name: NARROW NO-BREAK SPACE/
CharWidths[0x202f]=1

#By definition, should be wide = width of 1 em space
#0x002001 ' ' category: Zs name: EM QUAD/
#0x002003 ' ' category: Zs name: EM SPACE/
CharWidths[0x2001]=2
CharWidths[0x2003]=2

#############################################################################
# Output (to a file or pipe) for processing by data_generator.rb
# ... don't bother to output zero widths since that will be the default.

firstc = 0x000000
lastv = 0
uhex(c) = uppercase(hex(c,4))
for c in 0x0000:0x110000
    v = get(CharWidths, c, 0)
    if v != lastv || c == 0x110000
        v < 4 || error("invalid charwidth $v for $c")
        if firstc+1 < c
            println(uhex(firstc), "..", uhex(c-1), "; ", lastv)
        else
            println(uhex(firstc), "; ", lastv)
        end
        firstc = c
        lastv = v
    end
end


text:
println("Loading Libraries ....", "\n")

## libraries
using DataFrames,CSV, Dates, BSON, GLM, LinearAlgebra

## pull in functions
include("discharge_functions.jl")

## Take time and site input input
site = get_site_input()
time = get_dt_input()

println("Pulling Data From NWIS ....", "\n")
## Pull Data From NWIS from Fort Wayne Detroit River 04168510
cont_dis = nwis_pull(time)
## Clean continous discharge data
clean_dis!(cont_dis)

if site == "04165705"
	println("Estimated Discharge at Detroit River 04165705:  ",get_lag(time,cont_dis))
else
	## Make a DataFrame for date time of prediciton
	bel = DataFrame(datetime = time)
	make_lead_lag!(bel)

	## insert lead and lag values of discharge, gage height, and velocity into DataFrame
	fill_lead_lag!(bel, cont_dis)
	impute_z!(bel)

	println("Predicting Discharge ....", "\n")

	## Load Linear Model
	BSON.@load "./Discharge_Predict_Detroit_LinearRegression__20220321.bson" bel_lm_total

	## predict discharge
	println("\n")
	print("Estimated Discharge at Detroit River 04168557:  ",  round( GLM.predict(bel_lm_total, select(bel, Not(:datetime)))[1], digits=0), " cfs")
	println("\n")
end

ext: jl

max_stars_repo_name: atott/Detroit-Discharge-Prediction

max_issues_repo_name: atott/Detroit-Discharge-Prediction

max_forks_repo_name: atott/Detroit-Discharge-Prediction


text:
module RunMyNotes

export folder, package

using Literate, Conda

"""
    folder(in, [out])
Runs all files `*.jl` found in folder `in`, generating a Jupyter `.ipynb` notebook for each,
saved in the same folder or in `out` if provided.

Keywords:
* `html=true` then converts this notebook to HTML too, using `jupyter nbconvert`.
* `throw=true` exits with an error if any errors were encountered .
* `all=true` runs everything, `all=false` will skip any `.jl` files which are older than their `.ipynb` versions.
* 'verbose=true` adds details about time & computer at the end.
* `credit=true` adds links.
"""
function folder(indir::String, outdir::String=indir;
        html::Bool=true, throw::Bool=true, all::Bool=true, credit::Bool=true, verbose::Bool=true)

    fulldir = normpath(expanduser(indir))
    fullout = normpath(expanduser(outdir))
    @assert isdir(fulldir) "can't find input folder $fulldir"
    @assert isdir(fullout) "can't find output folder $fullout"

    olddir = pwd() # should just fix run command instead

    pre = verbose ? appendinfo : identity
    post = credit ? appendcredit : identity

    list = filter(s->endswith(s,".jl"), readdir(fulldir))

    errs = []

    for name in list
        jlname = joinpath(fulldir, name)
        nbname = joinpath(fullout, name)[1:end-2] * "ipynb"

        if isfile(nbname) && mtime(nbname) > mtime(jlname)
            all || continue # all=false means skip files for which .jl hasn't changed
        end

        ran = true
        nb = ""
        try
            nb = Literate.notebook(jlname, fullout; credit=credit, preprocess=pre, postprocess=post)
        catch err
            @error "error in creating or running notebook from $jlname \n" err
            push!(errs, err)
            ran = false
        end

        if html && ran
            try
                cd(Conda.PYTHONDIR)
                run(`./jupyter nbconvert $nb`) # writes to same folder
            catch err
                @error "error in converting notebook $nb \n" err
                push!(errs, err)
            end
        end
    end

    cd(olddir)

    if throw && length(errs)>0
        error("encountered $(length(errs)) errors executing & converting notebooks, details printed above")
    end

    return length(errs)==0 # for @test
end

"""
    package(ModuleName)
This essentially runs `RunMyNotes.foler("~/.julia/dev/ModuleName/notes/")`.
Change keyword `sub="notes"` to select a different folder.
"""
function package(mod::Module, outdir=nothing; sub::String="notes", kw...)

    indir = joinpath(dirname(pathof(mod)), "..", sub) # pathof gives "...dev/src/RunMyNotes.jl"
    if outdir==nothing
        outdir = indir # these will get expanduser |> normpath
    end

    return folder(indir, outdir; kw...)
end

using Dates, InteractiveUtils

function infostring()
    io = IOBuffer();
    InteractiveUtils.versioninfo(io)
    s = String(take!(io))
    vinfo = join(split(s, "\n")[[1,4,5]],",") # Julia Version, OS, CPU
    return string(DateTime(now()), ",  ", vinfo)
end

appendinfo(str) = str * string("\n # ----- \n # *", infostring(), "* ") # not entirely happy with \n here

function appendcredit(dict::Dict)
    str = dict["cells"][end]["source"][end]
    rep = replace(str, "/fredrikekre/Literate.jl)" =>
        "/fredrikekre/Literate.jl), called by [RunMyNotes](http://github.com/mcabbott/RunMyNotes.jl)")
    dict["cells"][end]["source"][end] = rep
    return dict
end

end # module


text:
struct RawSMPS{T <: AbstractFloat}
    name::String
    tim::RawTime
    cor::RawCor{T}
    sto::RawStoch

    function RawSMPS(tim::RawTime, cor::RawCor{T}, sto::RawStoch) where T <: AbstractFloat
        tim.name == cor.name && tim.name == sto.name || error("Inconsistent names of SMPS files.")
        return new{T}(tim.name, tim, cor, sto)
    end
end

function Base.read(io::IO, ::Type{RawSMPS})
    # Get filepath
    m = match(r"<file (.*).smps>", io.name)
    m === nothing && error("SMPS specification is malformed. Correct filename usage: /path/to/smps/files.smps")
    path = only(m.captures)
    # Parse tim file
    timfile = path * ".tim"
    tim = parse_tim(timfile)
    # Parse cor file
    corfile = path * ".cor"
    cor = parse_cor(corfile)
    # Parse sto file
    stofile = path * ".sto"
    sto = parse_sto(tim, cor, stofile)
    return RawSMPS(tim, cor, sto)
end

function Base.read(io::IO, ::Type{RawSMPS{T}}) where T <: AbstractFloat
    # Get filepath
    m = match(r"<file (.*).smps>", io.name)
    m === nothing && error("SMPS specification is malformed. Correct filename usage: /path/to/smps/files.smps")
    path = only(m.captures)
    # Parse tim file
    timfile = path * ".tim"
    tim = parse_tim(T, timfile)
    # Parse cor file
    corfile = path * ".cor"
    cor = parse_cor(T, corfile)
    # Parse sto file
    stofile = path * ".sto"
    sto = parse_sto(T, tim, cor, io)
    return RawSMPS(tim, cor, sto)
end


text:
module AllSubsetRegression

    using Printf

    using Distributed, Distributions, SharedArrays, LinearAlgebra
    using ..ModelSelection

    export to_string, AllSubsetRegressionResult, ALLSUBSETREGRESSION_EXTRAKEY
    export ols, ols!

    include("const.jl")
    include("strings.jl")
    include("utils.jl")
    include("structs/result.jl")
    include("core.jl")
    include("estimators/ols.jl")
end


text:
function _lll_gram(A::NfOrdIdl)
  K = nf(order(A))
  @assert istotally_real(K)
  g = trace_matrix(A)
  @hassert :ClassGroup 1 !iszero(det(g))
  @hassert :ClassGroup 1 isposdef(g)
  l, t = lll_gram_with_transform(g)
  return FakeFmpqMat(l, fmpz(1)), t::fmpz_mat
end

function _lll_quad(A::NfOrdIdl)
  K = nf(order(A))
  @assert degree(K) == 2 && discriminant(order(A)) < 0
  b = basis(A)
  a1 = 2*numerator(norm(b[1]))
  a2 = 2*numerator(norm(b[2]))
  a12 = numerator(trace(b[1] * conjugate_quad(K(b[2]))))
  g = matrix(FlintZZ, 2, 2, [a1, a12, a12, a2])
  @hassert :ClassGroup 1 isposdef(g)
  l, t = lll_gram_with_transform(g)
  return FakeFmpqMat(l, fmpz(1)), t::fmpz_mat
end

function lll_basis(A::NfOrdIdl, v::fmpz_mat = zero_matrix(FlintZZ, 1, 1); prec::Int = 100)
  L, T = lll(A, v, prec=prec)
  S = FakeFmpqMat(T)*basis_matrix(A)*basis_matrix(order(A))
  K = nf(order(A))
  q = nf_elem[elem_from_mat_row(K, numerator(S), i, denominator(S)) for i=1:degree(K)]
  return q
end

function cmpindex(A::fmpz_mat, i::Int, j::Int, b::fmpz)
  a = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), A, i-1, j-1)
  return ccall((:fmpz_cmp, :libflint), Int32, (Ptr{fmpz}, Ref{fmpz}), a, b)
end

function prod_diagonal(A::fmpz_mat)
  a = one(fmpz)
  GC.@preserve a begin
    for i=1:nrows(A)
      b = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), A, i - 1, i - 1)
      ccall((:fmpz_mul, :libflint), Nothing, (Ref{fmpz}, Ref{fmpz}, Ptr{fmpz}), a, a, b)
    end
  end
  return a
end

function lll(A::NfOrdIdl, v::fmpz_mat = zero_matrix(FlintZZ, 1, 1); prec::Int = 100)

  K = nf(order(A))

  if iszero(v) && istotally_real(K)
    #in this case the gram-matrix of the minkowski lattice is the trace-matrix
    #which is exact.
    return _lll_gram(A)
  end

  if iszero(v) && degree(K) == 2 && discriminant(order(A)) < 0
    #in this case the gram-matrix of the minkowski lattice is related to the
    #trace-matrix which is exact.
    #could be extended to CM-fields
    return _lll_quad(A)
  end

  n = degree(order(A))
  prec = max(prec, 4*n)

  l, t1 = lll_with_transform(basis_matrix(A, copy = false))

  if iszero(v)
    d = minkowski_gram_mat_scaled(order(A), prec)
    ccall((:fmpz_mat_mul, :libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz_mat}), d, d, l')
    ccall((:fmpz_mat_mul, :libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz_mat}), d, l, d)
    g = zero_matrix(FlintZZ, n, n)
    den = fmpz(1)
    sv = fmpz(0)
  else
    c = minkowski_matrix(nf(order(A)), prec) ## careful: current iteration
                                          ## c is NOT a copy, so don't change.
    b = FakeFmpqMat(l)*basis_matrix(order(A), copy = false)

    rt_c = roots_ctx(K)
    if !isdefined(rt_c, :cache_z1)
      rt_c.cache_z1 = zero_matrix(FlintZZ, n, n)
      rt_c.cache_z2 = zero_matrix(FlintZZ, n, n)
    end

    d::fmpz_mat = rt_c.cache_z1
    g::fmpz_mat = rt_c.cache_z2

    round_scale!(g, c, prec)

    sv = fmpz(0)
    if !iszero(v)
      @v_do :ClassGroup 2 println("using inf val", v)
      c = deepcopy(c)
      mult_by_2pow_diag!(c, v)
      sv = max(fmpz(0), sum(v[1,i] for i=1:ncols(l)))
    end

    round_scale!(d, c, prec)
    ccall((:fmpz_mat_mul, :libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat},  Ref{fmpz_mat}), g, (b.num), d)
    den = b.den

    ccall((:fmpz_mat_gram, :libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat}), d, g)
    shift!(d, -prec)
  end

  prec = div(prec, 2)
  shift!(d, -prec)  #TODO: remove?

  for i=1:n
    fmpz_mat_entry_add_ui!(d, i, i, UInt(nrows(d)))
  end

  ctx = Nemo.lll_ctx(0.99, 0.51, :gram)

  ccall((:fmpz_mat_one, :libflint), Nothing, (Ref{fmpz_mat}, ), g)
  ccall((:fmpz_lll, :libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{Nemo.lll_ctx}), d, g, ctx)

  l, t = d, g
  ## test if entries in l are small enough, if not: increase precision
  ## or signal that prec was too low
  @v_do :ClassGroup 2 printstyled("lll basis length profile\n", color=:green);
  @v_do :ClassGroup 2 for i=1:nrows(l)
    print(Float64(div(l[i,i], fmpz(2)^prec*den*den)*1.0), " : ")
  end
  @v_do :ClassGroup 2 println("")
  if nbits(maximum(abs, t)) >  div(prec, 2)
    @v_do :ClassGroup 2 printstyled("lll trafo too large\n", color=:red);
    throw(LowPrecisionLLL())
  end
  ## lattice has lattice disc = order_disc * norm^2
  ## lll needs to yield a basis sth
  ## l[1,1] = |b_i|^2 <= 2^((n-1)/2) disc^(1/n)
  ## and prod(l[i,i]) <= 2^(n(n-1)/2) disc
  n = nrows(l)
  disc = abs(discriminant(order(A)))*norm(A)^2 * den^(2*n) * fmpz(2)^(2*sv)
  di = root(disc, n)+1
  di *= fmpz(2)^(div(n+1,2)) * fmpz(2)^prec

  if cmpindex(l, 1, 1, di) > 0
    @v_do :ClassGroup 2 printstyled("LLL basis too large\n", color = :red);
    @v_do :ClassGroup 3 println("bound is ", di, " value at ", 1, " is ", l[1,1]);
    throw(LowPrecisionLLL())
  end
  pr = prod_diagonal(l)
  if pr > fmpz(2)^(div(n*(n-1), 2)) * disc * fmpz(2)^(n*prec)
    @v_do :ClassGroup 2 printstyled("LLL basis too large\n", color = :red);
    @v_do :ClassGroup 2 println("prod too large: ", pr, " > 2^(n(n-1)/2) disc = ", fmpz(2)^(div(n*(n-1), 2)) * disc * fmpz(2)^(n*prec));
    throw(LowPrecisionLLL())
  end

  return FakeFmpqMat(deepcopy(l), fmpz(2)^prec), t*t1
end

function short_elem(A::NfOrdFracIdl,
                v::fmpz_mat=zero_matrix(FlintZZ, 1,1); prec::Int = 100)
  return divexact(short_elem(A.num, v, prec = prec), A.den)
end
function _short_elem(A::NfOrdFracIdl,
                v::fmpz_mat=zero_matrix(FlintZZ, 1, 1))
  return divexact(_short_elem(A.num, v), A.den)
end

function short_elem(A::NfOrdIdl,
                v::fmpz_mat = zero_matrix(FlintZZ, 1,1); prec::Int = 100)
  K = nf(order(A))
  temp = basis_matrix(A, copy = false)*basis_matrix(order(A), copy = false)
  b = temp.num
  b_den = temp.den
  local t
  while true
    if prec > 2^18
      error("Something wrong in short_elem")
    end
    try
      l, t = lll(A, v, prec = prec)
      break
    catch e
      if !(e isa LowPrecisionLLL || e isa InexactError)
        rethrow(e)
      end
    end
    prec = 2 * prec
  end

  w = view(t, 1:1, 1:ncols(t))
  c = w*b
  q = elem_from_mat_row(K, c, 1, b_den)
  return q
end

function _short_elem(A::NfOrdIdl,
                v::fmpz_mat = zero_matrix(FlintZZ, 1, 1))
  p = 64
  while true
    if p > 40000
      error("Something wrong in reduce_ideal")
    end
    try
      b = short_elem(A, v, prec = p)
      return b
    catch e
      if e isa LowPrecisionLLL || e isa InexactError
        p = 2*p
      else
        rethrow(e)
      end
    end
  end
end

function reduce_ideal(A::NfOrdIdl)
  B = inv(A)
  success = false
  b = _short_elem(B)
  C = b*A
  simplify(C)
  @assert C.den == 1
  return C.num
end

function reduce_ideal(A::NfOrdFracIdl)
  B = inv(A)
  b = _short_elem(B.num)
  C = divexact(b, B.den)*A
  simplify(C)
  @assert C.den == 1
  return C.num
end

function reduce_ideal2(A::NfOrdIdl)
  B = inv(A)
  b = _short_elem(B)
  C = b*A
  simplify(C)
  @assert C.den == 1
  return C.num, b
end


text:
using Test
using JustSayIt
using PyCall
import JustSayIt: DEFAULT_MODEL_NAME, TYPE_MODEL_NAME, MODELDIR_PREFIX, DEFAULT_NOISES, DEFAULT_READER_ID, AUDIO_ELTYPE
import JustSayIt: init_jsi, finalize_jsi, reader, start_reading, stop_reading, set_default_streamer, default_streamer

# Test setup
const SAMPLEDIR_CMD = joinpath("samples", "commands")

commands = Dict("help"      => Help.help,
                "type"      => Keyboard.type)
modeldirs = Dict(DEFAULT_MODEL_NAME => joinpath(MODELDIR_PREFIX, "vosk-model-small-en-us-0.15"),
                 TYPE_MODEL_NAME    => joinpath(MODELDIR_PREFIX, "vosk-model-small-en-us-0.15"))

init_jsi(commands, modeldirs, DEFAULT_NOISES)
start_reading(joinpath(SAMPLEDIR_CMD, "help.wav"))

@testset "$(basename(@__FILE__))" begin
    @testset "1. readbytes! (from wave reader)" begin
        wav_reader = reader()
        sample = zeros(UInt8, wav_reader.getnframes()*sizeof(AUDIO_ELTYPE))
        bytes_read = readbytes!(wav_reader, sample)
        @test bytes_read == 68938
    end;
    @testset "2. set default streamer" begin
        @test_throws KeyError default_streamer()
        set_default_streamer(reader, DEFAULT_READER_ID)
        @test isa(default_streamer(), PyObject)
    end
end;

# Test tear down
stop_reading()
finalize_jsi()


text:
# This file was generated, do not modify it. # hide
figure(figsize=(8,6))
plot(curve.parameter_values, curve.measurements)
xlabel("Number of rounds", fontsize=14)
ylabel("HingeLoss", fontsize=14)
xticks([10, 100, 200, 500], fontsize=12)

savefig(joinpath(@OUTPUT, "EX-crabs-xgb-curve1.svg")) # hide


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Cylinder

"""
    Cylinder <: SimpleShape

A [`SimpleShape`](@ref) representing a cylinder of unitary height and diameter.

# Members

- `transformation::Transformation`: the `Transformation` associated with the cylinder.
- `material::Material`: the [`Material`](@ref) of the cylinder.
"""
Base.@kwdef struct Cylinder <: SimpleShape
    transformation::Transformation = Transformation()
    material::Material = Material()
end

@doc """
    Cylinder(transformation::Transformation, material::Material)

Constructor for a [`Cylinder`](@ref) instance.
""" Cylinder(::Transformation, ::Material)

@doc """
    Cylinder(transformation::Transformation = Transformation(),
           material::Material = Material())

Constructor for a [`Cylinder`](@ref) instance.
""" Cylinder(; ::Transformation, ::Material)

function get_t(::Type{Cylinder}, ray::Ray)
    sray = scaling(2) * ray
    ox, oy, oz = sray.origin.v
    dx, dy, dz = sray.dir

    # check if side is hit
    a = dx^2 + dy^2
    b = 2 * (ox * dx + oy * dy)
    c = ox^2 + oy^2 - 1
    Δ = b^2 - 4f0 * a * c
    if !(iszero(a) && iszero(b) && iszero(Δ)) && Δ >= 0
        sqrt_Δ = sqrt(Δ)
        t_1 = (-b - sqrt_Δ) / (2f0 * a)
        t_2 = (-b + sqrt_Δ) / (2f0 * a)
        # nearest point
        # @assert !isnan(t_1)
        # @assert !isnan(t_2)
        ray.tmin < t_1 < ray.tmax && abs(oz + t_1 * dz) <= 1f0 && return t_1
        ray.tmin < t_2 < ray.tmax && abs(oz + t_2 * dz) <= 1f0 && abs(oz) <= 1f0 && return t_2
    end

    # check if caps are hit
    tz1, tz2 = minmax(( 1f0 - oz) / dz, (-1f0 - oz) / dz)

    ray.tmin < tz1 < ray.tmax && (ox + tz1 * dx)^2 + (oy + tz1 * dy)^2 <= 1f0 && return tz1
    ray.tmin < tz2 < ray.tmax && (ox + tz2 * dx)^2 + (oy + tz2 * dy)^2 <= 1f0 && return tz2
    return Inf32
end

function get_all_ts(::Type{Cylinder}, ray::Ray)
    res = Vector{Float32}()
    sizehint!(res, 2)
    sray = scaling(2) * ray
    ox, oy, oz = sray.origin.v
    dx, dy, dz = sray.dir

    # check if side is hit
    a = dx^2 + dy^2
    b = 2 * (ox * dx + oy * dy)
    c = ox^2 + oy^2 - 1
    Δ = b^2 - 4f0 * a * c
    if !(iszero(a) && iszero(b) && iszero(Δ)) && Δ >= 0
        sqrt_Δ = sqrt(Δ)
        t_1 = (-b - sqrt_Δ) / (2f0 * a)
        t_2 = (-b + sqrt_Δ) / (2f0 * a)
        # nearest point
        # @assert !isnan(t_1)
        # @assert !isnan(t_2)
        abs(oz + t_1 * dz) <= 1f0 && push!(res, t_1)
        abs(oz + t_2 * dz) <= 1f0 && push!(res, t_2)
        length(res) == 2 && return res
    end

    # check if caps are hit
    tz1, tz2 = minmax(( 1f0 - oz) / dz, (-1f0 - oz) / dz)

    # ⪅(x::Number, y::Number) = x < y || x ≈ y
    # (ox + tz1 * dx)^2 + (oy + tz1 * dy)^2 ⪅ 1f0 && push!(res, tz1)
    # length(res) == 2 && return res
    # (ox + tz2 * dx)^2 + (oy + tz2 * dy)^2 ⪅ 1f0 && push!(res, tz2)
    # @assert (length(res) == 2 || length(res) == 0) "This cylinder does not have an entrance and an exit!\n\tres: $res\n\ttz1: $tz1\ttz2: $tz2"
    proj_height1, proj_height2 = (ox .+ (tz1, tz2) .* dx) .^ 2 .+ (oy .+ (tz1, tz2) .* dy) .^ 2 .|> x -> (x < 1f0 || x ≈ 1f0)
    if length(res) == 1
        if proj_height1
            push!(res, tz1)
        else
            @assert proj_height2
            push!(res, tz2)
        end
    elseif proj_height1 && proj_height2
        res = [tz1, tz2]
    end
    # @assert (length(res) == 2 || length(res) == 0) "This cylinder does not have an entrance and an exit!\n\tres: $res\n\ttz1: $tz1\ttz2: $tz2"
    return res
end

function get_normal(::Type{Cylinder}, point::Point, ray::Ray)
    z = point.v[3]
    # if it comes from the caps
    (abs(z) ≈ 0.5f0) && return -sign(z) * sign(ray.dir.z) * NORMAL_Z

    # if it comes from the side
    x, y = normalize(point.v[1:2])
    normal = Normal{true}(x, y, 0f0)
    (normal ⋅ ray.dir < 0) ? normal : -normal
end

function get_uv(::Type{Cylinder}, point::Point)
    x, y, z = point.v
    z ≈  0.5f0 && return Vec2D(3 - 2x, 3 - 2y) * 0.25f0
    z ≈ -0.5f0 && return Vec2D(3 - 2x, 1 + 2y) * 0.25f0
    (clamp(z + 0.5f0, 0, 1), (atan(y, x)/2π + 1) * 0.5f0)
end


text:
using Gurobi

using Compat
using Compat.Test, Compat.SparseArrays, Compat.Random

@testset "C API" begin
    include("c_wrapper.jl")
end

@testset "MathProgBase Tests" begin
    @testset for file in ["env", "large_coefficients", "range_constraints",]
        evalfile(joinpath("MathProgBase", "$(file).jl"))
    end
    @testset "MathProgJuMP" begin
        evalfile(joinpath("MathProgBase", "mathprog.jl"))
    end
end

@testset "MathOptInterface Tests" begin
    include("MOI_wrapper.jl")
end


text:
# SCIP enums
{% for enum_name, enum_vals in parser.enums.items() %}typealias {{ enum_name }} Int8
{% for val_name, val_init in enum_vals.items() %}const {{ val_name }} = int8({{ val_init }})
{% endfor %}
{% endfor %}


text:
# only Convex.jl tests
# using Convex, ProxSDP, Test
# using Convex.ProblemDepot: run_tests
# @testset "Convex Problem Depot tests" begin
#     run_tests(;  exclude=[r"mip", r"exp"]) do problem
#         solve!(problem, () -> ProxSDP.Optimizer(
#             log_freq = 1_000_000, log_verbose = true,
#             tol_gap = 5e-8, tol_feasibility = 1e-7,
#             max_iter = 10_000_000, time_limit = 30.)
#             )
#     end
# end

# Convex.jl and SumOfSquares.jl tests
using ConvexTests, ProxSDP
@info "Starting ProxSDP tests"
do_tests("ProxSDP", () -> ProxSDP.Optimizer(
    log_freq = 1_000_000, log_verbose = true,
    tol_gap = 5e-8, tol_feasibility = 1e-7,
    max_iter = 100_000_000, time_limit = 4 * 30.
); exclude = [r"mip", r"exp"])


text:
"""
Multivariate Gaussian distribution, pdf
X[i,:] is i'th case
"""
function normpdf{T<:Real,U<:Real,V<:Real}(μ::AbstractArray{U,1}, Σ::AbstractArray{V,2}, X::AbstractArray{T,2})
    @assert size(Σ, 1)==size(Σ, 2) string("Σ should be a square matrix, but is a ", size(Σ, 1), "x", size(Σ, 2), " matrix")
    d = size(Σ, 2)

    @assert size(X, 2)==d string("X μst be a NxD matrix, but is Nx", size(X, 2))
    @assert length(μ)==d string("μ μst be a D=", d, " but is a ", length(μ), " matrix")

    Xc = broadcast(-, μ', X)
    logP = -0.5*sum((Xc/Σ).*Xc, 2)
    logP = reshape(logP, size(logP, 1))
    logZ = (d/2)*log(2*pi) + 0.5*logdet(Σ)

    return exp(logP - logZ)
end


text:
# TODO: Do backup system
module RemAssistant

import InteractiveUtils
import InteractiveUtils: clipboard

import TOML
import BibTeX
import JSON
import ArgParse
import ZipFile

using OrderedCollections
using GitLinks
using Dates

using Reexport
@reexport using StringRepFilter

include("core/Ref.jl")
include("core/bib2rem.jl")
include("core/bibref.jl")
include("core/configfile.jl")
include("core/knowledgebase.jl")
include("core/remapi.jl")
include("core/remdepot.jl")
include("core/remote_backup.jl")
include("core/utils.jl")

export config_file

include("api/bibrefs.jl")
include("api/update_remote_backup.jl")

include("cli/cli_os.jl")
include("cli/copy_bibref.jl")
include("cli/filter_bibref.jl")
include("cli/update_remote_backup.jl")

export bibrefs
export filter_bibref, copy_bibref

end # module


text:
using PyPlot
using LatinHypercubeSampling
using LinearAlgebra
using Random
using Base.Sort

function get_range(ary::Array{Float64,1})
    return minimum(ary), maximum(ary)
end

function get_range(ary::Array{Float64,2})
    ranges = Array{Tuple}(undef,size(ary,2))
    for i = 1:size(ary,2)
        ranges[i] = get_range(ary[:,i])
    end
    return ranges
end

function make_spatial_bins(min::Float64,max::Float64,bins::Int64)
    return collect(range(min,stop=max,length=bins))
end

function make_spatial_bins(ranges::Array{Tuple};bins::Int64=10)
    bin_ary = Array{Array{Float64}}(undef,length(ranges))
    for i = 1:length(bin_ary)
        bin_ary[i] = make_spatial_bins(ranges[i][1],ranges[i][2],bins+1)
    end
    return bin_ary
end

function return_bin_center(bins::Array{Array{Float64}})
    bin_c = Array{Float64}(undef,length(bins[1])-1,length(bins))
    for i = 1:length(bins)
        bin_c[:,i] = diff(bins[i])./2.0 .+ bins[i][1:end-1]
    end
    return bin_c
end

function vnormalize(X::Array;dims::Int64=1)
    minx = minimum(X,dims=dims)
    maxx = maximum(X,dims=dims)

    minx_mat = repeat(minx,size(X,1),size(X,2))
    maxx_mat = repeat(maxx,size(X,1),size(X,2))

    X_scaled = (X.-minx_mat)./(maxx_mat.-minx_mat)

    return X_scaled
end

function vnorm(A::Array;dims::Int64=1,normalize::Bool=false)
    if normalize
        A = vnormalize(A,dims=dims)
    end

    return sqrt.(sum(abs2,A,dims=dims))[:]
end

function set_distance(X::Array{Float64})
    sum_dists = Array{Float64}(undef,0)
    for i = 1:size(X,1)
        not_i = [x for x = 1:size(X,1) if x != i]
        dtr = distance_to_reference(X[i,:],X[not_i,:])
        push!(sum_dists,sum(dtr))
    end
    return sum_dists
end

function distance_to_reference(ref::Array{Float64},X::Array{Float64})
    X_scaled = ones(size(X))*Inf
    ref_scaled = ones(size(ref))*Inf
    for i = 1:size(X,2)
        xci = X[:,i]
        min_xci = minimum(xci)
        max_xci = maximum(xci)

        ref_scaled[i] = (ref[i].-min_xci)./(max_xci-min_xci)
        X_scaled[:,i] = (xci.-min_xci)./(max_xci-min_xci)
    end
    X_norms = vnorm(X_scaled,dims=2)
    return X_norms
end

function find_closest_neighbor(bin_center::Array{Float64},X::Array{Float64};exclude=Array{Float64}=[],deterministic::Bool=true)
    # scale so the two dimensions have equal weight
    X_norms = distance_to_reference(bin_center,X)
    ndist = abs.(X_norms .- vnorm(bin_center))

    if !isempty(exclude)
        idx = collect(1:length(ndist))
        idx_mask = [x for x in idx if x ∉ exclude]
        #  println("exclude: $exclude")
        #  println("idx: $idx")
        #  println("ndist: $ndist")
        #  println("idx_mask: $idx_mask")
        indexed_ndist = [ndist idx]
        sorted_ndist = sortslices(indexed_ndist[idx_mask,:],dims=1)
        if deterministic
            neighbor = sorted_ndist[1,2]
        else
            num_neighbors = Int(round(0.5 * size(sorted_ndist,1)))
            rand_neighbor = randperm(num_neighbors)[1]
            neighbor = sorted_ndist[rand_neighbor,2]
            #  @info num_neighbors rand_neighbor neighbor size(sorted_ndist,1)
            #  readline()
        end
        nn = Int(neighbor)
        #  neighbor = argmin(indexed_ndist[idx_mask,:])[1][1]
        #  nn = Int(indexed_ndist[idx_mask,:][neighbor,2])
        #  @info "closest neighbor info:" bin_center nn ndist[nn] ndist X[nn,:]
        #  readline()
    else
        nn = argmin(ndist)
        #  @info nn
        #  readline()
    end

    return nn
end

function nearest_neighbor_set(X::Array{Float64},lhc_centers::Array{Float64})
    sample_set = []
    for i = 1:size(lhc_centers,1)
        bin_center = lhc_centers[i,:]
        #  @info "lhc info" lhc_centers[i,:] bin_center
        # elements of each dimension in bins indicated by lhc_design
        # if a point doesn't exist there then find closest one
        closest_neighbor = find_closest_neighbor(bin_center,X,exclude=sample_set)
        push!(sample_set,closest_neighbor)
    end
    #  @info sample_set
    return sample_set
end

function random_neighbor_set(X::Array{Float64},lhc_centers::Array{Float64})
    sample_set = []
    for i = 1:size(lhc_centers,1)
        bin_center = lhc_centers[i,:]
        #  @info "lhc info" lhc_centers[i,:] bin_center
        # elements of each dimension in bins indicated by lhc_design
        # if a point doesn't exist there then find closest one
        closest_neighbor = find_closest_neighbor(bin_center,X,exclude=sample_set,deterministic=false)
        push!(sample_set,closest_neighbor)
    end
    #  @info sample_set
    return sample_set
end

function return_lhc_bin_centers(lhc::Array{Int64},bins::Array{Float64})
    lhc_centers = zeros(size(lhc))
    #  @info lhc bins
    for i = 1:size(lhc,1)
        for j = 1:size(lhc,2)
            lhc_bin = lhc[i,j]
            lhc_centers[i,j] = bins[lhc_bin,j]
        end
    end
    return lhc_centers
end

function max_rand_distance(X::Array,n::Int64)
    i = 0
    set_history = Array{Array{Int64}}(undef,0)
    dist_history = []
    dist = 0.0
    while i < 100000
        candidate_set = randperm(size(X,1))[1:n]

        if n == 2
            dists = vnorm(diff(X[candidate_set,:],dims=2),dims=1,normalize=true)
        else
            dists = set_distance(X[candidate_set,:])
            #  dists = distance_to_reference(X[candidate_set[1],:],X[candidate_set[2:end],:])
        end
        #  @info dists sum(dists) dist candidate_set
        #  readline()

        if sum(dists) > dist
            push!(dist_history,sum(dists))
            push!(set_history,candidate_set)
            dist = sum(dists)
        end

        i += 1
    end
    #  @info set_history[end] minimum(dist_history) dist_history[end]
    readline()
    return set_history[end]
end

function bin_sample(X::Array{Float64},n::Int64)
    # select n samples from X spread out 'spatially'
    @assert length(size(X)) == 2
    if size(X,1) <= n
        # need to return everything to have enough
        return X
    end

    # which bins do we want to take data from?
    if n == 1
        return rand(1:size(X,1),1)
    else
        lhc_design = randomLHC(n,size(X,2))
    end
    #  @info "LHC Size" size(lhc_design)

    r = get_range(X)
    b = make_spatial_bins(r,bins=n)
    bc = return_bin_center(b)
    lhc_centers = return_lhc_bin_centers(lhc_design,bc)

    #  sample_set = max_rand_distance(X,n)
    #  sample_set = nearest_neighbor_set(X,lhc_centers)
    sample_set = random_neighbor_set(X,lhc_centers)
    return sample_set
end

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
# import ../shared.styl
#=

use https://github.com/Keno/julia-wasm

=#

function fib(n)
	if n <= 1 return n end
	return fib(n - 1) + fib(n - 2)
end

fibprint(n) = println(fib(n))

A = collect(0:9)  # create an array with each integer in range
fibprint.(A)      # apply fibprint to each array item

nothing           # return nothing; supress extra output

max_stars_repo_name: fiugd/plugins

max_issues_repo_name: fiugd/plugins

max_forks_repo_name: fiugd/plugins


text:
module Each
using Base.Meta

export each, @each

"""
    each(Type, collection)

Return the collection unchanged.

    each(Type, element::Type)

Return the element wrapped in a tuple.

Useful for functions where the arguments can be a list or a single element which should be
iterated over, but the elements could be iterable themselves.
Extra care should be taken when the list and the elements are similar types: When working
with tuples of tuples of integers, use `each(Tuple{Vararg{Integer}},collection)` rather than
`each(Tuple,collection)` as the latter will always use the entire list as one element.

# Examples
```julia-repl
julia>  function f(a::Union{String,AbstractArray{String}})
            # This will fail for "scalar" Strings
            for s in a
                println(s)
            end
        end
julia>  f(["apples","oranges"])
apples
oranges
julia>  f("apples")
a
p
p
l
e
s
julia>  function f(a::Union{String,AbstractArray{String}})
            # This handles strings correctly
            for s in each(String,a)
                println(s)
            end
        end
julia>  f(["apples","oranges"])
apples
oranges
julia>  f("apples")
apples

julia>  function f(a::Tuple)
            for t in each(Tuple,a)
                println(t)
            end
        end
julia>  f( (1, 2) ) # argument is a tuple, wrap in tuple
(1, 2)
julia>  f( ( (1, 2), (3, 4) ) ) # same thing
((1, 2), (3, 4))
julia>  function f(a::Tuple)
            for t in each(Tuple{Vararg{Number}},a)
                println(t)
            end
        end
julia>  f( (1, 2) )
(1, 2)
julia>  f( ( (1,2), (3, 4) ) )
(1, 2)
(3, 4)
```
"""
function each(::Type{T}, collection) where T
    eltype(collection) <:T && return collection
    error(typeof(collection)," does not contain elements of type ",T)
end

each(::Type{T}, element::T) where T = (element,)

"""
## @each exp::Expr -> ::Bool
Determines whether each element of an iterable meets a certain condition. Returns a boolean, true if all of the elements meet a condition, false if otherwise.
Used in conditional contexts.
### example
```
x = [5, 10, 15, 20]
if @each x % 5 == 0
   println("They are!")
end

They are!

if @each x > 25
    println("They are!")
end

if @each x < 25
    println("They are!")
end

They are!
```
"""
macro each(exp::Expr)
    strexpr = string(exp)
    comp = split(strexpr, ' ')
    for value in eval(Meta.parse(string(comp[1])))
        to_parse = string(value, " ", comp[2], " ", comp[3])
        state = eval(Meta.parse(to_parse))
        if state != true
            return(false)
        end
    end
    return(true)
end

end # Module


text:
# D2Q21
# See lectures in LBM

# Same as D2Q17, but vertical and horizontal vectors are extended by one

struct D2Q21{
    T <: Real,
    Abscissaes <: AbstractMatrix{<:Integer},
    Weights <: AbstractVector{<:Real},
} <: Quadrature
    abscissae::Abscissaes
    weights::Weights
    speed_of_sound_squared::T
end
function D2Q21()
    cs = sqrt(2 / 3)
    w_0 = 91 / 324
    w_1 = 1 / 12
    w_2 = 2 / 27
    w_3 = 7 / 360
    w_4 = 1 / 432
    w_5 = 0 / 36
    w_6 = 1 / 1620

    # https://strathprints.strath.ac.uk/29195/1/Zhang_YH_Pure_Accuracy_analysis_of_high_order_lattice_Boltzmann_models_for_rarefied_gas_flows_21_Jan_2011.pdf
    w_0 = 91 / 324
    w_1 = 1 / 12
    w_2 = 2 / 27
    w_3 = 7 / 360
    w_4 = 1 / 432
    w_5 = 1 / 1620
    w_6 = 0
    cs = 3 / 2

    return D2Q21(
        [   #A #B          #C          #D          #E          #F
            0 1 -1 0 0 1 -1 1 -1 2 -2 0 0 2 -2 2 -2 3 -3 0 0 3 -3 -3 3
            0 0 0 1 -1 1 -1 -1 1 0 0 2 -2 2 -2 -2 2 0 0 3 -3 3 -3 3 -3
        ],
        [
            w_0,
            w_1,
            w_1,
            w_1,
            w_1,
            w_2,
            w_2,
            w_2,
            w_2,
            w_3,
            w_3,
            w_3,
            w_3,
            w_4,
            w_4,
            w_4,
            w_4,
            w_5,
            w_5,
            w_5,
            w_5,
            w_6,
            w_6,
            w_6,
            w_6,
        ],
        cs,
    )
end

order(q::D2Q21) = 7
Base.show(io::IO, q::D2Q21) = show(io, "D2Q21")
Base.string(q::D2Q21) = "D2Q21"


text:
# Autogenerated wrapper script for Modflow6_jll for i686-w64-mingw32-libgfortran5
export libmf6, mf5to6, mf6, zbud6

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Modflow6")
JLLWrappers.@declare_library_product(libmf6, "libmf6.dll")
JLLWrappers.@declare_executable_product(mf5to6)
JLLWrappers.@declare_executable_product(mf6)
JLLWrappers.@declare_executable_product(zbud6)
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libmf6,
        "bin\\libmf6.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        mf5to6,
        "bin\\mf5to6.exe",
    )

    JLLWrappers.@init_executable_product(
        mf6,
        "bin\\mf6.exe",
    )

    JLLWrappers.@init_executable_product(
        zbud6,
        "bin\\zbud6.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Documenter
makedocs(;
    sitename="Demo"
)


text:
import Base: factor, +, -, *, size, copy, norm

export MNLdata,
	negloglik, grad_negloglik,
	nucnorm, grad_nucnorm, min_lin_st_nucnorm, prox_nucnorm,
	RectangularParam, FactoredParam, rectangular_part, factor,
	initialize_dropping_convexity,
	+, *, -, size, copy

### Data

type MNLdata<:AbstractData
	nobs::Int
	nrow::Int
	ncol::Int
	sets::Array{Array{Int,1},1}
	rows::Array{Int,1}
	cols::Array{Int,1}
end

### Statistical parameter to be estimated

# The idea is that Theta::RectangularParam, U::FactoredParam,
# and [W1 Theta; Theta' W2] = U U' \succeq 0 at the solution
# The minimizer Theta of L(Theta; data) + lambda*\|Theta\|_*
# and the minimizer U of L((U U')[1:m, m+1:m+n]) + lambda/2*\|U\|_F^2
# satisfy (U U')[1:m, m+1:m+n] = Theta

typealias RectangularParam AbstractArray

type FactoredParam
	asarray
	nrow
	ncol
end
top(x::FactoredParam) = x.asarray[1:x.nrow,:]
bottom(x::FactoredParam) = x.asarray[x.nrow+1:x.nrow+x.ncol,:]
rectangular_part(x::FactoredParam) = top(x)*bottom(x)'
copy(x::FactoredParam) = FactoredParam(copy(x.asarray), x.nrow, x.ncol)
function factor(Theta::RectangularParam; k = min(size(Theta)...))
	m, n = size(Theta)
	u,s,v = svd(Theta, k)
	k = min(k, sum(s .> 1e-12))
	return FactoredParam([u[:,1:k]*diagm(sqrt(s[1:k]));
					   v[:,1:k]*diagm(sqrt(s[1:k]))],
					   m, n)
end

function +(x::FactoredParam, y::FactoredParam)
	(x.nrow==y.nrow && x.ncol==y.ncol) || error("Incompatible sizes")
	FactoredParam(x.asarray+y.asarray, x.nrow, x.ncol)
end
+(x::FactoredParam, y) = FactoredParam(x.asarray+y, x.nrow, x.ncol)
+(y, x::FactoredParam) = FactoredParam(x.asarray+y, x.nrow, x.ncol)
*(x::FactoredParam, y::FactoredParam) = NotImplementedError()
*(x::FactoredParam, y::Number) = FactoredParam(x.asarray*y, x.nrow, x.ncol)
*(y::Number, x::FactoredParam) = FactoredParam(x.asarray*y, x.nrow, x.ncol)
norm(x::FactoredParam) = norm(x.asarray)
vecnorm(x::FactoredParam) = vecnorm(x.asarray)
-(x::FactoredParam) = FactoredParam(-x.asarray, x.nrow, x.ncol)
-(x, y::FactoredParam) = x + -y
size(x::FactoredParam) = (x.nrow, x.ncol)

### The function we'll be minimizing

function negloglik(Theta::RectangularParam, data::MNLdata)
	l = 0
	for t=1:data.nobs
		it = data.rows[t]
		jt = data.cols[t]
		St = data.sets[t]
		invlik = 0 # inverse likelihood of observation t
		# computing soft max directly is numerically unstable
		# instead note logsumexp(a_j) = logsumexp(a_j - M) + M
		# and we'll pick a good big (but not too big) M
		M = 0 #Theta[it,jt] - minimum(Float64[Theta[it,j] for j in St])
		for j in St
			invlik += exp(Theta[it,jt] - Theta[it,j] - M)
		end
		l += log(invlik) + M
	end
	return l / data.nobs
end

function grad_negloglik(Theta::RectangularParam, data::MNLdata)
	DTheta = zeros(size(Theta))
	for t=1:data.nobs
		it = data.rows[t]
		jt = data.cols[t]
		St = data.sets[t]
		# Using some nice algebra, you can show
		DTheta[it,jt] += 1
		# and for j \in S, j \ne jt,
		# DTheta -= 1/sum_{j' \in S} exp(Theta[it,j] - Theta[it,j'])
		# it's ok if this over/underflows, I think:
		# the contribution of one observation to one entry of the gradient
		# is always between -1 and 0
		sumexp = sum(map(j->exp(-Theta[it,j]), St))
		for j in St
			DTheta[it,j] -= exp(-Theta[it,j])/sumexp
		end
	end
	return DTheta / data.nobs
end

function nucnorm(x::FactoredParam)
	return 1/2*vecnorm(x.asarray)^2
end

function grad_nucnorm(U::FactoredParam)
	return U
end

# prox_nucnorm(U, alpha) solves minimize_X (1/2||X||^2 + 1/(2*alpha)||X-U||^2)
function prox_nucnorm(U::FactoredParam, alpha::AbstractFloat)
	return 1/(1+alpha)*U
end

function nucnorm(Theta::RectangularParam)
	u,s,v = svd(Theta)
	return sum(s)
end

function grad_nucnorm(Theta::RectangularParam)
	u,s,v = svd(Theta)
	return u*spdiagm(sign(s))*v'
end

# prox_nucnorm(T, alpha) solves minimize_X (||X||_* + 1/(2*alpha)||X-T||^2)
# via soft-thresholding the singular values
function prox_nucnorm(Theta::RectangularParam,
	alpha::AbstractFloat;
	k = 10, # k estimates the number of singular values that will be greater than alpha
	kpp = 10) # how much to increment k if we're wrong
	u,s,v = svd(Theta, k)
	while s[end] > alpha
		k += kpp
		u,s,v = svd(Theta, k)
	end
	return u*spdiagm(max(s - alpha, 0))*v'
end

# solves min_{nucnorm(x)<=delta} G \dot x
function min_lin_st_nucnorm(G, delta)
	# solution to min_{nucnorm(x)<=delta} G \dot x
	u,s,v = svd(G, 1)
	x = -delta*u*v'
	return x
end

function negloglik(U::FactoredParam, data::MNLdata)
	return negloglik(rectangular_part(U), data)
end

function grad_negloglik(U::FactoredParam, data::MNLdata)
	G = grad_negloglik(rectangular_part(U), data)
	# gradient is [0_m G; G' 0_n] * U
	# return FactoredParam([G*top(U); G'*bottom(U)], size(U)...)
	return FactoredParam([G*bottom(U); G'*top(U)], size(U)...)
end

### The initialization (aka very rough algorithm)

# this function isn't type stable (so compiles poorly)
# but we're only calling it once so we don't care
function initialize_dropping_convexity(grad_objective, zeroparam, k::Int=0)
	grad_0 = grad_objective(zeroparam)
	zeroparam[1,1] = 1
	grad_1 = grad_objective(zeroparam)
	normalization = vecnorm(grad_0 - grad_1)
	u,s,v = svd(-grad_0, k)

	if k==0 # return initialization for Theta
		psd_grad_0 = u*spdiagm(max(s,0))*v'
		return psd_grad_0 / normalization
	else # return initialization for U
		# Theta = u[:,1:k]*spdiagm(max(s[1:k],0))*v[:,1:k]'
		U = [u[:,1:k]*spdiagm(sqrt(max(s[1:k],0)));
			 v[:,1:k]*spdiagm(sqrt(max(s[1:k],0)))]
		return FactoredParam(U / sqrt(normalization), size(zeroparam)...)
	end
end

max_stars_repo_name: madeleineudell/FirstOrderOptimization


text:
#base.jl: base type, constant & function definitions
#==Comments
Main algorithm always tries to keep existing files to avoid copying for no
reason. This also allows successive sync operations to progressively overwrite
corrupt files (there seems to be a file corruption issue).

Algorithm also renames old "destination folder" with temporary name, and
rebuilds a new folder (moving files into that folder). This is done in an
attempt to satisfy media players that sort by creation order instead of file
names.
===============================================================================#

#==Useful constants
===============================================================================#

#==Main types
===============================================================================#
struct Playlist
	origin::String
	filelist::Vector{String}
end
Playlist(origin::String) = Playlist(origin, [])

#==Accessors & helper functions
===============================================================================#
Base.length(p::Playlist) = length(p.filelist)
Base.push!(p::Playlist, args...) = push!(p.filelist, args...)
Base.iterate(p::Playlist, args...) = iterate(p.filelist, args...)
Base.in(v::String, p::Playlist) = Base.in(v, p.filelist)
tmp_folder_name(destfld::String) = destfld * "_AStmpfolder"

#==Read/write functions
===============================================================================#
function read_m3u(path::String)
	path = abspath(path)
	prefix = dirname(path)
	result = Playlist(path)

	open(path, "r") do f
		for line in eachline(f)
			line = strip(line)
			if length(line) > 0 && line[1] != '#'
				push!(result, joinpath(prefix, line))
			end
		end
	end

	return result
end

#==File/folder manipulations
===============================================================================#
function clean_playlistfolder(tgtlist::Playlist)
	destfld = tgtlist.origin
	n_delete = 0
	if !isdir(destfld); return 0; end
	filelist = readdir(destfld)
	for filename in filelist
		if !(basename(filename) in tgtlist)
			@info("Deleting \"$filename\"...")
			rm(joinpath(destfld, filename))
			n_delete += 1
		end
	end
	return n_delete
end

#==Main algorithm
===============================================================================#
#Generate filesystem-based playlist (names for the output files)
function gen_fsplaylist(src::Playlist, destfld::String)
	result = Playlist(destfld)
	for (i, srcpath) in enumerate(src.filelist)
		push!(result,  @sprintf("%03d-%s", i, cleannamme(basename(srcpath))))
	end
	return result
end

"""`synchronize(src::Playlist, destfld::String)`

Synchronize all files in playlist to the destination directory.
# Inputs
 - src: Source Playlist
 - destfld: Destination folder (Path of filesystem-based playlist)
"""
function synchronize(src::Playlist, destfld::String)
	#Select file matching algorithm:
	filesmatch = filesmatch_SHA
	#filesmatch = filesmatch_fast

	@info("Synchronizing file-based playlist:\n$destfld...")
	println()

	destfld = abspath(destfld)
	tmpfld = tmp_folder_name(destfld)

	#Delete any files in playlist folder that are not in playlist.
	tgtlist = gen_fsplaylist(src, destfld)
	n_delete = clean_playlistfolder(tgtlist)

	#Move old playlist folder to tmpfld
	if isdir(tmpfld)
		rm(tmpfld)
	end
	if isdir(destfld)
		mv(destfld, tmpfld)
	else
		mkpath(tmpfld)
	end
	try
		mkpath(destfld)
	catch #Throw more useful error:
		error("Could not create destination folder: $destfld.")
	end

	n_missing = 0
	n_sync = 0
	n_sync_failed = 0
	for (srcpath, destfile) in zip(src, tgtlist)
		@info("Synchronizing \"$destfile\"\nto $srcpath")
		tmppath = joinpath(tmpfld, destfile) #Full path
		destpath = joinpath(destfld, destfile) #Full path

		if (!isfile(srcpath))
			@warn("Missing file: $srcpath")
			n_missing += 1
			continue
		end

		if !filesmatch(srcpath, tmppath)
			cp(srcpath, tmppath, force=true) #Won't preserve modified timestamps.
			#run(`cp -p $srcpath $destpath`) #No longer preserves timestamps either.
		end

		mv(tmppath, destpath)
		if !filesmatch(srcpath, destpath)
			@warn("Synchronization failed.")
			n_sync_failed += 1
		else
			n_sync += 1
		end
	end

	rm(tmpfld)
	missing = length(tgtlist) - n_sync
	print("\nDone. ")
	print("$n_sync files synchronized")
	print(", $n_sync_failed failed, $n_missing missing")
	print(", $n_delete files deleted.\n")
end

#Last line

max_stars_repo_name: ma-laforge/AudioStick

max_issues_repo_name: ma-laforge/AudioStick

max_forks_repo_name: ma-laforge/AudioStick


text:
# Swap time and space dimensions in netcdfs with forcing data

path_ts = "/data04/jmg/fsm_simulations/netcdf/forcings_ts"

path_st = "/data04/jmg/fsm_simulations/netcdf/forcings_st"

files = readdir(path_ts)

for file in files

    println("Running $(file)")

    file_ts = joinpath(path_ts, file)

    file_st = joinpath(path_st, file)

    run(`ncpdq -a dim_time,dim_space $(file_ts) $(file_st)`)

end


text:
# Copyright (c) 2020: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# DynamicBoundspODEIneq.jl
# See https://github.com/PSORLab/DynamicBoundspODEIneq.jl
#############################################################################
# example/sample_trajectory.jl
# A simple example illustrating how the underlying DifferentialInequality
# integrator may be used.
#############################################################################

# Imports Sundials and DifferentialEquations to access the CVODE solvers
# Imports DynamicBoundspODEsIneq to use DifferentialInequality integrator
# Imports DynamicBoundsBase for API and problem description functions.
using DynamicBoundspODEsIneq, DynamicBoundsBase, Sundials, DifferentialEquations

# Create problem
x0(p) = [34.0; 20.0; 0.0; 0.0; 16.0; 0.0]
function f!(du, u, p, t)
    du[1] = -p[1]*u[1]*u[2] + p[2]*u[3] + p[6]*u[6]
    du[2] = -p[1]*u[1]*u[2] + p[2]*u[3] + p[3]*u[3]
    du[3] =  p[1]*u[1]*u[2] - p[2]*u[3] - p[3]*u[3]
    du[4] =  p[3]*u[3] - p[4]*u[4]*u[5] + p[5]*u[6]
    du[5] = -p[4]*u[4]*u[5] + p[5]*u[6] + p[6]*u[6]
    du[6] =  p[4]*u[4]*u[5] - p[5]*u[6] - p[6]*u[6]
    return
end

tspan = (0.0, 0.1)
pL = [0.1; 0.033; 16.0; 5.0; 0.5; 0.3]
pU = 10.0*pL

prob = DynamicBoundsBase.ODERelaxProb(f!, tspan, x0, pL, pU)
set!(prob, SupportSet([i for i in 0.0:0.005:0.1]))

# Define a polyhedral constraint on state variables.
# In this case it arise from the stiochiometry of the
# underlying reacting system.
A = [0.0 -1.0 -1.0  0.0  0.0  0.0;
     0.0  0.0  0.0  0.0 -1.0 -1.0;
     1.0 -1.0 0.0 1.0 -1.0 0.0]
b = [-20.0; -16.0; -2.0]
set!(prob, PolyhedralConstraint(A, b))

# Provides apriori state bounds
xL = zeros(6)
xU = [34.0; 20.0; 20.0; 34.0; 16.0; 16.0]
set!(prob, ConstantStateBounds(xL, xU))

# Creates the integrator specifying that relaxations and subgradients
# thereof should be computed.
integrator = DifferentialInequality(prob, calculate_relax = false,
                                    calculate_subgradient = false,
                                    relax_ode_integrator = CVODE_Adams(),
                                    local_ode_integrator = CVODE_Adams())

# defines a point in p at which to compute relaxations
ratio = rand(6)
pstar = pL.*ratio .+ pU.*(1.0 .- ratio)
setall!(integrator, ParameterValue(), pstar)

# Computes the relaxation and subgradients
relax!(integrator)

# Computes a trajectory of the pODEs at pstar
integrate!(integrator)


text:
#=##############################################################################################
Filename: solve_wing_CF_CM.jl
Author: Ryan Anderson
Contact: rymanderson@gmail.com
README: define an `Action` object to solve a CCBlade rotor
=###############################################################################################

"""
    solve_wing_CF_CM(aircraft, parameters, freestream, environment, step_range, stepi, step_symbol)

Action function.

# Arguments:

* `aircraft::Aircraft`: `Aircraft` system object
* `parameters<:Parameters`: `Parameters` struct
* `freestream::Freestream`: `Freestream` object
* `environment::Environment` `Environment` object
* `step_range::AbstractArray`: array of times for which the simulation is run
* `stepi::Int`: index of the current step
* `step_symbol::String`: defines the step, e.g. `alpha` or `time`

`parameters <: Parameters` requires the following elements:

* `wake_function::Vector{Function}`: [i]th element is a function vwake(X::Vector{Float64}) describing the wake induced velocity at `X` at the ith step
* `CFs::Array{Float64,2}` : [:,j]th element is the [CD,CY,CL] force coefficients of the aircraft at the jth step
* `CMs::Array{Float64,2}` : [:,j]th element is the [CMx,CMy,CMz] moment coefficients of the aircraft at the jth step

"""
function solve_wing_CF_CM(aircraft, parameters, freestream, environment, step_range, stepi, step_symbol)

    # interpret freestream
    vlmfreestream = VL.Freestream(freestream)

    # get reference
    reference = aircraft.wing_system.system.reference[1]

    # solve vortex lattice
    VL.steady_analysis!(aircraft.wing_system.system, aircraft.wing_system.system.surfaces, aircraft.wing_system.system.reference[1], vlmfreestream; symmetric=true, additional_velocity = parameters.wake_function[stepi])

    # extract forces and moments
    CF, CM = VL.body_forces(aircraft.wing_system.system; frame=VL.Wind())

    # store to `parameters`
    parameters.CFs[:,stepi] .= CF
    parameters.CMs[:,stepi] .= CM

    return false
end

"""
solve_wing_CF_CM(system, step_range)

Method returns initialized elements required for the `parameters <: Parameters` struct during simulation.

# Arguments:

* `aircraft::Aircraft` : aircraft system to be simulated
* `step_range::AbstractArray` : defines each step of the simulation

# Returns:

* `wake_function::Vector{Function}`: [i]th element is a function vwake(X::Vector{Float64}) describing the wake induced velocity at `X` at the ith step
* `CFs::Array{Float64,2}` : [:,j]th element is the [CD,CY,CL] force coefficients of the aircraft at the jth step
* `CMs::Array{Float64,2}` : [:,j]th element is the [CMx,CMy,CMz] moment coefficients of the aircraft at the jth step

"""
function solve_wing_CF_CM(aircraft, step_range)

    wake_function = Function[(x) -> [0.0, 0.0, 0.0] for i in 1:length(step_range)]
    CFs = zeros(3, length(step_range))
    CMs = zeros(3, length(step_range))

    return wake_function, CFs, CMs
end

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
import CLBLAS
using OpenCL: cl

const clblas = CLBLAS
clblas.setup()
alpha = cl.cl_float(5)
beta = cl.cl_float(2)

device, context, queue = clblas.get_next_compute_context()
A = clblas.randBuf(8, 12; ctx=context, queue=queue)
B = clblas.randBuf(12, 24; ctx=context, queue=queue)
C = clblas.randBuf(8, 24; ctx=context, queue=queue)

matA = clblas.fetch(A)
matB = clblas.fetch(B)
matC = clblas.fetch(C)

data = (5.0*matA*matB) + (2.0*matC)
future = clblas.clblasSgemm(UInt32(0), UInt32(0), alpha, A, B, beta, C)

result = fetch(future)
println(result)
println("------------")
println(data)
for i in 1:length(data)
  Test.@test_approx_eq_eps result[i] data[i]  0.00001
end

info("success!")
clblas.teardown()


text:
# using Base.Tests
using Test
using SemanticModels
using SemanticModels.Unitful
import SemanticModels.Unitful: DimensionError, uconvert, NoUnits, s
using DifferentialEquations
using Distributions: Uniform
using GLM
using DataFrames
using Plots

include("parse.jl")
include("cassette.jl")
include("transform/ode.jl")

examples = ["agentbased.jl",
            "agentgraft.jl",
            "modelmacro.jl",
            "monomial_regression.jl",
            "multivariate_regression.jl",
            "polynomial_regression.jl",
            "workflow.jl",
            "pseudo_polynomial_regression.jl",
            "odegraft.jl",
            ]
for ex in examples
    @info "Running example: " file=ex
    try
      include(ex)
    catch err
      println(err)
      @warn "Error running: " file=ex
    end
end


text:
module GPUChecker

# package code goes here
function gpu_usages()
	rawtext = readall(pipeline(`nvidia-smi`, `grep |`, `grep N/A`))
	str_usages = split(rawtext, '\n')
	usages = []
	for usg in str_usages
		if usg != ""
			usage, capacity = map(x->parse(Int, match(r"(\d+)", x).captures[1]),split(strip(split(usg, '|')[3]), '/'))
			push!(usages, (usage, capacity))
		end
	end
	return usages
end

function first_min_used_gpu()
	usages = gpu_usages()
	return indmin(usages) - 1
end

export gpu_usages
export first_min_used_gpu

end # module


text:
abstract type DAEInitializationAlgorithm end

struct BrownFullBasicInit{T} <: DAEInitializationAlgorithm
	abstol::T
end
BrownFullBasicInit() = BrownFullBasicInit(1e-10)

function initialize_dae!(integrator, u, du, differential_vars, alg::BrownFullBasicInit, ::Val{true})
	@unpack p, t, f = integrator

	tmp = get_tmp_cache(integrator)[1]
	f(tmp, du, u, p, t)

	if integrator.opts.internalnorm(tmp,t) <= alg.abstol
		return
	elseif differential_vars === nothing
		error("differential_vars must be set for DAE initialization to occur. Either set consistent initial conditions or differential_vars")
	end

	nlequation = (out, x) -> begin
		for i in 1:length(x)
			if differential_vars[i]
				du[i] = x[i]
			else
				u[i] = x[i]
			end
		end
		f(out, du, u, p, t)
	end

	r = nlsolve(nlequation, zero(u))

	for i in 1:length(u)
		if differential_vars[i]
			du[i] = r.zero[i]
		else
			u[i] = r.zero[i]
		end
	end

	return
end

function initialize_dae!(integrator, _u, _du, differential_vars, alg::BrownFullBasicInit, ::Val{false})
	@unpack p, t, f = integrator

	if integrator.opts.internalnorm(f(_du, _u, p, t),t) <= alg.abstol
		return
	elseif differential_vars === nothing
		error("differential_vars must be set for DAE initialization to occur. Either set consistent initial conditions or differential_vars")
	end

	if _u isa Number && _du isa Number
		# This doesn't fix static arrays!
		u = [_u]
		du = [_du]
	else
		u = _u
		du = _du
	end

	nlequation = (out,x) -> begin
		for i in 1:length(x)
			if differential_vars[i]
				du[i] = x[i]
			else
				u[i] = x[i]
			end
		end
		out .= f(du, u, p, t)
	end

	r = nlsolve(nlequation, zero(u))

	for i in 1:length(u)
		if differential_vars[i]
			du[i] = r.zero[i]
		else
			u[i] = r.zero[i]
		end
	end

	if _u isa Number && _du isa Number
		# This doesn't fix static arrays!
		integrator.u = first(u)
		integrator.cache.du = first(du)
	else
		integrator.u = u
		integrator.cache.du = du
	end

	return
end


text:
abstract type MessagePassing <: Meta end

adjlist(m::T) where {T<:MessagePassing} = m.adjlist
message(m::T; kwargs...) where {T<:MessagePassing} = identity(; kwargs...)
update(m::T; kwargs...) where {T<:MessagePassing} = identity(; kwargs...)

function update_edge(m::T; gi::GraphInfo, kwargs...) where {T<:MessagePassing}
    adj = gi.adj
    edge_idx = gi.edge_idx
    M = message(m; adjacent_vertices_data(kwargs, 1, adj[1])...,
                   incident_edges_data(kwargs, 1, adj[1])...)
    Y = similar(M, size(M, 1), gi.E)
    Y[:, 1:edge_idx[2]] = M
    @inbounds Threads.@threads for i = 2:gi.V
        j = edge_idx[i]
        k = edge_idx[i+1]
        Y[:, j+1:k] = message(m; adjacent_vertices_data(kwargs, i, adj[i])...,
                                 incident_edges_data(kwargs, i, adj[i])...)
    end
    Y
end

update_vertex(m::T; kwargs...) where {T<:MessagePassing} = update(m; kwargs...)

aggregate_neighbors(m::T, aggr::Symbol; kwargs...) where {T<:MessagePassing} =
    pool(aggr, kwargs[:cluster], kwargs[:M])

function propagate(mp::T; aggr::Symbol=:add, adjl=adjlist(mp), kwargs...) where {T<:MessagePassing}
    gi = GraphInfo(adjl)

    # message function
    M = update_edge(mp; gi=gi, kwargs...)

    # aggregate function
    M = aggregate_neighbors(mp, aggr; M=M, cluster=generate_cluster(M, gi))

    # update function
    Y = update_vertex(mp; M=M, kwargs...)
    return Y
end


text:
module ClozeTest
using Fetch
using DataDeps

using ..Datasets: Dataset
import ..Datasets: testfile, trainfile, get_labels

export StoryCloze

function __init__()
    storycloze_init()
end

include("./storycloze.jl")

end


text:
struct EstSamplePath{T,R<:AbstractVector,S}
    tt::R
    yy::Vector{T}
    vv::Vector{S}
end

function piecewise(X::EstSamplePath, tend = X.tt[end])
    EstSamplePath(piecewise(SamplePath(X.tt, X.yy), tend)...,
        piecewise(SamplePath(X.tt, X.vv), tend)[2])
end

struct SamplePathBand{T,R<:AbstractVector}
    tt::R
    ll::Vector{T}
    uu::Vector{T}
end

function piecewise(X::SamplePathBand, tend = X.tt[end])
    SamplePathBand(piecewise(SamplePath(X.tt, X.ll), tend)...,
        piecewise(SamplePath(X.tt, X.uu), tend)[2])
end


text:
using GenericInstruments
GI = GenericInstruments
####################################
resmgr = GI.viOpenDefaultRM()
# Instantiate obj
dmm1 = GI.DMM.INSTR(:KE2000, "GPIB0::4::INSTR")
GI.DMM.connect!(resmgr, dmm1)
# Reset
GI.DMM.reset_instr(dmm1)
# Clear
GI.DMM.clear_instr(dmm1)
# IDN
GI.DMM.get_idn(dmm1)
# Function
GI.DMM.set_fc(dmm1,fc="VOLT:DC")
# Function and range
GI.DMM.set_fc_range(dmm1,fc="VOLT",rn=10)
GI.DMM.set_fc_range_auto(dmm1,fc="VOLT")
# Function and NPLC
GI.DMM.set_fc_nplc(dmm1,fc="VOLT",nplc=2)
# Init once
GI.DMM.init_meas(dmm1)
# Init continuous
GI.DMM.init_meas_continuous(dmm1,st="ON")
# Abort
#GI.DMM.abort_meas(dmm1)
# Read measurement
y = GI.DMM.read_meas(dmm1)
# Disconnect
GI.DMM.disconnect!(dmm1)


text:
#sigma x (NOT gate)
function sigma_x(o::Op, mode::Int64) #c is control, t is target

    base = basis(o)
    d = dim(o)
    l = 2^d

    if mode > d
        throw(ArgumentError("Your mode must be within your dimensions!"))
    end

    row = spzeros(l)
    col = spzeros(l)
    data = spzeros(l)
    for k in 1:l
        row[k] = k
        if base[k,mode] != 0
            col[k] = k - 2^(d-mode)
            data[k] = 1.0
        else
            col[k] = k + 2^(d-mode)
            data[k] = 1.0
        end
    end
    sigma_x = sparse(row, col, data)
    return sigma_x
end

#sigma y
function sigma_y(o::Op, mode::Int64) #c is control, t is target

    base = basis(o)
    d = dim(o)
    l = 2^d

    if mode > d
        throw(ArgumentError("Your mode must be within your dimensions!"))
    end

    row = spzeros(l)
    col = spzeros(l)
    data = spzeros(Complex{Float64},l)
    for k in 1:l
        row[k] = k
        if base[k,mode] != 0
            col[k] = k - 2^(d-mode)
            data[k] = 1.0*im
        else
            col[k] = k + 2^(d-mode)
            data[k] = -1.0*im
        end
    end
    sigma_y = sparse(row, col, data)
    return sigma_y
end

#sigma z
function sigma_z(o::Op, mode::Int64) #c is control, t is target

    base = basis(o)
    d = dim(o)
    l = 2^d

    if mode > d
        throw(ArgumentError("Your mode must be within your dimensions!"))
    end

    row = spzeros(l)
    col = spzeros(l)
    data = spzeros(l)
    for k in 1:l
        row[k] = k
        if base[k,mode] != 0
            col[k] = k
            data[k] = -1.0
        else
            col[k] = k
            data[k] = 1.0
        end
    end
    sigma_z = sparse(row, col, data)
    return sigma_z
end

#phase shift
function phase(o::Op, mode::Int64, phi::Real) #c is control, t is target

    base = basis(o)
    d = dim(o)
    l = 2^d

    if mode > d
        throw(ArgumentError("Your mode must be within your dimensions!"))
    end

    row = spzeros(l)
    col = spzeros(l)
    data = spzeros(Complex{Float64},l)
    for k in 1:l
        row[k] = k
        col[k] = k
        if base[k,mode] != 0
            data[k] = round(exp(im*phi), digits = 15)
        else
            data[k] = 1.0
        end
    end
    phase = sparse(row, col, data)
    return phase
end

#Hadamard
function hadamard(o::Op, mode1::Int64, mode2::Int64)
    #we first check that the input modes are different
    if mode1 == mode2
        throw(ArgumentError("Modes must be differents"))
    end

    i = min(mode1, mode2)
    j = max(mode1, mode2)
    base = basis(o)
    d = dim(o)
    l = 2^d

    if mode1 > d || mode2 > d
        throw(ArgumentError("Your modes must be within your dimensions!"))
    end

    row = spzeros(l+2^(d-1))
    col = spzeros(l+2^(d-1))
    data = spzeros(l+2^(d-1))
    counter = 1
    for k in 1:l
        if base[k,i] != 0
            if base[k,j] != 0
                row[counter] = k
                col[counter] = k
                data[counter] = 1.0
                counter = counter + 1
            else
                row[counter] = k
                col[counter] = k
                data[counter] = 1.0/sqrt(2)
                counter = counter + 1
                row[counter] = k
                col[counter] = k - 2^(d-i) + 2^(d-j)
                data[counter] = 1.0/sqrt(2)
                counter = counter + 1
            end

        elseif base[k,j] != 0
            row[counter] = k
            col[counter] = k
            data[counter] = 1.0/sqrt(2)
            counter = counter + 1
            row[counter] = k
            col[counter] = k + 2^(d-i) - 2^(d-j)
            data[counter] = 1.0/sqrt(2)
            counter = counter + 1
        else
            row[counter] = k
            col[counter] = k
            data[counter] = 1.0
            counter = counter + 1
        end
    end
    had = sparse(row,col,data)
    return had
end

#U control not
function ucnot(o::Op, control::Int64, target::Int64) #c is control, t is target

    c = control
    t = target
    base = basis(o)
    d = dim(o)

    #we first check basic conditions
    if c > d || t > d
        throw(ArgumentError("Your modes must be within your dimensions!"))
    end

    if c == t
        throw(ArgumentError("Control must be different from target"))
    end

    l = 2^d
    row = spzeros(l)
    col = spzeros(l)
    data = spzeros(l)
    for k in 1:l
        row[k] = k
        if base[k,c] != 0
            if base[k,t] != 0
                col[k] = k - 2^(d-t)
                data[k] = 1.0
            else
                col[k] = k + 2^(d-t)
                data[k] = 1.0
            end
        else
            col[k] = k
            data[k] = 1.0
        end
    end
    u_cnot = sparse(row, col, data)
    return u_cnot
end

#Swap:
function swap(o::Op, mode1::Int64, mode2::Int64) #c is control, t is target

    base = basis(o)
    d = dim(o)
    l = 2^d

    if mode1 > d || mode2 > d
        throw(ArgumentError("Your modes must be within your dimensions!"))
    end

    if mode1 == mode2
        throw(ArgumentError("Modes must be differents"))
    end

    row = spzeros(l)
    col = spzeros(l)
    data = spzeros(l)
    for k in 1:l
        row[k] = k
        if base[k, mode1] != 0
            if base[k, mode2] != 0
                col[k] = k
                data[k] = -1.0
            else
                col[k] = k - 2^(d-mode1) + 2^(d-mode2)
                data[k] = 1.0
            end
        else
            if base[k, mode2] != 0
                col[k] = k + 2^(d-mode1) - 2^(d-mode2)
                data[k] = 1.0
            else
                col[k] = k
                data[k] = 1.0
            end
        end
    end
    swap = sparse(row, col, data)
    return swap
end


text:
# --------------------------------------------------------------------------
# ACE1.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

module PairPotentials

import JuLIP
using JuLIP: JVec, JMat, Atoms, AtomicNumber
using JuLIP.MLIPs: IPBasis
using LinearAlgebra: norm, dot
using JuLIP.Potentials: ZList, SZList, @pot, @D,
                        PairPotential, SimplePairPotential
using StaticArrays: SMatrix

using ACE1: ScalarBasis, _allfieldsequal

import JuLIP: evaluate!, evaluate_d!, cutoff,
              evaluate, evaluate_d,
              read_dict, write_dict,
              energy, forces, virial,
              alloc_temp, alloc_temp_d,
              z2i, i2z, numz,
              fltype, rfltype

import JuLIP.Potentials: zlist

import ACE1: scaling

import JuLIP.MLIPs: alloc_B, alloc_dB

import Base: ==, length

include("pair_basis.jl")

include("pair_pot.jl")

include("repulsion.jl")

end


text:
module ITensorLocalMPO

using
  ITensors,
  LinearAlgebra

export
  combine_and_transform,
  inv_transform_and_uncombine,
  scale_bases

include("scale_bases.jl")

end

max_stars_repo_name: mtfishman/ITensorLocalMPO

max_issues_repo_name: mtfishman/ITensorLocalMPO

max_forks_repo_name: mtfishman/ITensorLocalMPO


text:
function SeisDecimate(in;mode="random",perc=50,incx1=1,incx2=1,incx3=1,incx4=1)

    if (mode=="random") # decimate data randomly
        out = copy(in)
	mask = rand(1,size(in,2)*size(in,3)*size(in,4)*size(in,5));
	mask[find(mask .< perc/100)] = 0;
	mask[find(mask .>= perc/100)] = 1;
	for it = 1 : size(in,1)
	    out[it,:] = out[it:it,:].*mask;
	end
    else # decimate data regularly with respect to 4 spatial dimensions
        out = zeros(in)
	if (size(in),5) > 1
	    out[:,1:incx1:end,1:incx2:end,1:incx3:end,1:incx4:end] = in[:,1:incx1:end,1:incx2:end,1:incx3:end,1:incx4:end]
	elseif (size(in),4) > 1
            out[:,1:incx1:end,1:incx2:end,1:incx3:end] = in[:,1:incx1:end,1:incx2:end,1:incx3:end]
        elseif (size(in),3) > 1
            out[:,1:incx1:end,1:incx2:end] = in[:,1:incx1:end,1:incx2:end]
        elseif (size(in),2) > 1
            out[:,1:incx1:end] = in[:,1:incx1:end]
        end

    end

    return out
end


text:
using FountainCodes, Test, Random
println("Starting tests")
@time @testset "GF256" begin include("GF256_test.jl") end
@time @testset "Gray sequence" begin include("Gray_test.jl") end
@time @testset "QMatrix" begin include("QMatrix_test.jl") end
@time @testset "Decoder" begin include("Decoder_test.jl") end
@time @testset "R10" begin include("R10_test.jl") end
@time @testset "RQ" begin include("RQ_test.jl") end
@time @testset "LT" begin include("LT_test.jl") end
@time @testset "Floating-point" begin include("Float_test.jl") end

max_stars_repo_name: severinson/RaptorCodes

max_issues_repo_name: severinson/RaptorCodes

max_forks_repo_name: severinson/RaptorCodes


text:
using MCVI
using Test
import MCVI: init_lower_action, lower_bound, upper_bound
using POMDPs
using Random
using POMDPModelTools
using POMDPSimulators
using POMDPModels # for LightDark1d
using POMDPLinter: @requirements_info, @show_requirements

# Bounds
mutable struct LightDark1DLowerBound
    rng::AbstractRNG
end

mutable struct LightDark1DUpperBound
    rng::AbstractRNG
end

function lower_bound(lb::LightDark1DLowerBound, p::LightDark1D, s::LightDark1DState)
    r = @gen(:r)(p, s, init_lower_action(p), lb.rng)
    return r * discount(p)
end

function upper_bound(ub::LightDark1DUpperBound, p::LightDark1D, s::LightDark1DState)
    steps = abs(s.y)/p.step_size + 1
    return p.correct_r*(discount(p)^steps)
end

function init_lower_action(p::LightDark1D)
    return 0 # Worst? This depends on the initial state? XXX
end

include("test_policy.jl")

include("test_updater.jl")

include("test_belief.jl")

@test test_dummy_graph()
@test test_dummy_graph2()

@test test_backup()

@test test_belief()

include("test_solve.jl")
@test test_solve()

include("test_simulation.jl")
test_simulation()

include("test_requirements.jl")
test_requirements()

# @test test_dummy_heuristics()


text:
@testset "PoissonLikelihood" begin
    lik = PoissonLikelihood()
    test_interface(lik, SqExponentialKernel(), rand(10))
end


text:
export play

using Colors

# coordinate transform for Makie.jl
transform(x::Int) = p -> CartesianIndex(p[2], x-p[1]+1)

using .Makie

function init_screen(env::Observable{<:AbstractGridWorld}; resolution=(1000,1000))
    scene = Scene(resolution = resolution, raw = true, camera = campixel!)

    area = scene.px_area
    grid_size = size(env[].world)[2:3]
    grid_inds = CartesianIndices(grid_size)
    tile_size = @lift((widths($area)[1] / size($env.world, 2), widths($area)[2] / size($env.world, 3)))
    T = transform(size(env[].world, 2))
    boxes(pos, s) = [FRect2D((T(p).I .- (1,1)) .* s, s) for p in pos]
    centers(pos, s) = [(T(p).I .- (0.5,0.5)) .* s for p in pos]

    # 1. paint background
    poly!(scene, area)

    # 2. paint each kind of object
    for o in get_object(env[])
        if o !== EMPTY
            scatter!(scene, @lift(centers(findall($env.world[o,:,:]), $tile_size)), color=get_color(o), marker=convert(Char, o), markersize=@lift(minimum($tile_size)))
        end
    end

    # 3. paint stroke
    poly!(scene, @lift(boxes(vec(grid_inds), $tile_size)), color=:transparent, strokecolor = :lightgray, strokewidth = 4)

    # 3. paint agent's view
    view_boxes = @lift boxes([p for p in get_agent_view_inds($env) if p ∈ grid_inds], $tile_size)
    poly!(scene, view_boxes, color="rgba(255,255,255,0.2)")

    # 4. paint agent
    agent = @lift(get_agent($env))
    agent_position = @lift((T(get_agent_pos($env)).I .- (0.5, 0.5)).* $tile_size)
    scatter!(scene, agent_position, color=@lift(get_color($agent)), marker=@lift(convert(Char, $agent)), markersize=@lift(minimum($tile_size)))

    display(scene)
    scene
end

function play(env::AbstractGridWorld;file_name=nothing,frame_rate=24)
    print("""
    Key bindings:
    ←: TurnLeft
    →: TurnRight
    ↑: MoveForward
    q: Quit
    """)
    env_node = Node(env)
    scene = init_screen(env_node)
    is_quit = Ref(false)

    if !isnothing(file_name)
        vs = VideoStream(scene)
        recordframe!(vs)
    end

    on(scene.events.keyboardbuttons) do b
        if ispressed(b, Keyboard.left)
            env(TURN_LEFT)
            env_node[] = env
        elseif ispressed(b, Keyboard.right)
            env(TURN_RIGHT)
            env_node[] = env
        elseif ispressed(b, Keyboard.up)
            env(MOVE_FORWARD)
            env_node[] = env
        elseif ispressed(b, Keyboard.q)
            is_quit[] = true
        end
        isnothing(file_name) || recordframe!(vs)
    end

    try
        while !is_quit[]
            sleep(0.5)
        end
    catch
    end
    isnothing(file_name) || save(file_name, vs;framerate=frame_rate)
    Makie.GLMakie.destroy!(Makie.GLMakie.global_gl_screen())
end


text:
module SparseMatrixCSCTests

    using SparseArrays
    using LinearAlgebra
    using Gridap.Algebra
    using Test

    maxnz=5
    maxrows=5
    maxcols=5

    for T in (Int32,Int,Float32,Float64)
        I = Vector{Int}()
        J = Vector{Int}()
        V = Vector{T}()
        for (ik, jk, vk) in zip(rand(1:maxrows, maxnz), rand(1:maxcols, maxnz), rand(1:T(maxnz), maxnz-1))
            push_coo!(I,J,V,ik,jk,vk)
        end
        push_coo!(I,J,V,maxrows,maxcols,maxnz)
        finalize_coo!(I,J,V,maxcols,maxrows)
        CSC = sparse(I, J, V, maxcols,maxrows)

        @test is_entry_stored(SparseMatrixCSC,1,1)
        @test is_entry_stored(SparseMatrixCSC,1,2)
        @test is_entry_stored(SparseMatrixCSC,2,1)

        _l = 10
        _I, _J, _V = allocate_coo_vectors(SparseMatrixCSC{T,Int},_l)
        @test length(_I) == _l
        @test length(_J) == _l
        @test length(_V) == _l
        @test eltype(_I) == Int
        @test eltype(_J) == Int
        @test eltype(_V) == T

        V2 = copy(V)
        V2 .= 0
        CSC2 = sparse(I, J, V2, maxcols,maxrows)

        copy_entries!(CSC2,CSC)
        @test CSC2 == CSC
        copy_entries!(CSC2,CSC2)
        @test CSC2 == CSC

        @test convert(SparseMatrixCSC, CSC) === CSC

        @test hasrowmajororder(CSC)  == false
        @test hascolmajororder(CSC) == true
        @test getptr(CSC)           == CSC.colptr
        @test getindices(CSC)       == rowvals(CSC)

        for col = 1:maxcols
            for j in nzrange(CSC, col)
                @test nonzeros(CSC)[j] == CSC[rowvals(CSC)[j], col]
            end
        end

        @test size(CSC) == (maxrows,maxcols)

        @test nnz(CSC) == count(!iszero, CSC) <= maxnz

        vold = getindex(CSC,maxrows,maxcols)
        add_entry!(CSC,1,maxrows,maxcols,+)
        @test getindex(CSC,maxrows,maxcols) == vold+1

        fill_entries!(CSC,0)
        @test all(x->x==0, nonzeros(CSC))

    end

end


text:
using ReadableNumbers
using Test

setprecision(BigFloat,128);

@test readable(BigFloat(pi)) == "3.14159_26535_89793_23846_26433_83279_50288_4195"
@test readable(factorial(BigInt(34))) == "295,232,799,039,604,140,847,618,609,643,520,000,000"
@test BigFloat(pi) == parse_readable(BigFloat, readable(BigFloat(pi)), '_')
@test factorial(BigInt(34)) == parse_readable(BigInt, readable(factorial(BigInt(34))), ',')


text:
cd(@__DIR__)
allowed = ["case", "space", "flux", "collision"]
D = KitBase.read_dict("config.txt", allowed)

D = KitBase.read_dict("config.txt")
ks, ctr, face, simTime = KitBase.initialize(D)

ks, ctr, face, simTime = KitBase.initialize("config.txt")
KitBase.write_sol(ks, ctr)
KitBase.write_sol(ks, ctr; mode = :jld)

KitBase.plot_line(ks, ctr)
plot(ks, ctr, legend = :none)

# tecplot writer
x = rand(5)
sol = randn(5, 2)
write_tec(x, sol)

x = zeros(4, 3)
y = zero(x)
for i in axes(x, 1), j in axes(y, 2)
    x[i, j] = 0.25 * (i - 1)
    y[i, j] = 0.5 * (j - 1)
end
sol = randn((size(x)..., 2))
write_tec(x, y, sol)


text:
"""
This struct acts as an infinity bus.
"""
struct StaticSource <: PSY.Injection
    number::Int64
    name::Symbol
    bus::PSY.Bus
    V_R::Float64
    V_I::Float64
    X_th::Float64
end