text:
#########################################################################
# Wrapper around Python exceptions

type PyError <: Exception
    msg::AbstractString # message string from Julia context, or "" if none

    # info returned by PyErr_Fetch/PyErr_Normalize
    T::PyObject
    val::PyObject
    traceback::PyObject

    # generate a PyError object.  Should normally only be called when
    # PyErr_Occurred returns non-NULL, and clears the Python error
    # indicator.
    function PyError(msg::AbstractString)
        exc = Array{PyPtr}(3)
        pexc = convert(UInt, pointer(exc))
        # equivalent of passing C pointers &exc[1], &exc[2], &exc[3]:
        ccall((@pysym :PyErr_Fetch), Void, (UInt,UInt,UInt),
              pexc, pexc + sizeof(PyPtr), pexc + 2*sizeof(PyPtr))
        ccall((@pysym :PyErr_NormalizeException), Void, (UInt,UInt,UInt),
              pexc, pexc + sizeof(PyPtr), pexc + 2*sizeof(PyPtr))
        new(msg, exc[1], exc[2], exc[3])
    end
end

function show(io::IO, e::PyError)
    print(io, "PyError",
          isempty(e.msg) ? e.msg : string(" (",e.msg,")"),
          " ")

    if e.T.o == C_NULL
        println(io, "None")
    else
        println(io, pystring(e.T), "\n", pystring(e.val))
    end

    if e.traceback.o != C_NULL
        o = pycall(format_traceback, PyObject, e.traceback)
        if o.o != C_NULL
            for s in PyVector{AbstractString}(o)
                print(io, s)
            end
        end
    end
end

#########################################################################
# Conversion of Python exceptions into Julia exceptions

# call to discard Python exceptions
pyerr_clear() = ccall((@pysym :PyErr_Clear), Void, ())

function pyerr_check(msg::AbstractString, val::Any)
    if ccall((@pysym :PyErr_Occurred), PyPtr, ()) != C_NULL
        throw(PyError(msg))
    end
    val # the val argument is there just to pass through to the return value
end

pyerr_check(msg::AbstractString) = pyerr_check(msg, nothing)
pyerr_check() = pyerr_check("")

# extract function name from ccall((@pysym :Foo)...) or ccall(:Foo,...) exprs
callsym(s::Symbol) = s
callsym(s::QuoteNode) = s.value
import Base.Meta.isexpr
callsym(ex::Expr) = isexpr(ex,:macrocall,2) ? callsym(ex.args[2]) : isexpr(ex,:ccall) ? callsym(ex.args[1]) : ex

# Macros for common pyerr_check("Foo", ccall((@pysym :Foo), ...)) pattern.
macro pycheck(ex)
    :(pyerr_check($(string(callsym(ex))), $(esc(ex))))
end

# Macros to check that ccall((@pysym :Foo), ...) returns value != bad
macro pycheckv(ex, bad)
    quote
        val = $(esc(ex))
        if val == $(esc(bad))
            # throw a PyError if available, otherwise throw ErrorException
            pyerr_check($(string(callsym(ex))))
            error($(string(callsym(ex))), " failed")
        end
        val
    end
end
macro pycheckn(ex)
    :(@pycheckv $(esc(ex)) C_NULL)
end
macro pycheckz(ex)
    :(@pycheckv $(esc(ex)) -1)
end

#########################################################################
# Mapping of Julia Exception types to Python exceptions

const pyexc = Dict{DataType, PyPtr}()
type PyIOError <: Exception end

function pyexc_initialize()
    pyexc[Exception] = @pyglobalobjptr :PyExc_RuntimeError
    pyexc[ErrorException] = @pyglobalobjptr :PyExc_RuntimeError
    pyexc[SystemError] = @pyglobalobjptr :PyExc_SystemError
    pyexc[TypeError] = @pyglobalobjptr :PyExc_TypeError
    pyexc[ParseError] = @pyglobalobjptr :PyExc_SyntaxError
    pyexc[ArgumentError] = @pyglobalobjptr :PyExc_ValueError
    pyexc[KeyError] = @pyglobalobjptr :PyExc_KeyError
    pyexc[LoadError] = @pyglobalobjptr :PyExc_ImportError
    pyexc[MethodError] = @pyglobalobjptr :PyExc_RuntimeError
    pyexc[EOFError] = @pyglobalobjptr :PyExc_EOFError
    pyexc[BoundsError] = @pyglobalobjptr :PyExc_IndexError
    pyexc[DivideError] = @pyglobalobjptr :PyExc_ZeroDivisionError
    pyexc[DomainError] = @pyglobalobjptr :PyExc_RuntimeError
    pyexc[OverflowError] = @pyglobalobjptr :PyExc_OverflowError
    pyexc[InexactError] = @pyglobalobjptr :PyExc_ArithmeticError
    pyexc[OutOfMemoryError] = @pyglobalobjptr :PyExc_MemoryError
    pyexc[StackOverflowError] = @pyglobalobjptr :PyExc_MemoryError
    pyexc[UndefRefError] = @pyglobalobjptr :PyExc_RuntimeError
    pyexc[InterruptException] = @pyglobalobjptr :PyExc_KeyboardInterrupt
    pyexc[PyIOError] = @pyglobalobjptr :PyExc_IOError
end

function pyraise(e)
    eT = typeof(e)
    pyeT = haskey(pyexc::Dict, eT) ? pyexc[eT] : pyexc[Exception]
    ccall((@pysym :PyErr_SetString), Void, (PyPtr, Cstring),
          pyeT, string("Julia exception: ", e))
end

function pyraise(e::PyError)
    ccall((@pysym :PyErr_Restore), Void, (PyPtr, PyPtr, PyPtr),
          e.T, e.val, e.traceback)
    e.T.o = e.val.o = e.traceback.o = C_NULL # refs were stolen
end


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

