text:
#=
Compute the condition number given the initial graph and the preconditioner.
Uses cholfact, so, if the initial graph is not sdd or if inverse has a lot of nonzeros this fails.
=#
function condNumber{Tv,Ti}(la::SparseMatrixCSC{Tv,Ti}, prec::SparseMatrixCSC{Tv,Ti}; tol = 1e-5)
	# cf = cholfact(prec)
	n = la.n
	ldl = ldltfact(prec)
	diag = 1 ./ (ldl[:D] \ ones(n))

	return computeCN(la, ldl[:UP], diag, tol=tol)
end

# M = Ut * d * U
function condNumber{Tv,Ti}(la::SparseMatrixCSC{Tv,Ti}, U::UpperTriangular{Tv,SparseMatrixCSC{Tv,Ti}}, d::Array{Tv,1}; tol = 1e-5)
	return computeCN(la, U, d, tol=tol)
end

# M = Ut * U
function condNumber{Tv,Ti}(la::SparseMatrixCSC{Tv,Ti}, U::UpperTriangular{Tv,SparseMatrixCSC{Tv,Ti}}; tol = 1e-5)
	return computeCN(la, U, ones(a.n), tol=tol)
end

function computeCN{Tv,Ti}(la::SparseMatrixCSC{Tv,Ti}, U, d::Array{Tv,1}; tol=1e-5)

	n = la.n

	# get the max eigenvalue of d^-1/2 * L^-1 * a * U^-1 * d^-1/2
	f = function(b::Array{Float64,1})
		res = copy(b)

		# diag sqrt inverse
		for i in 1:(n - 1)
            res[i] = res[i] * d[i]^(-1/2)
        end
		res[n] = 0

		res = U \ res
		res = la * res
		res = U' \ res

		# diag sqrt inverse
        for i in 1:(n - 1)
            res[i] = res[i] * d[i]^(-1/2)
        end
        res[n] = 0

		return res
	end
	fOp = SqLinOp(true,1.0,n,f)
	lmax = eigs(fOp;nev=1,which=:LM,tol=0.0)[1][1]

	if lmax < 1e-20
		# If lambdamax is too small, the computation of 1 - 1/K won't be possible.
		# So, just return a huge condition number (which should be true anyway)
		return 1,1e20
	end

	# get the max eigenvalue of 1 - (M^-1/2 L M^-1/2) / lambdaMax
	g = function(b::Array{Float64,1})
		res2 = copy(b)
		res = copy(b) / lmax

		# diag sqrt inverse
		for i in 1:(n - 1)
            res[i] = res[i] * d[i]^(-1/2)
        end
		res[n] = 0

		res = U \ res
		res = la * res
		res = U' \ res

		# diag sqrt inverse
		for i in 1:(n - 1)
            res[i] = res[i] * d[i]^(-1/2)
        end
		res[n] = res2[n] = 0

		return res2 - res
	end
	gOp = SqLinOp(true,1.0,n,g)

    #R = eigs(gOp;nev=1,which=:LM,tol=tol)[1][1]
	R = powerIteration(gOp, tol=tol, maxit=100)

    Kmin = 1 / (1 - R)
    Kmax = 1 / (1 - R - tol)

    if Kmax < 0
    	while Kmax < 0
    		tol = tol / 10

    		R = powerIteration(gOp, tol=tol)

		    Kmin = 1 / (1 - R)
		    Kmax = 1 / (1 - R - tol)
    	end

    	println("Initial tolerance too large. Recomputed with tolerance = ", tol)
    end

    return (Kmin, Kmax)

end

# function computeCN{Tv,Ti}(la::SparseMatrixCSC{Tv,Ti}, U::UpperTriangular{Tv,SparseMatrixCSC{Tv,Ti}}, d::Array{Tv,1}; tol=1e-3)

# 	n = la.n

# 	# get the max eigenvalue of d^-1/2 * L^-1 * a * U^-1 * d^-1/2
# 	f = function(b::Array{Float64,1})
# 		res = copy(b)

# 		# diag sqrt inverse
# 		res = res .* d.^-0.5
# 		res[n] = 0

# 		res = U \ res
# 		res = la * res
# 		res = U' \ res

# 		# diag sqrt inverse
#         res = res .* d.^-0.5
#         res[n] = 0

# 		return res
# 	end
# 	fOp = SqLinOp(true,1.0,n,f)
# 	lmax = eigs(fOp;nev=1,which=:LM,tol=0.0)[1][1]

# 	# get the max eigenvalue of 1 - (M^-1/2 L M^-1/2) / lambdaMax
# 	g = function(b::Array{Float64,1})
# 		res2 = copy(b)
# 		res = copy(b) / lmax

# 		# diag sqrt inverse
# 		res = res .* d.^-0.5
# 		res[n] = 0

# 		res = U \ res
# 		res = la * res
# 		res = U' \ res

# 		# diag sqrt inverse
# 		res = res .* d.^-0.5
# 		res[n] = res2[n] = 0

# 		return res2 - res
# 	end
# 	gOp = SqLinOp(true,1.0,n,g)

# 	R = eigs(gOp;nev=1,which=:LM,tol=tol)[1][1]

#     Kmin = 1 / (1 - R)
#     Kmax = 1 / (1 - R - tol)

#     return (Kmin, Kmax)

# end


text:
module TestSetup

using GLFW

include("init_test.jl")

end

ext: jl


text:
##########
# TArray #
##########

"""
    TArray{T}(dims, ...)

Implementation of data structures that automatically perform copy-on-write after task copying.

If current_task is an existing key in `s`, then return `s[current_task]`. Otherwise, return `s[current_task] = s[last_task]`.

Usage:

```julia
TArray(dim)
```

Example:

```julia
ta = TArray(4)              # init
for i in 1:4 ta[i] = i end  # assign
Array(ta)                   # convert to 4-element Array{Int64,1}: [1, 2, 3, 4]
```
"""
struct TArray{T,N} <: AbstractArray{T,N}
    ref :: Symbol  # object_id
    orig_task :: Task
    TArray{T,N}() where {T,N} = new(gensym(), current_task())
end

TArray{T,1}(d::Integer) where T = TArray(T,  d)
TArray{T}(d::Integer...) where T = TArray(T, d)
TArray{T}(UndefInitializer, d::Integer...) where T = TArray(T, d)
TArray{T,N}(d::Integer...) where {T,N} = length(d)==N ? TArray(T,d) : error("Malformed dims")
TArray{T,N}(UndefInitializer, d::Integer...) where {T,N} = length(d)==N ? TArray(T,d) : error("Malformed dims")
TArray{T,N}(dim::NTuple{N,Int}) where {T,N} = TArray(T, dim)

function TArray(T::Type, dim)
    res = TArray{T,length(dim)}();
    n = n_copies()
    d = Array{T}(undef, dim)
    task_local_storage(res.ref, (n,d))
    res
end

#
# Indexing Interface Implementation
#

function Base.getindex(S::TArray{T, N}, I::Vararg{Int,N}) where {T, N}
    t, d = task_local_storage(S.ref)
    return d[I...]
end

function Base.setindex!(S::TArray{T, N}, x, I::Vararg{Int,N}) where {T, N}
    n, d = task_local_storage(S.ref)
    cn   = n_copies()
    newd = d
    if cn > n
        # println("[setindex!]: $(S.ref) copying data")
        newd = deepcopy(d)
        task_local_storage(S.ref, (cn, newd))
    end
    newd[I...] = x
end

function Base.push!(S::TArray{T}, x) where T
    n, d = task_local_storage(S.ref)
    cn   = n_copies()
    newd = d
    if cn > n
        newd = deepcopy(d)
        task_local_storage(S.ref, (cn, newd))
    end
    push!(newd, x)
end

function Base.pop!(S::TArray)
    n, d = task_local_storage(S.ref)
    cn   = n_copies()
    newd = d
    if cn > n
        newd = deepcopy(d)
        task_local_storage(S.ref, (cn, newd))
    end
    pop!(d)
end

function Base.convert(::Type{TArray}, x::Array)
    res = TArray{typeof(x[1]),ndims(x)}();
    n   = n_copies()
    task_local_storage(res.ref, (n,x))
    return res
end

function Base.convert(::Array, x::Type{TArray})
    n,d = task_local_storage(S.ref)
    c = deepcopy(d)
    return c
end

function Base.display(S::TArray)
    arr = S.orig_task.storage[S.ref][2]
    @warn "display(::TArray) prints the originating task's storage, " *
        "not the current task's storage. " *
        "Please use show(::TArray) to display the current task's version of a TArray."
    display(arr)
end

Base.show(io::IO, S::TArray) = Base.show(io::IO, task_local_storage(S.ref)[2])

# Base.get(t::Task, S) = S
# Base.get(t::Task, S::TArray) = (t.storage[S.ref][2])
Base.get(S::TArray) = (current_task().storage[S.ref][2])

##
# Iterator Interface
IteratorSize(::Type{TArray{T, N}}) where {T, N} = HasShape{N}()
IteratorEltype(::Type{TArray}) = HasEltype()

# Implements iterate, eltype, length, and size functions,
# as well as firstindex, lastindex, ndims, and axes
for F in (:iterate, :eltype, :length, :size,
          :firstindex, :lastindex, :ndims, :axes)
    @eval Base.$F(a::TArray, args...) = $F(get(a), args...)
end

#
# Similarity implementation
#

Base.similar(S::TArray) = tzeros(eltype(S), size(S))
Base.similar(S::TArray, ::Type{T}) where {T} = tzeros(T, size(S))
Base.similar(S::TArray, dims::Dims) = tzeros(eltype(S), dims)

##########
# tzeros #
##########

"""
     tzeros(dims, ...)

Construct a distributed array of zeros.
Trailing arguments are the same as those accepted by `TArray`.

```julia
tzeros(dim)
```

Example:

```julia
tz = tzeros(4)              # construct
Array(tz)                   # convert to 4-element Array{Int64,1}: [0, 0, 0, 0]
```
"""
function tzeros(T::Type, dim)
    res = TArray{T,length(dim)}();
    n = n_copies()
    d = zeros(T,dim)
    task_local_storage(res.ref, (n,d))
    return res
end

tzeros(::Type{T}, d1::Integer, drest::Integer...) where T = tzeros(T, convert(Dims, tuple(d1, drest...)))
tzeros(d1::Integer, drest::Integer...) = tzeros(Float64, convert(Dims, tuple(d1, drest...)))
tzeros(d::Dims) = tzeros(Float64, d)

"""
     tfill(val, dim, ...)

Construct a TArray of a specified value.

```julia
tfill(val, dim)
```

Example:

```julia
tz = tfill(9.0, 4)            # construct
Array(tz)                     # convert to 4-element Array{Float64,1}:  [9.0  9.0  9.0  9.0]
```
"""
function tfill(val::Real, dim)
    res = TArray{typeof(val),length(dim)}();
    n = n_copies()
    d = fill(val,dim)
    task_local_storage(res.ref, (n,d))
    return res
end


text:
module DictWrappers

export DictWrapper,
    ShallowDictWrapper,
    RecursiveDictWrapper

const DICT_WRAPPER_ERROR = ArgumentError("DictWrappers can only wrap composite types.")

iscomposite(t::DataType) = Base.isstructtype(t)
iscomposite{T<:Tuple}(t::Type{T}) = false

abstract DictWrapper <: Associative

immutable ShallowDictWrapper <: DictWrapper
    wrapped

    function ShallowDictWrapper(wrapped)
        if iscomposite(typeof(wrapped))
            return new(wrapped)
        else
            throw(DICT_WRAPPER_ERROR)
        end
    end
end

immutable RecursiveDictWrapper <: DictWrapper
    wrapped

    function RecursiveDictWrapper(wrapped)
        if iscomposite(typeof(wrapped))
            return new(wrapped)
        else
            throw(DICT_WRAPPER_ERROR)
        end
    end
