text:
using DynamicalSystemsBase: DEFAULT_DIFFEQ_KWARGS, _get_solver
using Roots: find_zero, A42
export poincaresos, produce_orbitdiagram, PlaneCrossing, poincaremap, PoincareMap

const ROOTS_ALG = A42()

#####################################################################################
#                               Hyperplane                                          #
#####################################################################################
"""
    PlaneCrossing(plane, dir) → z
Create a struct that can be called as a function `z(u)` that returns the signed distance
of state `u` from the hyperplane `plane` (positive means in front of the hyperplane).
See [`poincaresos`](@ref) for what `plane` can be (tuple or vector).
"""
struct PlaneCrossing{P, D, T}
    plane::P
    dir::Bool
    n::SVector{D, T}  # normal vector
    p₀::SVector{D, T} # arbitrary point on plane
end
PlaneCrossing(plane::Tuple, dir) = PlaneCrossing(plane, dir, SVector(true), SVector(true))
function PlaneCrossing(plane::AbstractVector, dir)
    n = plane[1:end-1] # normal vector to hyperplane
    i = findfirst(!iszero, plane)
    D = length(plane)-1; T = eltype(plane)
    p₀ = zeros(D)
    p₀[i] = plane[end]/plane[i] # p₀ is an arbitrary point on the plane.
    PlaneCrossing(plane, dir, SVector{D, T}(n), SVector{D, T}(p₀))
end

# Definition of functional behavior
function (hp::PlaneCrossing{P})(u::AbstractVector) where {P<:Tuple}
    @inbounds x = u[hp.plane[1]] - hp.plane[2]
    hp.dir ? x : -x
end
function (hp::PlaneCrossing{P})(u::AbstractVector) where {P<:AbstractVector}
    x = zero(eltype(u))
    D = length(u)
    @inbounds for i in 1:D
        x += u[i]*hp.plane[i]
    end
    @inbounds x -= hp.plane[D+1]
    hp.dir ? x : -x
end

#####################################################################################
#                               Poincare Section                                    #
#####################################################################################
"""
    poincaresos(ds::ContinuousDynamicalSystem, plane, tfinal = 1000.0; kwargs...)
Calculate the Poincaré surface of section[^Tabor1989]
of the given system with the given `plane`.
The system is evolved for total time of `tfinal`.
Return a [`Dataset`](@ref) of the points that are on the surface of section.
See also [`poincaremap`](@ref) for the map version.

If the state of the system is ``\\mathbf{u} = (u_1, \\ldots, u_D)`` then the
equation defining a hyperplane is
```math
a_1u_1 + \\dots + a_Du_D = \\mathbf{a}\\cdot\\mathbf{u}=b
```
where ``\\mathbf{a}, b`` are the parameters of the hyperplane.

In code, `plane` can be either:

* A `Tuple{Int, <: Number}`, like `(j, r)` : the plane is defined
  as when the `j` variable of the system equals the value `r`.
* A vector of length `D+1`. The first `D` elements of the
  vector correspond to ``\\mathbf{a}`` while the last element is ``b``.

This function uses `ds` and higher order interpolation from DifferentialEquations.jl
to create a high accuracy estimate of the section.
See also [`produce_orbitdiagram`](@ref).

## Keyword Arguments
* `direction = -1` : Only crossings with `sign(direction)` are considered to belong to
  the surface of section. Positive direction means going from less than ``b``
  to greater than ``b``.
* `idxs = 1:dimension(ds)` : Optionally you can choose which variables to save.
  Defaults to the entire state.
* `Ttr = 0.0` : Transient time to evolve the system before starting
  to compute the PSOS.
* `u0 = get_state(ds)` : Specify an initial state.
* `warning = true` : Throw a warning if the Poincaré section was empty.
* `rootkw = (xrtol = 1e-6, atol = 1e-6)` : A `NamedTuple` of keyword arguments
  passed to `find_zero` from [Roots.jl](https://github.com/JuliaMath/Roots.jl).
* `diffeq` is a `NamedTuple` (or `Dict`) of keyword arguments propagated into
  `init` of DifferentialEquations.jl.
  See [`trajectory`](@ref) for examples. Only valid for continuous systems.

[^Tabor1989]: M. Tabor, *Chaos and Integrability in Nonlinear Dynamics: An Introduction*, §4.1, in pp. 118-126, New York: Wiley (1989)
"""
function poincaresos(
		ds::CDS{IIP, S, D}, plane, tfinal = 1000.0;
	    direction = -1, Ttr::Real = 0.0, warning = true, idxs = 1:D, u0 = get_state(ds),
	    rootkw = (xrtol = 1e-6, atol = 1e-6), diffeq = NamedTuple(), kwargs...
	) where {IIP, S, D}

    if !isempty(kwargs)
        @warn DIFFEQ_DEP_WARN
        diffeq = NamedTuple(kwargs)
    end

    _check_plane(plane, D)
    integ = integrator(ds, u0; diffeq)
    i = typeof(idxs) <: Int ? idxs : SVector{length(idxs), Int}(idxs...)
    planecrossing = PlaneCrossing(plane, direction > 0)
	Ttr ≠ 0 && step!(integ, Ttr)
	plane_distance = (t) -> planecrossing(integ(t))

	data = _poincaresos(integ, plane_distance, planecrossing, tfinal+Ttr, i, rootkw)
    warning && length(data) == 0 && @warn PSOS_ERROR
    return Dataset(data)
end

# The separation into two functions here exists only to introduce a function barrier
# for the low level method, to ensure optimization on argments of `poincaremap!`.
function _poincaresos(integ, plane_distance, planecrossing, tfinal, i, rootkw)
	data = _initialize_output(integ.u, i)
	while integ.t < tfinal
		out = poincaremap!(integ, plane_distance, planecrossing, tfinal, rootkw)
		if !isnothing(out)
            push!(data, out[i])
        else
            break # if we evolved for more than tfinal, we should break anyways.
        end
	end
	return data
end

_initialize_output(::S, ::Int) where {S} = eltype(S)[]
_initialize_output(u::S, i::SVector{N, Int}) where {N, S} = typeof(u[i])[]
function _initialize_output(u, i)
    error("The variable index when producing the PSOS must be Int or SVector{Int}")
end

const PSOS_ERROR = "The Poincaré surface of section did not have any points!"

"""
	poincaremap!(integ, plane_distance, planecrossing, Tmax, rootkw)
Low level function that actual performs the algorithm of finding the next crossing
of the Poincaré surface of section. Return the state at the section or `nothing` if
evolved for more than `Tmax` without any crossing.
"""
function poincaremap!(integ, plane_distance, planecrossing, Tmax, rootkw)
    t0 = integ.t
    # Check if initial condition is already on the plane
    side = planecrossing(integ.u)
    if side == 0
		dat = integ.u
        step!(integ)
		return dat
    end
    # Otherwise evolve until juuuuuust crossing the plane
    while side < 0
        (integ.t - t0) > Tmax && break
        step!(integ)
        side = planecrossing(integ.u)
    end
    while side ≥ 0
        (integ.t - t0) > Tmax && break
        step!(integ)
        side = planecrossing(integ.u)
    end
    (integ.t - t0) > Tmax && return nothing # we evolved too long and no crossing, return nothing
    # Else, we're guaranteed to have `t` after plane and `tprev` before plane
    tcross = Roots.find_zero(plane_distance, (integ.tprev, integ.t), Roots.A42(); rootkw...)
    ucross = integ(tcross)
    return ucross
end

function _check_plane(plane, D)
    P = typeof(plane)
    L = length(plane)
    if P <: AbstractVector
        if L != D + 1
            throw(ArgumentError(
            "The plane for the `poincaresos` must be either a 2-Tuple or a vector of "*
            "length D+1 with D the dimension of the system."
            ))
        end
    elseif P <: Tuple
        if !(P <: Tuple{Int, Number})
            throw(ArgumentError(
            "If the plane for the `poincaresos` is a 2-Tuple then "*
            "it must be subtype of `Tuple{Int, Number}`."
            ))
        end
    else
        throw(ArgumentError(
        "Unrecognized type for the `plane` argument."
        ))
    end
end

#####################################################################################
#                               Poincare Map                                        #
#####################################################################################
"""
	poincaremap(ds::ContinuousDynamicalSystem, plane, Tmax=1e6; kwargs...) → pmap

Return a map (integrator) that produces iterations over the Poincaré map of the `ds`.
This map is defined as the sequence of points on the Poincaré surface of section.
See [`poincaresos`](@ref) for details on `plane` and all other `kwargs`.
Keyword `idxs` does not apply to `poincaremap`, as it doesn't save any states.

You can progress the map one step on the section by calling `step!(pmap)`,
which also returns the next state crossing the hyperplane.
You can also set the integrator to start from a new
state `u`, which doesn't have to be on the hyperplane, by using `reinit!(pmap, u)`
and then calling `step!` as normally.

**Notice**: The argument `Tmax` exists so that the integrator can terminate instead
of being evolved for infinite time, to avoid cases where iteration would continue
forever for ill-defined hyperplanes or for convergence to fixed points.
If during `step!` the system has been evolved for more than `Tmax`,
then `step!(pmap)` will return `nothing`, while `get_state(pmap)` will return a vector
of `NaN`.

## Example
```julia
ds = Systems.rikitake([0.,0.,0.], μ = 0.47, α = 1.0)
pmap = poincaremap(ds, (3,0.), Tmax=200.)
next_state_on_psos = step!(pmap)
# Change initial condition
reinit!(pmap, [1.0, 0, 0])
next_state_on_psos = step!(pmap)
```
"""
function poincaremap(
		ds::CDS{IIP, S, D}, plane, Tmax = 1e6;
	    direction = -1, u0 = get_state(ds),
	    rootkw = (xrtol = 1e-6, atol = 1e-6), diffeq = NamedTuple(), kwargs...
	) where {IIP, S, D}

    if !isempty(kwargs)
        @warn DIFFEQ_DEP_WARN
        diffeq = NamedTuple(kwargs)
    end

	_check_plane(plane, D)
    integ = integrator(ds, u0; diffeq)
	planecrossing = PlaneCrossing(plane, direction > 0)
	plane_distance = (t) -> planecrossing(integ(t))
    v = SVector{length(u0), eltype(u0)}(u0)
	return PoincareMap(integ, plane_distance, planecrossing, Tmax, rootkw, v)
end

mutable struct PoincareMap{I, F, P, R, V}
	integ::I
	f::F
 	planecrossing::P
	Tmax::Float64
	rootkw::R
	proj_state::V
end

function DynamicalSystemsBase.step!(pmap::PoincareMap)
	u = poincaremap!(pmap.integ, pmap.f, pmap.planecrossing, pmap.Tmax, pmap.rootkw)
	if isnothing(u)
		return nothing
	else
		pmap.proj_state = u
		return pmap.proj_state
	end
end
function DynamicalSystemsBase.reinit!(pmap::PoincareMap, u0)
	reinit!(pmap.integ, u0)
	return
end
function DynamicalSystemsBase.get_state(pmap::PoincareMap)
	if pmap.integ.t < pmap.Tmax
		return pmap.proj_state
	else
		return pmap.proj_state*NaN
	end
end

function Base.show(io::IO, pmap::PoincareMap)
    println(io, "Iterator of the Poincaré map")
    println(io,  rpad(" rule f: ", 14),     DynamicalSystemsBase.eomstring(pmap.integ.f.f))
    println(io,  rpad(" hyperplane: ", 14),     pmap.planecrossing.plane)
end

#####################################################################################
# Poincare Section for Datasets (trajectories)                                      #
#####################################################################################
# TODO: Nice improvement would be to use cubic interpolation instead of linear,
# using points i-2, i-1, i, i+1
"""
    poincaresos(A::Dataset, plane; kwargs...)
Calculate the Poincaré surface of section of the given dataset with the given `plane`
by performing linear interpolation betweeen points that sandwich the hyperplane.

Argument `plane` and keywords `direction, warning, idxs` are the same as above.
"""
function poincaresos(A::Dataset, plane; direction = -1, warning = true, idxs = 1:size(A, 2))
    _check_plane(plane, size(A, 2))
    i = typeof(idxs) <: Int ? idxs : SVector{length(idxs), Int}(idxs...)
    planecrossing = PlaneCrossing(plane, direction > 0)
    data = poincaresos(A, planecrossing, i)
    warning && length(data) == 0 && @warn PSOS_ERROR
    return Dataset(data)
end
function poincaresos(A::Dataset, planecrossing::PlaneCrossing, j)
    i, L = 1, length(A)
    data = _initialize_output(A[1], j)
    # Check if initial condition is already on the plane
    planecrossing(A[i]) == 0 && push!(data, A[i][j])
    i += 1
    side = planecrossing(A[i])

    while i ≤ L # We always check point i vs point i-1
        while side < 0 # bring trajectory infront of hyperplane
            i == L && break
            i += 1
            side = planecrossing(A[i])
        end
        while side ≥ 0 # iterate until behind the hyperplane
            i == L && break
            i += 1
            side = planecrossing(A[i])
        end
        i == L && break
        # It is now guaranteed that A crosses hyperplane between i-1 and i
        ucross = interpolate_crossing(A[i-1], A[i], planecrossing)
        push!(data, ucross[j])
    end
    return data
end

function interpolate_crossing(A, B, pc::PlaneCrossing{<:AbstractVector})
    # https://en.wikipedia.org/wiki/Line%E2%80%93plane_intersection
    t = LinearAlgebra.dot(pc.n, (pc.p₀ .- A))/LinearAlgebra.dot((B .- A), pc.n)
    return A .+ (B .- A) .* t
end

function interpolate_crossing(A, B, pc::PlaneCrossing{<:Tuple})
    # https://en.wikipedia.org/wiki/Linear_interpolation
    y₀ = A[pc.plane[1]]; y₁ = B[pc.plane[1]]; y = pc.plane[2]
    t = (y - y₀) / (y₁ - y₀) # linear interpolation with t₀ = 0, t₁ = 1
    return A .+ (B .- A) .* t
end


text:
using DiffEqBase