mutable struct DiskProperties <: SwaggerModel
    timeCreated::Any # spec type: Union{ Nothing, DateTime } # spec name: timeCreated
    osType::Any # spec type: Union{ Nothing, String } # spec name: osType
    creationData::Any # spec type: Union{ Nothing, CreationData } # spec name: creationData
    diskSizeGB::Any # spec type: Union{ Nothing, Int32 } # spec name: diskSizeGB
    encryptionSettings::Any # spec type: Union{ Nothing, EncryptionSettings } # spec name: encryptionSettings
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState

    function DiskProperties(;timeCreated=nothing, osType=nothing, creationData=nothing, diskSizeGB=nothing, encryptionSettings=nothing, provisioningState=nothing)
        o = new()
        validate_property(DiskProperties, Symbol("timeCreated"), timeCreated)
        setfield!(o, Symbol("timeCreated"), timeCreated)
        validate_property(DiskProperties, Symbol("osType"), osType)
        setfield!(o, Symbol("osType"), osType)
        validate_property(DiskProperties, Symbol("creationData"), creationData)
        setfield!(o, Symbol("creationData"), creationData)
        validate_property(DiskProperties, Symbol("diskSizeGB"), diskSizeGB)
        setfield!(o, Symbol("diskSizeGB"), diskSizeGB)
        validate_property(DiskProperties, Symbol("encryptionSettings"), encryptionSettings)
        setfield!(o, Symbol("encryptionSettings"), encryptionSettings)
        validate_property(DiskProperties, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type DiskProperties

const _property_map_DiskProperties = Dict{Symbol,Symbol}(Symbol("timeCreated")=>Symbol("timeCreated"), Symbol("osType")=>Symbol("osType"), Symbol("creationData")=>Symbol("creationData"), Symbol("diskSizeGB")=>Symbol("diskSizeGB"), Symbol("encryptionSettings")=>Symbol("encryptionSettings"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_DiskProperties = Dict{Symbol,String}(Symbol("timeCreated")=>"DateTime", Symbol("osType")=>"String", Symbol("creationData")=>"CreationData", Symbol("diskSizeGB")=>"Int32", Symbol("encryptionSettings")=>"EncryptionSettings", Symbol("provisioningState")=>"String")
Base.propertynames(::Type{ DiskProperties }) = collect(keys(_property_map_DiskProperties))
Swagger.property_type(::Type{ DiskProperties }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_DiskProperties[name]))}
Swagger.field_name(::Type{ DiskProperties }, property_name::Symbol) =  _property_map_DiskProperties[property_name]

const _allowed_DiskProperties_osType = ["Windows", "Linux"]

function check_required(o::DiskProperties)
    (getproperty(o, Symbol("creationData")) === nothing) && (return false)
    true
end

function validate_property(::Type{ DiskProperties }, name::Symbol, val)
    if name === Symbol("osType")
        Swagger.validate_param(name, "DiskProperties", :enum, val, _allowed_DiskProperties_osType)
    end
end

ext: jl


text:
# Jonas August, 2/2022

#=
This problem was asked by Two Sigma.

Using a function rand7() that returns an integer from 1 to 7 (inclusive) with uniform probability, implement a function rand5() that returns an integer from 1 to 5 (inclusive).
=#

rand7() = rand(1:7)

function rand5()
    x = 6
    while x > 5  # Keep sampling uniformly on 1:7 until in range 1:5
        x = rand7()
    end
    x
end

sample = [rand5() for _ in 1:50000]

" Calculate 5 bin histogram. "
function hist5(sample::Vector{Int})
    h = zeros(Int, 5)
    for x in sample
        h[x] += 1
    end
    h
end

# Test uniformity of probability: All bins should have similar counts.
println(hist5(sample))

max_stars_repo_path: Random from random - Problem 71.jl

max_issues_repo_path: Random from random - Problem 71.jl

max_forks_repo_path: Random from random - Problem 71.jl


text:
### A Pluto.jl notebook ###
# v0.12.15

using Markdown
using InteractiveUtils

# ╔═╡ 695f3a2c-2640-11eb-3c75-93d1ed0e053a
html"""

"""

# ╔═╡ bf750d4c-25b9-11eb-17ca-dd85d876688d
md"""
# Statistics With Julia Example
"""

# ╔═╡ 30719ba6-25d4-11eb-38b9-11df9c03dfa8
md"""
## Example
### Bubble Sort

A basic sorting algorithm.

Takes an input array, index $1,...,n$, then sorts the elemnts smallesr to largest by allowing the larger elemnts, or "bubble".

```julia

function bubbleSort!(a)
    n = length(a)
    for i in 1:n-1
        for j in 1: n-i
            if a[j] > a[j+1]
                a[j], a[j+1] = a[j+1], a[j]
            end
        end
    end
    return a
end
```
"""

# ╔═╡ ac77aa56-2640-11eb-3a03-fdec7724034b
md"""
### Roots of a Polynomial

Consider the polynomial

$f(x) = a_nx^n + a_{n-1}x^{n-1} + ...+ a_1x + a_0$

with real-valued coefficients $a_0, ...a_n$.find all $x$ values that solve the equation $f(x)=0$.

```julia
using Roots
function polynomialGenerator(a...)
    n = length(a) - 1
    poly = function(x)
				return sum([a[i+1] * x^i for i in 0:n])
			end
	return poly
end

polynomial = polynomialGenerator(1,3,-10)
zeroVals = find_zeros(polynomial,-10,10)
println("Zeros of the function f(x): ", zeroVals)
```

"""

# ╔═╡ d15f3bc6-2641-11eb-2160-87492eb2595f
md"""
### Markov Chain

1. By raising the Matrix $P$ to a high power, the limiting distribution is obtained in any row.$\pi_i = lim_{n \to \infty}[P^n]_{i,j}$ for any index, $j$.

```julia
P = [0.5 0.4 0.1;
	 0.3 0.2 0.5;
	 0.5 0.3 0.2]

res1 = (P^100)[1, :]
```

2. We solve the linear systerm of equations,$\pi P = \pi$ and $\sum_{i=1}^{3} \pi_i = 1$

```julia
A = vcat((P' - I)[1:2, :], ones(3)')
b = [0 0 1]'
res2 = A\b
```
"""

# ╔═╡ 7b005a6e-2644-11eb-3f9a-0523247ad483
md"""
### Web,JSON and String Process

The **JSON** format uses '{}' characters to enclose a hierarchical nested structure of key value pairs.
```julia
using HTTP, JSON
data = HTTP.request("GET", "https://ocw.mit.edu/ans7870/6/6.006/s08/lecturenotes/files/t8.shakespeare.txt")
shakespeare = String(data.body)
shakespeareWords = split(shakespeare)

jsonWords = HTTP.request("GET", "https://raw.githubusercontent.com/"*
"h-Klok/StatsWithJuliaBook/master/1_chapter/jsonCode.json")
parsedJsonDict = JSON.parse(String(jsonWords.body))

keywords = Array{STring}(parsedJsonDict["words"])

numberToSHow = parsedJsonDict["numToShow"])

wordCount = Dict([(x, count(w -> lowrcase(w) == lowercase(x), shakespeareWords)) for x in keywords])
sortedWordCount = sort(collect(wordCount),by=last,rev=true)
sortedWordCount[1:numberToShow]

```

"""

# ╔═╡ 5f3a8406-2646-11eb-050b-47d517ceadf8
md"""
### Basic Plotting

you can visit the website[http://docs.juliaplots.org/](http://docs.juliaplots.org/), to get more information about plots

```julia
using Plots, LaTeXStrings, Measures
f(x,y) = x^2 + y^2
f0(x) = f(x,0)
f2(x) = f(x,2)
xVals, yVals = -5:0.1:5 , -5:0.1:5
plot(xVals, [f0.(xVals), f2.(xVals)],c=[:blue :red], xlims=(-5,5), legend=:top,
    ylims=(-5,25), ylabel=L"f(x,\cdot)", label=[L"f(x,0)" L"f(x,2)"])
p1 = annotate!(0, -0.2, text("(0,0) The minimum\n of f(x,0)", :left, :top, 10))
z=[f(x,y) for y in yVals, x in xVals]
p2 = surface(xVals, yVals, z, c=cgrad([:blue, :red]),legend=:none,
    ylabel="y", zlabel=L"f(x,y)")
M = z[1:10,1:10]
p3 = heatmap(M, c=cgrad([:blue, :red]), yflip=true, ylabel="y",
    xticks=([1:10;], xVals), yticks=([1:10;], yVals))
plot(p1, p2, p3, layout=(3,1), size=(1000,1200), xlabel="x", margin=5mm)
```

"""

# ╔═╡ Cell order:
# ╟─695f3a2c-2640-11eb-3c75-93d1ed0e053a
# ╟─bf750d4c-25b9-11eb-17ca-dd85d876688d
# ╟─30719ba6-25d4-11eb-38b9-11df9c03dfa8
# ╟─ac77aa56-2640-11eb-3a03-fdec7724034b
# ╟─d15f3bc6-2641-11eb-2160-87492eb2595f
# ╟─7b005a6e-2644-11eb-3f9a-0523247ad483
# ╟─5f3a8406-2646-11eb-050b-47d517ceadf8

max_stars_repo_name: ZhouZhuofei/juliaDL

max_issues_repo_name: ZhouZhuofei/juliaDL

max_forks_repo_name: ZhouZhuofei/juliaDL


text:
########################## non-stationary kernels ##############################
abstract type DotProductKernel{T} <: MercerKernel{T} end
(k::DotProductKernel)(x, y) = k(dot(x, y))

########################### dot product kernel #################################
struct Dot{T} <: DotProductKernel{T} end
Dot() = Dot{Union{}}()
@functor Dot
@inline (k::Dot)(d::Number) = d

# all of the kernels below can be defined as a composition of Dot(), with other kernel modifications
Line(σ::Real = 0.) = Dot() + σ
Polynomial(d::Int, σ::Real = 0.) = Line(σ)^d
const Poly = Polynomial

# exponential inner product kernel
# can be interpreted as infinite weighted combination of polynomial kernels
# WARNING: not well behaved for large inner products
struct ExponentialDot{T} <: DotProductKernel{T} end
ExponentialDot() = ExponentialDot{Union{}}()
@functor ExponentialDot
(k::ExponentialDot)(d::Number) = exp(d)

############################# Matrix kernel ####################################
# this could be discrete input kernel, as opposed to a matrix valued kernel
struct MatrixKernel{T, AT<:AbstractMatrix{T}} <: MercerKernel{T}
    A::AT
end
@functor MatrixKernel
(k::MatrixKernel)(x::Integer, y::Integer) = k.A[i,j]

############################ Brownian kernel ###################################
# has stationary increments
struct Brownian{T} <: MercerKernel{T} end
@functor Brownian
Brownian() = Brownian{Union{}}()
(k::Brownian)(x::Real, y::Real) = min(x, y)

######################### Finite Basis kernel ##################################
# finite basis function (linear regression) kernel,
struct FiniteBasis{T, B<:Union{AbstractVector, Tuple}} <: MercerKernel{T}
    basis::B # tuple of vector of functions
    function FiniteBasis{T}(basis) where {T}
        length(basis) ≥ 1 || throw("basis is empty: length(basis) = $(length(basis))")
        new{T, typeof(basis)}(basis)
    end
end
@functor FiniteBasis
FiniteBasis(basis) = FiniteBasis{Float64}(basis)
(k::FiniteBasis)(x, y) = sum(b->b(x)*b(y), k.basis)

# TODO: have in-place versions?
function basis(k::FiniteBasis, x::AbstractVector)
    U = zeros(eltype(k), (length(x), length(k.basis)))
    for (i, b) in enumerate(k.basis) # could use threads here
        @. U[:, i] = b(x)
    end
    U
end

function gramian(k::FiniteBasis, x::AbstractVector, y::AbstractVector)
    r = length(k.basis)
    if length(x) > r && length(y) > r # condition on low-rankness
        U = basis(k, x)
        V = x === y ? U : basis(k, y)
        return LazyMatrixProduct(U, V')
    else
        return Gramian(k, x, y) # should these be dense?
    end
end

########################### neural network kernel ##############################
struct NeuralNetwork{T} <: MercerKernel{T}
    σ::T # have to restrict to be positive
    # TODO: inner constructor
end
const NN = NeuralNetwork
@functor NN
# NN{T}(σ::T = one(T)) where {T} = NN(σ)
NN() = NN{Float64}(0.)

function (k::NN)(x, y)
    l = Line(k.σ)
    2/π * asin(l(x, y) / sqrt((1 + l(x, x)) * (1 + l(y, y))))
end
# NN(σ) = π/2 * (asin ∘ (Line(σ) / x->√(1+dot(x,x))))
# written with less sleek syntax: linear + normalized + asin (transform)
# NN(σ) = Transform(VerticalRescaling(Line(σ), x->1/√(1+dot(x,x))), asin)

########################## more future additions ###############################
# IDEA: Deep neural network kernel, convolutional NN kernel,
# Gibbs non-stationary construction etc.


text:
# * BitIntegers

module BitIntegers

import Base: &, *, +, -, <, <<, <=, ==, >>, >>>, |, ~, AbstractFloat, add_with_overflow,
             bswap, checked_abs, count_ones, div, flipsign, leading_zeros, mod,
             mul_with_overflow, ndigits0zpb, promote_rule, rem, sub_with_overflow,
             trailing_zeros, typemax, typemin, unsigned, xor

using Base: add_int, and_int, ashr_int, bswap_int, checked_sadd_int, checked_sdiv_int,
            checked_smul_int, checked_srem_int, checked_ssub_int, checked_uadd_int,
            checked_udiv_int, checked_umul_int, checked_urem_int, checked_usub_int, ctlz_int,
            ctpop_int, cttz_int, flipsign_int, lshr_int, mul_int, ndigits0z, ndigits0znb,
            neg_int, not_int, or_int, shl_int, sitofp, sle_int, slt_int, sub_int, uinttype,
            uitofp, ule_int, ult_int, xor_int

using Base.GMP: ispos, Limb

using Core: bitcast, check_top_bit, checked_trunc_sint, checked_trunc_uint, sext_int,
            trunc_int, zext_int

import Random: rand, Sampler
using Random: AbstractRNG, Repetition, SamplerType, LessThan, Masked

# * types definition & aliases

abstract type AbstractBitUnsigned <: Unsigned end
abstract type AbstractBitSigned   <: Signed   end

# ** @define_integers

macro define_integers(n::Int, SI=nothing, UI=nothing)
    if SI === UI === nothing
        SI = Symbol(:Int, n)
        UI = Symbol(:UInt, n)
    elseif SI === nothing || UI === nothing
        throw(ArgumentError(
            "@define_integers requires exactly 1 or 3 arguments, got 2"))
    end

    WSI = Symbol(:Int, 2n) # Wide
    WUI = Symbol(:UInt, 2n)

    WSI = isdefined(__module__, WSI) ? WSI : BigInt
    WUI = isdefined(__module__, WUI) ? WUI : BigInt

    quote
        # `esc` is necessary only on versions < 1.1
        primitive type $(esc(SI)) <: AbstractBitSigned   $n end
        primitive type $(esc(UI)) <: AbstractBitUnsigned $n end

        Base.Signed(x::$(esc(UI)))   = $(esc(SI))(x)
        Base.Unsigned(x::$(esc(SI))) = $(esc(UI))(x)
        Base.uinttype(::Type{$(esc(SI))}) = $(esc(UI))
        Base.uinttype(::Type{$(esc(UI))}) = $(esc(UI))

        Base.widen(::Type{$(esc(SI))}) = $(esc(WSI))
        Base.widen(::Type{$(esc(UI))}) = $(esc(WUI))
    end
end

# ** instanciation of types

const _DEFINED_SIZES = (256, 512, 1024)

# reverse so that widen knows about bigger types first:
for n = reverse(_DEFINED_SIZES)
    SI = Symbol(:Int, n)
    UI = Symbol(:UInt, n)
    sistr = Symbol(:int, n, :_str)
    uistr = Symbol(:uint, n, :_str)
    @eval begin
        @define_integers $n

        macro $sistr(s)
            return parse($SI, s)
        end
        macro $uistr(s)
            return parse($UI, s)
        end

        export $SI, $UI, $(Symbol("@", sistr)), $(Symbol("@", uistr))
    end
end

# ** type lists & aliases

const BitUnsigned_types = # (UInt256, ...)
    tuple((getfield(@__MODULE__, (Symbol(:UInt, n))) for n in _DEFINED_SIZES)...)

const BitSigned_types = # (Int256, ...)
    tuple((getfield(@__MODULE__, (Symbol(:Int, n))) for n in _DEFINED_SIZES)...)

const BitInteger_types = (BitSigned_types..., BitUnsigned_types...)

const BitUnsigned = Union{BitUnsigned_types...}
const BitSigned   = Union{BitSigned_types...}
const BitInteger  = Union{BitInteger_types...}

## convenient abbreviations, only for internal use

const BBS = Base.BitSigned
const BBU = Base.BitUnsigned
const BBI = Base.BitInteger

# eXtended
const XBS = AbstractBitSigned
const XBU = AbstractBitUnsigned
const XBI = Union{XBS,XBU}

# Union, Unified
const UBS = Union{BBS,XBS}
const UBU = Union{BBU,XBU}
const UBI = Union{BBI,XBI}

# ** typemin, typemax

typemin(::Type{T}) where {T<:XBU} = convert(T, 0)
typemax(::Type{T}) where {T<:XBU} = ~convert(T, 0)

typemin(::Type{T}) where {T<:XBS} = convert(T, 1) << (sizeof(T)*8-1)
typemax(::Type{T}) where {T<:XBS} = bitcast(T, typemax(uinttype(T)) >> 1)

# * conversions, promotions

# ** unsigned

unsigned(x::XBS) = reinterpret(typeof(convert(Unsigned, zero(x))), x)

# ** integers

# U -> X
(::Type{T})(x::Union{UBI,Bool}) where {T<:XBI} = convertto(T, x)::T

# X -> B
for T in Base.BitInteger_types
    @eval Core.$(Symbol(:to, T))(x::XBI) = convertto($T, x)
    @eval Base.$(Symbol(T))(x::XBI) = convertto($T, x)::$T
end

@generated function convertto(::Type{T}, x::Union{UBI,Bool}) where {T<:UBI}
    x === Bool && return :(and_int(zext_int(T, x), T(1)))
    if T <: Unsigned
        if x <: Signed
            if sizeof(x) < sizeof(T)
                :(sext_int(T, check_top_bit(x)))
            elseif sizeof(x) == sizeof(T)
                :(bitcast(T, check_top_bit(x)))
            else
                :(checked_trunc_uint(T, x))
            end
        else # x <: Unsigned
            if sizeof(x) < sizeof(T)
                :(zext_int(T, x))
            elseif sizeof(x) == sizeof(T)
                x === T ?
                    :x  :
                    :(reinterpret(T, x))
            else
                :(checked_trunc_uint(T, x))
            end
        end
    else # T <: Signed
        if x <: Signed
            if sizeof(x) < sizeof(T)
                :(sext_int(T, x))
            elseif sizeof(x) == sizeof(T)
                x === T ?
                    :x  :
                    :(reinterpret(T, x))
            else
                :(checked_trunc_sint(T, x))
            end
        else # x <: Unsigned
            if sizeof(x) < sizeof(T)
                :(zext_int(T, x))
            elseif sizeof(x) == sizeof(T)
                :(bitcast(T, check_top_bit(x)))
            else
                :(checked_trunc_sint(T, check_top_bit(x)))
            end
        end
    end
end

@generated function _rem(x::Union{UBI,Bool}, ::Type{to}) where {to<:UBI}
    from = x
    to === from && return :x # this replaces Base's method for BBI
    if to.size < from.size
        :(trunc_int(to, x))
    elseif from === Bool
        :(convert(to, x))
    elseif from.size < to.size
        if from <: Signed
            :(sext_int(to, x))
        else
            :(convert(to, x))
        end
    else
        :(bitcast(to, x))
    end
end

rem(x::Union{UBI,Bool}, ::Type{to}) where {to<:XBI} = _rem(x, to)
rem(x::XBI, ::Type{to}) where {to<:UBI} = _rem(x, to)
# to disambiguate
rem(x::XBI, ::Type{to}) where {to<:XBI} = _rem(x, to)

@generated function promote_rule(::Type{X}, ::Type{Y}) where {X<:XBI,Y<:UBI}
    if X.size > Y.size
        X
    elseif X.size == Y.size
        X <: Unsigned ?
            X :
            Y
    else
        Y
    end
end

# ** BigInt

function rem(x::BigInt, ::Type{T}) where T<:XBI
    if sizeof(T) <= sizeof(Limb)
        iszero(x) ? zero(T) : flipsign(unsafe_load(x.d) % T, x.size)
    else
        u = zero(T)
        for l = 1:min(abs(x.size), cld(sizeof(T), sizeof(Limb)))
            u += (unsafe_load(x.d, l) % T) << ((sizeof(Limb)<<3)*(l-1))
        end
        flipsign(u, x.size)
    end
end

function (::Type{T})(x::BigInt) where T<:XBU
    if sizeof(T) < sizeof(Limb)
        convert(T, convert(Limb,x))
    else
        0 <= x.size <= cld(sizeof(T),sizeof(Limb)) || throw(InexactError(Symbol(string(T)), T, x))
        x % T
    end
end

function (::Type{T})(x::BigInt) where T<:XBS
    n = abs(x.size)
    if sizeof(T) < sizeof(Limb)
        SLimb = typeof(Signed(one(Limb)))
        convert(T, convert(SLimb, x))
    else
        0 <= n <= cld(sizeof(T),sizeof(Limb)) || throw(InexactError(Symbol(string(T)), T, x))
        y = x % T
        ispos(x) ⊻ (y > 0) && throw(InexactError(Symbol(string(T)), T, x)) # catch overflow
        y
    end
end

# ** floats

AbstractFloat(x::XBI) = Float64(x)

for T in (Float32, Float64)
    @eval begin
        (::Type{$T})(x::XBS) = sizeof(x) > 16 ? $T(big(x)) : sitofp($T, x)
        (::Type{$T})(x::XBU) = sizeof(x) > 16 ? $T(big(x)) : uitofp($T, x)
        promote_rule(::Type{$T}, ::Type{<:XBI}) = $T
    end
end

promote_rule(::Type{Float16}, ::Type{<:XBI}) = Float16

# TODO: avoid conversion to BigInt
(::Type{T})(x::AbstractFloat) where {T<:XBI} = T(BigInt(x))::T
# to disambiguate:
(::Type{T})(x::Float16) where {T<:XBI} = T(BigInt(x))::T

# * comparisons

(<)(x::T, y::T) where {T<:XBU} = ult_int(x, y)
(<)(x::T, y::T) where {T<:XBS} = slt_int(x, y)

(<=)(x::T, y::T) where {T<:XBU} = ule_int(x, y)
(<=)(x::T, y::T) where {T<:XBS} = sle_int(x, y)

==(x::UBS, y::UBU) = (x >= 0) & (unsigned(x) == y)
==(x::UBU, y::UBS) = (y >= 0) & (x == unsigned(y))
<( x::UBS, y::UBU) = (x <  0) | (unsigned(x) <  y)
<( x::UBU, y::UBS) = (y >= 0) & (x <  unsigned(y))
<=(x::UBS, y::UBU) = (x <  0) | (unsigned(x) <= y)
<=(x::UBU, y::UBS) = (y >= 0) & (x <= unsigned(y))

# * bit operations

(~)(x::XBI) = not_int(x)
(&)(x::T, y::T) where {T<:XBI} = and_int(x, y)
(|)(x::T, y::T) where {T<:XBI} = or_int(x, y)
xor(x::T, y::T) where {T<:XBI} = xor_int(x, y)

>>( x::UBS, y::UBU) = ashr_int(x, y)
>>( x::UBU, y::UBU) = lshr_int(x, y)
>>>(x::UBI, y::UBU) = lshr_int(x, y)
<<( x::UBI, y::UBU) = shl_int(x, y)

>>( x::UBI, y::Int) = ifelse(0 <= y, x >> unsigned(y),  x << unsigned(-y))
<<( x::UBI, y::Int) = ifelse(0 <= y, x << unsigned(y),  x >> unsigned(-y))
>>>(x::UBI, y::Int) = ifelse(0 <= y, x >>> unsigned(y), x << unsigned(-y))

count_ones(    x::XBI) = Int(ctpop_int(x))
leading_zeros( x::XBI) = Int(ctlz_int(x))
trailing_zeros(x::XBI) = Int(cttz_int(x))

function bswap(x::XBI)
    if sizeof(x) % 2 != 0
        # llvm instruction is invalid
        error("unimplemented")
    else
        bswap_int(x)
    end
end

flipsign(x::T, y::T) where {T<:XBS} = flipsign_int(x, y)

# this doesn't catch flipsign(x::BBS, y::BBS), which is more specific in Base
flipsign(x::UBS, y::UBS) = flipsign_int(promote(x, y)...) % typeof(x)

# * arithmetic operations

(-)(x::XBI)                    = neg_int(x)
(-)(x::T, y::T) where {T<:XBI} = sub_int(x, y)
(+)(x::T, y::T) where {T<:XBI} = add_int(x, y)
(*)(x::T, y::T) where {T<:XBI} = mul_int(x, y)

div(x::XBS, y::Unsigned) = flipsign(signed(div(unsigned(abs(x)), y)), x)
div(x::Unsigned, y::XBS) = unsigned(flipsign(signed(div(x, unsigned(abs(y)))), y))

rem(x::XBS, y::Unsigned) = flipsign(signed(rem(unsigned(abs(x)), y)), x)
rem(x::Unsigned, y::XBS) = rem(x, unsigned(abs(y)))

mod(x::XBS, y::Unsigned) = rem(y + unsigned(rem(x, y)), y)

# these operations fail LLVM for bigger types than UInt128
div(x::T, y::T) where {T<:XBS} = sizeof(T) > 16 ? T(div(big(x), big(y))) : checked_sdiv_int(x, y)
rem(x::T, y::T) where {T<:XBS} = sizeof(T) > 16 ? T(rem(big(x), big(y))) : checked_srem_int(x, y)
div(x::T, y::T) where {T<:XBU} = sizeof(T) > 16 ? T(div(big(x), big(y))) : checked_udiv_int(x, y)
rem(x::T, y::T) where {T<:XBU} = sizeof(T) > 16 ? T(rem(big(x), big(y))) : checked_urem_int(x, y)

# ** checked operations

add_with_overflow(x::T, y::T) where {T<:XBS} = checked_sadd_int(x, y)
add_with_overflow(x::T, y::T) where {T<:XBU} = checked_uadd_int(x, y)

sub_with_overflow(x::T, y::T) where {T<:XBS} = checked_ssub_int(x, y)
sub_with_overflow(x::T, y::T) where {T<:XBU} = checked_usub_int(x, y)

mul_with_overflow(x::T, y::T) where {T<:XBS} = sizeof(T) >= 16 ? broken_mul_with_overflow(x, y) : checked_smul_int(x, y)
mul_with_overflow(x::T, y::T) where {T<:XBU} = sizeof(T) >= 16 ? broken_mul_with_overflow(x, y) : checked_umul_int(x, y)

# cf. base/checked.jl
# TODO: check whether the specific implementation for [U]Int128 is better suited here
function broken_mul_with_overflow(x::T, y::T) where T<:XBS
    r = widemul(x, y)
    f = r % T != r
    r % T, f
end

#= broken
function broken_mul_with_overflow(x::T, y::T) where T<:XBU
    r = widemul(x, y)
    f = r % T != r
    r % T, f
end
=#

function checked_abs(x::XBS)
    r = ifelse(x<0, -x, x)
    r<0 && throw(OverflowError(string("checked arithmetic: cannot compute |x| for x = ", x, "::", typeof(x))))
    r
end

checked_abs(x::XBU) = x

# * misc

function ndigits0zpb(x::XBU, b::Int)
    # precondition: b > 1
    x == 0 && return 0
    b < 0   && return ndigits0znb(signed(x), b)
    b == 2  && return sizeof(x)<<3 - leading_zeros(x)
    b == 8  && return (sizeof(x)<<3 - leading_zeros(x) + 2) ÷ 3
    b == 16 && return sizeof(x)<<1 - leading_zeros(x)>>2
    # b == 10 && return ndigits0z(x) # TODO: implement ndigits0z(x)

    d = 0
    while x > typemax(Int)
        x = div(x,b)
        d += 1
    end
    x = div(x,b)
    d += 1

    m = 1
    while m <= x
        m *= b
        d += 1
    end
    return d
end

ndigits0zpb(x::XBS, b::Integer) = ndigits0zpb(unsigned(abs(x)), Int(b))
ndigits0zpb(x::XBU, b::Integer) = ndigits0zpb(x, Int(b))

# * rand

# ** scalar

rand(rng::AbstractRNG, ::SamplerType{T}) where {T<:XBI} = rand(rng, NBits{T}(sizeof(T) << 3))

# sampler produce a T with at least n random bits (from low to high bits)
struct NBits{T<:XBI} <: Sampler{T}
    n::Int
end

function rand(rng::AbstractRNG, sp::NBits{T}) where {T<:XBI}
    n = (sp.n + 7) >> 3 # bytes
    if n <= 16 # 128 bits
        n > 8 ? rand(rng, UInt128) % T :
        n > 4 ? rand(rng, UInt64)  % T :
        n > 2 ? rand(rng, UInt32)  % T :
        n > 1 ? rand(rng, UInt16)  % T :
                rand(rng, UInt8)   % T
    else
        u = rand(rng, UInt128) % T
        while n > 16
            u = (u << 128) | rand(rng, UInt128)
            n -= 16
        end
        u
    end
end

# ** ranges

# we use SamplerRangeFast as divisions are not yet efficient (for SamplerRangeInt)

Sampler(::Type{<:AbstractRNG}, r::AbstractUnitRange{T}, ::Repetition) where {T<:XBI} =
    SamplerRangeFast(r)

# have to redefine SamplerRangeFast, as it's defined only for Base types in Random
struct SamplerRangeFast{U<:XBU,T<:XBI} <: Sampler{T}
    a::T      # first element of the range
    bw::UInt  # bit width
    m::U      # range length - 1
    mask::U   # mask generated values before threshold rejection
end

SamplerRangeFast(r::AbstractUnitRange{T}) where T<:BitInteger =
    SamplerRangeFast(r, uinttype(T))

function SamplerRangeFast(r::AbstractUnitRange{T}, ::Type{U}) where {T,U}
    isempty(r) && throw(ArgumentError("range must be non-empty"))
    m = (last(r) - first(r)) % uinttype(T) % U # % uinttype(T) to not propagate sign bit
    bw = (sizeof(U) << 3 - leading_zeros(m)) % UInt # bit-width
    mask = ((1 % U) << bw) - (1 % U)
    SamplerRangeFast{U,T}(first(r), bw, m, mask)
end

rand(rng::AbstractRNG, sp::SamplerRangeFast{<:XBI,T}) where {T} =
    rand(rng, LessThan(sp.m, Masked(sp.mask, NBits{T}(sp.bw)))) % T + sp.a

end # module


text:
# Predefined functions and values
const ϵ = 0.5
const ξ = 0.01

invfun(x) = 1/x

function rademacherDistribution!(v)
    o = one(eltype(v))
    v .= Base.rand.(Ref(-o:2*o:o))
end


text:
#=
Julia sample software pipeline usinn coroutine.
https://docs.julialang.org/en/v1.1/manual/parallel-computing/#Parallel-Computing-1
=#

const loader = Channel{Array{Float64,2}}(5)
const producer = Channel{Array{Float64,2}}(5)

function make_jobs(n)
   while true
       sleep(rand())
       put!(loader, rand(3, 3))
   end
end

function comm()
    println("start")
    while true
        data = take!(loader)
        sleep(rand())
        put!(producer, data)
    end
    println("end")
end

function use_pipeline(n)
    @async make_jobs(n)
    @async comm()
    println("process start")
    for i in 1:n
        @show take!(producer)
    end
end

function sequential(n)
    println("process start")
    for i in 1:n
        #do load something e.g. FileIO
        sleep(rand())
        # do some calculation
        sleep(rand())
        data = rand(3, 3)
        @show data
    end
end

function main()
    n=10
    @time sequential(n)
    @time use_pipeline(n)
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
@testset "Test issue #83" begin
    nxrec = 120
    nyrec = 100
    xrec = range(50f0, stop=1150f0, length=nxrec)
    yrec = range(100f0, stop=900f0, length=nyrec)
    zrec = 50f0

    # Construct 3D grid from basis vectors
    (x3d, y3d, z3d) = setup_3D_grid(xrec, yrec, zrec)

    for k=1:nyrec
        s = (k - 1) * nxrec + 1
        e = s + nxrec - 1
        @test all(x3d[k:nxrec:end] .== xrec[k])
        @test all(y3d[s:e] .== yrec[k])
    end
    @test all(z3d .== zrec)
end

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
T_u_inv(T) = [   0 1/T  0;
              -1/T   0  0;
                 0   0 -1]
T_u_inv_dot(T, Ṫ) = [    0 -Ṫ/T^2 0;
                     Ṫ/T^2      0 0;
                         0      0 0]
T_ω(T) = [0 -T  0;
          T  0  0;
          0  0  0]
skew(x) = [    0  -x[3]   x[2];
            x[3]      0  -x[1];
           -x[2]   x[1]     0]


text:
abstract type AbstractExecution end
abstract type Execution end
# Penalization API #
# ------------------------------------------------------------
init_subproblems!(ph::AbstractProgressiveHedgingSolver, subsolver::QPSolver) = init_subproblems!(ph, subsolver, ph.execution)
iterate!(ph::AbstractProgressiveHedgingSolver) = iterate!(ph, ph.execution)
init_workers!(ph::AbstractProgressiveHedgingSolver) = init_workers!(ph, ph.execution)
close_workers!(ph::AbstractProgressiveHedgingSolver) = close_workers!(ph, ph.execution)
resolve_subproblems!(ph::AbstractProgressiveHedgingSolver) = resolve_subproblems!(ph, ph.execution)
update_iterate!(ph::AbstractProgressiveHedgingSolver) = update_iterate!(ph, ph.execution)
update_subproblems!(ph::AbstractProgressiveHedgingSolver) = update_subproblems!(ph, ph.execution)
update_dual_gap!(ph::AbstractProgressiveHedgingSolver) = update_dual_gap!(ph, ph.execution)
calculate_objective_value(ph::AbstractProgressiveHedgingSolver) = calculate_objective_value(ph, ph.execution)
fill_submodels!(ph::AbstractProgressiveHedgingSolver, scenarioproblems) = fill_submodels!(ph, scenarioproblems, ph.execution)
# ------------------------------------------------------------
include("serial.jl")
include("distributed.jl")
include("channels.jl")
include("worker.jl")
include("synchronous.jl")
include("asynchronous.jl")


text:
module QuantumLattices

include("Interfaces.jl")
include("Prerequisites/Prerequisites.jl")
include("Essentials/Essentials.jl")

using .Interfaces
using .Prerequisites
using .Essentials

# Interfaces
export id, value, rank, dimension
export ⊕, ⊗, ⋅, add!, sub!, mul!, div!
export expand, expand!, decompose, decompose!, permute

# Essentials
export dtype, kind, update, update!, reset!

# Essentials.QuantumOperators
export QuantumOperator, OperatorUnit, ID, OperatorPack, OperatorProd, OperatorSum, Operator, Operators, LaTeX
export Transformation, Identity, Numericalization, MatrixRepresentation, Permutation, AbstractSubstitution, AbstractUnitSubstitution, UnitSubstitution, RankFilter
export ishermitian, idtype, optype, sequence, latexname, latexformat, matrix, matrix!

# Essentials.QuantumNumbers
export AbelianNumber, AbelianNumbers, @abeliannumber, SpinZ, ParticleNumber, SpinfulParticle, Momentum, Momentum₁, Momentum₂, Momentum₃
export periods, spinzs, particlenumbers, spinfulparticles

# Essentials.Spatials
export distance, azimuthd, azimuth, polard, polar, volume, isparallel, isonline, isintratriangle, issubordinate, reciprocals, translate, rotate, tile
export Translations, AbstractPID, PID, CPID, AbstractBond, Point, Bond, Lattice, SuperLattice, Cylinder, Bonds
export Segment, ReciprocalSpace, BrillouinZone, ReciprocalZone, ReciprocalPath
export pidtype, rcoord, icoord, isintracell, bonds!, bonds, @line_str, @rectangle_str, @hexagon_str, @translations_str
export allbonds, zerothbonds, insidebonds, acrossbonds, intrabonds, interbonds

# Essentials.DegreesOfFreedom
export IID, SimpleIID, CompositeIID, Internal, SimpleInternal, CompositeInternal, AbstractOID, Index, CompositeOID, OID, Hilbert
export IIDSpace, Subscript, Subscripts, AbstractCoupling, Coupling, Couplings, Metric, OIDToTuple, Table, Term
export statistics, iidtype, indextype, ismodulatable, abbr, @subscript_str, @couplings

# Essentials.QuantumSystems
## Canonical fermionic/bosonic systems
export majorana, annihilation, creation, flatex, blatex
export FID, Fock, FockCoupling, Onsite, Hopping, Pairing, Hubbard, InterOrbitalInterSpin, InterOrbitalIntraSpin, SpinFlip, PairHopping, Coulomb, FockTerm
export isnormalordered, @σ⁰_str, @σˣ_str, @σʸ_str, @σᶻ_str, @σ⁺_str, @σ⁻_str, @fc_str

## SU(2) spin systems
export slatex
export SID, Spin, SpinCoupling, SpinTerm, totalspin
export @heisenberg_str, @ising_str, @gamma_str, @dm_str, @sˣ_str, @sʸ_str, @sᶻ_str, @sc_str

## Phononic systems
export nlatex
export NID, Phonon, PhononCoupling, PhononKinetic, PhononPotential, PhononTerm

# Essentials.Frameworks
export Parameters, Boundary, Engine, AbstractGenerator, Formulation, Entry, CompositeGenerator, Generator, Image, Action, Assignment, Algorithm
export prepare!, run!, rundependences!, plain

end


text:
function LShape(domain::LShape, parent::ArbField)
    LShape(parent)
end

function Base.show(io::IO, domain::LShape)
    print(io, "L-shaped domain")
end

function area(domain::LShape)
    domain.parent(3)
end

function boundary_parameterization(t, domain::LShape, i::Integer)
    1 <= i <= 4 || throw(ErrorException("give 1 ≤ i ≤ 4 not i = $i"))

    if i == 1
        v = SVector(1, 0)
        w = SVector(1, 1)
    elseif i == 2
        v = SVector(1, 1)
        w = SVector(-1, 1)
    elseif i == 3
        v = SVector(-1, 1)
        w = SVector(-1, -1)
    elseif i == 4
        v = SVector(-1, -1)
        w = SVector(0, -1)
    end

    xy = v .+ t .* (w .- v)

    θ = atan(xy[2], xy[1])
    if typeof(t) == arb_series
        θ = ifelse(θ[0] < 0, θ + 2domain.parent(π), θ)
    else
        θ = ifelse(θ < 0, θ + 2domain.parent(π), θ)
    end
    (sqrt(xy[1]^2 + xy[2]^2), θ)
end

function boundary_points(
    domain::LShape,
    u::AbstractEigenfunction,
    ::Integer,
    n::Integer;
    distribution = :linear,
)
    points = Vector{Tuple{arb,arb}}(undef, n)

    for i = 1:n
        θ = 1.5 * domain.parent(π) * i / (n + 1)
        r = 1 / max(abs(sin(θ)), abs(cos(θ)))

        points[i] = (r, θ)
    end

    return points, fill(nothing, n)
end

function interior_points(domain::LShape, n::Integer; rng = MersenneTwister(42))
    # The points are generated by randomly picking a point in the
    # square [-1, 1]×[-1, 1] and keeping it if it is inside the
    # domain.
    points = Vector{Tuple{arb,arb}}(undef, n)

    for i = 1:n
        x = 2rand(rng) - 1
        y = 2rand(rng) - 1

        while x >= 0 && y <= 0
            x = 2rand(rng) - 1
            y = 2rand(rng) - 1
        end

        x, y = domain.parent(x), domain.parent(y)
        θ = atan(y, x)
        θ = ifelse(θ < 0, θ + 2domain.parent(π), θ)
        points[i] = (hypot(x, y), θ)
    end

    points
end


text:
using Distances

function get_3dtorus_distance_from_cols(X)

	dismat 		=	pairwise(Euclidean(), X)

	for i = [-1, 0, 1]
		for j = [-1, 0, 1]
			for k = [-1, 0, 1]
				Y 				=	X .+ [i,j,k]
				dismat_new		=	pairwise(Euclidean(), Y)
				dismat 			=	min.(dismat, dismat_new)
			end
		end
	end

	return dismat

end


text:
###############################################################################
#
# Function to deploy artifact
#
###############################################################################
"""
`PkgUtility` provides a convenient wrapper to deploy the artifact from a given list of files, for exmaple, `deploy_artifact!` function is used to deploy
    the [`GirddingMachine`](https://github.com/CliMA/GriddingMachine.jl) artifacts.

$(METHODLIST)

"""
function deploy_artifact! end

"""
What `deploy_artifact!` function does are
- determine if the artifact already exists in the `art_toml` file
- if `true`, skip the deployment
- if `false`
  - copy the file(s) to `~/.julia/artifacts/ARTIFACT_SHA/`
  - compress the artifact file(s) to a `.tar.gz` file
  - calculate the hash value of the compressed `tar.gz` file
  - bind the artifact file to the `.toml` file

Method for this deployment is

    deploy_artifact!(
                art_toml::String,
                art_name::String,
                art_locf::String,
                art_file::Vector{String},
                art_tarf::String,
                art_urls::Vector{String};
                new_file::Vector{String} = art_file)

Deploy the artifact, given
- `art_toml` Artifact `.toml` file location
- `art_name` Artifact name identitfier
- `art_locf` Local folder that stores the source files
- `art_file` Vector of the source file names
- `art_tarf` Folder location to store the compressed `.tar.gz` file
- `art_urls` Vector of public urls, where the compressed files are to be uploaded (user need to upload the file manually)
- `new_file` Optional. New file names of the copied files (same as `art_file` by default)

---
Examples
```julia
# deploy art_1.txt and art_2.txt as test_art artifact
deploy_artifact!("Artifacts.toml", "test_art", "./", ["art_1.txt", "art_2.txt], "./", ["https://public.server.url"]);

# deploy art_1.txt and art_2.txt as test_art artifact with new names
deploy_artifact!("Artifacts.toml", "test_art", "./", ["art_1.txt", "art_2.txt], "./", ["https://public.server.url"]; new_files=["new_1.txt", "new_2.txt"]);
```
"""
deploy_artifact!(
            art_toml::String,
            art_name::String,
            art_locf::String,
            art_file::Vector{String},
            art_tarf::String,
            art_urls::Vector{String};
            new_file::Vector{String} = art_file) =
(
    # querry whether the artifact exists
    _art_hash = artifact_hash(art_name, art_toml);

    # if artifact exists already skip
    if !isnothing(_art_hash) && artifact_exists(_art_hash)
        @info tinfo("Artifact $(art_name) already exists, skip it");
        return nothing;
    end;

    # create artifact
    @info tinfo("Artifact $(art_name) not found, deploy it now...");
    @info tinfo("Copying files into artifact folder...");
    _art_hash = create_artifact() do artifact_dir
        for i in eachindex(art_file)
            _in   = art_file[i];
            _out  = new_file[i];
            _path = joinpath(art_locf, _in);
            @info tinfo("Copying file $(_in)...");
            cp(_path, joinpath(artifact_dir, _out));
        end;
    end;

    # compress artifact
    @info tinfo("Compressing artifact $(art_name)...");
    _tar_loc  = "$(art_tarf)/$(art_name).tar.gz";
    _tar_hash = archive_artifact(_art_hash, _tar_loc);

    # bind artifact to download information
    _download_info = [("$(_url)/$(art_name).tar.gz", _tar_hash) for _url in art_urls];
    bind_artifact!(art_toml, art_name, _art_hash; download_info=_download_info, lazy=true, force=true);

    return nothing;
)

"""
In many cases, one might want to copy all the files in a folder to the target artifact, and iterate the file names is not convenient at all. Thus, a
    readily usable method is provided for this purpose:

    deploy_artifact!(
                art_toml::String,
                art_name::String,
                art_locf::String,
                art_tarf::String,
                art_urls::Vector{String})

Deploy the artifact, given
- `art_toml` Artifact `.toml` file location
- `art_name` Artifact name identitfier
- `art_locf` Local folder that stores the source files (all files will be copied into the artifact)
- `art_tarf` Folder location to store the compressed `.tar.gz` file
- `art_urls` Vector of public urls, where the compressed files are to be uploaded (user need to upload the file manually)

---
Examples
```julia
# deploy all files in target folder
deploy_artifact!("Artifacts.toml", "test_art", "./folder", "./", ["https://public.server.url"]);
```
"""
deploy_artifact!(
            art_toml::String,
            art_name::String,
            art_locf::String,
            art_tarf::String,
            art_urls::Vector{String}) =
(
    # querry all the files in the folder
    _art_files = String[_file for _file in readdir(art_locf)];

    # deploy the artifact
    deploy_artifact!(art_toml, art_name, art_locf, _art_files, art_tarf, art_urls);

    return nothing;
)


text:
"""
plotmap(dataset; args...)

dataset is a Cube from YAXArrays

variable :: the variable to plot

x        :: horizontal axis values, i.e. longitude

y        :: vertical axis values, i.e. latitude

reversey  :: if the data is reverse in the y direction

cmap     :: colormap

rmvalue  :: remove an specific value

utime    :: your time selection, format: Date("xxxx-xx-xx")

Returns :

fig, ax, pltobj = plot(...)

"""
function plotmap(dataset::YAXArray; x = :lon, y = :lat, variable = nothing, legcbar=" ",
    cmap = :plasma, figsize = (1240,650), utime = nothing, rmvalue =nothing, reversey = true)

    #steplonticks = nothing, steplatticks = nothing

    lon = collect(getAxis(string(x), dataset).values)
    lat = collect(getAxis(string(y), dataset).values)

    timestamp = getAxis("time", dataset)
    # drop menu with all variable symbols and timestamps if available. [Use Nodes and Menu functions]
    if timestamp == nothing
        data = variable == nothing ? dataset.data[:,:] : dataset[var = string(variable)].data[:,:]
        label = variable == nothing ? legcbar : string(variable)
    else
        timerange = getAxis("time", dataset).values
        if variable != nothing
            if utime != nothing
                data = dataset[var = string(variable), time = utime].data[:,:]
                label = string(variable) * "  $(utime)"
            else
                data = dataset[var = string(variable), time = Date(timerange[1])].data[:,:] # add keyword for user selection of time
                label = string(variable) * "  $(Date(timerange[1]))"
            end
        else
            if utime != nothing
                data = dataset[time = utime].data[:,:]
                label = legcbar * "  $(utime)"
            else
                data = dataset[time = Date(timerange[1])].data[:,:] # add keyword for user selection of time
                label = legcbar * "  $(Date(timerange[1]))"
            end
        end
    end

    minlon  = round(minimum(lon), RoundDown)
    maxlon  = round(maximum(lon), RoundDown)
    #steplon = (maxlon - minlon)/6
    #steplon = steplonticks == nothing ? steplon : steplonticks

    minlat  = round(minimum(lat), RoundDown)
    maxlat  = round(maximum(lat), RoundDown)
    #steplat = (maxlat - minlat)/6
    #steplat = steplatticks == nothing ? steplat : steplatticks

    data =  replace(data, missing=>NaN)
    if rmvalue != nothing
        data = rmvalue in data ? replace(data, rmvalue=>NaN) : data
    end
    data = reversey == true ? data[:,end:-1:1] : data

    fig = Figure(resolution = figsize, fontsize = 20)
    ax = Axis(fig, xlabel = string(x), ylabel = string(y),
        xtickalign = 1, ytickalign = 1, spinewidth = 0.75)

    pltobj = heatmap!(ax, lon, lat, data, colormap = cmap)
    cbar = Colorbar(fig, pltobj, label = label, width = 11, ticklabelsize =20,
        tickalign = 1, spinewidth =0.5)
    limits!(ax, minlon, maxlon, minlat, maxlat)
    hidespines!(ax, :t, :r)
    #ax.xticks = minlon:steplon:maxlon
    #ax.yticks = minlat:steplat:maxlat
    ax.xtickformat = "{:d}ᵒ"
    ax.ytickformat = "{:d}ᵒ"
    fig[1,1] = ax
    fig[1,2] = cbar
    fig, ax, pltobj #, cbar
end


text:
struct Coordinate{T, S}
    x::T
    y::S
end
function show(io::IO, c::Coordinate)
    print(io, "(")
    show(io, c.x)
    print(io, ", ")
    show(io, c.y)
    print(io, ")")
end

function compute_num_variables(hp::Hypothesis)
    total_free = 0
    total_dep = 0
    for (P, s) in hp.points
        total_free += numfree(s)
        total_dep += numdep(s)
    end
    return total_free, total_dep
end

numfree(::Free) = 2
numfree(::SemiFree) = 1
numfree(::Dependent) = 0
numdep(::Free) = 0
numdep(::SemiFree) = 1
numdep(::Dependent) = 2

function get_variables(hp::Hypothesis)
    nu, nx = compute_num_variables(hp)
    u, x = @polyvar u[1:nu-3] x[1:nx]
    return u, x
end

function assign_variables(hp::Hypothesis)
    u, x = get_variables(hp)
    point2var = OrderedDict{Symbol, Coordinate}()

    uidx = 1
    xidx = 1
    origin = true
    xaxis = true
    for (P, s) in hp.points
        if s isa Free
            if origin
                point2var[P] = Coordinate(0, 0)
                origin = false
            elseif xaxis
                point2var[P] = Coordinate(u[uidx], 0)
                uidx += 1
                xaxis = false
            else
                point2var[P] = Coordinate(u[uidx], u[uidx+1])
                uidx += 2
            end
        elseif s isa SemiFree
            point2var[P] = Coordinate(u[uidx], x[xidx])
            uidx += 1
            xidx += 1
        elseif s isa Dependent
            point2var[P] = Coordinate(x[xidx], x[xidx+1])
            xidx += 2
        end
    end
    for (P, s) in hp.points
        if s isa SemiFree
            eq = hp.constraints[s.eq1]
            pol = algebry(eq, point2var)[1]
            if maxdegree(pol, point2var[P].y) == 0
                point2var[P] = Coordinate(point2var[P].y, point2var[P].x)
            end
        end
    end
    return point2var, x, u
end


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

using Base.Test

# Tests for @__LINE__ inside and outside of macros
@test (@__LINE__) == 6

macro macro_caller_lineno()
    @test 9 == (@__LINE__) != __source__.line > 12
    return __source__.line
end

@test @macro_caller_lineno() == (@__LINE__) > 12

# @__LINE__ in a macro expands to the location of the macrocall in the source
# while __source__.line is the location of the macro caller
macro nested_LINE_expansion()
    return quote
        return (@emit_LINE, $(__source__.line))
    end
end
macro nested_LINE_expansion2()
    return :((@emit_LINE, $(__source__.line)))
end
macro emit_LINE()
    return quote
        (@__LINE__, $(__source__.line))
    end
end
@test (@emit_LINE) == ((@__LINE__) - 3, @__LINE__)
@test @nested_LINE_expansion() == ((@__LINE__() - 4, @__LINE__() - 12), @__LINE__())
@test @nested_LINE_expansion2() == ((@__LINE__() - 5, @__LINE__() - 9), @__LINE__())

include("test_sourcepath.jl")
thefname = "the fname!//\\&\1*"
include_string_test_func = include_string(@__MODULE__, "include_string_test() = @__FILE__", thefname)
@test include_string_test_func() == thefname
@test include_string(@__MODULE__, "Base.source_path()", thefname) == Base.source_path()
@test basename(@__FILE__) == "loading.jl"
@test isabspath(@__FILE__)

@test isdir(@__DIR__)
@test @__DIR__() == dirname(@__FILE__)
let exename = `$(Base.julia_cmd()) --precompiled=yes --startup-file=no`,
    wd = sprint(show, abspath(pwd(), "")),
    s_dir = sprint(show, joinpath(realpath(tempdir()), ""))
    @test wd != s_dir
    @test readchomp(`$exename -E "@__DIR__" -i`) == wd
    @test readchomp(`$exename -E "cd(()->eval(:(@__DIR__)), $s_dir)" -i`) == s_dir
    @test readchomp(`$exename -E "@__DIR__"`) == wd # non-interactive
end

# Issue #5789 and PR #13542:
mktempdir() do dir
    cd(dir) do
        let true_filename = "cAsEtEsT.jl", lowered_filename="casetest.jl"
            touch(true_filename)
            @test Base.isfile_casesensitive(true_filename)
            @test !Base.isfile_casesensitive(lowered_filename)

            # check that case-sensitivity only applies to basename of a path:
            if isfile(lowered_filename) # case-insensitive filesystem
                mkdir("cAsEtEsT")
                touch(joinpath("cAsEtEsT", true_filename))
                @test Base.isfile_casesensitive(joinpath("casetest", true_filename))
                @test !Base.isfile_casesensitive(joinpath("casetest", lowered_filename))
            end
        end

        # Test Unicode normalization; pertinent for OS X
        let nfc_name = "\U00F4.jl"
            touch(nfc_name)
            @test Base.isfile_casesensitive(nfc_name)
        end
    end
end

SAVED_LOAD_PATH = copy(LOAD_PATH)
empty!(LOAD_PATH)
dir = abspath(@__DIR__)
push!(LOAD_PATH, dir)

@test Base.find_in_path("test_sourcepath") == joinpath(dir, "test_sourcepath.jl")
@test Base.find_in_path(GenericString("test_sourcepath")) == joinpath(dir, "test_sourcepath.jl")
LOAD_PATH[end] = GenericString(LOAD_PATH[end])
@test Base.find_in_path("test_sourcepath") == joinpath(dir, "test_sourcepath.jl")

struct CustomLoader
    path::String
end
push!(LOAD_PATH, CustomLoader("abc"))
let name = randstring(20)
    @test_throws ArgumentError Base.find_in_path(name, nothing)
    Base.load_hook(prefix::CustomLoader, name::String, found) = joinpath(prefix.path, name)
    @test Base.find_in_path(name, nothing) == joinpath("abc", name)
end
@test Base.find_in_path("test_sourcepath", nothing) == joinpath("abc", "test_sourcepath")
Base.load_hook(prefix::CustomLoader, name::String, found::String) = found
@test Base.find_in_path("test_sourcepath", nothing) == joinpath(dir, "test_sourcepath.jl")

empty!(LOAD_PATH)
append!(LOAD_PATH, SAVED_LOAD_PATH)
@test LOAD_PATH == SAVED_LOAD_PATH

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
mutable struct IOState
	file::String
	f::IOStream
	tokens::Vector{String}
	function IOState(file)
		return new(file,open(file),String[])
	end
end

function tokenizer(line)
	tagopening = findall("<",line)
	tagclosing = findall(">",line)
	ntag = length(tagopening)
	@assert length(tagopening) == length(tagclosing)
	token = String[]
	if ntag <= 1
		push!(token,line)
		return token
	else
		push!(token,line[tagopening[1].start:tagclosing[1].start])
		append!(token,tokenizer(line[tagclosing[1].start+1:tagopening[end].start-1]))
		push!(token,line[tagopening[end].start:tagclosing[end].start])
	end
end

function comment_handler(line,state)
	comment_start = findall("<!--",line)
	comment_end = findall("-->",line)
	if length(comment_end) == length(comment_start) == 0
		if !isempty(strip(line))
			append!(state.tokens,tokenizer(line))
		end
		return nexttoken(state)
	elseif isempty(comment_end)
		line = line *" " * strip(readline(state.f))
		return comment_handler(line,state)
	elseif length(comment_end) != length(comment_start)
		error("Nested comments not allowed.")
	else
		tmp = line[1:comment_start[1].start-1]
		for i = 2:length(comment_start)
			tmp *= line[comment_end[i-1].stop+1:comment_start[i].start-1]
		end
		line = tmp * line[comment_end[end].stop+1:end]
		return comment_handler(line,state)
	end
end

function nexttoken(state)
	if !isempty(state.tokens)
		return popfirst!(state.tokens)
	elseif !eof(state.f)
		line = strip(readline(state.f))
		if isempty(line)
			return nexttoken(state)
		end
		if contains(line,"<!--")
			return comment_handler(line,state)
		end
		append!(state.tokens,tokenizer(line))
		return popfirst!(state.tokens)
	else
		return ""
	end
end

hastokens(state) = !isempty(state.tokens) || !eof(state.f)
istag(token) = !isempty(token) && token[1] == '<' && token[end] == '>' && token[end-1] != '/'
isemptytag(token) = !isempty(token) && token[1] == '<' && token[end] == '>' && token[end-1] == '/'

function readXMLTag(token,_empty=false)
	token = replace(token,"<"=>"")
	token = replace(token,">"=>"")
	if _empty
		if token[end] == '/'
			token = token[1:end-1]
		else
			error()
		end
	end
	token = replace(token, "="=>" ")
	pieces = split(token," ")
	filter!(!isempty,pieces)
	name = pieces[1]
	attributes = XMLAttribute[]
	@assert mod(length(pieces)-1,2)==0
	for i = 2:2:length(pieces)
		key = pieces[i]
		val = pieces[i+1]
		push!(attributes,XMLAttribute(key,val))
	end
	if _empty == false
		return XMLTag(name,attributes)
	else
		return XMLEmptyTag(name,attributes)
	end
end

function readXMLInclude(file)
	state = IOState(file)
	elements = Vector{XMLElement}()
	while hastokens(state)
		element = readXMLElement(state)
		if isdefined(element,:tag)
			push!(elements,element)
		end
	end
	close(state.f)
	return elements
end

function readXMLElement(state)
	element = XMLElement()
	while hastokens(state)
		token = nexttoken(state)
		#println(token)
		if istag(token)
			tag = readXMLTag(token)
			if !isdefined(element,:tag)
				element.tag = tag
			elseif tag.name[1]=='/' && element.tag.name == replace(tag.name,"/"=>"")
				return element
			else
				pushfirst!(state.tokens,token)
				if tag.name == "include" && tag.attributes[1].key == "file"
					file = replace(tag.attributes[1].val,"\""=>"")
					elements = readXMLInclude(joinpath(split(state.file,"/")[1:end-1]...,file)) #only linux
					append!(element.content,elements)
				else
					push!(element.content,readXMLElement(state))
				end
			end
		elseif isemptytag(token)
			@assert isdefined(element,:tag) "Empty-tag cannot be root"
			tag = readXMLTag(token,true)
			if tag.name == "include"
				@assert length(tag.attributes) == 1 "<include> does not support multiple attributes"
				if tag.attributes[1].key == "file"
					file = replace(tag.attributes[1].val,"\""=>"")
					elements = readXMLInclude(joinpath(split(state.file,"/")[1:end-1]...,file)) #only linux
					append!(element.content,elements)
				else
					error("Undefined attribute key $(tag.attributes)")
				end
			else
				push!(element.content,XMLElement(tag,Any[]))
			end
		else
			push!(element.content,token)
		end
	end
	return element
end

function writeAttribute(f::IOStream,attr)
	key,val = attr.key,attr.val
	write(f," $key=$val")
end

function writeTag(f::IOStream,tag::XMLTag,_sec=false)
	if !_sec
		write(f,"<$(tag.name)")
		for attr in tag.attributes
			writeAttribute(f,attr)
		end
		write(f,">\n")
	else
		write(f,"</$(tag.name)>\n")
	end
end

function writeTag(f::IOStream,tag::XMLEmptyTag,_sec=false)
	if !_sec
		write(f,"<$(tag.name)")
		for attr in tag.attributes
			writeAttribute(f,attr)
		end
		write(f,"/>\n")
	end
end

"""
`writeXMLElement(f::IOStream, el::XMLElement)`

Writes a `XMLElement` to an `IOStream`.
"""
function writeXMLElement(f::IOStream, el::XMLElement)
	writeTag(f,el.tag)
	for con in el.content
		if typeof(con) == XMLElement
			writeXMLElement(f,con)
		else
			write(f,con)
			write(f,"\n")
		end
	end
	el.tag.name[1] != '?' ? writeTag(f,el.tag,true) : nothing
end

"""
`Base.read(::Type{XMLElement}, file::String)`

Reads a XML file in location `file`.
Returns a `XMLElement`.
"""
function Base.read(::Type{XMLElement}, file::String)
	state = IOState(file)
	element = readXMLElement(state)
	close(state.f)
	return element
end

max_stars_repo_name: baxmittens/XMLParser

max_issues_repo_name: baxmittens/XMLParser

max_forks_repo_name: baxmittens/XMLParser


text:
using Revise,Zygote,LinearAlgebra,Random, ChainRules,BenchmarkTools
Random.seed!(1111)

#value of the function-used for testing purposes

#D(sum(A)) = [matrix of ones]
f_sA(A) = sum(A)
D_sA(A::TM, ::Type{T}=eltype(TM)) where {TM<:AbstractMatrix,T} = ones(T, size(A))

#D(At) = (noop- doesn't affect the input)
f_At(A) = A'
D_sAt(A::TM, ::Type{T}=eltype(TM)) where {TM,T} = D_A(A,T)

#D(AtA): let b=ones(K). Then D(AtA)=2A*bb'
f_AtA(A) = A'*A
function D_AtA(A::TM, ::Type{T}=eltype(TM)) where {TM,T}
  N,K = size(A)
  postmult = fill(T(2.0), K,K)
  #DA = D_sA(A)
  return A*postmult
end

#why doesn't this work?
function ChainRulesCore.rrule(::typeof(f_AtA), A::TM, ::Type{T}) where{TM,T}
  K::Int = size(A,2)
  f_AtA(A), (A)->A*fill(T(2.0), K,K)
end

function C_sAtA(A::TM, ::Type{T}=eltype(TM)) where {TM,T}
  N,K = size(A)
  AtA,pAtA  = rrule(*, A',A)
  sAtA,psAtA = rrule(sum, AtA)

  _, _psAtA = psAtA(T(1))
  _,_,_pAtA2 = pAtA(_psAtA)
  #_,_pAt = pAt(extern(_pAtA1))

  return T(2) .* _pAtA2
end

f_Aty(A,y) = A'*y
D_Aty(A::TM, y::AbstractVector{T}, ::Type{T}=eltype(A)) where {TM,T} =(
  reduce(hcat,[y for i in 1:size(A,2)]))

#D(AtAinv) = -AtAinv * D(AtA) * AtAinv
AtAinv(A) = cholesky(A'*A)\I
function D_AtAinv(A::TM, AtAinv::TM=cholesky(A'*A)\I, ::Type{T}=eltype(TM)) where {TM,T}
  N,K = size(A)
  return T(-1) .* A*AtAinv*fill(T(2.0),K,K)*AtAinv
  #return T(-1) .* DA*AtAinv*DAtA*AtAinv
end

#function rrule(::typeof(AtAinv), A::TM where TM, ::Type{T})
#  AtAinv(A),

f_AtAinvAty(A,y) = cholesky(X'*X)\y
#=D((X'X)^-1*X'y) wrt X:
  =D((X'X)^-1)*X'y+(X'X)^-1*D(X'y)
  =#

#=function D_AtAinvAty(A::TM, y::AbstractVector{T}, ::Type{T}=eltype(A)) where {TM,T}

  #build up least squares
  At, pAt = rrule(transpose, A)
  AtA,pAtA = rrule(*, At,A)
  cholAtA, pcholAtA = rrule(cholesky, AtA)
  Aty, pAty

  AtAinv = cholesky(A'*A)\I
  DAtAinv = D_AtAinv(A)
  Aty =A'*y
  DAty = D_Aty(A,y)

  term1 = DAtAinv*Aty
  term2 = AtAinv*DAty
  return  term1 .+ term2
end=#

function testD(f,D,args...; testname="", erroronfail::Bool = true)
  testval = D(args...)

  ∇(args) = gradient((args)->sum(f(args...)), args)
  verval = ∇(args)[1][1]

  local ε::Float64
  try
    ε =sum(abs.(testval .- verval))
  catch err
    msg = "Test $testname: failed probably due to dimensional mismatch: $err
      #args: $args
      testval: $testval
      verval: $verval"
    error(msg)
  end
  if args==verval
    @assert εrms==0
    @info "Test $testname: passed with a difference of 0"
  elseif testval≈verval
    @info "Test $testname: passed with absolute error of $(ε)"
  else
    msg = "Test $testname: failed with absolute error of $(ε)
      #args: $args
      testval: $testval
      verval: $verval"

    erroronfail && error(msg)
    @warn msg
  end

  return nothing
end

function testderivs(N=3,K=2)
  A= rand(N,K)
  y=rand(N)
  testD(f_sA, D_sA, A, testname="noop")
  testD(f_At, D_At, A, testname="transpose")
  testD(f_AtA, D_AtA, A, testname="AtA")
  testD(f_AtA, C_sAtA, A, testname="CAtA")
  #testD(f_AtA, C2_sAtA, A, testname="CAtA")

  testD(f_Aty, D_Aty, A, y, testname="Aty")

  g(A) = gradient((A)->sum(f_AtA(A)), A)
  @btime $g($A)[1][1]
  @btime C_sAtA($A)
  #@btime D_AtA($A)
  #testD(f_AtAinv, D_AtAinv, A, testname="AtAinv")
  #testD(f_AtAinvAty, D_AtAinvAty, A, y, testname="AtAinvAty")

end

testderivs(3,2)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
predl1_train_ens=Array{Float64}(undef,i_train,0)
predl1_test_ens=Array{Float64}(undef,i_test,0)
maetrs=Vector{Float64}(undef,0)
maetes=Vector{Float64}(undef,0)
for i=1:5
    global predl1_train_ens,predl1_test_ens,maetrs,maetes
    inds_all=1:i_train
    neurons=7000
    items_per_neuron=(Int64(floor(i_train/(neurons))))*ones(Int64,neurons);n_per_part=[0;cumsum(items_per_neuron)];n_per_part[end]=i_train
    @time include("~\\mnist_distort.jl")
    @time a_all,a_layer1,layer1=ANNBN.train_layer_1_sigmoid_fast(neurons,vars,i_train,n_per_part,inds_all,xx_train2,yy_train)
    predl1=[layer1 ones(i_train)]*a_layer1
    maetr=mean(abs.(yy_train-predl1))
    maximum(abs.(yy_train-predl1))
    predl1_test,layer1_test=ANNBN.predict_new(a_all,a_layer1,xx_test2,i_test,neurons);
    maete=mean(abs.(yy_test-predl1_test))
    maximum(abs.(yy_test-predl1_test))

    predl1_train_ens=[predl1_train_ens predl1]
    predl1_test_ens=[predl1_test_ens predl1_test]

    i1=predl1.<=0.5;predl1[i1].=tol1;i1=predl1.>0.5;predl1[i1].=1.0-tol1;i2=abs.(predl1.-yy_train).<1e-10;maetr=100sum(i2)/length(i2)
    i1=predl1_test.<=0.5;predl1_test[i1].=tol1;i1=predl1_test.>0.5;predl1_test[i1].=1.0-tol1;i2=abs.(predl1_test.-yy_test).<1e-10;
    maete=100sum(i2)/length(i2)

    println(i," ",maetr," ",maete)

    predl1=mean(predl1_train_ens,dims=2)
    maetr=mean(abs.(yy_train-predl1))
    maximum(abs.(yy_train-predl1))
    predl1_test=mean(predl1_test_ens,dims=2)
    maete=mean(abs.(yy_test-predl1_test))
    maximum(abs.(yy_test-predl1_test))
    i1=predl1.<=0.5;predl1[i1].=tol1;i1=predl1.>0.5;predl1[i1].=1.0-tol1;i2=abs.(predl1.-yy_train).<1e-10;maetr=100sum(i2)/length(i2)
    i1=predl1_test.<=0.5;predl1_test[i1].=tol1;i1=predl1_test.>0.5;predl1_test[i1].=1.0-tol1;i2=abs.(predl1_test.-yy_test).<1e-10;
    maete=100sum(i2)/length(i2)
    println(i," ",maetr," ",maete)
    push!(maetrs,maetr)
    push!(maetes,maete)
    # display(plot(maetrs))
    # display(plot!(maetes))

end

predl1=mean(predl1_train_ens,dims=2)
maetr=mean(abs.(yy_train-predl1))
maximum(abs.(yy_train-predl1))
predl1_test=mean(predl1_test_ens,dims=2)
maete=mean(abs.(yy_test-predl1_test))
maximum(abs.(yy_test-predl1_test))
i1=predl1.<=0.5;predl1[i1].=tol1;i1=predl1.>0.5;predl1[i1].=1.0-tol1;i2=abs.(predl1.-yy_train).<1e-10;maetr=100sum(i2)/length(i2)
i1=predl1_test.<=0.5;predl1_test[i1].=tol1;i1=predl1_test.>0.5;predl1_test[i1].=1.0-tol1;i2=abs.(predl1_test.-yy_test).<1e-10;
maete=100sum(i2)/length(i2)


text:
# set up

include("set_up.jl")

using JLD
using HDF5

################################################################################
##      parameters                                      					            ##
################################################################################

# set parameters for all jobs

# nbr parameters
set_nbr_params = 7

# nbr cores
nbr_of_cores = 4 # was 10

# length burn-in
burn_in = 1

# nbr iterations
nbr_iterations = 10000 #10000 # should be 20000

# length training data
length_training_data = 5000 # thid should ne 5000

# log-scale priors
log_scale_prior = false

# algorithm
mcmc_alg = "MCWM"  # set MCWM or PMCMC

# prob run MH update
beta_MH = 0.15 # should be 0.1

# load training data
load_tranining_data = true

# type of job
job = ARGS[1] # set work to simdata or new_data

# set jod dep. parameters
if job == "simdata"

	# jobname
	global_jobname = "est7"*job

	# nbr particels
	nbr_particels = 3*400

	# use simulated data
	sim_data = true # set to true to use sim data

	# data set
	data_set = "old" # was "old"

	# dt
	dt = 0.035 # new = 0.35 old = 0.035

	# dt_U
	dt_U = 1. # new = 1 old = 1

elseif job == "new_data"

	# jobname
	global_jobname = "est7"*job

	# nbr particels
	nbr_particels = 2000

	# use simulated data
	sim_data = false # set to true to use sim data

	# data set
	data_set = "new" # was "old"

	# dt
	dt = 0.35 # new = 0.35 old = 0.035

	# dt_U
	dt_U = 1. # new = 1 old = 1

end

################################################################################
##                         training data                                      ##
################################################################################

# set parameters
burn_in_tranining = 10000 # this should be 2000 when estimating 2 parameters

if set_nbr_params == 2
	burn_in_tranining = 1000 # this should be 2000
	length_training_data = 2000
end

nbr_iterations_tranining = burn_in_tranining +length_training_data
nbr_particels_tranining = 25 # should be 200
nbr_of_cores_tranining = 2 # should be 8

################################################################################
##      create (A)DA problem                                     					    ##
################################################################################

problem = set_up_gp_problem(nbr_of_unknown_parameters=set_nbr_params, use_sim_data = sim_data, data_set = data_set)
problem.alg_param.R = nbr_iterations
problem.alg_param.N = nbr_particels
problem.alg_param.burn_in = 1
problem.alg_param.nbr_of_cores = nbr_of_cores
problem.alg_param.length_training_data = length_training_data

problem.adaptive_update =  AMUpdate_gen(eye(set_nbr_params), 1/sqrt(set_nbr_params), 0.2, 1, 0.8, 25) # was 0.3

problem.alg_param.dt = dt
problem.alg_param.dt_U = dt_U
problem.alg_param.alg = mcmc_alg
problem.alg_param.compare_GP_and_PF = false
problem.alg_param.noisy_est = false
problem.alg_param.pred_method = "sample"
problem.alg_param.length_training_data = length_training_data
problem.alg_param.beta_MH = beta_MH

################################################################################
##                         set model parameters  for training                               ##
################################################################################

problem_training = set_up_problem(nbr_of_unknown_parameters=set_nbr_params, use_sim_data = sim_data, data_set = data_set)
problem_training.alg_param.alg = mcmc_alg
problem_training.alg_param.R = nbr_iterations_tranining
problem_training.alg_param.N = nbr_particels_tranining
problem_training.alg_param.burn_in = burn_in_tranining
problem_training.alg_param.nbr_of_cores = nbr_of_cores_tranining
problem_training.alg_param.dt = dt
problem_training.alg_param.dt_U = dt_U
problem_training.alg_param.alg = mcmc_alg
problem_training.adaptive_update =  AMUpdate_gen(eye(set_nbr_params), 1/sqrt(set_nbr_params), 0.15, 1, 0.8, 25)

################################################################################
##                generate training data                                     ###
################################################################################

if !load_tranining_data

	if !log_scale_prior
	  res_training, theta_training, loglik_training, cov_matrix = @time MCMC(problem_training, true, true)
	  #export_parameters(res_problem_normal_prior_est_AM_gen[2],jobname)
	else
	  res_training, theta_training, loglik_training, cov_matrix  = @time MCMC(problem_training_nonlog, true, true)
	end

	export_data(problem_training, res_training[1],"da_ada_gpMCMC_training_data"*global_jobname)

	save("gp_training_$(set_nbr_params)_par_test_new_code.jld", "res_training", res_training, "theta_training", theta_training, "loglik_training", loglik_training,"cov_matrix",cov_matrix)

else

	if job == "simdata"
		@load "gp_training_7_par_training_and_testsimdatalunarc_new.jld"
	elseif job == "new_data"
		@load "gp_training_7_par_training_and_testnew_datalunarc_new.jld"
	end

end

################################################################################
###            plot training data                                            ###
################################################################################

#export_data(problem_training, res_training[1],"da_ada_gpMCMC_training_data"*global_jobname)

#=
using PyPlot

for i = 1:set_nbr_params
  PyPlot.figure()
  PyPlot.plot(res_training[1].Theta_est[i,:])
end

for i = 1:set_nbr_params
  PyPlot.figure()
  PyPlot.plt[:hist](theta_training[i,:],100)
end

PyPlot.figure()
PyPlot.plt[:hist](loglik_training,100)
=#

println("Mean:")
println(mean(theta_training,2))
println("Std:")
println(std(theta_training,2))

################################################################################
###          Fit GP model                                                     ##
################################################################################

data_training = [theta_training; loglik_training']
data_test = [theta_test; loglik_test']

show(data_training[:,1:10])

if true

	# create gp object
	gp = GPModel("est_method",zeros(6), zeros(4),
	eye(length_training_data-20), zeros(length_training_data-20),zeros(2,length_training_data-20),
	collect(1:10))

	# fit gp model

	tic()

	# fit GP model
	if true #problem.alg_param.est_method == "ml"
	  # fit GP model using ml
	  #perc_outlier = 0.0 # used when using PMCMC for trainig data 0.05
	  #tail_rm = "left"

	  perc_outlier = 0.01 # should be 0.05
	  tail_rm = "left"
	  lasso = false # should be true

	  ml_est(gp, data_training,"SE", lasso,perc_outlier,tail_rm)
	else
	  error("The two stage estimation method is not in use")
	  #two_stage_est(gp, data_training)
	end

	time_fit_gp = toc()

end

################################################################################
##               Run DA-GP-MCMC                                              ###
################################################################################

accelerated_da = false

jobname = global_jobname*"da_gp_mcmc"

problem.model_param.theta_0 = mean(theta_training,2)

if !log_scale_prior

  res = dagpmcmc(problem_training, problem, gp, cov_matrix)

  mcmc_results = Result(res[1].Theta_est, res[1].loglik_est, res[1].accept_vec, res[1].prior_vec)

  # plot results
  export_data(problem, mcmc_results,jobname)
  #export_parameters(mcmc_results[2],jobname)
else
  # run adaptive PMCMC
  res, res_traning, theta_training, loglik_training, assumption_list, loglik_list = dagpMCMC(problem_training, problem, gp, cov_matrix)

  mcmc_results = Result(res[1][1].Theta_est, res[1][1].loglik_est, res[1][1].accept_vec, res[1][1].prior_vec)

  # plot results
  export_data(problem_nonlog, mcmc_results,jobname)
  #export_parameters(mcmc_results[2],jobname)
end

################################################################################
##                         set A-DA problem                               ##
################################################################################

accelerated_da = true

jobname = global_jobname*"ada_gp_mcmc_dt"

problem.model_param.theta_0 = mean(res_training[1].Theta_est[:,end-size(data_training,2)-size(data_test,2):end],2)[:]

################################################################################
##  Create features for classification models                                                                            ##
################################################################################

n = size(data_training,2)
n_burn_in = length(res_training[2]) -  length(loglik_training) - length(loglik_test) #problem_traning.alg_param.burn_in

idx_training_start = n_burn_in+1
idx_training_end = idx_training_start+n-1

idx_test_start = idx_training_end+1
idx_test_end = idx_test_start+n-1

idx_training_old = idx_training_start-1:idx_training_end-1
idx_test_old = idx_test_start-1:idx_test_end-1

loglik_training_old = res_training[1].loglik_est[idx_training_old]
loglik_test_old = res_training[1].loglik_est[idx_test_old]

dim = length(problem.model_param.theta_true)
data_signs = zeros(4,n)
data_signs[3,:] = data_training[dim+1,:]
data_signs[4,:] = loglik_training_old

std_pred_gp_star = zeros(n)

noisy_pred = problem.alg_param.noisy_est

for i = 1:n
  (loglik_est_star, var_pred_ml_star, prediction_sample_ml_star) = predict(data_training[1:dim,i],gp,noisy_pred)
  (loglik_est_old, var_pred_ml, prediction_sample_ml_old) = predict(res_training[1].Theta_est[:,idx_training_old[i]],gp,noisy_pred)
  data_signs[1,i] = prediction_sample_ml_star[1]
  data_signs[2,i] = prediction_sample_ml_old[1]
  std_pred_gp_star[i] = sqrt(var_pred_ml_star[1])
end

nbr_GP_star_geq_GP_old = zero(Int64)
nbr_case_1 = zero(Int64)
nbr_case_4 = zero(Int64)

targets_case_1_and_3 = []
data_case_1_and_3 = []

targets_case_2_and_4 = []
data_case_2_and_4 = []

for i = 1:n
  if data_signs[1,i] > data_signs[2,i]
    nbr_GP_star_geq_GP_old += 1
    data_case_1_and_3 = vcat(data_case_1_and_3, [data_training[1:dim,i]; data_signs[1,i]/data_signs[2,i]; std_pred_gp_star[i]])
    if data_signs[3,i] > data_signs[4,i]
      append!(targets_case_1_and_3, 1)
      nbr_case_1 += 1
    else
      append!(targets_case_1_and_3, 0)
    end
  elseif data_signs[3,i] > data_signs[4,i]
    data_case_2_and_4 = vcat(data_case_2_and_4, [data_training[1:dim,i]; data_signs[1,i]/data_signs[2,i]; std_pred_gp_star[i]])
    append!(targets_case_2_and_4, 0)
    nbr_case_4 += 1
  else
    data_case_2_and_4 = vcat(data_case_2_and_4, [data_training[1:dim,i]; data_signs[1,i]/data_signs[2,i]; std_pred_gp_star[i]])
    append!(targets_case_2_and_4, 1)
  end
end

# tansform features and set input data

# convert matricies to floats
data_case_1_and_3 = convert(Array{Float64,2},reshape(data_case_1_and_3, (dim+2, length(targets_case_1_and_3))))
data_case_2_and_4 = convert(Array{Float64,2},reshape(data_case_2_and_4, (dim+2, length(targets_case_2_and_4))))

targets_case_1_and_3 = convert(Array{Float64,1}, targets_case_1_and_3)
targets_case_2_and_4 = convert(Array{Float64,1}, targets_case_2_and_4)

################################################################################
##   set case model                                                          ###
################################################################################

select_case_model = ARGS[2] # logisticregression or dt

# fit model, i.e. est probabilities
nbr_GP_star_led_GP_old = n-nbr_GP_star_geq_GP_old

prob_case_1 = nbr_case_1/nbr_GP_star_geq_GP_old
prob_case_2 = (nbr_GP_star_led_GP_old-nbr_case_4)/nbr_GP_star_led_GP_old
prob_case_3 = 1-prob_case_1
prob_case_4 = nbr_case_4/nbr_GP_star_led_GP_old
prob_cases = [prob_case_1;prob_case_2;prob_case_3;prob_case_4]

println("Est prob:")
println(prob_cases)

if select_case_model == "biasedcoin"

  casemodel = BiasedCoin(prob_cases)

elseif select_case_model == "logisticregression"

  # tansformed data for logistic model
  mean_posterior = mean(theta_training,2)[:]

  input_data_case_1_and_3 = zeros(length(targets_case_1_and_3), dim+3)
  input_data_case_1_and_3[:,1] = sqrt((mean_posterior[1] - data_case_1_and_3[1,:]).^2)
  input_data_case_1_and_3[:,2] = sqrt((mean_posterior[2] - data_case_1_and_3[2,:]).^2)
  input_data_case_1_and_3[:,3] = sqrt((mean_posterior[3] - data_case_1_and_3[3,:]).^2)
  input_data_case_1_and_3[:,4] = sqrt((mean_posterior[4] - data_case_1_and_3[4,:]).^2)
  input_data_case_1_and_3[:,5] = sqrt((mean_posterior[5] - data_case_1_and_3[5,:]).^2)
  input_data_case_1_and_3[:,6] = sqrt((mean_posterior[6] - data_case_1_and_3[6,:]).^2)
  input_data_case_1_and_3[:,7] = sqrt((mean_posterior[7] - data_case_1_and_3[7,:]).^2)
  input_data_case_1_and_3[:,8] = sqrt(sum((repmat(mean_posterior', size(data_case_1_and_3,2))'-data_case_1_and_3[1:7,:]).^2,1))
  input_data_case_1_and_3[:,8] = data_case_1_and_3[8,:]
  input_data_case_1_and_3[:,end] = targets_case_1_and_3

  input_data_case_2_and_4 = zeros(length(targets_case_2_and_4), dim+3)
  input_data_case_2_and_4[:,1] = sqrt((mean_posterior[1] - data_case_2_and_4[1,:]).^2)
  input_data_case_2_and_4[:,2] = sqrt((mean_posterior[2] - data_case_2_and_4[2,:]).^2)
  input_data_case_2_and_4[:,3] = sqrt((mean_posterior[3] - data_case_2_and_4[3,:]).^2)
  input_data_case_2_and_4[:,4] = sqrt((mean_posterior[4] - data_case_2_and_4[4,:]).^2)
  input_data_case_2_and_4[:,5] = sqrt((mean_posterior[5] - data_case_2_and_4[5,:]).^2)
  input_data_case_2_and_4[:,6] = sqrt((mean_posterior[6] - data_case_2_and_4[6,:]).^2)
  input_data_case_2_and_4[:,7] = sqrt((mean_posterior[7] - data_case_2_and_4[7,:]).^2)
  input_data_case_2_and_4[:,8] = sqrt(sum((repmat(mean_posterior', size(data_case_2_and_4,2))'-data_case_2_and_4[1:7,:]).^2,1))
  input_data_case_1_and_3[:,9] = data_case_1_and_3[8,:]
  input_data_case_2_and_4[:,end] = targets_case_2_and_4

  input_data_case_1_and_3 = DataFrame(input_data_case_1_and_3)

  log_reg_model_case_1_and_3 = glm(@formula(x10 ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8  + x9), input_data_case_1_and_3, Binomial(), LogitLink())

  # fit model for cases 2 and 4

  input_data_case_2_and_4 = DataFrame(input_data_case_2_and_4)

  log_reg_model_case_2_and_4 = glm(@formula(x10 ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8  + x9), input_data_case_2_and_4, Binomial(), LogitLink())

  β_for_model1or3 = coef(log_reg_model_case_1_and_3)
  β_for_model2or4 = coef(log_reg_model_case_2_and_4)

  casemodel = LogisticRegression(β_for_model1or3, β_for_model2or4, mean_posterior)

elseif select_case_model == "dt"

  #standardization!(data_case_1_and_3)

  input_data_case_1_and_3 = zeros(length(targets_case_1_and_3), dim+2)

  input_data_case_1_and_3[:,1] =  data_case_1_and_3[1,:]
  input_data_case_1_and_3[:,2] =  data_case_1_and_3[2,:]
  input_data_case_1_and_3[:,3] =  data_case_1_and_3[3,:]
  input_data_case_1_and_3[:,4] =  data_case_1_and_3[4,:]
  input_data_case_1_and_3[:,5] =  data_case_1_and_3[5,:]
  input_data_case_1_and_3[:,6] =  data_case_1_and_3[6,:]
  input_data_case_1_and_3[:,7] =  data_case_1_and_3[7,:]
  input_data_case_1_and_3[:,8] =  data_case_1_and_3[8,:]
  input_data_case_1_and_3[:,end] = targets_case_1_and_3

  #standardization!(data_case_2_and_4)

  input_data_case_2_and_4 = zeros(length(targets_case_2_and_4), dim+2)

  input_data_case_2_and_4[:,1] =  data_case_2_and_4[1,:]
  input_data_case_2_and_4[:,2] =  data_case_2_and_4[2,:]
  input_data_case_2_and_4[:,3] =  data_case_2_and_4[3,:]
  input_data_case_2_and_4[:,4] =  data_case_2_and_4[4,:]'
  input_data_case_2_and_4[:,5] =  data_case_2_and_4[5,:]
  input_data_case_2_and_4[:,6] =  data_case_2_and_4[6,:]
  input_data_case_2_and_4[:,7] =  data_case_2_and_4[7,:]
  input_data_case_2_and_4[:,8] =  data_case_2_and_4[8,:]
  input_data_case_2_and_4[:,end] = targets_case_2_and_4

  # 4 features model
  features_1_and_3 = convert(Array, input_data_case_1_and_3[:, 1:dim+1])

  labels_1_and_3 = convert(Array, input_data_case_1_and_3[:, end])

  labels_1_and_3 = Array{String}(size(features_1_and_3,1))

  for i = 1:length(labels_1_and_3)
    if input_data_case_1_and_3[i,end] == 0
      labels_1_and_3[i] = "case 3"
    else
      labels_1_and_3[i] = "case 1"
    end
  end

  decisiontree1or3 = build_tree(labels_1_and_3, features_1_and_3)

  decisiontree1or3 = prune_tree(decisiontree1or3, 0.9)

  # tree based model for case 2 and 4

  # 4 features model
  features_case_2_and_4 = convert(Array, input_data_case_2_and_4[:, 1:dim+1])

  labels_case_2_and_4 = Array{String}(size(features_case_2_and_4,1))

  for i = 1:length(labels_case_2_and_4)
    if input_data_case_2_and_4[i,end] == 0
      labels_case_2_and_4[i] = "case 4"
    else
      labels_case_2_and_4[i] = "case 2"
    end
  end

  # train full-tree classifier
  decisiontree2or4 = build_tree(labels_case_2_and_4, features_case_2_and_4)

  decisiontree2or4 = prune_tree(decisiontree2or4, 0.9)

  casemodel = DT(decisiontree1or3, decisiontree2or4)

end

if !log_scale_prior
  # run adaptive PMCMC
  res = adagpmcmc(problem_training, problem, gp, casemodel, cov_matrix)

  mcmc_results = Result(res[1].Theta_est, res[1].loglik_est, res[1].accept_vec, res[1].prior_vec)

  # plot results
  export_data(problem, mcmc_results,jobname)
  #export_parameters(mcmc_results[2],jobname)
else
  # run adaptive PMCMC
  res, res_traning, theta_training, loglik_training, assumption_list, loglik_list = dagpMCMC(problem_training, problem, gp, cov_matrix)

  mcmc_results = Result(res[1][1].Theta_est, res[1][1].loglik_est, res[1][1].accept_vec, res[1][1].prior_vec)

  # plot results
  export_data(problem_nonlog, mcmc_results,jobname)
  #export_parameters(mcmc_results[2],jobname)
end

max_stars_repo_name: SamuelWiqvist/adamcmcpaper

max_issues_repo_name: SamuelWiqvist/adamcmcpaper

max_forks_repo_name: SamuelWiqvist/adamcmcpaper


text:
using Plots
gr(titlefontsize=12,size=(500,300),colorbar=false)

using FourierGPE

# Units: oscillator
L = (20.,20.,20.)
N = (64,64,64)
sim = Sim(L,N)
@unpack_Sim sim

import FourierGPE.V
V(x,y,z,t) = 0.5*(x^2 + y^2 + 4*z^2)

## set simulation parameters ====
γ = 0.5
μ = 15.0
tf = 0.3 # 1.5/γ
Nt = 150
t = LinRange(0.,tf,Nt)

## state functions
ψ0(x,y,z,μ,g) = sqrt(μ/g)*sqrt(max(1.0-V(x,y,z,0.0)/μ,0.0)+im*0.0)
healing(x,y,z,μ,g) = 1/sqrt(g*abs2(ψ0(x,y,z,μ,g)))

## initial state
x,y,z = X
y = y'; z = reshape(z,1,1,length(z))
ψi = ψ0.(x,y,z,μ,g)
ϕi = kspace(ψi,sim)

@pack_Sim! sim

plot(x,abs2.(ψ0.(x,0.,0.,μ,g)))
plot!(x,abs2.(ψi[:,32,32]))

## evolve in k space
sol = runsim(sim)

## ground state
ϕg = sol[end]
ψg = xspace(ϕg,sim)
showpsi(x,y',ψg[:,:,32])

R(w) = sqrt(2*μ/w^2)
R(1)
R(2)

plot(x,abs2.(ψ0.(x,0.,0.,μ,g)))
plot!(x,abs2.(ψg[:,32,32]))
μ/g

abs2.(ψg[32,32,32])

# animate a slice
anim = @animate for i=1:Nt
    ϕ = sol[i]
    zmid = findfirst(z .≈ 0)
    ψ = xspace(ϕ,sim)
    showpsi(x,y,ψ[:,:,zmid])
end
gif(anim,"./examples/3dquench.gif",fps=30)
#TODO debug
## animate isosurface in Makie
using Makie, AbstractPlotting

function dense(phi)
    ψm = xspace(phi,sim)
    density = abs2.(ψm)
    pmax = maximum(density)
    return density/pmax
end

function densityfilm(sol,Nt;file="3dquench.gif")
    saveto="examples/"*file
    scene = Scene()
    tindex = Node(1)
    scene = volume(lift(i -> dense(sol[i]), tindex),
    algorithm = :iso,
    color = (c3,0.25),
    show_axis=false,
    isovalue=3f0(.15))

    R = 70
    eyepos = Vec3f0(R,R,R)
    lookat = Vec3f0(18,18,0)

    record(scene, saveto, 1:Nt) do i
        update_cam!(scene, eyepos, lookat)
        rotate_cam!(scene, 0., -0.25, 0.)
        tindex[] = i
    end
    p = scene[end];
    return
end

p = densityfilm(sol,Nt,file="3dquench1.gif")

# p.attributes.attribute


text:
@testset "EigenAngles.jl" begin
    @info "Testing EigenAngles"

    @test isbits(EigenAngle(deg2rad(complex(85.0, -1.0))))
    @test EigenAngle(deg2rad(80-0.5im)) > EigenAngle(deg2rad(75-0.3im))

    @test_logs (:warn, "θ > 2π. Make sure θ is in radians.") EigenAngle(complex(85.0, 0.31))
end


text:
# Load D4M Module
include("loadD4M.jl")

# Run Assoc Intro Example
println("Running AI1_SetupTEST")
include("1Intro/1AssocIntro/AI1_SetupTEST.jl")
println("Running AI2_SubsrefTEST")
include("1Intro/1AssocIntro/AI2_SubsrefTEST.jl")
println("Running AI3_MathTEST")
include("1Intro/1AssocIntro/AI3_MathTEST.jl")
println("Running AI4_AdvConstructTEST")
include("1Intro/1AssocIntro/AI4_AdvConstructTEST.jl")

# Run Edge Art Example
println("Running EA1_GraphTEST")
include("1Intro/2EdgeArt/EA1_GraphTEST.jl")
println("Running EA2_SubsrefTEST")
include("1Intro/2EdgeArt/EA2_SubsrefTEST.jl")
println("Running EA3_SubGraphTEST")
include("1Intro/2EdgeArt/EA3_SubGraphTEST.jl")

# Run Entity Analysis Example
println("Running EA1_ReadTEST")
include("2Apps/1EntityAnalysis/EA1_ReadTEST.jl")
println("Running EA2_StatTEST")
include("2Apps/1EntityAnalysis/EA2_StatTEST.jl")
println("Running EA3_FacetTEST")
include("2Apps/1EntityAnalysis/EA3_FacetTEST.jl")
println("Running EA4_GraphTEST")
include("2Apps/1EntityAnalysis/EA4_GraphTEST.jl")
println("Running EA5_GraphQueryTEST")
include("2Apps/1EntityAnalysis/EA5_GraphQueryTEST.jl")
include("2Apps/2TrackAnalysis/findtrackgraph.jl")
include("2Apps/2TrackAnalysis/findtracks.jl")
println("Running TA1_BuildTEST")
include("2Apps/2TrackAnalysis/TA1_BuildTEST.jl")
println("Running TA2_QueryTEST")
include("2Apps/2TrackAnalysis/TA2_QueryTEST.jl")
println("Running TA3_GraphTEST")
include("2Apps/2TrackAnalysis/TA3_GraphTEST.jl")

# Run Parallel Database Example
println("Running pDB01_DataTEST")
include("3Scaling/2ParallelDatabase/pDB01_DataTEST.jl")
println("Running pDB02_FileTEST")
include("3Scaling/2ParallelDatabase/pDB02_FileTEST.jl")
println("Running pDB03_AssocTEST")
include("3Scaling/2ParallelDatabase/pDB03_AssocTEST.jl")
println("Running pDB04_DegreeTEST")
include("3Scaling/2ParallelDatabase/pDB04_DegreeTEST.jl")
println("Running pDB05_DBsetupTEST")
include("3Scaling/2ParallelDatabase/pDB05_DBsetupTEST.jl")
println("Running pDB06_AdjInsertTEST")
include("3Scaling/2ParallelDatabase/pDB06_AdjInsertTEST.jl")
println("Running pDB07_AdjQueryTEST")
include("3Scaling/2ParallelDatabase/pDB07_AdjQueryTEST.jl")
println("Running pDB08_AdjQueryItTEST")
include("3Scaling/2ParallelDatabase/pDB08_AdjQueryItTEST.jl")
println("Running pDB09_AdjJoinTEST")
include("3Scaling/2ParallelDatabase/pDB09_AdjJoinTEST.jl")
# Run Matrix Performance Example
println("Running MP1_DenseTEST")
include("3Scaling/3MatrixPerformance/MP1_DenseTEST.jl")
println("Running MP2_SparseTEST")
include("3Scaling/3MatrixPerformance/MP2_SparseTEST.jl")
println("Running MP3_AssocTEST")
include("3Scaling/3MatrixPerformance/MP3_AssocTEST.jl")
println("Running MP4_AssocCatKeyTEST")
include("3Scaling/3MatrixPerformance/MP4_AssocCatKeyTEST.jl")
println("Running MP5_AssocCatValKeyTEST")
include("3Scaling/3MatrixPerformance/MP5_AssocCatValKeyTEST.jl")
println("Running MP6_AssocPlus")
include("3Scaling/3MatrixPerformance/MP6_AssocPlus.jl")


text:
using Plots
using SpikingNeuralNetworks
SNN.@load_units

RS = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -65, d = 8))
IB = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -55, d = 4))
CH = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -50, d = 2))
FS = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.1, b = 0.2, c = -65, d = 2))
TC1 = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.25, c = -65, d = 0.05))
TC2 = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.25, c = -65, d = 0.05))
RZ = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.1, b = 0.26, c = -65, d = 2))
LTS = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.1, b = 0.25, c = -65, d = 2))
P = [RS, IB, CH, FS, TC1, TC2, RZ, LTS]