end

Base.getindex(dw::ShallowDictWrapper, key::Symbol) = getfield(dw.wrapped, key)

function Base.getindex(dw::RecursiveDictWrapper, key::Symbol)
    val = getfield(dw.wrapped, key)

    if iscomposite(typeof(val))
        return RecursiveDictWrapper(val)
    else
        return val
    end
end

Base.length(dw::DictWrapper) = length(fieldnames(dw.wrapped))

Base.keys(dw::DictWrapper) = fieldnames(dw.wrapped)

immutable DictWrapperValueIterator{T<:DictWrapper}
    dw::T
    keys::Vector{Symbol}
end

DictWrapperValueIterator(dw::DictWrapper) = DictWrapperValueIterator(dw, keys(dw))

Base.values(dw::DictWrapper) = DictWrapperValueIterator(dw)

Base.start(dwvi::DictWrapperValueIterator) = 1
Base.next(dwvi::DictWrapperValueIterator, state) = (dwvi.dw[dwvi.keys[state]], state + 1)
Base.done(dwvi::DictWrapperValueIterator, state) = state > length(dwvi.keys)

function Base.start(dw::DictWrapper)
    return (fieldnames(dw.wrapped), 1)
end

function Base.next(dw::DictWrapper, state)
    fields, idx = state
    curr_field = fields[idx]
    return (curr_field=>dw[curr_field], (fields, idx + 1))
end

function Base.done(dw::DictWrapper, state)
    fields, idx = state
    return idx > length(fields)
end

function Base.in{T}(item::Pair{Symbol, T}, dw::DictWrapper)
    return item.first in keys(dw) && item.second == dw[item.first]
end

function Base.get(dw::DictWrapper, key::Symbol, default)
    if key in keys(dw)
        return dw[key]
    else
        return default
    end
end

end # module


text:
using GPUArrays
using GPUArrays.TestSuite
using Base.Test

function run_fft(Typ)
    T = Typ{Complex64}
    for n = 1:3
        @testset "FFT with ND = $n" begin
            dims = ntuple(i-> 40, n)
            against_base(fft!, T, dims)
            against_base(ifft!, T, dims)

            against_base(fft, T, dims)
            against_base(ifft, T, dims)
        end
    end
end


text:
struct LHSMatrix

    M :: SparseMatrixCSC{Float64, Int64}
    Mdist :: DArray{Float64, 2, SparseMatrixCSC{Float64, Int64}}

end

function LHSMatrix(A::SparseMatrixCSC{Float64, Int64}; prols=workers())

    Mdist = DArray(size(A), workers(), (1,length(workers()))) do inds
        j = inds[2]
        A[:,j]
    end

    return LHSMatrix(A, Mdist)
end

function solve_local(M::DArray{Float64, 2, SparseMatrixCSC{Float64, Int64}}, r::Vector{Float64})

    delta = M[:L] \ r
    inds = localindices(M)[2]

    return delta, inds
end

function dlinsolve(L::LHSMatrix, x0::Vector{Float64}, b::Vector{Float64}, Np)

    n = length(x0)
    pids = procs(L.Mdist)
    nworkers = length(pids)
    chunks = collect(Iterators.partition(1:nworkers, Np))
    slices = fetch.([@spawnat pid localindices(L.Mdist)[2] for pid in pids])

    r = b - L.M*x0
    x = x0[:]

    bnorm = norm(b)
    err = norm(r) / bnorm
    while err > 1e-10
        for chunk in chunks

            res = Future[]
            for pid in pids[chunk]
                fut = @spawnat pid solve_local(L.Mdist, r)
                append!(res, [fut])
            end

            for f in res
                delta, inds = fetch(f)
                x[inds] += delta
            end

            r[:] = b - L.M*x
        end

        err = norm(r) / bnorm
        println(err)
    end

    return x
end


text:
#=Copyright (c) 2017 Gabriel Goh

Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in
the Software without restriction, including without limitation the rights to
use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
the Software, and to permit persons to whom the Software is furnished to do so,
subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.=#

# α = 0.076
# β = 0.513

α = 0.2
β = 0.99
x = 2
y = 0
X = zeros(100)
for i = 1:100
    y  = β*y - x
    x  = x - α*y
    X[i] = x
end
figure(); plot(X,"k"); plot(X,"k."); xlabel(L"k"); ylabel(L"f(x)")

max_stars_repo_name: sbeleidy/post--momentum

max_issues_repo_name: sbeleidy/post--momentum

max_forks_repo_name: sbeleidy/post--momentum


text:
Adapt.adapt_storage(to::PyAdaptor{T}, x::CUDA.AnyCuArray) where {T} = CUDA.CuArray(x)
Adapt.adapt_storage(to::PyAdaptor{T}, x::CUDA.StridedCuArray) where {T} = x
Adapt.adapt_storage(to::PyAdaptor{<:CUDA.CuArray}, x::FillArrays.AbstractFill) = CUDA.cu(collect(x))
Adapt.adapt_structure(to::PyAdaptor{<:CUDA.CuArray}, x::SubArray{F,I,G,H,L}) where {F,I,G <: CUDA.CuArray,H,L}  = copy(x)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "lttng_ust"
version = v"2.12.2"

