text:
module Nondeterministics

using Random
using StatsFuns: poisinvcdf
using StaticArrays
using CUDAnative
using CuArrays

include("specfuns.jl")

const _distributions = (
    :Categorical, :Poisson,
    :Uniform, :Normal, :Exponential, :Gamma,
    :Dirichlet
)

export Distribution
export params, logpdf
for s in _distributions
    @eval export $s
end

############################
# Distribution generalities
############################

primitive type Distribution 8 end

const Categorical = reinterpret(Distribution, 0x00)
const Poisson     = reinterpret(Distribution, 0x01)
const Uniform     = reinterpret(Distribution, 0x10)
const Normal      = reinterpret(Distribution, 0x11)
const Exponential = reinterpret(Distribution, 0x12)
const Gamma       = reinterpret(Distribution, 0x13)
const Dirichlet   = reinterpret(Distribution, 0x20)

function Base.show(io::IO, d::Distribution)
    for s in _distributions
        d == eval(s) && return print(io, string(s))
    end
    print(io, d)
end

Base.length(::Distribution) = 1
Base.iterate(d::Distribution) = (d, nothing)
Base.iterate(::Distribution, ::Nothing) = nothing

@inline function params(d::Distribution)
    d == Categorical && return Tuple{Vararg{Real}}
    d == Poisson     && return Tuple{Real}
    d == Uniform     && return Tuple{Real,Real}
    d == Normal      && return Tuple{Real,Real}
    d == Exponential && return Tuple{Real}
    d == Gamma       && return Tuple{Real,Real}
    d == Dirichlet   && return Tuple{Vararg{Real}}
end

###############################
# Distribution implementations
###############################

# Base.@irrational log2π 1.8378770664093454836 log(big(2.)*π)

function homogenize(t::Tuple)
    T = promote_type(typeof(t).parameters...)
    Tuple{Vararg{T}}(t)
end

include("scalars.jl")
include("arrays.jl")

@inline function logpdf(d::Distribution)
    for s in _distributions
        d == eval(s) && return eval(Symbol(:logpdf,s))
    end
end

@inline function random(d::Distribution, params...)
    for s in _distributions
        d == eval(s) && return eval(Symbol(:rand,s))(params...)
    end
end
(d::Distribution)(params...) = random(d, params...)

for s in _distributions
    @eval CuArrays.@cufunc $(Symbol(:logpdf,s))(args...) = $(Symbol(:_logpdf,s))(args...)
end

end # module


text:
include("fixtures.jl")

tests = ["message", "appender", "logconfig", "config", "logger_params", "level", "event", "filters"]

for t in tests
    fp = "$t.jl"
    include(fp)
end

ext: jl


text:
function test_adaptive_sketch(A, SketchType, r_true;
    tol=1e-6,
    ρ=1e-4,
    k_factor=0.9,
    condition_number=false,
)
    @testset "Adaptive Sketch Building" begin
        @testset "Norm" begin
            Anys_adapt = RP.adaptive_sketch(A, 2, SketchType; tol=10tol)
            @test round(r_true - rank(Anys_adapt) * 1/k_factor) == 0
            @test opnorm(A - Matrix(Anys_adapt)) <= tol
        end

        if condition_number
            @testset "Condition Number" begin
                Anys_adapt = RP.adaptive_sketch(A, 2, SketchType; tol=100tol, ρ=ρ, condition_number=true)
                @test round(r_true - rank(Anys_adapt) * 1/k_factor) == 0
                U, Λ = Anys_adapt.U, Anys_adapt.Λ
                M = Symmetric( (Λ[end] + ρ)*(U*inv(Λ + ρ*I)*U' + (I - U*U')) * (A + ρ*I) )
                @test (eigmax(M))/(eigmin(M) + ρ) - 1 <= tol
            end
        end
    end

    return nothing
end

@testset "Nystrom Sketch" begin
    # Data
    Random.seed!(0)
    n, r_true = 500, 100
    A = randn(n, r_true)
    A = A*A'

    r = round(Int, r_true * 1.2)
    k = round(Int, 0.9*r)
    λ = eigvals(A; sortby=x->-x)
    Anys = RP.NystromSketch(A, k, r)
    @test size(Anys, 1) == n && size(Anys, 2) == n
    @test rank(Anys) == k
    @test λ[1:k] ≈ svdvals(Anys)

    @test λ[1:k] ≈ eigvals(Anys)
    @test ≈(RP.estimate_norm_E(A, Anys; q=20), opnorm(A - Matrix(Anys)); rtol=1e-2)

    k, r = k ÷ 2, r ÷ 2
    Anys = RP.NystromSketch(A, k, r)
    @test ≈(RP.estimate_norm_E(A, Anys; q=50), opnorm(A - Matrix(Anys)); rtol=2e-1)

    x = randn(n)
    y = Anys * x
    @test y ≈ Matrix(Anys) * x
    z = zeros(n)

    test_adaptive_sketch(A, RP.NystromSketch, 128; k_factor=0.9, ρ=1e-4, tol=1e-6, condition_number=true)
end

@testset "EigenSketch" begin
    # Data
    Random.seed!(0)
    n, r1, r2 = 500, 60, 40
    r_true = r1+r2
    AAT(n, r) = randn(n, r) |> A->A*A'
    A = AAT(n, r1) - AAT(n, r2)

    r = round(Int, r_true * 1.2)
    k = round(Int, 0.9*r)
    A_sketch = RP.EigenSketch(A, k, r; check=true, q=5)
    @test size(A_sketch, 1) == n && size(A_sketch, 2) == n
    @test rank(A_sketch) == k
    λ = eigvals(A; sortby=x->-x) |> x->x[abs.(x) .≥ 1e-8]
    λhat =  eigvals(A_sketch) |> x->x[abs.(x) .≥ 1e-8]

    @test λ ≈ λhat
    @test RP.estimate_norm_E(A, A_sketch; q=20) <= 1e-8 && opnorm(A - Matrix(A_sketch)) <= 1e-8

    k, r = k ÷ 2, r ÷ 2
    A_sketch = RP.EigenSketch(A, k, r; q=100)
    @test ≈(RP.estimate_norm_E(A, A_sketch; q=100), opnorm(A - Matrix(A_sketch)); rtol=2e-1)

    x = randn(n)
    y = A_sketch * x
    @test y ≈ Matrix(A_sketch) * x
    z = zeros(n)

    test_adaptive_sketch(A, RP.EigenSketch, 128; k_factor=0.9, tol=1e-6)
end

@testset "Randomized SVD" begin
    # Data
    Random.seed!(0)
    n, r1, r2 = 500, 60, 40
    r = r1+r2
    AAT(n, r) = randn(n, r) |> A->A*A'
    A = AAT(n, r1) - AAT(n, r2)
    B = randn(n, n÷2) |> svd |> x->(x.U[:, 1:r] * Diagonal(x.S[1:r]) * x.Vt[1:r, :])
    B = B*B'*B

    for M in [A, B]
        r = round(Int, r1 + r2 * 1.2)
        k = round(Int, 0.9*r)
        Mhat = RP.RandomizedSVD(M, k, r; q=5)
        @test size(Mhat, 1) == size(M, 1) && size(Mhat, 2) == size(M, 2)
        @test rank(Mhat) == k
        @test svdvals(M)[1:k] ≈ svdvals(Mhat)
        @test ≈(RP.estimate_norm_E(M, Mhat; q=20), opnorm(M - Matrix(Mhat)); rtol=1e-1)

        k, r = k ÷ 2, r ÷ 2
        Mhat = RP.RandomizedSVD(M, k, r; q=5)
        @test ≈(RP.estimate_norm_E(M, Mhat; q=50), opnorm(M - Matrix(Mhat)); rtol=1e-1)

        x = randn(size(Mhat, 2))
        y = Mhat * x
        @test y ≈ Matrix(Mhat) * x
    end

    test_adaptive_sketch(A, RP.RandomizedSVD, 128; k_factor=0.9, tol=1e-6)

end


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

JLLWrappers.@generate_main_file_header("TempestRemap")
JLLWrappers.@generate_main_file("TempestRemap", UUID("8573a8c5-1df0-515e-a024-abad257ee284"))
end  # module TempestRemap_jll


text:
using Discord
c = Client("BOT-TOKEN-HERE")

function handler(c::Client, e::MessageCreate)
    if e.message.content == "ping--"
        reply(c, e.message, "im basically the same as php ***-julia***")
    end
    if e.message.content == "kys--"
        reply(c, e.message, "Bye! ***-julia***")
        close(c)
    end
    if e.message.content == "stopjl--"
        reply(c, e.message, "**WHEEEE**")
        close(c)
    end
    if e.message.content == "chat--asshole"
        reply(c, e.message, "thats right")
    end
    if e.message.content == "chat--you need help"
        reply(c, e.message, "so do you")
    end
    if e.message.content == "chat--bye"
        reply(c, e.message, "__**GREETINGS!!!**__")
    end
end
function ready(c::Client, e::Ready)
    update_status(
    c,
    Nothing(),
    Nothing(),
    "idle",
    false,
)
end

add_handler!(c, MessageCreate, handler)
add_handler!(c, Ready, ready)
open(c)
wait(c)

max_stars_repo_name: justletterh/variousdiscbotsig

max_issues_repo_name: justletterh/variousdiscbotsig

max_forks_repo_name: justletterh/variousdiscbotsig


text:
using Plots; pyplot()

xVals = [-2,3,5,6,12,14]
yVals = [7,2,9,3,12,3]
n = length(xVals)

V = [xVals[i+1]^(j) for i in 0:n-1, j in 0:n-1]
c = V \ yVals
xGrid = -5:0.01:20
f1(x) = c'*[x^i for i in 0:n-1]

beta0, beta1 = 4.58, 0.17
f2(x) = beta0 + beta1*x

plot(xGrid,f1.(xGrid), c=:blue, label="Polynomial 5th order")
plot!(xGrid,f2.(xGrid),c=:red, label="Linear model")
scatter!(xVals,yVals,
	c=:black, shape=:xcross, ms=8,
	label="Data points", xlims=(-5,20), ylims=(-50,50),
    	xlabel = "x", ylabel = "y")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using BinaryBuilder

name = "PETSc"
version = v"3.15.2"