SNN.monitor(P, [:v])
T = 2second
for t = 0:T
    for p in [RS, IB, CH, FS, LTS]
        p.I = [10]
    end
    TC1.I = [(t < 0.2T) ? 0mV : 2mV]
    TC2.I = [(t < 0.2T) ? -30mV : 0mV]
    RZ.I =  [(0.5T < t < 0.6T) ? 10mV : 0mV]
    SNN.sim!(P, [], 0.1ms)
end
SNN.vecplot(P, :v)


text:
# Packages used in Script
using Colors
using Interpolations
using LazySets
using Luxor

# Auxiliary drawing functions
function make_drawing(width, height, img_path, bkg_color, origin_p)
    d = Drawing(width, height, img_path)
    background(bkg_color)
    origin(origin_p)
    return d
end

map2luxor(p) = Point(p.x, -p.y)
map2luxor(p) = Point(p...)

# Constants
## Canvas constants
width = 500
height = 500
path = "voronoi.png"
color = RGBA(0.5, 0, 0.5)
op = Point(width / 2, height / 2)

# Setting up drawing
my_draw = make_drawing(width, height, path, color, op)

# Mathematical formulae for calculations
midpoint(x₁, y₁, x₂, y₂) = (x = (x₁ + x₂) / 2, y = (y₁ + y₂) / 2)
parametric_point(t, x, y) = x + y * t
d(x1, y1, x2, y2, x, y) = (x - x1) * (y2 - y1) - (y - y1) * (x2 - x1)