# Collection of sources required to complete build
sources = [
    ArchiveSource("http://lttng.org/files/lttng-ust/lttng-ust-2.12.2.tar.bz2", "bcd0f064b6ca88c72d84e760eac3472ae5c828411c634435922bee9fce359fc7")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/lttng-ust*
export CPPFLAGS="-I${includedir}"
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --disable-man-pages --disable-examples
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(Sys.islinux, supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("liblttng-ust-fd", :liblttng_ust_fd),
    LibraryProduct("liblttng-ust-dl", :liblttng_ust_dl),
    LibraryProduct("liblttng-ust-fork", :liblttng_ust_fork),
    LibraryProduct("liblttng-ust-libc-wrapper", :liblttng_ust_libc_wrapper),
    LibraryProduct("liblttng-ust-cyg-profile-fast", :liblttng_ust_cyg_profile_fast),
    LibraryProduct("liblttng-ust", :liblttng_ust),
    LibraryProduct("liblttng-ust-tracepoint", :liblttng_ust_tracepoint),
    LibraryProduct("liblttng-ust-pthread-wrapper", :liblttng_ust_pthread_wrapper),
    LibraryProduct("liblttng-ust-cyg-profile", :liblttng_ust_cyg_profile),
    LibraryProduct("liblttng-ust-ctl", :liblttng_ust_ctl)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="NUMA_jll", uuid="7f51dc2b-bb24-59f8-b771-bb1490e4195d"))
    Dependency(PackageSpec(name="URCU_jll", uuid="aa747835-a391-587f-982f-064ff03f7d29"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"5")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using BeamFEA
using Test

@testset "Simple Beam Deflection Tests" begin

    # --- Test data from "Finite Element Structural Analysis", Yang, pg. 145 ----
    # cantilever_deflection

    # arbitrary values
    E = 2.0
    I = 3.0
    L = 4.0
    p0 = 5.0

    x = [0, L]
    EIy = [E*I, E*I]
    EIz = [E*I, E*I]
    EA = [E*1.0, E*1.0]
    GJ = [1.0, 1.0]
    rhoA = [1.0, 1.0]
    rhoJ = [1.0, 1.0]
    beam = BeamProperties(x, EIy, EIz, EA, GJ, rhoA, rhoJ)

    Px = [0.0, 0.0]
    Py = [0.0, -p0]
    Pz = [0.0, 0.0]
    Fx = [0.0, 0.0]
    Fy = [0.0, 0.0]
    Fz = [0.0, 0.0]
    Mx = [0.0, 0.0]
    My = [0.0, 0.0]
    Mz = [0.0, 0.0]
    loads = Loads(x, Px, Py, Pz, Fx, Fy, Fz, Mx, My, Mz)

    kx = [0.0, Inf]
    ky = [0.0, Inf]
    kz = [0.0, Inf]
    kthetax = [0.0, Inf]
    kthetay = [0.0, Inf]
    kthetaz = [0.0, Inf]
    stiffness = Stiffness(kx, ky, kz, kthetax, kthetay, kthetaz)

    defl, freq = fea(beam, loads, stiffness)

    @test isapprox(defl.y[1], -p0*L^4/(30*E*I), atol=1e-15)
    @test isapprox(defl.thetaz[1], -p0*L^3/(24*E*I), atol=1e-15)
    @test defl.x[1] == 0.0
    @test defl.z[1] == 0.0
    @test defl.thetax[1] == 0.0
    @test defl.thetay[1] == 0.0

    # --- Test data from "Finite Element Structural Analysis", Yang, pg. 180 ----
    # tapered beam

    # arbitrary values
    E = 2.0
    I = 3.0
    L = 4.0
    P = 5.0

    x = [0, L/2, L]
    EIy = [E*I, 5*E*I, 9*E*I]
    EIz = [E*I, 5*E*I, 9*E*I]
    EA = [1.0, 1.0, 1.0]
    GJ = [1.0, 1.0, 1.0]
    rhoA = [1.0, 1.0, 1.0]
    rhoJ = [1.0, 1.0, 1.0]
    beam = BeamProperties(x, EIy, EIz, EA, GJ, rhoA, rhoJ)

    Px = [0.0, 0.0, 0.0]
    Py = [0.0, 0.0, 0.0]
    Pz = [0.0, 0.0, 0.0]
    Fx = [0.0, 0.0, 0.0]
    Fy = [-P, 0.0, 0.0]
    Fz = [0.0, 0.0, 0.0]
    Mx = [0.0, 0.0, 0.0]
    My = [0.0, 0.0, 0.0]
    Mz = [0.0, 0.0, 0.0]
    loads = Loads(x, Px, Py, Pz, Fx, Fy, Fz, Mx, My, Mz)

    kx = [0.0, 0.0, Inf]
    ky = [0.0, 0.0, Inf]
    kz = [0.0, 0.0, Inf]
    kthetax = [0.0, 0.0, Inf]
    kthetay = [0.0, 0.0, Inf]
    kthetaz = [0.0, 0.0, Inf]
    stiffness = Stiffness(kx, ky, kz, kthetax, kthetay, kthetaz)

    defl, freq = fea(beam, loads, stiffness)

    @test isapprox(defl.y[1], -0.051166*P*L^3/(E*I), atol=1e-2)
    @test isapprox(defl.thetaz[1], -0.090668*P*L^2/(E*I), atol=1e-2)
    @test defl.x[1] == 0.0
    @test defl.z[1] == 0.0
    @test defl.thetax[1] == 0.0
    @test defl.thetay[1] == 0.0

end  # end unit test

@testset "Beam Frequency Tests" begin

    # Test data from "Consistent Mass Matrix for Distributed Mass Systmes", John Archer,
    # Journal of the Structural Division Proceedings of the American Society of Civil Engineers,
    # pg. 168
    E = 2.0
    I = 3.0
    L = 4.0
    A = 5.0
    rho = 6.0

    n = 1
    nodes = n + 1
    x = range(0, n*L, length=nodes)
    EIy = E*I*ones(nodes)
    EIz = E*I*ones(nodes)
    EA = ones(nodes)
    GJ = ones(nodes)
    rhoA = rho*A*ones(nodes)
    rhoJ = ones(nodes)
    beam = BeamProperties(x, EIy, EIz, EA, GJ, rhoA, rhoJ)

    loads = NoLoads(x)

    kx = zeros(nodes)
    ky = zeros(nodes)
    kz = zeros(nodes)
    kthetax = zeros(nodes)
    kthetay = zeros(nodes)
    kthetaz = zeros(nodes)
    stiffness = Stiffness(kx, ky, kz, kthetax, kthetay, kthetaz)

    delta, freq = fea(beam, loads, stiffness)

    alpha = rho * A * (n*L)^4 / (840.0 * E * I)

    @test isapprox(freq[2], sqrt(0.85714 / alpha) / (2*pi), atol=1e-6)
    @test isapprox(freq[3], sqrt(0.85714 / alpha) / (2*pi), atol=1e-6)
    @test isapprox(freq[5], sqrt(10.0 / alpha) / (2*pi), atol=1e-6)
    @test isapprox(freq[6], sqrt(10.0 / alpha) / (2*pi), atol=1e-6)

    # simply supported
    kx[1] = Inf
    ky[1] = Inf
    kz[1] = Inf
    ky[end] = Inf
    stiffness = Stiffness(kx, ky, kz, kthetax, kthetay, kthetaz)

    delta, freq = fea(beam, loads, stiffness)

    @test isapprox(freq[2], sqrt(0.14286 / alpha) / (2*pi), atol=1e-6)
    @test isapprox(freq[5], sqrt(3.0 / alpha) / (2*pi), atol=1e-6)

    # second test from same paper (n = 2)

    n = 2
    nodes = n + 1
    x = range(0, n*L, length=nodes)
    EIy = E*I*ones(nodes)
    EIz = E*I*ones(nodes)
    EA = ones(nodes)
    GJ = ones(nodes)
    rhoA = rho*A*ones(nodes)
    rhoJ = ones(nodes)
    beam = BeamProperties(x, EIy, EIz, EA, GJ, rhoA, rhoJ)

    loads = NoLoads(x)

    kx = zeros(nodes)
    ky = zeros(nodes)
    kz = zeros(nodes)
    kthetax = zeros(nodes)
    kthetay = zeros(nodes)
    kthetaz = zeros(nodes)
    stiffness = Stiffness(kx, ky, kz, kthetax, kthetay, kthetaz)

    delta, freq = fea(beam, loads, stiffness)

    alpha = rho * A * (n*L)^4 / (840.0 * E * I)

    @test isapprox(freq[2], sqrt(0.59858 / alpha) / (2*pi), atol=1e-6)
    @test isapprox(freq[3], sqrt(0.59858 / alpha) / (2*pi), atol=1e-6)
    @test isapprox(freq[6], sqrt(5.8629 / alpha) / (2*pi), atol=1e-6)
    @test isapprox(freq[7], sqrt(5.8629 / alpha) / (2*pi), atol=1e-6)
    @test isapprox(freq[9], sqrt(36.659 / alpha) / (2*pi), atol=2e-6)
    @test isapprox(freq[10], sqrt(36.659 / alpha) / (2*pi), atol=2e-6)
    @test isapprox(freq[11], sqrt(93.566 / alpha) / (2*pi), atol=1e-6)
    @test isapprox(freq[12], sqrt(93.566 / alpha) / (2*pi), atol=1e-6)

    # simply supported
    kx[1] = Inf
    ky[1] = Inf
    kz[1] = Inf
    ky[end] = Inf
    stiffness = Stiffness(kx, ky, kz, kthetax, kthetay, kthetaz)
    delta, freq = fea(beam, loads, stiffness)

    @test isapprox(freq[2], sqrt(0.11688 / alpha) / (2*pi), atol=1e-6)
    @test isapprox(freq[5], sqrt(2.2858 / alpha) / (2*pi), atol=1e-6)
    @test isapprox(freq[8], sqrt(14.441 / alpha) / (2*pi), atol=1e-5)
    @test isapprox(freq[11], sqrt(48 / alpha) / (2*pi), atol=1e-6)

    # third test from same paper (n = 3)

    n = 3
    nodes = n + 1
    x = range(0, n*L, length=nodes)
    EIy = E*I*ones(nodes)
    EIz = E*I*ones(nodes)
    EA = ones(nodes)
    GJ = ones(nodes)
    rhoA = rho*A*ones(nodes)
    rhoJ = ones(nodes)
    beam = BeamProperties(x, EIy, EIz, EA, GJ, rhoA, rhoJ)

    loads = NoLoads(x)

    kx = zeros(nodes)
    ky = zeros(nodes)
    kz = zeros(nodes)
    kthetax = zeros(nodes)
    kthetay = zeros(nodes)
    kthetaz = zeros(nodes)
    stiffness = Stiffness(kx, ky, kz, kthetax, kthetay, kthetaz)

    delta, freq = fea(beam, loads, stiffness)

    alpha = rho * A * (n*L)^4 / (840.0 * E * I)

    @test isapprox(freq[2], sqrt(0.59919 / alpha)/(2*pi), atol=1e-6)
    @test isapprox(freq[3], sqrt(0.59919 / alpha)/(2*pi), atol=1e-6)

    @test isapprox(freq[6], sqrt(4.5750 / alpha)/(2*pi), atol=1e-6)
    @test isapprox(freq[7], sqrt(4.5750 / alpha)/(2*pi), atol=1e-6)

    @test isapprox(freq[9], sqrt(22.010 / alpha)/(2*pi), atol=1e-6)
    @test isapprox(freq[10], sqrt(22.010 / alpha)/(2*pi), atol=1e-6)

    @test isapprox(freq[12], sqrt(70.920 / alpha)/(2*pi), atol=1e-6)
    @test isapprox(freq[13], sqrt(70.920 / alpha)/(2*pi), atol=1e-6)

    @test isapprox(freq[15], sqrt(265.91 / alpha)/(2*pi), atol=1e-5)
    @test isapprox(freq[16], sqrt(265.91 / alpha)/(2*pi), atol=1e-5)

    @test isapprox(freq[17], sqrt(402.40 / alpha)/(2*pi), atol=1e-6)
    @test isapprox(freq[18], sqrt(402.40 / alpha)/(2*pi), atol=1e-6)

end

@testset "shear and bending" begin

    # Test data from "Mechanical of Materials", Gere, 6th ed., pg. 273
    # cantilevered beam with linear distributed load

    L = 10.0
    q0 = 3.0

    n = 1
    nodes = n+1

    x = range(0, L, length=nodes)
    EIy = ones(nodes)
    EIz = ones(nodes)
    EA = ones(nodes)

    beam = BeamProperties(x, EIy, EIz, EA, zeros(nodes), zeros(nodes), zeros(nodes))

    Px = [0.0, 0.0]
    Py = [0.0, -q0]
    Pz = [0.0, 0.0]
    Fx = zeros(nodes)
    Fy = zeros(nodes)
    Fz = zeros(nodes)
    Mx = zeros(nodes)
    My = zeros(nodes)
    Mz = zeros(nodes)
    loads = Loads(x, Px, Py, Pz, Fx, Fy, Fz, Mx, My, Mz)

    y = zeros(nodes)
    z = zeros(nodes)

    epsilon, sb = strain(y, z, beam, loads)

    @test isapprox(sb.Vy[1], 0.0, atol=1e-6)
    @test isapprox(sb.Vy[2], -q0*L/2.0, atol=1e-6)
    @test isapprox(sb.Mz[1], 0.0, atol=1e-6)
    @test isapprox(sb.Mz[2], -q0*L^2/6, atol=1e-6)

end


text:
using ReinforcementLearningBase
using ReinforcementLearningEnvironments
using PyCall
using UnicodePlots

@testset "d4rl_policies" begin
    model = d4rl_policy("ant", "online", 10)

    @test typeof(model) <: D4RLGaussianNetwork

    env = GymEnv("ant-medium-v0")

    a = state(env) |> model

    @test action_space(env) |> size == size(a[1])
end

@testset "d4rl_policy_evaluate" begin
    plt = deep_ope_d4rl_evaluate("halfcheetah", "online", 10; num_evaluations=100)
    @test typeof(plt) <: UnicodePlots.Plot
end


text:
function to_simple_graph(network::PowerFlowNetwork)
    if !has_continuous_index(network)
        network = deepcopy(network)
        normalize_index!(network)
    end
    g = SimpleGraph(nbus(network))
    for (src, dst) in eachrow(network.branch[!, [:SRC, :DST]])
        add_edge!(g, src, dst)
    end
    return g
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
#@everywhere include("count_heads.jl")
#a = @spawn count_heads(100000000)
#b = @spawn count_heads(100000000)
#println(fetch(a) + fetch(b))

function par(n::Int64, x::SharedArray{Float64}, y::SharedArray{Float64})
    fe = @inbounds @parallel  for i in 1:n
       y[i] = 0.5*x[i] + y[i]
    end
    return fe
end

function ser(n::Int64, x::SharedArray{Float64}, y::SharedArray{Float64})
    @inbounds for i in 1:n
       y[i] = 0.5*x[i] + y[i]
    end
end

function init(n::Int64, x::SharedArray{Float64}, y::SharedArray{Float64})
    for i in 1:n
        x[i] = i
        y[i] = 1.0
    end
end
count = 1000
n = 1000000
y = SharedArray{Float64}(n)
x = SharedArray{Float64}(n)

println(procs(y))

init(n, x, y)

@time for j in 1:count
    @sync begin
        fe = par(n, x, y)
    end
    # for f in fe
    #     fetch(f)
    # end
end

println(y[n])

# for i = 1:n
#    c[i] = 0.
# end
init(n, x, y)

@time for j in 1:count
    ser(n, x, y)
end
println(y[n])


text:
"""
Create a new feature definition object to hold field definitions.

The OGRFeatureDefn maintains a reference count, but this starts at zero, and
should normally be incremented by the owner.
"""
unsafe_createfeaturedefn(name::AbstractString) =
    FeatureDefn(GDAL.fd_create(name))

"""
Increments the reference count by one.

The reference count is used keep track of the number of OGRFeature objects
referencing this definition.

### Returns
the updated reference count.
"""
reference(fd::FeatureDefn) = GDAL.reference(fd.ptr)

"Decrements the reference count by one, and returns the updated count."
dereference(fd::FeatureDefn) = GDAL.dereference(fd.ptr)

"Fetch the current reference count."
nreference(fd::FeatureDefn) = GDAL.getreferencecount(fd.ptr)

"Destroy a feature definition object and release all memory associated with it"
destroy(fd::FeatureDefn) = (GDAL.destroy(fd.ptr); fd.ptr = C_NULL)

"Drop a reference, and destroy if unreferenced."
release(fd::FeatureDefn) = GDAL.release(fd.ptr)

"Get name of the OGRFeatureDefn passed as an argument."
getname(fd::FeatureDefn) = GDAL.getname(fd.ptr)

"Fetch number of fields on the passed feature definition."
nfield(fd::FeatureDefn) = GDAL.getfieldcount(fd.ptr)

"""
Fetch field definition of the passed feature definition.

### Parameters
* `fd` the feature definition to get the field definition from.
* `i`           the field to fetch, between `0` and `nfield(fd)-1`.

### Returns
an handle to an internal field definition object or NULL if invalid index. This
object should not be modified or freed by the application.
"""
getfielddefn(fd::FeatureDefn, i::Integer) =
    FieldDefn(GDAL.getfielddefn(fd.ptr, i))

"""
Find field by name.

### Returns
the field index, or -1 if no match found.
"""
getfieldindex(fd::FeatureDefn, name::AbstractString) =
    GDAL.getfieldindex(fd.ptr, name)

"""
Add a new field definition to the passed feature definition.

To add a new field definition to a layer definition, do not use this function
directly, but use OGR_L_CreateField() instead.

This function should only be called while there are no OGRFeature objects in
existence based on this OGRFeatureDefn. The OGRFieldDefn passed in is copied,
and remains the responsibility of the caller.
"""
addfielddefn!(fd::FeatureDefn, fielddefn::FieldDefn) =
    (GDAL.addfielddefn(fd.ptr, fielddefn.ptr); fd)

"""
Delete an existing field definition.

To delete an existing field definition from a layer definition, do not use this
function directly, but use `OGR_L_DeleteField()` instead.

This method should only be called while there are no OGRFeature objects in
existence based on this OGRFeatureDefn.
"""
function deletefielddefn!(fd::FeatureDefn, i::Integer)
    result = GDAL.deletefielddefn(fd.ptr, i)
    @ogrerr result "Failed to delete field $i in the feature definition"
    fd
end

"""
Reorder the field definitions in the array of the feature definition.

To reorder the field definitions in a layer definition, do not use this function
directly, but use `OGR_L_ReorderFields()` instead.

This method should only be called while there are no OGRFeature objects in
existence based on this OGRFeatureDefn.

### Parameters
* **fd**: handle to the feature definition.
* **indices**: an array of `GetFieldCount()` elements which is a permutation of
    `[0, GetFieldCount()-1]`. `indices` is such that, for each field definition
    at position `i` after reordering, its position before reordering was
    `indices[i]`.
"""
function reorderfielddefns!(fd::FeatureDefn, indices::Vector{Cint})
    result = GDAL.reorderfielddefns(fd.ptr, indices)
    @ogrerr result "Failed to reorder $indices in the feature definition"
    fd
end

"""
Fetch the geometry base type of the passed feature definition.

For layers without any geometry field, this method returns `wkbNone`.

This returns the same result as `OGR_FD_GetGeomType(OGR_L_GetLayerDefn(hLayer))`
but for a few drivers, calling `OGR_L_GetGeomType()` directly can avoid lengthy
layer definition initialization.

For layers with multiple geometry fields, this method only returns the geometry
type of the first geometry column. For other columns, use
    `OGR_GFld_GetType(OGR_FD_GetGeomFieldDefn(OGR_L_GetLayerDefn(hLayer), i))`.
"""
getgeomtype(fd::FeatureDefn) = GDAL.getgeomtype(fd.ptr)

"""
Assign the base geometry type for the passed layer (same as the fd).

All geometry objects using this type must be of the defined type or a derived
type. The default upon creation is `wkbUnknown` which allows for any geometry
type. The geometry type should generally not be changed after any OGRFeatures
have been created against this definition.
"""
function setgeomtype!(fd::FeatureDefn, etype::OGRwkbGeometryType)
    GDAL.setgeomtype(fd.ptr, etype)
    fd
end

"Determine whether the geometry can be omitted when fetching features."
isgeomignored(fd::FeatureDefn) = Bool(GDAL.isgeometryignored(fd.ptr))

"Set whether the geometry can be omitted when fetching features."
setgeomignored!(fd::FeatureDefn, ignore::Bool) =
    (GDAL.setgeometryignored(fd.ptr, ignore); fd)

"Determine whether the style can be omitted when fetching features."
isstyleignored(fd::FeatureDefn) = Bool(GDAL.isstyleignored(fd.ptr))

"Set whether the style can be omitted when fetching features."
setstyleignored!(fd::FeatureDefn, ignore::Bool) =
    (GDAL.setstyleignored(fd.ptr, ignore); fd)

"Fetch number of geometry fields on the passed feature definition."
ngeomfield(fd::FeatureDefn) = GDAL.getgeomfieldcount(fd.ptr)

"""
Fetch geometry field definition of the passed feature definition.

### Parameters
* `i`  geometry field to fetch, between `0` (default) and `ngeomfield(fd)-1`.

### Returns
an internal field definition object or `NULL` if invalid index. This object
should not be modified or freed by the application.
"""
getgeomfielddefn(fd::FeatureDefn, i::Integer = 0) =
    GeomFieldDefn(GDAL.getgeomfielddefn(fd.ptr, i))

"""
Find geometry field by name.

The geometry field index of the first geometry field matching the passed field
name (case insensitively) is returned.

### Returns
the geometry field index, or -1 if no match found.
"""
getgeomfieldindex(fd::FeatureDefn, name::AbstractString = "") =
    GDAL.getgeomfieldindex(fd.ptr, name)

"""
Add a new field definition to the passed feature definition.

To add a new geometry field definition to a layer definition, do not use this
function directly, but use OGRLayer::CreateGeomField() instead.

This method does an internal copy of the passed geometry field definition,
unless bCopy is set to FALSE (in which case it takes ownership of the field
definition.

This method should only be called while there are no OGRFeature objects in
existence based on this OGRFeatureDefn. The OGRGeomFieldDefn passed in is
copied, and remains the responsibility of the caller.
"""
addgeomfielddefn!(fd::FeatureDefn, geomfielddefn::GeomFieldDefn) =
    (GDAL.addgeomfielddefn(fd.ptr, geomfielddefn.ptr); fd)

"""
Delete an existing geometry field definition.

To delete an existing field definition from a layer definition, do not use this
function directly, but use OGRLayer::DeleteGeomField() instead.

This method should only be called while there are no OGRFeature objects in
existence based on this OGRFeatureDefn.
"""
function deletegeomfielddefn!(fd::FeatureDefn, i::Integer)
    result = GDAL.deletegeomfielddefn(fd.ptr, i)
    @ogrerr result "Failed to delete geom field $i in the feature definition"
    fd
end

"Test if the feature definition is identical to the other one."
issame(fd1::FeatureDefn, fd2::FeatureDefn) = Bool(GDAL.issame(fd1.ptr, fd2.ptr))

"""Returns the new feature object with null fields and no geometry

Note that the OGRFeature will increment the reference count of it's defining
OGRFeatureDefn. Destruction of the OGRFeatureDefn before destruction of all
OGRFeatures that depend on it is likely to result in a crash.

Starting with GDAL 2.1, returns NULL in case out of memory situation.
"""
unsafe_createfeature(fd::FeatureDefn) =
    Feature(GDAL.f_create(fd.ptr))

"Fetch feature definition."
getfeaturedefn(feature::Feature) = FeatureDefn(GDAL.getdefnref(feature.ptr))


text:
using .MDP
using Printf

# Only needed for rle function
using StatsBase

L = 100
alpha = [0.75 0.75]

P = [ alpha[1]		1-alpha[1]
	  1 - alpha[2]  alpha[2]  ]

# Initialize xi (could be optimized)
xi = [ P^n for n = 1:L+1]

r = 1.0

for c = [0.1 0.2 0.3 0.4 0.5]
    for p = [0.1 0.2 0.3 0.4]
        q = 1 - p

        function bellmanUpdate(v; discount=1.0)
            # Assume v is 2 x 2 x (L+1)
            v_next = zeros(Float64, size(v))
            g_next = zeros(Int32,   size(v))

            # v_next(n, s, l)
            # n = 1 == 0 packets in buffer
            # n = 2 == 1 packet in buffer

            # s is the last channel state observed
            # s = 1 == Idle
            # s = 2 == Busy

            # ell is the time when the last observation was made
            for s = 1:2
                for ell = 1:L+1
                    next_ell = min(ell+1, L+1)
                    v_next[1, s, ell] = discount * (q * v[1, s, next_ell] + p * v[2, s, next_ell])
                    g_next[1, s, ell] = 1

                    next = [ discount * v[2, s, next_ell],
                             xi[ell][s,1]*r - c + discount * ( xi[ell][s,1] * ( q * v[1, 1, 1] + p * v[2, 1, 1])
                                                             + xi[ell][s,2] * v[2, 2, 1]) ]
                    idx = g_next[2, s, ell] = indmax(next)
                    v_next[2, s, ell] = next[ idx ]
                end
            end

            return (v_next, g_next)
        end

        model = DynamicModel(bellmanUpdate; objective=:Max)
        v_initial = zeros(2,2,L+1)

        (v,g) = valueIteration(model, v_initial; discount=0.9)

        g1 = vec(g[2,1,:])
        g2 = vec(g[2,2,:])

        a1,b1 = rle(g1[1:end-1])
        a2,b2 = rle(g2[1:end-1])

        k1 = (a1[1] == 1.0) ? (b1[1]+1) : 1
        k2 = (a2[1] == 1.0) ? (b2[1]+1) : 1

        @printf("c=%0.2f\tp=%0.2f\t(%d,%d)\n", c,p, k1,k2)
  end
end

# v1 = vec(v[2,1,:])
# v2 = vec(v[2,2,:])
#
# println(v1)
# println(v2)


text:
using VariationalOnlineStats
using Test, Statistics

for n = 1:2
    for i = 1:100
        x = rand(n)
        @test Float64(MedianOnline(x)) ≈ median(x)
    end
end

@test Float64(MedianOnline(Float64[1, 3, 2])) ≈ 2.0
@test Float64(MedianOnline(Float64[3, 1, 2])) ≈ 2.0
@test abs(Float64(MedianOnline(Float64[1, 2, 3])) - 2) < 0.1
@test abs(Float64(MedianOnline(Float64[3, 2, 1])) - 2) < 0.1

v = [1, 11, 10, 2, 3, 9, 8, 4, 5, 7, 6]
@test abs(Float64(MedianOnline{Float64}(v)) - 6) < 0.1
@test Float64(MedianOnline{Float64}(1:11)) > 4.5
@test Float64(MedianOnline{Float64}(11:-1:1)) < 7.5


text:
using SpikingNN
using Plots

# SRM0 params
η₀ = 5.0
τᵣ = 1.0

# Input spike train params
rate = 0.01
T = 20
∂t = 0.01
n = convert(Int, ceil(T / ∂t))

srm = Neuron(QueuedSynapse(Synapse.Delta()), SRM0(η₀, τᵣ), Threshold.Poisson(5.0, 0.5, 0.1))
input = ConstantRate(rate)
spikes = excite!(srm, input, n)

# callback to record voltages
voltages = Float64[]
record = function ()
    push!(voltages, getvoltage(srm))
end

# simulate
@time output = simulate!(srm, n; dt = ∂t, cb = record, dense = true)

# plot raster plot
raster_plot = rasterplot(∂t .* spikes, ∂t .* output, label = ["Input", "Output"], xlabel = "Time (sec)",
                title = "Raster Plot (δ response)")
xlims!(0, T)

# plot dense voltage recording
plot(∂t .* collect(1:n), voltages,
    title = "SRM Membrane Potential with Varying Presynaptic Responses", xlabel = "Time (sec)", ylabel = "Potential (V)", label = "δ response")

# resimulate using presynaptic response
srm = Neuron(QueuedSynapse(Synapse.Alpha()), SRM0(η₀, τᵣ), Threshold.Poisson(5.0, 0.5, 0.1))
voltages = Float64[]
excite!(srm, spikes)
@time simulate!(srm, n; dt = ∂t, cb = record, dense = true)

# plot voltages with response function
voltage_plot = plot!(∂t .* collect(1:n), voltages, label = "α response")
xlims!(0, T)

plot(raster_plot, voltage_plot, layout = grid(2, 1), xticks = 0:T)


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

using  Base.MultiplicativeInverses: SignedMultiplicativeInverse

immutable ReshapedArray{T,N,P<:AbstractArray,MI<:Tuple{Vararg{SignedMultiplicativeInverse{Int}}}} <: AbstractArray{T,N}
    parent::P
    dims::NTuple{N,Int}
    mi::MI
end
ReshapedArray{T,N}(parent::AbstractArray{T}, dims::NTuple{N,Int}, mi) = ReshapedArray{T,N,typeof(parent),typeof(mi)}(parent, dims, mi)

# LinearFast ReshapedArray
typealias ReshapedArrayLF{T,N,P<:AbstractArray} ReshapedArray{T,N,P,Tuple{}}

# Fast iteration on ReshapedArrays: use the parent iterator
immutable ReshapedArrayIterator{I,M}
    iter::I
    mi::NTuple{M,SignedMultiplicativeInverse{Int}}
end
ReshapedArrayIterator(A::ReshapedArray) = _rs_iterator(parent(A), A.mi)
function _rs_iterator{M}(P, mi::NTuple{M})
    iter = eachindex(P)
    ReshapedArrayIterator{typeof(iter),M}(iter, mi)
end

immutable ReshapedIndex{T}
    parentindex::T
end

# eachindex(A::ReshapedArray) = ReshapedArrayIterator(A)  # TODO: uncomment this line
start(R::ReshapedArrayIterator) = start(R.iter)
@inline done(R::ReshapedArrayIterator, i) = done(R.iter, i)
@inline function next(R::ReshapedArrayIterator, i)
    item, inext = next(R.iter, i)
    ReshapedIndex(item), inext
end
length(R::ReshapedArrayIterator) = length(R.iter)

reshape(parent::AbstractArray, dims::IntOrInd...) = reshape(parent, dims)
reshape(parent::AbstractArray, shp::NeedsShaping) = reshape(parent, to_shape(shp))
reshape(parent::AbstractArray, dims::Dims)        = _reshape(parent, dims)

reshape{T,N}(parent::AbstractArray{T,N}, ndims::Type{Val{N}}) = parent
function reshape{T,AN,N}(parent::AbstractArray{T,AN}, ndims::Type{Val{N}})
    reshape(parent, rdims((), indices(parent), Val{N}))
end
# Move elements from inds to out until out reaches the desired
# dimensionality N, either filling with OneTo(1) or collapsing the
# product of trailing dims into the last element
@pure rdims{N}(out::NTuple{N}, inds::Tuple{}, ::Type{Val{N}}) = out
@pure function rdims{N}(out::NTuple{N}, inds::Tuple{Any, Vararg{Any}}, ::Type{Val{N}})
    l = length(last(out)) * prod(map(length, inds))
    (front(out)..., OneTo(l))
end
@pure rdims{N}(out::Tuple, inds::Tuple{}, ::Type{Val{N}}) = rdims((out..., OneTo(1)), (), Val{N})
@pure rdims{N}(out::Tuple, inds::Tuple{Any, Vararg{Any}}, ::Type{Val{N}}) = rdims((out..., first(inds)), tail(inds), Val{N})

# _reshape on Array returns an Array
_reshape(parent::Vector, dims::Dims{1}) = parent
_reshape(parent::Array, dims::Dims{1}) = reshape(parent, dims)
_reshape(parent::Array, dims::Dims) = reshape(parent, dims)

# When reshaping Vector->Vector, don't wrap with a ReshapedArray
function _reshape(v::AbstractVector, dims::Dims{1})
    len = dims[1]
    len == length(v) || throw(DimensionMismatch("parent has $(length(v)) elements, which is incompatible with length $len"))
    v
end
# General reshape
function _reshape(parent::AbstractArray, dims::Dims)
    n = _length(parent)
    prod(dims) == n || throw(DimensionMismatch("parent has $n elements, which is incompatible with size $dims"))
    __reshape((parent, linearindexing(parent)), dims)
end

# Reshaping a ReshapedArray
_reshape{T}(v::ReshapedArray{T,1}, dims::Dims{1}) = _reshape(v.parent, dims)
_reshape(R::ReshapedArray, dims::Dims) = _reshape(R.parent, dims)

function __reshape(p::Tuple{AbstractArray,LinearSlow}, dims::Dims)
    parent = p[1]
    strds = front(size_strides(parent))
    strds1 = map(s->max(1,s), strds)  # for resizing empty arrays
    mi = map(SignedMultiplicativeInverse, strds1)
    ReshapedArray(parent, dims, reverse(mi))
end

function __reshape(p::Tuple{AbstractArray,LinearFast}, dims::Dims)
    parent = p[1]
    ReshapedArray(parent, dims, ())
end

@inline size_strides(A::AbstractArray) = tail(size_strides((1,), size(A)...))
size_strides(out::Tuple) = out
@inline size_strides(out, s, sz...) = size_strides((out..., out[end]*s), sz...)

size(A::ReshapedArray) = A.dims
similar(A::ReshapedArray, eltype::Type, dims::Dims) = similar(parent(A), eltype, dims)
linearindexing{R<:ReshapedArrayLF}(::Type{R}) = LinearFast()
parent(A::ReshapedArray) = A.parent
parentindexes(A::ReshapedArray) = map(s->1:s, size(parent(A)))
reinterpret{T}(::Type{T}, A::ReshapedArray, dims::Dims) = reinterpret(T, parent(A), dims)

@inline ind2sub_rs(::Tuple{}, i::Int) = i
@inline ind2sub_rs(strds, i) = ind2sub_rs((), strds, i-1)
@inline ind2sub_rs(out, ::Tuple{}, ind) = (ind+1, out...)
@inline function ind2sub_rs(out, strds, ind)
    d, r = divrem(ind, strds[1])
    ind2sub_rs((d+1, out...), tail(strds), r)
end

@inline function getindex(A::ReshapedArrayLF, index::Int)
    @boundscheck checkbounds(A, index)
    @inbounds ret = parent(A)[index]
    ret
end
@inline function getindex(A::ReshapedArray, indexes::Int...)
    @boundscheck checkbounds(A, indexes...)
    _unsafe_getindex(A, indexes...)
end
@inline function getindex(A::ReshapedArray, index::ReshapedIndex)
    @boundscheck checkbounds(parent(A), index.parentindex)
    @inbounds ret = parent(A)[index.parentindex]
    ret
end

@inline function _unsafe_getindex(A::ReshapedArray, indexes::Int...)
    @inbounds ret = parent(A)[ind2sub_rs(A.mi, sub2ind(size(A), indexes...))...]
    ret
end
@inline function _unsafe_getindex(A::ReshapedArrayLF, indexes::Int...)
    @inbounds ret = parent(A)[sub2ind(size(A), indexes...)]
    ret
end

@inline function setindex!(A::ReshapedArrayLF, val, index::Int)
    @boundscheck checkbounds(A, index)
    @inbounds parent(A)[index] = val
    val
end
@inline function setindex!(A::ReshapedArray, val, indexes::Int...)
    @boundscheck checkbounds(A, indexes...)
    _unsafe_setindex!(A, val, indexes...)
end
@inline function setindex!(A::ReshapedArray, val, index::ReshapedIndex)
    @boundscheck checkbounds(parent(A), index.parentindex)
    @inbounds parent(A)[index.parentindex] = val
    val
end

@inline function _unsafe_setindex!(A::ReshapedArray, val, indexes::Int...)
    @inbounds parent(A)[ind2sub_rs(A.mi, sub2ind(size(A), indexes...))...] = val
    val
end
@inline function _unsafe_setindex!(A::ReshapedArrayLF, val, indexes::Int...)
    @inbounds parent(A)[sub2ind(size(A), indexes...)] = val
    val
end

# helpful error message for a common failure case
typealias ReshapedRange{T,N,A<:Range} ReshapedArray{T,N,A,Tuple{}}
setindex!(A::ReshapedRange, val, index::Int) = _rs_setindex!_err()
setindex!(A::ReshapedRange, val, indexes::Int...) = _rs_setindex!_err()
setindex!(A::ReshapedRange, val, index::ReshapedIndex) = _rs_setindex!_err()

_rs_setindex!_err() = error("indexed assignment fails for a reshaped range; consider calling collect")

unsafe_convert{T}(::Type{Ptr{T}}, a::ReshapedArray{T}) = unsafe_convert(Ptr{T}, parent(a))

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
# https://docs.julialang.org/en/stable/manual/getting-started/
println(PROGRAM_FILE)

for x in ARGS
   println(x)
end

println("julia -p --procs N|auto")

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
#base.jl: ExploreCMEO base structures and functionnality
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#
const LIST_GRADEID = ["1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12"]
#const MAXCOL_GRADEID = 7

#==Types
===============================================================================#
mutable struct ExploreSelection
	subject::String
	#zero-based indices:
	grade_idx::Int
	domain_idx::Int
	attente_idx::Int
	content_idx::Int
end
ExploreSelection(subject::String="", grade=0, domain=0, attente=0, content=0) =
	ExploreSelection(subject, grade, domain, attente, content)

#Different views made up of one or more fields:
abstract type AbstractFieldView; end

struct FViewNone<:AbstractFieldView; end
struct FViewSrcDoc<:AbstractFieldView; end
struct FViewDomain<:AbstractFieldView; end
struct FViewAttente<:AbstractFieldView; end
struct FViewAttente_LongDescr<:AbstractFieldView; end
struct FViewContent<:AbstractFieldView; end

struct FieldViewOpts
	include_id::Bool #Ex: A1.1
	include_shortdescr::Bool
	include_longdescr::Bool
end
FieldViewOpts() = FieldViewOpts(true, true, true)

#==Validation/Helper functions
===============================================================================#
#AbstractString: might have used strip() on string.
isempty_subject(subject::AbstractString) = (strip(subject) == "")
isvalid_subject(subject::AbstractString) = !isempty_subject(subject)

#==Create identifiers for different curriculum elements (ex: A1.1)
===============================================================================#
function _valididx(args...)
	for v in args
		if v < 1; return false; end
	end
	return true
end
function getid_grade(grade_idx::Int)
	if grade_idx in (1:length(LIST_GRADEID))
		return LIST_GRADEID[grade_idx]
	else
		return ""
	end
end
_getid_domain(domain_idx::Int) = string('A'+(domain_idx-1))
function _getid_attente(domain_idx::Int, attente_idx::Int)
	prefix = getid_domain(domain_idx)
	return "$prefix$attente_idx"
end
function _getid_content(domain_idx::Int, attente_idx::Int, content_idx::Int)
	prefix = _getid_attente(domain_idx, attente_idx)
	return "$prefix.$content_idx"
end
getid_domain(domain_idx::Int) = (_valididx(domain_idx) ? _getid_domain(domain_idx) : "")
getid_attente(domain_idx::Int, attente_idx::Int) =
	(_valididx(domain_idx, attente_idx) ? _getid_attente(domain_idx, attente_idx) : "")
getid_content(domain_idx::Int, attente_idx::Int, content_idx::Int) =
	(_valididx(domain_idx, attente_idx, content_idx) ? _getid_content(domain_idx, attente_idx, content_idx) : "")

#==Get identifiers for different field views
===============================================================================#
getlabel(::AbstractFieldView) = ""
getlabel(::FViewSrcDoc) = "Document source"
getlabel(::FViewDomain) = "Domaine"
getlabel(::FViewAttente) = "Attente"
getlabel(::FViewAttente_LongDescr) = getlabel(FViewAttente())
getlabel(::FViewContent) = "Contenu d'Apprentissage"

getid(::AbstractFieldView, sel::ExploreSelection) = ""
getid(fv::FViewDomain, sel::ExploreSelection) =
	getid_domain(sel.domain_idx)
getid(fv::FViewAttente, sel::ExploreSelection) =
	getid_attente(sel.domain_idx, sel.attente_idx)
getid(fv::FViewAttente_LongDescr, sel::ExploreSelection) =
	getid(FViewAttente(), sel)
getid(fv::FViewContent, sel::ExploreSelection) =
	getid_content(sel.domain_idx, sel.attente_idx, sel.content_idx)

getlabelledid(fv::AbstractFieldView, sel::ExploreSelection) =
	string(getlabel(fv), " ", getid(fv, sel))

#==Creating ExploreCMEO launcher
===============================================================================#
function get_launcher_cmd()
	launcherpath = abspath(joinpath(@__DIR__, "..", "lancer_explorecmeo.jl"))
	juliacmd = joinpath(Sys.BINDIR::String, Base.julia_exename())
	return string("\"", juliacmd, "\" -L \"", launcherpath, "\"")
end
function show_launcher_info()
	cmd = get_launcher_cmd()
	@info("Cible:\n$cmd")
	return
end

#==Accessing database file
===============================================================================#
function download_database(destfilepath::String)
	@info("Téléchargement de données:\n$URL_PUBLISHEDDATA")
	run(BinDeps.download_cmd(URL_PUBLISHEDDATA, destfilepath))
	return
end
function open_database(filepath::String, mode::String; log::Bool=true)
	filepath = abspath(filepath)
	db = HDF5.h5open(filepath, mode)
	@info("Données source:\n$(db.filename)")
	return db
end

#Last line


text:
using Plots
using Random
using DifferentialEquations
using LinearAlgebra
using DynamicalSystems
using Parameters
using BasicReservoirComputing

include("utilities.jl")

function test_rc(system::DynamicalSystem, params::NamedTuple; Δt = 0.01)
    D=length(system.u0)
    utrain, utest_arr = make_data(system; train_time = 100, test_time = 25, n_test = 100)

    @unpack N, SR, ρA, α, σ, β, σb, nspin = params

    # Build and Train RC
    reservoir = rc(D, N, ρA, Δt, SR, α, σ, σb, β)
    train_RC(reservoir, utrain, spinup=nspin)

    # Forecast RC
    utrue = utest_arr[1][:, nspin:end]
    nsteps = size(utrue)[2]
    upred = forecast_RC(reservoir, nsteps, uspin = utest_arr[1][:, 1:nspin])
    t = collect(0:nsteps-1)*Δt

    p1 = plot_prediction(t, upred, utrue)

    # Test on many points
    vpt = test_RC(reservoir, utest_arr; ϵ=0.3, spinup=nspin)
    p2 = histogram(vpt, legend=false,
                   overwrite_figure=false)
    xlabel!(p2, "Valid Prediction Time")
    pfinal = plot(p1, p2, layout = (1, 2), show=true)
    display(pfinal)
    return reservoir
end

params = (N=200, σ=0.084, α=0.6, SR=0.8, ρA=0.02, β=8.5e-8, σb=1.6, nspin=200)
system = Systems.lorenz()
reservoir = test_rc(system, params);

max_stars_repo_name: japlatt/BasicReservoirComputing

max_issues_repo_name: japlatt/BasicReservoirComputing

max_forks_repo_name: japlatt/BasicReservoirComputing


text:
"""
    barchart(values;
            boundingbox = BoundingBox(O + (-250, -120), O + (250, 120)),
            bargap=10,
            margin = 5,
            border=false,
            labels=false,
            labelfunction = (values, i, lowpos, highpos, barwidth, scaledvalue) -> begin
                    label(string(values[i]), :n, highpos, offset=10)
              end,
            barfunction =  (values, i, lowpos, highpos, barwidth, scaledvalue) -> begin
                @layer begin
                    setline(barwidth)
                    line(lowpos, highpos, :stroke)
                end
              end)

Draw a barchart where each bar is the height of a value in the `values` array. The bars
will be scaled to fit in a bounding box.

Text labels are drawn if the keyword `labels=true`.

# Extended help

The function returns a vector of points; each is the bottom center of a bar.

Draw a Fibonacci sequence as a barchart:

```
fib(n) = n > 2 ? fib(n - 1) + fib(n - 2) : 1
fibs = fib.(1:15)
@draw begin
    fontsize(12)
    barchart(fibs, labels=true)
end
```
To control the drawing of the text and bars, define functions that process the
end points:

`mybarfunction(values, i, lowpos, highpos, barwidth, scaledvalue)`

`mylabelfunction(values, i, lowpos, highpos, barwidth, scaledvalue)`

and pass them like this:

```julia
barchart(vals, barfunction=mybarfunction)
barchart(vals, labelfunction=mylabelfunction)
```

```
function myprologfunction(values, basepoint, minbarrange, maxbarrange, barchartheight)
    @layer begin
        setline(0.2)
        for i in 0:10:maximum(values)
            rule(boxbottomcenter(basepoint) + (0, -(rescale(i, minbarrange, maxbarrange) * barchartheight)))
        end
    end
end
```
"""
function barchart(values;
    boundingbox = BoundingBox(O + (-250, -120), O + (250, 120)),
    bargap=10,
    margin = 5,
    border=false,
    labels=false,
    labelfunction = (values, i, lowpos, highpos, barwidth, scaledvalue) -> begin
        label(string(values[i]), :n, highpos, offset=10)
    end,
    barfunction =  (values, i, lowpos, highpos, barwidth, scaledvalue) -> begin
        @layer begin
            setline(barwidth)
            line(lowpos, highpos, :stroke)
        end
    end,
    prologfunction = (values, basepoint, minbarrange, maxbarrange, barchartheight) -> ()
    )
    minvalue, maxvalue = extrema(values)
    barchartwidth  = boxwidth(boundingbox)  - 2bargap - 2margin
    barchartheight = boxheight(boundingbox) - 2margin
    barwidth = (barchartwidth - 2bargap)/length(values)
    if barwidth < 0.1
        throw(error("barchart() - bars are too small (< 0.1) at $(barwidth)"))
    end
    # if all bars are equal height, this will force a range
    minbarrange = minvalue - abs(minvalue)
    maxbarrange = maxvalue + abs(maxvalue)
    basepoint = boundingbox - (0, margin)
    hpositions = between.(
        boxbottomleft(basepoint),
        boxbottomright(basepoint),
        # skip first and last, then take every other one, which is at halfway
        range(0.0, stop=1.0, length=2length(values) + 1))[2:2:end-1]
    @layer begin
        if border
            box(boundingbox, :stroke)
        end
        prologfunction(values, basepoint, minbarrange, maxbarrange, barchartheight)
        for i in 1:length(values)
            scaledvalue = rescale(values[i], minbarrange, maxbarrange) * barchartheight
            lowposition = hpositions[i]
            highposition = lowposition - (0, scaledvalue) # -y coord
            barfunction(values, i, lowposition, highposition, barwidth, scaledvalue)
            labels && labelfunction(values, i, lowposition, highposition, barwidth, scaledvalue)
        end
    end
    return (positions = hpositions)
end


text:
"""
NueralNetwork Functions module

MIT License
Copyright (c) 2019 amenoyoya
"""

"""
活性化関数: 入力値に何らかの処理を施し出力する
    入力次元 = 出力次元となる
"""

# 恒等関数
identity_function(x::Float64)::Float64 = x

# ステップ関数
step_function(x::Float64)::Float64 = x > 0 ? x : 0

# シグモイド関数
sigmoid(x::Array{Float64,1})::Array{Float64,1} = 1 ./ (1 .+ exp.(-x))

# シグモイド関数（バッチ対応）
sigmoid(x::Array{Float64,2})::Array{Float64,2} = hcat([sigmoid(x[row, :]) for row in 1:size(x, 1)]...)'

# ReLU関数
relu(x::Array{Float64,1})::Array{Float64,1} = maximum.(0.0, x)

# ReLU関数（バッチ対応）
relu(x::Array{Float64,2})::Array{Float64,2} = hcat([relu(x[row, :]) for row in 1:size(x, 1)]...)'

# ソフトマックス関数
softmax(x::Array{Float64,1})::Array{Float64,1} = begin
    x = x .- maximum(x) # オーバーフロー対策
    exp.(x) ./ sum(exp.(x))
end

# ソフトマックス関数（バッチ対応）
softmax(x::Array{Float64,2})::Array{Float64,2} = hcat([softmax(x[row, :]) for row in 1:size(x, 1)]...)'

"""
損失関数: 得た予測値と正解値の誤差をFloat64型で出力
    予測次元 = 正解次元 ≠ 出力次元 となる
"""

# 二乗和誤差
## 予測値, 正解値 -> 誤差
mean_squared_error(y::Array{Float64,2}, t::Array{Float64,2})::Float64 = 0.5 * sum((y .- t).^2)

# 交差エントロピー誤差
## 予測値, 正解値 -> 誤差
cross_entropy_error(y::Array{Float64,1}, t::Array{Float64,1})::Float64 = begin
    # 教師データ: one-hot-vectorを想定 => 正解ラベルのindexのみ計算
    val, index = findmax(t) # 正解ラベル = 最大値(=1.0)のindex
    # ln(0) = -Inf が発生するのを防ぐため、予測値に微小な値（10^-7）を加算して計算する
    -log(y[index] .+ 1e-7)
end

# 交差エントロピー誤差（バッチ対応）
## 予測値, 正解値 -> 誤差
cross_entropy_error(y::Array{Float64,2}, t::Array{Float64,2})::Float64 = begin
    batch_size = size(y, 1)
    sum([cross_entropy_error(y[row, :], t[row, :]) for row in 1:batch_size]) / batch_size
end

max_forks_repo_name: amenoyoya/scala-ml


text:
EIN = (
    ("transformer.X.ix", [0.0, 0.0]),
    ("synUnit.SM.Psi", [1.0, 0.0]),
    ("synUnit.SM.Psi1", [1.0, 0.0]),
    ("synUnit.SM.phiKR", 0.0),
    ("synUnit.SM.Psifd", 1.20),
    ("synUnit.exciter.transducer.x", 1.0),
    ("synUnit.exciter.PI1.Integr.x", 0.0),
    ("synUnit.exciter.PI1f.x", 0.0),
    ("synUnit.exciter.PT1a.Integr.x", 0.0),
    (:VEXC, 1.01),
    (:LOAD, 0.5),
)

AUS = (
#     ("der(transformer.X.ix)", [0.0, 0.0]),
    ("der(synUnit.SM.Psi)", [0.0, 0.0]),
    ("der(synUnit.SM.Psi1)", [0.0, 0.0]),
    ("der(synUnit.SM.Psifd)", 0.0),
    ("der(synUnit.exciter.transducer.x)", 0.0),
    ("der(synUnit.exciter.PI1.Integr.x)", 0.0),
    ("der(synUnit.exciter.PI1f.x)", 0.0),
    ("der(synUnit.exciter.PT1a.Integr.x)", 0.0),
    ("der(J.w)", 0.0),
    ("meterGEN.powerRe", 0.8),
    ("meterGEN.powerIm", 0.3),
    ("transformer.X.derINI1", 0.0),
    ("transformer.X.derINI2", 0.0),

)

max_stars_repo_name: johhell/EPSphasor

max_issues_repo_name: johhell/EPSphasor

max_forks_repo_name: johhell/EPSphasor


text:
mutable struct Network
    input_size::Int64
    output_size::Int64
    hidden_neurons::Int64
    runtime::Int64

    weights::Array{Float32}
    voltage::Array{Float32}
    output::Array{Bool}

    fire_rec::Array{Int64}

    Pᵢ::Array{Float32}
    M::Array{Float32}

    construct
    activate

    function Network(;input_size::Int64, output_size::Int64, hidden_neurons::Int64, runtime::Int64)
        weights = rand(Float32, hidden_neurons+output_size, input_size+hidden_neurons)
        for i in 1:hidden_neurons
            weights[i, input_size+i] = 0.0f0
        end
        voltage = zeros(Float32, hidden_neurons+output_size, 1)
        output = zeros(Bool, input_size+hidden_neurons, 1)
        fire_rec = zeros(Int64, hidden_neurons+output_size, 1)

        Pᵢ = zeros(Float32, hidden_neurons+output_size, input_size+hidden_neurons)
        M = zeros(Float32, hidden_neurons+output_size, 1)

        new(input_size, output_size, hidden_neurons, runtime, weights, voltage, output, fire_rec, Pᵢ, M, construct, activate)
    end

    # n₂ => n₁
    function construct(self::Network; n₁::Int64, n₂::Int64)
        self.weights[n₁, n₂+self.input_size] = (-1)^rand(Bool)*rand(Float32)
        self.weights[n₂, n₁+self.input_size] = 0.0f0
    end

    # θ : threshold, r : refractory
    function activate(self::Network, spike_train::Array{Bool}; θ, r, η, γᵥ, γₚ, γₘ, A₊, A₋, w₊, w₋)
        self.fire_rec = fill(-r, self.hidden_neurons+self.output_size, 1)
        for t in 1:self.runtime
            # load input data
            self.output[1:self.input_size, 1] .= spike_train[:,t]

            # load result
            self.output[self.input_size+1:end, 1] .= self.voltage[1:self.hidden_neurons, 1] .>= θ

            # update the fire records
            self.fire_rec .= floor.(Int, self.fire_rec .* (t ./ self.fire_rec) .^ (self.voltage .>= θ))

            # reset voltage after firing
            self.voltage -= (self.voltage .>= θ) .* self.voltage

            # update pre-synaptic trace
            self.Pᵢ .*= γₚ
            for i in 1:self.hidden_neurons+self.output_size
                if t-self.fire_rec[i,1] > r
                    self.Pᵢ[i,:] .+= self.output[:,1] .* (self.weights[i,:] .!= 0.0f0) * A₊
                    self.weights[i,:] .+= η * (self.output[:,1] .* (self.weights[i,:] .!= 0.0f0))' * self.M[i,1] * (w₋ .- self.weights[i,:])
                end
            end

            self.voltage .*= γᵥ
            self.voltage .= self.weights * self.output .* ((t .- self.fire_rec) .> r) + self.voltage
            println(self.voltage[self.hidden_neurons+1:end])
            self.output[self.input_size+1:end, 1] .= self.voltage[1:self.hidden_neurons, 1] .>= θ
            for i in axes(self.weights, 1)
                self.weights[i,:] .+= η * (self.output[i,1] .* (self.weights[i,:] .!= 0.0f0))' .* self.Pᵢ[i,:] * (w₊ .- self.weights[i,:])
            end

            # update post-synaptic trace
            self.M .*= γₘ
            #self.M .-= (self.voltage[1:self.hidden_neurons, 1] .>= θ) * A₋
            self.M .-= (self.voltage .>= θ) * A₋
        end
    end
end

max_stars_repo_name: Louisenlund-plusMINT/NeuronSim

max_issues_repo_name: Louisenlund-plusMINT/NeuronSim

max_forks_repo_name: Louisenlund-plusMINT/NeuronSim


text:
using SequentialMonteCarlo
using RNGPool
import SMCExamples.StochVolatility.defaultSVModel

include("test.jl")

setRNGs(0)

model, guidedmodel, theta, ys = defaultSVModel(10)

numParticles = 1024*1024
numTrials = 2

testSMCParallel(model, numParticles, numTrials, false)
testSMCParallel(model, numParticles, numTrials, false, 0.5)
testSMCParallel(guidedmodel, numParticles, numTrials, false)
testSMCParallel(guidedmodel, numParticles, numTrials, false, 0.5)


text:
# Wrapper to CARET library.
module CaretWrapper

importall Orchestra.Types
importall Orchestra.Util

using PyCall
@pyimport rpy2.robjects as RO
@pyimport rpy2.robjects.packages as RP
@pyimport rpy2.robjects.numpy2ri as N2R
N2R.activate()
RP.importr("caret")

export CRTLearner,
       fit!,
       transform!

# Convert vector to R equivalent.
vector_to_r{T<:Int}(vector::Vector{T}) = RO.IntVector(vector)
vector_to_r{T<:Real}(vector::Vector{T}) = RO.FloatVector(vector)
vector_to_r{T<:Bool}(vector::Vector{T}) = RO.BoolVector(vector)
vector_to_r{T<:String}(vector::Vector{T}) = RO.StrVector(vector)
function vector_to_r(vector::Vector{Any})
  vec_eltype = infer_eltype(vector)
  if vec_eltype == Any
    error("Cannot handle R conversion for vector with differing element types.")
  end

  return vector_to_r(convert(Vector{vec_eltype}, vector))
end
vector_to_r(vector::Vector) = error(
  "Cannot handle R conversion for $(typeof(vector))."
)

# Builds R dataframe out of dataset.
# Returns (dataframe, label_factor_levels).
function dataset_to_r_dataframe(
  instances::Matrix, labels=nothing)

  # Build dataframe
  df_dict = Dict()
  for col in 1:size(instances, 2)
    df_dict["X$col"] = vector_to_r(instances[:, col])
  end

  if labels != nothing
    r_labels = RO.FactorVector(labels)
    df_dict["Y"] = r_labels
    return (RO.DataFrame(df_dict), r_labels[:levels])
  else
    return (RO.DataFrame(df_dict), nothing)
  end
end

# CARET wrapper that provides access to all learners.
#
# Options for the specific CARET learner is to be passed
# in `options[:impl_options]` dictionary.
type CRTLearner <: Learner
  model
  options

  function CRTLearner(options=Dict())
    default_options = {
      # Output to train against
      # (:class).
      :output => :class,
      :learner => "svmLinear",
      :impl_options => Dict()
    }
    new(nothing, nested_dict_merge(default_options, options))
  end
end

function fit!(crtw::CRTLearner, instances::Matrix, labels::Vector)
  impl_options = crtw.options[:impl_options]
  crtw.model = Dict()
  crtw.model[:learner] = crtw.options[:learner]

  # Build R dataframe out of dataset
  (r_dataset_df, label_factors) = dataset_to_r_dataframe(instances, labels)

  # Assign label factors
  crtw.model[:label_factors] = collect(label_factors)

  # Train model
  caret_formula = RO.Formula("Y ~ .")
  r_fit_control = pycall(RO.r[:trainControl], PyObject,
    method = "none"
  )
  if isempty(impl_options)
    r_model = pycall(RO.r[:train], PyObject,
      caret_formula,
      method = crtw.model[:learner],
      data = r_dataset_df,
      trControl = r_fit_control,
      tuneLength = 1
    )
  else
    r_model = pycall(RO.r[:train], PyObject,
      caret_formula,
      method = crtw.model[:learner],
      data = r_dataset_df,
      trControl = r_fit_control,
      tuneGrid = RO.DataFrame(impl_options)
    )
  end
  crtw.model[:r_model] = r_model
end

function transform!(crtw::CRTLearner, instances::Matrix)
  (r_instance_df, _) = dataset_to_r_dataframe(instances)
  predictions = collect(pycall(RO.r[:predict], PyObject,
    crtw.model[:r_model],
    newdata = r_instance_df
  ))
  label_factors = crtw.model[:label_factors]
  predictions = map(x -> label_factors[x], predictions)
  return predictions
end

end # module


text:
module PkgDeps

using Pkg.Types: VersionNumber, VersionRange
using REPL
using TOML: parsefile
using UUIDs

export PkgEntry, RegistryInstance
export NoUUIDMatch, PackageNotInRegistry
export users, reachable_registries

include("pkg_entry.jl")
include("registry_instance.jl")
include("exceptions.jl")
include("utilities.jl")

const GENERAL_REGISTRY = "General"

"""
    reachable_registries(registry_names::Array)
    reachable_registries(registry_name::String; depots::Union{String, Vector{String}}=Base.DEPOT_PATH)
    reachable_registries(; depots::Union{String, Vector{String}}=Base.DEPOT_PATH)

Get an array of found registries.

# Arguments
- `registry_names::Array`: List of registries to retrieve

# Keywords
- `depots::Union{String, Vector{String}}=Base.DEPOT_PATH`: Depots to look in for registries

# Returns
- `Array{RegistryInstance}`: List of all found registries
"""
function reachable_registries(
    registry_names::Array;
    depots::Union{String, Vector{String}}=Base.DEPOT_PATH,
    kwargs...
)
    registries = RegistryInstance[]

    if depots isa String
        depots = [depots]
    end

    for d in depots
        isdir(d) || continue
        reg_dir = joinpath(d, "registries")
        isdir(reg_dir) || continue

        for name in readdir(reg_dir)
            if isempty(registry_names) || name in registry_names
                file = joinpath(reg_dir, name, "Registry.toml")
                isfile(file) || continue
                push!(registries, RegistryInstance(joinpath(reg_dir, name)))
            end
        end
    end

    return registries
end
reachable_registries(registry_name::String; depots::Union{String, Vector{String}}=Base.DEPOT_PATH, kwargs...) = reachable_registries([registry_name]; depots=depots, kwargs...)
reachable_registries(; depots::Union{String, Vector{String}}=Base.DEPOT_PATH, kwargs...) = reachable_registries([]; depots=depots, kwargs...)

"""
    users(uuid::UUID; registries::Array{RegistryInstance}=reachable_registries(), depots::Union{String, Vector{String}}=Base.DEPOT_PATH)
    users(pkg_name::String, pkg_registry_name::String=GENERAL_REGISTRY; kwargs...)
    users(pkg_name::String, pkg_registry::RegistryInstance; kwargs...))

Find the users of a given package.

# Arguments
- `uuid::UUID`: UUID of the package.
- `pkg_name::String`: Find users of this package.
- `pkg_registry_name::String="$GENERAL_REGISTRY"`: Name of registry where `pkg_name` is
  registered. This is used to look up the UUID of `pkg_name`.

# Keywords
- `registries::Array{RegistryInstance}=reachable_registries()`: Registries to search for users.

# Returns
- `Array{String}`: All packages which are dependent on the given package.
"""
function users(
    uuid::UUID;
    registries::Array{RegistryInstance}=reachable_registries(),
    depots::Union{String, Vector{String}}=Base.DEPOT_PATH
)
    pkg_name = _get_pkg_name(uuid, reachable_registries(; depots=depots))
    downstream_dependencies = String[]

    for rego in registries
        for (pkg, v) in rego.pkgs
            base_path = joinpath(v.registry_path, v.path)
            deps_file_path = joinpath(base_path, "Deps.toml")

            if isfile(deps_file_path)
                deps_content = parsefile(deps_file_path)
                dependency_versions = collect(keys(deps_content))
                latest_version = _get_latest_version(base_path)

                # Use the latest_version of pkg, and check to see if pkg_name is in its dependents
                for version_range in dependency_versions
                    if in(latest_version, VersionRange(version_range))
                        dependencies = collect(keys(deps_content[version_range]))

                        # Check if pkg_name is used in the latest version of pkg
                        if pkg_name in dependencies
                            push!(downstream_dependencies, pkg)
                        end
                    end
                end
            end
        end
    end

    return downstream_dependencies
end

function users(pkg_name::String, pkg_registry_name::String=GENERAL_REGISTRY; kwargs...)
    uuid = _get_pkg_uuid(pkg_name, pkg_registry_name)
    return users(uuid; kwargs...)
end

function users(pkg_name::String, pkg_registry::RegistryInstance; kwargs...)
    uuid = _get_pkg_uuid(pkg_name, pkg_registry)
    return users(uuid; kwargs...)
end

end  # module


text:
using Documenter

makedocs(
    sitename = "FastRationals.jl",
    authors = "Jeffrey Sarnoff",
    pages = Any[
        "The BPP formula for PI" => "bpp.md",
        "Finding the Range" => "findingtherange.md",
        "The Stateless Way" => "thestatelessway.md",
        "What cannot overflow?" => "mayoverflow.md",
        "Rational Magnitude In Action" => "rationalmagnitude.md",
        "A metaphor that illuminates" => "metaphoricalflashlight.md",
    ]
)

deploydocs(
    repo = "github.com/JeffreySarnoff/FastRationals.jl.git",
    target = "build"
)


text:
struct Permutation{T <: Integer}
    M::Matrix{T}
    function Permutation(M::AbstractMatrix{<:Integer})
        row = [i for i = 1:size(M, 1) if M[i,1] ≠ M[i,2]]
        dup = [i for i = 2:length(row) if M[row[i-1], 1] == M[row[i], 1]]
        deleteat!(row, dup)
        new{eltype(M)}(Array(M[row, :]))
    end
end

struct PermutationGroup <: AbstractFiniteGroup
    name::String
    multab::Matrix{Int64}
    inv::Vector{Int64}
    cls::Vector{Vector{Int64}}
    clsv::Vector{Int64}
    mult::Vector{Int64}
    operations::Vector{Permutation}
end

function Base.display(g::PointGroup)
    println("Permutation group : $(name(g))")
    println("Group order       : $(order(g))")
    println("Classes           : $(length(class(g)))")
end

operation(g::PermutationGroup) = g.operations
operation(g::PermutationGroup, i) = g.operations[i]

export permutationgroup
function permutationgroup(gens::AbstractVector{<:Permutation}; name::String="Unnamed Group")
    multab, eles = generate_group(gens)
    ginv = group_inverse(multab)
    cls, clsv = conjugate_class(multab, ginv)
    mult = length.(cls)
    PermutationGroup(name, multab, ginv, cls, clsv, mult, eles)
end

permutationgroup(n::Integer) = permutationgroup([cycles(1:2), cycles(1:n)], name="S$n")

function Base.string(p::Permutation)
    cyc = tocycles(p)
    str = "Cycle: "
    for c in cyc
        n = length(c)
        str *= "("
        for i = 1:n-1
            str *= "$(c[i]), "
        end
        str *= "$(c[n]))"
    end
    str
end
Base.repr(p::Permutation) = string(p)
Base.display(p::Permutation) = println(string(p))

Base.isequal(p1::Permutation, p2::Permutation) = isequal(p1.M, p2.M)
Base.isless(p1::Permutation, p2::Permutation) = isless(reshape(p1.M, :), reshape(p2.M, :))

export permutation
function permutation(M::AbstractMatrix{<:Integer})
    p = sortperm(view(M, :, 1))
    M = M[p, 1:2]
    row = [i for i = 1:size(M, 1) if M[i,1] ≠ M[i,2]]
    dup = [i for i = 2:length(row) if M[row[i-1], 1] == M[row[i], 1]]
    deleteat!(row, dup)
    Permutation(Array(M[row, :]))
end

export cycles
function cycles(cyc::Union{AbstractVector{<:Integer}, Tuple}...)
    n = sum(length(c) for c in cyc)
    iszero(n) && (return Permutation(zeros(Int64, 0, 2)))
    M = Matrix{Int64}(undef, n, 2)
    i = 1
    for c in cyc
        j = i + length(c) - 1
        M[i:j, 1] .= c
        M[j, 2] = c[1]
        M[i:j-1, 2] .= c[2:end]
        i = j + 1
    end
    permutation(M)
end

export permutation
function tocycles(p::Permutation)
    cyc = []
    ind = p.M[:, 1]
    n = length(ind)
    used = fill(false, n)
    i = 1
    while i <= n
        used[i] ? (i += 1; continue) : (used[i] = true)
        j = ind[i]
        c = [j]
        k = p(j)
        while k ≠ j
            used[binary_search(ind, k)[2]] = true
            push!(c, k)
            k = p(k)
        end
        isone(length(c)) || push!(cyc, c)
    end
    cyc
end

function (p::Permutation)(i::Integer)
    q, j = binary_search(view(p.M, :, 1), i)
    q ? p.M[j, 2] : i
end

function *(p1::Permutation, p2::Permutation)
    i = sort([view(p1.M, :, 1); view(p2.M, :, 1)]) |> delete_duplicate!
    j = p1.(p2.(i))
    Permutation([i j])
end

function generate_group(gens::AbstractVector)
    N = 0
    eles = gens
    while length(eles) > N
        N = length(eles)
        eles = vcat(eles, ([gen * ele for ele in eles] for gen in gens)...) |> sort
        delete_duplicate!(eles)
    end
    tb = Matrix{Int64}(undef, N, N)
    @threads for i = 1:N
        for j=1:N
            _, k = binary_search(eles, eles[i] * eles[j])
            tb[i, j] = k
        end
    end
    tb, eles
end

#-------------------------------------------------------------------------------
# Helpers
#-------------------------------------------------------------------------------
"""
    binary_search
"""
function binary_search(list::AbstractVector, i)
    l, r = 1, length(list)
    if iszero(r) || i < list[l]
        return false, 1
    elseif i > list[r]
        return false, r+1
    elseif isequal(i, list[l])
        return true, l
    elseif isequal(i, list[r])
        return true, r
    end
    while true
        c = (l + r) ÷ 2
        t = list[c]
        if i < t
            r = c
        elseif i > t
            l = c
        else
            return true, c
        end
        isone(r - l) && (return false, r)
    end
end

function delete_duplicate!(v)
    N = length(v)
    i = 2
    while i<=N
        if isequal(v[i], v[i-1])
            deleteat!(v, i)
            N -= 1
        else
            i += 1
        end
    end
    v
end


text:
import LatentGaussianMixtureModel
using Distributions, StatsBase
using Yeppp, Random
b=1
Ctrue=3

nF = 282
Random.seed!(100)
n_ij = round.(Int64, rand(Poisson(5), 282).+rand(Exponential(45), 282))
N = sum(n_ij)

groupindex = inverse_rle(1:nF, n_ij)
J=2  #42 #beta dimension
Random.seed!(100)
X = rand(Normal(), N, J)
betas_true=ones(J) #rand(Normal(0,1), J)
if Ctrue == 1
    mu_true = [-log(1/0.1 - 1)]
    wi_true = [1.0]
    sigmas_true = [1.2]
elseif Ctrue == 2
    mu_true = [-log(1/0.1 - 1) - 3.0, -log(1/0.1 - 1) + 1.0]
    wi_true = [.5, .5]
    sigmas_true = [1.2, .8]
elseif Ctrue == 3
    mu_true = [-log(1/0.1 - 1) - 6, -log(1/0.1 - 1) - 2, -log(1/0.1 - 1) + 2;]
    wi_true = [.3, .4, .3]
    sigmas_true = [1.2, .8, .9]
end

#Randomly data generation based on the setting on datagen.jl
Random.seed!(b * 100)
m = MixtureModel(map((u, v) -> Normal(u, v), mu_true, sigmas_true), wi_true)
gamma_true = rand(m, nF)

prob = exp.(gamma_true[groupindex] .+ X*betas_true)
prob= prob ./ (1 .+ prob)
Y = Array(Bool, N)
Random.seed!(b * 100)
for i in 1:N
    if rand(Binomial(1, prob[i])) == 1
        Y[i] = true
    else
        Y[i] = false
    end
end
#lr=loglikelihoodratio(X, Y, groupindex, Calternative, ntrials=25, debuginfo=debuginfo)

ncomponent1=Ctrue+1
C0 = ncomponent1 - 1
C1 = ncomponent1
nF = maximum(groupindex)
an1 = 0.0 # 1/nF
gamma_init, beta_init, sigmas_tmp = LatentGaussianMixtureModel.maxposterior(X, Y, groupindex)
wi_init, mu_init, sigmas_init, ml_tmp = LatentGaussianMixtureModel.gmm(gamma_init, C0, ones(C0)/C0, quantile(gamma_init, linspace(0, 1, C0+2)[2:end-1]), ones(C0), an=an1, maxiter=1)

 Random.seed!(100);@time wi_init, mu_init, sigmas_init, betas_init, ml_tmp, gamma_mat =LatentGaussianMixtureModel.latentgmm(X, Y, groupindex, C0, beta_init, wi_init, mu_init, sigmas_init, Mmax=5000, initial_iteration=10, maxiteration=30, an=an1, sn=std(gamma_init).*ones(C0))

trand=LatentGaussianMixtureModel.asymptoticdistribution(X, Y, groupindex, wi_init, mu_init, sigmas_init, betas_init)

vtau=[.5, .3, .1]
ngh=1000
Mctau=1000
gamma0 = vec(mean(gamma_mat, 2))
mingamma = minimum(gamma0)
maxgamma = maximum(gamma0)

lr = zeros(length(vtau), C0)
or = sortperm(mu_init)
wi0 = wi_init[or]
mu0 = mu_init[or]
sigmas0 = sigmas_init[or]
betas0 = betas_init
an = 0.0 # decidepenalty(wi0, mu0, sigmas0, nF)

N,J=size(X)
sample_gamma_mat = zeros(nF, Mctau)
sumlogmat = zeros(nF, ngh*ncomponent1)
llvec = zeros(N)
llvecnew = zeros(N)
xb = zeros(N)

whichtosplit=2
ind = [1:whichtosplit, whichtosplit:C0;]
mu_lb = [mingamma, (mu0[1:(C0-1)] .+ mu0[2:C0])./2;]
mu_ub = [(mu0[1:(C0-1)] .+ mu0[2:C0])./2, maxgamma;]
mu_lb = mu_lb[ind]
mu_ub = mu_ub[ind]
sigmas_lb = 0.25 .* sigmas0[ind]
sigmas_ub = 2 .* sigmas0[ind]

tau=.1

ind = [1:whichtosplit, whichtosplit:C0;]
using FastGaussQuadrature
ngh=1000
ghx, ghw = gausshermite(ngh)
ntrials=3
wi_C1 = wi0[ind]
wi_C1[whichtosplit] = wi_C1[whichtosplit]*tau
wi_C1[whichtosplit+1] = wi_C1[whichtosplit+1]*(1-tau)

wi = repmat(wi_C1, 1, 4*ntrials)
mu = zeros(ncomponent1, 4*ntrials)
sigmas = ones(ncomponent1, 4*ntrials)
betas = repmat(betas0, 1, 4*ntrials)
ml = -Inf .* ones(4*ntrials)

i=4
Random.seed!(i)
mu[:, i] = rand(ncomponent1) .* (mu_ub .- mu_lb) .+ mu_lb
sigmas[:, i] = rand(ncomponent1) .* (sigmas_ub .- sigmas_lb) .+ sigmas_lb

Random.seed!(100);LatentGaussianMixtureModel.latentgmm_ctau(X, Y, groupindex, ncomponent1, betas0, wi[:, i], mu[:, i], sigmas[:, i], whichtosplit, tau, ghx, ghw, mu_lb=mu_lb,mu_ub=mu_ub, maxiteration=30, Mmax=50, M_discard=500, an=0., Q_maxiter=2, debuginfo=true)

Random.seed!(100)
for i in 1:4*ntrials
    mu[:, i] = rand(ncomponent1) .* (mu_ub .- mu_lb) .+ mu_lb
    sigmas[:, i] = rand(ncomponent1) .* (sigmas_ub .- sigmas_lb) .+ sigmas_lb
    # if ncomponent1 != 2
    #     #fit gmm on gamma_hat with the starting points, to accelerate the latentgmm_ctau
    #     wi[:, i], mu[:, i], sigmas[:, i], tmp = gmm(gamma0, ncomponent1, wi[:, i], mu[:, i], sigmas[:, i], whichtosplit=whichtosplit, tau=tau, mu_lb=mu_lb,mu_ub=mu_ub, maxiter=1, wifixed=true, sn=sn, an=an)
    # end
    wi[:, i], mu[:, i], sigmas[:, i], betas[:, i], ml[i] = LatentGaussianMixtureModel.latentgmm_ctau(X, Y, groupindex, ncomponent1, betas0, wi[:, i], mu[:, i], sigmas[:, i], whichtosplit, tau, ghx, ghw, mu_lb=mu_lb,mu_ub=mu_ub, maxiteration=8, Mmax=50, M_discard=500, an=0., Q_maxiter=2, debuginfo=false)
end

mlperm = sortperm(ml)
for j in 1:ntrials
    i = mlperm[4*ntrials+1 - j] # start from largest ml
    wi[:, i], mu[:, i], sigmas[:, i], betas[:, i], ml[i] = LatentGaussianMixtureModel.latentgmm_ctau(X, Y, groupindex, ncomponent1, betas0, wi[:, i], mu[:, i], sigmas[:, i], whichtosplit, tau, ghx, ghw, mu_lb=mu_lb,mu_ub=mu_ub, maxiteration=50, Mmax=50, M_discard=500, an=0., Q_maxiter=2, debuginfo=true)
end

mlmax, imax = findmax(ml[mlperm[(3*ntrials+1):4*ntrials]])
imax = mlperm[3*ntrials+imax]

re=LatentGaussianMixtureModel.latentgmm(X, Y, groupindex, ncomponent1, betas[:, imax], wi[:, imax], mu[:, imax], sigmas[:, imax], Mmax=5000, maxiteration=3, initial_iteration=0, an=an)
2(re[5] - ml_C0)