text:
# core formatting functions

### auxiliary functions

function _repwrite(out::IO, c::Char, n::Int)
    while n > 0
        write(out, c)
        n -= 1
    end
end

### print string or char

@compat function _pfmt_s(out::IO, fs::FormatSpec, s::Union{AbstractString,Char})
    wid = fs.width
    slen = length(s)
    if wid <= slen
        write(out, s)
    else
        a = fs.align
        if a == '<'
            write(out, s)
            _repwrite(out, fs.fill, wid-slen)
        else
            _repwrite(out, fs.fill, wid-slen)
            write(out, s)
        end
    end
end

### print integers

_mul(x::Integer, ::_Dec) = x * 10
_mul(x::Integer, ::_Bin) = x << 1
_mul(x::Integer, ::_Oct) = x << 3
@compat _mul(x::Integer, ::Union{_Hex, _HEX}) = x << 4

_div(x::Integer, ::_Dec) = div(x, 10)
_div(x::Integer, ::_Bin) = x >> 1
_div(x::Integer, ::_Oct) = x >> 3
@compat _div(x::Integer, ::Union{_Hex, _HEX}) = x >> 4

function _ndigits(x::Integer, op)  # suppose x is non-negative
    m = 1
    q = _div(x, op)
    while q > 0
        m += 1
        q = _div(q, op)
    end
    return m
end

_ipre(op) = ""
@compat _ipre(::Union{_Hex, _HEX}) = "0x"
_ipre(::_Oct) = "0o"
_ipre(::_Bin) = "0b"

_digitchar(x::Integer, ::_Bin) = @compat Char(x == 0 ? '0' : '1')
_digitchar(x::Integer, ::_Dec) = @compat Char('0' + x)
_digitchar(x::Integer, ::_Oct) = @compat Char('0' + x)
_digitchar(x::Integer, ::_Hex) = @compat Char(x < 10 ? '0' + x : 'a' + (x - 10))
_digitchar(x::Integer, ::_HEX) = @compat Char(x < 10 ? '0' + x : 'A' + (x - 10))

_signchar(x::Number, s::Char) = x < 0 ? '-' :
                                s == '+' ? '+' :
                                s == ' ' ? ' ' : '\0'

function _pfmt_int{Op}(out::IO, sch::Char, ip::Compat.ASCIIString, zs::Integer, ax::Integer, op::Op)
    # print sign
    if sch != '\0'
        write(out, sch)
    end
    # print prefix
    if !isempty(ip)
        write(out, ip)
    end
    # print padding zeros
    if zs > 0
        _repwrite(out, '0', zs)
    end
    # print actual digits
    if ax == 0
        write(out, '0')
    else
        _pfmt_intdigits(out, ax, op)
    end
end

function _pfmt_intdigits{Op,T<:Integer}(out::IO, ax::T, op::Op)
    b_lb = _div(ax, op)
    b = one(T)
    while b <= b_lb
        b = _mul(b, op)
    end
    r = ax
    while b > 0
        (q, r) = divrem(r, b)
        write(out, _digitchar(q, op))
        b = _div(b, op)
    end
end

function _pfmt_i{Op}(out::IO, fs::FormatSpec, x::Integer, op::Op)
    # calculate actual length
    ax = abs(x)
    xlen = _ndigits(abs(x), op)
    # sign char
    sch = _signchar(x, fs.sign)
    if sch != '\0'
        xlen += 1
    end
    # prefix (e.g. 0x, 0b, 0o)
    ip = ""
    if fs.ipre
        ip = _ipre(op)
        xlen += length(ip)
    end

    # printing
    wid = fs.width
    if wid <= xlen
        _pfmt_int(out, sch, ip, 0, ax, op)
    elseif fs.zpad
        _pfmt_int(out, sch, ip, wid-xlen, ax, op)
    else
        a = fs.align
        if a == '<'
            _pfmt_int(out, sch, ip, 0, ax, op)
            _repwrite(out, fs.fill, wid-xlen)
        else
            _repwrite(out, fs.fill, wid-xlen)
            _pfmt_int(out, sch, ip, 0, ax, op)
        end
    end
end

### print floating point numbers

function _pfmt_float(out::IO, sch::Char, zs::Integer, intv::Real, decv::Real, prec::Int)
    # print sign
    if sch != '\0'
        write(out, sch)
    end
    # print padding zeros
    if zs > 0
        _repwrite(out, '0', zs)
    end
    idecv = round(Integer, decv * exp10(prec))
    if idecv == exp10(prec)
        intv += 1
        idecv = 0
    end
    # print integer part
    if intv == 0
        write(out, '0')
    else
        _pfmt_intdigits(out, intv, _Dec())
    end
    # print decimal point
    write(out, '.')
    # print decimal part
    if prec > 0
        nd = _ndigits(idecv, _Dec())
        if nd < prec
            _repwrite(out, '0', prec - nd)
        end
        _pfmt_intdigits(out, idecv, _Dec())
    end
end

function _pfmt_f(out::IO, fs::FormatSpec, x::AbstractFloat)
    # separate sign, integer, and decimal part
    ax = abs(x)
    sch = _signchar(x, fs.sign)
    intv = trunc(Integer, ax)
    decv = ax - intv

    # calculate length
    xlen = _ndigits(intv, _Dec()) + 1 + fs.prec
    if sch != '\0'
        xlen += 1
    end

    # print
    wid = fs.width
    if wid <= xlen
        _pfmt_float(out, sch, 0, intv, decv, fs.prec)
    elseif fs.zpad
        _pfmt_float(out, sch, wid-xlen, intv, decv, fs.prec)
    else
        a = fs.align
        if a == '<'
            _pfmt_float(out, sch, 0, intv, decv, fs.prec)
            _repwrite(out, fs.fill, wid-xlen)
        else
            _repwrite(out, fs.fill, wid-xlen)
            _pfmt_float(out, sch, 0, intv, decv, fs.prec)
        end
    end
end

function _pfmt_floate(out::IO, sch::Char, zs::Integer, u::Real, prec::Int, e::Int, ec::Char)
    intv = trunc(Integer,u)
    decv = u - intv
    if round(Integer, decv * exp10(prec)) == exp10(prec)
        intv += 1
        if intv == 10
            intv = 1
            e += 1
        end
        decv = 0.
    end
    _pfmt_float(out, sch, zs, intv, decv, prec)
    write(out, ec)
    if e >= 0
        write(out, '+')
    else
        write(out, '-')
        e = -e
    end
    (e1, e2) = divrem(e, 10)
    write(out, @compat Char('0' + e1))
    write(out, @compat Char('0' + e2))
end

function _pfmt_e(out::IO, fs::FormatSpec, x::AbstractFloat)
    # extract sign, significand, and exponent
    ax = abs(x)
    sch = _signchar(x, fs.sign)
    if ax == 0.0
        e = 0
        u = zero(x)
    else
        e = floor(Integer,log10(ax))  # exponent
        u = ax / exp10(e)  # significand
    end

    # calculate length
    xlen = 6 + fs.prec
    if sch != '\0'
        xlen += 1
    end

    # print
    ec = isupper(fs.typ) ? 'E' : 'e'
    wid = fs.width
    if wid <= xlen
        _pfmt_floate(out, sch, 0, u, fs.prec, e, ec)
    elseif fs.zpad
        _pfmt_floate(out, sch, wid-xlen, u, fs.prec, e, ec)
    else
        a = fs.align
        if a == '<'
            _pfmt_floate(out, sch, 0, u, fs.prec, e, ec)
            _repwrite(out, fs.fill, wid-xlen)
        else
            _repwrite(out, fs.fill, wid-xlen)
            _pfmt_floate(out, sch, 0, u, fs.prec, e, ec)
        end
    end
end

function _pfmt_g(out::IO, fs::FormatSpec, x::AbstractFloat)
    # number decomposition
    ax = abs(x)
    if 1.0e-4 <= ax < 1.0e6
        _pfmt_f(out, fs, x)
    else
        _pfmt_e(out, fs, x)
    end
end

function _pfmt_specialf(out::IO, fs::FormatSpec, x::AbstractFloat)
    if isinf(x)
        if x > 0
            _pfmt_s(out, fs, "Inf")
        else
            _pfmt_s(out, fs, "-Inf")
        end
    else
        @assert isnan(x)
        _pfmt_s(out, fs, "NaN")
    end
end


text:
include("in_source.jl")
include("out_recording.jl")
include("out_recwf.jl")
include("out_wavefield.jl")

function runsimu(model::Union{elmod2d,acmod2d},
    sou::Array{MTsource};
    rec = nothing,
    showevery = 100,
    RecDataPath = nothing,
    WFDataPath = nothing,
    Slices = nothing,
    OptCpnt = nothing)

    message(rec,RecDataPath,Slices,WFDataPath,OptCpnt)
#=== No output ===#
    if rec == nothing && Slices == nothing
        for it = 1:model.medium.nT
            for isn in 1:length(sou)
                addmt!(model.wf, sou[isn], it)
            end
            run!(model)
        end

#=== Output Recordings only ===#
  elseif rec != nothing && Slices == nothing
      simurec(model,sou,rec,OptCpnt,RecDataPath)

#=== Output Wavefields only ===#
  elseif rec == nothing && Slices != nothing
      simuwf(model,sou,Slices,OptCpnt,WFDataPath)

#=== Output both Recordings and Wavefields ===#
  elseif rec != nothing && Slices != nothing
      simurecwf(model,sou,rec,RecDataPath,Slices,WFDataPath,OptCpnt)
  end
end

ext: jl


text:
# Julia code
# Copyright (c) 2022 the SLICOTMath.jl developers
# Portions extracted from SLICOT-Reference distribution:
# Copyright (c) 2002-2020 NICONET e.V.