## Point constants
p1 = [125; 125]
p2 = [170; 10]
p3 = [-200; 0]
p4 = [-30; -140]
p5 = [0; -5]

sethue(RGBA(1, 1, 0))
circle(map2luxor(p1), 15; action = :fill)
circle(map2luxor(p2), 15; action = :fill)
circle(map2luxor(p3), 15; action = :fill)
circle(map2luxor(p4), 15; action = :fill)
circle(map2luxor(p5), 15; action = :fill)

pts_list = [p1, p2, p3, p4, p5]
border_seeds = convex_hull(pts_list)
inner_seeds = setdiff(pts_list, border_seeds)

inner_hulls = []

for pt in inner_seeds
    edges = []
    sub_pts_list = [x for x in pts_list if x != pt]
    for s in sub_pts_list
        p_vec = pt - s
        mp = midpoint(s[1], s[2], pt[1], pt[2])

        x = 1
        y = 1
        if (-1 * p_vec[1] * x) / p_vec[2] |> isinf
            x = (p_vec[2] * y) / (-1 * p_vec[1])
        else
            y = (-1 * p_vec[1] * x) / p_vec[2]
        end

        p = [parametric_point(i, [mp.x; mp.y], [x; y]) for i = -10:10]
        x = [i[1] for i in p]
        y = [i[2] for i in p]
        ext =
            interpolate((x,), y, Gridded(Linear())) |>
            foo -> extrapolate(foo, Interpolations.Line())

        sethue(RGBA(1, 1, 1))
        push!(edges, [map2luxor([-250, ext(-250)]), map2luxor([250, ext(250)])])
        # Luxor.line(
        # map2luxor([-250, ext(-250)]),
        # map2luxor([250, ext(250)]);
        # action = :stroke,
        # )
    end

    vtx = []
    for e1 in edges
        for e2 in edges
            if e1 != e2
                push!(vtx, (intersectionlines(e1[1], e1[2], e2[1], e2[2]))[2])
            end
        end
    end

    # TODO: Add line check here

    vtx = vtx |> unique
    distances = Luxor.distance.(vtx, Point(pt...))
    v_vtx = vtx[sortperm(distances)][1:4] |> x -> convert(Vector{Point}, x)
    hull = [[p.x, p.y] for p in v_vtx] |> convex_hull
    lux_hull = hull |> x -> [Point(p[1], p[2]) for p in x]

    sethue(RGBA(1.0, 0.7529411764705882, 0.796078431372549))
    poly(lux_hull; close = true, action = :fill)

    sethue(RGBA(0, 0, 0))
    poly(lux_hull; close = true, action = :stroke)

    push!(inner_hulls, hull)

    # TODO: Determine how to get canvas edge points