function test_num_args()
  f(x) = 2x
  f(x,y) = 2xy

  numpar = numargs(f) # Should be 2
  g = (x,y) -> x^2
  numpar2 = numargs(g)
  numpar == 2 && numpar2 == 2
end
@test test_num_args()

ext: jl


text:
module MDSim

using Random: AbstractRNG, MersenneTwister
using Distributions: Normal
using LinearAlgebra: diagind
using Statistics: mean

export

MDSimulation,
velocity_verlet!

include("./MDSimulation.jl")
include("./update_methods.jl")
include("./initialization_methods.jl")
include("./computation_methods.jl")

end # module


text:
# implicit scheme for advection-diffusion

using Revise
using AdFem
using ADCME
using PyPlot
using Statistics
using LinearAlgebra
using ADCMEKit

m = 40
n = 20
h = 1/n
NT = 100
Δt = 1/NT
Ra = 1.0

function solve_stokes(T, η)
    η = convert_to_tensor(η)
    z = constant(zeros(4m*n))
    hmat = reshape([η z z z η z z z η], (4m*n,3,3))
    K = compute_fem_stiffness_matrix(hmat, m, n, h)
    B = constant(compute_interaction_matrix(m, n, h))
    Z = [K -B'
    B spzero(size(B,1))]
    bd = bcnode("lower", m, n, h)
    Z, _ = fem_impose_Dirichlet_boundary_condition(Z, bd, m, n, h)
    T = reshape(repeat(T, 1, 4), (-1,))
    rhs = Ra * compute_fem_source_term(constant(zeros(4*m*n)), T, m, n, h)
    rhs = [rhs;constant(zeros(m*n))]
    rhs = scatter_update(rhs, [bd; bd .+ (m+1)*(n+1)], zeros(2length(bd)))
    sol = Z\rhs
    sol[1:2*(m+1)*(n+1)], sol[2*(m+1)*(n+1)+1:end]
end

function solve_heat_eq(u, T)
    up_and_down = bcedge("upper|lower", m, n, h)
    M = compute_fvm_mass_matrix(m, n, h)
    u_on_fvm_grid = [fem_to_fvm(u[1:(m+1)*(n+1)], m, n, h); fem_to_fvm(u[(m+1)*(n+1)+1:end], m, n, h)]
    K, rhs1 = compute_fvm_advection_matrix(u_on_fvm_grid, up_and_down, zeros(size(up_and_down, 1)), m, n, h)
    S, rhs2 = compute_fvm_tpfa_matrix(missing, up_and_down, zeros(size(up_and_down, 1)), m, n, h)
    A = M/Δt + K - S
    T_new = A\(M*T/Δt - rhs1 + rhs2)
end

function compute_viscosity_parameter(u, T; kwargs...)
    # η = mantle_viscosity(u, T, m, n, h; kwargs...)
    η = ones(4*m*n)
end

xy = fvm_nodes(m, n, h)
T0 = @. exp(-10*((xy[:,1]-1.25)^2 + (xy[:,2]-0.5)^2))

u_arr = TensorArray(NT+1)
T_arr = TensorArray(NT+1)
p_arr = TensorArray(NT+1)
i = constant(1, dtype = Int32)
u_arr = write(u_arr, 1, zeros(2*(m+1)*(n+1)))
T_arr = write(T_arr, 1, T0)
p_arr = write(p_arr, 1, zeros(m*n))

function condition(i, tas...)
    i<=NT
end
function body(i, tas...)
    u_arr, T_arr, p_arr = tas
    u, T, p = read(u_arr, i), read(T_arr, i), read(p_arr, i)
    T = solve_heat_eq(u, T)
    η = compute_viscosity_parameter(u, T)
    u, p = solve_stokes(T, η)
    i+1, write(u_arr, i+1, u), write(T_arr, i+1, T), write(p_arr, i+1, p)
end

_, u, T, p = while_loop(condition, body, [i, u_arr, T_arr, p_arr])
u = set_shape(stack(u), (NT+1, 2(m+1)*(n+1)))
T = set_shape(stack(T), (NT+1, m*n))
p = set_shape(stack(p), (NT+1, m*n))

sess = Session(); init(sess)
u_, T_, p_ = run(sess, [u, T, p])

p = visualize_displacement(u_, m, n, h)
saveanim(p, "displacement.gif")
# visualize_scalar_on_fvm_points(T_, m, n, h)
# visualize_scalar_on_fvm_points(p_, m, n, h)

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
@doc raw"""
    calculate_mutual_capacitance(sim::Simulation, ij::Tuple{Int, Int}; consider_multiplicity::Bool = true)

Returns the mutual capacitance between the contacts with ID `i = ij[1]` and `j = ij[2]`.
It is calculated via the weighting potentials of the contacts, ``\Phi_i^w(\vec{r})`` and ``\Phi_j^w(\vec{r})``:
```math
c_{ij} = \epsilon_0 \int_{World} \nabla \Phi_i^w(\vec{r}) ϵ_r(\vec{r}) \nabla \Phi_j^w(\vec{r}) d\vec{r}
```

!!! note
    These are elements of the Mawell Capcitance Matrix. Look up [Capacitances](@ref) for more information.

!!! note
    The electric potential as well as the two weighting potentials of both contacts have to be calculated.

## Arguments
* `sim::Simulation`: [`Simulation`](@ref) for which the capacitance matrix is calculated.
* `ij::Tuple{Int,Int}`: Tuple of indices of the contacts for which the capacitance should be calculated.

## Keywords
* `consider_multiplicity::Bool = true`: Whether symmetries of the system should be taken into account.
    For example, in case of true coaxial detector center around the origin and calculated on a cartesian grid
    with the `x-axis` going from `[0, x_max]` and the `y-axis` going from `[0, y_max]` the multiplicity is 4
    and, if `consider_multiplicity == true`, the returned value is already multiplied by 4.
"""
function calculate_mutual_capacitance(sim::Simulation, ij::Tuple{Int, Int}; consider_multiplicity::Bool = true)
    _calculate_mutual_capacitance(sim.weighting_potentials[ij[1]], sim.weighting_potentials[ij[2]], sim.ϵ_r; consider_multiplicity)
end

function _calculate_mutual_capacitance(
            wpi::WeightingPotential{T,3,CS}, wpj::WeightingPotential{T,3,CS},
            ϵ_r::DielectricDistribution{T,3,CS};
            consider_multiplicity::Bool = true ) where {T, CS}

    int_p1 = interpolated_scalarfield(wpi)
    int_p2 = interpolated_scalarfield(wpj)
    int_ϵ_r = interpolated_scalarfield(ϵ_r)
    cylindrical = CS == Cylindrical
    phi_2D = cylindrical && size(wpi, 2) == size(wpj, 2) == 1
    grid = phi_2D ? get_2π_potential(wpi, n_points_in_φ = 2).grid : _get_closed_potential(wpi).grid
    grid_mps = get_extended_midpoints_grid(grid)

    c_ij::T = _calculate_mutual_capacitance(grid, grid_mps, int_ϵ_r, int_p1, int_p2)

    phi_2D && (c_ij *= 2)
    consider_multiplicity && (c_ij *= multiplicity(grid))
    return uconvert(u"pF", c_ij*u"m" * ϵ0*u"F/m" )
end

function _calculate_mutual_capacitance(grid::Grid{T, 3, CS}, grid_mps, int_ϵ_r, int_p1, int_p2) where {T, CS}
    c_ij_vector = zeros(Float64, size(grid, 3)-1)
    # In some cases, the sum of many Float32's (of different order of magnitudes)
    # can lead to large errors on the sum => Use Float64 as the datatype for the sum.
    @inbounds Base.Threads.@threads for i3 in 1:size(grid, 3)-1
        for i2 in 1:size(grid, 2)-1
            for i1 in 1:size(grid, 1)-1
                w1, w2, w3 = voxel_widths(grid, i1, i2, i3)
                dV = voxel_volume(grid, i1, i2, i3, w1, w2, w3)

                pt_voxel_mid = GridPoint(grid_mps, (i1 + 1, i2 + 1, i3 + 1))
                ϵ_r_voxel = get_interpolation(int_ϵ_r, pt_voxel_mid, CS)

                efs_1 = _approximate_potential_gradient(int_p1, grid, i1, i2, i3, w1, w2, w3)
                efs_2 = _approximate_potential_gradient(int_p2, grid, i1, i2, i3, w1, w2, w3)

                c_ij_vector[i3] += Float64(sum(efs_1 .* efs_2) * dV * ϵ_r_voxel)
            end
        end
    end
    return T(sum(c_ij_vector))
end

function _approximate_potential_gradient(int_p, grid::Grid{T, 3, CS}, i1, i2, i3, w1, w2, w3) where {T, CS}
    p000 = get_interpolation(int_p, GridPoint(grid, (i1    , i2    , i3    )), CS)
    p100 = get_interpolation(int_p, GridPoint(grid, (i1 + 1, i2    , i3    )), CS)
    p010 = get_interpolation(int_p, GridPoint(grid, (i1    , i2 + 1, i3    )), CS)
    p110 = get_interpolation(int_p, GridPoint(grid, (i1 + 1, i2 + 1, i3    )), CS)
    p001 = get_interpolation(int_p, GridPoint(grid, (i1    , i2    , i3 + 1)), CS)
    p101 = get_interpolation(int_p, GridPoint(grid, (i1 + 1, i2    , i3 + 1)), CS)
    p011 = get_interpolation(int_p, GridPoint(grid, (i1    , i2 + 1, i3 + 1)), CS)
    p111 = get_interpolation(int_p, GridPoint(grid, (i1 + 1, i2 + 1, i3 + 1)), CS)
    efv1 = ( (p100 - p000) + (p110 - p010) + (p101 - p001) + (p111 - p011) ) / (4 * w1)
    efv2 = if CS == Cylindrical
        _w2 = (grid.axes[2].ticks[i2 + 1] - grid.axes[2].ticks[i2])
        if i1 == 1
            ((p110 - p100)/(_w2*grid.axes[1].ticks[i1+1]) +
                (p111 - p101)/(_w2*grid.axes[1].ticks[i1+1]) ) / 2
        else
            ((p010 - p000)/(_w2*grid.axes[1].ticks[i1]) +
                (p110 - p100)/(_w2*grid.axes[1].ticks[i1+1]) +
                (p011 - p001)/(_w2*grid.axes[1].ticks[i1]) +
                (p111 - p101)/(_w2*grid.axes[1].ticks[i1+1])) / 4
        end
    else
        ( (p010 - p000) + (p110 - p100) + (p011 - p001) + (p111 - p101) ) / (4 * w2)
    end
    efv3 = ( (p001 - p000) + (p101 - p100) + (p011 - p010) + (p111 - p110) ) / (4 * w3)
    return (efv1, efv2, efv3)
end

@doc raw"""
    calculate_capacitance_matrix(sim::Simulation{T}; consider_multiplicity::Bool = true) where {T}

Calculates the Maxwell Capacitance `N×N`-Matrix in units of pF,
where `N` is the number of contacts of `sim.detector`.
The individual elements, ``c_{i,j}``, are calculated via
[`calculate_mutual_capacitance(sim::Simulation, (i,j)::Tuple{Int,Int})`](@ref).
The matrix should be symmetric. The difference of `C[i,j]` and `C[j,i]` are due
to numerical precision in the integration due to the different grids of the two weighting potentials.

## Arguments
* `sim::Simulation`: [`Simulation`](@ref) for which the capacitance matrix is calculated.

## Keywords
* `consider_multiplicity::Bool = true`: Whether symmetries of the system should be taken into account.
    For example, in case of true coaxial detector center around the origin and calculated on a cartesian grid
    with the `x-axis` going from `[0, x_max]` and the `y-axis` going from `[0, y_max]` the multiplicity is 4
    and, if `consider_multiplicity == true`, the returned value is already multiplied by 4.
"""
function calculate_capacitance_matrix(sim::Simulation{T}; consider_multiplicity::Bool = true) where {T}
    @assert !ismissing(sim.ϵ_r) "The electric potential needs to be calculated first."
    @assert !ismissing(sim.weighting_potentials) "The weighting_potentials needs to be calculated first."
    n = length(sim.weighting_potentials)
    C = zeros(typeof(one(T) * u"pF"), (n, n))
    for i in 1:n
        for j in 1:n
            C[j, i] = if !ismissing(sim.weighting_potentials[i]) && !ismissing(sim.weighting_potentials[j])
                calculate_mutual_capacitance(sim, (i, j); consider_multiplicity)
            else
                missing
            end
        end
    end
    return C
end


text:
module ITensors

using Random,
      Printf,
      LinearAlgebra,
      StaticArrays,
      TimerOutputs

# TODO: move imports to individual files
import Base.adjoint,
       Base.conj,
       Base.convert,
       Base.copy,
       Base.deepcopy,
       Base.copyto!,
       Base.eltype,
       Base.fill!,
       Base.getindex,
       Base.in,
       Base.isapprox,
       Base.isless,
       Base.iterate,
       Base.length,
       Base.push!,
       Base.setindex!,
       Base.eachindex,
       Base.show,
       Base.sum,
       Base.summary,
       Base.similar,
       Base.size,
       Base.ndims,
       Base.!=,
       Base.+,
       Base.-,
       Base.*,
       Base./,
       Base.^,
       Base.setdiff,  # Since setdiff doesn't
                      # work with IndexSet, overload it
       Base.lastindex,
       LinearAlgebra.axpby!,
       LinearAlgebra.axpy!,
       LinearAlgebra.dot,
       LinearAlgebra.norm,
       LinearAlgebra.mul!,
       LinearAlgebra.rmul!,
       LinearAlgebra.normalize!,
       Random.randn!

#####################################
# Global Variables
#
const warnTensorOrder = 10
const GLOBAL_TIMER = TimerOutput()

#####################################
# Index and IndexSet
#
# TODO: load these after Tensor functions
include("smallstring.jl")
include("readwrite.jl")
include("tagset.jl")
include("index.jl")
include("indexset.jl")

#####################################
# Tensor
#
include("tensor/tensor.jl")
include("tensor/tensorstorage.jl")
include("tensor/contraction_logic.jl")
include("tensor/dense.jl")
include("tensor/linearalgebra.jl")
include("tensor/diag.jl")
include("tensor/combiner.jl")
include("tensor/truncate.jl")
include("tensor/svd.jl")