# Collection of sources required to build PETSc. Avoid using the git repository, it will
# require building SOWING which fails in all non-linux platforms.
sources = [
    ArchiveSource("https://www.mcs.anl.gov/petsc/mirror/release-snapshots/petsc-3.15.2.tar.gz",
    "3b10c19c69fc42e01a38132668724a01f1da56f5c353105cd28f1120cc9041d8"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/petsc*
atomic_patch -p1 $WORKSPACE/srcdir/patches/petsc_name_mangle.patch

BLAS_LAPACK_LIB="${libdir}/libopenblas.${dlext}"

if [[ "${target}" == *-mingw* ]]; then
    #atomic_patch -p1 $WORKSPACE/srcdir/patches/fix-header-cases.patch
    MPI_LIBS="${libdir}/msmpi.${dlext}"
else
    MPI_LIBS="[${libdir}/libmpifort.${dlext},${libdir}/libmpi.${dlext}]"
fi

build_petsc()
{

    if [[ "${3}" == "Int64" ]]; then
        USE_INT64=1
    else
        USE_INT64=0
    fi

    ./configure --prefix=${prefix} \
        CC=${CC} \
        FC=${FC} \
        CXX=${CXX} \
        COPTFLAGS='-O3' \
        CXXOPTFLAGS='-O3' \
        CFLAGS='-fno-stack-protector' \
        FOPTFLAGS='-O3' \
        --with-64-bit-indices=${USE_INT64} \
        --with-debugging=0 \
        --with-batch \
        --PETSC_ARCH=${target}_${1}_${2}_${3} \
        --with-blaslapack-lib=$BLAS_LAPACK_LIB \
        --with-blaslapack-suffix="" \
        --known-64-bit-blas-indices=0 \
        --with-mpi-lib="${MPI_LIBS}" \
        --known-mpi-int64_t=0 \
        --with-mpi-include="${includedir}" \
        --with-sowing=0 \
        --with-precision=${1} \
        --with-scalar-type=${2}

    if [[ "${target}" == *-mingw* ]]; then
        export CPPFLAGS="-Dpetsc_EXPORTS"
    elif [[ "${target}" == powerpc64le-* ]]; then
        export CFLAGS="-fPIC"
        export FFLAGS="-fPIC"
    fi

    make -j${nproc} \
        PETSC_DIR="${PWD}" \
        PETSC_ARCH="${target}_${1}_${2}_${3}" \
        CPPFLAGS="${CPPFLAGS}" \
        CFLAGS="${CFLAGS}" \
        FFLAGS="${FFLAGS}" \
        DEST_DIR="${prefix}" \
        all

    make PETSC_DIR=$PWD PETSC_ARCH=${target}_${1}_${2}_${3} DEST_DIR=$prefix install

    # add suffix to library name
    if [[ "${target}" == *-mingw* ]]; then
        # changing the extension from so to dll.
        mv ${prefix}/lib/libpetsc.so.*.*.* "${libdir}/libpetsc_${1}_${2}_${3}.${dlext}"
    elif [[ "${target}" == *-apple* ]]; then
        mv ${prefix}/lib/libpetsc.*.*.*.${dlext} "${libdir}/libpetsc_${1}_${2}_${3}.${dlext}"
    else
        mv ${prefix}/lib/libpetsc.${dlext}.*.*.* "${libdir}/libpetsc_${1}_${2}_${3}.${dlext}"
    fi
    # Remove useless links
    rm ${prefix}/lib/libpetsc.*
    # Remove duplicated file
    rm ${prefix}/lib/pkgconfig/PETSc.pc
}

build_petsc double real Int32
build_petsc single real Int32
build_petsc double complex Int32
build_petsc single complex Int32
build_petsc double real Int64
build_petsc single real Int64
build_petsc double complex Int64
build_petsc single complex Int64
"""

# We attempt to build for all defined platforms
platforms = expand_gfortran_versions(supported_platforms(exclude=[Platform("i686", "windows")]))

products = [
    LibraryProduct("libpetsc_double_real_Int32", :libpetsc), # Current default build
    LibraryProduct("libpetsc_double_real_Int64", :libpetsc_Float64_Real_Int64),
    LibraryProduct("libpetsc_single_real_Int64", :libpetsc_Float32_Real_Int64),
    LibraryProduct("libpetsc_double_complex_Int64", :libpetsc_Float64_Complex_Int64),
    LibraryProduct("libpetsc_single_complex_Int64", :libpetsc_Float32_Complex_Int64),
    LibraryProduct("libpetsc_double_real_Int32", :libpetsc_Float64_Real_Int32),
    LibraryProduct("libpetsc_single_real_Int32", :libpetsc_Float32_Real_Int32),
    LibraryProduct("libpetsc_double_complex_Int32", :libpetsc_Float64_Complex_Int32),
    LibraryProduct("libpetsc_single_complex_Int32", :libpetsc_Float32_Complex_Int32),
]

dependencies = [
    Dependency("OpenBLAS32_jll"),
    Dependency("MPICH_jll"),
    Dependency("MicrosoftMPI_jll"),
    Dependency("CompilerSupportLibraries_jll"),
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"9")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module Softmax_CEL
    using LoopVectorization

    function opt_func(values)
        return @avxt exp.(values.-maximum(values))/sum(exp.(values.-maximum(values)))
    end

    function func!(output_matrix::Array{Float32}, value_matrix::Array{Float32})
        @avxt for i in axes(value_matrix, 1), j in axes(value_matrix, 2)
            value_matrix[i,j] = ifelse(value_matrix[i,j]>3.0f38, 3.0f38, value_matrix[i,j])
            value_matrix[i,j] = ifelse(value_matrix[i,j]<-3.0f38, -3.0f38, value_matrix[i,j])
        end
        for b in axes(value_matrix, 2)
            output_matrix[:,b] = opt_func(value_matrix[:,b])
        end
    end

    function pre_diff(inputs::Array{Float32}, position::Int64)
        outputs = opt_func(inputs)
        derivative_vector = zeros(Float32, size(inputs))
        @avxt for i in 1:length(derivative_vector)
            derivative_vector[i] = ifelse(i!=position, outputs[i], outputs[i]-1)
        end
        return derivative_vector
    end

    function get_∇biases!(∇biases::Array{Float32}, input_matrix::Array{Float32}, δ::Array{Float32})
        @avxt for i in axes(∇biases, 1), j in axes(∇biases, 2)
            ∇biases[i,j] = 0.0f0
        end

        for b in axes(input_matrix, 2)
            for l in axes(input_matrix, 1)
                if δ[l,b]!=0
                    ∇biases[:,b] += pre_diff(input_matrix[:,b], l)
                end
            end
        end
    end

    function get_name()
        return "Softmax_CEL"
    end
end

module Softmax
    using LoopVectorization

    function opt_func(values)
        return @avxt exp.(values.-maximum(values))/sum(exp.(values.-maximum(values)))
    end

    function func!(output_matrix::Array{Float32}, value_matrix::Array{Float32})
        @avxt for i in axes(value_matrix, 1), j in axes(value_matrix, 2)
            value_matrix[i,j] = ifelse(value_matrix[i,j]>3.0f38, 3.0f38, value_matrix[i,j])
            value_matrix[i,j] = ifelse(value_matrix[i,j]<-3.0f38, -3.0f38, value_matrix[i,j])
        end
        for b in axes(value_matrix, 2)
            output_matrix[:,b] = opt_func(value_matrix[:,b])
        end
    end

    function pre_diff(inputs::Array{Float32}, position::Int64)
        outputs = opt_func(inputs)
        derivative_vector = zeros(Float32, size(inputs))
        @avxt for i in 1:length(derivative_vector)
            derivative_vector[i] = ifelse(i!=position, -outputs[i]*outputs[position], outputs[i]*(1-outputs[i]))
        end
        return derivative_vector
    end

    function get_∇biases!(∇biases::Array{Float32}, input_matrix::Array{Float32}, δ::Array{Float32})
        @avxt for i in axes(∇biases, 1), j in axes(∇biases, 2)
            ∇biases[i,j] = 0.0f0
        end

        for b in axes(input_matrix, 2)
            for l in axes(input_matrix, 1)
                if δ[l,b]!=0
                    ∇biases[:,b] += pre_diff(input_matrix[:,b], l)
                end
            end
        end

        @avxt for i in axes(∇biases, 1), j in axes(∇biases, 2)
            ∇biases[i,j] *= δ[i,j]
        end
    end

    function get_name()
        return "Softmax"
    end
end


text:
## (De)compose a polynomial into even and odd degrees
function evenodd_coeffs(p)
    l = length(p); ## degree + 1
    ## For the zero polynomial, 'odd' would be a 0-element array
    if l < 1
        even = [0]
        odd  = [0]
    else
        even = get_coeffs(p, 0:2:l)
        odd  = get_coeffs(p, 1:2:l)
        # even = reverse([coeff(p, i) for i in 0:2:l])
        # odd  = reverse([coeff(p, i) for i in 1:2:l])
    end
    return even, odd
end

evenodd_polys(p) = parent(p).(evenodd_coeffs(p))

function mix_evenodd(even, odd)
    n = length(even) + length(odd)
    out = Vector{eltype(even)}(undef, n)
    out[1:2:end] .= even
    out[2:2:end] .= odd
    return out
end

mix_evenodd(pe::PolyElem, po::PolyElem) = mix_evenodd(get_coeffs(pe), get_coeffs(po))

principaldeterminantsof(M) = i -> det(M[1:i,1:i])

# Structure for the Hurwitz matrix of a polynomial.
# Inherits matrix structure from MatElem general type.
struct HurwitzMatrix{T} <: MatElem{T}
    p::PolyElem{T}
    # HurwitzMatrix(p::PolyElem{T}) where T = new{T}(parent(p)(normalize_Hurwitz(p)))
end

AA.nrows(p::HurwitzMatrix) = degree(p.p)
AA.ncols(p::HurwitzMatrix) = AA.nrows(p)
AA.length(p::HurwitzMatrix) = AA.nrows(p)^2
AA.isempty(p::HurwitzMatrix) = false
AA.base_ring(p::HurwitzMatrix) = AA.base_ring(p.p)

Base.show(io::IO, p::HurwitzMatrix) = print(io, p.p)
Base.show(io::IO, ::MIME"text/plain", p::HurwitzMatrix) =
    print(io, "HurwitzMatrix for poly (coeffs $(base_ring(p.p))) :\n", p)

Base.convert(::Type{AbstractMatrix}, p::HurwitzMatrix) = Matrix(p)

Base.IndexStyle(::Type{<:HurwitzMatrix}) = IndexCartesian()
function Base.getindex(p::HurwitzMatrix, i::Int, j::Int)
    n = degree(p.p)
    # For i > deg(p), we already have `coeff(p,i)=0`
    return -1 < (n - 2 * i + j) ? coeff(p.p, n - 2 * i + j) : zero(base_ring(p.p))
end

# Functions computing Hurwitz something.
function normalize_Hurwitz(p)
    even, odd = normalize.(evenodd_polys(p))
    return parent(p)(mix_evenodd(even, odd))
end

function Hurwitz_genericmatrix(p)
    n = degree(p)
    R, vars = PolynomialRing(base_ring(p), ["a$i" for i in 0:(n - 1)])
    S, t = PolynomialRing(R, "t")
    Q = S([vars; one(R)])
    H = []
    for i in Hurewitz_odd_indices(p)
        D = Hurwitzdeterminants(Q)
        subs(D, vars .=> coeffs(p))
    end
end

Hurwitzmatrix(p) = HurwitzMatrix(p)[1:degree(p), 1:degree(p)]
Hurwitzmatrix_sylvester(p) = sylvester_matrix(evenodd_polys(p)...)

Hurwitzdeterminants(p) = principaldeterminantsof(HurwitzMatrix(p))

Hurwitz_odd_indices(p) = (degree(p) - 1):-2:2

Hurwitz_odd_determinants(p) = reverse(Hurwitzdeterminants(p).(Hurwitz_odd_indices(p)))

Hurwitz_subresultants(p) = Hurwitz_subresultants_signs(p) .* AA.subresultants_ducos(evenodd_polys(p)...)
Hurwitzdeterminants_subresultant(p) = AA.leading_coefficient.(Hurwitz_subresultants(p))

## The signs of the antidiagonal.
function Hurwitz_subresultants_signs(p)
    n = length(Hurwitz_odd_indices(p))
    out = ones(Int, n)
    goes_neg = iseven(degree(p)) ? [1,2] : [2,3]
    neg_indices = filter(i -> i % 4 in goes_neg, 1:n)
    out[neg_indices] .= -1
    # for i in .%(1:n,4)
    #     if i goes_neg
    #         out[i] = -1
    #     end
    # end
    return out
end

## Numerical results form test_subresultant_ducos
## Deg | length | signs
##  22 |     10 | [-1, -1, 1, 1, -1, -1, 1, 1, -1, -1]
##  20 |      9 | [-1, -1, 1, 1, -1, -1, 1, 1, -1]
##  18 |      8 | [-1, -1, 1, 1, -1, -1, 1, 1]
##
##  21 |     10 | [1, -1, -1, 1, 1, -1, -1, 1, 1, -1]
##  19 |      9 | [1, -1, -1, 1, 1, -1, -1, 1, 1]
##  17 |      8 | [1, -1, -1, 1, 1, -1, -1, 1]

#############################################
#############################################
################################### Stability

## See https://sci-hub.se/https://doi.org/10.1017/S0305004100049872
## Barnett 1971, Proceedings of the Cambridge Philosophical Society, 70(02), 269.
##
##

# Check the signs of the Hurwitz determinants for the Ruth-Hurwitz criteria
# getsigns(D) = unique(sign.(reshape(hcat, coeffs.(coeffs(D)))))
function signs_stabilitymatrices(io,
                                 J;
                                 mindeg::Integer=size(J, 1) - rank(J),
                                 unique_signs::Function=x -> unique(sign.(coeffs(x))))
    P, x = PolynomialRing(base_ring(J), "x")
    p = divexact(charpoly(P, J), x^mindeg)
    L = []
    # println(p)
    range = Hurwitz_odd_indices(p)
    # println(collect(range))
    write(io, "==============================================\n")
    write(io, "==============================================\n")
    write(io, "$(length(range)) expressions will be studied\n")
    for i in reverse(range)
        # This poly belongs to base_ring(J), which is defined in the scope calling this function
        # @computing H = Hurwitzdeterminants(p)(i) "$i-th Hurwitz determinant"
        H = Hurwitzdeterminants(p)(i)
        write(io, "\n-- $i th determinant --\n")
        signs = unique_signs(H)
        write(io, "$(signs)\n")
        if -1 in signs
            push!(L, H)
        end
    end
    write(io, "==============================================\n")
    write(io, "==============================================\n")
    return L
end

function signs_stabilitymatrices(J;
                                 mindeg::Integer=size(J, 1) - rank(J),
                                 unique_signs::Function=x -> unique(sign.(coeffs(x))))
    signs_stabilitymatrices(IOBuffer(), J, mindeg=mindeg, unique_signs=unique_signs)
end

function signs_stabilitymatrices(file::String,
                                 J;
                                 netname=nothing,
                                 mode="a", # append
                                 mindeg::Integer=size(J, 1) - rank(J),
                                 unique_signs::Function=x -> unique(sign.(coeffs(x))))
    open(file, mode) do io
        if !(isnothing(netname))
            write(io, "==============================================\n")
            write(io, "==============================================\n")
            write(io, "$netname\n")
            write(io, "==============================================\n")
            write(io, "==============================================\n")
        end
        signs_stabilitymatrices(io, J, mindeg=mindeg, unique_signs=unique_signs)
    end
end


text:
@testset "HubbardAtom" begin

# test two-site Hubbard model
t, U, μ, β = 0.1, 1.0, 0.5, 0.1
m = Hubbard.hubbardAtom(:fermi, U, μ, β)
@test minimum(m.E) ≈ -μ

m = Hubbard.hubbardAtom2(:fermi, t, U, μ, β)
@test minimum(m.E) ≈ -1.0385164807134504

m = Hubbard.hubbardAtom4(:fermi, t, U, μ, β)
@test minimum(m.E) ≈ -2.109987777274974
end


text:
using Test
using TestSetExtensions
using LinearAlgebra
using Qaintessent

##==----------------------------------------------------------------------------------------------------------------------

# adapted from https://github.com/FluxML/Zygote.jl/blob/master/test/gradcheck.jl
function ngradient(f, xs::AbstractArray...)
    grads = zero.(xs)
    for (x, Δ) in zip(xs, grads), i in 1:length(x)
        δ = sqrt(eps())
        tmp = x[i]
        x[i] = tmp - δ/2
        y1 = f(xs...)
        x[i] = tmp + δ/2
        y2 = f(xs...)
        x[i] = tmp
        Δ[i] = (y2-y1)/δ
        if eltype(x) <: Complex
            # derivative with respect to imaginary part
            x[i] = tmp - im*δ/2
            y1 = f(xs...)
            x[i] = tmp + im*δ/2
            y2 = f(xs...)
            x[i] = tmp
            Δ[i] += im*(y2-y1)/δ
        end
    end
    return grads
end

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "statevector circuit gradients" begin

    # construct parametrized circuit
    N = 4
    A = rand(ComplexF64, 2 ,2)
    U, R = qr(A)
    U = Array(U)
    i = rand(1:N)

    cgc(θ, ϕ, χ, κ, ωn, η, ι, λ) = CircuitGate[
        circuit_gate(3, HadamardGate()),
        circuit_gate(2, RzGate(θ), (1, 4)),
        circuit_gate(2, TGate(), 4),
        circuit_gate(2, 3, SwapGate()),
        circuit_gate(2, SdagGate(), 1),
        circuit_gate(3, PhaseShiftGate(ϕ)),
        circuit_gate(1, RyGate(χ)),
        circuit_gate(4, HadamardGate(), 2),
        circuit_gate(3, X, 1),
        circuit_gate(2, Z, 4),
        circuit_gate(4, Y, 1),
        circuit_gate(3, 1, EntanglementYYGate(κ)),
        circuit_gate(3, RotationGate(ωn)),
        circuit_gate(3, RxGate(η), (1, 4)),
        circuit_gate(2, 3, EntanglementXXGate(ι)),
        circuit_gate(4, 1, EntanglementZZGate(λ)),
        circuit_gate(4, SGate()),
        circuit_gate(2, TdagGate(), 3),
        circuit_gate(i, MatrixGate(U)),
        circuit_gate(2, Z),
        circuit_gate(3, X),
        circuit_gate(2, Y, 4),
        circuit_gate(1, SGate()),
        circuit_gate(2, 3, SwapGate(), 4),
        circuit_gate(3, SdagGate()),
        circuit_gate(2, TGate()),
    ]

    # measurement operators
    meas(M) = MeasurementOperator.([Matrix{Float64}(I, 2^N, 2^N), Hermitian(M)], (Tuple(1:N),))

    # parameter values
    θ = 1.5π
    ϕ = 0.3
    χ = √2
    κ = exp(-0.4)
    η = √3
    ι = 3.33π
    λ = exp(-0.11)
    n = randn(Float64, 3)
    n /= norm(n)
    ωn = 0.2π * n
    M = randn(ComplexF64, 2^N, 2^N)
    M = 0.5*(M + adjoint(M))

    c = Circuit{N}(cgc(θ, ϕ, χ, κ, ωn, η, ι, λ), meas(M))
    # input quantum state
    ψ = Statevector(randn(ComplexF64, 2^N))
    ψl = deepcopy(ψ)
    # fictitious gradients of cost function with respect to circuit output
    Δ = [0.3, -1.2]

    dc = Qaintessent.gradients(c, ψ, Δ)[1]

    f(rθ, rϕ, rχ, rκ, ωn, rη, rι, rλ, M) = dot(Δ, apply!(copy!(ψl, ψ), Circuit{N}(cgc(rθ[], rϕ[], rχ[], rκ[], ωn, rη[], rι[], rλ[]), meas(M))))
    # numeric gradients
    ngrad = ngradient(f, [θ], [ϕ], [χ], [κ], ωn, [η], [ι], [λ], M)
    # symmetrize gradient with respect to measurement operator

    ngrad[end][:] = 0.5*(ngrad[end] + adjoint(ngrad[end]))
    @test all(isapprox.(ngrad,
        (dc.moments[1][2].gate.U.θ,
         dc.moments[4][2].gate.ϕ,
         dc.moments[5][1].gate.θ,
         dc.moments[8][1].gate.θ,
         dc.moments[9][1].gate.nθ,
         dc.moments[10][1].gate.U.θ,
         dc.moments[11][1].gate.θ,
         dc.moments[11][2].gate.θ,
         sparse_matrix(dc.meas[2])), rtol=1e-5, atol=1e-5))
end

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "circuit gradients with moments" begin

   # construct parametrized circuit
   N = 4
    cgc(θ, ϕ, χ, ωn) = Moment[
        Moment([
            circuit_gate(3, HadamardGate()),
            circuit_gate(2, RzGate(θ), (1, 4)),
        ]),
        Moment(circuit_gate(2, 3, SwapGate())),
        Moment(circuit_gate(3, PhaseShiftGate(ϕ))),
        Moment([
            circuit_gate(3, RotationGate(ωn)),
            circuit_gate(1, RyGate(χ)),
        ]),
    ]

   # measurement operators
   meas(M) = MeasurementOperator.([Matrix{Float64}(I, 2^N, 2^N), Hermitian(M)], (Tuple(1:N),))

   # parameter values
   θ = 1.5π
   ϕ = 0.3
   χ = √2
   n = randn(Float64, 3)
   n /= norm(n)
   ωn = 0.2π * n
   M = randn(ComplexF64, 2^N, 2^N)
   M = 0.5*(M + adjoint(M))

   c = Circuit{N}(cgc(θ, ϕ, χ, ωn), meas(M))

   # input quantum state
   ψ = Statevector(randn(ComplexF64, 2^N))
   ψl = deepcopy(ψ)
   # fictitious gradients of cost function with respect to circuit output
   Δ = [0.3, -1.2]

   dc = Qaintessent.gradients(c, ψ, Δ)[1]

   f(rθ, rϕ, rχ, ωn, M) = dot(Δ, apply!(copy!(ψl, ψ), Circuit{N}(cgc(rθ[], rϕ[], rχ[], ωn), meas(M))))
   # numeric gradients
   ngrad = ngradient(f, [θ], [ϕ], [χ], ωn, M)
   # symmetrize gradient with respect to measurement operator
   ngrad[end][:] = 0.5*(ngrad[end] + adjoint(ngrad[end]))

   @test all(isapprox.(ngrad,
       (dc.moments[1][2].gate.U.θ,
        dc.moments[3][1].gate.ϕ,
        dc.moments[4][2].gate.θ,
        dc.moments[4][1].gate.nθ,
        sparse_matrix(dc.meas[2])), rtol=1e-5, atol=1e-5))
end


text:
using BinaryBuilder

name = "OpenLibm"
version = v"0.7.4"
sources = [
    ArchiveSource("https://github.com/JuliaMath/openlibm/archive/v$(version).tar.gz",
                  "eb585204d6c995691e5545ea4ffcb887c159bf05cc5ee1bbf92a6fcdde2fccae"),
]

script = raw"""
# Enter the funzone
cd ${WORKSPACE}/srcdir/openlibm*

# Install into output
flags=("prefix=${prefix}")

# Build ARCH from ${target}
flags+=("ARCH=${target%-*-*}")

# OpenLibm build system doesn't recognize our windows cross compilers properly
if [[ ${target} == *mingw* ]]; then
    flags+=("OS=WINNT")
fi

# Add `CC` override, since OpenLibm seems to think it knows best:
flags+=("CC=$CC")

# Build the library
make "${flags[@]}" -j${nproc}

# Install the library
make "${flags[@]}" install

install_license ./LICENSE.md
"""

# We enable experimental platforms as this is a core Julia dependency
platforms = supported_platforms(;experimental=true)

products = [
    LibraryProduct("libopenlibm", :libopenlibm),
]

dependencies = [
]

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

max_stars_repo_name: ericphanson/Yggdrasil

max_issues_repo_name: ericphanson/Yggdrasil

max_forks_repo_name: ericphanson/Yggdrasil


text:
export get_model_action, get_indiv_action

function get_model_action(features::Array{Float64}, model)
    inputs = get_inputs(features)
    rewards = zeros(Float32, 30)
    for i in eachindex(rewards)
        actions = 0:29 .== i
        rewards[i] = model([inputs; actions])[1]
    end
    if maximum(rewards) > minimum(rewards)
        return argmax(rewards)-1
    end
    rand(0:29)
end

function get_indiv_action(features::Array{Float64}, indiv::CartesianGeneticProgramming.CGPInd)
    inputs = get_inputs(features)
    outputs = CartesianGeneticProgramming.process(indiv, inputs)
    argmax(outputs) - 1
end

function get_indiv_action(features::Array{Float64}, indiv::NEAT.NEATInd)
    inputs = get_inputs(features)
    outputs = NEAT.process(indiv, inputs)
    argmax(outputs) - 1
end


text:
module DashBase
import JSON2
include("components.jl")
include("registry.jl")
export Component, push_prop!, get_name, get_type, get_namespace,
ResourcePkg, Resource, ResourcesRegistry, isdynamic, register_package!, main_registry,
get_dash_dependencies, get_dash_renderer_pkg, get_componens_pkgs,
has_relative_path, has_dev_path, has_external_url, get_type,
get_external_url, get_dev_path, get_relative_path

end # module

max_stars_repo_name: waralex/DashBase

max_issues_repo_name: waralex/DashBase

max_forks_repo_name: waralex/DashBase


text:
function generatemonomials(x,m,;option=0)
    mons = [];
    if option == 0
        mons = monomials(x,[1:m;], mon -> (degree(mon) == m));
    elseif 1 <= option <= n
        mons = monomials(x,[1:m;], mon -> ((degree(mon) == m)&&(degree(mon,x[option])<=1)));
    else
        throw(DomainError(option,"option must be non-negative and <=$(n)"))
    end
    return mons
end

function differentiatepow(u,var::MP.AbstractVariable;power=1)
    res = u;
    if power>1
        res = differentiate(differentiatepow(u,var,power=power-1),var);
    elseif power==1
        res = differentiate(u,var);
    elseif power < 0
        throw(DomainError(power,"power must be non-negative"));
    end
    return res;
end

function dalpha(u,alpha,vars)
    if length(alpha) != length(vars)
        throw(DomainError(alpha,"alpha must be a valid multi-index for the variables x"))
    end
    res = u;
    for i in 1:1:length(alpha)
        res = differentiatepow(res,vars[i],power=alpha[i]);
    end
    return res;
end

function polyderivative(u,poly,vars)
    pterms = terms(poly);
    pcoefs = map(coefficient,pterms);
    pmons = map(monomial,pterms);
    fixedexp = z -> exponents(prod(vars.^0)*z);
    ppows = map(fixedexp,pmons);
    Dvec = map(z->dalpha(u,z,vars),ppows);
    if length(pcoefs)!= length(Dvec)
        throw(ErrorException("poly was not well defined"))
    end
    return sum(pcoefs[i]*Dvec[i] for i in 1:1:length(Dvec));
end

function laplacian(u,vars,;power=1)
    res = u;
    if power > 1
        res = polyderivative(laplacian(u,vars,power=power-1),sum(vars.*vars),vars);
    elseif power == 1
        res = polyderivative(u,sum(vars .* vars),vars);
    elseif power < 0
        throw(DomainError(power,"power must be non-negative"));
    end
    return res;
end

function normpartialder(hvec,var::MP.AbstractVariable,vars,;power=1)
    res = hvec;
    if power>1
        der = normpartialder(hvec,var,vars,power=power-1)
        res = (der[1]-2,der[1]*var*der[2]+sum(vars.*vars)*differentiate(der[2],var));
    elseif power==1
        res = (hvec[1]-2,hvec[1]*var*hvec[2]+sum(vars.*vars)*differentiate(hvec[2],var));
    elseif power < 0
        throw(DomainError(power,"power must be non-negative"));
    end
    return res;
    return
end

function normdalpha(hvec,alpha,vars)
    if length(alpha) != length(vars)
        throw(DomainError(alpha,"alpha must be a valid multi-index for the variables x"))
    end
    res = hvec;
    for i in 1:1:length(alpha)
        res = normpartialder(res,vars[i],vars,power=alpha[i]);
    end
    return res;
end

function normpolyderivative(hvec,poly,vars)
    pterms = terms(poly);
    pcoefs = map(coefficient,pterms);
    pmons = map(monomial,pterms);
    fixedexp = z -> exponents(prod(vars.^0)*z);
    ppows = map(fixedexp,pmons);
    Dvec = map(z->normdalpha(hvec,z,x),ppows);
    if length(pcoefs)!= length(Dvec)
        throw(ErrorException("poly was not well defined"))
    end
    return Dvec
end

function generatebasissphere(m, vars)
    n = length(vars)
    genpoly = (2-n,1);
    monos = generatemonomials(vars,m,option=1);
    fixedexp = z -> exponents(prod(vars.^0)*z);
    mpows = map(fixedexp,monos)
    return map(m->normdalpha(genpoly,m,vars)[2],mpows)
end

function dcoefficient(n,m,i,j)
    d = 0;
    if i == j
        d = 1.0/(2.0^j*factorial(j)*prod([(2*m + n - 2*j - 2*l) for l in 1:j]));
    elseif i < j
        d = -dcoefficient(n,m,i,j-1)/(2*(j-i)*(2*m+n-2-2*i-2*j));
    end
    return d;
end

function harmonicdecomposition(u,vars)
    norm = sum(vars .* vars);
    n =length(vars);
    m = maxdegree(u);
    if mindegree(u) != m
        throw(DomainError(u,"polynomial must be homogeneous"));
    end
    return [ (2*i,sum([dcoefficient(n,m,i,j)*norm^(j-i)*laplacian(u,vars,power=j) for j in i:floor(Int,m/2)])) for i in 0:floor(Int,m/2) ];
end


text:
using Iterators, Base.Test
using Compat

# count
# -----

i = 0
for j = countfrom(0, 2)
	@test j == i*2
	i += 1
	i <= 10 || break
end

# take
# ----

t = take(0:2:8, 10)

@test length(collect(t)) == 5

i = 0
for j = t
	@test j == i*2
	i += 1
end
@test i == 5

i = 0
for j = take(0:2:100, 10)
	@test j == i*2
	i += 1
end
@test i == 10

# drop
# ----

i = 0
for j = drop(0:2:10, 2)
	@test j == (i+2)*2
	i += 1
end
@test i == 4

# cycle
# -----

i = 0
for j = cycle(0:3)
	@test j == i % 4
	i += 1
	i <= 10 || break
end

# repeated
# --------

i = 0
for j = repeated(1, 10)
	@test j == 1
	i += 1
end
@test i == 10
i = 0
for j = repeated(1)
	@test j == 1
	i += 1
	i <= 10 || break
end

# repeatedly
# ----------

i = 0
for j = repeatedly(() -> 1, 10)
	@test j == 1
	i += 1
end
@test i == 10
for j = repeatedly(() -> 1)
	@test j == 1
	i += 1
	i <= 10 || break
end

# chain
# -----

@test collect(chain(1:2:5, 0.2:0.1:1.6)) == [1:2:5; 0.2:0.1:1.6]
@test collect(chain(1:0, 1:2:5, 0.2:0.1:1.6)) == [1:2:5; 0.2:0.1:1.6]

# product
# -------

x1 = 1:2:10
x2 = 1:5
@test collect(product(x1, x2)) == vec([(y1, y2) for y1 in x1, y2 in x2])

@test length(product()) == 1
@test collect(product()) == [()]

# distinct
# --------

x = [5, 2, 2, 1, 2, 1, 1, 2, 4, 2]
@test collect(distinct(x)) == unique(x)

# partition
# ---------

@test collect(partition(take(countfrom(1), 6), 2)) == [(1,2), (3,4), (5,6)]
@test collect(partition(take(countfrom(1), 4), 2, 1)) == [(1,2), (2,3), (3,4)]
@test collect(partition(take(countfrom(1), 8), 2, 3)) == [(1,2), (4,5), (7,8)]
@test collect(partition(take(countfrom(1), 0), 2, 1)) == []

@test_throws ArgumentError partition(take(countfrom(1), 8), 2, 0)

# imap
# ----

function test_imap(expected, input...)
  result = collect(imap(+, input...))
  @test result == expected
end

# Empty arrays
test_imap(
  Any[],
  []
)

test_imap(
  Any[],
  Union{}[]
)

# Simple operation
test_imap(
  Any[1,2,3],
  [1,2,3]
)

# Multiple arguments
test_imap(
  Any[5,7,9],
  [1,2,3],
  [4,5,6]
)

# Different-length arguments
test_imap(
  Any[2,4,6],
  [1,2,3],
  countfrom(1)
)

# groupby
# -------

function test_groupby(input, expected)
  result = collect(groupby(x -> x[1], input))
  @test result == expected
end

# Empty arrays
test_groupby(
  [],
  Any[]
)

test_groupby(
  Union{}[],
  Any[]
)

# Singletons
test_groupby(
  ["xxx"],
  Any[["xxx"]]
)

# Typical operation
test_groupby(
  ["face", "foo", "bar", "book", "baz"],
  Any[["face", "foo"], ["bar", "book", "baz"]]
)

# Trailing singletons
test_groupby(
  ["face", "foo", "bar", "book", "baz", "xxx"],
  Any[["face", "foo"], ["bar", "book", "baz"], ["xxx"]]
)

# Leading singletons
test_groupby(
  ["xxx", "face", "foo", "bar", "book", "baz"],
  Any[["xxx"], ["face", "foo"], ["bar", "book", "baz"]]
)

# Middle singletons
test_groupby(
  ["face", "foo", "xxx", "bar", "book", "baz"],
  Any[["face", "foo"], ["xxx"], ["bar", "book", "baz"]]
)

# subsets
# -------

@test collect(subsets(Any[])) == Any[Any[]]

@test collect(subsets([:a])) == Any[Symbol[], Symbol[:a]]

@test collect(subsets([:a, :b, :c])) ==
      Any[Symbol[], Symbol[:a], Symbol[:b], Symbol[:a, :b], Symbol[:c],
          Symbol[:a, :c], Symbol[:b, :c], Symbol[:a, :b, :c]]

# subsets of size k
# -----------------

@test collect(subsets(Any[],0)) == Any[Any[]]
@test collect(subsets([:a, :b, :c],1)) == Any[Symbol[:a], Symbol[:b], Symbol[:c]]
@test collect(subsets([:a, :b, :c],2)) == Any[Symbol[:a,:b], Symbol[:a,:c], Symbol[:b,:c]]
@test collect(subsets([:a, :b, :c],3)) == Any[Symbol[:a,:b,:c]]
@test collect(subsets([:a, :b, :c],4)) == Any[]
@test length(collect(subsets(collect(1:4),1))) == binomial(4,1)
@test length(collect(subsets(collect(1:4),2))) == binomial(4,2)
@test length(collect(subsets(collect(1:4),3))) == binomial(4,3)
@test length(collect(subsets(collect(1:4),4))) == binomial(4,4)
@test length(collect(subsets(collect(1:4),5))) == binomial(4,5)

# nth element
# -----------
for xs in Any[[1, 2, 3], 1:3, reshape(1:3, 3, 1)]
    @test nth(xs, 3) == 3
    @test_throws BoundsError nth(xs, 0)
    @test_throws BoundsError nth(xs, 4)
end

for xs in Any[take(1:3, 3), drop(-1:3, 2)]
    @test nth(xs, 3) == 3
    @test_throws BoundsError nth(xs, 0)
end

s = subsets([1, 2, 3])
@test_throws BoundsError nth(s, 0)
@test_throws BoundsError nth(s, length(s) + 1)

# Every nth
# ---------

@test collect(takenth([], 10)) == []
@test_throws ArgumentError takenth([], 0)
@test collect(takenth(10:20, 3)) == [12,15,18]
@test collect(takenth(10:20, 1)) == collect(10:20)

## @itr
## ====

# @zip
# ----

macro test_zip(input...)
    n = length(input)
    x = Expr(:tuple, ntuple(i->gensym(), n)...)
    v = Expr(:tuple, map(esc, input)...)
    w = :(zip($(map(esc, input)...)))
    quote
	br = Any[]
	for $x in zip($v...)
	    push!(br, $x)
	end
	mr = Any[]
	@itr for $x in $w
	    push!(mr, $x)
	end
	@test br == mr
    end
end

@test_zip [1,2,3] [:a, :b, :c] ['x', 'y', 'z']
@test_zip [1,2,3] [:a, :b] ['w', 'x', 'y', 'z']
@test_zip [1,2,3] Any[] ['w', 'x', 'y', 'z']
@test_zip [1,2,3] Union{}[] ['w', 'x', 'y', 'z']

# @enumerate
# ----------

macro test_enumerate(input)
    i = gensym()
    x = gensym()
    v = esc(input)
    quote
	br = Any[]
	for ($i,$x) in enumerate($v)
	    push!(br, ($i,$x))
	end
	mr = Any[]
	@itr for ($i,$x) in enumerate($v)
	    push!(mr, ($i,$x))
	end
	@test br == mr
    end
end

@test_enumerate [:a, :b, :c]
@test_enumerate Union{}[]
@test_enumerate Any[]

# @take
# -----

macro test_take(input, n)
    x = gensym()
    v = esc(input)
    quote
	br = Any[]
	for $x in take($v, $n)
	    push!(br, $x)
	end
	mr = Any[]
	@itr for $x in take($v, $n)
	    push!(mr, $x)
	end
	@test br == mr
    end
end

@test_take [:a, :b, :c] 2
@test_take [:a, :b, :c] 5
@test_take [:a, :b, :c] 0
@test_take Any[] 2
@test_take Union{}[] 2
@test_take Any[] 0
@test_take [(:a,1), (:b,2), (:c,3)] 2

# @takestrict
# -----

macro test_takestrict(input, n)
    x = gensym()
    v = esc(input)
    quote
	br = Any[]
	bfailed = false
	try
	    for $x in takestrict($v, $n)
		push!(br, $x)
	    end
	catch
	    bfailed = true
	end

	mr = Any[]
	mfailed = false
	try
	    @itr for $x in takestrict($v, $n)
		push!(mr, $x)
	    end
	catch
	    mfailed = true
	end
	@test br == mr
	@test bfailed == mfailed
    end
end

@test_takestrict [:a, :b, :c] 2
@test_takestrict [:a, :b, :c] 3
@test_takestrict [:a, :b, :c] 5
@test_takestrict [:a, :b, :c] 0
@test_takestrict Any[] 2
@test_takestrict Union{}[] 2
@test_takestrict Any[] 0
@test_takestrict [(:a,1), (:b,2), (:c,3)] 2
@test_takestrict [(:a,1), (:b,2), (:c,3)] 3
@test_takestrict [(:a,1), (:b,2), (:c,3)] 4

# @drop
# -----

macro test_drop(input, n)
    x = gensym()
    v = esc(input)
    quote
	br = Any[]
	for $x in drop($v, $n)
	    push!(br, $x)
	end
	mr = Any[]
	@itr for $x in drop($v, $n)
	    push!(mr, $x)
	end
	@test br == mr
    end
end

@test_drop [:a, :b, :c] 2
@test_drop [:a, :b, :c] 5
@test_drop [:a, :b, :c] 0
@test_drop Any[] 2
@test_drop Union{}[] 2
@test_drop Any[] 0
@test_drop [(:a,1), (:b,2), (:c,3)] 2

# @chain
# -----

macro test_chain(input...)
    x = gensym()
    v = Expr(:tuple, map(esc, input)...)
    w = :(chain($(map(esc, input)...)))
    quote
	br = Any[]
	for $x in chain($v...)
	    push!(br, $x)
	end
	mr = Any[]
	@itr for $x in $w
	    push!(mr, $x)
	end
	@test br == mr
    end
end

@test_chain [1,2,3] [:a, :b, :c] ['x', 'y', 'z']
@test_chain [1,2,3] [:a, :b] ['w', 'x', 'y', 'z']
@test_chain [1,2,3] Any[] ['w', 'x', 'y', 'z']
@test_chain [1,2,3] Union{}[] ['w', 'x', 'y', 'z']
@test_chain [1,2,3] 4 [('w',3), ('x',2), ('y',1), ('z',0)]


text:
module TNTRun

    include("Builder.jl")
    include("RunMode.jl")

    export buildArena
    export fillArena
    export runMode

end # module


text:
using InvariantDistribution

@testset  "Markov matrix: discrete, dense" begin
    @test all(sum(test_mm_discrete_bullmap_cube(n_randpts = 10), 2) .≈ 1)
end

#
# using BenchmarkTools
# all(sum(test_mm_discrete_bullmap_cube(n_randpts = 5, prefilter = true, use_orientations = false, split_factor = 2), 2) .≈ 1)
# all(sum(test_mm_discrete_bullmap_cube(n_randpts = 5, prefilter = true, use_orientations = true, split_factor = 2), 2) .≈ 1)
# all(sum(test_mm_discrete_bullmap_cube(n_randpts = 5, prefilter = false, use_orientations = false, split_factor = 2), 2) .≈ 1)
# all(sum(test_mm_discrete_bullmap_cube(n_randpts = 5, prefilter = false, use_orientations = true, split_factor = 2), 2) .≈ 1)
#
# split = 6
# n_pts = 1000
# t1 = @btime all(sum(test_mm_discrete_bullmap_cube(n_randpts = n_pts, prefilter = true, use_orientations = false, split_factor = split), 2) .≈ 1)
# t2 = @btime all(sum(test_mm_discrete_bullmap_cube(n_randpts = n_pts, prefilter = true, use_orientations = true, split_factor = split), 2) .≈ 1)
# #t3 = @btime all(sum(test_mm_discrete_bullmap_cube(n_randpts = n_pts, prefilter = false, use_orientations = false, split_factor = split), 2) .≈ 1)
# #t4 = @btime all(sum(test_mm_discrete_bullmap_cube(n_randpts = n_pts, prefilter = false, use_orientations = true, split_factor = split), 2) .≈ 1)
#
# size(invariant_bullmap_on_cube(split_factor = split).simplex_inds, 1)
# @time all(sum(test_mm_discrete_bullmap_cube(n_randpts = n_pts, prefilter = true, use_orientations = false, split_factor = 7), 2) .≈ 1)

max_stars_repo_name: kahaaga/InvariantDistribution

max_issues_repo_name: kahaaga/InvariantDistribution

max_forks_repo_name: kahaaga/InvariantDistribution


text:
## Exercise 7-4
## The mathematician Srinivasa Ramanujan found an infinite series that can be used to generate a numerical approximation of 1/π:

## 1/π = 2√2 /9801 ∑_{k=0}^{∞} ((4k)!(1103+26390k))/(k!)^4 * 396^{4k}

## Write a function called estimatepi that uses this formula to compute and return an estimate of π. It should use a while loop to compute terms of the summation until the last term is smaller than 1e-15 (which is Julia notation for 10^{−15}). You can check the result by comparing it to π.
println("Ans: ")

function last_part(k::Int)
    numerator = factorial(4k) * (1103 + 26390 * k)
    denominator = factorial(k)^4 * 396^(4k)

    return numerator/denominator
end

function estimatepi()
    k = 0
    sum = 0.0

    while true
        last_term = last_part(k)

        if last_term < 1e-15
            break
        end

        k += 1
        sum += last_term
    end

    one_over_pi = ((2 * √2) / 9801) * sum

    return 1 / one_over_pi
end

println(estimatepi())

println("End.")


text:
using KnetLayers, BytePairEncoding
import KnetLayers: IndexedDict
import KnetLayers.AutoGrad: full

ENV["DATADEPS_ALWAYS_ACCEPT"] = true

const startsym = "_start_"
const delisym = "_deli_"
const clfsym = "_clf_"
const unksym = "<unk>"
const labels = ["1", "2"]

const vocab_out = IndexedDict(labels)

gpt_model, bpe, vocab_in, tokenizer = GPT.load_gpt_pretrain(12,
                                                    startsym=startsym,
                                                    delisym=delisym,
                                                    clfsym=clfsym,
                                                    unksym=unksym)

set_tokenizer(tokenizer)

gpt = Chain(gpt_model, Dense(input=768,output=1,activation=nothing), Dropout(0.1), CrossEntropyLoss())
embed = gpt[1].embed.tok

function loss(gpt, x1, x2, y, x1_mask, x2_mask, c1_index, c2_index)
    e1 = gpt[1].embed(x1)
    e2 = gpt[1].embed(x2)

    t1 = gpt[1].transformers(e1, x1_mask)
    t2 = gpt[1].transformers(e2, x2_mask)
    lm = GPT.lmloss(embed, x1, t1, x1_mask) + GPT.lmloss(embed, x2, t2, x2_mask)

    c1 = hcat((t1[:,v,i] for (i,v) in enumerate(c1_index))...)
    c2 = hcat((t2[:,v,i] for (i,v) in enumerate(c2_index))...)
    p1 = gpt[2](gpt[3](c1))
    p2 = gpt[2](gpt[3](c2))
    p = vcat(p1, p2)

    cl = gpt[4](p, y)
    #unstable type will cause performance issue
    0.5f0 * lm + cl
end

function predict(gpt, x1, x2, y, x1_mask, x2_mask, c1_index, c2_index)
    e1 = gpt[1].embed(x1)
    e2 = gpt[1].embed(x2)
    t1 = gpt[1].transformers(e1, x1_mask)
    t2 = gpt[1].transformers(e2, x2_mask)
    lm = GPT.lmloss(embed, x1, t1, x1_mask) + GPT.lmloss(embed, x2, t2, x2_mask)
    c1 = hcat((t1[:,v,i] for (i,v) in enumerate(c1_index))...)
    c2 = hcat((t2[:,v,i] for (i,v) in enumerate(c2_index))...)
    p1 = gpt[2](gpt[3](c1))
    p2 = gpt[2](gpt[3](c2))
    p  = vcat(p1, p2)
    cl = gpt[4](p, y)
    #unstable type will cause performance issue
    p
end

const rocs = Datasets.StoryCloze()
const ps = params(gpt)
const opt = Adam(lr=6.25e-5)
setoptim!(M, optimizer) = for p in params(M); p.opt = deepcopy(optimizer); end
lrdecay!(M, decay::Real) = for p in params(M); p.opt.lr = p.opt.lr*decay; end

function eval(gpt,data,vocab_in, vocab_out; Batch=8)
    println("eval:")
    datas = dataset(Datasets.Test, data)
    al,i = 0.0, 0
    while (batch = get_batch(datas, Batch)) !== nothing
        b1, b2, c1i, c2i, y, b1_mask, b2_mask = preprocess(batch, vocab_in, vocab_out)
        p = predict(gpt, b1, b2, y, b1_mask, b2_mask, c1i, c2i)
        a = accuracy(p, y)
        al += a
        i  += 1
    end
    total =  al /= i
    @show total
    return total
end

function train1!(gpt, data, vocab_in, vocab_out; opt=nothing, epoch=10, Batch=8)
    #global Batch, rocs, opt, ps
    ps  = params(gpt)
    if opt !== Nothing
        setoptim!(gpt,opt)
    end
    lss = typemax(Float64)
    for e = 1:epoch
        println("start training: $e")
        datas  = dataset(Datasets.Train, data)
        i, al  = 0, 0.0
        grads = nothing

        while (batch = get_batch(datas,Batch)) !== nothing
            b1, b2, c1i, c2i, y, b1_mask, b2_mask, = preprocess(batch, vocab_in, vocab_out)
            J =  @diff loss(gpt, b1, b2, y, b1_mask, b2_mask, c1i, c2i)
            if grads == nothing
                grads = map(w->full(grad(J,w)),params(J))
            else
                for (k,w) in enumerate(ps)
                    grads[k] += grad(J,w)
                end
            end

            if i % 4 == 1
                for (k,w) in enumerate(ps)
                    update!(value(w),grads[k]/4.0f0,w.opt)
                end
                for g in grads
                    fill!(g,0.0f0)
                end
            end

            i  += 1
            al += value(J)
            if i%100==1
                println(al/i)
            end
        end
        nwlss = al/i
        if nwlss < lss
            lss=nwlss
        else
            lrdecay!(gpt,0.002)
            println("lr decay!")
        end
        total = eval(gpt, data, vocab_in, vocab_out)
    end
end

# function accuracy(p, y)
#     pred = onecold(collect(p))
#     label = onecold(collect(y))
#     sum(pred .== label) / length(label)
# end

#
#
# const embed = gpt.embed.embeddings.tok
# const clf = arrtype(Dense(input=768, output=1))
# const ansdrop = Dropout(0.1)
# const loss = CrossEntropyLoss()
#
#
# function loss(x1, x2, y, x1_mask, x2_mask, c1_index, c2_index; train=false)
#     e1 = gpt.embed(x1)
#     e2 = gpt.embed(x2)
#     t1 = gpt.transformers(e1, x1_mask)
#     t2 = gpt.transformers(e2, x2_mask)
#     lm = lmloss(embed, vocab[x1], t1, x1_mask) + lmloss(embed, vocab[x2], t2, x2_mask)
#     c1 = t1[:,c1_index]
#     c2 = t2[:,c2_index]
#
#     p1 = clf(c1)
#     p2 = clf(c2)
#     p = vcat(p1, p2)
#     p = ansdrop(p, 1)
#     cl = loss(p, y)
#     if train
#         convert(eltype(arrtype), 0.5) * lm
#     else
#         convert(eltype(arrtype), 0.5) * lm + cl,p
#     end
# end
#
# const rocs = StoryCloze()
# setoptim!(params([gpt, clf]); lr=ADAM(6.25e-5))
#
# train!(args["epoch"])


text:
using Documenter, PackageCompilerX

makedocs(
    sitename = "PackageCompilerX",
    pages = Any[
        "Home" => "index.md",
        "Examples" => Any[
            "examples/ohmyrepl.md"
        ]
    ]
)

deploydocs(
    repo = "github.com/KristofferC/PackageCompilerX.jl.git",
)


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Julia 1.2.0
#     language: julia
#     name: julia-1.2
# ---

# ## 1D Vlasov–Ampere system
#
# $$
# \frac{\partial f}{\partial t} + \upsilon \frac{\partial f}{\partial x}
# - E(t,x) \frac{\partial f}{\partial \upsilon} = 0
# $$
#
# $$
# \frac{\partial E}{\partial t} = - J = \int f\upsilon \; d\upsilon
# $$

# ## Algorithm
#
# - For each $j$ compute discrete Fourier transform in $x$ of $f^n(x_i,\upsilon_j)$ yielding $f_k^n(\upsilon_j)$,
#
# - For $ k \neq 0 $
#
#     - Compute
#
#     $$f^{n+1}_k(\upsilon_j) = e^{−2i\pi k \upsilon
#     \Delta t/L} f_n^k(\upsilon_j),$$
#
#     - Compute
#
#     $$\rho_k^{n+1} = \Delta \upsilon \sum_j􏰄 f^{n+1}_k(\upsilon_j),$$
#
#     - Compute
#
#     $$E^{n+1}_k = \rho^{n+1}_k L/(2i\pi k \epsilon_0),$$
#
# - For $k = 0$ do nothing:
#
# $$f_{n+1}(\upsilon_j) = f^n_k(\upsilon_j), E^{n+1}_k = E^n_k$$.
#
# - Perform inverse discrete Fourier transform of $E^{n+1}_k$ and for each $j$ of $f^{n+1}_k (\upsilon_j)$.

using ProgressMeter, FFTW, Plots, LinearAlgebra
using BenchmarkTools, Statistics

# +
"""

    UniformMesh(start, stop, length)

    1D uniform mesh data.

    length   : number of points
    length-1 : number of cells

    To remove the last point, set endpoint=false

"""
struct UniformMesh

   start    :: Float64
   stop     :: Float64
   length   :: Int
   step     :: Float64
   points   :: Vector{Float64}
   endpoint :: Bool

   function UniformMesh(start, stop, length::Int; endpoint=true)

       if (endpoint)
           points = range(start, stop=stop, length=length)
       else
           points = range(start, stop=stop, length=length+1)[1:end-1]
       end
       step = points[2]-points[1]

       new( start, stop, length, step, points, endpoint)

   end

end

# +
using FFTW, LinearAlgebra

abstract type AbstractAdvection end
abstract type AdvectionType end

# +
import FFTW: fft!, ifft!
import LinearAlgebra: transpose

export Ampere

struct Ampere end

"""

    advection!( fᵀ, mesh1, mesh2, E, dt, type, axis )

    ∂f/∂t − υ∂f/∂x  = 0
    ∂E/∂t = −J = ∫ fυ dυ
    ∂f/∂t − E(x) ∂f/∂υ  = 0

"""
struct AmpereAdvection <: AbstractAdvection

    mesh :: UniformMesh
    kx   :: Vector{Float64}

    function AmpereAdvection( mesh )

        nx  = mesh.length
        dx  = mesh.step
        Lx  = mesh.stop - mesh.start
        kx  = zeros(Float64, nx)
        kx .= 2π/Lx * [0:nx÷2-1;-nx÷2:-1]
        new( mesh, kx)

    end

end

function (adv :: AmpereAdvection)( fᵗ  :: Array{ComplexF64,2},
		                           e   :: Vector{ComplexF64},
		                           dt  :: Float64 )
    fft!(fᵗ, 1)
    fᵗ .= fᵗ .* exp.(-1im * dt * adv.kx * transpose(e))
    ifft!(fᵗ, 1)

end

function (adv :: AmpereAdvection)( f   :: Array{ComplexF64,2},
		                           e   :: Vector{ComplexF64},
		                           v   :: Vector{Float64},
		                           dt  :: Float64 )

    ev = exp.(-1im*dt * adv.kx * transpose(v))

    fft!(f,1)
    f .= f .* ev
    dv = v[2]-v[1]
    ρ = dv * vec(sum(f,dims=2))
    for i in 2:length(e)
        e[i] = -1im * ρ[i] ./ adv.kx[i]
    end
    e[1] = 0.0
    ifft!(f,1)
    ifft!(e)
    e .= real(e)
end

# +
using FFTW, LinearAlgebra

"""
    bspline(p, j, x)

Return the value at x in [0,1[ of the B-spline with integer nodes of degree p with support starting at j.
Implemented recursively using the [De Boor's Algorithm](https://en.wikipedia.org/wiki/De_Boor%27s_algorithm)

```math
B_{i,0}(x) := \\left\\{
\\begin{matrix}
1 & \\mathrm{if}  \\quad t_i ≤ x < t_{i+1} \\\\
0 & \\mathrm{otherwise}
\\end{matrix}
\\right.
```

```math
B_{i,p}(x) := \\frac{x - t_i}{t_{i+p} - t_i} B_{i,p-1}(x)
+ \\frac{t_{i+p+1} - x}{t_{i+p+1} - t_{i+1}} B_{i+1,p-1}(x).
```
"""
function bspline(p::Int, j::Int, x::Float64)
   if p == 0
       if j == 0
           return 1.0
       else
           return 0.0
       end
   else
       w = (x - j) / p
       w1 = (x - j - 1) / p
   end
   return (w * bspline(p - 1, j, x) + (1 - w1) * bspline(p - 1, j + 1, x))
end

export BSpline

struct BSpline <: AdvectionType
    p :: Int64
end

"""
    BsplineAdvection(p, mesh)

Advection type

"""
mutable struct BSplineAdvection <: AbstractAdvection

    p        :: Int64
    mesh     :: UniformMesh
    modes    :: Vector{Float64}
    eig_bspl :: Vector{Float64}
    eigalpha :: Vector{Complex{Float64}}

    function BSplineAdvection( p, mesh )

        println(" Create BSL advection with bspline of degree $p ")
        @show nx  = mesh.length

        modes     = zeros(Float64, nx)
        modes    .= [2π * i / nx for i in 0:nx-1]
        eig_bspl  = zeros(Float64, nx)
        eig_bspl  = zeros(Float64, nx)
        eig_bspl .= bspline(p, -div(p+1,2), 0.0)
        for i in 1:div(p+1,2)-1
            eig_bspl .+= bspline(p, i-(p+1)÷2, 0.0) * 2 .* cos.(i * modes)
        end
        eigalpha  = zeros(Complex{Float64}, nx)
        new( p, mesh, modes, eig_bspl, eigalpha )
    end

end

"""

    advection( f, v, dt)

Instantiate an advection type

```julia
advection! = Advection( p, mesh)

advection!(f, v, dt)
```

"""
function (adv :: BSplineAdvection)(f    :: Array{Complex{Float64},2},
                                   v    :: Vector{Float64},
                                   dt   :: Float64)

   nx = adv.mesh.length
   nv = length(v)
   dx = adv.mesh.step

   fft!(f,1)

   for j in 1:nv

      @inbounds alpha = dt * v[j] / dx
      # compute eigenvalues of cubic splines evaluated at displaced points
      ishift = floor(-alpha)
      beta   = -ishift - alpha
      fill!(adv.eigalpha,0.0im)
      for i in -div(adv.p-1,2):div(adv.p+1,2)
         adv.eigalpha .+= (bspline(adv.p, i-div(adv.p+1,2), beta)
                        .* exp.((ishift+i) * 1im .* adv.modes))
      end

      # compute interpolating spline using fft and properties of circulant matrices

      @inbounds f[:,j] .*= adv.eigalpha ./ adv.eig_bspl

   end

   ifft!(f,1)

end

# +
function Advection( b::BSpline, mesh:: UniformMesh )

    BSplineAdvection( b.p, mesh)

end

function Advection( b::Ampere, mesh:: UniformMesh )

    AmpereAdvection(mesh)

end

# +
"""

    compute_rho( mesh, f)

    Compute charge density

    ρ(x,t) = ∫ f(x,v,t) delta2

    return ρ - ρ̄ if neutralized=true

"""
function compute_rho(meshv::UniformMesh, f, neutralized=true)

   local dv  = meshv.step
   ρ = dv * vec(sum(real(f), dims=2))
   if (neutralized)
       ρ .- mean(ρ)
   else
       ρ
   end

end

# +
function compute_e(mesh::UniformMesh, ρ)

   local n = mesh.length
   local k =  2π / (mesh.stop - mesh.start)
   local modes  = zeros(Float64, n)
   modes .= k * vcat(0:n÷2-1,-n÷2:-1)
   modes[1] = 1.0
   ρ̂ = fft(ρ)./modes
   vec(real(ifft(-1im*ρ̂)))

end

# +
function landau( ϵ, kx, x, v )

    (1.0.+ϵ*cos.(kx*x))/sqrt(2π) .* transpose(exp.(-0.5*v.*v))

end
# -

function vlasov_ampere( nx, nv, xmin, xmax, vmin, vmax , tf, nt)

    meshx = UniformMesh(xmin, xmax, nx)
    meshv = UniformMesh(vmin, vmax, nv)

    # Initialize distribution function
    x = meshx.points
    v = meshv.points
    ϵ, kx = 0.001, 0.5

    f = zeros(Complex{Float64},(nx,nv))
    fᵀ= zeros(Complex{Float64},(nv,nx))

    f .= landau( ϵ, kx, x, v)

    transpose!(fᵀ,f)

    ρ  = compute_rho(meshv, f)
    e  = zeros(ComplexF64, nx)
    e .= compute_e(meshx, ρ)

    nrj = Float64[]

    dt = tf / nt

    advection_x! = Advection( Ampere(), meshx )
    advection_v! = Advection( Ampere(), meshv )

    @showprogress 1 for i in 1:nt
        advection_v!(fᵀ, e,  0.5*dt)
        transpose!(f,fᵀ)
        advection_x!( f, e, v, dt)
        push!(nrj, log(sqrt((sum(e.^2))*meshx.step)))
        transpose!(fᵀ,f)
        advection_v!(fᵀ, e,  0.5*dt)
    end
    nrj
end

function vlasov_poisson( nx, nv, xmin, xmax, vmin, vmax , tf, nt)

    meshx = UniformMesh(xmin, xmax, nx)
    meshv = UniformMesh(vmin, vmax, nv)

    # Initialize distribution function
    x = meshx.points
    v = meshv.points
    ϵ, kx = 0.001, 0.5

    e = zeros(Float64, nx)
    ρ = zeros(Float64, nx)
    f = zeros(Complex{Float64},(nx,nv))
    fᵀ= zeros(Complex{Float64},(nv,nx))

    f .= landau( ϵ, kx, x, v)

    nrj = Float64[]

    dt = tf / nt

    advection_x! = Advection( BSpline(5), meshx )
    advection_v! = Advection( BSpline(5), meshv )

    @showprogress 1 for i in 1:nt
        advection_x!(f, v,  0.5dt)
        ρ  .= compute_rho(meshv, f)
        e  .= compute_e(meshx, ρ)
        push!(nrj, log(sqrt((sum(e.^2))*meshx.step)))
        transpose!(fᵀ, f)
        advection_v!(fᵀ, e, dt)
        transpose!(f, fᵀ)
        advection_x!( f, v, 0.5dt)
    end
    nrj
end

nx, nv = 256, 256
xmin, xmax =  0., 4*π
vmin, vmax = -6., 6.
tf = 60
nt = 600

t =  range(0,stop=tf,length=nt)
plot(t, -0.1533*t.-5.48)
plot!(t, vlasov_ampere(nx, nv, xmin, xmax, vmin, vmax, tf, nt), label=:ampere )
plot!(t, vlasov_poisson(nx, nv, xmin, xmax, vmin, vmax, tf, nt), label=:poisson )

max_stars_repo_name: JuliaVlasov/VlasovNotebooks

max_issues_repo_name: JuliaVlasov/VlasovNotebooks

max_forks_repo_name: JuliaVlasov/VlasovNotebooks


text:
# using Ipopt
include(joinpath(pwd(),"src/Ipopt.jl"))
using Test, ForwardDiff

function row_col!(row,col,r,c)
    for cc in c
        for rr in r
            push!(row,convert(Int,rr))
            push!(col,convert(Int,cc))
        end
    end
    return row, col
end

obj(z) = (z[1] - 5)^2 + (2*z[2] + 1)^2 + λ_global'*z[9:11] + 0.5*mu_prev*z[9:11]'*z[9:11]

function con!(g,z)
	g[1] = 2*(z[2] - 1) - 1.5*z[2] + z[3] - 0.5*z[4] + z[5]
	g[2] = 3*z[1] - z[2] - 3.0 - z[6]
	g[3] = -z[1] + 0.5*z[2] + 4.0 - z[7]
	g[4] = -z[1] - z[2] + 7.0 - z[8]
	g[5] = z[9] - z[3]*z[6]
	g[6] = z[10] - z[4]*z[7]
	g[7] = z[11] - z[5]*z[8]

	return g
end

function eval_f(x::Vector{Float64})
  x_global .= x
  return obj(x)
end

function eval_g(x::Vector{Float64}, g::Vector{Float64})
	con!(g,x)

	return g
end

function eval_grad_f(x::Vector{Float64}, grad_f::Vector{Float64})
	grad_f .= ForwardDiff.gradient(obj,x)
end

function eval_jac_g(x::Vector{Float64}, mode, rows::Vector{Int32}, cols::Vector{Int32}, values::Vector{Float64})
  if mode == :Structure
		row = []
		col = []
    row_col!(row,col,(1:m),(1:n))
		rows .= vec(row)
		cols .= vec(col)
  else

  	values .= vec(ForwardDiff.jacobian(con!,zeros(7),x))
  end
end

# This tests callbacks.
function intermediate(alg_mod::Int, iter_count::Int,
  obj_value::Float64, inf_pr::Float64, inf_du::Float64, mu::Float64,
  d_norm::Float64, regularization_size::Float64, alpha_du::Float64, alpha_pr::Float64,
  ls_trials::Int)

	println("CALLBACK")
  println("   mu_prev: $(mu_prev)")
  println("   mu: $(mu)")
  println("   x: $(x_global)")

  if mu != mu_prev
    println("-mu difference-")
    println("-AL update-")
    λ_global .+= 1/mu_prev*x_global[9:11]
    global mu_prev = mu
  end

  return true#!(norm(x_global[9:11],Inf) < 1.0e-3 && norm(eval_g(x_global,zeros(7))) < 1.0e-3) #iter_count < 1  # Interrupts after one iteration.
end

n = 11
m = 7
x_L = zeros(n)
x_U = 1.0e8*ones(n)

g_L = zeros(m)
g_U = zeros(m)
# g_U[5:7] .= 1.0e8

prob = createProblem(n, x_L, x_U, m, g_L, g_U, m*n, n*n,
                     eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h)
addOption(prob,"hessian_approximation","limited-memory")
addOption(prob,"mu_strategy","monotone")

addOption(prob,"constr_viol_tol",1.0e-3)
addOption(prob,"tol",1.0e-3)# addOption(prob,"")
setIntermediateCallback(prob, intermediate)

global x_global = zeros(n)
global λ_global = zeros(3)
global mu_prev = 0.1

prob.x .= 1.0e-5*rand(n)
x_global .= prob.x
@time solvestat = solveProblem(prob)

λ_global
mu_prev
# @test Ipopt.ApplicationReturnStatus[solvestat] == :User_Requested_Stop

prob.x[9:11]


text:
using Libdl, Base.Sys

depsfile = joinpath(dirname(@__FILE__), "deps.jl")

if isfile(depsfile)
    rm(depsfile)
end

function write_depsfile(knpath, libpath)
    open(depsfile,"w") do f
        print(f,"const libknitro = ")
        show(f, libpath)
        println(f)
        print(f,"const amplexe = ")
        show(f, joinpath(knpath, "knitroampl", "knitroampl"))
        println(f)
    end
end

paths_to_try = String[]

libname = string(Sys.iswindows() ? "" : "lib", "knitro", ".", Libdl.dlext)

# try to load absolute path before
if haskey(ENV, "KNITRODIR")
    push!(paths_to_try, ENV["KNITRODIR"])
end

push!(paths_to_try, libname)

global found_knitro = false
for path in paths_to_try
    l = joinpath(path, "lib", libname)
    d = Libdl.dlopen_e(l)
    if d != C_NULL
        global found_knitro = true
        write_depsfile(path, l)
        break
    end
end

if !found_knitro
    error("Unable to locate KNITRO installation, " *
          "please check your enviroment variable KNITRODIR.")
end


text:
# binary operations that can not be simplified by a unitary operations
struct SimpleBinaryRule{ix1,ix2,iy} <: EinRule{2} end
function SimpleBinaryRule(code::EinCode)
    ixs, iy = getixs(code), getiy(code)
    @assert length(ixs)==2 "fail to construct simple binary rule from $code"
    SimpleBinaryRule{ixs[1], ixs[2], iy}()
end
SimpleBinaryRule(ix1, ix2, iy) = SimpleBinaryRule{ix1, ix2, iy}()

@inline function _add_patch(::SimpleBinaryRule{ix1,ix2,iy}) where {ix1,ix2,iy}
    SimpleBinaryRule{(ix1...,'l'), (ix2...,'l'), (iy...,'l')}()
end
@inline _add_patch(::DefaultRule) = DefaultRule()

function match_rule_binary(ix1, ix2, iy)
    Nx1, Nx2, Ny = length(ix1), length(ix2), length(iy)
    if (Nx1 + Nx2 + Ny) % 2 == 0 # no batch
        _match_simple2(ix1,ix2,iy,Nx1,Nx2,Ny)
    elseif Nx1>0 && Nx2>0 && Ny>0 && ix1[Nx1]==ix2[Nx2]==iy[Ny]
        rule = _match_simple2(ix1,ix2,iy,Nx1-1,Nx2-1,Ny-1)
        _add_patch(rule)
    else
        DefaultRule()
    end
end

function _match_simple2(ix1, ix2, iy, Nx1, Nx2, Ny)
    if Nx1==0
        if (Ny==Nx2==0)
            return SimpleBinaryRule((), (), ())
        elseif (Ny==Nx2==1 && ix2[1] == iy[1])
            return SimpleBinaryRule((), ('k',), ('k',))
        end
    elseif Nx1==1
        if (Nx2==0 && Ny==1 && iy[1]==ix1[1])
            return SimpleBinaryRule(('i',), (), ('i',))
        elseif (Nx2==1 && Ny==0 && ix1[1]==ix2[1])
            return SimpleBinaryRule(('j',), ('j',), ())
        elseif Nx2==1 && Ny==2
            if (iy[1]==ix1[1] && iy[2]==ix2[1])
                return SimpleBinaryRule(('i',), ('k',), ('i','k'))
            elseif iy[1]==ix2[1] && iy[2]==ix1[1]
                return SimpleBinaryRule(('i',), ('k',), ('k','i'))
            end
        elseif Nx2==2 && Ny==1
            if ix2[1]==ix1[1] && ix2[2]==iy[1]
                return SimpleBinaryRule(('j',), ('j','k'), ('k',))
            elseif ix2[1]==iy[1] && ix2[2]==ix1[1]
                return SimpleBinaryRule(('j',), ('k','j'), ('k',))
            end
        end
    elseif Nx1==2
        if Nx2==1 && Ny==1
            if ix1[1]==ix2[1] && ix1[2]==iy[1]
                return SimpleBinaryRule(('j','i'), ('j',), ('i',))
            elseif ix1[1]==iy[1] && ix1[2]==ix2[1]
                return SimpleBinaryRule(('i','j'), ('j',), ('i',))
            end
        elseif (Nx2==2 && Ny==2)
            if ix1[1]==ix2[1] && ix1[2]==iy[1] && ix2[2]==iy[2]
                return SimpleBinaryRule(('j','i'), ('j','k'), ('i','k'))
            elseif ix1[1]==ix2[2] && ix1[2]==iy[1] && ix2[1]==iy[2]
                return SimpleBinaryRule(('j','i'), ('k','j'), ('i','k'))
            elseif ix1[1]==ix2[1] && ix1[2]==iy[2] && ix2[2]==iy[1]
                return SimpleBinaryRule(('j','i'), ('j','k'), ('k','i'))
            elseif ix1[1]==ix2[2] && ix1[2]==iy[2] && ix2[1]==iy[1]
                return SimpleBinaryRule(('j','i'), ('k','j'), ('k','i'))
            elseif ix1[2]==ix2[1] && ix1[1]==iy[1] && ix2[2]==iy[2]
                return SimpleBinaryRule(('i','j'), ('j','k'), ('i','k'))
            elseif ix1[2]==ix2[2] && ix1[1]==iy[1] && ix2[1]==iy[2]
                return SimpleBinaryRule(('i','j'), ('k','j'), ('i','k'))
            elseif ix1[2]==ix2[1] && ix1[1]==iy[2] && ix2[2]==iy[1]
                return SimpleBinaryRule(('i','j'), ('j','k'), ('k','i'))
            elseif ix1[2]==ix2[2] && ix1[1]==iy[2] && ix2[1]==iy[1]
                return SimpleBinaryRule(('i','j'), ('k','j'), ('k','i'))
            end
        end
    end
    return DefaultRule()
end

function einsum(rule::SimpleBinaryRule, ixs, iy, xs::NTuple{2, Any}, size_dict)
    @debug rule size.(xs)
    einsum(rule, xs)
end
# Code is a binary representation of `(O1,I,O2,B)`.
# Because the time complexity of `GEMM` and `BatchedGEMM` are higher than space complexity, we allow `permutedims`.
# We reduce the contraction to these basic forms through `permutedims` and reshape,
# because we assume most using cases have both inner and outer degrees on freedom.

# ,-> : 000
# S = 1
# T = 1
function einsum(::SimpleBinaryRule{(),(), ()}, xs::NTuple{2, Any})
    asarray(asscalar(xs[1]) * xs[2], xs[1])
end

# i,->i : 100
# S = N
# T = N
function einsum(::SimpleBinaryRule{('i',),(), ('i',)}, xs::NTuple{2, Any})
    xs[1] .* asscalar(xs[2])
end

# j,j-> : 010
# S = N
# T = N
function einsum(::SimpleBinaryRule{('j',), ('j',), ()}, xs::NTuple{2, Any})
    asarray(transpose(xs[1]) * xs[2], xs[1])
end

# ,k->k : 001
# S = N
# T = N
@inline function einsum(::SimpleBinaryRule{(), ('k',), ('k',)}, xs::NTuple{2, Any})
    einsum(SimpleBinaryRule{('i',),(),('i',)}(), (xs[2], xs[1]))
end

# j,jk->k : 011
# S = N^2
# T = N^2
function einsum(::SimpleBinaryRule{('j',), ('j','k'), ('k',)}, xs::NTuple{2, Any})
    vec(transpose(xs[1]) * xs[2])
end
function einsum(::SimpleBinaryRule{('j',), ('k','j'), ('k',)}, xs::NTuple{2, Any})
    xs[2] * xs[1]
end

# ij,j->i : 110
# S = N^2
# T = N^2
@inline function einsum(::SimpleBinaryRule{('i','j'),('j',), ('i',)}, xs::NTuple{2, Any})
    einsum(SimpleBinaryRule{('j',),('k','j'), ('k',)}(), (xs[2], xs[1]))
end
@inline function einsum(::SimpleBinaryRule{('j','i'),('j',), ('i',)}, xs::NTuple{2, Any})
    einsum(SimpleBinaryRule{('j',),('j','k'), ('k',)}(), (xs[2], xs[1]))
end

# i,k->ik : 101
# S = N^2
# T = N^2
function einsum(::SimpleBinaryRule{('i',), ('k',), ('i','k')}, xs::NTuple{2, Any})
    xs[1] * transpose(xs[2])
end
@inline function einsum(::SimpleBinaryRule{('i',), ('k',),('k','i')}, xs::NTuple{2, Any})
    einsum(SimpleBinaryRule{('i',),('k',),('i','k')}(), (xs[2], xs[1]))
end

# 000
function einsum(::SimpleBinaryRule{('l',),('l',), ('l',)}, xs::NTuple{2, Any})
    xs[1] .* xs[2]
end

# 100
function einsum(::SimpleBinaryRule{('i','l'),('l',), ('i','l')}, xs::NTuple{2, Any})
    xs[1] .* transpose(xs[2])
end

# 001
@inline function einsum(::SimpleBinaryRule{('l',), ('k','l'), ('k','l')}, xs::NTuple{2, Any})
    einsum(SimpleBinaryRule{('i','l'),('l',),('i','l')}(), (xs[2], xs[1]))
end

# 010
function einsum(::SimpleBinaryRule{('j','l'), ('j','l'), ('l',)}, xs::NTuple{2, Any})
    a, b = xs
    dropdims(mapreduce(*, +, a, b; dims=1); dims=1)
end

# 101
function einsum(::SimpleBinaryRule{('i','l'), ('k','l'), ('i','k','l')}, xs::NTuple{2, Any})
    a, b = xs
    _batched_gemm('N', 'N', reshape(a, size(a, 1), 1, size(a, 2)), reshape(b, 1, size(b, 1), size(b, 2)))
end
@inline function einsum(::SimpleBinaryRule{('i','l'), ('k','l'), ('k','i','l')}, xs::NTuple{2, Any})
    einsum(SimpleBinaryRule{('i','l'),('k','l'), ('i','k','l')}(), (xs[2], xs[1]))
end

# 011
function einsum(::SimpleBinaryRule{('j','l'), ('j','k','l'), ('k','l')}, xs::NTuple{2, Any})
    reshape(_batched_gemm('N', 'N', reshape(xs[1], 1, size(xs[1],1), size(xs[1],2)), xs[2]), size(xs[2],2), size(xs[2],3))
end
function einsum(::SimpleBinaryRule{('j','l'), ('k','j','l'), ('k','l')}, xs::NTuple{2, Any})
    reshape(_batched_gemm('N', 'T', reshape(xs[1], 1, size(xs[1],1), size(xs[1],2)), xs[2]), size(xs[2],1), size(xs[2],3))
end

# 110
function einsum(::SimpleBinaryRule{('i','j','l'), ('j','l'), ('i','l')}, xs::NTuple{2, Any})
    reshape(_batched_gemm('N', 'N', xs[1], reshape(xs[2], size(xs[2],1), 1, size(xs[2],2))), size(xs[1],1), size(xs[1],3))
end
function einsum(::SimpleBinaryRule{('j','i','l'), ('j','l'), ('i','l')}, xs::NTuple{2, Any})
    reshape(_batched_gemm('T', 'N', xs[1], reshape(xs[2], size(xs[2],1), 1, size(xs[2],2))), size(xs[1],2), size(xs[1],3))
end

# ij,jk->ik : 111
# S = N^2
# T = N^3
for (i1, X1) in enumerate([('i', 'j'), ('j', 'i')])
    for (i2, X2) in enumerate([('j', 'k'), ('k', 'j')])
        for (i3, X3) in enumerate([('i', 'k'), ('k', 'i')])
            A1 = i1==i3 ? :(xs[1]) : :(transpose(xs[1]))
            A2 = i2==i3 ? :(xs[2]) : :(transpose(xs[2]))
            @eval function einsum(::SimpleBinaryRule{$X1,$X2, $X3}, xs::NTuple{2, Any})
                $(i3==1 ? :($A1*$A2) : :($A2*$A1))
            end
            X1B = (X1...,'l')
            X2B = (X2...,'l')
            X3B = (X3...,'l')
            C1 = i1==i3 ? 'N' : 'T'
            C2 = i2==i3 ? 'N' : 'T'
            @eval function einsum(::SimpleBinaryRule{$X1B,$X2B,$X3B}, xs::NTuple{2, Any})
                $(i3==1 ? :(_batched_gemm($C1, $C2, xs[1], xs[2])) : :(_batched_gemm($C2, $C1, xs[2], xs[1])))
            end
        end
    end
end

# there are too many combination in the binary case, so nospecialize
function einsum(::DefaultRule, ixs, iy, @nospecialize(xs::NTuple{2, Any}), size_dict::Dict{LT}) where LT
    @debug "DefaultRule binary" ixs => iy size.(xs)
    ix1, ix2 = ixs
    x1, x2 = xs
    c1, c2, cy, s1, s2, i1, i2, iyb = analyze_binary(_collect(LT,ix1), _collect(LT,ix2), _collect(LT,iy), size_dict)
    rule = SimpleBinaryRule{(i1...,), (i2...,), (iyb...,)}()
    x1 = simplify_unary(_collect(LT,ix1), c1, x1, size_dict)
    x2 = simplify_unary(_collect(LT,ix2), c2, x2, size_dict)
    x1_ = reshape(x1, s1...)
    x2_ = reshape(x2, s2...)
    @debug rule size.((x1_, x2_))
    y_ = reshape(einsum(rule, (x1_, x2_)), [size_dict[x] for x in cy]...)
    return expand_unary(cy, _collect(LT,iy), y_, size_dict)
end

function simplify_unary(ix::Vector{T}, iy::Vector{T}, x, size_dict::Dict{T}) where T
    if ix == iy
        return x
    elseif length(ix) == length(iy) # permutation
        return einsum(Permutedims(), (ix,), iy, (x,), size_dict)
    else
        # diag
        ix_ = unique(ix)
        x_ = length(ix_) != length(ix) ? einsum(Diag(), (ix,), ix_, (x,), size_dict) : x
        # sum
        if length(ix_) != length(iy)
            return einsum(Sum(), (ix_,), iy, (x_,), size_dict)
        elseif ix_ != iy
            return einsum(Permutedims(), (ix_,), iy, (x_,), size_dict)
        else
            return x_
        end
    end
end

function expand_unary(ix::Vector{T}, iy::Vector{T}, x::AbstractArray, size_dict::Dict{T}) where T
    iy_b = unique(iy)
    iy_a = filter(i->i ∈ ix, iy_b)
    y_a = if ix != iy_a
        einsum(Permutedims(), (ix,), iy_a, (x,), size_dict)
    else
        x
    end
    # repeat
    y_b = length(iy_a) != length(iy_b) ? einsum(Repeat(), (iy_a,), iy_b, (y_a,), size_dict) : y_a
    # duplicate
    length(iy_b) != length(iy) ? einsum(Duplicate(), (iy_b,), iy, (y_b,), size_dict) : y_b
end

"""
Get the expected labels.
"""
function analyze_binary(ix1::Vector{T}, ix2::Vector{T}, iy::Vector{T}, size_dict::Dict{T,Int}) where T
    ix_inner, ix1_outer, ix2_outer, batch = _analyze_binary_input(ix1, ix2, iy)
    c1 = vcat(ix1_outer, ix_inner, batch)
    c2 = vcat(ix_inner, ix2_outer, batch)
    cy = vcat(ix1_outer, ix2_outer, batch)
    si = prod(map(x->size_dict[x], ix1_outer))
    sj = prod(map(x->size_dict[x], ix_inner))
    sk = prod(map(x->size_dict[x], ix2_outer))
    sl = prod(map(x->size_dict[x], batch))
    has_i = !isempty(ix1_outer)
    has_j = !isempty(ix_inner)
    has_k = !isempty(ix2_outer)
    has_l = !isempty(batch)
    i1 = Char[]
    i2 = Char[]
    iyb = Char[]
    s1 = Int[]
    s2 = Int[]
    if has_i
        push!(i1, 'i')
        push!(iyb, 'i')
        push!(s1, si)
    end
    if has_j
        push!(i1, 'j')
        push!(i2, 'j')
        push!(s1, sj)
        push!(s2, sj)
    end
    if has_k
        push!(i2, 'k')
        push!(iyb, 'k')
        push!(s2, sk)
    end
    if has_l
        push!(i1, 'l')
        push!(i2, 'l')
        push!(iyb, 'l')
        push!(s1, sl)
        push!(s2, sl)
    end
    return c1, c2, cy, s1, s2, i1, i2, iyb
end

function _analyze_binary_input(ix1::Vector{T}, ix2::Vector{T}, iy::Vector{T}) where T
    ix1_batch = T[]
    ix1_inner = T[]
    ix1_outer = T[]
    for l1 in ix1
        if l1 ∈ ix2
            if l1 ∈ iy  # batch
                l1 ∉ ix1_batch && push!(ix1_batch, l1)
            else        # inner
                l1 ∉ ix1_inner && push!(ix1_inner, l1)
            end
        elseif l1 ∈ iy  # outer
            l1 ∉ ix1_outer && push!(ix1_outer, l1)
        else
            # dangling
        end
    end
    ix2_outer = T[]     # outer dimension of x2
    for l2 in ix2
        if l2 ∉ ix1 && l2 ∈ iy && l2 ∉ ix2_outer
            push!(ix2_outer, l2)
        end
    end
    ix1_inner, ix1_outer, ix2_outer, ix1_batch
end


text:
# LabJackM constants and types

const LJM_StreamReadCallback = Ptr{Cvoid}
const LJM_DeviceReconnectCallback = Ptr{Cvoid}

@cenum IODirection::Int32 begin
    READ = 0
    WRITE = 1
end

@cenum LJMDataType::Int32 begin
    UINT16  = 0
    UINT32  = 1
    INT32   = 2
    FLOAT32 = 3
    BYTE    = 99
    STRING  = 98
end

@cenum LJMStringConstant::Int32 begin
    ALLOCATION_SIZE     = 50 # Buffer size include null-terminator
    MAX_SIZE            = 49 # String length
    MAX_NAME_SIZE       = 256
    MAC_STRING_SIZE     = 18
    IPv4_STRING_SIZE    = 16 # for INE_ADDRSTRLEN
end

const BYTES_PER_REGISTER = 2

"NamesToAddresses uses this when a register name is not found"
const INVALID_NAME_ADDRESS = -1

@cenum LJDeviceType::Int32 begin
    dtANY   = 0
    T4      = 4
    T7      = 7
    DIGIT   = 200
    TSERIES = 84
end

@cenum LJConnectionType::Int32 begin
    ctANY           = 0
    ctANY_TCP       = 0
    USB             = 1
    TCP             = 2
    NETWORK_TCP     = 2
    ETHERNET        = 3
    ETHERNET_TCP    = 3
    WIFI            = 4
    WIFI_TCP        = 4
    NETWORK_UDP     = 4
    ETHERNET_UDP    = 6
    WIFI_UDP        = 7
    NETWORK_ANY     = 8
    ETHERNET_ANY    = 9
    WIFI_ANY        = 10
end

@cenum LJNetworkConstant::Int32 begin
    TCP_PORT            = 502
    ETHERNET_UDP_PORT   = 52362
    WIFI_UDP_PORT       = 502
    NO_IP_ADDRESS       = 0
    NO_PORT             = 0
end

const LIST_ALL_SIZE = 128

@cenum LJTimeouts::Int32 begin
    NO_TIMEOUT                                  = 0
    DEFAULT_USB_SEND_RECEIVE_TIMEOUT_MS         = 2600
    DEFAULT_ETHERNET_OPEN_TIMEOUT_MS            = 1000
    DEFAULT_ETHERNET_SEND_RECEIVE_TIMEOUT_MS    = 2600
    DEFAULT_WIFI_OPEN_TIMEOUT_MS                = 1000
    DEFAULT_WIFI_SEND_RECEIVE_TIMEOUT_MS        = 4000
end

@cenum LJThermocouple::Clong begin
    ttB = 6001
    ttE = 6002
    ttJ = 6003
    ttK = 6004
    ttN = 6005
    ttR = 6006
    ttS = 6007
    ttT = 6008
    ttC = 6009
end


text:
using Simulate, Random, Unitful, Test
import Unitful: Time, ms, s, minute, hr

println(".... testing Simulate.jl .....")
@testset "clock.jl" begin
    include("test_clock.jl")
end

@testset "channel 1 example" begin
    include("test_channels1.jl")
end

@testset "process.jl" begin
    include("test_process.jl")
end

@testset "logger.jl" begin
        include("test_logger.jl")
end


text:
using ArrayViews

type PulseSummaries
    pretrig_mean      ::Vector{Float64}
    pretrig_rms       ::Vector{Float64}
    pulse_average     ::Vector{Float64}
    pulse_rms         ::Vector{Float64}
    rise_time         ::Vector{Float64}
    postpeak_deriv    ::Vector{Float64}
    timestamp         ::Vector{Float64}
    peak_index        ::Vector{Int16}
    peak_value        ::Vector{UInt16}
    min_value         ::Vector{UInt16}

    function PulseSummaries(n::Integer)
        pretrig_mean = Array(Float64, n)
        pretrig_rms = Array(Float64, n)
        pulse_average = Array(Float64, n)
        pulse_rms = Array(Float64, n)
        rise_time = Array(Float64, n)
        postpeak_deriv = Array(Float64, n)
        timestamp = Array(Float64, n)
        peak_index = Array(Int16, n)
        peak_value = Array(UInt16, n)
        min_value = Array(UInt16, n)

        new(pretrig_mean, pretrig_rms, pulse_average, pulse_rms, rise_time,
            postpeak_deriv, timestamp, peak_index, peak_value, min_value)
    end
end

function compute_summary(pulseiterator, Npre, frame_time)
    summary = PulseSummaries(length(pulseiterator))
    Nsamp = length(first(pulseiterator))
    Npost = Nsamp-Npre
    post_peak_deriv_vect = zeros(Float64, Npost)
    for (p, data) in enumerate(pulseiterator)
        # Pretrigger computation first
        s = s2 = 0
        min_idx = 0
        peak_idx = 0
        peak_val = UInt16(0)
        min_idx = 0
        min_val = typemax(UInt16)
        for j = 1:Npre
            d=Int(data[j])
            if d > peak_val
                peak_idx, peak_val = j, d
            elseif d < min_val
                min_idx, min_val = j,d
            end
            s+=d
            s2+=d*d
        end
        ptm = s/Npre
        summary.pretrig_mean[p] = ptm
        summary.pretrig_rms[p] = sqrt(abs(s2/Npre - ptm*ptm))

        # Now post-trigger calculations
        s = s2 = 0
        for j = Npre+1:Nsamp
            d=Int(data[j])
            if d > peak_val
                peak_idx, peak_val = j, d
            elseif d < min_val
                min_idx, min_val = j,d
            end
            # d2=d-ptm
            s+=d
            s2+=d*d
        end

        rise_time::Float64 = estimate_rise_time(data, Npre+1:peak_idx,
                                       peak_val, ptm, frame_time)

        postpeak_data = view(data, peak_idx:Nsamp)
        postpeak_deriv = max_timeseries_deriv_simple(postpeak_data)

        # Copy results into the PulseSummaries object
        summary.pulse_average[p] = s/Npost-ptm
        summary.pulse_rms[p] = sqrt(abs(s2/Npost - ptm*(ptm+2*summary.pulse_average[p])))
        summary.pulse_rms[p] = sqrt(abs(s2/Npost - ptm*(ptm+2*summary.pulse_average[p])))
        summary.rise_time[p] = rise_time
        summary.postpeak_deriv[p] = postpeak_deriv
        summary.peak_index[p] = peak_idx
        summary.min_value[p] = min_val
        if ptm < 0
            summary.peak_value[p] = peak_val
        elseif peak_val > ptm
            summary.peak_value[p] = peak_val - round(UInt16,ptm)
        else
            summary.peak_value[p] = UInt16(0)
        end
    end
    (summary.pretrig_mean, summary.pretrig_rms, summary.pulse_average, summary.pulse_rms, summary.rise_time, summary.postpeak_deriv,
	summary.peak_index, summary.peak_value, summary.min_value)
end

"estimate_rise_time(pulserecord, searchrange, peakval, ptm, frametime)
pulserecord = Vector
searchrange = range to look for risetime in, should go from base of pulse to peak, last(searchrange) is peakindex
peakval = maximum(pulserecord[postrig_region], has not had ptm subtracted off)
ptm = pretrigger mean
framtime = time spacing between points"
function estimate_rise_time(pulserecord, searchrange, peakval,ptm,frametime)
    idx10 = first(searchrange)
    peakindex = last(searchrange)
    (peakindex > length(pulserecord) || length(searchrange)==0) && (return length(pulserecord))

    idx90 = peakindex
    thresh10 = 0.1*(peakval-ptm)+ptm
    thresh90 = 0.9*(peakval-ptm)+ptm
    j=0 # to make j exist after the for loop
    for j = 1:peakindex
        if pulserecord[j] > thresh10
            idx10 = j-1
            break
        end
    end
    for j = j+1:peakindex
        if pulserecord[j] > thresh90
            idx90 = j-1
            break
        end
    end
    rise_time = (idx90-idx10)*frametime
    rise_time
end

# Estimate the derivative (units of arbs / sample) for a pulse record or other timeseries.
# This version uses the default kernel of [-2,-1,0,1,2]/10.0
#
max_timeseries_deriv!(deriv, pulserecord, reject_spikes::Bool) =
    max_timeseries_deriv!(deriv, pulserecord, convert(Vector{eltype(deriv)},[.2 : -.1 : -.2]), reject_spikes)

# Post-peak derivative computed using Savitzky-Golay filter of order 3
# and fitting 1 point before...3 points after.
#
max_timeseries_deriv_SG!(deriv, pulserecord, reject_spikes::Bool) =
    max_timeseries_deriv!(deriv, pulserecord, [-0.11905, .30952, .28572, -.02381, -.45238],
                            reject_spikes)

# Estimate the derivative (units of arbs / sample) for a pulse record or other timeseries.
# Caller pre-allocates the full derivative array, which is available as deriv.
# Returns the maximum value of the derivative.
# The kernel should be a short *convolution* (not correlation) kernel to be convolved
# against the input pulserecord.
# If reject_spikes is true, then the max value at sample i is changed to equal the minimum
# of the values at (i-2, i, i+2). Note that this test only makes sense for kernels of length
# 5 (or less), because only there can it be guaranteed insensitive to unit-length spikes of
# arbitrary amplitude.
#
function max_timeseries_deriv!{T}(
        deriv::Vector{T},       # Modified! Pre-allocate an array of sufficient length
        pulserecord, # The pulse record (presumably starting at the pulse peak)
        kernel::Vector{T},      # The convolution kernel that estimates derivatives
        reject_spikes::Bool  # Whether to employ the spike-rejection test
        )
    N = length(pulserecord)
    Nk = length(kernel)
    @assert length(deriv) >= N+1-Nk
    if Nk > N
        return 0.0
    end
    if Nk+4 > N
        reject_spikes = false
    end
    fill!(deriv, zero(eltype(deriv)))
    for i=1:N-Nk+1
        for j=1:Nk
            deriv[i] += pulserecord[i+Nk-j]*kernel[j] #float
        end
    end
    for i=N-Nk+2:length(deriv)
        deriv[i]=deriv[N-Nk+1]
    end
    if reject_spikes
        for i=3:N-Nk-2
            if deriv[i] > deriv[i+2]
                deriv[i] = deriv[i+2]
            end
            if deriv[i] > deriv[i-2]
                deriv[i] = deriv[i-2]
            end
        end
    end
    maximum(deriv)
end

function max_timeseries_deriv_simple(pulserecord)
    max_deriv = realmin(Float64)
    for j = 1:length(pulserecord)-1
        deriv = float(pulserecord[j+1])-float(pulserecord[j])
        deriv > max_deriv && (max_deriv = deriv)
    end
    max_deriv
end


text:
#reload("Nodes.jl")
#reload("Expressions.jl")
#reload("Parsimonious.jl")
using Base.Test

using Parsimonious
# include("Util.jl")

# BootstrappingGrammarTests
#Tests for the expressions in the grammar that parses the grammar
#definition syntax

# test_quantifier(self):
text = "*"
@show parse(rule_grammar.exprs["quantifier"], text)
@test parse(rule_grammar.exprs["quantifier"], text) == Node("quantifier", text, 1, 1, (Node("", text, 1, 1), Node("_", text, 2, 1)))

# TODO: Whether nodeA == nodeB regardless of "match" to make tests less verbose
# After all, if they match the same text what's the difference?
text = "?"
@test parse(rule_grammar.exprs["quantifier"], text) == Node("quantifier", text, 1, 1, (Node("", text, 1, 1), Node("_", text, 2, 1)))

text = "+"
@test parse(rule_grammar.exprs["quantifier"], text) == Node("quantifier", text, 1, 1, (Node("", text, 1, 1), Node("_", text, 2, 1)))

eq_(a::AnyNode, b::AnyNode) = a == b

# test spaceless_literal(self):
text = dq"anything but quotes#$*&^"
@test eq_(parse(rule_grammar.exprs["spaceless_literal"], text) , Node("spaceless_literal", text, 1, length(text), (Node("", text, 1, length(text)),)))

# TODO: benchmark visitors with array, tuple and vararg paramters
# TODO: generic visit throw error when called for a node for which visit(::MyVisitor, ::Node{:mynode}) exists

# text = r'''r"\""'''
# text = p"""'''r"\""'''"""
# TODO: Bug in spaceless_literal on this text. Or maybe not. At least it probably reports the wrong character for unconsumed text
# TODO: Cool colorized, unicodized underlining of errors like clang
# text = s""" "\\"" """
text = unescape_string(s""" "\\"" """)
@show text
@show length(text)
@show length(unescape_string(text))
println(text)
println(unescape_string(text))
@test eq_(parse(rule_grammar.exprs["spaceless_literal"], text), Node("spaceless_literal", text, 1, 4, ( Node("", text, 1, 4),)))

# test regex(self):
text = p"""~"[a-zA-Z_][a-zA-Z_0-9]*"LI"""

@test eq_(parse(rule_grammar.exprs["regex"], text),
    Node("regex", text, 1, length(text), (
         Node("", text, 1, 1),
         Node("spaceless_literal", text, 2, 25, (
             Node("", text, 2, 25),)),
         Node("", text, 26, 27),
         Node("_", text, 28, 27),)))

ok_(n::MatchNode) = true

# test successes(self):
# Make sure the PEG recognition grammar succeeds on various inputs.
@test ok_(parse(rule_grammar.exprs["label"], "_"))
@test ok_(parse(rule_grammar.exprs["label"], "jeff"))
@test ok_(parse(rule_grammar.exprs["label"], "_THIS_THING"))

@test ok_(parse(rule_grammar.exprs["atom"], s""" "some_label" """))
@test ok_(parse(rule_grammar.exprs["atom"], s""" "some literal" """))
@test ok_(parse(rule_grammar.exprs["atom"], s""" ~"some regex"i """))

@test ok_(parse(rule_grammar.exprs["quantified"], s""" ~"some regex"i* """))
@test ok_(parse(rule_grammar.exprs["quantified"], "thing+"))
@test ok_(parse(rule_grammar.exprs["quantified"], s""" "hi"? """))

@test ok_(parse(rule_grammar.exprs["term"], "this"))
@test ok_(parse(rule_grammar.exprs["term"], "that+"))

@test ok_(parse(rule_grammar.exprs["sequence"], "this that? other"))

@test ok_(parse(rule_grammar.exprs["ored"], s"""this / that+ / "other" """))

# + is higher precedence than &, so 'anded' should match the whole
# thing:
@test ok_(parse(rule_grammar.exprs["lookahead_term"], "&this+"))

@test ok_(parse(rule_grammar.exprs["expression"], "this"))
@test ok_(parse(rule_grammar.exprs["expression"], "this? that other*"))

@test ok_(parse(rule_grammar.exprs["expression"], s""" &this / that+ / "other" """))
@test ok_(parse(rule_grammar.exprs["expression"], s""" this / that? / "other"+ """))
@test ok_(parse(rule_grammar.exprs["expression"], "this? that other*"))

@test ok_(parse(rule_grammar.exprs["rule"], "this = that\r"))
@test ok_(parse(rule_grammar.exprs["rule"], "this = the? that other* \t\r"))
@test ok_(parse(rule_grammar.exprs["rule"], """the=~"hi*"\n"""))

@test ok_(parse(rule_grammar, """
    this = the? that other*
    that = "thing"
    the=~"hi*"
    other = "ahoy hoy"
    """))

"""
class RuleVisitorTests(TestCase):
Tests for ``RuleVisitor``

As I write these, Grammar is not yet fully implemented. Normally, there'd
be no reason to use ``RuleVisitor`` directly.

# test round_trip(self):
Test a simple round trip.

Parse a simple grammar, turn the parse tree into a map of expressions,
and use that to parse another piece of text.

Not everything was implemented yet, but it was a big milestone and a
proof of concept.
"""

import Parsimonious.RuleVisitor
import Parsimonious.DEFUNCT_visit

# TODO: This causes an error because no visit(v, ::Node{:number}, ...) implemented but the error message is messed up.
txt = """number = ~"[0-9]+"\n"""
tree = parse(rule_grammar, txt)
rules, default_rule = DEFUNCT_visit(RuleVisitor(), tree)

text = "98"
@test eq_(parse(default_rule, text), Node("number", text, 1, 2))

# test undefined_rule(self):
# Make sure we throw the right exception on undefined rules.

txt = "boy = howdy\n"
tree = parse(rule_grammar, txt)

# TODO: fails to fail:
# Also fails to remove the LazyReferences
# Also, the references are wrong: Where's "howdy"?
#  ({"boy"=>LazyReference("boy")},LazyReference("boy"))
# @test_throws visit(RuleVisitor(), tree)

@test_throws DEFUNCT_visit(RuleVisitor(), tree)  # ; debug=true)
try
    DEFUNCT_visit(RuleVisitor(), tree)
catch e
    @test isa(e.exc, Parsimonious.UndefinedLabel)
end

# test optional(self):
txt = """boy = "howdy"?\n"""
tree = parse(rule_grammar, txt)
rules, default_rule = DEFUNCT_visit(RuleVisitor(), tree)

howdy = "howdy"

# It should turn into a Node from the Optional and another from the
# Literal within.
@test eq_(parse(default_rule, howdy), Node("boy", howdy, 1, 5, (Node("", howdy, 1, 5),)))

"""
class GrammarTests(TestCase):
Integration-test ``Grammar``: feed it a PEG and see if it works.

# test expressions_from_rules(self):
Test the ``Grammar`` base class's ability to compile an expression
tree from rules.

That the correct ``Expression`` tree is built is already tested in
``RuleGrammarTests``. This tests only that the ``Grammar`` base class's
``_expressions_from_rules`` works.

"""

greeting_grammar = Grammar(s"""greeting = "hi" / "howdy" """)
tree = parse(greeting_grammar, "hi")
@test eq_(tree, Node("greeting", "hi", 1, 2, (Node("", "hi", 1, 2),)))

# test unicode(self):
# Assert that a ``Grammar`` can convert into a string-formatted series of rules.

grammar = Grammar("""
                  bold_text  = bold_open text bold_close
                  text       = ~"[A-Z 0-9]*"is
                  bold_open  = "(("
                  bold_close = "))"
                  """)
lines = split(unicode(grammar), "\n")
function showln(x...)
    show(x...)
    println()
end
[showln(l) for l in lines]

@test lines[1] == "bold_text = bold_open text bold_close"
@show """text = ~"[A-Z 0-9]*"is"""
@show lines
@test in("""text = ~"[A-Z 0-9]*"is""", lines)
@test in(s"""bold_open = "((" """, lines)
@test in(s"""bold_close = "))" """, lines)
@test length(lines) == 4

# test match(self):
# Make sure partial-matching (with pos) works.
grammar = Grammar("""
                  bold_text  = bold_open text bold_close
                  text       = ~"[A-Z 0-9]*"i
                  bold_open  = "(("
                  bold_close = "))"
                  """)
s = " ((boo))yah"
@test eq_(match(grammar, s, 2), Node("bold_text", s, 2, 8, (
                                 Node("bold_open", s, 2, 3),
                                 Node("text", s, 4, 6),
                                 Node("bold_close", s, 7, 8))))

# test bad_grammar(self):
# Constructing a Grammar with bad rules should raise ParseError.

@assert_raises ParseError Grammar("just a bunch of junk")

# test comments(self):
# Test tolerance of comments and blank lines in and around rules.
grammar = Grammar("""# This is a grammar.

                  # It sure is.
                  bold_text  = stars text stars  # nice
                  text       = ~"[A-Z 0-9]*"i #dude

                  stars      = "**"
                  # Pretty good
                  #Oh yeah.#""")  # Make sure a comment doesn't need a
                                  # \n or \r to end.

# TODO: Maybe re-implement show(io::IO, g::Grammar) to do what str and unicode do. But not right now.
# eq_(list(sorted(str(grammar).splitlines())),

@test sort(split(unicode(grammar), "\n")) ==
    ["""bold_text = stars text stars""",
     s"""stars = "**" """,
     """text = ~"[A-Z 0-9]*"i"""]

# test multi_line(self):
# Make sure we tolerate all sorts of crazy line breaks and comments in
# the middle of rules.
grammar = Grammar("""
    bold_text  = bold_open  # commenty comment
                 text  # more comment
                 bold_close
    text       = ~"[A-Z 0-9]*"i
    bold_open  = "((" bold_close =  "))"
    """)
ok_(parse(grammar, "((booyah))"))

# test not(self):
# Make sure "not" predicates get parsed and work properly.
grammar = Grammar(s""" not_arp = !"arp" ~"[a-z]+" """)
@assert_raises ParseError parse(grammar, "arp")
@test ok_(parse(grammar, "argle"))

# test lookahead(self):
grammar = Grammar(s"""starts_with_a = &"a" ~"[a-z]+" """)
@assert_raises ParseError parse(grammar, "burp")

s = "arp"
@test eq_(parse(grammar, "arp"), Node("starts_with_a", s, 1, 3, (
                              Node("", s, 1, 0),
                              Node("", s, 1, 3))))

# test parens(self):
grammar = Grammar("""sequence = "chitty" (" " "bang")+""")
# Make sure it's not as if the parens aren't there:
@assert_raises ParseError parse(grammar, "chitty bangbang")

# test lookahead(self):
grammar = Grammar(s"""starts_with_a = &"a" ~"[a-z]+" """)
@assert_raises ParseError parse(grammar, "burp")

s = "arp"
eq_(parse(grammar, "arp"), Node("starts_with_a", s, 1, 3, (
                              Node("", s, 1, 0),
                              Node("", s, 1, 3))))

# test parens(self):
grammar = Grammar("""sequence = "chitty" (" " "bang")+""")
# Make sure it's not as if the parens aren't there:
@assert_raises ParseError parse(grammar, "chitty bangbang")

s = "chitty bang bang"
# TODO: println(NodeText(parse(grammar, s), s))  -- rather akward syntax
println(string(parse(grammar, s)))
# TODO: give a shit
#@test string(parse(grammar, s)) ==
#"""<Node called 'sequence' matching 'chitty bang bang'>
#    <Node called '' matching 'chitty'>
#    <Node called '' matching ' bang bang'>
#        <Node called '' matching ' bang'>
#            <Node called '' matching ' '>
#            <Node called '' matching 'bang'>
#        <Node called '' matching ' bang'>
#            <Node called '' matching ' '>
#            <Node called '' matching 'bang'>
#"""

# test resolve_refs_order(self):
# Smoke-test a circumstance where lazy references don't get resolved.
grammar = Grammar("""
    expression = "(" terms ")"
    terms = term+
    term = number
    number = ~"[0-9]+"
    """)
parse(grammar, "(34)")

# test infinite_loop(self):
# Smoke-test a grammar that was causing infinite loops while building.
# This was going awry because the "int" rule was never getting marked as
# resolved, so it would just keep trying to resolve it over and over.

Grammar("""
    digits = digit+
    int = digits
    digit = ~"[0-9]"
    number = int
    main = number
    """)

# test right_recursive(self):
# Right-recursive refs should resolve.
grammar = Grammar("""
    digits = digit digits?
    digit = ~"[0-9]"
    """)
ok_(parse(grammar, "12"))

# test badly_circular(self):
# Uselessly circular references should be detected by the grammar
# compiler.
# raise SkipTest('We have yet to make the grammar compiler detect these.')

# TODO: this should at least raise an error because not all of the references have been resolved
grammar = Grammar("""
    foo = bar
    bar = foo
    """)

# test parens_with_leading_whitespace(self):
# Make sure a parenthesized expression is allowed to have leading
# whitespace when nested directly inside another.
parse(Grammar("""foo = ( ("c") )"""), "c")

# test single_quoted_literals(self):
parse(Grammar("""foo = 'a' '"'"""), s"""a" """)


text:
module NewtonsMethod

using ForwardDiff, LinearAlgebra

D(g) = x -> ForwardDiff.derivative(g, x)

function newtonroot(f,fprime, x_0=0.8, tolerance=1.0E-7, maxiter=1000)
    # setup the algorithm
    x_old = x_0
    normdiff = Inf
    iter = 1
    while normdiff > tolerance && iter <= maxiter
        x_new = x_old-f(x_old)/fprime(x_old) # use the passed in map
        normdiff = norm(x_new - x_old)
        x_old = x_new
        iter = iter + 1
    end
    return (x_old, normdiff, iter)
end

newtonroot(f,x_0=0.8, tolerance=1.0E-7, maxiter=1000)=newtonroot(f,D(f), x_0, tolerance, maxiter)

export newtonroot

end # module