end

for pt in border_seeds
    edges = [
        [map2luxor([-250, -250]), map2luxor([250, -250])],
        [map2luxor([-250, 250]), map2luxor([250, 250])],
        [map2luxor([-250, -250]), map2luxor([-250, 250])],
        [map2luxor([250, -250]), map2luxor([250, 250])],
    ]

    sethue(RGBA(1, 1, 1))
    # Top edge
    Luxor.line(map2luxor([-250, -250]), map2luxor([250, -250]); action = :stroke)

    # Bottom edge
    Luxor.line(map2luxor([-250, 250]), map2luxor([250, 250]); action = :stroke)

    # Left edge
    Luxor.line(map2luxor([-250, -250]), map2luxor([-250, 250]); action = :stroke)

    # Right edge
    Luxor.line(map2luxor([250, -250]), map2luxor([250, 250]); action = :stroke)

    sub_pts_list = [x for x in pts_list if x != pt]
    for s in sub_pts_list
        p_vec = pt - s
        mp = midpoint(s[1], s[2], pt[1], pt[2])

        x = 1
        y = 1
        if (-1 * p_vec[1] * x) / p_vec[2] |> isinf
            x = (p_vec[2] * y) / (-1 * p_vec[1])
        else
            y = (-1 * p_vec[1] * x) / p_vec[2]
        end

        p = [parametric_point(i, [mp.x; mp.y], [x; y]) for i = -10:10]
        x = [i[1] for i in p]
        y = [i[2] for i in p]
        ext =
            interpolate((x,), y, Gridded(Linear())) |>
            foo -> extrapolate(foo, Interpolations.Line())

        sethue(RGBA(1, 1, 1))
        push!(edges, [map2luxor([-250, ext(-250)]), map2luxor([250, ext(250)])])
        Luxor.line(
            map2luxor([-250, ext(-250)]),
            map2luxor([250, ext(250)]);
            action = :stroke,
        )
    end

    vtx = []
    inner_vtx = []
    for e1 in edges
        for e2 in edges
            if e1 != e2
                v = (intersectionlines(e1[1], e1[2], e2[1], e2[2]))[2]
                push!(vtx, v)
                if v in inner_hulls[1] .|> map2luxor
                    push!(inner_vtx, v)
                end
            end
        end
    end

    # TODO: Add line check here

    sethue(RGBA(1.0, 0.7529411764705882, 0.796078431372549))

    vtx = vtx |> unique |> x -> [[p.x, p.y] for p in x]
    v_vtx = []
    for p in vtx
        if !(element(Singleton(p)) ∈ VPolygon(inner_hulls[1]))
            push!(v_vtx, p)
        end
    end

    v_vtx = vcat(v_vtx, inner_vtx |> unique)
    filt_v_vtx = []
    for p in v_vtx
    println(p)
	if (-width <= p[1] <= width) && (-height <= p[2] <= height)
		push!(filt_v_vtx, p)
	end
    end

    distances = filt_v_vtx .|> x -> Point(x...) |> x -> Luxor.distance.(x, Point(pt...))
    filt_v_vtx = filt_v_vtx .|> x -> Point(x...)
    clipboard(filt_v_vtx)
    filt_v_vtx = filt_v_vtx[sortperm(distances)][1:5]
    hull = [[p.x, p.y] for p in filt_v_vtx] |> convex_hull |> x -> [Point(p[1], p[2]) for p in x]

    sethue(RGBA(1.0, 0.7529411764705882, 0.796078431372549))
    poly(hull; close = true, action = :fill)

    sethue(RGBA(0, 0, 0))
    poly(hull; close = true, action = :stroke)