#####################################
# ITensor
#
include("itensor.jl")
include("decomp.jl")
include("iterativesolvers.jl")

###########################################################
# MPS/MPO
#
include("mps/mps.jl")
include("mps/mpo.jl")
include("mps/sweeps.jl")
include("mps/projmpo.jl")
include("mps/observer.jl")
include("mps/dmrg.jl")

###########################################################
# Physics
#
include("physics/tag_types.jl")
include("physics/lattices.jl")
include("physics/site_types/spinhalf.jl")
include("physics/site_types/spinone.jl")
include("physics/site_types/electron.jl")
include("physics/site_types/tj.jl")
include("physics/autompo.jl")

end # module ITensors


text:
# following
#  https://www.mrao.cam.ac.uk/~dag/CUBEHELIX/
#
# published in
#  Green, D. A., 2011, `A colour scheme for the display of astronomical intensity images', Bulletin of the Astronomical Society of India, 39, 289.
#    (2011BASI...39..289G at ADS.)

import Colors: RGB

# parameters: N number of entries in result colormap (btw: Array(RGB{Float64},N))
#  rotations: rotations -> added term
#  start: startpoint in H -> 0.0 -> Red
#  hue: saturation of hue -> 0.0 -> grayscale only, >1.0 -> clipping in some areas
#  gamma: gamma factor

function cubehelix(N::Integer; rotations::Real=-1.5,start::Real=0.5,hue::Real=1.0,gamma::Real=1.0)
    colormap = Array(RGB{Float64},N)
    fract_all = linspace(0,1,N)
    coeff_cos = [+0.14861, -0.29227, +1.97294]
    coeff_sin = [+1.78277, -0.90649, 0.0]
    c = zeros(3)
    for (i,fract) in enumerate(fract_all)
        c[1]=c[2]=c[3]=0

        angle = 2*pi*(((start/3.0)+rotations)*fract)
        fract = fract ^ gamma
        amp = hue*fract*(1-fract)/2.0
        for j in 1:3
            c[j] = clamp(
                fract + amp * ((coeff_cos[j]*cos(angle)) + (coeff_sin[j] * sin(angle))),
                0.0,1.0)
        end
        colormap[i] = RGB{Float64}(c[1],c[2],c[3])
    end
    colormap
end

cubehelix() = cubehelix(64)
access_dict["cubehelix"] = cubehelix

max_stars_repo_name: lobingera/ColorMaps


text:
##############################################################################
#
# Fallback methods, usually overridden for specific distributions
#
##############################################################################

# generic function to get number of samples

nsamples{D<:UnivariateDistribution}(dt::Type{D}, x::Array) = length(x)
nsamples{D<:MultivariateDistribution}(dt::Type{D}, x::Matrix) = size(x, 2)
nsamples{D<:MatrixDistribution,T}(dt::Type{D}, x::Array{Matrix{T}}) = length(x)

#### Statistics ####
mean(d::Distribution) = throw(MethodError(mean,(d,)))
std(d::Distribution) = sqrt(var(d))

# What's the purpose for this function?
function var{M <: Real}(d::UnivariateDistribution, mu::AbstractArray{M})
    V = similar(mu, Float64)
    for i in 1:length(mu)
        V[i] = var(d, mu[i])
    end
    return V
end

function cor(d::MultivariateDistribution)
    R = cov(d)
    m, n = size(R)
    for j in 1:n
        for i in 1:n
            R[i, j] = d.cov[i, j] / sqrt(d.cov[i, i] * d.cov[j, j])
        end
    end
    return R
end

binaryentropy(d::Distribution) = entropy(d) / log(2)

# kurtosis returns excess kurtosis by default
# proper kurtosis can be returned with correction = false
function kurtosis(d::Distribution, correction::Bool)
    if correction
        return kurtosis(d)
    else
        return kurtosis(d) + 3.0
    end
end

excess(d::Distribution) = kurtosis(d)
excess_kurtosis(d::Distribution) = kurtosis(d)
proper_kurtosis(d::Distribution) = kurtosis(d, false)

isplatykurtic(d::Distribution) = kurtosis(d) > 0.0
isleptokurtic(d::Distribution) = kurtosis(d) < 0.0
ismesokurtic(d::Distribution) = kurtosis(d) == 0.0

median(d::UnivariateDistribution) = quantile(d, 0.5)

#### pdf, cdf, and quantile ####

logpdf(d::UnivariateDistribution, x::Real) = log(pdf(d, x))
pdf(d::MultivariateDistribution, x::Vector) = exp(logpdf(d, x))

ccdf(d::UnivariateDistribution, q::Real) = 1.0 - cdf(d, q)
cquantile(d::UnivariateDistribution, p::Real) = quantile(d, 1.0 - p)

logcdf(d::Distribution, q::Real) = log(cdf(d,q))
logccdf(d::Distribution, q::Real) = log(ccdf(d,q))
invlogccdf(d::Distribution, lp::Real) = quantile(d, -expm1(lp))
invlogcdf(d::Distribution, lp::Real) = quantile(d, exp(lp))

#### insupport ####

insupport(d::Distribution, x) = false

function insupport(d::UnivariateDistribution, X::Array)
    for x in X; insupport(d, x) || return false; end
    true
end

function insupport(t::DataType, X::Array)
    for x in X; insupport(t, x) || return false; end
    true
end

function insupport(d::MultivariateDistribution, X::Matrix)
    for i in 1 : size(X, 2)
        if !insupport(d, X[:, i])  # short-circuit is generally faster
            return false
        end
    end
    return true
end

function insupport(d::MatrixDistribution, X::Array)
    for i in 1 : size(X, 3)
        if !insupport(d, X[:, :, i])
            return false
        end
    end
    return true
end

#### log likelihood ####

function loglikelihood(d::UnivariateDistribution, X::Array)
    ll = 0.0
    for i in 1:length(X)
        ll += logpdf(d, X[i])
    end
    return ll
end

function loglikelihood(d::MultivariateDistribution, X::Matrix)
    ll = 0.0
    for i in 1:size(X, 2)
        ll += logpdf(d, X[:, i])
    end
    return ll
end

#### Vectorized functions for univariate distributions ####

for fun in [:pdf, :logpdf, :cdf, :logcdf, :ccdf, :logccdf, :invlogcdf, :invlogccdf, :quantile, :cquantile]
    fun! = symbol(string(fun, '!'))

    @eval begin
        function ($fun!)(r::AbstractArray, d::UnivariateDistribution, X::AbstractArray)
            if length(r) != length(X)
                throw(ArgumentError("Inconsistent array dimensions."))
            end
            for i in 1 : length(X)
                r[i] = ($fun)(d, X[i])
            end
            r
        end

        function ($fun)(d::UnivariateDistribution, X::AbstractArray)
            ($fun!)(Array(Float64, size(X)), d, X)
        end
    end
end

#### Vectorized functions for multivariate distributions ####

function logpdf!(r::AbstractArray, d::MultivariateDistribution, x::AbstractMatrix)
    n::Int = size(x, 2)
    if length(r) != n
        throw(ArgumentError("Inconsistent array dimensions."))
    end
    for i in 1 : n
        r[i] = logpdf(d, x[:, i])
    end
    r
end

function logpdf(d::MultivariateDistribution, X::AbstractMatrix)
    logpdf!(Array(Float64, size(X, 2)), d, X)
end

function pdf!(r::AbstractArray, d::MultivariateDistribution, X::AbstractMatrix)
    logpdf!(r, d, X)  # size checking is done by logpdf!
    for i in 1 : size(X, 2)
        r[i] = exp(r[i])
    end
    r
end

function pdf(d::MultivariateDistribution, X::AbstractMatrix)
    pdf!(Array(Float64, size(X, 2)), d, X)
end

#### logpmf & pmf for discrete distributions ####

logpmf(d::DiscreteDistribution, args::Any...) = logpdf(d, args...)
logpmf!(r::AbstractArray, d::DiscreteDistribution, args::Any...) = logpdf!(r, d, args...)
pmf(d::DiscreteDistribution, args::Any...) = pdf(d, args...)

#### Sampling: rand & rand! ####

function rand!(d::UnivariateDistribution, A::Array)
    for i in 1:length(A)
        A[i] = rand(d)
    end
    return A
end

function rand(d::ContinuousUnivariateDistribution, dims::Dims)
    return rand!(d, Array(Float64, dims))
end

function rand(d::DiscreteUnivariateDistribution, dims::Dims)
    return rand!(d, Array(Int, dims))
end

function rand(d::UnivariateDistribution, dim1::Integer, dims::Integer...)
    return rand(d, map(int, tuple(dim1,dims...)))
end

function rand(d::ContinuousMultivariateDistribution)
    return rand!(d, Array(Float64, dim(d)))
end

function rand(d::DiscreteMultivariateDistribution)
    return rand!(d, Array(Int, dim(d)))
end

function rand(d::ContinuousMultivariateDistribution, n::Integer)
    return rand!(d, Array(Float64, dim(d), n))
end

function rand(d::DiscreteMultivariateDistribution, n::Integer)
    return rand!(d, Array(Int, dim(d), n))
end

function rand(d::MatrixDistribution, n::Integer)
    return rand!(d, Array(Matrix{Float64}, n))
end

function rand!(d::MultivariateDistribution, X::Matrix)
    if size(X, 1) != dim(d)
        error("Inconsistent argument dimensions")
    end
    for i in 1 : size(X, 2)
        X[:, i] = rand(d)
    end
    X
end

function rand!(d::MatrixDistribution, X::Array{Matrix{Float64}})
    for i in 1:length(X)
        X[i] = rand(d)
    end
    return X
end

function sprand(m::Integer, n::Integer, density::Real, d::Distribution)
    return sprand(m, n, density, n -> rand(d, n))
end

# Fitting

function suffstats{D<:Distribution}(dt::Type{D}, xs...)
    argtypes = tuple(D, map(typeof, xs)...)
    error("suffstats is not implemented for $argtypes.")
end

fit_mle{D<:UnivariateDistribution}(dt::Type{D}, x::Array) = fit_mle(D, suffstats(D, x))
fit_mle{D<:UnivariateDistribution}(dt::Type{D}, x::Array, w::Array) = fit_mle(D, suffstats(D, x, w))

fit_mle{D<:MultivariateDistribution}(dt::Type{D}, x::Matrix) = fit_mle(D, suffstats(D, x))
fit_mle{D<:MultivariateDistribution}(dt::Type{D}, x::Matrix, w::Array) = fit_mle(D, suffstats(D, x, w))

fit{D <: Distribution}(dt::Type{D}, x::Array) = fit_mle(D, x)

# Conjugates

posterior{D<:Distribution}(pri::Distribution, G::Type{D}, x) = posterior(pri, suffstats(G, x))
posterior(pri::Distribution, G::GenerativeFormulation, x) = posterior(pri, suffstats(G, x))

posterior{D<:Distribution}(pri::Distribution, G::Type{D}, x, w) = posterior(pri, suffstats(G, x, w))
posterior(pri::Distribution, G::GenerativeFormulation, x, w) = posterior(pri, suffstats(G, x, w))

posterior_rand(pri::Distribution, s::SufficientStats) = rand(posterior(pri, s))
posterior_rand{D<:Distribution}(pri::Distribution, G::Type{D}, x) = rand(posterior(pri, G, x))
posterior_rand{D<:Distribution}(pri::Distribution, G::Type{D}, x, w) = rand(posterior(pri, G, x, w))
posterior_rand(pri::Distribution, G::GenerativeFormulation, x) = rand(posterior(pri, G, x))
posterior_rand(pri::Distribution, G::GenerativeFormulation, x, w) = rand(posterior(pri, G, x, w))

posterior_rand!(r::Array, pri::Distribution, s::SufficientStats) = rand!(posterior(pri, s), r)
posterior_rand!{D<:Distribution}(r::Array, pri::Distribution, G::Type{D}, x) = rand!(posterior(pri, G, x), r)
posterior_rand!{D<:Distribution}(r::Array, pri::Distribution, G::Type{D}, x, w) = rand!(posterior(pri, G, x, w), r)
posterior_rand!(r::Array, pri::Distribution, G::GenerativeFormulation, x) = rand!(posterior(pri, G, x), r)
posterior_rand!(r::Array, pri::Distribution, G::GenerativeFormulation, x, w) = rand!(posterior(pri, G, x, w), r)

posterior_mode(pri::Distribution, s::SufficientStats) = mode(posterior(pri, s))
posterior_mode{D<:Distribution}(pri::Distribution, G::Type{D}, x) = mode(posterior(pri, G, x))
posterior_mode{D<:Distribution}(pri::Distribution, G::Type{D}, x, w) = mode(posterior(pri, G, x, w))
posterior_mode(pri::Distribution, G::GenerativeFormulation, x) = mode(posterior(pri, G, x))
posterior_mode(pri::Distribution, G::GenerativeFormulation, x, w) = mode(posterior(pri, G, x, w))