function run_mb04fd(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    NMAX = 50
    LDA  = NMAX÷2
    LDB = NMAX÷2
    LDDE = NMAX÷2
    LDFG = NMAX÷2
    LDQ = NMAX
    LDWORK = 3*NMAX*NMAX÷4
    A = Array{Float64,2}(undef, LDA,NMAX÷2)
    DE = Array{Float64,2}(undef, LDDE,NMAX÷2+1)
    B = Array{Float64,2}(undef, LDB,NMAX÷2)
    FG = Array{Float64,2}(undef, LDFG,NMAX÷2+1)
    Q = Array{Float64,2}(undef, LDQ,NMAX)
    ALPHAR = Array{Float64,1}(undef, NMAX÷2)
    ALPHAI = Array{Float64,1}(undef, NMAX÷2)
    BETA = Array{Float64,1}(undef, NMAX÷2)
    DWORK = Array{Float64,1}(undef, LDWORK)
    IWORK = Array{BlasInt,1}(undef, NMAX÷2+1)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    JOB =  vs[1][1]
    COMPQ =  vs[2][1]
    N = parse(BlasInt, vs[3])

    vs = String[]
    _isz,_jsz = (N÷2,N÷2)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (N÷2,N÷2+1)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       DE[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (N÷2,N÷2)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       B[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (N÷2,N÷2+1)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       FG[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
if ( N<0 || N>NMAX || mod( N, 2 )!=0 )
else
# interp call 1

    INFO = SLICOT.mb04fd!(JOB, COMPQ, N, A, DE, B, FG, Q, ALPHAR, ALPHAI, BETA)
    @test INFO == 0
    INFO == 0 || return

if ( INFO!=0 )
else
if ( LSAME( JOB, 'T' ) )
# interp output 1
    println(io,"A:")
    _nc = N÷2
    _nr = N÷2
    show(io,"text/plain",A[1:_nr,1:_nc])
    println(io,)

end # if
# interp output 2
    println(io,"DE:")
    _nc = N÷2+1
    _nr = N÷2
    show(io,"text/plain",DE[1:_nr,1:_nc])
    println(io,)

if ( LSAME( JOB, 'T' ) )
# interp output 3
    println(io,"B:")
    _nc = N÷2
    _nr = N÷2
    show(io,"text/plain",B[1:_nr,1:_nc])
    println(io,)

end # if
# interp output 4
    println(io,"FG:")
    _nc = N÷2+1
    _nr = N÷2
    show(io,"text/plain",FG[1:_nr,1:_nc])
    println(io,)

if ( !LSAME( COMPQ, 'N' ) )
# interp output 5
    println(io,"Q:")
    _nc = N
    _nr = N
    show(io,"text/plain",Q[1:_nr,1:_nc])
    println(io,)

end # if
# interp output 6
    println(io,"ALPHAR:")
    _nr = N÷2
    show(io,"text/plain",ALPHAR[1:_nr])
    println(io,)

# interp output 7
    println(io,"ALPHAI:")
    _nr = N÷2
    show(io,"text/plain",ALPHAI[1:_nr])
    println(io,)

# interp output 8
    println(io,"BETA:")
    _nr = N÷2
    show(io,"text/plain",BETA[1:_nr])
    println(io,)

end # if
end # if
    close(f)
end # run_X()


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

mutable struct ConnectionMonitorQueryResult <: SwaggerModel
    sourceStatus::Any # spec type: Union{ Nothing, String } # spec name: sourceStatus
    states::Any # spec type: Union{ Nothing, Vector{ConnectionStateSnapshot} } # spec name: states

    function ConnectionMonitorQueryResult(;sourceStatus=nothing, states=nothing)
        o = new()
        validate_property(ConnectionMonitorQueryResult, Symbol("sourceStatus"), sourceStatus)
        setfield!(o, Symbol("sourceStatus"), sourceStatus)
        validate_property(ConnectionMonitorQueryResult, Symbol("states"), states)
        setfield!(o, Symbol("states"), states)
        o
    end
end # type ConnectionMonitorQueryResult

const _property_map_ConnectionMonitorQueryResult = Dict{Symbol,Symbol}(Symbol("sourceStatus")=>Symbol("sourceStatus"), Symbol("states")=>Symbol("states"))
const _property_types_ConnectionMonitorQueryResult = Dict{Symbol,String}(Symbol("sourceStatus")=>"String", Symbol("states")=>"Vector{ConnectionStateSnapshot}")
Base.propertynames(::Type{ ConnectionMonitorQueryResult }) = collect(keys(_property_map_ConnectionMonitorQueryResult))
Swagger.property_type(::Type{ ConnectionMonitorQueryResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ConnectionMonitorQueryResult[name]))}
Swagger.field_name(::Type{ ConnectionMonitorQueryResult }, property_name::Symbol) =  _property_map_ConnectionMonitorQueryResult[property_name]

const _allowed_ConnectionMonitorQueryResult_sourceStatus = ["Unknown", "Active", "Inactive"]

function check_required(o::ConnectionMonitorQueryResult)
    true
end

function validate_property(::Type{ ConnectionMonitorQueryResult }, name::Symbol, val)
    if name === Symbol("sourceStatus")
        Swagger.validate_param(name, "ConnectionMonitorQueryResult", :enum, val, _allowed_ConnectionMonitorQueryResult_sourceStatus)
    end
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

using Unitful

const TEMP_REF = 20.0
const PRESSURE_REF = 1.0
const TEMP_REF_UNITFUL = TEMP_REF * u"°C"
export TEMP_REF, PRESSURE_REF

const DISPFORM_NAMES = ["Schott", "Sellmeier1", "Herzberger", "Sellmeier2", "Conrady", "Sellmeier3", "HandbookOfOptics1", "HandbookOfOptics2", "Sellmeier4", "Extended1", "Sellmeier5", "Extended2", "Extended3"]
const STATUS = ["Standard", "Preferred", "Obsolete", "Special", "Melt"]

#TODO need to fix this to look up the path in the Artifacts.toml file
# const AGFGLASSCAT_PATH = joinpath(JL_DIR, "AGFGlassCat.jl")

# NOTE: if you change JL_DIR or AGFGLASSCAT_PATH, you also need to change the include statement in AGFFileReader.jl

max_stars_repo_name: rambunctiousapple/GlassCat

max_issues_repo_name: rambunctiousapple/GlassCat

max_forks_repo_name: rambunctiousapple/GlassCat


text:
using Omega
using Test

function testcond()
  x = normal(0.0, 1.0)
  x_ = cond(x, x < 1.0)
  x_samples = rand(x_, 10000, alg = RejectionSample)
  x_samples = filter(s -> s != nothing, x_samples)
  @test maximum(x_samples) < 1.0
  x__ = cond(x_, x > -1.0)
  x__samples = rand(x__, 10000, alg = RejectionSample)
  x__samples = filter(s -> s != nothing, x__samples)
  @test minimum(x__samples) > -1.0
end

testcond()

function testcond2()
  x = cond(poisson(2.0), iseven)
  xsamples = rand(x, 100, alg = RejectionSample)
  @test all(iseven, xsamples)
end

testcond2()


text:
import FFTW
using LinearAlgebra

function setup_threading(;n_fft=1, n_blas=Threads.nthreads())
    n_julia = Threads.nthreads()
    FFTW.set_num_threads(n_fft)
    BLAS.set_num_threads(n_blas)
    mpi_master() && @info "Threading setup:" n_fft n_blas n_julia
end

"""
Convenience function to disable all threading in DFTK and assert that Julia threading
is off as well.
"""
function disable_threading()
    n_julia = Threads.nthreads()
    n_julia > 1 && mpi_master() && error(
        "Julia currently uses $n_julia threads. Ensure that the environment variable " *
        "JULIA_NUM_THREADS is unset and julia does not get the `-t` flag passed."
    )
    @assert n_julia == 1  # To exit in non-master MPI nodes
    setup_threading(;n_fft=1, n_blas=1)
end


text:
using ITensors
using Test

@testset "Sweeps constructor" begin

  #Sweeps
  #1 cutoff=1.0E-12, maxdim=50, mindim=10, noise=1.0E-07
  #2 cutoff=1.0E-12, maxdim=100, mindim=20, noise=1.0E-08
  #3 cutoff=1.0E-12, maxdim=200, mindim=20, noise=1.0E-10
  #4 cutoff=1.0E-12, maxdim=400, mindim=20, noise=0.0E+00
  #5 cutoff=1.0E-12, maxdim=800, mindim=20, noise=1.0E-11
  #6 cutoff=1.0E-12, maxdim=800, mindim=20, noise=0.0E+00

  sweep_args =
    [
     "maxdim" "mindim" "cutoff" "noise"
      50       10       1e-12    1e-7
      100      20       1e-12    1e-8
      200      20       1e-12    1e-10
      400      20       1e-12    0
      800      20       1e-12    1e-11
      800      20       1e-12    0
     ]

  @testset "Don't specify nsweep" begin
    nsw = size(sweep_args, 1) - 1
    sw = Sweeps(sweep_args)

    @test nsweep(sw) == nsw

    @test maxdim(sw, 1) == 50
    @test maxdim(sw, 2) == 100
    @test maxdim(sw, 3) == 200
    @test maxdim(sw, 4) == 400
    for n in 5:nsw
      @test maxdim(sw, n) == 800
    end

    @test mindim(sw, 1) == 10
    for n in 2:nsw
      @test mindim(sw, n) == 20
    end

    for n in 1:nsw
      @test cutoff(sw, n) == 1e-12
    end

    @test noise(sw, 1) == 1e-7
    @test noise(sw, 2) == 1e-8
    @test noise(sw, 3) == 1e-10
    @test noise(sw, 4) == 0
    @test noise(sw, 5) == 1e-11
    @test noise(sw, 6) == 0
  end

  @testset "Specify nsweep, more than data" begin
    nsw = 7
    sw = Sweeps(nsw, sweep_args)

    @test nsweep(sw) == nsw

    @test maxdim(sw, 1) == 50
    @test maxdim(sw, 2) == 100
    @test maxdim(sw, 3) == 200
    @test maxdim(sw, 4) == 400
    for n in 5:nsw
      @test maxdim(sw, n) == 800
    end

    @test mindim(sw, 1) == 10
    for n in 2:nsw
      @test mindim(sw, n) == 20
    end

    for n in 1:nsw
      @test cutoff(sw, n) == 1e-12
    end

    @test noise(sw, 1) == 1e-7
    @test noise(sw, 2) == 1e-8
    @test noise(sw, 3) == 1e-10
    @test noise(sw, 4) == 0
    @test noise(sw, 5) == 1e-11
    @test noise(sw, 6) == 0
    @test noise(sw, 7) == 0
  end

  @testset "Specify nsweep, less than data" begin
    nsw = 5
    sw = Sweeps(nsw, sweep_args)

    @test nsweep(sw) == nsw

    @test maxdim(sw, 1) == 50
    @test maxdim(sw, 2) == 100
    @test maxdim(sw, 3) == 200
    @test maxdim(sw, 4) == 400
    for n in 5:nsw
      @test maxdim(sw, n) == 800
    end

    @test mindim(sw, 1) == 10
    for n in 2:nsw
      @test mindim(sw, n) == 20
    end

    for n in 1:nsw
      @test cutoff(sw, n) == 1e-12
    end

    @test noise(sw, 1) == 1e-7
    @test noise(sw, 2) == 1e-8
    @test noise(sw, 3) == 1e-10
    @test noise(sw, 4) == 0
    @test noise(sw, 5) == 1e-11
  end

  @testset "Variable types of input" begin
    sw = Sweeps(5)
    setnoise!(sw,1E-8,0)
    @test noise(sw,1) ≈ 1E-8
    @test noise(sw,2) ≈ 0.0
    @test noise(sw,3) ≈ 0.0
    setcutoff!(sw,0,1E-8,0,1E-12)
    @test cutoff(sw,1) ≈ 0.0
    @test cutoff(sw,2) ≈ 1E-8
    @test cutoff(sw,3) ≈ 0.0
    @test cutoff(sw,4) ≈ 1E-12
  end
end

nothing


text:
# @testset "Savitzky-Golay definition" begin
#     sg=SavitzkyGolay_Filter(rand(11))
#     @test offset(sg) == 5
#     @test range(sg) == -5:5
#     @test length(sg) == 11
# end

@testset "Savitzky-Golay" begin

    m=SG_Filter(Float64,halfWidth=2,degree=0)

    @test filter(m,derivativeOrder=0) ≈ [1/5; 1/5; 1/5; 1/5; 1/5]

    m=SG_Filter(Float64,halfWidth=3,degree=2)

    @test filter(m,derivativeOrder=1) ≈ [-(3/28); -(1/14); -(1/28); 0;  (1/28); (1/14); (3/28)]
    @test filter(m,derivativeOrder=2) ≈ [5/42; 0; -(1/14); -(2/21); -(1/14); 0; 5/42]

    @test maxDerivativeOrder(m) == 2
    @test length(m) == 2*3+1
    @test polynomialOrder(m) == 2
end

max_stars_repo_name: vincent-picaud/DirectConvolution

max_issues_repo_name: vincent-picaud/DirectConvolution

max_forks_repo_name: vincent-picaud/DirectConvolution


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -17.39516631142554 <= q <= 12.050938789121952)

                     add_NL_constraint(m, :(swish(0.1892165831624366 + -0.6797303891528506*swish(0.23382152176968107 + 0.970568112692761*$(x[1]) + 0.9147961292868287*$(x[2]) + 0.25099364821370385*$(x[3])) + -0.5563706493225857*swish(0.1824369939040409 + -0.636682724457545*$(x[1]) + 0.684450108697257*$(x[2]) + -0.3896097552087272*$(x[3])) + -0.5881166064998418*swish(-0.7017137443246066 + -0.579761469455264*$(x[1]) + -0.38090529589576505*$(x[2]) + -0.30474275998045197*$(x[3])) + -0.5324749908347566*swish(0.495407339331571 + -0.11045663155169017*$(x[1]) + -0.19556503402268088*$(x[2]) + 0.8751710578699647*$(x[3]))) + swish(-0.7087227741366906 + -0.3339546766337098*swish(0.23382152176968107 + 0.970568112692761*$(x[1]) + 0.9147961292868287*$(x[2]) + 0.25099364821370385*$(x[3])) + -0.8987724666577948*swish(0.1824369939040409 + -0.636682724457545*$(x[1]) + 0.684450108697257*$(x[2]) + -0.3896097552087272*$(x[3])) + 0.3907505378456695*swish(-0.7017137443246066 + -0.579761469455264*$(x[1]) + -0.38090529589576505*$(x[2]) + -0.30474275998045197*$(x[3])) + -0.7983638213649877*swish(0.495407339331571 + -0.11045663155169017*$(x[1]) + -0.19556503402268088*$(x[2]) + 0.8751710578699647*$(x[3]))) + swish(0.16171979073463483 + 0.31171990946083294*swish(0.23382152176968107 + 0.970568112692761*$(x[1]) + 0.9147961292868287*$(x[2]) + 0.25099364821370385*$(x[3])) + 0.2757921127232579*swish(0.1824369939040409 + -0.636682724457545*$(x[1]) + 0.684450108697257*$(x[2]) + -0.3896097552087272*$(x[3])) + -0.732639957443523*swish(-0.7017137443246066 + -0.579761469455264*$(x[1]) + -0.38090529589576505*$(x[2]) + -0.30474275998045197*$(x[3])) + -0.947011017876227*swish(0.495407339331571 + -0.11045663155169017*$(x[1]) + -0.19556503402268088*$(x[2]) + 0.8751710578699647*$(x[3]))) + swish(-0.17479526548317725 + -0.7121643152938684*swish(0.23382152176968107 + 0.970568112692761*$(x[1]) + 0.9147961292868287*$(x[2]) + 0.25099364821370385*$(x[3])) + -0.20573796119512977*swish(0.1824369939040409 + -0.636682724457545*$(x[1]) + 0.684450108697257*$(x[2]) + -0.3896097552087272*$(x[3])) + 0.920351043144791*swish(-0.7017137443246066 + -0.579761469455264*$(x[1]) + -0.38090529589576505*$(x[2]) + -0.30474275998045197*$(x[3])) + -0.8770524767395065*swish(0.495407339331571 + -0.11045663155169017*$(x[1]) + -0.19556503402268088*$(x[2]) + 0.8751710578699647*$(x[3]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Py2JlRuntime
# runtests/benchmark.jl, line 1
function var".sum_by_2"(_f′, _xs′)
    f = _f′
    xs = _xs′
    local f
    local xs
    local s
    local e
    # runtests/benchmark.jl, line 2
    s = jpy_literal(0)
    # runtests/benchmark.jl, line 3
    var".iter_1" = jpy_getiter(xs)
    @noscope while jpy_movenext(var".iter_1")
        e = jpy_getcurrent(var".iter_1")
        # runtests/benchmark.jl, line 4
        s = jpy_add(s, jpy_call(f, (e, ), (;)))
        jpy_none;
    end
    # runtests/benchmark.jl, line 5
    return s
    jpy_none;
end
const sum_by = var".sum_by_2"
jpy_none;


text:
"""

# cmdline

Recursively parse the model to construct command line.

### Method
```julia
cmdline(m)
```

### Required arguments
```julia
* `m::CmdStanSampleModel`                : CmdStanSampleModel
```

### Related help
```julia
?CmdStanSampleModel                      : Create a CmdStanSampleModel
```
"""
function cmdline(m::Union{VariationalModel, Variational}, id)

  #=
  `/Users/rob/.julia/dev/StanVariational/examples/Bernoulli/tmp/bernoulli
  variational algorithm=meanfield grad_samples=1 elbo_samples=100
  iter=10000 tol_rel_obj=0.01 eval_elbo=100 output_samples=10000
  random seed=-1 init=2 id=1
  data file=/Users/rob/.julia/dev/StanVariational/examples/Bernoulli/tmp/bernoulli_data_1.R
  output file=/Users/rob/.julia/dev/StanVariational/examples/Bernoulli/tmp/bernoulli_chain_1.csv
  refresh=100`
  =#

  cmd = ``
  if isa(m, VariationalModel)
    # Handle the model name field for unix and windows
    cmd = `$(m.exec_path)`

    # Sample() specific portion of the model
    cmd = `$cmd $(cmdline(getfield(m, :method), id))`

    # Common to all models
    cmd = `$cmd random seed=$(getfield(m, :seed).seed)`

    # Init file required?
    if length(m.init_file) > 0 && isfile(m.init_file[id])
      cmd = `$cmd init=$(m.init_file[id])`
    else
      cmd = `$cmd init=$(m.init.bound)`
    end

    # Data file required?
    if length(m.data_file) > 0 && isfile(m.data_file[id])
      cmd = `$cmd id=$(id) data file=$(m.data_file[id])`
    end

    # Output options
    cmd = `$cmd output`
    if length(getfield(m, :output).file) > 0
      cmd = `$cmd file=$(string(getfield(m, :output).file))`
    end
    if length(m.diagnostic_file) > 0
      cmd = `$cmd diagnostic_file=$(string(getfield(m, :output).diagnostic_file))`
    end
    cmd = `$cmd refresh=$(string(getfield(m, :output).refresh))`

  else

    # The 'recursive' part
    cmd = `$cmd $(split(lowercase(string(typeof(m))), '.')[end])`
    for name in fieldnames(typeof(m))
      cmd = `$cmd $(name)=$(getfield(m, name))`
    end
  end

  cmd

end


text:
function preconditioner(d, q)
    P_preallocated = Diagonal(zeros(d*(q+1)))

    @fastmath @inbounds function P(h)
        val = h^(-q-1/2)
        @simd for j in 0:q
            @simd for i in 1:d
                # P_preallocated[j*d + i,j*d + i] = h^(j-q-1/2)
                P_preallocated[j*d + i,j*d + i] = val
            end
            val *= h
        end
        return P_preallocated
    end

    return P
end


text:
import Primes: isprime
import Base: +,-,*,/,^,==,inv,show

# Field has only one attribute, p
# Implicitly Fp^2 initialized as Fp/<x^2+1>
# Requires p = 3 (mod 4) for polynomial to be irredusible
# The parameter T determines the type of integer coefficients (i.e. Int64, Int32, BigInt etc...)
struct Fp2{T <: Integer}
    p::T
end

# FieldElements point to their parent field
# Consists of integers a,b representing the integer a + b*i
struct Fp2Elem{T <: Integer}
    a::T
    b::T
    field::Fp2{T}
end

# Field-Constructor.
function Fp2(p::T) where T <: Integer
    !isprime(p) && error("p must prime")
    p%4 != 3 && error("Prime p must be 3 mod 4")
    return Fp2{T}(p)
end

# Element-Constructor
# ASSUMES a,b are both between 0 and p for speed.
# Arithmetic functions will make sure of this
function (F::Fp2{T})(a::T,b::T) where T <: Integer
    return Fp2Elem{T}(a, b, F)
end

############################
#
# Printing
#
############################

function show(io::IO, F::Fp2{T}) where {T <: Integer}
    print(io, "F_$(F.p)[X]/<X^2+1>")
end

function show(io::IO, F::Fp2Elem{T}) where {T <: Integer}
    print(io, "$(F.a) + $(F.b)*i")
end

############################
#
# Some helper functions
#
############################

function zero(F::Fp2{T}) where T <: Integer
    return F(T(0),T(0))
end

function one(F::Fp2{T}) where T <: Integer
    return F(T(1),T(0))
end

function iszero(x::Fp2Elem{T}) where T <: Integer
    return x.a == 0 && x.b == 0
end

function isone(x::Fp2Elem{T}) where T <: Integer
    return x.a == 1 && x.b == 0
end

function ==(x::Fp2Elem{T}, y::Fp2Elem{T}) where T <: Integer
    return (x.a == y.a) && (x.b == y.b)
end

############################
#
# Basic operators
#
############################
# All operators trusts that elements are from the same ring
# Additive inverse
function -(x::Fp2Elem{T}) where T <: Integer
    a_new = mod(-x.a, x.field.p)
    b_new = mod(-x.b, x.field.p)
    return x.field(a_new, b_new)
end

# Addition
function +(x::Fp2Elem{T}, y::Fp2Elem{T}) where T <: Integer
    a_new = mod((x.a + y.a), x.field.p)
    b_new = mod((x.b + y.b), x.field.p)
    return x.field(a_new, b_new)
end

# Subtraction
function -(x::Fp2Elem{T}, y::Fp2Elem{T}) where T <: Integer
    return x+(-y)
end

# Multiplication
function *(x::Fp2Elem{T}, y::Fp2Elem{T}) where T <: Integer
    a_new = mod((x.a*y.a - x.b*y.b), x.field.p)
    b_new = mod((x.a*y.b + y.a*x.b), x.field.p)
    return x.field(a_new, b_new)
end

# Scalar multiplication on left
function *(x::Integer, y::Fp2Elem{T}) where T <: Integer
    return y.field(mod(T(x), y.field.p),T(0))*y
end

# Scalar multiplication on right
function *(x::Fp2Elem{T}, y::Integer) where T <: Integer
    return x*x.field(mod(T(y), x.field.p),T(0))
end

# Multiplicative inverse
function inv(x::Fp2Elem{T}) where T <: Integer
    x.a == 0 && x.b == 0 && error("0 is not invertible...")
    #x^-1 = x.a*(x.a^2+x.b^2)^-1 + (-x.b)(x.a^2+x.b^2)^-1 * i
    d = invmod((x.a*x.a + x.b*x.b),x.field.p)
    a_new = mod((x.a*d), x.field.p)
    b_new = mod((-x.b*d), x.field.p)
    return x.field(a_new, b_new)
end

# Exponensiation
# Uses the simple square and multiply algorithm
function ^(x::Fp2Elem{T}, y::Integer) where T <: Integer
    if iszero(x) || y == 1
        return x
    end
    if y == 0
        return one(x.field)
    end
    if y < 0
        x = inv(x)
        y = -y
    end
    bit = T(1) << (ndigits(y, base = 2) - 1)
    z = x
    bit >>= 1
    while bit != 0
        z = z*z
        if (bit & y) != 0
            z *= x
        end
        bit >>=1
    end
    return z
end


text:
import Base: convert, promote_rule, show, real, imag, conj, abs, abs2, inv, +, -, /, *

immutable Quaternion{T<:Real} <: Number
    q0::T
    q1::T
    q2::T
    q3::T
end

Quaternion(q0::Real,q1::Real,q2::Real,q3::Real) = Quaternion(promote(q0,q1,q2,q3)...)

convert{T}(::Type{Quaternion{T}}, x::Real) =
    Quaternion(convert(T,x), zero(T), zero(T), zero(T))
convert{T}(::Type{Quaternion{T}}, z::Complex) =
    Quaternion(convert(T,real(z)), convert(T,imag(z)), zero(T), zero(T))
convert{T}(::Type{Quaternion{T}}, z::Quaternion) =
    Quaternion(convert(T,z.q0), convert(T,z.q1), convert(T,z.q2), convert(T,z.q3))

promote_rule{T,S}(::Type{Complex{T}}, ::Type{Quaternion{S}}) = Quaternion{promote_type(T,S)}
promote_rule{S}(::Type{Bool}, ::Type{Quaternion{S}}) = Quaternion{S}
promote_rule{T<:Real,S}(::Type{T}, ::Type{Quaternion{S}}) = Quaternion{promote_type(T,S)}
promote_rule{T,S}(::Type{Quaternion{T}}, ::Type{Quaternion{S}}) = Quaternion{promote_type(T,S)}

function show(io::IO, z::Quaternion)
    pm(x) = x < 0 ? " - $(-x)" : " + $x"
    print(io, z.q0, pm(z.q1), "i", pm(z.q2), "j", pm(z.q3), "k")
end

real(z::Quaternion) = z.q0
imag(z::Quaternion) = z.q1

conj(z::Quaternion) = Quaternion(z.q0, -z.q1, -z.q2, -z.q3)
abs(z::Quaternion) = sqrt(z.q0*z.q0 + z.q1*z.q1 + z.q2*z.q2 + z.q3*z.q3)
abs2(z::Quaternion) = z.q0*z.q0 + z.q1*z.q1 + z.q2*z.q2 + z.q3*z.q3
inv(z::Quaternion) = conj(z)/abs2(z)

(-)(z::Quaternion) = Quaternion(-z.q0, -z.q1, -z.q2, -z.q3)

(/)(z::Quaternion, x::Real) = Quaternion(z.q0/x, z.q1/x, z.q2/x, z.q3/x)

(+)(z::Quaternion, w::Quaternion) = Quaternion(z.q0 + w.q0, z.q1 + w.q1,
                                               z.q2 + w.q2, z.q3 + w.q3)
(-)(z::Quaternion, w::Quaternion) = Quaternion(z.q0 - w.q0, z.q1 - w.q1,
                                               z.q2 - w.q2, z.q3 - w.q3)
(*)(z::Quaternion, w::Quaternion) = Quaternion(z.q0*w.q0 - z.q1*w.q1 - z.q2*w.q2 - z.q3*w.q3,
                                               z.q0*w.q1 + z.q1*w.q0 + z.q2*w.q3 - z.q3*w.q2,
                                               z.q0*w.q2 - z.q1*w.q3 + z.q2*w.q0 + z.q3*w.q1,
                                               z.q0*w.q3 + z.q1*w.q2 - z.q2*w.q1 + z.q3*w.q0)
(/)(z::Quaternion, w::Quaternion) = z*inv(w)

q = Quaternion(1,0,0,0)
x = Quaternion(0,1,1,1)

println("q = $q")
println("q*2.0+2 = $(q*2.0+2)")
println("abs((-q+x*2)/4) = ", abs((-q+x*2)/4))

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
# This file is a part of SimilaritySearch.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

using SimilaritySearch
using Dates

import SimilaritySearch:
    search, fit, push!, optimize!

export Knr

mutable struct Knr{T} <: Index
    db::Vector{T}
    refs::Vector{T}
    k::Int
    ksearch::Int
    minmatches::Int
    invindex::Vector{Vector{Int32}}
    verbose::Bool
end

"""
    fit(::Type{Knr}, dist::Function, db::AbstractVector{T}, refs::AbstractVector{T}, k::Int, minmatches::Int=1, verbose=false) where T
    fit(::Type{Knr}, dist::Function, db::AbstractVector{T}; numrefs::Int=1024, k::Int=3, minmatches::Int=1, tournamentsize::Int=3, verbose=false) where T

Creates a `Knr` index using the given references or the number of references to be used.

- `dist`: the distance function to be used
- `db`: the dataset to bbe indexed
- `refs`: the references to be used
- `k`: the number of references for representing objects inside the index
- `ksearch`: the same than `k` but used while searching
- `minmatches`: at query time, it determines the minimum number of references matched to allow a distance evaluation
- `verbose`: controls if the index must have a verbose output for its operations
- `numrefs`: if `refs` is not given, then `numrefs` objects are selected from `db` as `refs`
- `tournamentsize`: `numrefs` is specified, an incremental construction of `refs` is performed, each reference is selected as the more distant (to already selected references) among `tournamentsize` candidates
"""
function fit(::Type{Knr}, dist::Function, db::AbstractVector{T}, refs::AbstractVector{T}, k::Int, minmatches::Int=1; verbose=false) where T
    verbose && println(stderr, "Knr> refs=$(typeof(db)), k=$(k), numrefs=$(length(refs)), dist=$(dist)")
    invindex = [Vector{Int32}(undef, 0) for i in 1:length(refs)]
    seqindex = fit(Sequential, refs)

    pc = round(Int, length(db) / 20)
    for i=1:length(db)
        if verbose && i % pc == 0
            println(stderr, "Knr> advance $(round(i/length(db), digits=4)), now: $(now())")
        end

        res = search(seqindex, dist, db[i], KnnResult(k))
        for p in res
            push!(invindex[p.objID], i)
        end
    end

    Knr(db, refs, k, k, minmatches, invindex, verbose)
end

function fit(::Type{Knr}, dist::Function, db::AbstractVector{T}; numrefs::Int=1024, k::Int=3, minmatches::Int=1, tournamentsize::Int=3, verbose=false) where T
    # refs = rand(db, numrefs)
    refs = [db[x] for x in select_tournament(dist, db, numrefs, tournamentsize)]
    fit(Knr, dist, db, refs, k, minmatches, verbose=verbose)
end

"""
    search(index::Knr, dist::Function, q, res::KnnResult)

Solves the query specified by `q` and `res` using the `Knr` index
"""
function search(index::Knr, dist::Function, q, res::KnnResult)
    dz = zeros(Int16, length(index.db))
    # M = BitArray(length(index.db))
    seqindex = fit(Sequential, index.refs)
    kres = search(seqindex, dist, q, KnnResult(index.ksearch))

    for p in kres
        @inbounds for objID in index.invindex[p.objID]
            c = dz[objID] + 1
            dz[objID] = c

            if c == index.minmatches
                d = dist(q, index.db[objID])
                push!(res, objID, d)
            end
        end
    end

    return res
end

"""
    push!(index::Knr, dist::Function, obj)

Inserts `obj` into the index
"""
function push!(index::Knr, dist::Function, obj)
    push!(index.db, obj)
    seqindex = fit(Sequential, index.refs)
    res = search(seqindex, dist, obj, KnnResult(index.k))
    for p in res
        push!(index.invindex[p.objID], length(index.db))
    end

    length(index.db)
end

"""
    optimize!(index::Knr, dist::Function; recall=0.9, k=10, num_queries=128, perf=nothing)

Optimizes the index to achieve the specified recall.
"""
function optimize!(index::Knr, dist::Function; recall=0.9, k=10, num_queries=128, perf=nothing)
    index.verbose && println(stderr, "Knr> optimizing index for recall=$(recall)")
    if perf == nothing
        perf = Performance(index.db, dist; expected_k=k, num_queries=num_queries)
    end
    index.minmatches = 1
    index.ksearch = 1
    p = probe(perf, index, dist)

    while p.recall < recall && index.ksearch < length(index.refs)
        index.ksearch += 1
        index.verbose && println(stderr, "Knr> opt step ksearch=$(index.ksearch), performance $(p)")
        p = probe(perf, index, dist)

    end

    index.verbose && println(stderr, "Knr> reached performance $(p)")
    index
end


text:
using Pkg
Pkg.add("JLD")

using JLD

x = collect(-3:0.1:3)
y = collect(-3:0.1:3)

xx = reshape([xi for xi in x for yj in y], length(y), length(x))
yy = reshape([yj for xi in x for yj in y], length(y), length(x))

z = sin.(xx .+ yy .^ 2)

data_dict = Dict("x" => x, "y" => y, "z" => z)

save("data_dict.jld", data_dict)

#%% restart the REPL
using JLD
data_dict2 = load("data_dict.jld")

x2 = data_dict2["x"]
y2 = data_dict2["y"]
z2 = data_dict2["z"]

using Plots
plotly()

plot(x2, y2, z2, st = :surface, color = :ice)
savefig("img1")

#%%
using JLD
struct Person
    height::Float64
    weight::Float64
end

bob = Person(1.84, 74)

dict_new = Dict("bob" => bob)
save("bob.jld", dict_new)

#%% please restart the REPL
using JLD
struct Person
    height::Float64
    weight::Float64
end
bob2 = load("bob.jld")

bob2["bob"]

#%% please restart the REPL
using JLD

bob3 = load("bob.jld")

bob3["bob"]

bob3["bob"].height

max_stars_repo_name: JeffreySarnoff/techytok-examples

max_issues_repo_name: JeffreySarnoff/techytok-examples

max_forks_repo_name: JeffreySarnoff/techytok-examples


text:
function _log_preamble(
    fm::DMIPFormulation,
    primal_bound::Float64,
    config::AlgorithmConfig,
)
    config.silent && return nothing
    m = num_constraints(fm)
    n = num_variables(fm)
    d = length(fm.disjunction_formulaters)
    @info "Cerberus: An (experimental) solver for disjunctive mixed-integer programming."
    @info "Optimizing a model with $m rows, $n columns, and $d disjunctive constraints."
    kinds = [get_variable_kind(fm, CVI(i)) for i in 1:n]
    n_c = count(v -> v === nothing, kinds)
    n_i = count(v -> v !== nothing, kinds)
    n_b = count(v -> v isa ZO, kinds)
    @assert n_i + n_c == n
    @info "Variables: $n_c continuous, $n_i integer ($n_b binary)."
    # TODO: Log problem size after formulation of disjunctions
    # TODO: Print table of coefficient statistics
    if primal_bound < Inf
        @info "Starting primal bound of $primal_bound"
    end
    @info "    Nodes    |    Current Node    |     Objective Bounds      |    Work"
    @info " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time"
    return nothing
end

function _log_postamble(result::Result, config::AlgorithmConfig)
    config.silent && return nothing
    # TODO: Summarize cutting planes added (none right now...)
    @info Printf.@sprintf(
        "Explored %u nodes (%u simplex iterations) in %.2f seconds.",
        result.total_node_count,
        result.total_simplex_iters,
        result.total_elapsed_time_sec,
    )
    if result.termination_status == MOI.INFEASIBLE
        @info "Model is infeasible."
    elseif result.termination_status == MOI.DUAL_INFEASIBLE
        @info "Model is unbounded."
    else
        if result.termination_status == MOI.OPTIMAL
            @info "Optimal solution found."
        elseif result.primal_bound < Inf && all(!isnan, result.best_solution)
            @info "Feasible solution found."
        end
        @info Printf.@sprintf(
            "Best objective %10e, best bound %10e, gap %.4f%%.",
            result.primal_bound,
            result.dual_bound,
            _optimality_gap(result.primal_bound, result.dual_bound),
        )
    end
    return nothing
end

function _log_node_update(state::CurrentState, node_result::NodeResult)
    cost = node_result.cost
    gap = _optimality_gap(state.primal_bound, state.dual_bound)
    @info Printf.@sprintf(
        "%5u %5u   %8s %4u %4s %8s %8.5f %8s  %5.1f %5us",
        state.total_node_count,
        length(state.tree),
        if cost == Inf
            "Inf      "
        else
            Printf.@sprintf("%8.5f", cost)
        end,
        node_result.depth,
        # If infeasible, don't report int infeasibility as 0.
        if cost == Inf
            "    -"
        else
            Printf.@sprintf("%5u", node_result.int_infeas)
        end,
        if state.primal_bound == Inf
            "         -"
        else
            Printf.@sprintf("%8.5f", state.primal_bound)
        end,
        state.dual_bound,
        # If we don't have a primal bound, just don't report a gap.
        if isnan(gap)
            "      - "
        else
            (isinf(gap) ? "      ∞ " : Printf.@sprintf("%7.2f%%", gap))
        end,
        state.polling_state.period_simplex_iters /
        state.polling_state.period_node_count,
        floor(UInt, state.total_elapsed_time_sec),
    )
end

const EARLY_POLLING_CUTOFF = 10.0
const EARLY_POLLING_CADENCE = 0.25
const NORMAL_POLLING_CADENCE = 5.0

# TODO: Unit test
function _log_if_necessary(
    state::CurrentState,
    result::NodeResult,
    config::AlgorithmConfig,
)
    config.silent && return nothing
    elapsed_time_sec = state.total_elapsed_time_sec
    if elapsed_time_sec > state.polling_state.next_polling_target_time_sec
        update_dual_bound!(state)
        _log_node_update(state, result)
        if elapsed_time_sec <= EARLY_POLLING_CUTOFF
            state.polling_state.next_polling_target_time_sec +=
                EARLY_POLLING_CADENCE
        else
            state.polling_state.next_polling_target_time_sec +=
                NORMAL_POLLING_CADENCE
        end
        state.polling_state.period_node_count = 0
        state.polling_state.period_simplex_iters = 0
    end
    return nothing
end


text:
#EasyPlotPlots AbstractBuilder interface implementation
#-------------------------------------------------------------------------------

#==Types
===============================================================================#
mutable struct Builder <: EasyPlot.AbstractBuilder
	backend::Symbol
	postproc::Optional{Function}
	args::Tuple
	kwargs::Base.Iterators.Pairs
end
EasyPlot.AbstractBuilder(::DS{:PlotsJl}) = Builder #Register builder

#==Constructor interface
===============================================================================#
#Use same builder irrespective of target application:
function EasyPlot.getbuilder(::Target{T}, ::Type{Builder}, args...;
		postproc=nothing, backend=nothing, kwargs...) where T
	if isnothing(backend); backend = defaults.backend; end
	return Builder(backend, postproc, args, kwargs)
end

function EasyPlot._show(io::IO, mime::MIME, opt::EasyPlot.ShowOptions, fig::FigureMulti)
	if !EasyPlot.isauto(opt.dim)
		@warn("Cannot currently specify plot size. Using defaults.")
	end
	show(io, mime, fig.p)
end

#Currently no support for non-FigureMulti plots... but could generate a single image...:
EasyPlot._show(io::IO, mime::MIME, opt::EasyPlot.ShowOptions, fig::Figure) =
	throw(MethodError(EasyPlot._show, (io, mime, opt, fig)))

#Module does not yet support other inline formats (must figure out how)
EasyPlot._showable(mime::MIME, b::Builder) = false
#Assume PNG is supported by all backends:
EasyPlot._showable(mime::MIME"image/png", b::Builder) = true

#==Display interface
===============================================================================#
EasyPlot.displaygui(fig::FigureMulti) = display(fig.p)
function EasyPlot.displaygui(fig::FigureSng)
	for s in fig.subplots
		display(s) #Defined in Plots.jl
	end
	nothing
end

#==Build interface
===============================================================================#
function _build(sp::Plots.Subplot, eplot::EasyPlot.Plot, theme::EasyPlot.Theme)
	Plots.title!(sp, eplot.title)
	Plots.plot!(sp, legend=eplot.legend)
	fold = isa(eplot.xaxis, EasyPlot.FoldedAxis) ? eplot.xaxis : nothing

	builder = WfrmBuilder(sp, theme, fold)
	for (i, wfrm) in enumerate(eplot.wfrmlist)
		EasyPlot.addwfrm(builder, wfrm, i)
	end

	#x-axis properties:
	xscale = Symbol(eplot.xaxis)
	xmin = eplot.xext.min; xmax = eplot.xext.max

	#y-axis properties:
	ylabel = ""
	yscale = :lin
	ymin, ymax = (NaN, NaN)
	if length(eplot.ystriplist) > 0
		strip = eplot.ystriplist[1]
		ylabel = strip.axislabel
		yscale = strip.scale
		ymin = strip.ext.min; ymax = strip.ext.max
	end

	#Apply x/y labels:
	Plots.xlabel!(sp, eplot.xlabel)
	Plots.ylabel!(sp, ylabel)

	#Apply x/y scales:
	Plots.plot!(sp, xscale=scalemap[xscale], yscale=scalemap[yscale])

	#Update axis limits:
	xlims!(sp, (xmin, xmax))
	ylims!(sp, (ymin, ymax))
	return sp
end

function _build(fig::Figure, ecoll::EasyPlot.PlotCollection)
	ecoll = EasyPlot.condxfrm_multistrip(ecoll, "EasyPlotPlots") #Emulate multi-strip plots
	ncols = ecoll.ncolumns
	nsubplots = length(ecoll.plotlist)
	plt = addsubplots(fig, ncols, nsubplots)

	for (i, plot) in enumerate(ecoll.plotlist)
		sp = getsubplot(fig, i)
		_build(sp, plot, ecoll.theme)
	end

	return fig
end

function EasyPlot.build(b::Builder, ecoll::EasyPlot.PlotCollection)
	local fig = nothing
	try
		bknd = Plots.backend(b.backend) #Activate backend
		fig = Figure(b.backend)
		_build(fig, ecoll)
	finally
		#TODO: Restore state
	end
	return fig
end

#Last line


text:
function solve(n::Integer, from::Integer, to::Integer, via::Integer)
  if n == 1
    println("Move disk from $from to $to")
  else
    solve(n - 1, from, via, to)
    solve(1, from, to, via)
    solve(n - 1, via, to, from)
  end
end

solve(4, 1, 2, 3)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export Diagm

"C[:,:,k]=diagm(A[:,k])"
function Diagm(A::Array{Float32,2})
  C=zeros(Float32,size(A,1),size(A,1),size(A,2))
  for k=1:1:size(A,2)
    @inbounds C[:,:,k]=diagm(A[:,k])
  end
  C
end


text:
function compute_density_profile(orbital::T, xgrid::Array{DType,1}, obdm::Array{DType,2}; kwargs...)::Array{DType,1} where T<:Orbital
    """ Takes the one-body density matrix and returns the spatial density-profile
        on the specified position grid (unnormalized).
    """
    density = zeros(DType, (length(xgrid),))
    for i=1:size(obdm)[1]
        di = conj.(orbital.(i, xgrid))
        for j=1:size(obdm)[2]
            density .+= di.* orbital.(j, xgrid) .* obdm[i,j]
        end
    end
    return -density
end


text:
using PyPlot, Seismic

dt = 0.002
fn = 1/(2*dt)

# Ricker
w1 = Ricker(dt=dt)
nw = length(w1)
nc = floor(Int, nw/2)
t1 = dt*collect(-nc:1:nc)
nf = 8*nextpow2(nw)
df = 1/(nf*dt)
f1 = df*collect(0:1:nf-1)
wpad = cat(1,w1,zeros(nf-nw))
W1 = abs(fft(wpad))
W1 = W1/maximum(W1)

# Ormsby
w2 = Ormsby(dt=dt, f=[5.0, 10.0, 30.0, 55.0])
nw = length(w2)
t2 = dt*collect(0:1:nw-1)
nf = 8*nextpow2(nw)
df = 1/(nf*dt)
f2 = df*collect(0:1:nf-1)
wpad = cat(1,w2,zeros(nf-nw))
W2 = abs(fft(wpad))
W2 = W2/maximum(W2)

# Berlage
w3 = Berlage(dt=dt)
nw = length(w3)
t3 = dt*collect(0:1:nw-1)
nf = 8*nextpow2(nw)
df = 1/(nf*dt)
f3 = df*collect(0:1:nf-1)
wpad = cat(1,w3,zeros(nf-nw))
W3 = abs(fft(wpad))
W3 = W3/maximum(W3)

figure("Suite of wavelets in time")

subplot(221)
plot(t1, w1)
axis("tight")
xlabel("Time (s)")
title("Ricker wavelet")

subplot(222)
plot(t2, w2)
axis("tight")
xlabel("Time (s)")
title("Ormsby wavelet")

subplot(223)
plot(t3, w3)
axis("tight")
xlabel("Time (s)")
title("Berlage wavelet")

tight_layout()

figure("Suite of amplitude spectrums")

subplot(221)
plot(f1, W1)
xlim([0, fn])
xlabel("Frequency (Hz)")
title("Ricker spectrum")

subplot(222)
plot(f2, W2)
xlim([0, fn])
xlabel("Frequency (Hz)")
title("Ormsby spectrum")

subplot(223)
plot(f3, W3)
xlim([0, fn])
xlabel("Frequency (Hz)")
title("Berlage spectrum")

tight_layout()


text:
using PyCall
using Tweaks

py"""
class X:
    def __str__(self): return 'Y'
"""

@testset "print" begin
    buf = IOBuffer()
    println(buf, py"X()")
    seek(buf, 0)
    @test read(buf, String) == "Y\n"
end

@testset "string" begin
    @test string(py"X()") == "Y"
end


text:
using IterableTables
using TimeSeries
using DataFrames
using TableTraits
using IteratorInterfaceExtensions
using Test

@testset "TimeSeries" begin

dates  = collect(Date(1999,1,1):Day(1):Date(1999,1,3))

source_ta1 = TimeArray(dates, collect(1:length(dates)))

@test IteratorInterfaceExtensions.isiterable(source_ta1) == true

source_it = IteratorInterfaceExtensions.getiterator(source_ta1)

@test length(source_it) == 3

df1 = DataFrame(source_ta1)
@test size(df1) == (3,2)
@test df1[:timestamp] == dates
@test df1[:A] == [1,2,3]

source_ta2 = TimeArray(dates, collect(1:length(dates)), [:a])
df2 = DataFrame(source_ta2)
@test size(df2) == (3,2)
@test df2[:timestamp] == dates
@test df2[:a] == [1,2,3]

source_ta3 = TimeArray(dates, hcat(collect(1:length(dates)), collect(length(dates):-1:1)), [:a, :b])
df3 = DataFrame(source_ta3)
@test size(df3) == (3,3)
@test df3[:timestamp] == dates
@test df3[:a] == [1,2,3]
@test df3[:b] == [3,2,1]

source_ta4 = TimeArray(dates, hcat(collect(1:length(dates)), collect(length(dates):-1:1)))
df4 = DataFrame(source_ta4)
@test size(df4) == (3,3)
@test df4[:timestamp] == dates
@test df4[:A] == [1,2,3]
@test df4[:B] == [3,2,1]

source_df = DataFrame(a=[4.,5.], timestamp=[Date(1999,1,1),Date(1999,1,2)], b=[6.,8.])
ta2 = TimeArray(source_df, timestamp=:timestamp)
@test length(ta2) == 2
@test TimeSeries.colnames(ta2) == [:a, :b]
@test TimeSeries.values(ta2) == [4. 6.;5. 8.]
@test TimeSeries.timestamp(ta2) == [Date(1999,1,1),Date(1999,1,2)]

end


text:
# main.jl - Code to generate the flashcards in newline-delimited format

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

DEFAULT_DIR = "..\\Julia-Associativity-Flashcard-Generation\\"

julia_ops_syntax_0 = [:(.), :(::)]
julia_ops_exponentiation = [:(^)]
julia_ops_bitshifts = [:(<<), :(>>), :(>>>)]
julia_ops_fractions = [:(//)]
julia_ops_multiplication = [:(*), :(/), :(%), :(&), :(\), :(÷)]
julia_ops_addition = [:(+), :(-), :(|), :(⊻)]
julia_ops_syntax_1 = [:(:), :(..)]
julia_ops_syntax_2 = [:(|>)]
julia_ops_syntax_3 = [:(<|)]
julia_ops_comparisons = [:(>), :(<), :(>=), :(<=), :(==), :(===), :(!=), :(!==), :(<:)]
julia_ops_control_flow = [:(&&), :(||)]
julia_ops_pairs = [:(=>)]
julia_ops_assignments = [:(=), :(+=), :(-=), :(*=), :(/=), :(//=), :(\=), :(^=), :(÷=), :(%=), :(|=), :(&=), :(⊻=), :(<<=), :(>>=), :(>>>=)]

# https://docs.julialang.org/en/v1/manual/mathematical-operations/#Operator-Precedence-and-Associativity-1
# Source of the operators list.

julia_common_operators = [julia_ops_syntax_0
                          julia_ops_exponentiation
                          julia_ops_bitshifts
                          julia_ops_fractions
                          julia_ops_multiplication
                          julia_ops_addition
                          julia_ops_syntax_1
                          julia_ops_syntax_2
                          julia_ops_syntax_3
                          julia_ops_comparisons
                          julia_ops_control_flow
                          julia_ops_pairs
                          julia_ops_assignments]

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

# Code to generate random flashcards for each operator.
function generate_flashcards(stream_out=stdout, ops=julia_common_operators)
    for operator1 in ops
        operator2 = ops[rand(1:end)]

        while operator2 == operator1
            operator2 = ops[rand(1:end)]
        end

        p1 = Base.operator_precedence(operator1)
        p2 = Base.operator_precedence(operator2)

        write(stream_out, "Which has higher precedence in Julia: $operator1, or $operator2?\n")

        if p1 == p2
            write(stream_out, "They have the same precedence. Be careful!\n\n")
        elseif p1 > p2
            write(stream_out, "$operator1 has higher precedence.\n\n")
        elseif p1 < p2
            write(stream_out, "$operator2 has higher precedence.\n\n")
        end
    end
end

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

function make_card_file(dir, filename, ops)
    try
        cd(dir)
        file = open(filename, "w")

        println("Generating and writing flashcards to $file ...")
        generate_flashcards(file, ops)
        println("Done. Closing $file .")
        close(file)
    catch
        println("Uh oh! Something went wrong...\n")
        println("Not to worry. We'll just generate and print the cards here.\n\n")
        try
            generate_flashcards(stdout)
        catch
            println("What kind of wacky setup are you running? I can't even get onto stdout!")
            println("Better luck next time, I guess ...", raw" ¯\_(ツ)_/¯")
        end
    end

    println("Have a nice day! :)\n")
end

# Everything looks good. 😊

make_card_file(DEFAULT_DIR,
               "generated_flashcards.txt",
               julia_common_operators)

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

# Some more for ya, because I need something a little more
# difficult!

julia_common_ops_sans_assignments = [
                          julia_ops_syntax_0
                          julia_ops_exponentiation
                          julia_ops_bitshifts
                          julia_ops_fractions
                          julia_ops_multiplication
                          julia_ops_addition
                          julia_ops_syntax_1
                          julia_ops_syntax_2
                          julia_ops_syntax_3
                          julia_ops_comparisons
                          julia_ops_control_flow
                          julia_ops_pairs
                        # julia_ops_assignments
                          ]

make_card_file(DEFAULT_DIR,
               "generated_flashcards-sans-assignments.txt",
               julia_common_ops_sans_assignments)


text:
"""
Module for performing Scalpels (Self-Correlation Analysis of Line Profiles for Extracting Low-amplitude Shifts)
based on a CCF timeseries.

For algorithm information, see Collier-Cameron, Ford, Shahaf et al. 2020

Author: Eric Ford
Date:   September 2020
"""
module Scalpels

using Statistics, StatsBase
using LinearAlgebra

default_max_num_basis = 16
default_num_basis = 4

include("scalpels_code.jl")
export clean_rvs_scalpels, calc_clean_rvs_scores_basis_scalpels, rms_clean_rvs_vs_num_basis_scalpels

end


text:
module NodeNets
include("DBFs.jl"); using .DBFs;
include("PointArrays.jl"); using .PointArrays;
using ..RealNeuralNetworks.SWCs
import LinearAlgebra: norm

export NodeNet

import LightGraphs
import BigArrays
using SparseArrays
import DataStructures: IntSet
using NearestNeighbors

const OFFSET = (zero(UInt32), zero(UInt32), zero(UInt32))
# rescale the skeleton
const EXPANSION = (one(UInt32), one(UInt32), one(UInt32))
# control the removing points around path based on DBF
const REMOVE_PATH_SCALE = 3
const REMOVE_PATH_CONST = 4

mutable struct NodeNet{T}
    # x,y,z,r
    nodeList            :: Vector{NTuple{4,T}}
    nodeClassList       :: Vector{UInt8}
    # connectivity matrix to represent edges
    # conn[2,3]=true means node 2 and 3 connected with each other
    # conn[3,2] should also be true since this is undirected graph
    connectivityMatrix  :: SparseMatrixCSC{Bool,UInt32}
end

@inline function NodeNet(nodeList::Vector{NTuple{4,T}},
                         connectivityMatrix::SparseMatrixCSC{Bool, UInt32}) where {T}
    nodeClassList = zeros(UInt8, length(nodeList))
    NodeNet{T}(nodeList, nodeClassList, connectivityMatrix)
end

"""
    NodeNet( seg, obj_id; penalty_fn=alexs_penalty)
Perform the teasar algorithm on the passed binary array.
"""
function NodeNet( seg::Array{T,3};
                     obj_id::T = convert(T,1),
                     expansion::NTuple{3, UInt32} = EXPANSION,
                     penalty_fn::Function = alexs_penalty ) where T
    # note that the object voxels are false and non-object voxels are true!
    # bin_im = DBFs.create_binary_image( seg, obj_id )
    points = PointArrays.from_seg(seg; obj_id=obj_id)
    NodeNet(points; expansion=expansion, penalty_fn=penalty_fn)
end

"""
    NodeNet(bin_im)
Parameters:
    bin_im: binary mask. the object voxel should be false, non-object voxel should be true
Return:
    nodeNet object
"""
function NodeNet(bin_im::Union{BitArray, Array{Bool,3}};
                 offset::NTuple{3, UInt32} = OFFSET,
                 expansion::NTuple{3, UInt32} = EXPANSION,
                 penalty_fn::Function = alexs_penalty)
        # transform segmentation to points
    points = PointArrays.from_binary_image(bin_im)

    println("computing DBF");
    # boundary_point_indexes = PointArrays.get_boundary_point_indexes(points, seg; obj_id=obj_id)
    #@time DBF = DBFs.compute_DBF( points, boundary_point_indexes );
    @time DBF = DBFs.compute_DBF(points)
    # @time dbf = DBFs.compute_DBF(points, bin_im)

    PointArrays.add_offset!(points, offset)
    NodeNet(points; dbf=dbf, penalty_fn=penalty_fn, expansion = expansion)
end

"""
    NodeNet( points; penalty_fn = alexs_penalty )

  Perform the teasar algorithm on the passed Nxd array of points
"""
function NodeNet( points::Matrix{T}; dbf::DBF=DBFs.compute_DBF(points),
                         penalty_fn::Function = alexs_penalty,
                         expansion::NTuple{3, UInt32} = EXPANSION) where T
    @assert length(dbf) == size(points, 1)
    println("total number of points: $(size(points,1))")
    points, bbox_offset = translate_to_origin!( points );
    # volumeIndex2NodeIndex represent points as a sparse vector containing the whole volume
    # in this way, we can fetch the node directly according to the coordinate.
    # the coordinate should be transfered to vector index though
    volumeIndex2NodeId, max_dims = create_node_lookup( points );
    max_dims_arr = [max_dims...];#use this for rm_nodes, but ideally wouldn't
    # transfer the coordinate to node
    # sub2node = x -> volumeIndex2NodeIndex[ sub2ind(max_dims, x[1],x[2],x[3]) ];#currently only used in line 48

    println("making graph (2 parts)");
    @time G, weights = make_neighbor_graph( points, volumeIndex2NodeId, max_dims;)
    println("build dbf weights from penalty function ...")
    @time dbf_weights = penalty_fn( weights, dbf, G )

    #init
    #nonzeros SHOULD remove duplicates, but it doesn't so
    # I have to do something a bit more complicated
    _,nonzero_vals = findnz(volumeIndex2NodeId);
    disconnectedNodeIdSet = IntSet( nonzero_vals );
    pathList = Vector(); # holds vector of nodeNet paths
    destinationNodeIdList = Vector{Int}(); #host dest node for each path
    # set of nodes for which we've "inspected" already
    # removing their neighbors based on DBF
    inspectedNodeIdList = Set{Int}();

    println("Finding paths")
    while length(disconnectedNodeIdSet) > 0
        rootNodeId = find_new_root_node_id( dbf, disconnectedNodeIdSet );
        @assert rootNodeId in disconnectedNodeIdSet

        #do a graph traversal to find farthest node and
        # find reachable nodes
        dsp_euclidean = LightGraphs.dijkstra_shortest_paths(G, rootNodeId, weights);
        #and another to find the min DBF-weighted paths to all nodes
        dsp_dbf       = LightGraphs.dijkstra_shortest_paths(G, rootNodeId, dbf_weights);

        #remove reachable nodes from the disconnected ones
        reachableNodeIdList = findall(.!(isinf.(dsp_euclidean.dists)));
        #another necessary precaution for duplicate nodes
        reachableNodeIdList = intersect(reachableNodeIdList, disconnectedNodeIdSet);
        setdiff!(disconnectedNodeIdSet, reachableNodeIdList);
        empty!(inspectedNodeIdList);

        while length(reachableNodeIdList) > 0

            #find the node farthest away from the root
            # by euc distance
            _, farthestNodeIndex = findmax( dsp_euclidean.dists[[reachableNodeIdList...]] );
            farthestNodeId = reachableNodeIdList[farthestNodeIndex];
            push!(destinationNodeIdList, farthestNodeId);
            println("dest node index: $(farthestNodeId)")

            if farthestNodeId == rootNodeId break end #this can happen apparently

            new_path = LightGraphs.enumerate_paths( dsp_dbf, farthestNodeId );

            push!(pathList, new_path)
            #can't do this in-place with arrays
            #this fn call is getting ridiculous
            @time reachableNodeIdList = remove_path_from_rns!( reachableNodeIdList,
                                                        new_path, points, volumeIndex2NodeId,
                                                        dbf, max_dims_arr,
                                                        inspectedNodeIdList );
        end #while reachable nodes from root
    end #while disconnected nodes

    println("Consolidating Paths...")
    path_nodes, path_edges = consolidate_paths( pathList );
    node_radii = dbf[path_nodes];

    # build a new graph containing only the nodeNet nodes and edges
    nodes, edges = distill!(points, path_nodes, path_edges)

    conn = get_connectivity_matrix(edges)
    nodeList = Vector{NTuple{4,Float32}}()
    sizehint!(nodeList, length(node_radii))
    for i in 1:length(node_radii)
        push!(nodeList, (map(Float32,nodes[i,:])..., node_radii[i]))
    end

    nodeNet = NodeNet(nodeList, conn)
    # add the offset from shift bounding box function
    bbox_offset = map(Float32, bbox_offset)
    @show bbox_offset
    add_offset!(nodeNet, bbox_offset)
    return nodeNet
end

function NodeNet(swc::SWC)
    nodeList = Vector{NTuple{4, Float32}}()
    nodeClassList = Vector{UInt8}()
    I = Vector{UInt32}()
    J = Vector{UInt32}()
    for (index, point) in enumerate(swc)
        push!(nodeList, (point.x, point.y, point.z, point.radius))
        push!(nodeClassList, point.class)
        if point.parent != -one(typeof(point.parent))
            # the connectivity matrix is symetric
            push!(I, index)
            push!(J, point.parent)
            push!(I, point.parent)
            push!(J, index)
        end
    end
    connectivityMatrix = sparse(I,J,true, length(swc), length(swc))
    NodeNet(nodeList, nodeClassList, connectivityMatrix)
end

function SWCs.SWC(nodeNet::NodeNet)
    edges = NodeNets.get_edges(nodeNet)
    swc = SWC()
    sizehint!(swc, NodeNets.get_node_num(nodeNet))
    nodeClassList = get_node_class_list(nodeNet)

    for (i, node) in NodeNets.get_node_list(nodeNet) |> enumerate
        point = SWCs.PointObj(nodeClassList[i], node[1], node[2], node[3], node[4], -1)
        push!(swc, point)
    end
    # assign parents according to edge
    for e in edges
        swc[e[2]].parent = e[1]
    end
    swc
end

##################### properties ###############################
@inline function get_node_list(self::NodeNet) self.nodeList end
@inline function get_connectivity_matrix(self::NodeNet) self.connectivityMatrix end
@inline function get_node_class_list(self::NodeNet) self.nodeClassList end
@inline function get_xyz(self::NodeNet) map(x->x[1:3], self.nodeList) end
@inline function get_radii(self::NodeNet) map(x->x[4],  self.nodeList) end
@inline function get_node_num(self::NodeNet; class::Union{Nothing, UInt8}=nothing)
    if class == nothing
        return length(self.nodeList)
    else
        return count(self.nodeClassList .== class)
    end
end

function node_list_to_array(nodeList::Vector{NTuple{4,T}}) where T
    N = length(nodeList)
    ret = Array{T,2}(undef, 4, N)
    for (i,node) in enumerate(nodeList)
        ret[:,i] = [node...]
    end
    ret
end

"""
    get_node_array(self::NodeNet)

Return:
    nodeArray::Array{T,2}: size is (nd, np), nd is the dimention of each node,
                np is the number of nodes
"""
@inline function get_node_array(self::NodeNet)
    nodeList = get_node_list(self)
    node_list_to_array(nodeList)
end

"""
the connectivity matrix is symmetric, so the connection is undirected
"""
@inline function get_edge_num(self::NodeNet) div(nnz(self.connectivityMatrix), 2) end

function get_edges(self::NodeNet)
    edges = Vector{Tuple{UInt32,UInt32}}()
    conn = get_connectivity_matrix(self)
    I,J,V = findnz(conn)
    # only record the triangular part of the connectivity matrix
    for index in 1:length(I)
        if I[index] > J[index]
            push!(edges, (I[index], J[index]) )
        end
    end
    edges
end

@inline function is_terminal_node(self::NodeNet, i::Integer)
    connectivityMatrix = get_connectivity_matrix(self)
    dropzeros!(connectivityMatrix)
    return nnz(connectivityMatrix[i,:]) < 2
end

function get_terminal_node_id_list(self::NodeNet)
    terminalNodeIdList = Vector{Int}()
    connectivityMatrix = get_connectivity_matrix(self)
    dropzeros!(connectivityMatrix)
    for i in 1:get_node_num(self)
        if nnz(connectivityMatrix[i,:]) < 2
            push!(terminalNodeIdList, i)
        end
    end
    return terminalNodeIdList
end

"""
    get_segment_point_num(self::NodeNet)

get number of branching points
"""
function get_num_segment_point(self::NodeNet)
    conn = get_connectivity_matrix(self)
    num_segment_point = 0
    for col in size(conn, 2)
        if nnz(conn[:,col]) > 2
            num_segment_point += 1
        end
    end
    return num_segment_point
end

"""
assume that the graph is acyclic, no loop.
"""
@inline function get_num_segmentes(self::NodeNet)
    get_num_segment_point(self) + 2
end

"""
    get_sholl_number(self::NodeNet, radius::AbstractFloat)
get the number of points which is in neurite and incounters with a sphere centered on root node
"""
function get_sholl_number(self::NodeNet, radius::AbstractFloat; rootNodeIndex::Integer=1)
    error("depracated, use the function in Neurons.")
    shollNum = 0
    nodeList = get_node_list(self)
    rootNode = nodeList[ rootNodeIndex ]
    for edge in get_edges(self)
        node1 = nodeList[ edge[1] ]
        node2 = nodeList[ edge[2] ]
        r1 = distance_from( rootNode, node1 )
        r2 = distance_from( rootNode, node2 )
        if r1>radius != r2>radius
            shollNum += 1
        end
    end
    shollNum
end

@inline function Base.length(self::NodeNet)
    length(self.nodeList)
end

@inline function Base.getindex(self::NodeNet, i::Integer)
    get_node_list(self)[i]
end

function Base.UnitRange(self::NodeNet)
    minCoordinates = [typemax(UInt32), typemax(UInt32), typemax(UInt32)]
    maxCoordinates = [zero(UInt32), zero(UInt32), zero(UInt32)]
    for node in get_node_list(self)
        minCoordinates = map(min, minCoordinates, node[1:3])
        maxCoordinates = map(max, maxCoordinates, node[1:3])
    end
    return [minCoordinates[1]:maxCoordinates[1],
            minCoordinates[2]:maxCoordinates[2],
            minCoordinates[3]:maxCoordinates[3]]
end

"""
    find_closest_node_id(self::NodeNet{T}, point::NTuple{3,T}) where T

look for the id of the closest node
"""
@inline function find_closest_node_id(self::NodeNet{T}, point::NTuple{N,T}) where {N,T}
    @assert N>=3 && N<=4
    find_closest_node_id(self, [point[1:3]...])
end

function find_closest_node_id(self::NodeNet{T}, point::Vector{T}) where T
    nodeArray = get_node_array(self)
    kdtree = KDTree(nodeArray; leafsize=10)
    idxs, _ = knn(kdtree, point, 1)
    return idxs[1]
end

##################### transformation ##########################
"""
get binary buffer formatted as neuroglancer nodeNet.

# Binary format
    UInt32: number of vertex
    UInt32: number of edges
    Array{Float32,2}: Nx3 array, xyz coordinates of vertex
    Array{UInt32,2}: Mx2 arrray, node index pair of edges
reference:
https://github.com/seung-lab/neuroglancer/wiki/Skeletons
"""
function get_neuroglancer_precomputed(self::NodeNet)
    @show get_node_num(self)
    @show get_edge_num(self)
    # total number of bytes
    num_bytes = 4 + 4 + 4*3*get_node_num(self) + 4*2*get_edge_num(self)
    buffer = IOBuffer( read=false, write=true, maxsize=num_bytes )
    # write the number of vertex, and edges
    write(buffer, UInt32(get_node_num(self)))
    write(buffer, UInt32(get_edge_num(self)))
    # write the node coordinates
    for node in get_node_list(self)
        write(buffer, [node[1:3]...])
    end
    # write the edges
    for edge in get_edges( self )
        # neuroglancer index is 0-based
        write(buffer, UInt32( edge[1] - one(UInt32) ))
        write(buffer, UInt32( edge[2] - one(UInt32) ))
    end
    bin = Vector{UInt8}(take!(buffer))
    close(buffer)
    return bin
end

"""
    get_connectivity_matrix(edges::Vector)

construct sparse connectivity matrix accordint to the edges
"""
@inline function get_connectivity_matrix(edges::Vector)
    I = [e[1] for e in edges]
    J = [e[2] for e in edges]
    return sparse([I...,J...],[J...,I...],ones(Bool, 2*length(edges)))
end

###################### IO #################################
function save(self::NodeNet, fileName::AbstractString)
    if endswith(fileName, ".swc") || endswith(fileName, ".swc.bin")
        swc = SWC(self)
        SWCs.save(swc, fileName)
    else
        error("only support format of swc and swc.bin, but getting: ", fileName)
    end
end

"""
    save(self::NodeNet, cellId::UInt32, d_json::Associative, d_bin::Associative)

save nodeNet in google cloud storage for neuroglancer visualization
the format is the same with meshes
"""
function save(self::NodeNet, cellId::Integer, d_bin::AbstractDict)
    # get the bounding box of nodeNet and transfer to string representation
    # example string: 1432-1944_1264-1776_16400-16912
    rangeString = BigArrays.Indexes.unit_range2string( UnitRange(self) )
    # write the binary representation of nodeNet
    d_bin["$(cellId)"] = get_neuroglancer_precomputed(self)
end

"""
save binary file of point pairs
used in neuroglancer python interface to visualize the nodeNet
"""
function save_edges(self::NodeNet, fileName::String)
    open(fileName, "w") do f
        nodeList = get_node_list(self)
        edges = get_edges(self)
        for edge in edges
            write(f, nodeList[ edge[1] ])
            write(f, nodeList[ edge[2] ])
        end
    end
end

##################### manipulate ############################
@inline function add_offset!(self::NodeNet{T}, offset::Union{Vector{T},NTuple{3,T}} ) where T
    @assert length(offset) == 3
    for i in 1:get_node_num(self)
        xyz = map((x,y)->x+y, self.nodeList[i][1:3],offset)
        self.nodeList[i] = (xyz..., self.nodeList[i][4])
    end
end

@inline function stretch_coordinates!(self::NodeNet, mip::Real)
    expansion = (2^(mip-1), 2^(mip-1), 1)
    stretch_coordinates!(self, expansion)
end
@inline function stretch_coordinates!(self::NodeNet, expansion::Union{Vector, Tuple})
    @assert length(expansion) == 3
    nodeList = get_node_list(self)
    for (nodeIndex, node) in enumerate(nodeList)
        self.nodeList[ nodeIndex ] = (map((x,y)->x*y, node[1:3], expansion)...,
                                      node[4]* (prod(expansion)^(1/3)))
    end
end

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

"""

    translate_to_origin!( points::Matrix{T} )

  Normalize the point dimensions by subtracting the min
  across each dimension. This step isn't extremely necessary,
  but might be useful for compatibility with the MATLAB code.
  record the offset and add it back after building the nodeNet
"""
function translate_to_origin!( points::Matrix{T} ) where T
    offset = minimum( points, dims=1 ) .- one(T) ;
    # transform to 1d vector
    offset = vec( offset )
    @assert length(offset) == 3
    @assert offset[3] < T(20000)
    points[:,1] .-= offset[1]
    points[:,2] .-= offset[2]
    points[:,3] .-= offset[3]

    points, offset
end

"""

    create_node_lookup( points )

  Abstractly represent the points as a volume by means of
  linear indexes into a sparse vector.
"""
function create_node_lookup( points::Matrix{T} ) where T

  #need the max in order to define the bounds of that volume
  # tuple allows the point to be passed to fns as dimensions
  max_dims = ( maximum( points, dims=1 )... ,);
  max_dims = map(Int, max_dims)
  num_points = size(points,1);
  #creating the sparse vector mapping indices -> node ids
  volumeIndex2NodeIndex = sparsevec( [
                        (LinearIndices(max_dims))[Vector{UInt32}(points[i,:])...,]
                            for i=1:num_points ], 1:num_points, prod(max_dims), (x,y) -> x )

  volumeIndex2NodeIndex, max_dims
end

"""

    make_neighbor_graph( points )

  Creates a LightGraphs graph from a point cloud, such that
  each point has a node (indexed in the same order as the rows),
  and each node is connected to its neighbors (by 26-connectivity).

  Also returns a sparse matrix of weights for these edges equal to the
  euclidean distance between the indices of each point. These can
  be weighted or modified easily upon return.
"""
function make_neighbor_graph( points::Array{T,2}, volumeIndex2NodeIndex=nothing, max_dims=nothing;
                             expansion::NTuple{3, UInt32} = EXPANSION ) where T

  if volumeIndex2NodeIndex == nothing volumeIndex2NodeIndex, max_dims = create_node_lookup(points) end

  #init
  num_points = size(points,1);
  g = LightGraphs.Graph(num_points);

  #26-connectivity neighborhood
  # weights computed by euc_dist to center voxel
  nhood = [[i,j,k] for i=1:3,j=1:3,k=1:3];

  map!( x-> (x .- [2,2,2]).*[expansion ...] , nhood, nhood)
  nhood_weights = map( norm, nhood );

  #only adding weights for non-duplicate nodes
  _,non_duplicates = findnz( volumeIndex2NodeIndex );

  #init - descriptors of the edges for graph
  # weights later
  is = Vector{Int}();
  js = Vector{Int}();
  ws = Vector{Float64}();

  dest_sub = Array{Int}(undef, 3)

  for i in eachindex(nhood)

  	#only need to explicitly tell LightGraphs to add the edges in
  	# one direction (undirected), though doing this properly
    # can be tricky
    #direction = ind2sub((3,3,3),i)
    direction = Tuple(CartesianIndices((3,3,3))[i])

    if (direction[1] == 3 ||
       (direction[1] == 2 && 2*direction[2] + direction[3] >= 6)) continue end

    for p in non_duplicates

      #location of potential node to connect
      for j=1:3 dest_sub[j] = nhood[i][j] + points[p,j]' end

      #bounds_checking
      if (dest_sub[1] < 1 ||
      	  dest_sub[2] < 1 ||
      	  dest_sub[3] < 1 ) continue end

      if (dest_sub[1] > max_dims[1] ||
      	  dest_sub[2] > max_dims[2] ||
      	  dest_sub[3] > max_dims[3] ) continue end

      #fetching the id of the proposed node location
      dest_ind = (LinearIndices(max_dims))[dest_sub[1:3]...,]
      dest_node = volumeIndex2NodeIndex[dest_ind];

      if dest_node != 0 #valid edge
        LightGraphs.add_edge!(g,p,dest_node)
        push!(is,p); push!(js,dest_node); push!(ws,nhood_weights[i]);
      end

    end #for p
  end #for i in nhood

  #creating weight matrix
  g_weights = sparse(vcat(is,js),vcat(js,is),vcat(ws,ws), num_points, num_points);

  g,g_weights
end

"""

    alexs_penalty( weights, dbf )

  Returns a version of the edge weights, modified by the DBF-based
  teasar penalty:

  w = w * 5000 .* (1 - dbf/maxDBF).^(16)

  The factor of 5000 shouldn't make a difference, but I've
  left it here
"""
function alexs_penalty( weights, dbf, G, mult=true )

  # *1.01 ensures later quotient is >0
  M = maximum(dbf, dims=1)*1.01

  p_v = (1 .- dbf./M ) .^ 16;

  #DBF penalty is defined by the destination node, only need js
  is, js = begin
      I = findall(!iszero, weights)
      (getindex.(I, 1), getindex.(I, 2))
  end
  ws = nonzeros(weights);

  if mult
    new_ws = ws .* p_v[js];
  else
    new_ws = ws .+ p_v[js];
  end

  sparse( is, js, new_ws, size(weights)... );
end

"""

    local_max_multiplicative_penalty( weights, dbf, G )

  Returns a version of the edge weights, modified by the DBF-based
  penalty

  w = w * (1 - DBFdest/DBFstar)

  Where DBFstar is the maximum DBF value of an outwards neighbor for
  each node
"""
function local_max_multiplicative_penalty( weights, dbf, G )

  ns = LightGraphs.vertices(G)
  dbf_star = Vector{Float64}( length(ns) );

  #this is ok since ns is always an int range 1:N
  for n in ns
    neighborhood = LightGraphs.out_neighbors(G, n);

    if length(neighborhood) > 0
      dbf_star[n] = maximum(dbf[neighborhood]);
    else
      dbf_star[n] = 1
    end
  end

  is, js = findn(weights)
  ws = nonzeros(weights)
  new_ws = Vector{Float64}(length(ws));

  for i in eachindex(is)
    dbf_dest = dbf[js[i]];
    dbf_s = dbf_star[is[i]];

    new_ws[i] = ws[i] * (1 - dbf_dest/dbf_s);
  end

  sparse( is, js, new_ws, size(weights)... );
end

"""

    find_new_root_node( dbf )

  Extracts the point with the largest dbf from the
  Array of passed points
"""
@inline function find_new_root_node_id(dbf::DBF, disconnectedNodeIdSet::IntSet)
    disconnectedNodeIndexList = [disconnectedNodeIdSet...]
    disconnectedNodeDBFList = dbf[disconnectedNodeIndexList]
    rootIndexInDisconnectedList = argmax( disconnectedNodeDBFList )
    rootIndex = disconnectedNodeIndexList[rootIndexInDisconnectedList]
    return rootIndex
end

"""

    find_new_root_node( points )

  Extracts the point with the lowest linear index from the
  Array of passed points
"""
@inline function find_new_root_node_V1( points::Array{T,2} ) where T

  res = points;
  for i in 3:-1:1
    res = filter_rows_by_min( res, i );
  end

  res
end

function filter_rows_by_min( arr::Array{T,2}, dim ) where T

  selected_rows = arr[:,dim] .== minimum(arr,1)[dim];

  arr[selected_rows,:];
end

function create_dummy_matrix( s )
  points = hcat( findn(ones(s,s,s))... );
  i = div(size(points,1),2)+1;

  points = points[[1:i-1;i+1:end],:];
  points;
end

"""

    remove_path_from_rns!( reachableNodeList, path, points, sub2node, dbf, max_dims, scale_param, const_param)

  Identifies the nodes to remove from the reachable nodes within the graph.
  Probably the ugliest function here, and should be optimized later

  TO OPTIMIZE
"""
function remove_path_from_rns!( reachableNodeList::Vector, path::Vector{Int},
                points, volumeIndex2NodeIndex, dbf, max_dims, inspectedNodeIdList,
                scale_param::Integer=REMOVE_PATH_SCALE,
                const_param::Integer=REMOVE_PATH_CONST )
  radii = dbf[path].*scale_param .+ const_param;

  to_remove = Set{Int}();
  for i in eachindex(path)
    nodeIndex = path[i];
    if !(nodeIndex in inspectedNodeIdList)
        push!(inspectedNodeIdList, nodeIndex);
	    union!( to_remove, nodes_within_radius(points[nodeIndex,:], volumeIndex2NodeIndex,
                                               radii[i], max_dims ) );
    end
  end

  return setdiff(reachableNodeList, to_remove);
end

"""

    nodes_within_radius( sub, sub2node, r, max_dims )

  Identifies the node indices within the subscript
"""
function nodes_within_radius( sub::Array{T,1}, volumeIndex2NodeIndex, r, max_dims::Vector ) where T;

  beginning = convert(Vector{Int}, ceil.(max.(sub[:] .- r,1)));
  ending    = convert(Vector{Int}, floor.(min.(sub[:] .+ r, max_dims)));
  ind::Int = convert(Int,0);
    max_dims = Vector{Int}(max_dims)

  nodes = Set{T}();
  for x in beginning[1]:ending[1]
    for y in beginning[2]:ending[2]
      for z in beginning[3]:ending[3]
        ind = x + (y-1)*max_dims[1] + (z-1)*max_dims[2]*max_dims[1];
        push!( nodes, volumeIndex2NodeIndex[ind] );
        #don't need to worry about the '0 node' here as
        # the only result that matters is the setdiff from
        # the reachable nodes
      end
    end
  end

  nodes
end

"""

    consolidate_paths( path_list )

  Extracts the unique nodes and edges among all paths in the list
"""
function consolidate_paths( path_list::Vector )

    nodes = Set{UInt32}();
    edges = Set{NTuple{2,UInt32}}();

    # estimate the total number of nodes and edges
    num_nodes = 0
    num_edges = 0
    for path in path_list
        path_length = length(path)
        num_nodes += path_length
        num_edges += path_length-1
    end
    # give some hints for memory allocation
    sizehint!(nodes, num_nodes)
    sizehint!(edges, num_edges)

    for path in path_list
        path_length = length(path);
        @assert path_length > 1

        for i in 1:(path_length-1)
            @assert path[i] != path[i+1]
            push!( edges, (path[i], path[i+1]) );
            push!( nodes, path[i] );
        end
        push!( nodes, path[path_length] );
    end

    collect(nodes), collect(edges)
end

"""
Parameters:
===========
    point_array: a Nx3 array recording all the voxel coordinates inside the object.
    path_nodes: the indexes of nodeNet voxels in the point_array
    path_edges: the index pairs of nodeNet in the point_array

In the path_nodes, the nodes were encoded as the index of the point_array. Since we do
not need point_array anymore, which could be pretty big, we'll only reserve the nodeNet
coordinates and let the gc release the point_array.
the same applys to path_edges
"""
function distill!(point_array::Array{T,2},
                  path_nodes::Vector, path_edges::Vector) where T
    num_nodes = length(path_nodes)
    num_edges = length(path_edges)
    nodes = Array{T}(undef, num_nodes, 3)
    # map the old path node id to new id
    id_map = Dict{T, T}()
    sizehint!(id_map, num_nodes)
    # build new nodes and the id map
    for i in 1:num_nodes
        nodes[i,:] = point_array[path_nodes[i], :]
        id_map[path_nodes[i]] = i
    end

    # rebuild the edges
    for i in 1:num_edges
        path_edges[i] = (id_map[path_edges[i][1]], id_map[path_edges[i][2]])
    end
    for edge in path_edges
        # no self connecting edge!
        @assert edge[1] != edge[2]
    end
    return nodes, path_edges
end

end # module


text:
@osx_only begin
    # OS X only for the moment, because TestImages is only
    # required for OS X (see test/REQUIRE).  That seems to cause
    # a testing failure on Travis for Linux.
    # I assume that at some point, streaming will expand and other tests
    # will be needed, at which point we can open up this test
    using Images, TestImages
    using Base.Test

    workdir = joinpath(tempdir(), "ImagesStreamTest")
    testfile = joinpath(workdir,"jigsaw_tmpl.png")
    if !isfile(testfile)
        mkpath(workdir)
        download("http://www.imagemagick.org/Usage/images/jigsaw_tmpl.png",testfile)
    end

    open(testfile,"r") do io
        img = imread(io,Images.OSXNative)
    end

    ## cleanup
    rm(workdir,recursive=true)
end


text:
import BasicModelInterface as BMI
import ModflowInterface as MF

mf6_modelname = "MF6LAKE"
cd("examples/mf6lake")

m = BMI.initialize(MF.ModflowModel)
BMI.get_component_name(m)

##

MF.prepare_time_step(m, 0.0)
MF.prepare_solve(m, 1)

##

function solve_to_convergence(model)
    converged = false
    iteration = 1
    while !converged
        @show iteration
        converged = MF.solve(model, 1)
        iteration += 1
    end
    println("converged")
    return
end

solve_to_convergence(m)

##

headtag = MF.get_var_address(m, "X", mf6_modelname)
head = BMI.get_value_ptr(m, headtag)

##

# Note, this de-allocates the head array!
BMI.finalize(m)


text:
using Turing, BridgeSampling

## Test on bivariate standard Normal distribution
samples = rand(MvNormal(ones(2)), 10_000)

log_posterior(x) = -0.5 * x' * x

Bridge = bridgesampling(samples, log_posterior, [-Inf, -Inf], [Inf, Inf])
Analytical = log(2*pi)

isapprox(value(Bridge), Analytical, atol = 1e-3)

## Hierarchical Normal example

# Generate data
true_μ = 0.0
true_τ² = 0.5
true_σ² = 1.0
n = 20
true_θ = rand(Normal(true_μ, sqrt(true_τ²)), n)
#y = rand.(Normal.(true_θ, sqrt(true_σ²)))
# Same values as in : https://cran.r-project.org/web/packages/bridgesampling/vignettes/bridgesampling_example_jags.html
y = vec([1.19365332  1.95745331 -0.72161754 -1.87380833 -1.16928239  0.51960853 -0.03610041  0.42508815  0.41119221 -0.81236980  0.72967357  3.48186722  2.31126381  2.00029422 -0.27643507  1.06882370 -0.95083599 -1.89651101  2.56019737  0.23703060])

# DynamicPPL models
@model H0(Y) = begin
    τ² ~ InverseGamma(1,1)
    θ ~ filldist(Normal(0.0, sqrt(τ²)), length(Y))
    Y ~ arraydist(Normal.(θ, 1.0))
end

@model H1(Y) = begin
    μ ~ Normal()
    τ² ~ InverseGamma(1,1)
    θ ~ filldist(Normal(μ, sqrt(τ²)), length(Y))
    Y ~ arraydist(Normal.(θ, 1.0))
end

M0 = H0(y)
M1 = H1(y)

# Sampling
chn0 = sample(M0, NUTS(), MCMCThreads(), 50_000, 3, burnin=2000)
chn1 = sample(M1, NUTS(), MCMCThreads(), 50_000, 3, burnin=2000)

# Estimate log marginal likelihood : directly with the DynamicPPL model (lower precision)
L0 = bridgesampling(chn0, M0)
L1 = bridgesampling(chn1, M1)

# Or with specified log posterior functions and boundaries
log_posterior_H0(x) = sum([logpdf(Normal(x.θ[i], 1.0), y[i]) for i = 1:n]) +
                sum([logpdf(Normal(0.0, sqrt(x.τ²)), x.θ[i]) for i = 1:n]) +
                logpdf(InverseGamma(1,1), x.τ²)

log_posterior_H1(x) = sum([logpdf(Normal(x.θ[i], 1.0), y[i]) for i = 1:n]) +
                sum([logpdf(Normal(x.μ, sqrt(x.τ²)), x.θ[i]) for i = 1:n]) +
                logpdf(InverseGamma(1,1), x.τ²) +
                logpdf(Normal(), x.μ)

samplesH0, names = BridgeSampling.extract_samples(chn0, M0)
L0 = bridgesampling(samplesH0, log_posterior_H0, [0.0, -Inf], [Inf, Inf]; names=names)

samplesH1, names = BridgeSampling.extract_samples(chn1, M1)
L1 = bridgesampling(samplesH1, log_posterior_H1, [-Inf, 0.0, -Inf], [Inf, Inf, Inf]; names=names)

# Bayes factor
bayes_factor(L0, L1)

# Posterior model probabilities
posterior_probabilities([L0, L1])


text:
"""
    unitRandom3Cartesian(rng = MersenneTwister())

Returns a unit vector pointing in a random direction in Cartesian coordinates.
"""
function unitRandom3Cartesian(rng = MersenneTwister())
    θ = π*rand(rng)
    ϕ = 2.0*π*rand(rng)
    return Vector([
        cos(ϕ)*sin(θ),
        sin(ϕ)*sin(θ),
        cos(θ)
    ])
end

export unitRandom3Cartesian

"""
    randbtw(a::T,b::T,rng=MersenneTwister()) where {T<:Real}

Generates a random number between a and b.
"""
function randbtw(a::T,b::T,rng=MersenneTwister()) where {T<:Real}
    if a > b
        return randbtw(b,a,rng)
    end
    return a+rand(rng)*(b-a)
end

export randbtw
"""
    dist(n,func,domain,range,rng=MersenneTwister())

Generates n values following the distribution func in a given domain and range
    - SLOW
    - Can be used on any distribution function
"""
function dist(n,func,domain,range,rng=MersenneTwister())
    out = zeros(n)
    for i=1:n
        val = 0.0
        tryagain = true
        while tryagain
            xtest = randbtw(domain...,rng)
            fxtest = abs(func(xtest))
            if abs(randbtw(range...,rng)) < fxtest
                val = xtest
                tryagain = false
            end
        end
        out[i] = val
    end
    return out
end

export dist

"""
    PowerLaw

Sampler for a power law distribution
"""
struct PowerLaw <: Sampleable{Univariate, Continuous}
    α::Real
    domain::Tuple{T,T} where {T<:Real}
end

export PowerLaw

"""
    rand(rng::AbstractRNG,s::PowerLaw)

Generate one sample from a power law distribution
"""
function Base.rand(rng::AbstractRNG,s::PowerLaw)
    ξ = rand(rng)
    γ = 1.0-s.α
    e0γ= s.domain[1]^(γ)
    efγ = s.domain[2]^(γ)
    out = (efγ-e0γ)*ξ+e0γ
    return out^(1/γ)
end

"""
    BrokenPowerLaw

Sampler for a broken power law distribution
"""
struct BrokenPowerLaw <: Sampleable{Univariate, Continuous}
    α::Array{T} where {T<:Real}
    domain::Array{T} where {T<:Real}
    function BrokenPowerLaw(α::Array{T}, domain::Array{T}) where {T<:Real}
        @assert length(α)+1 == length(domain) "Requires length(α) == length(domain)-1"
        bOrder = true
        for i=1:(length(domain)-1)
            bOrder *= domain[i] < domain[i+1]
        end

        @assert bOrder "domain breaks must be in order"
        new(α, domain)
    end
end

export BrokenPowerLaw

function afterIndex(val,vec)
	for i=1:(length(vec)-1)
		check = vec[i] < val
		check *= val < vec[i+1]
		if check
			return i
		end
	end
	return 0
end

"""
    rand(rng::AbstractRNG,s::BrokenPowerLaw)

Generate one sample from a broken power law distribution.
"""
function Base.rand(rng::AbstractRNG,s::BrokenPowerLaw)
	ξ = rand(rng)
	γ = 1 .- s.α

	cx = ones(length(s.α))
	for i = 1:(length(cx)-1)
		β = s.α[i]/s.α[i+1]
		cx[i+1] = cx[i]^β*s.domain[i+1]^(β-1)
	end
	cx = map(/,cx,γ)

	arts = [
		begin
			(s.domain[i+1]^γ[i] - s.domain[i]^γ[i])*cx[i]
		end for i=1:length(γ)
			]

	parts = vcat([0],arts)
	sumParts = cumsum(parts)

	norm = sumParts[end]
	cnorm = ξ*norm

	supIndex = afterIndex(cnorm,sumParts)

	out = cnorm - sumParts[supIndex]
	out /= cx[supIndex]
	out += s.domain[supIndex]^γ[supIndex]
	out = out^(1/γ[supIndex])
	return out
end


text:
@newimmutable NewFoo{T} <: Integer
@newtype NewBar{T <: AbstractString} <: AbstractString

@newimmutable NewBaz <: UTF8String
@stringinterface NewBaz

@newtype NewQux <: UTF8String
@stringinterface NewBaz

@newtype NewFnord <: Int64

facts("New types") do
    @fact UserId("abc") --> UserId("abc")
    @fact UserId("abc") --> not(UserId("foo"))

    takes_id(a::UserId) = a
    @fact_throws MethodError takes_id(SlackName("abc"))

    newfoo = NewFoo{Int}(1)
    @fact typeof(newfoo.v) --> Int
    @fact newfoo.v         --> 1

    newbar = NewBar{UTF8String}(utf8("bar"))
    @fact typeof(newbar.v) --> UTF8String
    @fact newbar.v         --> utf8("bar")

    newbaz = NewBaz(utf8("baz"))
    @fact typeof(newbaz.v) --> UTF8String
    @fact_throws MethodError NewBaz(1)

    newqux = NewQux(utf8("qux"))
    @fact typeof(newqux.v) --> UTF8String
    @fact newqux.v --> utf8("qux")
    newqux.v = "qux2"
    @fact newqux.v --> utf8("qux2")
    @fact_throws MethodError NewBaz(1)

    @fact NewFnord(123).v --> 123
end


text:
# =============================================================================
# Bouncing ball with air friction
#
# system type: hybrid system with polynomial dynamics
# state dimension: 2
# modes: 2
# transitions: 2
#
# See Example 4.1.2 (page 98f) in [1].
#
# [1] Xin Chen: Reachability analysis of non-linear hybrid systems using Taylor
# models. PhD thesis, 2015.
# =============================================================================
using HybridSystems, MathematicalSystems, LazySets, TaylorIntegration

@taylorize function flow_down!(dx, x, params, t)
    dx[1] = x[2]
    dx[2] = -9.8 + 0.1*(x[2])^2
    return dx
end

@taylorize function flow_up!(dx, x, params, t)
    dx[1] = x[2]
    dx[2] = -9.8 - 0.1*(x[2])^2
    return dx
end

function bouncingBallNonlinear_model()
    # hybrid automaton with state variables x, v
    HA = LightAutomaton(2)

    # mode 1 ("down")
    X = HPolyhedron([HalfSpace([-1.0, 0.0], 0.0),   # x ≥ 0
                     HalfSpace([ 0.0, 1.0], 0.0)])  # v ≤ 0
    m1 = @system(x' = flow_down!(x), dim: 2, x ∈ X)

    # mode 2 ("up")
    X = HPolyhedron([HalfSpace([-1.0, 0.0], 0.0),   # x ≥ 0
                     HalfSpace([0.0, -1.0], 0.0)])  # v ≥ 0
    m2 = @system(x' = flow_up!(x), dim: 2, x ∈ X)

    # α transition down → up
    add_transition!(HA, 1, 2, 1)
    G = Hyperplane([1.0, 0.0], 0.0)  # x = 0
    A = [1.0 0.0; 0.0 -0.8]
    Rα = ConstrainedLinearMap(A, G)  # v := -0.8v

    # β transition up → down
    add_transition!(HA, 2, 1, 2)
    G = Hyperplane([0.0, 1.0], 0.0)  # v = 0
    Rβ = ConstrainedIdentityMap(2, G)

    # hybrid system
    S = HybridSystem(HA, [m1, m2], [Rα, Rβ], fill(AutonomousSwitching(), 2))

    return S
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
function bitArray2Str(x::Union{BitArray{1},Array{Bool,1}})
    parseString(b::Bool) = b == 1 ? "1" : "0"
    number = ""
    for i = 1:1:length(x)
        number = number * parseString(x[i])
    end
    return number
end

function intArray2Str(x::Array{Int64,1})
    parseString(n::Int64) = n == 1 ? "1" : "0"
    number = ""
    for i = 1:1:length(x)
        number = number * parseString(x[i])
    end
    return number
end

function bin2dec(buf)
    num = 0
    for i in 1:length(buf)
        if buf[length(buf) - i + 1]
            num += 2^(i - 1)
        end
    end
    return num
end

function bin2dec_twoscomp(x::Union{BitArray{1},Array{Bool,1}})
    if x[1] == false
        num = bin2dec(x)
    else
        _xor = Array{Bool,1}(trues(length(x)))
        _sum1 = [falses(length(x) - 1);true]
        x = map(⊻, x, _xor) # invert Bits
        x = map(+, x, _sum1) # sum 1. At this point x is no longer a BitArray but an Array{Int64,1}

        for i = length(x):-1:2
            if x[i] == 2
                x[i] = 0
                x[i - 1] = x[i - 1] + 1
            end
        end
        x[1] == 2 ? x[1] = 0 : x[1] = x[1]

        number = intArray2Str(x)
        num = -parse(Int, number, base = 2)
    end
    return num
end


text:
using Revise
using ADCME
using NNFEM
using JLD2
using PyCall
using LinearAlgebra
using Random, Distributions

reset_default_graph()

E = 200e3

nntype = "ae_scaled"

#nntype = "piecewise"

idx = 1
H_function = "None"
include("nnutil.jl")

n = 50
X_1d = LinRange(-2.0, 2.0, n)
X = zeros(Float64, n, 1)
for i = 1:n
    X[i, 1] = X_1d[i]
end

ε_μ = 0.0; ε_σ = 0.00
noise = rand(Normal(ε_μ, ε_σ), (n,1))
Y = X * E + noise

loss = constant(0.0)
x = constant(X)
dummy = x

#ε = x; ε0 = x; σ0 = x

y = nn(x, dummy, dummy)

loss += mean((y-Y)^2)

sess = Session(); init(sess)
@show run(sess, loss)
for i = 1:10
    BFGS!(sess, loss, 1000)
    ADCME.save(sess, "Data/NNLearn_$(nntype)_ite$(i).mat")
end


text:
module RemoteFilesystem

"""
    RemoteFS(host, path=""; ro::Bool=false)

Mount a remote filesystem located at `host:path` to the local filesystem.

## Examples
```julia
# Mount remote user home directory
RemoteFs("fredrik@remotehost")

# Mount remote directory /home/fredrik/files in read-only mode:
RemoteFS("fredrik@remotehost", "/home/fredrik/files"; ro=true)
"""
mutable struct RemoteFS
    host::String
    remotepath::String
    localpath::String
    function RemoteFS(host, path=""; ro::Bool=false)
        sshfs = Sys.which("sshfs")
        sshfs === nothing && error("sshfs not found")
        fusermount = Sys.which("fusermount")
        fusermount === nothing && error("fusermount not found")
        localpath = mktempdir(; cleanup=false)
        cmd = `$sshfs $host:$path $localpath`
        ro && (cmd = `$cmd -o ro`)
        run(cmd)
        rfs = new(host, path, localpath)
        function unmount_fs(rfs::RemoteFS)
            run(`$fusermount -u $(rfs.localpath)`)
            rm(rfs.localpath; recursive=true, force=true)
        end
        finalizer(unmount_fs, rfs)
        return rfs
    end
end

Base.joinpath(rfs::RemoteFS, path::String, paths::String...) = joinpath(rfs.localpath, path, paths...)
Base.readdir(rfs::RemoteFS) = readdir(rfs.localpath)

end # module


text:
function onedElliptic(nElements=40,order=2)
#  Solves a linear elliptic PDE in 1D

  f(x)=x           # right-hand side function definition.

  x,eConn,iU,iC = onedMesh( [0.0,1.0], order, nElements )
  nNodes = size(x,1)

  ide = [1:nNodes;]     # index into unknowns  (set b.c. later in 1D)

  r,w = onedQuadratureRule(5)
  o   = ones(size(r))

  nElDOF = order+1
  nEntries = nElements*nElDOF^2
  II  = zeros(Int64,nEntries) #Array{Int32,1}(undef,nEntries)     # zeros(Int32,nEntries)
  JJ  = zeros(Int64,nEntries) #Array{Int32,1}(undef,nEntries)     # zeros(Int32,nEntries)
  AA  = zeros(Float64,nEntries) #Array{Float64,1}(undef,nEntries)   # zeros(Float64,nEntries)
  b   = zeros(Float64,nNodes) #Array{Float64,1}(undef,nNodes)     # zeros(Float64,nNodes)

  index = 0
  for k=1:nElements
    nLocal = eConn[k,:]
    xLocal = view(x,nLocal)

    xg,wg,phi,p_x,p_xx = onedShape( xLocal, r, w )
    fg     = f(xg)             # forcing function evaluated at quadrature points

    ALoc   = onedBilinear(  o, p_x, p_x, wg )
    bLoc   = onedLinForm(  fg,      phi, wg )

    lDOF = ide[nLocal]
    for nt = 1:nElDOF
      for nu = 1:nElDOF
        index = index + 1
        II[index]   = lDOF[nu]
        JJ[index]   = lDOF[nt]
        AA[index]   = ALoc[nt,nu]
      end
      b[lDOF[nt]] = b[lDOF[nt]] + bLoc[nt]
    end
  end

  A = sparse(II,JJ,AA)

  u = A[iU,iU]\b[iU]

  #  Calculate the exact solution for unit test
  xi = x[iU];

  uExact = -(xi.^3-xi)/6

  return norm(u-uExact)
# using Gaston
# figure(1)
# Gaston.set_terminal("x11");
# Gaston.closeall();
# c = Gaston.CurveConf();
# c.color = "blue";
# c.legend = "FEM";
# Gaston.addcoords(xi,u,c);

# figure(2)
# Gaston.set_terminal("x11");
# c = Gaston.CurveConf();
# c.color = "red";
# c.legend = "Exact";
# c.marker = "fdmd";
# c.plotstyle = "linespoints";
# Gaston.addcoords(xi,-(xi.^3-xi)/6,c);
# Gaston.llplot();

# set_filename("test.png")
# printfigure("png")

#  using PyPlot
#  plot(xi,u,"r",label="FEM solution")
#  plot(xi,-(xi.^3-xi)/6,"go",label="Exact at nodes")
#  legend()
end


text:
let
	f5(x) = (x[1] - 5.0)^2 + (x[2] - 11.0)^2
	function g5!(x, storage)
		storage[1] = 2.0 * (x[1] - 5.0)
		storage[2] = 2.0 * (x[2] - 11.0)
	end
	d = DifferentiableFunction(f5, g5!)

	x = [0.0, 0.0]
	x_new = [0.0, 0.0]
	gr_new = [0.0, 0.0]
	phi0 = d.fg!(x, gr_new)
	p = -gr_new
	dphi0 = dot(p, gr_new)

	lsr = Optim.LineSearchResults(eltype(x))
	push!(lsr, 0.0, phi0, dphi0)

	alpha = 1.0
	mayterminate = false
	alpha, f_update, g_update = Optim.backtracking_linesearch!(d, x, p, x_new, gr_new, lsr, alpha, mayterminate)
	alpha, f_update, g_update = Optim.interpolating_linesearch!(d, x, p, x_new, gr_new, lsr, alpha, mayterminate)
	alpha, f_update, g_update = Optim.mt_linesearch!(d, x, p, x_new, gr_new, lsr, alpha, mayterminate)
	alpha, f_update, g_update = Optim.hz_linesearch!(d, x, p, x_new, gr_new, lsr, alpha, mayterminate)

	Optim.optimize(d, [0.0, 0.0], method=LBFGS())
end