text:
export BrickletDistanceIRDistanceCallbackThreshold
struct BrickletDistanceIRDistanceCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletDistanceIRAnalogValueCallbackThreshold
struct BrickletDistanceIRAnalogValueCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletDistanceIRIdentity
struct BrickletDistanceIRIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletDistanceIR
"""
Measures distance up to 150cm with infrared light
"""
mutable struct BrickletDistanceIR <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletDistanceIR(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_distance_ir")
        deviceInternal = package.BrickletDistanceIR(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_distance
"""
    $(SIGNATURES)

Returns the distance measured by the sensor. Possible
distance ranges are 40 to 300, 100 to 800 and 200 to 1500, depending on the
selected IR sensor.

If you want to get the distance periodically, it is recommended to use the
:cb:`Distance` callback and set the period with
:func:`Set Distance Callback Period`.
"""
function get_distance(device::BrickletDistanceIR)
    return device.deviceInternal.get_distance()
end

export get_analog_value
"""
    $(SIGNATURES)

Returns the value as read by a 12-bit analog-to-digital converter.

.. note::
 The value returned by :func:`Get Distance` is averaged over several samples
 to yield less noise, while :func:`Get Analog Value` gives back raw
 unfiltered analog values. The only reason to use :func:`Get Analog Value` is,
 if you need the full resolution of the analog-to-digital converter.

If you want the analog value periodically, it is recommended to use the
:cb:`Analog Value` callback and set the period with
:func:`Set Analog Value Callback Period`.
"""
function get_analog_value(device::BrickletDistanceIR)
    return device.deviceInternal.get_analog_value()
end

export set_sampling_point
"""
    $(SIGNATURES)

Sets a sampling point value to a specific position of the lookup table.
The lookup table comprises 128 equidistant analog values with
corresponding distances.

If you measure a distance of 50cm at the analog value 2048, you
should call this function with (64, 5000). The utilized analog-to-digital
converter has a resolution of 12 bit. With 128 sampling points on the
whole range, this means that every sampling point has a size of 32
analog values. Thus the analog value 2048 has the corresponding sampling
point 64 = 2048/32.

Sampling points are saved on the EEPROM of the Distance IR Bricklet and
loaded again on startup.

.. note::
 An easy way to calibrate the sampling points of the Distance IR Bricklet is
 implemented in the Brick Viewer. If you want to calibrate your Bricklet it is
 highly recommended to use this implementation.
"""
function set_sampling_point(device::BrickletDistanceIR, position, distance)
    device.deviceInternal.set_sampling_point(position, distance)
end

export get_sampling_point
"""
    $(SIGNATURES)

Returns the distance to a sampling point position as set by
:func:`Set Sampling Point`.
"""
function get_sampling_point(device::BrickletDistanceIR, position)
    return device.deviceInternal.get_sampling_point(position)
end

export set_distance_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Distance` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Distance` callback is only triggered if the distance has changed since the
last triggering.
"""
function set_distance_callback_period(device::BrickletDistanceIR, period)
    device.deviceInternal.set_distance_callback_period(period)
end

export get_distance_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Distance Callback Period`.
"""
function get_distance_callback_period(device::BrickletDistanceIR)
    return device.deviceInternal.get_distance_callback_period()
end

export set_analog_value_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Analog Value` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Analog Value` callback is only triggered if the analog value has
changed since the last triggering.
"""
function set_analog_value_callback_period(device::BrickletDistanceIR, period)
    device.deviceInternal.set_analog_value_callback_period(period)
end

export get_analog_value_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Analog Value Callback Period`.
"""
function get_analog_value_callback_period(device::BrickletDistanceIR)
    return device.deviceInternal.get_analog_value_callback_period()
end

export set_distance_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Distance Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the distance is *outside* the min and max values"
 "'i'",    "Callback is triggered when the distance is *inside* the min and max values"
 "'<'",    "Callback is triggered when the distance is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the distance is greater than the min value (max is ignored)"
"""
function set_distance_callback_threshold(device::BrickletDistanceIR, option, min, max)
    device.deviceInternal.set_distance_callback_threshold(option, min, max)
end

export get_distance_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Distance Callback Threshold`.
"""
function get_distance_callback_threshold(device::BrickletDistanceIR)
    return device.deviceInternal.get_distance_callback_threshold()
end

export set_analog_value_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Analog Value Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the analog value is *outside* the min and max values"
 "'i'",    "Callback is triggered when the analog value is *inside* the min and max values"
 "'<'",    "Callback is triggered when the analog value is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the analog value is greater than the min value (max is ignored)"
"""
function set_analog_value_callback_threshold(device::BrickletDistanceIR, option, min, max)
    device.deviceInternal.set_analog_value_callback_threshold(option, min, max)
end

export get_analog_value_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Analog Value Callback Threshold`.
"""
function get_analog_value_callback_threshold(device::BrickletDistanceIR)
    return device.deviceInternal.get_analog_value_callback_threshold()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callbacks

* :cb:`Distance Reached`,
* :cb:`Analog Value Reached`

are triggered, if the thresholds

* :func:`Set Distance Callback Threshold`,
* :func:`Set Analog Value Callback Threshold`

keep being reached.
"""
function set_debounce_period(device::BrickletDistanceIR, debounce)
    device.deviceInternal.set_debounce_period(debounce)
end

export get_debounce_period
"""
    $(SIGNATURES)

Returns the debounce period as set by :func:`Set Debounce Period`.
"""
function get_debounce_period(device::BrickletDistanceIR)
    return device.deviceInternal.get_debounce_period()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletDistanceIR)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletDistanceIR, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
module TestDicts

using Base.Test
using DataFramesMeta
using Compat

y = 3
@compat d = Dict(:s => 3, :y => 44, :d => 5, :e => :(a + b))
@test @with(d, :s + :y) == d[:s] + d[:y]
@test @with(d, :s + y)  == d[:s] + y
@test @with(d, d)  == d
@test @with(d, :s + d[^(:y)])  == d[:s] + d[:y]
@test @with(d, :e.head) == d[:e].head
@test @compat @with(Dict(:s => 3, :y => 44, :d => 5, :e => :(a + b)), :e.head) == d[:e].head

x = @with d begin
    z = y + :y - 1
    :s + z
end
@test x == y + d[:y] - 1 + d[:s]

fun = d -> @with d begin
    z = y + :y - 1
    :s + z
end
@test fun(d) == y + d[:y] - 1 + d[:s]

d2 = @transform(d, z = :y + :s)
@test d2[:z] == 47

d2 = @select(d, :y, z = :y + :s, :e)
@test d2[:z] == 47

end # module

ext: jl


text:
using LinearOperators
using MatrixDepot
using KrylovMethods
using Base.Test

println("=== Testing CLGS ===")

# test with sparse matrix
A   = sprandn(100,10,.2)
Af = (x,flag) -> ((flag=='F') ? A*x : A'*x)
rhs = randn(100)

# test early stopping behaviour
tt = cgls(Af,rhs,tol=1e-20,maxIter=2,out=2,storeInterm=true)
@test tt[2]==-1

# test behaviour for zero rhs
tt = cgls(Af,zeros(size(A,1)),tol=1e-20,maxIter=2,out=2)
@test tt[2]==-9
@test all(tt[1].==0)

xgt = full(A)\rhs
xt  = cgls(LinearOperator(A),rhs,tol=1e-20)
xt2 = cgls(Af,rhs,tol=1e-20,maxIter=100,x=randn(size(xt[1])))
Xt  = cgls(A,rhs,tol=1e-20,maxIter=100,storeInterm=true,out=1)

@test norm(xgt-xt[1])/norm(xgt) < 1e-6
@test norm(xgt-xt2[1])/norm(xgt) < 1e-6
@test norm(Xt[1][:,end]-xt[1])/norm(xt[1]) < 1e-15

# test with dense hilbert matrix
A   = matrixdepot("hilb",10,3)
rhs = randn(size(A,1))
xgt = A\rhs

x1  = cgls(LinearOperator(A),rhs,tol=1e-10,maxIter=10)
x2  = cgls(A,rhs,tol=1e-10,maxIter=10)

@test norm(xgt-x1[1])/norm(xgt) < 1e-6
@test norm(xgt-x2[1])/norm(xgt) < 1e-6
@test norm(x1[1]-x2[1])/norm(x1[1]) < 1e-12

println("=== CGLS: All tests passed. ===")


text:
function Rasters.RasterStack(output::DeterministicOutput)
    Rasters.RasterStack(output.name)
end

abstract type AbstractMixingRatio end
struct MixingRatio{T, N} <: AbstractMixingRatio
    raster::AbstractRaster{T, N}
    relpoints::AbstractArray
end
# MixingRatio(path::String) = MixingRatio(Raster(path, name= :spec001_mr))
function MixingRatio(output::AbstractOutputFile)
    rasterstack = RasterStack(output)
    rellat, rellon = rasterstack[:RELLAT1], rasterstack[:RELLNG1]
    [[lon, lat] for (lon, lat) in zip(rellon, rellat)]
    MixingRatio(rasterstack[:spec001_mr], [[lon, lat] for (lon, lat) in zip(rellat, rellon)])
end
Base.show(io::IO, mime::MIME"text/plain", mr::AbstractMixingRatio) = show(io, mime, mr.raster)

struct CubicMixingRatio
    raster::AbstractRaster
    relpoints::AbstractArray
end


text:
module DtCalc

using Reexport
@reexport using DomTerm
@reexport using FixedPointNumbers
@reexport using Colors
@reexport using FileIO
@reexport using ImageShow
@reexport using ImageTransformations
@reexport using SVGDraw
@reexport using MathMLRepr
@reexport using Mxma

include("display.jl")
include("htmlio.jl")
include("utils.jl")

export
    # basics
    display, make_newline, make_div, make_img,
    # stream io
    hprint, hprintln,
    # utilities
    showIEEE754, ls_imginfo, preview_dir,
    # setings
    typeset_matrix
end


text:
function cddual_shrinking(X::SparseOrFullMat,
	            Y::Vector;
	            C::Real = 1.0,
	            norm::Integer = 2,
	            randomized::Bool = true,
	            maxpasses::Integer = 1000,
	            eps::Real = 0.1)
	const PG_EPS = 1.0e-12

	# l: # of samples
	# n: # of features
	n, l = size(X)
	alpha = zeros(l)
	w = zeros(n)

	# Set U and D
	#  * L1-SVM: U = C, D[i] = 0
	#  * L2-SVM: U = Inf, D[i] = 1 / (2C)
	U = 0.0
	D = Array(Float64, l)
	if norm == 1
		U = C
		for i in 1:l
			D[i] = 0.0
		end
	elseif norm	== 2
		U = Inf
		for i in 1:l
			D[i] = 1.0 / (2.0 * C)
		end
	else
		DomainError("Only L1-SVM and L2-SVM are supported")
	end

	# Set Qbar
	Qbar = Array(Float64, l)
	for i in 1:l
		Qbar[i] = D[i] + vecnorm(X[:, i], 2)^2
	end

	# Loop over examples
	converged = false
	pass = 0

	pg_max_old = Inf
	pg_min_old = -Inf
	active_size = l

	indices = randomized ? randperm(l) : [1:l]

	while !converged
		pass += 1
		if pass > maxpasses
			break
		end

		pg_max_new = -Inf
		pg_min_new = Inf

		if randomized && active_size > 1
			for i = active_size:-1:2
				j = rand(1:i)
				indices[i], indices[j] = indices[j], indices[i]
			end
		end

		k = 1
		while k <= active_size
			decrement_active = false

			i = indices[k]
			g = Y[i] * inner_prod(w, X[:, i]) - 1.0 + D[i] * alpha[i]

			if alpha[i] == 0.0
				if g > pg_max_old
					decrement_active = true
				end
				pg = min(g, 0.0)
			elseif alpha[i] == U
				if g < pg_min_old
					decrement_active = true
				end
				pg = max(g, 0.0)
			else
				pg = g
			end

			pg_max_new = max(pg_max_new, pg)
			pg_min_new = min(pg_min_new, pg)

			if abs(pg) > PG_EPS
				alphabar = alpha[i]
				alpha[i] = min(max(alpha[i] - g / Qbar[i], 0.0), U)
				add_vec!(w, X[:, i], e -> (alpha[i] - alphabar) * Y[i] * e)
			end

			if decrement_active
				indices[k], indices[active_size] = indices[active_size], indices[k]
				active_size -= 1
			else
				k += 1
			end
		end

		if (pg_max_new - pg_min_new) <= eps
			if active_size == l
				converged = true
				break
			else
				active_size = l
				pg_max_old = Inf
				pg_min_old = -Inf
			end
		else
			pg_max_old = pg_max_new > 0.0 ? pg_max_new : Inf
			pg_min_old = pg_min_new < 0.0 ? pg_min_new : -Inf
		end
	end
	return SVMFit(w, pass, converged)
end


text:
#export

## experimental and not exported or documented!
function _adjacency_matrix(data::PopData)
    data_loci = groupby(data.genodata, :locus)
    out_vec = Vector{Array{Int8,2}}(undef, data.metadata.loci)
    for (j,i) in enumerate(data_loci)
        uniq = unique(skipmissing(i.genotype))
        adj_mat = fill(Int8(0), data.metadata.samples, length(uniq))
        for (j,k) in zip(i.genotype, eachrow(adj_mat))
            k .= Ref(j) .=== uniq
        end
        out_vec[j] = adj_mat
    end
    return out_vec
end

"""
    _p_adjust(pvals::Vector{T}, method::String) where T <: Union{Missing, <:AbstractFloat}
Modification to `MultipleTesting.adjust` to include `missing` values in the
returned array. See MultipleTesting.jl docs for full more detailed information.

**Example**
```
julia> _p_adjust([0.1, 0.01, 0.005, 0.3], "bh")

```
### `correction` methods (case insensitive)
- `"bonferroni"` : Bonferroni adjustment
- `"holm"` : Holm adjustment
- `"hochberg"` : Hochberg adjustment
- `"bh"` : Benjamini-Hochberg adjustment
- `"by"` : Benjamini-Yekutieli adjustment
- `"bl"` : Benjamini-Liu adjustment
- `"hommel"` : Hommel adjustment
- `"sidak"` : Šidák adjustment
- `"forwardstop"` or `"fs"` : Forward-Stop adjustment
- `"bc"` : Barber-Candès adjustment
"""
@inline function _p_adjust(pvals::Vector{T}, method::String) where T <: Union{Missing, <:AbstractFloat}
    # make a dict of all possible tests and their respective functions
    d = Dict(
        "bonferroni" => Bonferroni(),
        "holm" => Holm(),
        "hochberg" => Hochberg(),
        "bh" => BenjaminiHochberg(),
        "by" => BenjaminiYekutieli(),
        "bl" => BenjaminiLiu(),
        "hommel" => Hommel(),
        "sidak" => Sidak(),
        "forwardstop" => ForwardStop(),
        "fs" => ForwardStop(),
        "bc" => BarberCandes(),
    )
    p_copy = copy(pvals)
    p_copy[.!ismissing.(p_copy)] .= adjust(p_copy[.!ismissing.(p_copy)] |> Vector{Float64}, d[lowercase(method)])
    return p_copy
end

feature_req() = "\nPlease open an Issue or Pull Request on https://www.github.com/biojulia/PopGen.jl if you would like this feature implemented"


text:
@info("--------- Testing OrthogonalPolynomials ----------")

##
using ACE, Test, ForwardDiff

using LinearAlgebra: norm, cond
using ACE.OrthPolys: OrthPolyBasis
using ACE: evaluate, evaluate_d, evaluate_ed
using ACEbase.Testing: print_tf

##
@info("de-dictionisation")

for ntest = 1:10
   local N, Nquad, tdf, ww, Jd
   N = 8
   Nquad = 1000
   tdf = rand(1000)
   ww = 1.0 .+ rand(1000)
   Jd = OrthPolyBasis(N, 2, 1.0, 1, -1.0, tdf, ww)
   print_tf(@test all(ACE.Testing.test_fio(Jd)))
end
println()

##
@info("Construction and FD vs Grad for randomly generated OrthPolyBasis")

N = 8
Nquad = 1000
tdf = rand(1000)
ww = 1.0 .+ rand(1000)
Jd = OrthPolyBasis(N, 2, 1.0, 2, -1.0, tdf, ww)

let errtol = 1e-12, ntest = 50
   for itest = 1:ntest
      x = 2*rand() - 1
      Jx = evaluate(Jd, x)
      dJx = evaluate_d(Jd, x)
      Jx1, dJx1 = evaluate_ed(Jd, x)
      print_tf(@test all( (Jx, dJx) .≈ (Jx1, dJx1) ) )
      adJx = ForwardDiff.derivative(x -> evaluate(Jd, x), x)
      err = maximum(abs.(dJx - adJx) ./ (1.0 .+ abs.(dJx)))
      print_tf(@test (err < errtol))
   end
   println()
end

##

# TODO: replace this test with an orthogonality test
# @info("Discretised Jacobi are close to the real Jacobi Poly's")
#
# N = 15
# Nquad = 1000
# dt =  2 / Nquad
# tdf = range(-1.0+dt/2, 1.0-dt/2, length=Nquad)
# Jd = OrthPolyBasis(N,  0, 1.0, 0, -1.0, tdf)
# J = Jacobi(0.0, 0.0, N-1, normalise=true)
#
# for ntest = 1:30
#    x = 2*rand() - 1
#    Jdx = evaluate(Jd, x)
#    Jx = evaluate(J, x)
#    Jx /= Jx[1]
#    print_tf((@test norm(Jx - Jdx, Inf) < 10/N^2))
# end
# println()

##

# TODO: convert into a quadrature orthogonality test
# @info("Testing the orthogonality (via A-basis)")
#
# spec = SparseSHIP(3, 10)
# P = ACE.OrthPolys.transformed_jacobi(ACE.maxK(spec)+1, trans, 2.0, 0.5;
#                                       pcut = 2, pin = 2)
# aceB = SHIPBasis(spec, P)
#
# function evalA(aceB, tmp, Rs)
#    Zs = zeros(Int16, length(Rs))
#    ACE.precompute_A!(tmp, aceB, Rs, Zs, 1)
#    return tmp.A[1]
# end
#
# function A_gramian(aceB, Nsamples = 100_000)
#    tmp = ACE.alloc_temp(aceB)
#    lenA = length(tmp.A[1])
#    G = zeros(ComplexF64, lenA, lenA)
#    for n = 1:Nsamples
#       R = ACE.rand_vec(aceB.J)
#       A = evalA(aceB, tmp, [R])
#       for i = 1:lenA, j = 1:lenA
#          G[i,j] +=  A[i] * A[j]'
#       end
#    end
#    return G
# end
#
# G = A_gramian(aceB)
# println(@test cond(G) < 1.2)

##

# ## Quick look at the basis
# using Plots
# N = 5
# Jd = ACE.OrthPolys.discrete_jacobi(N; pcut = 3, pin = 2)
# tp = range(-1, 1, length=100)
# Jp = zeros(length(tp), N)
# for (i,t) in enumerate(tp)
#    Jp[i, :] = evaluate(Jd, t)
# end
# plot(tp, Jp)


text:
module KeplerFlow

	using LinearAlgebra

	include("cCode.jl")
	include("coordinateTransform.jl")
	include("energy.jl")
	include("inputOutput.jl")
	include("interaction.jl")
	include("kahanSummationAlgorithm.jl")
	include("keplerianOrbits.jl")
	include("masses.jl")
	include("simulation.jl")

	export cKeplerFlow!, cartesian2jacobi!, jacobi2cartesian!, calculateEnergy, saveOutput, showOutput, firstInteractionStep!, secondInteractionStep!, kahanSumOneValue, kahanSum!, stumpff!, keplerSolve!, keplerFlow, keplerFlowWithKahan, keplerStep!, keplerStepWithKahan!, calculateMassMu!, simulation

end # module

max_stars_repo_name: salanueva/KeplerFlow

max_issues_repo_name: salanueva/KeplerFlow

max_forks_repo_name: salanueva/KeplerFlow


text:
using BernsteinExpansions, Compat, BenchmarkTools

SUITE = BenchmarkGroup()  # parent BenchmarkGroup to contain our suite

SUITE["Univariate"] = SUni = BenchmarkGroup()
SUni["Univariate monomial, floating point input, degree 3"] = @benchmarkable univariate(3, 3, 1.0, 2.0)
SUni["Univariate monomial, floating point input, degree 5"] = @benchmarkable univariate(1, 5, 0.0, 1.0)

SUITE["Multivariate"] = SMulti = BenchmarkGroup()
k = [1,2,3,4,5,2,3,4]
l = [5,5,5,5,5,5,5,5]

low = [1//1,1//1,1//1,1//1,1//1,1//1,1//1,1//1]
high = [2//1,2//1,2//1,2//1,2//1,2//1,2//1,2//1]
m = multivariate(k, l, low, high)
SMulti["Rational multivariate implicit expansion"] = @benchmarkable multivariate($k, $l, $low, $high)
SMulti["Rational multivariate full expansion (tensor algorithm)"] = @benchmarkable generate_tensor_form($m, algorithm="tensor_generated")
SMulti["Rational multivariate full expansion (kron algorithm)"] = @benchmarkable generate_tensor_form($m, algorithm="linear_kron")

low = convert(Vector{Float64}, low)
high = convert(Vector{Float64}, high)
m = multivariate(k, l, low, high)
SMulti["Float64 multivariate implicit expansion"] = @benchmarkable multivariate($k, $l, $low, $high)
SMulti["Float64 multivariate full expansion (tensor algorithm)"] = @benchmarkable generate_tensor_form($m, algorithm="tensor_generated")
SMulti["Float64 multivariate full expansion (kron algorithm)"] = @benchmarkable generate_tensor_form($m, algorithm="linear_kron")

max_stars_repo_name: JuliaReach/BernsteinExpansions

max_issues_repo_name: JuliaReach/BernsteinExpansions

max_forks_repo_name: JuliaReach/BernsteinExpansions


text:
import Reachability.ReachSets.discretize

# ===================================================================
# Discretization of a continuous-time system without input (ZeroSet)
# ===================================================================
A = sparse([1, 1, 2, 3, 4], [1, 2, 2, 4, 3], [1., 2., 3., 4., 5.], 4, 4)
X0 = BallInf(zeros(4), 0.1)
#U = ZeroSet(size(A, 1))
cont_sys_homog = ContinuousSystem(A, X0)
δ = 0.01

# no bloating, do not use Pade approximation
discr_sys_homog = discretize(cont_sys_homog, δ, approx_model="nobloating", pade_expm=false)
@test inputdim(discr_sys_homog) == 0 # there is no input set!

# no bloating, use Pade approximation
discr_sys_homog = discretize(cont_sys_homog, δ, approx_model="nobloating", pade_expm=true)
@test inputdim(discr_sys_homog) == 0

# bloating, do not use Pade approximation
discr_sys_homog = discretize(cont_sys_homog, δ, pade_expm=false)
@test inputdim(discr_sys_homog) == 0

# bloating, first order
discr_sys_homog = discretize(cont_sys_homog, δ, approx_model="firstorder")
@test inputdim(discr_sys_homog) == 0

# ===============================================================
# Discretization of a continuous-time system with constant input
# ===============================================================
U = Ball2(ones(4), 0.5)
cont_sys = ContinuousSystem(A, X0, U)

# no bloating, do not use Pade approximation
discr_sys = discretize(cont_sys, δ, approx_model="nobloating", pade_expm=false)
@test inputdim(discr_sys) == 4
#@test length(discr_sys.U) == 1 # a constant input has no length method!
inputs = next_set(inputset(discr_sys)) # getter: inputset(discr_sys) = discr_sys.s.U
@test dim(inputs) == 4
@test isa(inputs, LinearMap)
@test isa(inputs.X, Ball2) && inputs.X.center == ones(4) && inputs.X.radius == 0.5

# no bloating, use Pade approximation
discr_sys = discretize(cont_sys, δ, approx_model="nobloating", pade_expm=true)
@test inputdim(discr_sys) == 4

# bloating, do not use Pade approximation
discr_sys = discretize(cont_sys, δ, pade_expm=false)
@test inputdim(discr_sys) == 4

#@test length(discr_sys.U) == 1
inputs = next_set(inputset(discr_sys))
@test dim(inputs) == 4
@test isa(inputs, MinkowskiSum)

# bloating, use Pade approximation
discr_sys = discretize(cont_sys, δ, pade_expm=true)
@test inputdim(discr_sys) == 4

discr_sys = discretize(cont_sys, δ, approx_model="firstorder")
@test inputdim(discr_sys) == 4

# ===================================================================
# Discretization of a continuous-time system with time-varying input
# ===================================================================
Ui = [Ball2(0.01*i*ones(4), i*0.2) for i in 1:3]
cont_sys = ContinuousSystem(A, X0, Ui)

# no bloating, do not use Pade approximation
discr_sys = discretize(cont_sys, δ, approx_model="nobloating", pade_expm=false)
Ui_d = inputset(discr_sys)
@test length(Ui_d) == 3

input_state = start(Ui_d)
inputs, input_state = next(Ui_d, input_state)
@test dim(inputs) == 4
@test isa(inputs, LinearMap)
@test isa(inputs.X, Ball2) && inputs.X.center == 0.01*ones(4) && inputs.X.radius == 0.2

inputs, input_state = next(Ui_d, input_state)
@test dim(inputs) == 4
@test isa(inputs, LinearMap)
@test isa(inputs.X, Ball2) && inputs.X.center == 0.01*2*ones(4) && inputs.X.radius == 0.2*2

inputs, input_state = next(Ui_d, input_state)
@test dim(inputs) == 4
@test isa(inputs, LinearMap)
@test isa(inputs.X, Ball2) && inputs.X.center == 0.01*3*ones(4) && inputs.X.radius == 0.2*3

# no bloating, use Pade approximation
discr_sys = discretize(cont_sys, δ, approx_model="nobloating", pade_expm=true)

# bloating, do not use Pade approximation
discr_sys = discretize(cont_sys, δ, pade_expm=false)
Ui_d = inputset(discr_sys)
@test length(Ui_d) == 3

input_state = start(Ui_d)
inputs, input_state = next(Ui_d, input_state)
@test dim(inputs) == 4
@test isa(inputs, MinkowskiSum)
X = inputs.X.X
@test isa(X, Ball2) && X.center == 0.01*ones(4) && X.radius == 0.2

inputs, input_state = next(Ui_d, input_state)
@test dim(inputs) == 4
@test isa(inputs, MinkowskiSum)
X = inputs.X.X
@test isa(X, Ball2) && X.center == 0.01*2*ones(4) && X.radius == 0.2*2

inputs, input_state = next(Ui_d, input_state)
@test dim(inputs) == 4
@test isa(inputs, MinkowskiSum)
X = inputs.X.X
@test isa(X, Ball2) && X.center == 0.01*3*ones(4) && X.radius == 0.2*3

# bloating, use Pade approximation
discr_sys = discretize(cont_sys, δ, pade_expm=true)

# lazy symmetric interval hull
discr_sys = discretize(cont_sys, δ, lazy_sih=true)
discr_sys = discretize(cont_sys, δ, lazy_sih=false)


text:
using Pkg
using SafeTestsets
using Test

function dev_subpkg(subpkg)
    subpkg_path = joinpath(dirname(@__DIR__), "lib", subpkg)
    Pkg.develop(PackageSpec(path=subpkg_path))
end

dev_subpkg("IntegralsCuba")
dev_subpkg("IntegralsCubature")

@time @safetestset "Interface Tests" begin include("interface_tests.jl") end
@time @safetestset "Derivative Tests" begin include("derivative_tests.jl") end
@time @safetestset "Infinite Integral Tests" begin include("inf_integral_tests.jl") end


text:
"""
    setentry(logger, name::String, type; kwargs...)

Used to add a new entry/field or modify an existing entry/field in the `logger`.

# Adding a new entry
Specify a unique `name` to add a new entry to the logger. The `type` is used to provide
reasonable formatting defaults and must be included. The keyword arguments control
the behavior/formatting of the field:
* `fmt` A C-style format string used to control the format of the field
* `index` Column for the entry in the output. Negative numbers insert from the end.
* `level` Verbosity level for the entry. A higher level will be printed less often.
        Level 0 will always be printed, unless the logger is disabled. Prefer to use
        a minimum level of 1.
* `width` Width of the column. Data is left-aligned to this width.
* `ccrayon` A [`ConditionalCrayon`](@ref) for conditional formatting.

# Modified an existing entry
This method can also modify an existing entry, if `name` is already a field int the logger.
The `type` can be omitted in this case. Simply specify any of the keyword arguments
with the new setting.
"""
function setentry(log::Logger, name::String, type::Type{T}=Float64;
        fmt::String=default_format(log, name, T),
        index::Integer=default_index(log, name),
        level=default_level(log, name),
        width::Integer=default_width(log, name, fmt),
        ccrayon=nothing
    ) where T
    @assert width > 0 "Field width must be positive"

    # Check if index is valid
    newentry = !haskey(log.fmt, name)
    index == 0 && error("Index can't be zero. Must be positive or negative")
    if length(log.data) == 0 && abs(index) == 1
        index = 1
    elseif index < 0
        if index >= -length(log.data)
            index = length(log.data) + index + 1 + newentry  # add one to the end
        else
            error("Invalid index. Negative indices must be greater than $(-length(log.data))")
        end
    elseif index > length(log.data)
        error("Invalid index. Must be less than $(length(log.data))")
    end

    # Check if the field already exists
    if haskey(log.fmt, name)
        espec = log.fmt[name]
        fid = espec.uid
        oldindex = log.idx[fid]

        # Shift data
        if oldindex != index
            shiftswap!(log.data, index, oldindex)
            shiftswap!(log.crayons, index, oldindex)
            shiftidx!(log.idx,  fid, index)
        end

        if isnothing(ccrayon)
            ccrayon = espec.ccrayon
        end

        if fmt != espec.fmt || level != espec.level || width != espec.width ||
            ccrayon != espec.ccrayon
            log.fmt[name] = EntrySpec(espec.type, fmt, fid, level, width, ccrayon)
        end
    else
        @assert type != Nothing "Must specify type for a new field"

        # Insert new field
        fid = length(log.idx) + 1
        insert!(log.data, index, " "^width)
        insert!(log.crayons, index, Crayon(reset=true))
        push!(log.idx, fid)
        shiftidx!(log.idx, fid, index)

        if isnothing(ccrayon)
            ccrayon = ConditionalCrayon()
        end

        # Set field format and index
        log.fmt[name] = EntrySpec(T, fmt, fid, level, width, ccrayon)
    end

    # Add formatter if it doesn't exist
    if !haskey(log.fmtfun, fmt)
        log.fmtfun[fmt] = generate_formatter(fmt)
    end
    return
end

function default_index(log::Logger, name::String)
    if haskey(log.fmt, name)
        uid = log.fmt[name].uid
        return log.idx[uid]::Int16
    end
    Int16(-1)
end

function default_level(log::Logger, name::String)
    if haskey(log.fmt, name)
        return log.fmt[name].level::UInt8
    end
    UInt8(DEFAULT_LEVEL)
end

function default_width(log::Logger, name::String, fmt::String)::UInt8
    if haskey(log.fmt, name)
        width = log.fmt[name].width::UInt8
    else
        width = getwidth(fmt)
        if width <= 0
            width = DEFAULT_WIDTH
        else
            width += 1
        end
    end
    return UInt8(max(length(name) + 1, width))
end

# Default Formatting
function default_format(log::Logger, name::String, ::Type{T}) where T
    if haskey(log.fmt, name)
        return log.fmt[name].fmt
    end
    _getformat(log.defaults, T)
end
default_format(log::Logger, ::Type{T}) where T = _getformat(log.defaults, T)

"""
    set_default_format(logger, type, fmt)

Set the default format for entries type type is a sub-type of `type`. For example:

    set_default_format(logger, AbstractFloat, "%0.2d")

Will print all floating point fields with 2 decimal places. The most format for the
type closest to the default is always chosen, such that if we specified

    set_default_format(logger, Float64, "%0.1d")

This would override the previous behavior for `Float64`s.

"""
function set_default_format(log::Logger, ::Type{T}, fmt::String) where T
    log.defaults[T] = fmt
end

function _getformat(fmt::Dict, ::Type{T}) where T
    if haskey(fmt, T)
        return fmt[T]
    else
        return _newdefault(fmt, T)
    end
end

function _newdefault(fmt::Dict, ::Type{T}) where T
    Tsuper = Any
    for (k,v) in pairs(fmt)
        if (T <: k) && (k <: Tsuper)
            Tsuper = k
        end
    end
    fmt[T] = fmt[Tsuper]
end


text:
export ObjectAgent_DEP, calcDEPForce, jac!

# ---------------------------- DEP -------------------------
const k_DEP = 3.47538687303371e-23; # pi*obj.e_0*obj.Medium.Permittivity*obj.Particle.Radius^3
const f_CM_R = -0.461847136009005;
const f_CM_I = -0.145447673066018;
const DEP_ampl = 20;

mutable struct ObjectAgent_DEP{T,U} <: ObjectAgent{T,U}
    name::String;
    pos::Tuple{T, T, T};
    Fdes::Tuple{T,T,T};
    actList::Array{Tuple{U, U}, 1};
    act_used::Array{U, 1};
    Γ::Array{T, 2};
    Λₐ::Tuple{Array{T, 2}, Array{T, 2}, Array{T, 2}}
    neighbors::Array{Tuple{ObjectAgent{T,U}, Array{U, 1}, Array{U, 1}},1};
    xk::Vector{T};
    zk::Vector{T};
    uk::Vector{T};
    rcvBuffer_x::Vector{T};
    rcvBuffer_u::Vector{T};
    rcvBuffer_N::Vector{T};
    J::Array{T, 2};
    lsd::LinSystemData{T};
    decDir::Vector{T};
    fvk::Vector{T};
    Cr::Tuple{Array{T, 2}, Array{T, 2}, Array{T, 2}};
    Ci::Tuple{Array{T, 2}, Array{T, 2}, Array{T, 2}};
    Cr_red::Tuple{Array{T, 2}, Array{T, 2}, Array{T, 2}};
    Ci_red::Tuple{Array{T, 2}, Array{T, 2}, Array{T, 2}};
    Fdes_sc::T;
    tmpA1::Array{T, 2};
    tmpA2::Array{T, 2};
    A_kkt::Array{T, 2};
    b_kkt::Vector{T};
    a_r::Vector{T};
    a_i::Vector{T};
    b_r::Vector{T};
    b_i::Vector{T};

    function ObjectAgent_DEP(name::String,
        pos::Tuple{T, T, T},
        Fdes::Tuple{T,T,T},
        aa::ActuatorArray{T,U},
        actList::Array{Tuple{U, U}, 1},
        act_used::Array{U, 1}) where {T<:Real, U<:Unsigned}

        neighbrs = Array{Tuple{ObjectAgent{T,U}, Array{U, 1}, Array{U, 1}},1}[];

        # Compute Γ and Λₐ matrices matrix
        Γ, Λₐ = genC_DEP(aa, pos, actList);

        # Normalize the optimization problem
        Fdes_sc = norm(Fdes);
        Fdes = Fdes./Fdes_sc;
        Γ /= Fdes_sc;
        Γ *= DEP_ampl^2;

        # Initialize LinSystemData = preallocate data for solving the optimization problem
        lsd = LinSystemData{T}(length(act_used), 3);

        ####################### Will be removed when optimized #################
        Cr = (zeros(T, length(actList), length(actList)), zeros(T, length(actList), length(actList)), zeros(T, length(actList), length(actList)));
        Ci = (zeros(T, length(actList), length(actList)), zeros(T, length(actList), length(actList)), zeros(T, length(actList), length(actList)));
        Cr_red = (zeros(T, length(act_used), length(actList)), zeros(T, length(act_used), length(actList)), zeros(T, length(act_used), length(actList)));
        Ci_red = (zeros(T, length(act_used), length(actList)), zeros(T, length(act_used), length(actList)), zeros(T, length(act_used), length(actList)));
        for i in 1:3
            Cr[i] .= k_DEP*f_CM_R*(Λₐ[i]*Γ' + Γ*Λₐ[i]');
            Ci[i] .= k_DEP*f_CM_I*(Λₐ[i]*Γ' - Γ*Λₐ[i]');
            Cr_red[i] .= Cr[i][act_used,:];
            Ci_red[i] .= Ci[i][act_used,:];
        end
        ########################################################################
        b_kkt = zeros(T, length(act_used)+4);
        b_kkt[end-3] = 1;

        # number of actuators
        n = length(actList);

        new{T,U}(name, pos, Fdes, actList, act_used, Γ, Λₐ, neighbrs,
        zeros(T, n), # xk
        zeros(T, n), # zk
        zeros(T, n), # uk
        zeros(T, n), # buffer xk
        zeros(T, n), # buffer uk
        zeros(T, n), # buffer N
        zeros(T, length(act_used), 3), # J
        lsd,         # LinSystemData
        lsd.x_sol,   # decDir
        zeros(T, 3), # fvk
        Cr,
        Ci,
        Cr_red,
        Ci_red,
        Fdes_sc,
        zeros(T, length(act_used), length(actList)), # tmpA1
        zeros(T, length(act_used), length(actList)),  # tmpA2
        zeros(T, length(act_used)+4, length(act_used)+4), # A_kkt
        b_kkt, # b_kkt
        zeros(T, 3), # a_r
        zeros(T, 3), # a_i
        zeros(T, 3), # b_r
        zeros(T, 3), # b_i
        );
    end
end

include("DEP_model.jl");
function genC_DEP(aa::ActuatorArray{T,U}, pos::Tuple{T, T, T},
    al::Array{Tuple{U, U}, 1} = vec([(i,j) for j in oneunit(aa.ny):aa.ny, i in oneunit(aa.nx):aa.nx]) ) where {T<:Real,U<:Unsigned}
    # p = Vector{ComplexF32}(undef, length(al));
    Γ = zeros(T, length(al), 3);
    Λₐ = (zeros(T, length(al), 3), # Lambda_x
                zeros(T, length(al), 3), # Lambda_y
                zeros(T, length(al), 3)); # Lambda_z

    # precomp_w = 1e-4.*[0.560343281319084;0.641589559368632;0.533388218135986;0.966133522019187;0.729809440489818;1.3223366991602;0.584185556810487;0.818917874471294;1.11114468489008;0.713304002034397;1.22069212004038;0.953479444255371;1.08520203870701;1.46115156961446;1.09120813884737;0.865971706264724;1.28472318044273;1.23782780533269;1.72344388564653;1.45405798718119;2.08577707140132;2.70136063519805;2.63772101370926;1.32738453095412;0.644683064696835];
    # precomp_l = 1e-4.*[0.55876527575103;0.645337660700216;0.945329383178398;0.534095933916202;0.726695701399419;0.630721025019882;1.25654560374921;0.825210823499908;0.776820179903889;1.10664391236893;0.950397116802921;1.03393240360361;1.07713933735852;1.1808595051766;1.40115863622489;1.34097784885123;1.32872852691204;1.47701744598014;1.44687111683753;1.71675951300292;2.02665619561559;1.32765680200996;0.644584541216948;2.78286596946413;2.70273226082876];
    # precomp_h = [0.100739487195064;0.140522401433007;0.0506445451505626;0.0667646496786449;0.159161476319683;0.0702663244060299;0.0590559460583131;0.0892352746802388;0.0470985196175195;0.0560667423873221;0.0301643706052842;0.0203694962578092;0.0110378724835092;0.0269432413334298;0.0218821277880894;0.0290242270512687;-0.0207383124040105;0.0218872886332642;0.0140689493436172;0.0130489329027799;0.0055462869079359;0.0314610233484279;-0.0368002751182956;0.0303826109876248;-0.0349416987084138];

    # Iterate over all the electrodes
    for (k, a_pos) in enumerate(al)
        # position of the k-th actuator
        xa, ya = actuatorPosition(aa, a_pos[1:2]);

        # x = (pos[1] - xa);
        # y = (pos[2] - ya);
        # z = pos[3];
        #
        # for i in 1:length(precomp_w)
        #     dVx, dVy, dVz, dVxx, dVxy, dVxz, dVyy, dVyz, dVzz = dV_opt( x/(precomp_w[i]/2), y/(precomp_l[i]/2), z);
        #
        #     Γ[k,1] -= precomp_h[i]*dVx/(precomp_w[i]/2);
        #     Γ[k,2] -= precomp_h[i]*dVy/(precomp_l[i]/2);
        #     Γ[k,3] -= precomp_h[i]*dVz;
        #
        #     # Lambda_x
        #     Λₐ[1][k,1] -= precomp_h[i]*dVxx/(precomp_w[i]/2)^2;
        #     Λₐ[1][k,2] -= precomp_h[i]*dVxy/(precomp_w[i]/2)/(precomp_l[i]/2);
        #     Λₐ[1][k,3] -= precomp_h[i]*dVxz/(precomp_w[i]/2);
        #
        #     # Lambda_y
        #     Λₐ[2][k,2] -= precomp_h[i]*dVyy/(precomp_l[i]/2)^2;
        #     Λₐ[2][k,3] -= precomp_h[i]*dVyz/(precomp_l[i]/2);
        #
        #     # Lambda_z
        #     Λₐ[3][k,3] -= precomp_h[i]*dVzz;
        # end

        sz = precomp_x[1,2]-precomp_x[1,1];

        x = (pos[1] - xa)/(sz/2);
        y = (pos[2] - ya)/(sz/2);
        z = (pos[3])/(sz/2);

        for i in 1:size(precomp_x)[1], j in 1:size(precomp_x)[2]
            dVx, dVy, dVz, dVxx, dVxy, dVxz, dVyy, dVyz, dVzz = dV_opt( x - precomp_x[i,j]/(sz/2), y - precomp_y[i,j]/(sz/2), z);

            Γ[k,1] -= precomp_z[i,j]*dVx/(sz/2);
            Γ[k,2] -= precomp_z[i,j]*dVy/(sz/2);
            Γ[k,3] -= precomp_z[i,j]*dVz/(sz/2);

            # Lambda_x
            Λₐ[1][k,1] -= precomp_z[i,j]*dVxx/(sz/2)^2;
            Λₐ[1][k,2] -= precomp_z[i,j]*dVxy/(sz/2)^2;
            Λₐ[1][k,3] -= precomp_z[i,j]*dVxz/(sz/2)^2;

            # Lambda_y
            Λₐ[2][k,2] -= precomp_z[i,j]*dVyy/(sz/2)^2;
            Λₐ[2][k,3] -= precomp_z[i,j]*dVyz/(sz/2)^2;

            # Lambda_x
            Λₐ[3][k,3] -= precomp_z[i,j]*dVzz/(sz/2)^2;
        end

        Λₐ[2][k,1] = Λₐ[1][k,2];
        Λₐ[3][k,1] = Λₐ[1][k,3];
        Λₐ[3][k,2] = Λₐ[2][k,3];

    end

    return Γ, Λₐ;
end

function calcDEPForce(aa::ActuatorArray{T,U}, pos::Tuple{T,T,T}, phases::Array{T,2}) where {T<:Real, U<:Unsigned}
    phsvec = vec(phases');
    ur = DEP_ampl*cos.(phsvec);
    ui = DEP_ampl*sin.(phsvec);
    ur[isnan.(real(ur))] .= 0;
    ui[isnan.(real(ui))] .= 0;
    Γ, Λₐ = genC_DEP(aa, pos);

    F_DEP = (
                2*k_DEP*f_CM_R*ur'*Λₐ[1]*Γ'*ur + 2*k_DEP*f_CM_R*ui'*Λₐ[1]*Γ'*ui + 2*k_DEP*f_CM_I*ur'*(Γ*Λₐ[1]' - Λₐ[1]*Γ')*ui,
                2*k_DEP*f_CM_R*ur'*Λₐ[2]*Γ'*ur + 2*k_DEP*f_CM_R*ui'*Λₐ[2]*Γ'*ui + 2*k_DEP*f_CM_I*ur'*(Γ*Λₐ[2]' - Λₐ[2]*Γ')*ui,
                2*k_DEP*f_CM_R*ur'*Λₐ[3]*Γ'*ur + 2*k_DEP*f_CM_R*ui'*Λₐ[3]*Γ'*ui + 2*k_DEP*f_CM_I*ur'*(Γ*Λₐ[3]' - Λₐ[3]*Γ')*ui
            );

    return F_DEP;
end

function updatex_dirFixed!(oa::ObjectAgent_DEP{T,U}, λ::T) where {T<:Real, U<:Unsigned}
    oa.xk .= oa.zk .- oa.uk;

    vk_r = cos.(oa.xk);
    vk_i = sin.(oa.xk);

    jac!(oa); # DEP ~40 us

    Fd = [oa.Fdes[1] oa.Fdes[2] oa.Fdes[3]];

    n_red = length(oa.act_used);
    oa.A_kkt .=  vcat([(1/λ)*I zeros(n_red,1) -oa.J],
                 [zeros(1,n_red) 1 Fd],
                 [-oa.J' Fd' zeros(3,3)]);

    for i in 1:3
        oa.b_kkt[end-3+i] = vk_r'*oa.Cr[i]*vk_r + vk_i'*oa.Cr[i]*vk_i - 2*vk_r'*oa.Ci[i]*vk_i;
    end

    sol = oa.A_kkt\oa.b_kkt;
    oa.xk[oa.act_used] .+= sol[1:end-4];

    nothing;
end

function costFun(oa::ObjectAgent_DEP{T, U}, xk::Vector{T}) where {T<:Real,U<:Unsigned}
    xk_r = cos.(xk);
    xk_i = sin.(xk);

    cost = 0;

    mul!(oa.b_r, oa.Γ', xk_r);
    mul!(oa.b_i, oa.Γ', xk_i);
    for k in 1:3
        mul!(oa.a_r, oa.Λₐ[k]', xk_r);
        mul!(oa.a_i, oa.Λₐ[k]', xk_i);
        cost += (2*k_DEP*f_CM_R*dot(oa.a_r,oa.b_r) + 2*k_DEP*f_CM_R*dot(oa.a_i,oa.b_i) - 2*k_DEP*f_CM_I*(dot(oa.a_r,oa.b_i) - dot(oa.b_r,oa.a_i)) - oa.Fdes[k])^2;
    end

    return cost
end

function costFun!(oa::ObjectAgent_DEP)
    vk_r = cos.(oa.xk);
    vk_i = sin.(oa.xk);
    # This function is equivalent to the following lines of code:
    # for i in 1:3
        # oa.fvk[i] = vk_r'*oa.Cr[i]*vk_r + vk_i'*oa.Cr[i]*vk_i - 2*vk_r'*oa.Ci[i]*vk_i - oa.Fdes[i];
    # end

    mul!(oa.b_r, oa.Γ', vk_r);
    mul!(oa.b_i, oa.Γ', vk_i);
    for k in 1:3
        mul!(oa.a_r, oa.Λₐ[k]', vk_r);
        mul!(oa.a_i, oa.Λₐ[k]', vk_i);
        oa.fvk[k] = 2*k_DEP*f_CM_R*dot(oa.a_r,oa.b_r) + 2*k_DEP*f_CM_R*dot(oa.a_i,oa.b_i) - 2*k_DEP*f_CM_I*(dot(oa.a_r,oa.b_i) - dot(oa.b_r,oa.a_i)) - oa.Fdes[k];
    end

    nothing
end

function jac!(oa::ObjectAgent_DEP)
    # This function is equivalent to the following lines of code:
    # for i in 1:3
        # oa.J[:,i] .= 2*(dur*oa.Ci_red[i] - dui*oa.Cr_red[i])*vk_r .+ 2*(dur*oa.Cr_red[i] .+ dui*oa.Ci_red[i])*vk_i;
    # end

    n = length(oa.act_used);
    m = length(oa.actList);

    vk_r = cos.(oa.xk);
    vk_i = sin.(oa.xk);

    @inbounds for k in 1:3
        for i in 1:n
            oa.J[i,k] = 0;
            @simd for j in 1:m
                oa.J[i,k] += 2*(vk_r[oa.act_used[i]]*oa.Ci_red[k][i,j] - vk_i[oa.act_used[i]]*oa.Cr_red[k][i,j])*vk_r[j];
                oa.J[i,k] += 2*(vk_r[oa.act_used[i]]*oa.Cr_red[k][i,j] + vk_i[oa.act_used[i]]*oa.Ci_red[k][i,j])*vk_i[j];
            end
        end
    end

    nothing;
end


text:
function greeting(user::AppUser)
    println("Hello [$(user.lastname)]")
end

function greeting(user::AppUser, date::Date)
    println("Hello[$(user.lastname)] on[$date]")
end

function function1()
    println("In function1")
end


text:
"""
    tfill(v, ::Val{D}) where D

Returns a tuple of length `D` that contains `D` times the object `v`.
In contrast to `tuple(fill(v,D)...)` which returns the same result, this function is type-stable.
"""
function tfill(v, ::Val{D}) where D
  t = tfill(v, Val{D-1}())
  (v,t...)
end

tfill(v,::Val{0}) = ()
tfill(v,::Val{1}) = (v,)
tfill(v,::Val{2}) = (v,v)
tfill(v,::Val{3}) = (v,v,v)

"""
    get_val_parameter(::Val{T}) where T
    get_val_parameter(::Type{Val{T}}) where T

Returns `T`.
"""
function get_val_parameter(::Val{T}) where T
  T
end

function get_val_parameter(::Type{Val{T}}) where T
  T
end


text:
using RuntimeGeneratedFunctions, BenchmarkTools
using Serialization
using Test

RuntimeGeneratedFunctions.init(@__MODULE__)

function f(_du,_u,_p,_t)
    @inbounds _du[1] = _u[1]
    @inbounds _du[2] = _u[2]
    nothing
end

ex1 = :((_du,_u,_p,_t) -> begin
    @inbounds _du[1] = _u[1]
    @inbounds _du[2] = _u[2]
    nothing
end)

ex2 = :(function f(_du,_u,_p,_t)
    @inbounds _du[1] = _u[1]
    @inbounds _du[2] = _u[2]
    nothing
end)

ex3 = :(function (_du::T,_u::Vector{E},_p::P,_t::Any) where {T<:Vector,E,P}
    @inbounds _du[1] = _u[1]
    @inbounds _du[2] = _u[2]
    nothing
end)

f1 = @RuntimeGeneratedFunction(ex1)
f2 = @RuntimeGeneratedFunction(ex2)
f3 = @RuntimeGeneratedFunction(ex3)

@test f1 isa Function

du = rand(2)
u = rand(2)
p = nothing
t = nothing

@test f1(du,u,p,t) === nothing
du == u
du = rand(2)
f2(du,u,p,t)
@test du == u
du = rand(2)
@test f3(du,u,p,t) === nothing
du == u

t1 = @belapsed $f($du,$u,$p,$t)
t2 = @belapsed $f1($du,$u,$p,$t)
t3 = @belapsed $f2($du,$u,$p,$t)
t4 = @belapsed $f3($du,$u,$p,$t)

@test t1 ≈ t2 atol = 3e-9
@test t1 ≈ t3 atol = 3e-9
@test t1 ≈ t4 atol = 3e-9

function no_worldage()
    ex = :(function f(_du,_u,_p,_t)
        @inbounds _du[1] = _u[1]
        @inbounds _du[2] = _u[2]
        nothing
    end)
    f1 = @RuntimeGeneratedFunction(ex)
    du = rand(2)
    u = rand(2)
    p = nothing
    t = nothing
    f1(du,u,p,t)
end
@test no_worldage() === nothing

# Test show()
@test sprint(show, MIME"text/plain"(), @RuntimeGeneratedFunction(Base.remove_linenums!(:((x,y)->x+y+1)))) ==
     """
     RuntimeGeneratedFunction(#=in $(@__MODULE__)=#, #=using $(@__MODULE__)=#, :((x, y)->begin
               x + y + 1
           end))"""

# Test with precompilation
push!(LOAD_PATH, joinpath(@__DIR__, "precomp"))
using RGFPrecompTest

@test RGFPrecompTest.f(1,2) == 3
@test RGFPrecompTest.g(40) == 42

# Test that RuntimeGeneratedFunction with identical body expressions (but
# allocated separately) don't clobber each other when one is GC'd.
f_gc = @RuntimeGeneratedFunction(Base.remove_linenums!(:((x,y)->x+y+100001)))
let
    @RuntimeGeneratedFunction(Base.remove_linenums!(:((x,y)->x+y+100001)))
end
GC.gc()
@test f_gc(1,-1) == 100001

# Test that threaded use works
tasks = []
for k=1:4
    let k=k
        t = Threads.@spawn begin
            r = Bool[]
            for i=1:100
                f = @RuntimeGeneratedFunction(Base.remove_linenums!(:((x,y)->x+y+$i*$k)))
                x = 1; y = 2;
                push!(r, f(x,y) == x + y + i*k)
            end
            r
        end
        push!(tasks, t)
    end
end
@test all(all.(fetch.(tasks)))

# Test that globals are resolved within the correct scope

module GlobalsTest
    using RuntimeGeneratedFunctions
    RuntimeGeneratedFunctions.init(@__MODULE__)

    y_in_GlobalsTest = 40
    f = @RuntimeGeneratedFunction(:(x->x + y_in_GlobalsTest))
end

@test GlobalsTest.f(2) == 42

f_outside = @RuntimeGeneratedFunction(GlobalsTest, :(x->x + y_in_GlobalsTest))
@test f_outside(2) == 42

@test_throws ErrorException @eval(module NotInitTest
    using RuntimeGeneratedFunctions
    # RuntimeGeneratedFunctions.init(@__MODULE__) # <-- missing
    f = @RuntimeGeneratedFunction(:(x->x+y))
end)

# closures
if VERSION >= v"1.7.0-DEV.351"
    ex = :(x -> (y -> x + y))
    @test @RuntimeGeneratedFunction(ex)(2)(3) === 5

    ex = :(x -> (f(y::Int)::Float64 = x + y; f))
    @test @RuntimeGeneratedFunction(ex)(2)(3) === 5.0

    ex = :(x -> function (y::Int)
                    return x + y
                end)
    @test @RuntimeGeneratedFunction(ex)(2)(3) === 5

    ex = :(x -> function f(y::Int)::UInt8
                    return x + y
                end)
    @test @RuntimeGeneratedFunction(ex)(2)(3) === 0x05

    ex = :(x -> sum(i^2 for i in 1:x))
    @test @RuntimeGeneratedFunction(ex)(3) === 14

    ex = :(x -> [2i for i in 1:x])
    @test @RuntimeGeneratedFunction(ex)(3) == [2, 4, 6]
end

# Serialization

buf = IOBuffer(read(`$(Base.julia_cmd()) "serialize_rgf.jl"`))
deserialized_f = deserialize(buf)
@test deserialized_f(11) == "Hi from a separate process. x=11"


text:
__precompile__(false)

module ProjectEuler

import Glob

function __init__()
    withenv("PROJECT_EULER_WITHOUT_RESULT"=>true) do
        for file in Glob.glob(["..",r"\d\..*",r".*\.jl"],@__DIR__)
            include(file)
        end
    end
end

end

max_stars_repo_name: miRoox/ProjectEuler

max_issues_repo_name: miRoox/ProjectEuler

max_forks_repo_name: miRoox/ProjectEuler


text:
using Test
using DiskArrays: eachchunk, haschunks, Chunked, GridChunks, readblock!
import ArchGDAL as AG

@testset "test_array.jl" begin
    @testset "RasterDataset Type" begin
        AG.readraster("ospy/data4/aster.img") do ds
            @testset "Test forwarded methods" begin
                @test ds isa AG.RasterDataset{UInt8}
                @test AG.getgeotransform(ds) ==
                      [419976.5, 15.0, 0.0, 4.6624225e6, 0.0, -15.0]
                @test AG.nraster(ds) == 3
                @test AG.getband(ds, 1) isa AG.AbstractRasterBand
                @test startswith(AG.getproj(ds), "PROJCS")
                @test AG.width(ds) == 5665
                @test AG.height(ds) == 5033
                @test AG.getdriver(ds) isa AG.Driver
                @test splitpath(AG.filelist(ds)[1]) ==
                      ["ospy", "data4", "aster.img"]
                @test splitpath(AG.filelist(ds)[2]) ==
                      ["ospy", "data4", "aster.rrd"]
                @test AG.listcapability(ds) isa Dict
                @test AG.ngcp(ds) == 0
                @test AG.write(ds, tempname()) == nothing
                @test AG.testcapability(ds, "ODsCCreateLayer") == false
            end
            @testset "DiskArray chunk interface" begin
                b = AG.getband(ds, 1)
                @test eachchunk(ds) == GridChunks(size(ds), (64, 64, 1))
                @test eachchunk(b) == GridChunks(size(b), (64, 64))
                @test haschunks(ds) == Chunked()
                @test haschunks(b) == Chunked()
            end
            @testset "Reading into non-arrays" begin
                data1 = view(zeros(3, 3, 3), 1:3, 1:3, 1:3)
                readblock!(ds, data1, 1:3, 1:3, 1:3)
                @test data1 == ds[1:3, 1:3, 1:3]
            end
        end
    end

    @testset "Test Array getindex" begin
        AG.readraster("ospy/data4/aster.img") do ds
            @testset "Dataset indexing" begin
                @testset "dims dropped correctly" begin
                    @test typeof(ds[:, :, :]) <: Array{UInt8,3}
                    @test typeof(ds[:, :, 1]) <: Array{UInt8,2}
                    @test typeof(ds[:, 1, 1]) <: Array{UInt8,1}
                    @test typeof(ds[1, 1, 1]) <: UInt8
                end
                @testset "range indexing" begin
                    buffer = ds[1:AG.width(ds), 1:AG.height(ds), 1:1]
                    @test typeof(buffer) <: Array{UInt8,3}
                    total = sum(buffer)
                    count = sum(buffer .> 0)
                    @test total / count ≈ 76.33891347095299
                    @test total / (AG.height(ds) * AG.width(ds)) ≈
                          47.55674749653172
                end
                @testset "colon indexing" begin
                    buffer = ds[:, :, 1]
                    total = sum(buffer)
                    count = sum(buffer .> 0)
                    @test total / count ≈ 76.33891347095299
                    @test total / (AG.height(ds) * AG.width(ds)) ≈
                          47.55674749653172
                end
                @testset "int indexing" begin
                    @test ds[755, 2107, 1] == 0xff
                end
            end
            @testset "RasterBand indexing" begin
                band = AG.getband(ds, 1)
                @testset "dims dropped correctly" begin
                    @test typeof(band[:, :]) <: Array{UInt8,2}
                    @test typeof(band[:, 1]) <: Array{UInt8,1}
                    @test typeof(band[1, 1]) <: UInt8
                end
                @testset "range indexing" begin
                    buffer = band[1:AG.width(band), 1:AG.height(band)]
                    @test typeof(buffer) <: Array{UInt8,2}
                    total = sum(buffer)
                    count = sum(buffer .> 0)
                    @test total / count ≈ 76.33891347095299
                    @test total / (AG.height(band) * AG.width(band)) ≈
                          47.55674749653172
                end
                @testset "colon indexing" begin
                    buffer = band[:, :]
                    total = sum(buffer)
                    count = sum(buffer .> 0)
                    @test total / count ≈ 76.33891347095299
                    @test total / (AG.height(band) * AG.width(band)) ≈
                          47.55674749653172
                end
                @testset "int indexing" begin
                    @test band[755, 2107] == 0xff
                end
            end
        end
    end

    cp("ospy/data4/aster.img", "ospy/data4/aster_write.img"; force = true)

    @testset "Test Array setindex" begin
        AG.readraster("ospy/data4/aster_write.img"; flags = AG.OF_UPDATE) do ds
            @testset "Dataset setindex" begin
                @test ds[755, 2107, 1] == 0xff
                ds[755, 2107, 1] = 0x00
                @test ds[755, 2107, 1] == 0x00
                ds[755:755, 2107:2107, 1:1] = reshape([0x01], 1, 1, 1)
                @test ds[755, 2107, 1] == 0x01
                ds[755:755, 2107:2107, 1] = reshape([0x02], 1, 1)
                @test ds[755, 2107, 1] == 0x02
                ds[755:755, 2107, 1] = [0x03]
                @test ds[755, 2107, 1] == 0x03
            end
            @testset "RasterBand setindex" begin
                band = AG.getband(ds, 1)
                band[755:755, 2107] = [0x00]
                @test band[755, 2107] == 0x00
                band[755, 2107] = 0x01
                @test band[755, 2107] == 0x01
                band[755:755, 2107:2107] = reshape([0xff], 1, 1)
                @test band[755, 2107] == 0xff
                buffer = band[:, :]
                band[:, :] = buffer .* 0x00
                @test sum(band[:, :]) == 0x00
                band[:, 1:500] = buffer[:, 1:500]
                band[:, 501:end] = buffer[:, 501:end]
                @test sum(buffer) / sum(buffer .> 0) ≈ 76.33891347095299
                @test_throws DimensionMismatch band[:, 501:end] = [1, 2, 3]
            end
        end
    end

    @testset "Test Array constructor" begin
        AG.readraster("ospy/data4/aster_write.img"; flags = AG.OF_UPDATE) do ds
            @test sprint(print, ds) == """
            GDAL Dataset (Driver: HFA/Erdas Imagine Images (.img))
            File(s):
              ospy/data4/aster_write.img

            Dataset (width x height): 5665 x 5033 (pixels)
            Number of raster bands: 3
              [GA_Update] Band 1 (Undefined): 5665 x 5033 (UInt8)
              [GA_Update] Band 2 (Undefined): 5665 x 5033 (UInt8)
              [GA_Update] Band 3 (Undefined): 5665 x 5033 (UInt8)
            """
            buffer = Array(ds)
            typeof(buffer) <: Array{UInt8,3}
            total = sum(buffer[:, :, 1:1])
            count = sum(buffer[:, :, 1:1] .> 0)
            @test total / count ≈ 76.33891347095299
            @test total / (AG.height(ds) * AG.width(ds)) ≈ 47.55674749653172

            AG.copy(ds) do copy1
                @test typeof(AG.copywholeraster!(ds, copy1)) == typeof(copy1)
                @test typeof(
                    AG.copywholeraster!(AG.RasterDataset(copy1), ds),
                ) == typeof(ds)
                @test typeof(AG.copywholeraster!(copy1, ds)) == typeof(ds)
            end
        end

        AG.create(
            AG.getdriver("MEM"),
            width = 2,
            height = 2,
            nbands = 0,
            dtype = UInt8,
        ) do dataset
            @test_throws Union{ArgumentError,MethodError} AG.RasterDataset(
                dataset,
            )
            @test_throws DimensionMismatch AG._common_size(dataset)
        end
    end
end


text:
"""
AMQP functionality for Genie.
"""
module AMQP
using AMQP
import Genie, Logging

function getconnection(port::Int = Genie.config.AMQP_PORT,  user::String = Genie.config.AMQP_USER, password::String = Genie.config.AMQP_PASSWORD, virtualhost::String = Genie.config.AMQP_VIRTUALHOST, host::String = Genie.config.AMQP_HOST)
    try
        auth_params = Dict{String,Any}("MECHANISM"=>"AMQPLAIN", "LOGIN"=>user, "PASSWORD"=>password)
        conn = Connection(;virtualhost=virtualhost, host=host,port=port,auth_params=auth_params)
        @info conn
    catch(e)
        @error "Failed parsing `server` parameter info."
        @error ex
    end
end


text:
include("const.jl")
include("dir_tree.jl")
include("get_files.jl")
include("NCI60_data.jl")


text:
using DelimitedFiles: readdlm
using ProgressBars

function load_edges(filename, delim, comments)
    return readdlm(filename, delim, Int, comments=comments)
end

function kcore(nodes::Vector{Int64}, edges::Matrix{Int64})
    if minimum(edges) == 0
        edges = edges .+ 1
        nodes = Vector(1:length(nodes)+1)
    end
    n = length(nodes)
    i = n
    c = 0
    # Compute degrees
    degrees = zeros(Int64, n)
    for edge in eachrow(edges)
        degrees[edge[1]] += 1
        degrees[edge[2]] += 1
    end

    for i in ProgressBar(1:n)
        index = argmin(degrees)
        node = nodes[index]
        c = max(c, degrees[index])
        # Update graph
        deleteat!(nodes, index)

        neighbor_edges_indices = Vector{Int64}()

        for (k, edge) in enumerate(eachrow(edges))
            if edge[1] == node || edge[2] == node
                push!(neighbor_edges_indices, k)
            end
        end

        edges = edges[findall(x->!(x in neighbor_edges_indices), 1:size(edges, 1)),:]

        degrees = zeros(Int64, length(nodes))
        for edge in eachrow(edges)
            index1 = findall(x->x==edge[1], nodes)[1]
            index2 = findall(x->x==edge[2], nodes)[1]
            degrees[index1] += 1
            degrees[index2] += 1
        end
        i = i - 1
    end

    return c
end

filename = "com-amazon.ungraph.txt"
edges = load_edges(filename, '\t', true)
nodes = Vector(1:maximum(edges))
core_value = kcore(nodes, edges)

max_stars_repo_name: quentin-burthier/GrandsGraphes

max_issues_repo_name: quentin-burthier/GrandsGraphes

max_forks_repo_name: quentin-burthier/GrandsGraphes


text:
const StringList = Ptr{Cstring}

"""
The value of `GDALDataType` could be different from `GDAL.GDALDataType`.

It maps correctly to `GDAL.GDALDataType` if you do e.g.

```jldoctest; output = false
convert(GDAL.GDALDataType, ArchGDAL.GDT_Unknown)

# output
GDT_Unknown::GDALDataType = 0x00000000
```
"""
@enum(
    GDALDataType,
    GDT_Unknown = 0,
    GDT_Byte = 1,
    GDT_UInt16 = 2,
    GDT_Int16 = 3,
    GDT_UInt32 = 4,
    GDT_Int32 = 5,
    GDT_Float32 = 6,
    GDT_Float64 = 7,
    GDT_CInt16 = 8,
    GDT_CInt32 = 9,
    GDT_CFloat32 = 10,
    GDT_CFloat64 = 11,
    GDT_TypeCount = 12,
)

"""
The value of `OGRFieldType` could be different from `GDAL.OGRFieldType`.

It maps correctly to `GDAL.OGRFieldType` if you do e.g.

```jldoctest; output = false
convert(GDAL.OGRFieldType, ArchGDAL.OFTInteger)

# output
OFTInteger::OGRFieldType = 0x00000000
```
"""
@enum(
    OGRFieldType,
    OFTInteger = 0,
    OFTIntegerList = 1,
    OFTReal = 2,
    OFTRealList = 3,
    OFTString = 4,
    OFTStringList = 5,
    OFTWideString = 6,
    OFTWideStringList = 7,
    OFTBinary = 8,
    OFTDate = 9,
    OFTTime = 10,
    OFTDateTime = 11,
    OFTInteger64 = 12,
    OFTInteger64List = 13,
    OFTMaxType = 14, # 13
)

"""
The value of `OGRFieldSubType` could be different from `GDAL.OGRFieldSubType`.

It maps correctly to `GDAL.OGRFieldSubType` if you do e.g.

```jldoctest; output = false
convert(GDAL.OGRFieldSubType, ArchGDAL.OFSTNone)

# output
OFSTNone::OGRFieldSubType = 0x00000000
```
"""
@enum(
    OGRFieldSubType,
    OFSTNone = 0,
    OFSTBoolean = 1,
    OFSTInt16 = 2,
    OFSTFloat32 = 3,
    OFSTJSON = 4,
    OFSTMaxSubType = 5, # 4
)

"""
The value of `OGRJustification` could be different from `GDAL.OGRJustification`.

It maps correctly to `GDAL.OGRJustification` if you do e.g.

```jldoctest; output = false
convert(GDAL.OGRJustification, ArchGDAL.OJUndefined)

# output
OJUndefined::OGRJustification = 0x00000000
```
"""
@enum(OGRJustification, OJUndefined = 0, OJLeft = 1, OJRight = 2,)

"""
The value of `GDALRATFieldType` could be different from `GDAL.GDALRATFieldType`.

It maps correctly to `GDAL.GDALRATFieldType` if you do e.g.

```jldoctest; output = false
convert(GDAL.GDALRATFieldType, ArchGDAL.GFT_Integer)

# output
GFT_Integer::GDALRATFieldType = 0x00000000
```
"""
@enum(GDALRATFieldType, GFT_Integer = 0, GFT_Real = 1, GFT_String = 2,)

"""
The value of `GDALRATFieldUsage` could be different from `GDAL.GDALRATFieldUsage`.

It maps correctly to `GDAL.GDALRATFieldUsage` if you do e.g.

```jldoctest; output = false
convert(GDAL.GDALRATFieldUsage, ArchGDAL.GFU_Generic)

# output
GFU_Generic::GDALRATFieldUsage = 0x00000000
```
"""
@enum(
    GDALRATFieldUsage,
    GFU_Generic = 0,
    GFU_PixelCount = 1,
    GFU_Name = 2,
    GFU_Min = 3,
    GFU_Max = 4,
    GFU_MinMax = 5,
    GFU_Red = 6,
    GFU_Green = 7,
    GFU_Blue = 8,
    GFU_Alpha = 9,
    GFU_RedMin = 10,
    GFU_GreenMin = 11,
    GFU_BlueMin = 12,
    GFU_AlphaMin = 13,
    GFU_RedMax = 14,
    GFU_GreenMax = 15,
    GFU_BlueMax = 16,
    GFU_AlphaMax = 17,
    GFU_MaxCount = 18,
)

"""
The value of `GDALAccess` could be different from `GDAL.GDALAccess`.

It maps correctly to `GDAL.GDALAccess` if you do e.g.

```jldoctest; output = false
convert(GDAL.GDALAccess, ArchGDAL.GA_ReadOnly)

# output
GA_ReadOnly::GDALAccess = 0x00000000
```
"""
@enum(GDALAccess, GA_ReadOnly = 0, GA_Update = 1,)

"""
The value of `GDALRWFlag` could be different from `GDAL.GDALRWFlag`.

It maps correctly to `GDAL.GDALRWFlag` if you do e.g.

```jldoctest; output = false
convert(GDAL.GDALRWFlag, ArchGDAL.GF_Read)

# output
GF_Read::GDALRWFlag = 0x00000000
```
"""
@enum(GDALRWFlag, GF_Read = 0, GF_Write = 1,)

"""
The value of `GDALPaletteInterp` could be different from `GDAL.GDALPaletteInterp`.

It maps correctly to `GDAL.GDALPaletteInterp` if you do e.g.

```jldoctest; output = false
convert(GDAL.GDALPaletteInterp, ArchGDAL.GPI_Gray)

# output
GPI_Gray::GDALPaletteInterp = 0x00000000
```
"""
@enum(GDALPaletteInterp, GPI_Gray = 0, GPI_RGB = 1, GPI_CMYK = 2, GPI_HLS = 3,)

"""
The value of `GDALColorInterp` could be different from `GDAL.GDALColorInterp`.

It maps correctly to `GDAL.GDALColorInterp` if you do e.g.

```jldoctest; output = false
convert(GDAL.GDALColorInterp, ArchGDAL.GCI_Undefined)

# output
GCI_Undefined::GDALColorInterp = 0x00000000
```
"""
@enum(
    GDALColorInterp,
    GCI_Undefined = 0,
    GCI_GrayIndex = 1,      # GreyScale
    GCI_PaletteIndex = 2,   # Paletted (see associated color table)
    GCI_RedBand = 3,        # Red band of RGBA image
    GCI_GreenBand = 4,      # Green band of RGBA image
    GCI_BlueBand = 5,       # Blue band of RGBA image
    GCI_AlphaBand = 6,      # Alpha (0=transparent, 255=opaque)
    GCI_HueBand = 7,        # Hue band of HLS image
    GCI_SaturationBand = 8, # Saturation band of HLS image
    GCI_LightnessBand = 9,  # Lightness band of HLS image
    GCI_CyanBand = 10,      # Cyan band of CMYK image
    GCI_MagentaBand = 11,   # Magenta band of CMYK image
    GCI_YellowBand = 12,    # Yellow band of CMYK image
    GCI_BlackBand = 13,     # Black band of CMYK image
    GCI_YCbCr_YBand = 14,   # Y Luminance
    GCI_YCbCr_CbBand = 15,  # Cb Chroma
    GCI_YCbCr_CrBand = 16,  # Cr Chroma
    GCI_Max = 17,           # Max current value = 16
)

"""
The value of `GDALAsyncStatusType` could be different from `GDAL.GDALAsyncStatusType`.

It maps correctly to `GDAL.GDALAsyncStatusType` if you do e.g.

```jldoctest; output = false
convert(GDAL.GDALAsyncStatusType, ArchGDAL.GARIO_PENDING)

# output
GARIO_PENDING::GDALAsyncStatusType = 0x00000000
```
"""
@enum(
    GDALAsyncStatusType,
    GARIO_PENDING = 0,
    GARIO_UPDATE = 1,
    GARIO_ERROR = 2,
    GARIO_COMPLETE = 3,
    GARIO_TypeCount = 4,
)

"""
The value of `OGRSTClassId` could be different from `GDAL.OGRSTClassId`.

It maps correctly to `GDAL.OGRSTClassId` if you do e.g.

```jldoctest; output = false
convert(GDAL.OGRSTClassId, ArchGDAL.OGRSTCNone)

# output
OGRSTCNone::ogr_style_tool_class_id = 0x00000000
```
"""
@enum(
    OGRSTClassId,
    OGRSTCNone = 0,
    OGRSTCPen = 1,
    OGRSTCBrush = 2,
    OGRSTCSymbol = 3,
    OGRSTCLabel = 4,
    OGRSTCVector = 5,
)

"""
The value of `OGRSTUnitId` could be different from `GDAL.OGRSTUnitId`.

It maps correctly to `GDAL.OGRSTUnitId` if you do e.g.

```jldoctest; output = false
convert(GDAL.OGRSTUnitId, ArchGDAL.OGRSTUGround)

# output
OGRSTUGround::ogr_style_tool_units_id = 0x00000000
```
"""
@enum(
    OGRSTUnitId,
    OGRSTUGround = 0,
    OGRSTUPixel = 1,
    OGRSTUPoints = 2,
    OGRSTUMM = 3,
    OGRSTUCM = 4,
    OGRSTUInches = 5,
)

"""
The value of `OGRwkbGeometryType` could be different from `GDAL.OGRwkbGeometryType`.

It maps correctly to `GDAL.OGRwkbGeometryType` if you do e.g.

```jldoctest; output = false
convert(GDAL.OGRwkbGeometryType, ArchGDAL.wkbUnknown)

# output
wkbUnknown::OGRwkbGeometryType = 0x00000000
```
"""
@enum(
    OGRwkbGeometryType,
    wkbUnknown = 0,
    wkbPoint = 1,
    wkbLineString = 2,
    wkbPolygon = 3,
    wkbMultiPoint = 4,
    wkbMultiLineString = 5,
    wkbMultiPolygon = 6,
    wkbGeometryCollection = 7,
    wkbCircularString = 8,
    wkbCompoundCurve = 9,
    wkbCurvePolygon = 10,
    wkbMultiCurve = 11,
    wkbMultiSurface = 12,
    wkbCurve = 13,
    wkbSurface = 14,
    wkbPolyhedralSurface = 15,
    wkbTIN = 16,
    wkbTriangle = 17,
    wkbNone = 18,
    wkbLinearRing = 19,
    wkbCircularStringZ = 20,
    wkbCompoundCurveZ = 21,
    wkbCurvePolygonZ = 22,
    wkbMultiCurveZ = 23,
    wkbMultiSurfaceZ = 24,
    wkbCurveZ = 25,
    wkbSurfaceZ = 26,
    wkbPolyhedralSurfaceZ = 27,
    wkbTINZ = 28,
    wkbTriangleZ = 29,
    wkbPointM = 30,
    wkbLineStringM = 31,
    wkbPolygonM = 32,
    wkbMultiPointM = 33,
    wkbMultiLineStringM = 34,
    wkbMultiPolygonM = 35,
    wkbGeometryCollectionM = 36,
    wkbCircularStringM = 37,
    wkbCompoundCurveM = 38,
    wkbCurvePolygonM = 39,
    wkbMultiCurveM = 40,
    wkbMultiSurfaceM = 41,
    wkbCurveM = 42,
    wkbSurfaceM = 43,
    wkbPolyhedralSurfaceM = 44,
    wkbTINM = 45,
    wkbTriangleM = 46,
    wkbPointZM = 47,
    wkbLineStringZM = 48,
    wkbPolygonZM = 49,
    wkbMultiPointZM = 50,
    wkbMultiLineStringZM = 51,
    wkbMultiPolygonZM = 52,
    wkbGeometryCollectionZM = 53,
    wkbCircularStringZM = 54,
    wkbCompoundCurveZM = 55,
    wkbCurvePolygonZM = 56,
    wkbMultiCurveZM = 57,
    wkbMultiSurfaceZM = 58,
    wkbCurveZM = 59,
    wkbSurfaceZM = 60,
    wkbPolyhedralSurfaceZM = 61,
    wkbTINZM = 62,
    wkbTriangleZM = 63,
    wkbPoint25D = 64,
    wkbLineString25D = 65,
    wkbPolygon25D = 66,
    wkbMultiPoint25D = 67,
    wkbMultiLineString25D = 68,
    wkbMultiPolygon25D = 69,
    wkbGeometryCollection25D = 70,
)

"""
The value of `OGRwkbByteOrder` could be different from `GDAL.OGRwkbByteOrder`.

It maps correctly to `GDAL.OGRwkbByteOrder` if you do e.g.

```jldoctest; output = false
convert(GDAL.OGRwkbByteOrder, ArchGDAL.wkbXDR)

# output
wkbXDR::OGRwkbByteOrder = 0x00000000
```
"""
@enum(OGRwkbByteOrder, wkbXDR = 0, wkbNDR = 1,)

@enum(
    GDALOpenFlag,
    OF_READONLY = GDAL.GDAL_OF_READONLY,                # 0x00
    OF_UPDATE = GDAL.GDAL_OF_UPDATE,                  # 0x01
    # OF_All                  = GDAL.GDAL_OF_ALL,                   # 0x00
    OF_RASTER = GDAL.GDAL_OF_RASTER,                  # 0x02
    OF_VECTOR = GDAL.GDAL_OF_VECTOR,                  # 0x04
    OF_GNM = GDAL.GDAL_OF_GNM,                     # 0x08
    OF_KIND_MASK = GDAL.GDAL_OF_KIND_MASK,               # 0x1e
    OF_SHARED = GDAL.GDAL_OF_SHARED,                  # 0x20
    OF_VERBOSE_ERROR = GDAL.GDAL_OF_VERBOSE_ERROR,           # 0x40
    OF_INTERNAL = GDAL.GDAL_OF_INTERNAL,                # 0x80
    # OF_DEFAULT_BLOCK_ACCESS = GDAL.GDAL_OF_DEFAULT_BLOCK_ACCESS,  # 0
    OF_ARRAY_BLOCK_ACCESS = GDAL.GDAL_OF_ARRAY_BLOCK_ACCESS,      # 0x0100
    OF_HASHSET_BLOCK_ACCESS = GDAL.GDAL_OF_HASHSET_BLOCK_ACCESS,    # 0x0200
    # OF_RESERVED_1           = GDAL.GDAL_OF_RESERVED_1,            # 0x0300
    OF_BLOCK_ACCESS_MASK = GDAL.GDAL_OF_BLOCK_ACCESS_MASK,       # 0x0300
)

@enum(
    FieldValidation,
    F_VAL_NULL = GDAL.OGR_F_VAL_NULL,                      # 0x0001
    F_VAL_GEOM_TYPE = GDAL.OGR_F_VAL_GEOM_TYPE,                 # 0x0002
    F_VAL_WIDTH = GDAL.OGR_F_VAL_WIDTH,                     # 0x0004
    F_VAL_ALLOW_NULL_WHEN_DEFAULT = GDAL.OGR_F_VAL_ALLOW_NULL_WHEN_DEFAULT,   # 0x0008
    F_VAL_ALLOW_DIFFERENT_GEOM_DIM = GDAL.OGR_F_VAL_ALLOW_DIFFERENT_GEOM_DIM,  # 0x0010
)


text:
module Mathematica

# This package is based on MathLink
import MathLink
import MathLink: WExpr, WSymbol
import IterTools: takewhile

export weval, @importsymbol

include("types.jl")
include("operators.jl")
include("show.jl")
include("eval.jl")

end # module


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

#using Compat

# begin enum curandStatus
typealias curandStatus UInt32
const CURAND_STATUS_SUCCESS = (UInt32)(0)
const CURAND_STATUS_VERSION_MISMATCH = (UInt32)(100)
const CURAND_STATUS_NOT_INITIALIZED = (UInt32)(101)
const CURAND_STATUS_ALLOCATION_FAILED = (UInt32)(102)
const CURAND_STATUS_TYPE_ERROR = (UInt32)(103)
const CURAND_STATUS_OUT_OF_RANGE = (UInt32)(104)
const CURAND_STATUS_LENGTH_NOT_MULTIPLE = (UInt32)(105)
const CURAND_STATUS_DOUBLE_PRECISION_REQUIRED = (UInt32)(106)
const CURAND_STATUS_LAUNCH_FAILURE = (UInt32)(201)
const CURAND_STATUS_PREEXISTING_FAILURE = (UInt32)(202)
const CURAND_STATUS_INITIALIZATION_FAILED = (UInt32)(203)
const CURAND_STATUS_ARCH_MISMATCH = (UInt32)(204)
const CURAND_STATUS_INTERNAL_ERROR = (UInt32)(999)
# end enum curandStatus

typealias curandStatus_t curandStatus

# begin enum curandRngType
typealias curandRngType UInt32
const CURAND_RNG_TEST = (UInt32)(0)
const CURAND_RNG_PSEUDO_DEFAULT = (UInt32)(100)
const CURAND_RNG_PSEUDO_XORWOW = (UInt32)(101)
const CURAND_RNG_PSEUDO_MRG32K3A = (UInt32)(121)
const CURAND_RNG_PSEUDO_MTGP32 = (UInt32)(141)
const CURAND_RNG_PSEUDO_MT19937 = (UInt32)(142)
const CURAND_RNG_PSEUDO_PHILOX4_32_10 = (UInt32)(161)
const CURAND_RNG_QUASI_DEFAULT = (UInt32)(200)
const CURAND_RNG_QUASI_SOBOL32 = (UInt32)(201)
const CURAND_RNG_QUASI_SCRAMBLED_SOBOL32 = (UInt32)(202)
const CURAND_RNG_QUASI_SOBOL64 = (UInt32)(203)
const CURAND_RNG_QUASI_SCRAMBLED_SOBOL64 = (UInt32)(204)
# end enum curandRngType

typealias curandRngType_t curandRngType

# begin enum curandOrdering
typealias curandOrdering UInt32
const CURAND_ORDERING_PSEUDO_BEST = (UInt32)(100)
const CURAND_ORDERING_PSEUDO_DEFAULT = (UInt32)(101)
const CURAND_ORDERING_PSEUDO_SEEDED = (UInt32)(102)
const CURAND_ORDERING_QUASI_DEFAULT = (UInt32)(201)
# end enum curandOrdering

typealias curandOrdering_t curandOrdering

# begin enum curandDirectionVectorSet
typealias curandDirectionVectorSet UInt32
const CURAND_DIRECTION_VECTORS_32_JOEKUO6 = (UInt32)(101)
const CURAND_SCRAMBLED_DIRECTION_VECTORS_32_JOEKUO6 = (UInt32)(102)
const CURAND_DIRECTION_VECTORS_64_JOEKUO6 = (UInt32)(103)
const CURAND_SCRAMBLED_DIRECTION_VECTORS_64_JOEKUO6 = (UInt32)(104)
# end enum curandDirectionVectorSet

typealias curandDirectionVectorSet_t curandDirectionVectorSet

immutable Array_32_UInt32
    d1::UInt32
    d2::UInt32
    d3::UInt32
    d4::UInt32
    d5::UInt32
    d6::UInt32
    d7::UInt32
    d8::UInt32
    d9::UInt32
    d10::UInt32
    d11::UInt32
    d12::UInt32
    d13::UInt32
    d14::UInt32
    d15::UInt32
    d16::UInt32
    d17::UInt32
    d18::UInt32
    d19::UInt32
    d20::UInt32
    d21::UInt32
    d22::UInt32
    d23::UInt32
    d24::UInt32
    d25::UInt32
    d26::UInt32
    d27::UInt32
    d28::UInt32
    d29::UInt32
    d30::UInt32
    d31::UInt32
    d32::UInt32
end

zero(::Type{Array_32_UInt32}) = begin  # /home/shindo/local-lemon/.julia/v0.4/Clang/src/wrap_c.jl, line 266:
        Array_32_UInt32(fill(zero(UInt32),32)...)
    end

typealias curandDirectionVectors32_t Array_32_UInt32

immutable Array_64_Culonglong
    d1::Culonglong
    d2::Culonglong
    d3::Culonglong
    d4::Culonglong
    d5::Culonglong
    d6::Culonglong
    d7::Culonglong
    d8::Culonglong
    d9::Culonglong
    d10::Culonglong
    d11::Culonglong
    d12::Culonglong
    d13::Culonglong
    d14::Culonglong
    d15::Culonglong
    d16::Culonglong
    d17::Culonglong
    d18::Culonglong
    d19::Culonglong
    d20::Culonglong
    d21::Culonglong
    d22::Culonglong
    d23::Culonglong
    d24::Culonglong
    d25::Culonglong
    d26::Culonglong
    d27::Culonglong
    d28::Culonglong
    d29::Culonglong
    d30::Culonglong
    d31::Culonglong
    d32::Culonglong
    d33::Culonglong
    d34::Culonglong
    d35::Culonglong
    d36::Culonglong
    d37::Culonglong
    d38::Culonglong
    d39::Culonglong
    d40::Culonglong
    d41::Culonglong
    d42::Culonglong
    d43::Culonglong
    d44::Culonglong
    d45::Culonglong
    d46::Culonglong
    d47::Culonglong
    d48::Culonglong
    d49::Culonglong
    d50::Culonglong
    d51::Culonglong
    d52::Culonglong
    d53::Culonglong
    d54::Culonglong
    d55::Culonglong
    d56::Culonglong
    d57::Culonglong
    d58::Culonglong
    d59::Culonglong
    d60::Culonglong
    d61::Culonglong
    d62::Culonglong
    d63::Culonglong
    d64::Culonglong
end

zero(::Type{Array_64_Culonglong}) = begin  # /home/shindo/local-lemon/.julia/v0.4/Clang/src/wrap_c.jl, line 266:
        Array_64_Culonglong(fill(zero(Culonglong),64)...)
    end

typealias curandDirectionVectors64_t Array_64_Culonglong
typealias curandGenerator_st Void
typealias curandGenerator_t Ptr{curandGenerator_st}
typealias curandDistribution_st Cdouble
typealias curandDistribution_t Ptr{curandDistribution_st}
typealias curandDistributionShift_st Void
typealias curandDistributionShift_t Ptr{curandDistributionShift_st}
typealias curandDistributionM2Shift_st Void
typealias curandDistributionM2Shift_t Ptr{curandDistributionM2Shift_st}
typealias curandHistogramM2_st Void
typealias curandHistogramM2_t Ptr{curandHistogramM2_st}
typealias curandHistogramM2K_st UInt32
typealias curandHistogramM2K_t Ptr{curandHistogramM2K_st}
typealias curandHistogramM2V_st curandDistribution_st
typealias curandHistogramM2V_t Ptr{curandHistogramM2V_st}
typealias curandDiscreteDistribution_st Void
typealias curandDiscreteDistribution_t Ptr{curandDiscreteDistribution_st}

# begin enum curandMethod
typealias curandMethod UInt32
const CURAND_CHOOSE_BEST = (UInt32)(0)
const CURAND_ITR = (UInt32)(1)
const CURAND_KNUTH = (UInt32)(2)
const CURAND_HITR = (UInt32)(3)
const CURAND_M1 = (UInt32)(4)
const CURAND_M2 = (UInt32)(5)
const CURAND_BINARY_SEARCH = (UInt32)(6)
const CURAND_DISCRETE_GAUSS = (UInt32)(7)
const CURAND_REJECTION = (UInt32)(8)
const CURAND_DEVICE_API = (UInt32)(9)
const CURAND_FAST_REJECTION = (UInt32)(10)
const CURAND_3RD = (UInt32)(11)
const CURAND_DEFINITION = (UInt32)(12)
const CURAND_POISSON = (UInt32)(13)
# end enum curandMethod

typealias curandMethod_t curandMethod


text:
using Plots, Colors, PyCall, Compose
pyplot()
include("graph_environment.jl")
# This code simulates a patrol of a grid where expected commute time
# averaged over all pairs of nodes is minimized.

"""
    find_path(time,startNode,mg)

Takes in amount of time to simulate, t, the initial node, node0, and the
MetaGraph, mg, and outputs the calculated travel path
"""
function find_path(time_end,startNode,mg)

    M = mg
    n = nv(M)
    k = sqrt(n)
    # path in node form
    Path = zeros(time_end,1);
    Path[1] = startNode;
    for t = 2:time_end
        x = rand();
        num_neighbors = length(neighbors(M,Int(Path[t-1])))
        prob_choices = zeros(num_neighbors,1)
        for neighbor in enumerate(neighbors(M,Int(Path[t-1])))
            prob_choices[neighbor[1]] = get_prop(M,Int(Path[t-1]),neighbor[2],:weight)
            if x < sum(prob_choices[1:neighbor[1]])
                Path[t] = neighbor[2]
                break
            end
        end
    end

    # path in x and y coordinates
    travel = zeros(time_end,2);
    for i = 1:length(Path)
        if mod(Path[i],k) != 0
            travel[i,1] = k*(Path[i]/k - floor(Path[i]/k));
            travel[i,2] = ceil(Path[i]/k)
        else
            travel[i,1]  = k;
            travel[i,2] = ceil(Path[i]/k)
        end
    end
    return travel
end
"""
    plot_path(travel,fn,mg)

Plots the travel path given (also needs the overall MetaGraph) and saves it to
the given file location
"""
function plot_path(travel,fn,mg)
    pyplot()
    M = mg
    n = nv(M)
    k = sqrt(n)
    time_end = size(travel,1)
    finalplot = plot(legend=false)
    title!("Patrolling of Robots")
    xlabel!("x position")
    ylabel!("y position")
    plot!(collect(1:(k-1)/99:k), ones(100,1), color = :black)
    plot!(collect(1:(k-1)/99:k), ones(100,1)*k, color = :black)
    plot!(ones(100,1), collect(1:(k-1)/99:k), color = :black)
    plot!(ones(100,1)*k, collect(1:(k-1)/99:k), color = :black)
    for ii = 1:time_end-1
        line_color = weighted_color_mean(ii/time_end, colorant"black", colorant"red")
        plot!(travel[ii:ii+1,1], travel[ii:ii+1,2],color = line_color)
    end
    # annotate!(k+5,0,text("this is #5", 16, :red, :center))
    png(finalplot,fn)
    return finalplot
end

"""
    main

I don't believe in the concept of "main" functions, but whatever.
"""
function main()
    time_end = 1000
    startNode = 32
    A = []
    path = find_path(time_end,startNode,A)
    plot_path(path,"ethantest.png",A)
end


text:
export MDString, MDNode, operands, Metadata

@checked struct MetadataAsValue <: Value
    ref::reftype(Value)
end
identify(::Type{Value}, ::Val{API.LLVMMetadataAsValueValueKind}) = MetadataAsValue

# NOTE: the C API doesn't allow us to differentiate between MD kinds,
#       all are wrapped by the opaque MetadataAsValue...

const MDString = MetadataAsValue

MDString(val::String) = MDString(API.LLVMMDString(val, length(val)))

MDString(val::String, ctx::Context) =
    MDString(API.LLVMMDStringInContext(ref(ctx), val, length(val)))

function Base.convert(::Type{String}, md::MDString)
    len = Ref{Cuint}()
    ptr = API.LLVMGetMDString(ref(md), len)
    ptr == C_NULL && throw(ArgumentError("invalid metadata, not a MDString?"))
    return unsafe_string(convert(Ptr{Int8}, ptr), len[])
end

const MDNode = MetadataAsValue

MDNode(vals::Vector{T}) where {T<:Value} =
    MDNode(API.LLVMMDNode(ref.(vals), length(vals)))

MDNode(vals::Vector{T}, ctx::Context) where {T<:Value} =
    MDNode(API.LLVMMDNodeInContext(ref(ctx), ref.(vals), length(vals)))

function operands(md::MDNode)
    nops = API.LLVMGetMDNodeNumOperands(ref(md))
    ops = Vector{API.LLVMValueRef}(undef, nops)
    API.LLVMGetMDNodeOperands(ref(md), ops)
    return Value[Value(op) for op in ops]
end

@checked struct Metadata
    ref::API.LLVMMetadataRef
end
reftype(::Type{Metadata}) = API.LLVMMetadataRef

function Metadata(val::Value)
    return Metadata(LLVM.API.LLVMValueAsMetadata(ref(val)))
end

function Value(md::Metadata, ctx::Context)
    return MetadataAsValue(API.LLVMMetadataAsValue(ref(ctx), ref(md)))
end


text:
"""
    function save_plot(;
        y_var::Vector{Float64},
        title::String,
        f_name::String,
        xlabel::String,
        ylabel::String,
        tournaments::Vector{Tournament}
    )

Creates a plot.

Arguments
=========
- `y_var`:
    y variable.
- `title`:
    plot title.
- `f_name`:
    name of plot file.
- `xlabel`:
    x variable label.
- `ylabel`:
    y variable label.
- `tournaments`:
    tournaments played.
"""

function save_plot(;
    y_var::Vector{Float64},
    title::String,
    f_name::String,
    xlabel::String,
    ylabel::String,
    tournaments::Vector{Tournament}
)
    tournament_types = permutedims(map(t -> t.type, tournaments))
    plot(
        y_var,
        title=title,
        seriestype=:scatter,
        xlabel=xlabel,
        ylabel=ylabel,
        labels=tournament_types,
        titlefont=font(10,"Computer Modern"),
        yguidefont=font(8,"Computer Modern"),
        xguidefont=font(8,"Computer Modern"),
    )
    savefig("result/$f_name.png")
end

max_stars_repo_name: Danial-Hussain/Tennis-Ranking

max_issues_repo_name: Danial-Hussain/Tennis-Ranking

max_forks_repo_name: Danial-Hussain/Tennis-Ranking


text:
using ICU
using Base.Test

# Tests for not overrunning buffer
str = "\u3b0"
upp = "\u3a5\u308\u301"

@test ICU.toupper(ICU.cvt_utf16(str^8)) == ICU.cvt_utf16(upp^8)

@testset "utext" begin
end

@testset "ustring (general string)" begin
end

@testset "ucasemap (Case mapping)" begin
    @test ICU.tolower("ABCDEF") == "abcdef"
    @test ICU.toupper("abcdef") == "ABCDEF"
    @test ICU.totitle("this") == "This"
end

@testset "ubrk (Word/Line/Sentence breaking)" begin
    str = ICU.cvt_utf16("This is a test of splitting sentences and words.  Scott is busy.\nHow does this\nsplit? Tell me.")
    for (sep, res) in ((UBRK.WORD,
                        [4,5,7,8,9,10,14,15,17,18,27,28,37,38,41,42,47,48,49,50,
                         55,56,58,59,63,64,65,68,69,73,74,78,79,84,85,86,90,91,93,94]),
                       (UBRK.SENTENCE, [50,65,79,86,94]),
                       (UBRK.LINE, [5,8,10,15,18,28,38,42,50,56,59,65,69,74,79,86,91,94]))
        brk = UBrk(sep, str)
        vec = Vector{Int}()
        while (v = ICU.next(brk)) != UBRK.DONE
            push!(vec, v)
        end
        close(brk)
        @test vec == res
    end
end

@testset "custom ubrk (Word/Line/Sentence breaking)" begin
    str = ICU.cvt_utf16("This is a test of splitting sentences and words with custom rules")
    rul = ICU.cvt_utf16(".+ {0};")
    brk = UBrk(rul, str)
    vec = Vector{Int}()
    while (v = ICU.next(brk)) != UBRK.DONE
            push!(vec, v)
    end
    close(brk)
    @test length(vec) == 1
end

@testset "ucnv (Conversion)" begin
    @test ICU.tolower(ICU.cvt_utf16("ABCDEF")) == "abcdef"
    @test ICU.toupper(ICU.cvt_utf16("abcdef")) == "ABCDEF"
    @test ICU.totitle(ICU.cvt_utf16("this")) == "This"
end

@testset "ucol (Collation)" begin
end

@testset "ucal (Calendar)" begin
    fields = [UCAL.YEAR, UCAL.MONTH, UCAL.DATE, UCAL.HOUR_OF_DAY, UCAL.MINUTE, UCAL.SECOND]
    cal = UCalendar("America/Los_Angeles")
    ICU.set_millis!(cal, 1486735888612.)
    @test ICU.get(cal,fields) == [2017, 1, 10, 6, 11, 28]
    cal = UCalendar("Belgium/Antwerp")
    ICU.set_millis!(cal, 1486735888612.)
    @test ICU.get(cal,fields) == [2017, 1, 10, 14, 11, 28]
    ICU.clear!(cal)
end

@testset "udat (Date)" begin
    df = UDateFormat(UDAT.SHORT, UDAT.SHORT, "en_US")
    @test ICU.format(df, 1486735888612.) == "2017-02-10 14:11"
end

@testset "ucsdet (Character Set Detection)" begin
    csd = UCharsetDetector()
    for (s, det) in ((b"This is ASCII test, let's see how it does", "ISO-8859-1"),
                     (b"See if it detects Windows-1252, with the Euro \x80 and other characters \x8c \x9c", "windows-1252"),
                     ("にほんでは、近頃ちかごろ多おおくの人ひとが保育園ほいくえん問題もんだいについて話はなしている。特とくに東京とうきょうでは十分じゅうぶんな施設しせつがないので、子こどもを保育園ほいくえんに入いれることがとても大変たいへんだ。今いま私わたしは東京とうきょうに住すんでいるので、息子むすこを保育園ほいくえんに入いれるのは不可能ふかのうだろうと思おもっていた。しかし驚おどろいたことに、息子むすこは受うけ入いれてもらえた", "UTF-8"))
        ICU.set!(csd, s)
        m = ICU.detect(csd)
        @test m != nothing
        if m != nothing
            @test ICU.get_name(m) == det
        end
    end
    ICU.close(csd)
end


text:
module SHIPs

const IntS = Int32

import JuLIP.Potentials: z2i

using Reexport
@reexport using JuLIP

include("aux.jl")
include("prototypes.jl")

# specification of the radial basis
include("jacobi.jl")
include("transforms.jl")

# specification of the angular basis
include("sphericalharmonics.jl")
include("rotations.jl")

# basis specification: subsets of the full expansion
include("basisspecs.jl")

# implements the A functions ∏A functions
include("Alist.jl")

# SHIPBasis definition
include("basis.jl")

# SHIP interatomic potential definition
include("fast.jl")

include("regularisers.jl")

include("descriptors.jl")

include("utils.jl")

end # module

max_stars_repo_name: cortner/SHIPs

max_issues_repo_name: cortner/SHIPs

max_forks_repo_name: cortner/SHIPs


text:
module Interpret2
using ..LambdaC: AExpr
export empty_env, Environment
import MLStyle

const Environment = NamedTuple
empty_env() = NamedTuple()

"Produces new environment Γ' s.t. `Γ(x) = v` (and everything else unchanged)"
update(Γ, x::Symbol, v) = merge(Γ, NamedTuple{(x,)}((v,)))

function interpret_assign(aex, Γ)
  x = aex.args[1] # Variable that we're updating
  v = aex.args[2] # Value that we're updating it to
  Γ_ = update(Γ, x, v)
end

function interpret_program(aex, Γ)
  aex.head == :program || error("Must be a program aex")
  for stmt in aex.args
    Γ = interpret_assign(stmt, Γ)
  end
  return Γ
end

prim_to_func = Dict(:+ => +,
                    :- => -,
                    :* => *,
                    :/ => /)

isprim(f) = f in keys(prim_to_func)
primapl(f, x...) = prim_to_func[f](x...)

function interpret(aex, Γ)
  arr = [aex.head, aex.args...]
  next(x) = interpret(x, Γ)
  isaexpr(x) = x isa AExpr
  MLStyle.@match arr begin
    [:program, args...]                            => (AExpr(:assign, x, interpret(v, Γ)), Γ)
    [:assign, x, v::AExpr]                         => (AExpr(:assign, x, interpret(v, Γ)), Γ)
    [:assign, x, v::Symbol]                        => (AExpr(:assign, x, interpret(v, Γ)), Γ)
    [:assign, x, v]                                => (AExpr(:assign, x, v), update(Γ, x, v))
    s::Symbol                                      => (Γ[x], Γ)
    [:if, true, t, e]                              => (t, Γ)
    [:if, false, t, e]                             => (e, Γ)
    [:if, i::AExpr, t, e]                          => (AExpr(:if, interpret(i, Γ), t, e), Γ)
    [:call, f, args...] && if !isprim(f) end       => (AExpr(:call, next(f), args...), Γ)
    [:call, f, args...] && if any(isaexpr, args) end  => (AExpr(:call, map((i, arg) -> i == shit ? next(arg) : arg)), Γ)
    [:call, f, args...] && isprim(f) end           => (primapl(f, args...), Γ)
    _                                              => error("Could not interpret $arr")
  end
end

# Interpret until termination
function interpret_terminate(aex, env)
  while !done(aex, env)
    aex, env = interpret(aex, env)
  end
  aex, env
end

function start(aex::AExpr)::Environment
  env = empty_env()
  aex_, env_ = interpret_terminate(aex, env)
  env_
end

function step(aex::AExpr, env::Environment)::Environment
  aex_, env_ = interpret_terminate(aex, env)
  env_
end

end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using Plots
using Plots.PlotMeasures
using DelimitedFiles
using Random, Distributions

"""
Call as a jl script on any function in this file.
e.g. julia -L plot_errors.jl -e 'plot_errors(args)'.

Author: Ariel Kellison, 12/2021
"""

fpath = @__DIR__

function plot_errors(err_files::Array{String,1}, steps::Array{Int64,1})
	"""
	Args: err_files is an array of the files to be read
	"""
	# read data as matrices from files and save max from each file

	max_errs = []

	for file in err_files
		out_xv   = readdlm(file, Float64)
		# data to vec
		vec_xv       = [x::Float64 for x in out_xv]
		len          = convert(Int64,length(vec_xv)/2)
		out_dat_x    = [x for x in out_xv[1:len]]
		out_dat_v    = [v for v in out_xv[len+1:len*2]]

		norms = []

		for i = 1:length(out_dat_x)
			x    = out_dat_x[i]
			v    = out_dat_v[i]
			append!(norms, sqrt(x^2 + v^2))
		end
		append!(max_errs, maximum(norms))
	end

	norm = ((1.235*(1/10^7))^2+(6.552*(1/10^8))^2)^0.5
	bnd_array = [0.0]
	sb = 1.000003814704543
	for n in steps
		if (n >0)
			append!(bnd_array, norm*(sb^n-1)/(sb-1))
		end
	end
	print(bnd_array)
	print(max_errs)

	scatter([n for n in steps], bnd_array, label="VCFloat norm bound",
	left_margin = 15mm, right_margin = 15mm)

	plotx = scatter!(
	[n for n in steps], max_errs, ylabel = "floating-point error",xlabel = "number of steps",
	size = (900, 500),bottom_margin = 10mm,left_margin = 10mm,
	yguidefontsize=12, xguidefontsize=12, xtickfontsize = 12, ytickfontsize = 12,
	label="empirical error",
	title = "Leapfrog integration normwise absolute error vs iteration number"
	)

	savefig(plotx,fpath*"/normwise_errors.png")

end


text:
#------------------------------
# functions for minibatch approach
#------------------------------

function batchloss_individual_parameters(bX, bY, bts, bweights, m::odevae)
    batchloss = 0.0
    for ind in 1:batchsize
        latentμ, latentlogσ = m.encodedμ(m.encoder(bX[ind])), m.encodedlogσ(m.encoder(bX[ind]))
        curts = Int.(bts[ind] .*(1 ./dt) .+1)
        curparams = m.paramNN(bY[ind])
        smoothμ = Array(solve(m.ODEprob, Tsit5(), u0 = [latentμ[1,1], latentμ[2,1]], p=curparams, saveat=dt))[:,curts]
        combinedμ = hcat(latentμ[:,1],smoothμ)
        combinedz = latentz.(combinedμ, latentlogσ)
        ELBO = 1.0 .* logp_x_z(m, bX[ind], combinedz) .- 0.5 .* kl_q_p(combinedμ, latentlogσ)
        lossval = sum(-ELBO) .+ 0.01*reg(m) .+ 0.5 .* sum((smoothμ .- latentμ[:,2:end]).^2)
        batchloss += lossval * bweights[ind]
    end
    return batchloss
end

function batchloss_wrapper(m::odevae)
    return function(bX, bY, bts, bweights) batchloss_individual_parameters(bX, bY, bts, bweights, m) end
end

function tricube(x)
    if abs(x) <= 1
        val = (1 .- abs.(x).^3).^3
    else
        val = 0
    end
    return val
end

function maxk(a, indrange) #inspired and adapted from https://discourse.julialang.org/t/what-is-julias-maxk-matlab-that-returns-the-indice-of-top-k-largest-values/14100/7
    maxkinds = partialsortperm(a, indrange)
    return [maxkinds, a[maxkinds]]
end

euclidean_dist(x,y) = length(x) == length(y) ? sqrt(sum((x .- y).^2)) : error("lengths of the two vectors don't match")

mean_dist(x,y) = length(x) == length(y) ? (1/length(x))*sum(abs.(x.-y)) : error("lengths of the two vectors don't match")

function randomminibatches(Y, batchsize)
    minibatches = []
    n=length(Y)
    for ind in 1:n
        distrange = collect(1:n)[1:end .!= ind]
        randombatch = cat(ind, shuffle(distrange)[1:(batchsize-1)], dims=1)
        push!(minibatches, randombatch)
    end
    randomweights = collect(fill(1/batchsize, batchsize) for i=1:n)
    minibatches, randomweights
end

function getdistmat_odesols_mean(m::odevae, xs, x_params; centralise=true)
    n = length(xs)
    odeparams = m.paramNN.(x_params)
    odeμ0s = collect(m.encodedμ.(m.encoder.(xs))[ind][:,1] for ind in 1:n)
    tgrid = collect(11:10:101)
    dims = length(odeμ0s[1])
    odesols = collect(Array(solve(m.ODEprob, Tsit5(), u0 = odeμ0s[ind], p=odeparams[ind], saveat=dt))[:,tgrid] for ind in 1:n)
    if centralise # additional centralising step: subtract mean of trajectory for all inds and dims
        for ind in 1:n
            for d in 1:dims
                mean_per_dim = mean(odesols[ind][d,j] for j in 1:length(tgrid))
                odesols[ind][d,:] .-= mean_per_dim
            end
        end
    end
    mean_distmat = fill(0.0,(n,n));     # für jeden Zeitpunkt eine Distanzmatrix, dann mitteln
    for i in 1:n
        for j in 1:i
            mean_distmat[i,j] = (1.0/sqrt(length(tgrid))) .* sum(euclidean_dist(odesols[i][:,tp], odesols[j][:,tp]) for tp in 1:length(tgrid))
            #mean_distmat[i,j] = (1.0/sqrt(length(tgrid))) .* mean(sqrt(sum((odesols[i][d,tp]-odesols[j][d,tp]).^2 for tp in 1:length(tgrid))) for d in 1:dims)
            mean_distmat[j,i] = mean_distmat[i,j]
        end
    end
    return mean_distmat
end

function findminibatches_distmat(distmat, batchsize, kernel; bandwidth=1.0)
    n = size(distmat,2)
    minibatches =  []
    allweights = []
    for ind in 1:n
        distrange = collect(1:n)[1:end .!= ind]
        distvec = collect(distmat[ind,j] for j in distrange)
        batch_dist = maxk(distvec,1:(batchsize-1))
        map!(x -> (x >= ind ? x+1 : x), batch_dist[1], batch_dist[1])
        new_batch_dist = cat(0, batch_dist[2]..., dims=1)
        minibatch = cat(ind, Int.(batch_dist[1]), dims=1)
        kernelsum = sum(kernel.(new_batch_dist ./ bandwidth))
        weights = kernel.(new_batch_dist ./ bandwidth) ./kernelsum
        push!(allweights, weights)
        push!(minibatches, minibatch)
    end
    minibatches, allweights
end

function findminibatches_distmat(distmat, batchsize)
    n = size(distmat,2)
    minibatches =  []
    for ind in 1:n
        distrange = collect(1:m)[1:end .!= ind]
        distvec = collect(distmat[ind,j] for j in distrange)
        batch_dist = maxk(distvec,1:(batchsize-1))
        map!(x -> (x >= ind ? x+1 : x), batch_dist[1], batch_dist[1])
        new_batch_dist = cat(0, batch_dist[2]..., dims=1)
        minibatch = cat(ind, Int.(batch_dist[1]), dims=1)
        push!(minibatches, minibatch)
    end
    allweights = fill(fill(1/batchsize, batchsize), n)
    minibatches, allweights
end

# count the number of individuals of the same group as the reference individual in each minibatch
function evaluate_minibatches(minibatches, group1, group2)
    correctly_classified = 0
    for ind in 1:length(minibatches)
        cur_minibatch = minibatches[ind]
        correct_group = cur_minibatch[1] ∈ group1 ? group1 : group2
        curnum = length(intersect(cur_minibatch[2:end],correct_group))
        correctly_classified += curnum
    end
    return correctly_classified
end

function evaluate_minibatches_pergroup(minibatches, group1, group2)
    correctly_classified_group1 = 0
    correctly_classified_group2 = 0
    for ind in 1:length(minibatches)
        cur_minibatch = minibatches[ind]
        correct_group = cur_minibatch[1]∈group1 ? group1 : group2
        curnum = length(intersect(cur_minibatch[2:end],correct_group))
        if correct_group == group1
            correctly_classified_group1 += curnum
        else
            correctly_classified_group2 += curnum
        end
    end
    total_gr1 = length(group1) .* (batchsize-1)
    total_gr2 = length(group2) .* (batchsize-1)
    total = length(minibatches) .* (batchsize-1)
    return (correctly_classified_group1 ./total_gr1 .*100),
    (correctly_classified_group2 ./total_gr2 .* 100),
    ((correctly_classified_group1 + correctly_classified_group2)./total .*100)
end


text:
@testset "Memory" begin

#= FIXME
let
    a,b = Mem.info()
    # NOTE: actually testing this is pretty fragile on CI
    #=@test a == =# Mem.free()
    #=@test b == =# Mem.total()
    #=@test b-a == =# Mem.used()
end
=#

# pointer-based
src = 42

buf1 = Mem.alloc(sizeof(src))

Mem.set!(buf1, UInt32(0), sizeof(Int)÷sizeof(UInt32))

Mem.upload!(buf1, Ref(src), sizeof(src))

dst1 = Ref(0)
Mem.download!(dst1, buf1, sizeof(src))
@test src == dst1[]

buf2 = Mem.alloc(sizeof(src))

Mem.transfer!(buf2, buf1, sizeof(src))

dst2 = Ref(0)
Mem.download!(dst2, buf2, sizeof(src))
@test src == dst2[]

Mem.free(buf2)
Mem.free(buf1)

# array-based
src = [42]

buf1 = Mem.alloc(src)

Mem.upload!(buf1, src)

dst1 = similar(src)
Mem.download!(dst1, buf1)
@test src == dst1

buf2 = Mem.upload(src)

dst2 = similar(src)
Mem.download!(dst2, buf2)
@test src == dst2

Mem.free(buf1)

# type-based
buf = Mem.alloc(Int)

# there's no type-based upload, duh
src = [42]
Mem.upload!(buf, src)

dst = Mem.download(eltype(src), buf)
@test src == dst

let
    @test_throws ArgumentError Mem.alloc(Function, 1)   # abstract
    @test_throws ArgumentError Mem.alloc(Array{Int}, 1) # UnionAll
    @test_throws ArgumentError Mem.alloc(Integer, 1)    # abstract
    # TODO: can we test for the third case?
    #       !abstract && leaftype seems to imply UnionAll nowadays...

    # zero-width allocations should be permitted
    null = Mem.alloc(Int, 0)
    Mem.free(null)

    # double-free should throw
    x = Mem.alloc(1)
    Mem.free(x)
    # FIXME: Segfaults... @test_throws HSAError Mem.free(x)
end

let
    @eval mutable struct MutablePtrFree
        foo::Int
        bar::Int
    end
    buf = Mem.alloc(MutablePtrFree)
    Mem.upload!(buf, [MutablePtrFree(0,0)])
    Mem.free(buf)
end

let
    @eval mutable struct MutableNonPtrFree
        foo::Int
        bar::String
    end
    @test_throws ArgumentError Mem.alloc(MutableNonPtrFree)
end

end


text:
include(joinpath(dirname(@__DIR__), "src", "FaceDetection.jl"))
using Documenter, .FaceDetection

Documenter.makedocs(
    clean = true,
    doctest = true,
    modules = Module[FaceDetection],
    repo = "",
    highlightsig = true,
    sitename = "FaceDetection Documentation",
    expandfirst = [],
    pages = [
        "Home" => "index.md",
        "Usage" => "usage.md",
        "Examples" => "examples.md",
        "Benchmarking Results" => "benchmarking.md",
        "Caveats" => "caveats.md",
        "Other Resources" => "resources.md",
        "A Few Acknowledgements" => "acknowledgements.md"
    ]
)

deploydocs(;
    repo  =  "github.com/jakewilliami/FaceDetection.jl.git",
)


text:
#Demo 17: Rise/fall tests
#-------------------------------------------------------------------------------

using CMDimCircuits
CMDimCircuits.@using_CData()

#Get a demo display:
include(CMDimCircuits.demoplotcfgscript)
#Normally use something like:
#CMDimData.@includepkg EasyPlotInspect

#==Constants
===============================================================================#
LBLAX_AMPV = "Amplitude (V)"
LBLAX_RFTIMEPS = "Rise/Fall (ps)"
LBLAX_TIME = "Time (s)"
delayattr = cons(:a,
	line = set(style=:solid, width=3),
	glyph = set(shape=:x, size=2),
)

#==Input data
===============================================================================#
tbit = 1e-9 #Bit period
osr = 20 #samples per bit
nsamples = 127
amp = 2

#==Computations
===============================================================================#
seq = 1.0*prbs(reglen=7, seed=1, nsamples=nsamples)
t = DataF1(0:(tbit/osr):(nsamples*tbit))
tmax = maximum(t)

#Generate parameter sweeps:
sweeplist = PSweep[
	PSweep("tau", tbit.*[1/2.5])
]

#Generate data:
Π = fill(DataHR, sweeplist) do tau
	return pulse(t, Pole(1/tau,:rad), tpw=tbit)
end

pat = (pattern(seq, Π, tbit=tbit)-0.5)*amp #Center data pattern

lthresh = -1+amp*.2 #20%
hthresh = -1+amp*.8 #80%
@show lthresh, hthresh
trise = measrise(pat, lthresh=lthresh, hthresh=hthresh)
tfall = measfall(pat, lthresh=lthresh, hthresh=hthresh)

#==Generate plot
===============================================================================#
plot = cons(:plot, nstrips=2, legend=false,
	ystrip1 = set(axislabel=LBLAX_AMPV, striplabel="Patterns"),
	ystrip2 = set(axislabel=LBLAX_RFTIMEPS, striplabel="20-80 Rise/Fall Times"),
	xaxis = set(label=LBLAX_TIME, max=tmax),
)
push!(plot,
	cons(:wfrm, pat, label="pat", strip=1),

	cons(:wfrm, trise/1e-12, delayattr, line=set(color=:blue), label="rise", strip=2),
	cons(:wfrm, tfall/1e-12, delayattr, line=set(color=:red), label="fall", strip=2),
)

pcoll = push!(cons(:plotcoll, title="Rise/Fall Tests"), plot)
	pcoll.ncolumns = 1

#==Display results in pcoll
===============================================================================#
EasyPlot.displaygui(pcoll)

#==Return pcoll to user (call evalfile(...))
===============================================================================#
pcoll