text:
function solve()
    v, t, s, d = [parse(Int, x) for x in split(readline())]
    if t ≤ d / v ≤ s
        "No"
    else
        "Yes"
    end
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
module SAMPLER
    include("io/util.jl")
    include("io/args.jl")
    include("io/netcdf.jl")
    include("rasterizer/rasterization.jl")

    using NetCDF
    using Printf
    using Base.Filesystem
    using Base
    using SeisSolXDMF
    using ..Args
    using ..Rasterization

    function main()
        println("Using $(Threads.nthreads()) threads.")

        ARGS = read_args()

        #============================================#
        # Read and preprocess args
        #============================================#

        mem_limit       = ARGS["memory-limit"]
        sampling_rate   = ARGS["sampling-rate"]
        domain          = ARGS["domain"]

        has_3d          = !isempty(ARGS["input-file-3d"])
        seafloor_only   = ARGS["seafloor-only"]

        seafloor_vars   = ARGS["seafloor-vars"]
        surface_vars    = ARGS["surface-vars"]
        volumetric_vars = has_3d ? ARGS["volumetric-vars"] : Dict{String,String}()

        water_height    = ARGS["water-height"]
        has_tanioka     = ARGS["tanioka"]

        out_filename    = ARGS["output-file"]
        endswith(out_filename, ".nc") || (out_filename = out_filename * ".nc")
        println("Will write output to \"$out_filename\".")

        surface_output  = !seafloor_only
        volume_output   = !seafloor_only && has_3d

        #============================================#
        # Compare timesteps of 2D and 3D files.
        # They must be equal.
        #============================================#

        xdmf2d = XDMFFile(ARGS["input-file-2d"])
        times  = timesteps_of(xdmf2d)

        if has_3d
            xdmf3d = XDMFFile(ARGS["input-file-3d"])
            times_3d = timesteps_of(xdmf3d)

            if (times_3d != times)
                throw(ArgumentError("Timesteps of 2D and 3D input files do not match!"))
            end

            xdmf3d = nothing
            times_3d = nothing
        end

        #============================================#
        # Compute timesteps to be processed from args
        #============================================#

        t_start = 1
        t_end   = length(times)

        if !isnothing(ARGS["output-time"])
            out_time = ARGS["output-time"]
            while t_start < length(times) && times[t_start + 1] ≤ out_time.t_start
                t_start += 1
            end

            while t_end > t_start && times[t_end - 1] ≥ out_time.t_end
                t_end -= 1
            end
        elseif !isnothing(ARGS["output-steps"])
            out_steps = ARGS["output-steps"]
            t_start   = min(1 + out_steps.t_start, length(times))
            t_end     = min(1 + out_steps.t_end, length(times))
        end

        #============================================#
        # Process 2D seafloor
        #============================================#

        empty_mapping    = Args.VarMapping()
        stat_var_mapping = Args.VarMapping()
        # "b" is the only static variable. If it is found in the mapping, it is separated into the stat_var_mapping.
        if haskey(seafloor_vars, "b");
            stat_var_mapping["b"] = seafloor_vars["b"]
            delete!(seafloor_vars, "b")
        end

        Rasterization.rasterize(xdmf2d, seafloor_vars, stat_var_mapping, times, sampling_rate, out_filename, mem_limit;
                                create_nc_dyn_vars=[seafloor_vars, surface_vars, volumetric_vars],
                                create_nc_stat_vars=[stat_var_mapping],
                                z_range=Rasterization.z_floor, t_begin=t_start, t_end=t_end, water_height=water_height,
                                tanioka=has_tanioka, domain=domain)

        GC.gc(true)

        #============================================#
        # Process 2D sea surface
        #============================================#

        if surface_output
            Rasterization.rasterize(xdmf2d, surface_vars, empty_mapping, times, sampling_rate, out_filename, mem_limit;
                                    z_range=Rasterization.z_surface, t_begin=t_start, t_end=t_end, water_height=water_height,
                                    domain=domain)
        end

        xdmf2d = nothing
        GC.gc(true)

        #============================================#
        # Process 3D mesh
        #============================================#

        if volume_output
            xdmf3d = XDMFFile(ARGS["input-file-3d"])

            Rasterization.rasterize(xdmf3d, volumetric_vars, empty_mapping, times, sampling_rate, out_filename, mem_limit;
                                    z_range=Rasterization.z_all, t_begin=t_start, t_end=t_end, water_height=water_height,
                                    domain=domain)
        end
    end
end

max_stars_repo_name: SeisSol/SAMPLER

max_issues_repo_name: SeisSol/SAMPLER

max_forks_repo_name: SeisSol/SAMPLER


text:
using Test

@testset "Finite FK Model tests" begin
  include("ffk_test.jl")
end

@testset "LGModel tests" begin
  include("lg_test.jl")
end

@testset "MVLGModel tests" begin
  include("mvlg_test.jl")
end

@testset "SMC Sampler tests" begin
  include("smcs_test.jl")
end

@testset "Lorenz96 tests" begin
  include("lorenz96_test.jl")
end

@testset "Netto tests" begin
  include("netto_test.jl")
end

@testset "PMMH tests" begin
  include("lgPMMH_test.jl")
end

@testset "Stochastic Volatility tests" begin
  include("sv_test.jl")
end


text:
using TopOpt, Test, LinearAlgebra

# write your own tests here
#include("metadata_tests.jl")
include("csimp.jl")


text:
#!/usr/bin/env julia

# CUDAnative.jl version

using CUDAdrv, CUDAnative

using Statistics
using Printf

function kernel_dummy(ptr)
    Base.pointerset(ptr, 0f0, Int(blockIdx().x), 8)
    return
end

const len = 1000
const ITERATIONS = 100

function benchmark(gpu_buf)
    @cuda threads=len kernel_dummy(Base.unsafe_convert(Ptr{Float32}, gpu_buf))
end

function main()
    cpu_time = Vector{Float64}(undef, ITERATIONS)
    gpu_time = Vector{Float64}(undef, ITERATIONS)

    gpu_buf = Mem.alloc(len*sizeof(Float32))
    for i in 1:ITERATIONS
        i == ITERATIONS-4 && CUDAdrv.Profile.start()

        gpu_tic, gpu_toc = CuEvent(), CuEvent()

        cpu_tic = time_ns()
        record(gpu_tic)
        benchmark(gpu_buf)
        record(gpu_toc)
        synchronize(gpu_toc)
        cpu_toc = time_ns()

        cpu_time[i] = (cpu_toc-cpu_tic)/1000
        gpu_time[i] = CUDAdrv.elapsed(gpu_tic, gpu_toc)*1000000
    end
    CUDAdrv.Profile.stop()
    Mem.free(gpu_buf)

    popfirst!(cpu_time)
    popfirst!(gpu_time)

    @printf("CPU time: %.2f ± %.2f us\n", mean(cpu_time), std(cpu_time))
    @printf("GPU time: %.2f ± %.2f us\n", mean(gpu_time), std(gpu_time))
end

main()


text:
a = [1 2  3  4
     5 6  7  8
     9 10 11 12]
da = DimensionalArray(a, (Y(10:30), Time((1:4)u"s")))

@testset "selectors with dim wrappers" begin
    @test da[Y<|At([10, 30]), Time<|At([1u"s", 4u"s"])] == [1 4; 9 12]
    @test_throws ArgumentError da[Y<|At([9, 30]), Time<|At([1u"s", 4u"s"])]
    @test view(da, Y<|At(20), Time<|At((3:4)u"s")) == [7, 8]
    @test view(da, Y<|Near(17), Time<|Near([1.3u"s", 3.3u"s"])) == [5, 7]
    @test view(da, Y<|Between(9, 31), Time<|At((3:4)u"s")) == [3 4; 7 8; 11 12]
end

@testset "selectors without dim wrappers" begin
    @test da[At(20:10:30), At(1u"s")] == [5, 9]
    @test view(da, Between(9, 31), Near((3:4)u"s")) == [3 4; 7 8; 11 12]
    @test view(da, Near(22), At(3u"s", 4u"s")) == [7, 8]
    @test view(da, At(20), At((2:3)u"s")) == [6, 7]
    @test view(da, Near<|13, Near<|[1.3u"s", 3.3u"s"]) == [1, 3]
    # Near works with a tuple input
    @test view(da, Near<|(13,), Near<|(1.3u"s", 3.3u"s")) == [1 3]
    @test view(da, Between(11, 20), At((2:3)u"s")) == [6 7]
    # Between also accepts a tuple input
    @test view(da, Between((11, 20)), Between((2u"s", 3u"s"))) == [6 7]
end

@testset "setindex! with selectors" begin
    da[Near(11), At(3u"s")] = 100
    @test a[1, 3] == 100
    da[Time<|Near(2.2u"s"), Y<|Between(10, 30)] = [200, 201, 202]
    @test a[1:3, 2] == [200, 201, 202]
end

a = [1 2  3  4
     5 6  7  8
     9 10 11 12]

@testset "more Unitful dims" begin
    dimz = Time<|1.0u"s":1.0u"s":3.0u"s", Y<|(1u"km", 4u"km")
    da = DimensionalArray(a, dimz)
    @test da[Y<|Between(2u"km", 3.9u"km"), Time<|At<|3.0u"s"] == [10, 11]
end

@testset "CategoricalGrid" begin
    dimz = Time([:one, :two, :three]; grid=CategoricalGrid()),
           Y([:a, :b, :c, :d]; grid=CategoricalGrid())
    da = DimensionalArray(a, dimz)
    @test da[Time<|At(:one, :two), Y<|At(:b)] == [2, 6]
    @test da[At([:one, :three]), At([:b, :c, :d])] == [2 3 4; 10 11 12]
    @test da[At(:two), Between(:b, :d)] == [6, 7, 8]
    # Near doesn't make sense for categories
    @test_throws ArgumentError da[Near(:two), At([:b, :c, :d])]
end

@testset "TranformedGrid " begin
    using CoordinateTransformations

    m = LinearMap([0.5 0.0; 0.0 0.5])

    dimz = Dim{:trans1}(m; grid=TransformedGrid(X())),
           Dim{:trans2}(m, grid=TransformedGrid(Y()))

    @testset "permutedims works on grid dimensions" begin
        @test permutedims((Y(), X()), dimz) == (X(), Y())
    end

    a = [1 2  3  4
         5 6  7  8
         9 10 11 12]
    da = DimensionalArray(a, dimz)

    @testset "Indexing with array dims indexes the array as usual" begin
        @test da[Dim{:trans1}(3), Dim{:trans2}(1)] == 9
        # Using selectors works the same as indexing with grid
        # dims - it applies the transform function.
        # It's not clear this should be allowed or makes sense,
        # but it works anyway because the permutatoin is correct either way.
        @test da[Dim{:trans1}(At(6)), Dim{:trans2}(At(2))] == 9
    end

    @testset "Indexing with grid dims uses the transformation" begin
        @test da[X(Near(6.1)), Y(Near(8.5))] == 12
        @test da[X(At(4.0)), Y(At(2.0))] == 5
        @test_throws InexactError da[X(At(6.1)), Y(At(8))]
        # Indexing directly with grid dims also just works, but maybe shouldn't?
        @test da[X(2), Y(2)] == 6
    end
end


text:
export PDNode

mutable struct PDNode
    data
    children::Vector{PDNode}
    parent
end

PDNode(data) = PDNode(data, PDNode[], nothing)
function PDNode(data, children::Vector{PDNode})
    node = PDNode(data, children, nothing)
    push!(node, children...)
    node
end

isroot(node::PDNode) = node.parent == nothing
Base.isempty(node::PDNode) = isempty(node.children)
Base.length(node::PDNode) = length(node.children)
Base.size(node::PDNode, i::Int) = size(node.children, i)
#Base.endof(node::PDNode) = endof(node.children)

function Base.setindex!(node::PDNode, child::PDNode, i::Int)
    node.children[i] = child
    child.parent = node
end

function Base.push!(node::PDNode, children::PDNode...)
    push!(node.children, children...)
    for i = length(children):-1:1
        c = children[i]
        delete!(c)
        c.parent = node
    end
end

function Base.insert!(node::PDNode, i::Int, child::PDNode)
    insert!(node.children, i, child)
    child.parent = node
end

function Base.delete!(node::PDNode)
    node.parent == nothing && return
    i = findfirst(x -> x == node, node.parent.children)
    deleteat!(node.parent, i)
end

function topdown(f, node::PDNode)
    f(node)
    for c in node.children
        topdown(f, c)
    end
end

function bottomup(f, node::PDNode)
    for c in node.children
        bottomup(f, c)
    end
    f(node)
end

function toconll(root::PDNode)
    dict = Dict()
    i = 1
    bottomup(root) do node
        if isempty(node)
            dict[node] = i:i
            i += 1
        else
            r1 = dict[node[1]]
            r2 = dict[node[end]]
            dict[node] = first(r1):last(r2)
        end
    end
end


text:
@testset "Flux utilities" begin
    model = load_model(StandardModel, model_paths["e_coli_core.json"])

    fluxes = flux_balance_analysis_dict(
        model,
        Tulip.Optimizer;
        modifications = [change_objective("BIOMASS_Ecoli_core_w_GAM")],
    )

    consuming, producing = metabolite_fluxes(model, fluxes)
    @test isapprox(consuming["atp_c"]["PFK"], -7.47738; atol = TEST_TOLERANCE)
    @test isapprox(producing["atp_c"]["PYK"], 1.75818; atol = TEST_TOLERANCE)

    # remove the biomass production from the fluxes, so that there's some atom
    # disbalance that can be measured
    delete!(fluxes, "BIOMASS_Ecoli_core_w_GAM")

    # atom tracker
    atom_fluxes_out = atom_fluxes(model, fluxes)
    @test isapprox(atom_fluxes_out["C"], 37.190166489763214; atol = TEST_TOLERANCE)
    @test isapprox(atom_fluxes_out["O"], 41.663071522672226; atol = TEST_TOLERANCE)
    @test isapprox(atom_fluxes_out["N"], 4.765319167566247; atol = TEST_TOLERANCE)
end


text:
WeightedToSimple(graph) = LightGraphs.Graph( SimpleWeightedGraphs.adjacency_matrix( graph ) )

function countitems(x)
    un = unique(x)
    counts = [ sum(x .== i) for i in un ]
    return un .=> counts
end

function EmpiricalFormula(data::Array{Element, 1} )
    countem = Dict(countitems( abbreviation.( data ) ))
    nhydrogens = sum(H, data)
    if nhydrogens > 0
        countem["H"] = get(countem, "H", 0) + nhydrogens
    end
    sortem = sort( collect( keys( countem ) ) )
    return join([ (countem[ele] == 1) ? ele : ele * string( countem[ele] ) for ele in sortem])
end