end

finish()

max_stars_repo_name: TheCedarPrince/DataExplorers

max_issues_repo_name: TheCedarPrince/DataExplorers

max_forks_repo_name: TheCedarPrince/DataExplorers


text:
function choropleth(; x::AbstractVector = Int[], y::AbstractVector = Int[], entity::Symbol = :uscounties)
    v = VegaVisualization(name = "choropleth", width = 960, height = 500)

    #Add numeric data
    add_data!(v, x = x, y = y)

    #Get unique table name
    table = v.data[1].name

    #Add geodata

    #US County level map
    if entity == :uscounties
      push!(v.data, VegaData(name = "geodata",
                             url = "http://trifacta.github.io/vega/data/us-10m.json",
                             format = VegaFormat(_type = "topojson", feature = "counties"),
                             transform = [VegaTransform(Dict{Any, Any}("type" => "geopath", "projection" => "albersUsa")),
                                          VegaTransform(Dict{Any, Any}("type" => "lookup", "on" => table, "onKey" => "x", "keys" => ["id"], "as" => ["table2"])),
                                          VegaTransform(Dict{Any, Any}("type" => "filter", "test" => "datum.layout_path!=null && datum.table2!=null"))]
                            )
           )

    #United States map
    elseif entity == :usstates

            push!(v.data, VegaData(name = "geodata",
                             url = "http://trifacta.github.io/vega/data/us-10m.json",
                             format = VegaFormat(_type = "topojson", feature = "states"),
                             transform = [VegaTransform(Dict{Any, Any}("type" => "geopath", "projection" => "albersUsa")),
                                          VegaTransform(Dict{Any, Any}("type" => "lookup", "on" => table, "onKey" => "x", "keys" => ["id"], "as" => ["table2"])),
                                          VegaTransform(Dict{Any, Any}("type" => "filter", "test" => "datum.layout_path!=null && datum.table2!=null"))]
                            )
           )

    end

    #Generically applicable to all data
    v.scales = [VegaScale(name = "group",
                         _type = "quantize",
                          domain = [minimum(y), maximum(y)],
                          range = ["#f7fbff", "#deebf7", "#c6dbef", "#9ecae1", "#6baed6", "#4292c6", "#2171b5", "#08519c", "#08306b"]
                        )]

    v.marks = [VegaMark(_type = "path",
                        from = VegaMarkFrom(data="geodata"),
                        properties = VegaMarkProperties(enter = VegaMarkPropertySet(path = VegaValueRef(field = "layout_path"), fill = VegaValueRef(scale = "group", field = "table2.y"))

                                                        )
            )]

    return v

end

max_stars_repo_name: jpthor/vega