"""
    gH, atomnodesH = instantiate_hydrogens(g, atomnodes)

Given a molecule, where `g` is the connectivity graph and `atomnodes` contains additional
information about the vertices, return a new representation `gH, atomnodesH` where all hydrogens
in `g`, `atomnodes` have been instantiated as full nodes in `gH, atomnodesH`.
"""
function instantiate_hydrogens(g, atomnodes::AbstractVector{Element})
    gH, atomnodesH = copy(g), [copy(atom) for atom in atomnodes]
    nhydrogens = sum(H, atomnodes)
    nhydrogens == 0 && return gH, atomnodesH
    lastH = SimpleWeightedGraphs.nv(g)
    SimpleWeightedGraphs.add_vertices!(gH, nhydrogens)
    for (i, atom) in enumerate(atomnodesH)
        for _ = 1:H(atom)
            SimpleWeightedGraphs.add_edge!(gH, i, lastH += 1)
        end
        atom.explicithydrogens = atom.implicithydrogens = 0
    end
    append!(atomnodesH, [Element("H") for _ in 1:nhydrogens])
    return gH, atomnodesH
end


text:
include("monkey/monkey.jl")
using .Monkey

data = String(read("program.monkey"))
lexer = Lexer(data)
parser = Parser(lexer)

prog = Monkey.parse(parser)
println(prog)

max_stars_repo_name: viert/monkeylang-julia

max_issues_repo_name: viert/monkeylang-julia

max_forks_repo_name: viert/monkeylang-julia


text:
# Further tests from Szabo and Ostlund, Table $3.10, Standard geoms used in calcs
# units are 'a.u.' ; so assume Bohr
CO = Quante.Molecule([  Atom(AtoZ["C"], 0, 0, 0),
                        Atom(AtoZ["O"], 2.132, 0, 0)])
N2 = Quante.Molecule([  Atom(AtoZ["N"], 0, 0, 0),
                        Atom(AtoZ["N"], 2.074, 0, 0)])

@testset "SzaboOstlundReference" begin

# Szabo and Ostlund
@time Energy,E,U = rhf(CO)
println("CO Energy: $Energy S&O STO-3G: -111.225")
@test Energy ≈ -111.225 atol=2e-3

@time Energy,E,U = rhf(N2)
println("N2 Energy: $Energy S&O STO-3G: -107.496")
#@test Energy ≈ -107.496 atol=2e-3

end #testset


text:
using FractionalDiffEq, Documenter, DemoCards, JSON

ChaosGallery, postprocess_cb, ChaosGallery_assets = makedemos("ChaosGallery")

makedocs(;
    modules=[FractionalDiffEq],
    authors="Qingyu Qu",
    repo="https://github.com/SciFracX/FractionalDiffEq.jl/blob/{commit}{path}#{line}",
    sitename="FractionalDiffEq.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://SciFracX.github.io/FractionalDiffEq.jl",
        assets=["assets/favicon.ico"],
    ),
    pages = [
        "FractionalDiffEq.jl" => "index.md",
        "Get Started" => "get_started.md",
        "Multi-term Fractional Differential Equations" => "multiterm.md",
        "System of Fractional Differential Equations" => "fdesystem.md",
        "Fractional Partial Differential Equations" => "fpde.md",
        "Fractional Delay Differential Equations" => "fdde.md",
        "Distributed Order Differential Equations" => "dode.md",
        "Fractional Differences Equations" => "differences.md",
        "Fractional Integral Equations" => "fie.md",
        "Ordinary Differential Equations" => "ode.md",
        ChaosGallery,
        "Problems" => "problems.md",
        "Algorithms" => "algorithms.md",
        "Mittag Leffler Function" => "mittagleffler.md",
        "Comparison" => "comparison.md"
    ],
)

postprocess_cb()

deploydocs(;
    repo="github.com/SciFracX/FractionalDiffEq.jl",
)


text:
#=
# Test out different node orderings
=#
if !@isdefined(Femshop)
    include("../Femshop.jl");
    using .Femshop
end
n = 6;
ord = 7;

init_femshop("efordering");
@useLog("eforderinglog")

@domain(3)
@functionSpace(LEGENDRE, ord)

@matrixFree(200,1e-6)

@mesh(HEXMESH, n)

@variable(u)

@testSymbol(v)

@boundary(u, 1, DIRICHLET, 0)

@coefficient(f, "-3*pi*pi*sin(pi*x)*sin(pi*y)*sin(pi*z)")
@weakForm(u, "-dot(grad(u), grad(v)) - f*v")

#
#using Plots
#pyplot();

# exact solution is sin(pi*x)*sin(pi*y)*sin(pi*z)
# check error
function check_error()
    maxerr = 0;
    exact(x,y,z) = sin(pi*x)*sin(pi*y)*sin(pi*z);

    for i=1:size(Femshop.grid_data.allnodes,1)
        x = Femshop.grid_data.allnodes[i,1];
        y = Femshop.grid_data.allnodes[i,2];
        z = Femshop.grid_data.allnodes[i,3];
        err = abs(u.values[i] - exact(x,y,z));
        tmp = maxerr;
        maxerr = max(err,maxerr);
    end
    println("max error for u = "*string(maxerr));
end

# First do u (one DOF)
times = 1; # average times times
println("One DOF per node, average of "*string(times)*" runs")
#warm up
solve(u);

# Lex. ordering
regtime = Base.Libc.time();
for iter=1:times
    solve(u);
end
regtime = Base.Libc.time() - regtime;
regtime /= times;
println("Lex. time = "*string(regtime)*"sec.");
check_error();
#display(spy(Femshop.CGSolver.Amat, legend=nothing, reuse=false));

# morton ordering
@mesh(HEXMESH, n)
hilbert_elements([n,n,n]);
ef_nodes();

#solve(u);
tiletime = Base.Libc.time();
for iter=1:times
    solve(u);
end
tiletime = Base.Libc.time() - tiletime;
tiletime /= times;

println("ef = "*string(tiletime)*"sec.");
check_error()
#display(spy(Femshop.CGSolver.Amat, legend=nothing, reuse=false));

@finalize()


text:
using Qsosed
using Test

@testset "Test Black Hole properties" begin
    black_hole = BlackHole(1e8 * M_SUN, 0.5, 0.5)
    @test compute_eddington_luminosity(black_hole) ≈ 1.25750e46 * 1.17 rtol = 1e-3 #erg/s
    @test compute_bolometric_luminosity(black_hole) ≈
          0.5 * compute_eddington_luminosity(black_hole) rtol = 1e-3
    #@test mass_accretion_rate(black_hole) ≈ 1.39916e26 rtol = 1e-4
    @test black_hole.Rg ≈ 14766250380501.246 rtol = 1e-4
    @test black_hole.isco ≈ 4.2330025 rtol = 1e-4
    black_hole2 = BlackHole(1e8 * M_SUN, 0.5, 0.0)
    @test black_hole2.isco ≈ 6.0 rtol = 1e-4
    @test black_hole2.efficiency ≈ 0.057190958417936644
    @test compute_mass_accretion_rate(black_hole2) ≈ 1.2228101097715833e+26 * 1.17
end

@testset "Initialise black holes" begin
    parameters = Parameters(M=1e9, mdot=0.2, spin=0.6)
    bh = BlackHole(parameters)
    @test bh.M == 1e9 * M_SUN
    @test bh.mdot == 0.2
    @test bh.spin == 0.6
end


text:
using FractionalSystems, ControlSystems
using Test

@testset "Test FOSS" begin
    s = foss([-5 0; 0 -5], [2; 2], [3 3],[0], 0.5, 2, [2], 2)

    @test foss2ss(s)  == ss([-5 0; 0 -5], [2; 2], [3 3], [0])
end


text:
"""
    AbstractPointProcess{T}

Abstract type encoding the notion of [point process](https://en.wikipedia.org/wiki/Point_process)
The type `T` corresponds to the elements' type of the point process (vectors, graph edges, etc.)

**See also**
- [`PRS.AbstractSpatialPointProcess`](@ref)
- [`PRS.AbstractGraphPointProcess`](@ref)
"""
abstract type AbstractPointProcess{T} end
Base.eltype(pp::AbstractPointProcess{T}) where {T} = T

@doc raw"""
    AbstractGraphPointProcess{T<:Vector{Float64}} <: AbstractPointProcess{T}

Abstract type encoding point processes defined on graphs
"""
abstract type AbstractGraphPointProcess{T} <: AbstractPointProcess{T} end

@doc raw"""
    AbstractSpatialPointProcess{T<:Vector{Float64}} <: AbstractPointProcess{T}

Abstract type encoding point processes defined on ``\mathbb{R}^d``.

Concrete instances must have a `window` field of type [`PRS.AbstractSpatialWindow`](@ref)
"""
abstract type AbstractSpatialPointProcess{T<:Vector{Float64}} <: AbstractPointProcess{T} end

"""
    window(pp::AbstractSpatialPointProcess)::AbstractSpatialWindow = pp.window
"""
function window(pp::AbstractSpatialPointProcess)::AbstractSpatialWindow
    return pp.window
end

"""
    dimension(pp::AbstractSpatialPointProcess) = dimension(window(pp))
"""
dimension(pp::AbstractSpatialPointProcess) = dimension(window(pp))

"""
Abstract type representing a window

**See also**

- [`PRS.AbstractDiscreteWindow`](@ref)
- [`PRS.AbstractSpatialWindow`](@ref)
"""
abstract type AbstractWindow end

@doc raw"""
    AbstractSpatialWindow{T<:Float64} <: AbstractWindow{T}

Abstract type representing a spatial window ``\subseteq \mathbb{R}^d``

**See also**

- [`PRS.AbstractRectangleWindow`](@ref)
    - [`PRS.RectangleWindow`](@ref)
    - [`PRS.SquareWindow`](@ref)
- [`PRS.BallWindow`](@ref)
"""
abstract type AbstractSpatialWindow{T<:Float64} <: AbstractWindow end

"""
Abstract type representing a window on a discrete state space

**See also**

- [`PRS.GraphNode`](@ref)
"""
abstract type AbstractDiscreteWindow{T} <: AbstractWindow end

function generate_sample(pp::AbstractPointProcess, args...; kwargs...)
    return generate_sample(Random.default_rng(), pp, args...; kwargs...)
end

function generate_sample_prs(pp::AbstractPointProcess, args...; kwargs...)
    return generate_sample_prs(Random.default_rng(), pp, args...; kwargs...)
end


text:
__precompile__()
module RocksDB

# The author(s) would like to acknowledge jerryzhenleicai and others who
# contributed to LevelDB.jl
# https://github.com/jerryzhenleicai/LevelDB.jl

using Pkg

#=
depsfile = joinpath(pathof(RocksDB), "deps", "deps.jl")
if isfile(depsfile)
    include(depsfile)
else
    error("RocksDB not properly installed. Please run Pkg.build(\"RocksDB\")")
end
=#

const depsfile = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(depsfile)
    include(depsfile)
else
    error("RocksDB not properly installed. Please run Pkg.build(\"RocksDB\") then restart Julia.")
end

export open_db, close_db
export create_write_batch, batch_put, write_batch
export db_put, db_get, db_delete
export db_range, range_close, db_delete_range, db_compact_range
export db_snap_key_range, db_create_snapshot, db_release_snapshot
export db_backup_open, db_backup_create, db_backup_close
export db_backup_purge, db_backup_restore
export db_create_checkpoint

struct RocksDBException <: Exception
    msg::String
end

include("serialize.jl")
include("interface.jl")
end # module


text:
using JuliaTutorial
using Test

# to group test together, to oranise serveral tests
# @testset "JuliaTutorial.jl" begin
#     # @test myaddition(1, 1) == 2
# end

sexonthebeach = Cocktail(("vodka", "schnaps", "orange juice", "cranberry juice"))
tipunch = Cocktail(("rhum", "sugar", "lime juice"))

olivierssecret = myaddition(sexonthebeach, tipunch)

@test olivierssecret isa Cocktail


text:
export evolve

"""
    *( U :: Unitary, ket :: Ket ) :: Ket
    *( bra :: Bra, U :: Unitary ) :: Bra
    *(U1 :: Unitary, ρ :: State, U2 ::  Unitary) :: Matrix
"""
*(U :: Unitary, ket :: Ket) :: Ket = Ket(U.M * ket.ψ, atol=max(U.atol, ket.atol))
*(bra :: Bra, U :: Unitary) :: Bra = Bra(bra.ψ * U.M, atol=max(U.atol, bra.atol))
*(U1 :: Unitary, ρ :: State, U2 ::  Unitary) :: Matrix =  U1.M*ρ.M*U2.M

"""
Apply a unitary evolution `U` to density matrix `ρ`: ``\\rho' = U\\rho U^{\\dagger}``.

    evolve(
        U :: Unitary,
        ρ :: State
    ) :: State
"""
evolve(U :: Unitary, ρ :: State) :: State = State(
    U*ρ*U', atol=max(U.atol, ρ.atol)
)
evolve(U :: AbstractMatrix, ρ :: AbstractMatrix; atol=ATOL :: Float64) :: State = State(U*ρ*U', atol=atol)
"""
Apply a unitary evolution `U` to a state ket `ψ`: ``|\\psi'\\rangle = U|\\psi\\rangle``.

    evolve(
        U ::  Unitary,
        ψ :: Ket
    ) :: Ket
"""
evolve(U::Unitary, ψ::Ket) :: Ket = U*ψ
evolve(U :: AbstractMatrix, ψ :: AbstractVector; atol=ATOL :: Float64) :: Ket = Ket(U*ψ, atol=atol)


text:
export get_size, set_size, get_edge, set_edge, has_directed_path, has_cycle
export is_DAG, to_DAG, set, to_network, has_incoming_edge, has_outgoing_edge
export print, get_adjacent_nodes, get_parents, get_children
export DSL_EDGETYPE_NONE, DSL_EDGETYPE_UNDIRECTED, DSL_EDGETYPE_DIRECTED

# ------

const DSL_EDGETYPE_NONE       = unsafe_load(cglobal((:EDSL_EDGETYPE_NONE,       LIB_SMILE), Cint))
const DSL_EDGETYPE_UNDIRECTED = unsafe_load(cglobal((:EDSL_EDGETYPE_UNDIRECTED, LIB_SMILE), Cint))
const DSL_EDGETYPE_DIRECTED   = unsafe_load(cglobal((:EDSL_EDGETYPE_DIRECTED,   LIB_SMILE), Cint))

# ------

get_size( pat::Pattern ) = ccall( (:pattern_getSize, LIB_SMILE), Cint, (Ptr{Void},), pat.ptr )
set_size( pat::Pattern, size::Integer ) = ccall( (:pattern_setSize, LIB_SMILE), Void, (Ptr{Void},Cint), pat.ptr, size )
get_edge( pat::Pattern, from::Integer, to::Integer ) = ccall( (:pattern_getEdge, LIB_SMILE), Cint, (Ptr{Void},Cint,Cint), pat.ptr, from, to )
set_edge( pat::Pattern, from::Integer, to::Integer, edgetype::Cint ) = ccall( (:pattern_setEdge, LIB_SMILE), Void, (Ptr{Void},Cint,Cint,Cint), pat.ptr, from, to, edgetype )
has_directed_path( pat::Pattern, from::Integer, to::Integer ) = ccall( (:pattern_hasDirectedPath, LIB_SMILE), Bool, (Ptr{Void},Cint,Cint), pat.ptr, from, to )
has_cycle( pat::Pattern ) = ccall( (:pattern_hasCycle, LIB_SMILE), Bool, (Ptr{Void},), pat.ptr )
is_DAG( pat::Pattern ) = ccall( (:pattern_isDAG, LIB_SMILE), Bool, (Ptr{Void},), pat.ptr )
to_DAG( pat::Pattern ) = ccall( (:pattern_toDAG, LIB_SMILE), Bool, (Ptr{Void},), pat.ptr )
set( pat::Pattern, net::Network ) = ccall( (:pattern_set, LIB_SMILE), Void, (Ptr{Void},Ptr{Void}), pat.ptr, net.ptr )
to_network( pat::Pattern, dset::Dataset, net::Network ) = ccall( (:pattern_toNetwork, LIB_SMILE), Bool, (Ptr{Void},Ptr{Void},Ptr{Void}), pat.ptr, dset.ptr, net.ptr )
has_incoming_edge( pat::Pattern, to::Integer ) = ccall( (:pattern_hasIncomingEdge, LIB_SMILE), Bool, (Ptr{Void},Cint), pat.ptr, to )
has_outgoing_edge( pat::Pattern, from::Integer ) = ccall( (:pattern_hasOutgoingEdge, LIB_SMILE), Bool, (Ptr{Void},Cint), pat.ptr, from )
print( pat::Pattern ) = ccall( (:pattern_print, LIB_SMILE), Void, (Ptr{Void},), pat.ptr )

function get_adjacent_nodes( pat::Pattern, node::Integer )
	size = get_size(pat)
	arr  = Array(Cint, size)
	len  = ccall( (:pattern_getAdjacentNodes, LIB_SMILE), Cint, (Ptr{Void},Cint,Ptr{Cint}), pat.ptr, node, arr )
	arr[1:len]
end
function get_parents( pat::Pattern, node::Integer )
	size = get_size(pat)
	arr  = Array(Cint, size)
	len  = ccall( (:pattern_getParents, LIB_SMILE), Cint, (Ptr{Void},Cint,Ptr{Cint}), pat.ptr, node, arr )
	arr[1:len]
end
function get_children( pat::Pattern, node::Integer)
	size = get_size(pat)
	arr  = Array(Cint, size)
	len  = ccall( (:pattern_getChildren, LIB_SMILE), Cint, (Ptr{Void},Cint,Ptr{Cint}), pat.ptr, node, arr )
	arr[1:len]
end


text:
using SimpleDirectMediaLayer

# import SimpleDirectMediaLayer:
#     Window, Renderer

using ..Ranger:
    World

const SDL2 = SimpleDirectMediaLayer

const SAMPLESIZE = 16

mutable struct SDL

    function SDL(world::World)
        # SDL2.GL_SetAttribute(SDL2.GL_MULTISAMPLEBUFFERS, SAMPLESIZE)
        # SDL2.GL_SetAttribute(SDL2.GL_MULTISAMPLESAMPLES, SAMPLESIZE)

        SDL2.Init(UInt32(SDL2.INIT_VIDEO))

        if world.window_centered
            world.window = SDL2.CreateWindow(world.title, Int32(SDL2.WINDOWPOS_CENTERED()), Int32(SDL2.WINDOWPOS_CENTERED()),
                Int32(world.window_width), Int32(world.window_height),
                UInt32(SDL2.WINDOW_SHOWN))
        else
            world.window = SDL2.CreateWindow(world.title, Int32(world.window_position_x), Int32(world.window_position_y),
                Int32(world.window_width), Int32(world.window_height),
                UInt32(SDL2.WINDOW_SHOWN))
        end

        SDL2.SetWindowResizable(world.window, false)

        world.renderer = SDL2.CreateRenderer(world.window, Int32(-1),
            UInt32(SDL2.RENDERER_ACCELERATED | SDL2.RENDERER_PRESENTVSYNC))

        new()
    end
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
using IntegratedOptics
using Test

@testset "IntegratedOptics.jl" begin
    # Write your tests here.
end


text:
"""
    web_socket_server(port::UInt16, io_port::Integer, start::Condition [; ipaddr::String])

Starts a bi-directional web socket server that uses a WEB-socket
at address `ipaddr::String` (defaults to `"127.0.0.1"`) and port
`port` and communicates with the search server through the
TCP port `io_port`. The server is started once the condition
`start` is triggered.
"""
function web_socket_server(port::UInt16,
                           io_port::Integer,
                           start::Condition;
                           ipaddr::String="127.0.0.1")
    #Checks
    if port <= 0
        @error "Please specify a WEB-socket port of positive integer value."
    end

    # Wait for search server to be ready
    wait(start)
    @info "Web-Socket server online @$ipaddr:$port..."

    # Start serving requests
    @async HTTP.WebSockets.listen(ipaddr, port, verbose=false) do ws
        while !eof(ws)
            # Read data
            data = readavailable(ws)
            # Convert to String
            request = join(Char.(data))
            if !isempty(request)
                # Send request to search server and get response
                ssconn = connect(Sockets.localhost, io_port)
                println(ssconn, request)                                 # writes a "\n" as well
                response = ifelse(isopen(ssconn), readline(ssconn), "")  # expects a "\n" terminator
                close(ssconn)
                # Return response
                write(ws, response)
            end
        end
    end
end


text:
abstract type BoundaryCondition end

"""
Apply a group of boundary conditions after streaming
"""
function apply!(
    boundary_conditions::BCs,
    q::Quadrature,
    f_new,
    f_old;
    time = 0.0,
) where {BCs <: AbstractVector{<:BoundaryCondition}}
    for boundary_condition in boundary_conditions
        apply!(boundary_condition, q, f_new, f_old, time = time)
    end
end

function apply!(bc::BoundaryCondition, q::Quadrature, f_new, f_old; time = 0.0)
    apply!(bc::BoundaryCondition, q::Quadrature, f_new, f_old)
end

"""
Often we need special rules when a boundary is placed in a specific direction

Currently we only support a few 2D boundary conditions
In the future we could add additional directions to support specific corners
"""
abstract type Direction end
struct North <: Direction end
struct East <: Direction end
struct South <: Direction end
struct West <: Direction end

include("boundary_conditions/bounce_back.jl")
include("boundary_conditions/moving_wall.jl")


text:
#=
renormalize!
=#

function renormalize!{T}(lattice::Classical2dFractalLattice{T}, dimM::Int)
	tensorT, tensorTtilde = getTensorT(lattice)

	tensorT, tensorTtilde, legTensorY, caretTensorY =
	constructHalf(tensorT, tensorTtilde, lattice, dimM, 2)

	tensorS, tensorStilde, legTensorX, caretTensorX =
	constructHalf(tensorT, tensorTtilde, lattice, dimM, 1)

	updateLocalTensors!(tensorS, tensorStilde, legTensorX, legTensorY, caretTensorX, caretTensorY, lattice)
end

function updateLocalTensors!{T}(tensorS::Array{T,4}, tensorStilde::Array{T,4}, legTensorX::Array{T,2}, legTensorY::Array{T,2}, caretTensorX::Array{T,1}, caretTensorY::Array{T,1}, lattice::Classical2dFractalLattice{T})
 	@tensor begin
 		dumA[a,b,c,yp] := tensorS[a,b,c,d] * legTensorX[d,yp]
 		dummy1[a,xp,c,yp] := dumA[a,b,c,yp] * legTensorY[b,xp]
 		dummy2[x,xp,c,yp] := dummy1[a,xp,c,yp] * legTensorY[x,a]
 		newTensorT[x,xp,y,yp] := dummy2[x,xp,c,yp] * legTensorX[y,c]
 		newTensorPy[x,xp,yp] := dummy2[x,xp,c,yp] * caretTensorX[c]
 		dummy3[xp,c,yp] := dummy1[a,xp,c,yp] * caretTensorY[a]
 		newTensorPx[y,yp,xp] := dummy3[xp,c,yp] * legTensorX[y,c]
 		newTensorQ[xp,yp] := dummy3[xp,c,yp] * caretTensorX[c]

 		dumAt[a,b,c,yp] := tensorStilde[a,b,c,d] * legTensorX[d,yp]
 		dummy1t[a,xp,c,yp] := dumAt[a,b,c,yp] * legTensorY[b,xp]
 		dummy2t[x,xp,c,yp] := dummy1t[a,xp,c,yp] * legTensorY[x,a]
 		newTensorTtilde[x,xp,y,yp] := dummy2t[x,xp,c,yp] * legTensorX[y,c]
 	end
	setTensorT!(lattice, newTensorT)
	setTensorT!(lattice, newTensorTtilde; tilde = true)
	setTensorP!(lattice, 1, newTensorPx)
	setTensorP!(lattice, 2, newTensorPy)
	setTensorQ!(lattice, newTensorQ)
end

function constructHalf{T}(tensorT::Array{T,4}, tensorTtilde::Array{T,4}, lattice::Classical2dFractalLattice{T}, dimM::Int, XorY::Int; writefile = true, filename::AbstractString = "singularvalues.txt")
	tenmatM = getTensorM(tensorT)
	matrixMMd = tenmatM.matrix * transpose(tenmatM.matrix)
	Ul, lambdaVector, Uld = svd(matrixMMd)
	writefile && writeVector(lambdaVector, filename)
	trunUl = truncMatrixU(Ul,dimM)
	tensorUl = matU2tenU(trunUl)

	newTensorT = getNewTensorT(tensorT,tensorUl)
	newRotTensorT = rotateTensorT90(newTensorT)

	newTensorTtilde = getNewTensorT(tensorT,tensorUl,tensorTtilde)
	newRotTensorTtilde = rotateTensorT90(newTensorTtilde)

	tensorP = getTensorP(lattice, XorY)
	legTensor = getNewLegTensor(tensorP, tensorUl)
	caretTensor = getNewCaretTensor(getTensorQ(lattice), tensorUl, XorY)

	legextension = getLegextension(lattice)
	if legextension == 0
		legTensor = legTensor^getLegextension(lattice)
		caretTensor = ones(size(caretTensor))
	elseif legextension > 0
		extraleg = legTensor^(getLegextension(lattice) - 1)
		legTensor = legTensor * extraleg
		caretTensor = extraleg * caretTensor
	else
		error("legextension should be not negative integer")
	end

	return newRotTensorT, newRotTensorTtilde, legTensor, caretTensor
end

function getNewLegTensor{T}(tensorP::Array{T,3}, tensorU::Array{T,3})
	@tensor begin
		#####? maybe tensorU[x2,x1,i]
		dummy1[i,x2,s,x1p] := tensorU[x1,x2,i] * tensorP[x1,x1p,s]
		dummy2[i,x1p,x2p] := dummy1[i,x2,s,x1p] * tensorP[x2p,x2,s]
		legTensor[i,j] :=  dummy2[i,x1p,x2p] * tensorU[x1p,x2p,j]
	end
	return legTensor
end

function getNewCaretTensor{T}(tensorQ::Array{T,2}, tensorU::Array{T,3}, XorY::Int)
	if XorY == 1 # C[X]
		@tensor begin
			dummy[y1,y2] := tensorQ[s,y1] * tensorQ[s,y2]
			caretTensor[j] := dummy[y1,y2] * tensorU[y1,y2,j]
		end
	else # C[Y]
		@tensor begin
			dummy[x1,x2] := tensorQ[x1,s] * tensorQ[x2,s]
			caretTensor[j] := dummy[x1,x2] * tensorU[x1,x2,j]
		end
	end
	return caretTensor
end


text:
# # Binary Basics
# This tutorial is about operations about basis, it is mainly designed for developers, but is also useful to users.
using Yao, Yao.Intrinsics
#------------------
# ## Table of Contents
# * Storage of Computing Bases
# * Binray Operations
# * Number Readouts
# * Iterating over Bases
#---------------------
# ## Storage of Computing Bases
# We use an `Int` type to store spin configurations, e.g. `0b011101` (`29`) represents qubit configuration
# ```math
# \sigma_1=1, \sigma_2=0, \sigma_3=1, \sigma_4=1, \sigma_5=1, \sigma_6=0
# ```
# so we relate the configurations $\vec σ$ with integer $b$ by $b = \sum\limits_i 2^{i-1}σ_i$.
#----------------------------
# related APIs are
# * `integer(s) |> bitarray(nbit)`, transform integers to bistrings of type `BitArray`.
# * `bitstring |> packabits`, transform bitstrings to integers.
# * `integer |> baddrs`, get the locations of nonzero qubits.
@show 4 |> bitarray(5)
@show [4, 5, 6] |> bitarray(5)
@show [1, 1 , 0] |> packbits
@show [4, 5, 6] |> bitarray(5) |> packbits
@show baddrs(0b011);

# ## Binray Operations
takebit(0b11100, 2, 3)

# Masking is an important concept for binary operations, to generate a mask with specific position masked, e.g. we want to mask qubits `1, 3, 4`
mask = bmask(UInt8, 1,3,4)
@assert mask == 0b1101;

# with this mask, we can
@show testall(0b1011, mask) # true if all masked positions are 1
@show testany(0b1011, mask) # true if any masked positions is 1
@show testval(0b1011, mask, 0b1001)  # true if mask outed position matches `0b1001`
@show flip(0b1011, mask)  # flip masked positions
@show swapbits(0b1011, 0b1100)  # swap masked positions
@show setbit(0b1011, 0b1100);  # set masked positions 1

# For more interesting bitwise operations, see manual page [Yao.Intrinsics](@ref Intrinsics).
#----------------------------------
# ## Number Readouts
# In phase estimation and HHL algorithms, we sometimes need to readouts qubits as integer or float point numbers.
# We can read the register in different ways, like
# * bint, the integer itself
# * bint_r, the integer with bits small-big end reflected.
# * bfloat, the float point number 0.σ₁σ₂...σₙ.
# * bfloat_r, the float point number 0.σₙ...σ₂σ₁.
@show bint(0b010101)
@show bint_r(0b010101, nbit=6)
@show bfloat(0b010101)
@show bfloat_r(0b010101, nbit=6);

# Notice here functions with `_r` ending always require `nbit` as an additional input parameter to help reading, which is regarded as less natural way of expressing numbers.
#----------------------
# ## Iterating over Bases
# Counting from `0` is very natural way of iterating quantum registers, very pity for `Julia`
@show basis(4);

# `itercontrol` is a complicated API, but it plays an fundamental role in high performance quantum simulation of `Yao`.
# It is used for iterating over basis in controlled way, its interface looks like
@doc itercontrol
# Here, `poss` is a vector of controled positions, `vals` is a vector of values in controled positions.
#-------------
# ### example
# In a 4 qubit system, find out basis with 1st and 3rd qubits in state `0` and `1` respectively.
ic = itercontrol(4, [1,3], [0,1])
for i in ic
    println(i |> bitarray(4) .|> Int)
end
# Here, we have 1st and 3rd bits controlled,
# only 2 qubits are free, so the size of phase space here is `4`.


text:
using Distributions
using Parameters
using HamiltonianABC
using ContinuousTransformations
using Plots
using StatsBase
using StatPlots                 # for kernel density
import HamiltonianABC: logdensity, simulate!
using Base.Test
using PlotlyJS
plotlyjs()

############################################################################
## Simple simultaneous equations
############################################################################

## True model:
##  C_t = β * Y_t + us_t where  us_t ∼ N(0, τ)   (1)
##  Y_t = C_t + X_t                              (2)

## C_t is the consumption at time t
## X_t is the non-consumption at time t
## Y_t is the output at time t
## So output is used for consumption and non-consumption

## We have simultaneous equations as c_t depends on Y_t
## And Y_t depends on C_t as well

## Only X_t is exogenous in this model, Y_t and C_t are endogenous

"""
    simulate_simultaneous(β, X, us)

Take in the prameter (β), X and errors us, give back the endogenous variables of the system (Y and C).
"""
function simulate_simultaneous(β, X, us)
    N = length(us)
    C = Vector{Float64}
    Y = (X .+ us) / (1 - β)
    C = β * Y .+ us
    return (C, Y)
end

struct ToySimultaneousModel
    "observed consumption"
    Cs::Vector{Float64}
    "observed output"
    Ys::Vector{Float64}
    "non-comsumption"
    Xs::Vector{Float64}
    "distribution of Xs"
    dist_x
    "prior for β"
    prior_β
    "Normal(0,τ) draws for simulation, where τ is fixed"
    us::Vector{Float64}
    "distribution of us"
    dist_us
end

"""
    ToySimultaneousModel(Cs, Ys, prior_β, dist_x, dist_us, M)

Convenience constructor for ToySimultaneousModel.
Take in the observed data, the prior, and number of simulations (M).
"""
function ToySimultaneousModel(Cs, Ys, prior_β, dist_x, dist_us, M)
    ToySimultaneousModel(Cs, Ys, rand(dist_x, M), dist_x, prior_β, rand(dist_us, M), dist_us)
end

"""
    OLS(y, x)

Take in the dependant variable (y) and the regressor (x), give back the estimated coefficients (β) and the variance (σ_2).
"""
function OLS(y, x)
    β = x \ y
    err = (y - x * β)
    σ_2 = mean(abs2, err)
    return β, σ_2
end

"""
    simulate!(pp::ToySimultaneousModel)

Updates the shock and the non-consumption variable of the model.
"""
function simulate!(pp::ToySimultaneousModel)
    @unpack Xs, us, dist_x, dist_us = pp
    rand!(dist_x, Xs)
    rand!(dist_us, us)
end

## logdensity uses the following auxiliary model:
## C_t ∼  N(β_1 + β_2 * X_t, √σ_2)
function logdensity(pp::ToySimultaneousModel, β)
    @unpack Cs, Ys, Xs, prior_β, us = pp
    logprior = logpdf(prior_β, β[1])

    ## Generating the data
    C, Y = simulate_simultaneous(β[1], Xs, us)
    # OLs estimatation, regressing C on [1 X]
    est, σ_2 = OLS(C, [ones(length(us)) (Y - C)])
    β_1, β_2 = est

    log_likelihood = 0
    for i in 1:length(Cs)
        log_likelihood += logpdf(Normal(β_1 + β_2 * (Ys[i] - Cs[i]), √(σ_2)), Cs[i])
    end

    return(logprior + log_likelihood)
end

# Trial
β = 0.9
C, Y = simulate_simultaneous(β, rand(Normal(100, 3), 100), rand(Normal(0, 5), 100))
pp = ToySimultaneousModel(C, Y, Uniform(0, 1), Normal(100, 3), Normal(0, 5), 1000)
# variance is hand-tuned to get a ≈ 0.6
chain, a = mcmc(RWMH(diagm([4e-7])), pp, [0.1], 5000)
result = vcat(chain[2500:end]'...)

# Analysis with plotting
plt = plot(density(result), label = "posterior", title = "β")
plot!(plt, a -> pdf(pp.prior_β, a), linspace(0, 1, 100), label = "prior")
vline!(plt, [β], label = "true value")

@testset "simultaneous equation" begin
    @test mean(result) ≈ β rtol = 0.05
end


text:
#
#  Split - Apply - Combine operations
#

##############################################################################
##
## GroupedDataFrame...
##
##############################################################################

"""
The result of a `groupby` operation on an AbstractDataFrame; a
view into the AbstractDataFrame grouped by rows.

Not meant to be constructed directly, see `groupby`.
"""
type GroupedDataFrame
    parent::AbstractDataFrame
    cols::Vector         # columns used for sorting
    idx::Vector{Int}     # indexing vector when sorted by the given columns
    starts::Vector{Int}  # starts of groups
    ends::Vector{Int}    # ends of groups
end

#
# Split
#

function groupsort_indexer(x::AbstractVector, ngroups::Integer, null_last::Bool=false)
    # translated from Wes McKinney's groupsort_indexer in pandas (file: src/groupby.pyx).

    # count group sizes, location 0 for NULL
    n = length(x)
    # counts = x.pool
    counts = fill(0, ngroups + 1)
    for i = 1:n
        counts[x[i] + 1] += 1
    end

    # mark the start of each contiguous group of like-indexed data
    where = fill(1, ngroups + 1)
    if null_last
        for i = 3:ngroups+1
            where[i] = where[i - 1] + counts[i - 1]
        end
        where[1] = where[end] + counts[end]
    else
        for i = 2:ngroups+1
            where[i] = where[i - 1] + counts[i - 1]
        end
    end

    # this is our indexer
    result = fill(0, n)
    for i = 1:n
        label = x[i] + 1
        result[where[label]] = i
        where[label] += 1
    end
    result, where, counts
end

"""
A view of an AbstractDataFrame split into row groups

```julia
groupby(d::AbstractDataFrame, cols)
groupby(cols)
```

### Arguments

* `d` : an AbstractDataFrame to split (optional, see [Returns](#returns))
* `cols` : data frame columns to group by

### Returns

* `::GroupedDataFrame` : a grouped view into `d`
* `::Function`: a function `x -> groupby(x, cols)` (if `d` is not specified)

### Details

An iterator over a `GroupedDataFrame` returns a `SubDataFrame` view
for each grouping into `d`. A `GroupedDataFrame` also supports
indexing by groups and `map`.

See the following for additional split-apply-combine operations:

* `by` : split-apply-combine using functions
* `aggregate` : split-apply-combine; applies functions in the form of a cross product
* `combine` : combine (obviously)
* `colwise` : apply a function to each column in an AbstractDataFrame or GroupedDataFrame

Piping methods `|>` are also provided.

See the
[DataFramesMeta](https://github.com/JuliaStats/DataFramesMeta.jl)
package for more operations on GroupedDataFrames.

### Examples

```julia
df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
               b = repeat([2, 1], outer=[4]),
               c = randn(8))
gd = groupby(df, :a)
gd[1]
last(gd)
vcat([g[:b] for g in gd]...)
for g in gd
    println(g)
end
map(d -> mean(dropnull(d[:c])), gd)   # returns a GroupApplied object
combine(map(d -> mean(dropnull(d[:c])), gd))
df |> groupby(:a) |> [sum, length]
df |> groupby([:a, :b]) |> [sum, length]
```

"""
function groupby{T}(d::AbstractDataFrame, cols::Vector{T})
    ## a subset of Wes McKinney's algorithm here:
    ##     http://wesmckinney.com/blog/?p=489

    ncols = length(cols)
    # use CategoricalArray to get a set of integer references for each unique item
    nv = NullableCategoricalArray(d[cols[ncols]])
    # if there are NULLs, add 1 to the refs to avoid underflows in x later
    anynulls = (findfirst(nv.refs, 0) > 0 ? 1 : 0)
    # use UInt32 instead of the original array's integer size since the number of levels can be high
    x = similar(nv.refs, UInt32)
    for i = 1:nrow(d)
        if nv.refs[i] == 0
            x[i] = 1
        else
            x[i] = CategoricalArrays.order(nv.pool)[nv.refs[i]] + anynulls
        end
    end
    # also compute the number of groups, which is the product of the set lengths
    ngroups = length(levels(nv)) + anynulls
    # if there's more than 1 column, do roughly the same thing repeatedly
    for j = (ncols - 1):-1:1
        nv = NullableCategoricalArray(d[cols[j]])
        anynulls = (findfirst(nv.refs, 0) > 0 ? 1 : 0)
        for i = 1:nrow(d)
            if nv.refs[i] != 0
                x[i] += (CategoricalArrays.order(nv.pool)[nv.refs[i]] + anynulls - 1) * ngroups
            end
        end
        ngroups = ngroups * (length(levels(nv)) + anynulls)
        # TODO if ngroups is really big, shrink it
    end
    (idx, starts) = groupsort_indexer(x, ngroups)
    # Remove zero-length groupings
    starts = _uniqueofsorted(starts)
    ends = starts[2:end] - 1
    GroupedDataFrame(d, cols, idx, starts[1:end-1], ends)
end
groupby(d::AbstractDataFrame, cols) = groupby(d, [cols])

# add a function curry
groupby{T}(cols::Vector{T}) = x -> groupby(x, cols)
groupby(cols) = x -> groupby(x, cols)

Base.start(gd::GroupedDataFrame) = 1
Base.next(gd::GroupedDataFrame, state::Int) =
    (sub(gd.parent, gd.idx[gd.starts[state]:gd.ends[state]]),
     state + 1)
Base.done(gd::GroupedDataFrame, state::Int) = state > length(gd.starts)
Base.length(gd::GroupedDataFrame) = length(gd.starts)
Base.endof(gd::GroupedDataFrame) = length(gd.starts)
Base.first(gd::GroupedDataFrame) = gd[1]
Base.last(gd::GroupedDataFrame) = gd[end]

Base.getindex(gd::GroupedDataFrame, idx::Int) =
    sub(gd.parent, gd.idx[gd.starts[idx]:gd.ends[idx]])
Base.getindex(gd::GroupedDataFrame, I::AbstractArray{Bool}) =
    GroupedDataFrame(gd.parent, gd.cols, gd.idx, gd.starts[I], gd.ends[I])

Base.names(gd::GroupedDataFrame) = names(gd.parent)
_names(gd::GroupedDataFrame) = _names(gd.parent)

##############################################################################
##
## GroupApplied...
##    the result of a split-apply operation
##    TODOs:
##      - better name?
##      - ref
##      - keys, vals
##      - length
##      - start, next, done -- should this return (k,v) or just v?
##      - make it a real associative type? Is there a need to look up key columns?
##
##############################################################################

"""
The result of a `map` operation on a GroupedDataFrame; mainly for use
with `combine`

Not meant to be constructed directly, see `groupby` abnd
`combine`. Minimal support is provided for this type. `map` is
provided for a GroupApplied object.

"""
immutable GroupApplied{T<:AbstractDataFrame}
    gd::GroupedDataFrame
    vals::Vector{T}

    @compat function (::Type{GroupApplied})(gd::GroupedDataFrame, vals::Vector)
        length(gd) == length(vals) ||
            throw(DimensionMismatch("GroupApplied requires keys and vals be of equal length (got $(length(gd)) and $(length(vals)))."))
        new{eltype(vals)}(gd, vals)
    end
end

#
# Apply / map
#

# map() sweeps along groups
function Base.map(f::Function, gd::GroupedDataFrame)
    GroupApplied(gd, [wrap(f(df)) for df in gd])
end
function Base.map(f::Function, ga::GroupApplied)
    GroupApplied(ga.gd, [wrap(f(df)) for df in ga.vals])
end

wrap(df::AbstractDataFrame) = df
wrap(A::Matrix) = convert(DataFrame, A)
wrap(s::Any) = DataFrame(x1 = s)

"""
Combine a GroupApplied object (rudimentary)

```julia
combine(ga::GroupApplied)
```

### Arguments

* `ga` : a GroupApplied

### Returns

* `::DataFrame`

### Examples

```julia
df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
               b = repeat([2, 1], outer=[4]),
               c = randn(8))
combine(map(d -> mean(dropnull(d[:c])), gd))
```

"""
function combine(ga::GroupApplied)
    gd, vals = ga.gd, ga.vals
    valscat = vcat(vals)
    idx = Vector{Int}(size(valscat, 1))
    j = 0
    @inbounds for (start, val) in zip(gd.starts, vals)
        n = size(val, 1)
        idx[j + (1:n)] = gd.idx[start]
        j += n
    end
    hcat!(gd.parent[idx, gd.cols], valscat)
end

"""
Apply a function to each column in an AbstractDataFrame or
GroupedDataFrame

```julia
colwise(f::Function, d)
colwise(d)
```

### Arguments

* `f` : a function or vector of functions
* `d` : an AbstractDataFrame of GroupedDataFrame

If `d` is not provided, a curried version of groupby is given.

### Returns

* various, depending on the call

### Examples

```julia
df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
               b = repeat([2, 1], outer=[4]),
               c = randn(8))
colwise(sum, df)
colwise(sum, groupby(df, :a))
```

"""
colwise(f::Function, d::AbstractDataFrame) = Any[vcat(f(d[idx])) for idx in 1:size(d, 2)]
colwise(f::Function, gd::GroupedDataFrame) = map(colwise(f), gd)
colwise(f::Function) = x -> colwise(f, x)
colwise(f) = x -> colwise(f, x)
# apply several functions to each column in a DataFrame
colwise{T<:Function}(fns::Vector{T}, d::AbstractDataFrame) =
    reshape(Any[vcat(f(d[idx])) for f in fns, idx in 1:size(d, 2)],
            length(fns)*size(d, 2))
colwise{T<:Function}(fns::Vector{T}, gd::GroupedDataFrame) = map(colwise(fns), gd)
colwise{T<:Function}(fns::Vector{T}) = x -> colwise(fns, x)

"""
Split-apply-combine in one step; apply `f` to each grouping in `d`
based on columns `col`

```julia
by(d::AbstractDataFrame, cols, f::Function)
by(f::Function, d::AbstractDataFrame, cols)
```

### Arguments

* `d` : an AbstractDataFrame
* `cols` : a column indicator (Symbol, Int, Vector{Symbol}, etc.)
* `f` : a function to be applied to groups; expects each argument to
  be an AbstractDataFrame

`f` can return a value, a vector, or a DataFrame. For a value or
vector, these are merged into a column along with the `cols` keys. For
a DataFrame, `cols` are combined along columns with the resulting
DataFrame. Returning a DataFrame is the clearest because it allows
column labeling.

A method is defined with `f` as the first argument, so do-block
notation can be used.

`by(d, cols, f)` is equivalent to `combine(map(f, groupby(d, cols)))`.

### Returns

* `::DataFrame`

### Examples

```julia
df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
               b = repeat([2, 1], outer=[4]),
               c = randn(8))
by(df, :a, d -> sum(d[:c]))
by(df, :a, d -> 2 * dropnull(d[:c]))
by(df, :a, d -> DataFrame(c_sum = sum(d[:c]), c_mean = mean(dropnull(d[:c]))))
by(df, :a, d -> DataFrame(c = d[:c], c_mean = mean(dropnull(d[:c]))))
by(df, [:a, :b]) do d
    DataFrame(m = mean(dropnull(d[:c])), v = var(dropnull(d[:c])))
end
```

"""
by(d::AbstractDataFrame, cols, f::Function) = combine(map(f, groupby(d, cols)))
by(f::Function, d::AbstractDataFrame, cols) = by(d, cols, f)

#
# Aggregate convenience functions
#

# Applies a set of functions over a DataFrame, in the from of a cross-product
"""
Split-apply-combine that applies a set of functions over columns of an
AbstractDataFrame or GroupedDataFrame

```julia
aggregate(d::AbstractDataFrame, cols, fs)
aggregate(gd::GroupedDataFrame, fs)
```

### Arguments

* `d` : an AbstractDataFrame
* `gd` : a GroupedDataFrame
* `cols` : a column indicator (Symbol, Int, Vector{Symbol}, etc.)
* `fs` : a function or vector of functions to be applied to vectors
  within groups; expects each argument to be a column vector

Each `fs` should return a value or vector. All returns must be the
same length.

### Returns

* `::DataFrame`

### Examples

```julia
df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
               b = repeat([2, 1], outer=[4]),
               c = randn(8))
aggregate(df, :a, sum)
aggregate(df, :a, [sum, x->mean(dropnull(x))])
aggregate(groupby(df, :a), [sum, x->mean(dropnull(x))])
df |> groupby(:a) |> [sum, x->mean(dropnull(x))]   # equivalent
```

"""
aggregate(d::AbstractDataFrame, fs::Function) = aggregate(d, [fs])
function aggregate{T<:Function}(d::AbstractDataFrame, fs::Vector{T})
    headers = _makeheaders(fs, _names(d))
    _aggregate(d, fs, headers)
end

# Applies aggregate to non-key cols of each SubDataFrame of a GroupedDataFrame
aggregate(gd::GroupedDataFrame, f::Function) = aggregate(gd, [f])
function aggregate{T<:Function}(gd::GroupedDataFrame, fs::Vector{T})
    headers = _makeheaders(fs, _setdiff(_names(gd), gd.cols))
    combine(map(x -> _aggregate(without(x, gd.cols), fs, headers), gd))
end
(|>)(gd::GroupedDataFrame, fs::Function) = aggregate(gd, fs)
(|>){T<:Function}(gd::GroupedDataFrame, fs::Vector{T}) = aggregate(gd, fs)

# Groups DataFrame by cols before applying aggregate
function aggregate{S <: ColumnIndex, T <:Function}(d::AbstractDataFrame,
                                     cols::@compat(Union{S, AbstractVector{S}}),
                                     fs::@compat(Union{T, Vector{T}}))
    aggregate(groupby(d, cols), fs)
end

function _makeheaders{T<:Function}(fs::Vector{T}, cn::Vector{Symbol})
    fnames = _fnames(fs) # see other/utils.jl
    reshape([Symbol(colname,'_',fname) for fname in fnames, colname in cn],
            length(fnames)*length(cn))
end

function _aggregate{T<:Function}(d::AbstractDataFrame, fs::Vector{T}, headers::Vector{Symbol})
    DataFrame(colwise(fs, d), headers)
end


text:
#=

  Euler #37 in Julia.

  Problem 37:
  """
  The number 3797 has an interesting property. Being prime itself, it is possible to
  continuously remove digits from left to right, and remain prime at each stage:
  3797, 797, 97, and 7. Similarly we can work from right to left: 3797, 379, 37, and 3.

  Find the sum of the only eleven primes that are both truncatable from left to right
  and right to left.

  NOTE: 2, 3, 5, and 7 are not considered to be truncatable primes.
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

# check is n is a truncated prime
function truncated_prime(n)
    len = nlen(n)
    for i in 1:len-1
        ii = 10^i
        if !isPrime(n % ii) || !isPrime(div(n,ii))
            return false
        end
    end
    return true
end

# 0.02413909s
function euler37a()
    p = 11
    sum = 0
    c = 0
    while c < 11
        if truncated_prime(p) && isPrime(p)
            c += 1
            sum += p
        end
        p += 2
    end
    return sum
end

run_euler(euler37a)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
using BinaryProvider
# This is where all binaries will get installed
const prefix = Prefix(!isempty(ARGS) ? ARGS[1] : joinpath(@__DIR__,"usr"))
# Instantiate products here.  Examples:
libstemmer = LibraryProduct(prefix, "libstemmer", :libstemmer)
# foo_executable = ExecutableProduct(prefix, "fooifier")
# libfoo_pc = FileProduct(joinpath(libdir(prefix), "pkgconfig", "libfoo.pc"))
# Assign products to `products`:
products = [libstemmer]
# Download binaries from hosted location
bin_prefix = "https://github.com/aviks/SnowballBuilder/releases/download/v0.1.0"
# Listing of files generated by BinaryBuilder:
download_info = Dict(
    BinaryProvider.Linux(:aarch64, :glibc) => ("$bin_prefix/LibStemmer.aarch64-linux-gnu.tar.gz", "9fc5727ee09ea42a9ea0234f3cb87973da06801c51d61a03340fc09d52bb9a9f"),
    BinaryProvider.Linux(:armv7l, :glibc) => ("$bin_prefix/LibStemmer.arm-linux-gnueabihf.tar.gz", "6741ccb3e07c1f0a437a681860bea4280cd26e45360456a594202b4af1e2525c"),
    BinaryProvider.Linux(:i686, :glibc) => ("$bin_prefix/LibStemmer.i686-linux-gnu.tar.gz", "635dfd5016439798d4ce8b6b816329fa27f61aedc29d146e6c27bed9686f9fa7"),
    BinaryProvider.Windows(:i686) => ("$bin_prefix/LibStemmer.i686-w64-mingw32.tar.gz", "c3818e8971124436afc588c3dbf93d42d5ebc23f6abd81e7c3fd3bc9a2531e39"),
    BinaryProvider.Linux(:powerpc64le, :glibc) => ("$bin_prefix/LibStemmer.powerpc64le-linux-gnu.tar.gz", "c44a22f4dad364b954ea0906f4de65ed3e32433facb7621c20be4c52f71ff510"),
    BinaryProvider.MacOS() => ("$bin_prefix/LibStemmer.x86_64-apple-darwin14.tar.gz", "991f333ee1c4cc1e87a07cd76028de1604181e0b5b20f71d0d460aab1f20424d"),
    BinaryProvider.Linux(:x86_64, :glibc) => ("$bin_prefix/LibStemmer.x86_64-linux-gnu.tar.gz", "9daf7ca5f2df44f8208c719dae99e6d52f42285e96a64140a0814b188704dcdc"),
    BinaryProvider.Windows(:x86_64) => ("$bin_prefix/LibStemmer.x86_64-w64-mingw32.tar.gz", "3fa2aeb67cfd8cd62ad67543d2cdbdd362e0bcaf804e33a1d90e0e4c983f7d6a"),
)
if platform_key() in keys(download_info)
    # First, check to see if we're all satisfied
    if any(!satisfied(p; verbose=true) for p in products)
        # Download and install binaries
        url, tarball_hash = download_info[platform_key()]
        install(url, tarball_hash; prefix=prefix, force=true, verbose=true)
    end
    # Finally, write out a deps.jl file that will contain mappings for each
    # named product here: (there will be a "libfoo" variable and a "fooifier"
    # variable, etc...)
    write_deps_file(joinpath(@__DIR__, "deps.jl"), products)
else
    error("Your platform $(Sys.MACHINE) is not supported by this package!")
end


text:
# functions that return the xml manipulation expressions for the each argument

################################################################
"""
Get a argument creator expression
"""
function get_arg_xmlcreator(argcustomcreator, has_arg_xmlchecker::Bool, argtype, esc_argvar, argname, argliteraltype, esc_argfunction, argsym, argvarcall)

    if !has_arg_xmlchecker
        arg_creator = quote
            $(esc(argcustomcreator))
            addnode!(aml, $argname, $esc_argvar, $argliteraltype)
        end
    else
        arg_creator=quote
            $(esc(argcustomcreator))
            if isnothing($esc_argvar) || ($esc_argfunction)($esc_argvar)
                addnode!(aml, $argname, $esc_argvar, $argliteraltype)
            else
                error("$($argname) doesn't meet criteria function")
            end
        end
    end
    return arg_creator
end
################################################################
"""
Get a argument extractor expression
"""
function get_arg_xmlextractor(argcustomextractor, has_arg_xmlchecker::Bool, esc_argtype, esc_argvar, argname, argliteraltype, esc_argfunction, argsym, argvarcall)

    if !has_arg_xmlchecker
        arg_extractor=quote
            $(esc(argcustomextractor))
            $esc_argvar = findcontent($esc_argtype, $argname, aml, $argliteraltype)
        end
    else
        arg_extractor=quote
            $(esc(argcustomextractor))

            $esc_argvar = findcontent($esc_argtype, $argname, aml, $argliteraltype)

            if !isnothing($esc_argvar) && !(($esc_argfunction)($esc_argvar))
                error("$($argname) doesn't meet criteria function")
            end
        end
    end
    return arg_extractor
end
################################################################
"""
Get a argument updater expression
"""
function get_arg_xmlupdater(argcustomupdater, has_arg_xmlchecker::Bool, esc_argtype, esc_argvar, argname, argliteraltype, esc_argfunction, argsym, argvarcall)

    if !has_arg_xmlchecker
        arg_updater = quote
            $(esc(argcustomupdater))

            if name == $argsym
                updatecontent!(value, $argname, str.aml, $argliteraltype)
            end
        end
    else
        arg_updater = quote
            $(esc(argcustomupdater))

            if name == $argsym
                if isnothing($(argvarcall)) || ($esc_argfunction)($(argvarcall))
                    updatecontent!(value, $argname, str.aml, $argliteraltype)
                else
                    error("$($argname) doesn't meet criteria function")
                end
            end
        end
    end
    return arg_updater
end

################################################################
# """
# Parses argument type for sentence creation. Only three output are considered:
# - Vector
# - Dict
# - Other (Any)
#
# Not used anymore, but can be used to determine the correct xml utils sets of functions
# """
# function parse_argtype(argtype)
#     # Vector
#     if (isa(argtype, Expr) && argtype.args[1] == :Vector) ||
#        (!isa(argtype, Union{Symbol,Expr}) && argtype <: AbstractVector)
#
#        arg_parsedtype = AbstractVector
#
#     # Dict
#     elseif (isa(argtype, Expr) && argtype.args[1] == :Dict) ||
#            (!isa(argtype, Union{Symbol,Expr}) && argtype <: AbstractDict)
#
#         arg_parsedtype = AbstractDict
#
#     # # Non Vector
#     # elseif isa(argtype, Symbol) ||
#     #        (isa(argtype, Expr) && argtype.args[1] == :Union) ||
#     #        (isa(argtype, Expr) && argtype.args[1] == :UN) ||
#     #        !(esc_argtype <: AbstractVector)
#     else
#         arg_parsedtype = Any
#     end
#
#     return arg_parsedtype
# end

max_stars_repo_name: aminya/AML

max_forks_repo_name: aminya/AML


text:
@testset "Polysegment Binarization" begin
    # Verify that the output type matches the input type.
    img = testimage("lena")
    for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64})
        img₁ = T.(img)
        img₂ = binarize(Polysegment(), img₁)
        @test eltype(img₁) == eltype(img₂)
    end
end