text:
"""
    RCutOut(m; rng = GLOBAL_RNG)
    RCutOut(m...; rng = GLOBAL_RNG)

This operation applies a mask of zeros to a random location within image. `m`
is a tuple or integers defining the mask size.

# Example
```julia
julia> augment(RCutOut((100, 100))) |> plot

julia> augment(RCutOut(100, 100)) |> plot
```
"""
struct RCutOut{N, I <: Tuple, R <: AbstractRNG} <: ImageOperation
    msze::I
    rng::R
end

RCutOut{N}(m::T, rng::R = GLOBAL_RNG) where {N, T <: NTuple{N, <:Integer}, R} =
    RCutOut{N, typeof(m), R}(m, rng)

RCutOut(m; rng = GLOBAL_RNG) = RCutOut{length(m)}(m, rng)

RCutOut(m::Integer...; rng = GLOBAL_RNG) = RCutOut{length(m)}(m, rng)

@inline supports_eager(::Type{<:RCutOut}) = false

@inline supports_affineview(::Type{<:RCutOut}) = false

@inline supports_view(::Type{<:RCutOut}) = false

@inline supports_stepview(::Type{<:RCutOut}) = true

applystepview(op::RCutOut, img::AbstractArray, pr) = applyeager(op, img, pr)

function applyeager(op::RCutOut{N}, img::AbstractArray, pr) where N
    st = rand.(Ref(op.rng), UnitRange.(1, size(img)[1:N] .- op.msze))
    ct = UnitRange.(st, st .+ op.msze .- 1)
    idx = repeat(Any[Colon()], ndims(img))
    for (i, c) in enumerate(ct)
        if !isempty(c)
            idx[i] = c
        end
    end
    nim = copy(img)
    z = img |> eltype |> zero
    nim[idx...] .= fill(z, op.msze)
    nim
end


text:
#This file is simply a small script to load the Utils files/types in the correct order

#import needed libraries
import Base: norm

#include MathUtils.jl first
include("MathUtils.jl")

#then Quaterion Library
include("QuatLib.jl")
include("State.jl")

ext: jl


text:
using ExpressCommands
using Documenter

DocMeta.setdocmeta!(ExpressCommands, :DocTestSetup, :(using ExpressCommands); recursive=true)

makedocs(;
    modules=[ExpressCommands],
    authors="Reno <singularitti@outlook.com>",
    repo="https://github.com/MineralsCloud/ExpressCommands.jl/blob/{commit}{path}#{line}",
    sitename="ExpressCommands.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://MineralsCloud.github.io/ExpressCommands.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Commands" => "commands.md",
    ],
)

deploydocs(;
    repo="github.com/MineralsCloud/ExpressCommands.jl",
)


text:
module GraphDatasets

import Base: getindex, length, eltype, iterate, show, firstindex, lastindex

import Graphs: loadgraphs, SimpleGraph

import SimpleValueGraphs:
    nv, has_edge, is_directed,
    get_vertexval, get_edgeval, get_graphval,
    outedgevals, outneighbors,
    ValGraph

export
    list_datasets,

    loadgraphs,
    loadreadme,

    ValGraphCollection,
    ValGraphCollectionView,
    ng,

    TUDatasets,

    # reexport from Base
    getindex,
    length,
    eltype,
    iterate,
    show,
    firstindex,
    lastindex,

    # reexport from Graphs.jl & SimpleValueGraphs
    nv,
    has_edge,

    get_vertexval,
    get_edgeval,
    get_graphval,

    outedgevals,
    outneighbors,

    SimpleGraph,
    ValGraph

include("utils.jl")
include("graphdataset.jl")
include("valgraphcollection.jl")
include("TUDatasets.jl")

end


text:
# this is a wrapper module for the benchmarks
module ARCH2018_RE
    nothing
end # module


text:
"""
Convert an array of tuples to a 2D array.
"""
# could restrict to number types, but allowing all types for now
#function array_from_tuples(T::Array{Tuple{Vararg{T,N}},1} where {T<:Number,N})
function array_from_tuples(arrtup::Array{Tuple{Vararg{T,N}},1} where {T,N})
    m = length(arrtup)
    n = length(arrtup[1])
    typename = typeof(arrtup[1][1])
    A = Array{typename}(undef, (m,n))
    for i in 1:m
        for j in 1:n
            A[i,j] = arrtup[i][j]
        end
    end
    return A
end
##


text:
@testset "resource.jl" begin
    shared = nothing
    nested = nothing
    sharedList = []
    nestedList = []
    @process control() begin
        @with_resource shared begin
            @with_resource nested begin
                push!(nestedList, nested.available.value)
                @schedule now nestedConsumer()
                work(2)
            end
            push!(nestedList, nested.available.value)
            push!(sharedList, shared.available.value)
            @schedule now Consumer()
            work(2)
        end
        # control has released its resource, went from 0 available to 1
        push!(sharedList, shared.available.value)

        request(shared)
        push!(sharedList, shared.available.value)
        release(shared)
        push!(sharedList, shared.available.value)

        # extra release
        release(shared)

    end

    @process Consumer() begin
        @with_resource shared begin
            push!(sharedList, shared.available.value)
            work(4)
        end
        # Consumer has released its resource, went from 1 available to 2
        push!(sharedList, shared.available.value)
    end

    @process nestedConsumer() begin
        @with_resource nested begin
            push!(nestedList, nested.available.value)
            work(4)
        end
        push!(nestedList, nested.available.value)
    end

    @simulation begin
        shared = Resource(2, "")
        nested = Resource(2, "")
        push!(nestedList, nested.available.value)
        push!(sharedList, shared.available.value)
        @schedule now control()
        start_simulation()
    end

    wantShared = [2,1,0,1,0,1,2]
    wantNested = [2,1,0,1,2]
    @test sharedList == wantShared
    @test nestedList == wantNested

end


text:
export mean_shift

using mlpack._Internal.io

import mlpack_jll
const mean_shiftLibrary = mlpack_jll.libmlpack_julia_mean_shift

# Call the C binding of the mlpack mean_shift binding.
function mean_shift_mlpackMain()
  success = ccall((:mean_shift, mean_shiftLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module mean_shift_internal
  import ..mean_shiftLibrary

end # module

"""
    mean_shift(input; [force_convergence, in_place, labels_only, max_iterations, radius, verbose])

This program performs mean shift clustering on the given dataset, storing the
learned cluster assignments either as a column of labels in the input dataset or
separately.

The input dataset should be specified with the `input` parameter, and the radius
used for search can be specified with the `radius` parameter.  The maximum
number of iterations before algorithm termination is controlled with the
`max_iterations` parameter.

The output labels may be saved with the `output` output parameter and the
centroids of each cluster may be saved with the `centroid` output parameter.

For example, to run mean shift clustering on the dataset `data` and store the
centroids to `centroids`, the following command may be used:

```julia
julia> using CSV
julia> data = CSV.read("data.csv")
julia> centroids, _ = mean_shift(data)
```

# Arguments

 - `input::Array{Float64, 2}`: Input dataset to perform clustering on.
 - `force_convergence::Bool`: If specified, the mean shift algorithm will
      continue running regardless of max_iterations until the clusters converge.
       Default value `false`.

 - `in_place::Bool`: If specified, a column containing the learned cluster
      assignments will be added to the input dataset file.  In this case,
      --output_file is overridden.  (Do not use with Python.)  Default value
      `false`.

 - `labels_only::Bool`: If specified, only the output labels will be
      written to the file specified by --output_file.  Default value `false`.

 - `max_iterations::Int`: Maximum number of iterations before mean shift
      terminates.  Default value `1000`.

 - `radius::Float64`: If the distance between two centroids is less than
      the given radius, one will be removed.  A radius of 0 or less means an
      estimate will be calculated and used for the radius.  Default value `0`.

 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `centroid::Array{Float64, 2}`: If specified, the centroids of each
      cluster will be written to the given matrix.
 - `output::Array{Float64, 2}`: Matrix to write output labels or labeled
      data to.

"""
function mean_shift(input;
                    force_convergence::Union{Bool, Missing} = missing,
                    in_place::Union{Bool, Missing} = missing,
                    labels_only::Union{Bool, Missing} = missing,
                    max_iterations::Union{Int, Missing} = missing,
                    radius::Union{Float64, Missing} = missing,
                    verbose::Union{Bool, Missing} = missing,
                    points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, mean_shiftLibrary), Nothing, ());

  IORestoreSettings("Mean Shift Clustering")

  # Process each input argument before calling mlpackMain().
  IOSetParamMat("input", input, points_are_rows)
  if !ismissing(force_convergence)
    IOSetParam("force_convergence", convert(Bool, force_convergence))
  end
  if !ismissing(in_place)
    IOSetParam("in_place", convert(Bool, in_place))
  end
  if !ismissing(labels_only)
    IOSetParam("labels_only", convert(Bool, labels_only))
  end
  if !ismissing(max_iterations)
    IOSetParam("max_iterations", convert(Int, max_iterations))
  end
  if !ismissing(radius)
    IOSetParam("radius", convert(Float64, radius))
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("centroid")
  IOSetPassed("output")
  # Call the program.
  mean_shift_mlpackMain()

  return IOGetParamMat("centroid", points_are_rows),
         IOGetParamMat("output", points_are_rows)
end


text:
# Function that updates partial force and energy vectors, used for parallel computation

function compute_uf_partials_QUAD!(it :: Int64, f :: Forces,
                                   i :: Int64, j :: Int64,
                                   xj :: Float64, yj :: Float64, r :: Float64,
                                   input :: MDInput)

   f.upartial[it] = f.upartial[it] + input.eps*(input.sig-r)^2

   drdxi = -xj/r
   drdyi = -yj/r

   fx = 2*input.eps*(input.sig-r)*drdxi
   fy = 2*input.eps*(input.sig-r)*drdyi

   f.fpartial[it][i,1] = f.fpartial[it][i,1] + fx
   f.fpartial[it][i,2] = f.fpartial[it][i,2] + fy
   f.fpartial[it][j,1] = f.fpartial[it][j,1] - fx
   f.fpartial[it][j,2] = f.fpartial[it][j,2] - fy

end


text:
function assert_input_parameters(
	agentcount::Int64,
	nettopology::String,
	networkprops::Dict,
	stubborncount::Int64
)
	if stubborncount > agentcount
		stubborncount = trunc(Int64, agentcount / 10)
		print_warning(
			"agentcount must be greater than stubborncount",
			"stubborncount", stubborncount
		)
	end
	if nettopology == "grid"
		handle_faulty_grid_params!(networkprops, agentcount)
	elseif nettopology == "watts_strogatz"
		handle_faulty_watts_strogatz_params!(networkprops, agentcount)
	elseif nettopology == "barabasi_albert"
		handle_faulty_barabasi_albert_params!(networkprops, agentcount)
	end
	return agentcount, nettopology, networkprops, stubborncount
end

function handle_faulty_grid_params!(networkprops, agentcount)
	if !("grid_height" in keys(networkprops))
		grid_height = 1
		print_warning(
			"parameter 'grid_height' is missing from networkprops",
			"grid_height", grid_height
		)
		networkprops["grid_height"] = grid_height
	else
		if (agentcount % networkprops["grid_height"]) != 0
			grid_height = 1
			print_warning(
				"agentcount must be divisible by grid_height",
				"grid_height", grid_height
			)
			networkprops["grid_height"] = grid_height
		end
	end
	return networkprops
end

function handle_faulty_watts_strogatz_params!(networkprops, agentcount)
	if !("k" in keys(networkprops))
		k = 2
		print_warning(
			"parameter 'k' is missing from networkprops",
			"k", k
		)
		networkprops["k"] = 2
	end
	if !("beta" in keys(networkprops))
		beta = 0.5
		print_warning(
			"parameter 'beta' is missing from networkprops",
			"beta", beta
		)
		networkprops["beta"] = beta
	end
	if "k" in keys(networkprops)
		if networkprops["k"] <= 0
			k = 2
			print_warning(
				"parameters 'k' must be greater than 0",
				"k", k
			)
			networkprops["k"] = k
		end
		if networkprops["k"] > (agentcount / 10)
			k = trunc(Int64, (agentcount / 10))
			print_warning(
				"Watts-Strogatz model requires n >> k",
				"k", k
			)
			networkprops["k"] = k
		end
	end
	if "beta" in keys(networkprops)
		if (networkprops["beta"] < 0) | (networkprops["beta"] > 1)
			beta = 0.5
			print_warning(
				"parameter 'beta' must be between 0 and 1",
				"beta", beta
			)
			networkprops["beta"] = 0.5
		end
	end
	return networkprops
end

function handle_faulty_barabasi_albert_params!(networkprops, agentcount)
	if !("m0" in keys(networkprops))
		m0 = trunc(Int64, (agentcount / 10))
		print_warning(
			"parameter 'm0' is missing from networkprops",
			"m0", m0
		)
		networkprops["m0"] = m0
	else
		if networkprops["m0"] > agentcount
			m0 = trunc(Int64, (agentcount / 10))
			print_warning(
				"Barabási-Albert model requires m0 < agentcount",
				"m0", m0
			)
			networkprops["m0"] = m0
		end
	end
	return networkprops
end

function print_warning(message::String, varname::String, default_value::Any)
	println("WARNING: " * message)
	println("         defaulted to " * varname * " = $default_value")
end

max_stars_repo_name: JohannesNakayama/axelrodnet

max_issues_repo_name: JohannesNakayama/axelrodnet

max_forks_repo_name: JohannesNakayama/axelrodnet


text:
using ProfileView
using MCTS
using POMDPToolbox
using POMDPs
using Multilane
using POMCPOW
using D3Trees

tii = true
wup = WeightUpdateParams(smoothing=0.0, wrong_lane_factor=0.5)

# solver = POMCPOWSolver(tree_queries=10_000_000,
#                                criterion=MaxUCB(10.0),
#                                max_depth=20,
#                                max_time=1.0,
#                                enable_action_pw=false,
#                                k_observation=4.0,
#                                alpha_observation=1/8,
#                                estimate_value=FORollout(SimpleSolver()),
#                                check_repeat_obs=false,
#                                node_sr_belief_updater=AggressivenessPOWFilter(wup),
#                                tree_in_info=tii
#                               )

dpws = DPWSolver(depth=20,
                 n_iterations=100_000,
                 max_time=1.0,
                 exploration_constant=2.0,
                 k_state=4.0,
                 alpha_state=1/8,
                 enable_action_pw=false,
                 check_repeat_state=false,
                 estimate_value=RolloutEstimator(SimpleSolver()),
                 tree_in_info=tii
                )
# solver = MLMPCSolver(dpws)
solver = QBSolver(dpws)

behaviors = standard_uniform(correlation=0.75)
pp = PhysicalParam(4, lane_length=100.0)
dmodel = NoCrashIDMMOBILModel(10, pp,
                              behaviors=behaviors,
                              p_appear=1.0,
                              lane_terminate=true,
                              max_dist=500.0
                             )
rmodel = SuccessReward(lambda=2,
                       target_lane=4,
                       brake_penalty_thresh=4.0
                      )
pomdp = NoCrashPOMDP{typeof(rmodel)}(dmodel, rmodel, 0.95, false)
mdp = NoCrashMDP{typeof(rmodel)}(dmodel, rmodel, 0.95, false)
is = initial_state(pomdp, Base.GLOBAL_RNG)
ips = MLPhysicalState(is)

agg_up = AggressivenessUpdater(pomdp, 500, 0.1, 0.1, wup, MersenneTwister(50000))
ib = initialize_belief(agg_up, ips)

planner = solve(solver, mdp)

@show collect(actions(pomdp, ips))
@show collect(actions(pomdp, ib))

action(planner, ib)
@time action(planner, ib)
Profile.clear()
@profile action(planner, ib)

a, i = action_info(planner, ib)
@show i[:tree_queries]

ProfileView.view()
inchrome(D3Tree(i[:tree]))


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

@doc raw"""Subject contains a reference to the object or user identities a role binding applies to.  This can either hold a direct API object reference, or a value for non-objects such as user and group names.

    IoK8sApiRbacV1beta1Subject(;
        apiGroup=nothing,
        kind=nothing,
        name=nothing,
        namespace=nothing,
    )

    - apiGroup::String : APIGroup holds the API group of the referenced subject. Defaults to \"\" for ServiceAccount subjects. Defaults to \"rbac.authorization.k8s.io\" for User and Group subjects.
    - kind::String : Kind of object being referenced. Values defined by this API group are \"User\", \"Group\", and \"ServiceAccount\". If the Authorizer does not recognized the kind value, the Authorizer should report an error.
    - name::String : Name of the object being referenced.
    - namespace::String : Namespace of the referenced object.  If the object kind is non-namespace, such as \"User\" or \"Group\", and this value is not empty the Authorizer should report an error.
"""
mutable struct IoK8sApiRbacV1beta1Subject <: SwaggerModel
    apiGroup::Any # spec type: Union{ Nothing, String } # spec name: apiGroup
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    namespace::Any # spec type: Union{ Nothing, String } # spec name: namespace

    function IoK8sApiRbacV1beta1Subject(;apiGroup=nothing, kind=nothing, name=nothing, namespace=nothing)
        o = new()
        validate_property(IoK8sApiRbacV1beta1Subject, Symbol("apiGroup"), apiGroup)
        setfield!(o, Symbol("apiGroup"), apiGroup)
        validate_property(IoK8sApiRbacV1beta1Subject, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiRbacV1beta1Subject, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiRbacV1beta1Subject, Symbol("namespace"), namespace)
        setfield!(o, Symbol("namespace"), namespace)
        o
    end
end # type IoK8sApiRbacV1beta1Subject

const _property_map_IoK8sApiRbacV1beta1Subject = Dict{Symbol,Symbol}(Symbol("apiGroup")=>Symbol("apiGroup"), Symbol("kind")=>Symbol("kind"), Symbol("name")=>Symbol("name"), Symbol("namespace")=>Symbol("namespace"))
const _property_types_IoK8sApiRbacV1beta1Subject = Dict{Symbol,String}(Symbol("apiGroup")=>"String", Symbol("kind")=>"String", Symbol("name")=>"String", Symbol("namespace")=>"String")
Base.propertynames(::Type{ IoK8sApiRbacV1beta1Subject }) = collect(keys(_property_map_IoK8sApiRbacV1beta1Subject))
Swagger.property_type(::Type{ IoK8sApiRbacV1beta1Subject }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiRbacV1beta1Subject[name]))}
Swagger.field_name(::Type{ IoK8sApiRbacV1beta1Subject }, property_name::Symbol) =  _property_map_IoK8sApiRbacV1beta1Subject[property_name]

function check_required(o::IoK8sApiRbacV1beta1Subject)
    (getproperty(o, Symbol("kind")) === nothing) && (return false)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiRbacV1beta1Subject }, name::Symbol, val)
end


text:
"""
This is the parent type for all surfaces. Surfaces
are purely a geometrical abstraction.
"""
abstract type AbstractSurface{T<:Real} end

abstract type AbstractRotationalSymmetricSurface{T} <: AbstractSurface{T} end

function sag(x, y, s::AbstractRotationalSymmetricSurface)
    radius = sqrt(x^2 + y^2)
    return (sag(radius, s))
end

function transfer_and_refract(ray, n1, t, s::AbstractSurface, n2, λ)
    transfered = transfer_to_intersection!(ray, t, s)
    if iserror(transfered)
        return transfered
    else
        refracted = refract!(transfered, n1, s, n2, λ)
        return refracted
    end
end

struct OpticalSurface{T<:Real}
    surface::AbstractSurface{T}
    aperture::Union{ClearDiameter{T},Nothing}
    n::AbstractMedium{T}
    t::T
    id::Union{Symbol,Nothing}
end


text:
module HeatEqCalc

function laplacian4!(v, u, dx)
    v[begin] = (u[end] + u[begin+1] - 2u[begin])/dx^2
    @. @views v[begin+1:end-1] = (u[begin:end-2] + u[begin+2:end] - 2u[begin+1:end-1])/dx^2
    v[end] = (u[end-1] + u[begin] - 2u[1end])/dx^2
    return
end

function heateq4(u0, dx, tmax, N=200)
    t = 0:dx:tmax
    dt = step(t)
    u = similar(u0, length(u0), length(t)+1)
    u[:, 1] = u0
    v = similar(u0)
    N = 100
    for i in 2:length(t)+1
        @. @views u[:, i] = u[:, i-1]
        for _ in 1:N
            @views laplacian4!(v, u[:, i], dx)
            @. @views u[:, i] += v*dt/N
        end
    end
    t, u
end

using StaticArrays

function calc_sol(
        n = 200,
        x = range(-π, π; length=n+1)[1:end-1],
        m = 200
    )
    dx = step(x)
    noise1, noise2 = 0.3randn(m), 0.3randn(m)
    f(x) = SVector{m}((1 .+ noise1)*sin(x) .- (1 .+ noise2)*cos(2x))
    u0 = f.(x)
    tmax = 1.0
    t, u = heateq4(u0, dx, tmax)
    (; t, x, u, m)
end

using DelimitedFiles

function save_sol(sol)
    t, x, u, m = sol
    U = vcat((getindex.(u, k) for k in 1:m)...)
    writedlm("sol_m.txt", m)
    writedlm("sol_t.txt", t)
    writedlm("sol_x.txt", x)
    writedlm("sol_U.txt", U)
end

function load_sol(name)
    m = parse(Int, read(name * "_m.txt", String))
    t = vec(readdlm(name * "_t.txt"))
    x = vec(readdlm(name * "_x.txt"))
    U = readdlm(name * "_U.txt")
    u = [SVector((U[k*length(x) + i, j] for k in 0:m-1)...) for i in 1:length(x), j in 1:length(t)]
    (; t, x, u, m)
end

if isfile(joinpath(pkgdir(@__MODULE__), "src", "trace.jl"))
    include("trace.jl")
end

end # module

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using DOIHS
using ValidatedNumerics
using Distributions
using Base.Test

@testset "CRRA utility" begin

    @test_throws DomainError crra_utility_function(0.1)
    @test_throws DomainError crra_utility(0.1, 0.1)
    @test_throws DomainError crra_utility(-1, 1)
    @test crra_utility(0, 1) ≡ -∞

    u = crra_utility_function(1)
    @test_throws DomainError u(-10)
    @test u(0) ≡ -∞
    @test u(1) == 0.0
    @test isa(u(1), Float64)

    u = crra_utility_function(2)
    @test_throws DomainError u(-10)
    @test u(0.0) ≡ -∞
    @test u(1) == 0.0
    @test isa(u(1), Float64)

    @test_throws DomainError crra_u′(-0.1, 0.1)
    @test crra_u′(2, 1) == 1/2

    @test crra_inv_u′(crra_u′(9.0, 3.0), 3.0) ≈ 9.0
end

@testset "function approximation" begin

    function approxerror(basis, f, n=100)
        α = basis_matrix(basis) \ f(collocation_points(basis))
        x = linspace(basis, n)
        y = basis_matrix(basis, x) * α
        maximum(abs(y-f(x)) for (x,y) in zip(x,y))
    end

    @test approxerror(ChebyshevBasis(10), exp) ≤ 1e-9

    @test approxerror(IntervalAB(0..pi, ChebyshevBasis(10)), sin) ≤ 1e-7

end

@testset "Chebyshev extrapolation" begin
    basis = IntervalAB(0..10, ChebyshevBasis(10; extrapolate = true))
    f(x) = x^2 - 5*x + 6
    g = basis \ f
    @test g(11) ≈ f(11)
    @test g(-1) ≈ f(-1)
end

@testset "linear interpolation" begin
    @test_throws Exception LinearInterpolation([2.0,1.0])
    li = LinearInterpolation(Float64[1,2,4])
    @test collocation_points(li) == [1.0,2.0,4.0]
    @test domain(li) == 1.0..4.0
    @test degf(li) == 3
    y = [5.0, 3.0, 7.0]
    lif = li \ y
    @test collocation_values(lif) == y
    @test lif(1.0) == 5.0
    @test lif(1.5) == 4.0
    @test lif(1.75) == 3.5
    @test lif(3.0) == 5.0
    @test lif(3.5) == 6.0
    @test lif(3.75) == 6.5
    @test lif(4.0) == 7.0
    @test_throws Exception lif(-1.0)
    @test_throws Exception lif(5.0)
end

@testset "extrapolate level" begin
    inner_basis = IntervalAB(1..10, ChebyshevBasis(3))
    basis = ExtrapolateLevel(inner_basis)
    @test collocation_points(inner_basis) == collocation_points(basis)
    @test degf(inner_basis) == degf(basis)
    y = [2.0, 3.0, 7.0]
    ipf = basis \ y
    @test collocation_values(ipf) ≈ y
    a = ipf(1.0)
    b = ipf(10.0)
    @test ipf(-1.0) ≈ a
    @test ipf(11.0) ≈ b
end

@testset "quadrature" begin

    d = Truncated(Normal(0,1), -1, 2)
    q = quadrature(10, d)
    m = mean(d)
    @test isapprox(q(identity), m; rtol=1e-6)
    @test isapprox(q(x->(x-m)^2), var(d); rtol=1e-6)

end

@testset "quadrature normal" begin

    μ = 1.0
    σ = 2.0
    q = quadrature_normal(10, μ, σ)

    n = Normal(μ,σ)

    @test mean(n) ≈ q(identity)
    @test var(n) ≈ q(x->(x-μ)^2)

end


text:
function filter_WF_adaptive_precomputation_keep_fixed_number(α, data, fixed_number::Int64; silence = false, trim0 = true)
    # println("filter_WF_mem2")

    function prune_keeping_fixed_number(Λ_of_t, wms_of_t)
        Λ_of_t_kept, wms_of_t_kept = keep_fixed_number_of_weights(Λ_of_t, wms_of_t, fixed_number)
        return Λ_of_t_kept, normalise(wms_of_t_kept)
    end

    filter_WF_adaptive_precomputation_ar(α, data, prune_keeping_fixed_number; silence = silence, trim0 = trim0)

end

function filter_WF_adaptive_precomputation_keep_above_threshold(α, data, ε::Float64; silence = false, trim0 = true)
    # println("filter_WF_mem2")

    function prune_keeping_above_threshold(Λ_of_t, wms_of_t)
        Λ_of_t_kept, wms_of_t_kept = keep_above_threshold(Λ_of_t, wms_of_t, ε)
        return Λ_of_t_kept, normalise(wms_of_t_kept)
    end

    filter_WF_adaptive_precomputation_ar(α, data, prune_keeping_above_threshold; silence = silence, trim0 = trim0)

end

function filter_WF_adaptive_precomputation_keep_fixed_fraction(α, data, fraction::Float64; silence = false, trim0 = true)
    # println("filter_WF_mem2")

    function prune_keeping_fixed_fraction(Λ_of_t, wms_of_t)
        Λ_of_t_kept, wms_of_t_kept = keep_fixed_fraction(Λ_of_t, wms_of_t, fraction)
        return Λ_of_t_kept, normalise(wms_of_t_kept)
    end

    filter_WF_adaptive_precomputation_ar(α, data, prune_keeping_fixed_fraction; silence = silence, trim0 = trim0)

end

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
# wrappers of low-level functionality

function curandCreateGenerator(typ)
  handle_ref = Ref{curandGenerator_t}()
  @check unsafe_curandCreateGenerator(handle_ref, typ) CURAND_STATUS_INITIALIZATION_FAILED
  handle_ref[]
end

@memoize function curandGetProperty(property::libraryPropertyType)
  value_ref = Ref{Cint}()
  curandGetProperty(property, value_ref)
  value_ref[]
end

@memoize version() = VersionNumber(curandGetProperty(CUDA.MAJOR_VERSION),
                                   curandGetProperty(CUDA.MINOR_VERSION),
                                   curandGetProperty(CUDA.PATCH_LEVEL))


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 5934e838-c887-11eb-3d34-5364828175f3
using JLD2, FFTW

# ╔═╡ 6490f3c9-828c-4a61-a21c-201ca2134772
using StatsBase

# ╔═╡ c9739713-12b6-43a7-afe2-6e4fea25b732
using DataFrames

# ╔═╡ 0ff69784-d617-476d-8793-7d82d8f1ef38
using Statistics

# ╔═╡ eb8fbe5e-edf9-452b-a83c-37b8d5258344
using CairoMakie

# ╔═╡ 5445a823-f610-47e0-81da-ee8a302271e1
using PlutoUI

# ╔═╡ 5dfb64de-c28c-4aa7-abca-80933117ce09
using Distributions

# ╔═╡ 23956119-4eca-490d-8f4e-bb001ecbdd7d
using Latexify

# ╔═╡ fd498b55-6d68-4418-8b50-de5cef64931a
using Clustering

# ╔═╡ b1b25bed-79ac-442a-8c10-985d8d40fdea
using CSV, DataFramesMeta

# ╔═╡ 46c9dbc9-6a9e-4393-a251-4b2bd1d7e638
using DelimitedFiles

# ╔═╡ 7338aa8b-63c8-440f-9e6b-172965b27c23
using HypothesisTests

# ╔═╡ b784697f-0644-40ef-97b7-e6fe4af4188e
resultdir = joinpath("..", "results")

# ╔═╡ 66571363-4357-4de9-b932-2731aa3a6b42
jlds = joinpath.(resultdir, filter(endswith("jld2"), readdir(resultdir)))

# ╔═╡ 909914e4-a7c1-4479-801d-9f92ccefafc4
function getlastrice(f)
	rices = jldopen(f) do io
		map(1:1000) do i
			df = io[string(i)]
			df.food[end] ≥ df.food[1] ÷ 2
		end
	end
	# convert.(Int32, rices)
end

# ╔═╡ 20b53174-e6ab-4035-a1aa-a1ad6a8b100b
readdir(resultdir)

# ╔═╡ a0d379a9-75be-4c18-87a3-7410d2bc912c
lastrices = getlastrice.(jlds)

# ╔═╡ ff5d4c71-b9dd-4028-8c92-3951a4ebcf3c
lastrices_ = filter(!isempty, convert.(Vector{Int32}, lastrices))

# ╔═╡ 011014c1-ac43-44a8-a358-c9d86ddeb6bf
stds = map(r -> std(r, corrected=false), lastrices_)

# ╔═╡ 26862094-b78f-4984-89ef-aa25177c0492
means = map(mean, lastrices_)

# ╔═╡ c09ebe2c-e31a-49fd-9b36-2b653f682db5
function cvs(X)
	map(100:1000) do n
		std(X[1:n]) / mean(X[1:n])
	end
end

# ╔═╡ da80c530-f5be-4a52-937b-abbf366fff67
function dcvs(X)
	C = cvs(X)
	[abs(C[m] - C[n]) for m in 1:length(C), n in 1:length(C) if n > m]
end

# ╔═╡ 6d8c724a-cf10-4938-982c-7bff67480260
findmax(dcvs(lastrices[1]))

# ╔═╡ 5f56093d-3c6f-46a0-85ce-e0f5a6fc792f
function sliding(z, w)
	((@view z[i:i+w-1]) for i in 1:length(z)-w+1)
end

# ╔═╡ 7c96d734-9876-4215-810a-854afd73d90c
map(sum, sliding(1:10, 3))

# ╔═╡ 9477c96b-8a6b-49c0-8a22-ed91cebf3e30
jlds[1]

# ╔═╡ 8f252bb5-ae7e-4de3-a956-94ebfc056e1e
jldopen(jlds[1]) do f
	f["500"][900:905, :]
end

# ╔═╡ 38fb200c-b683-4b02-9341-ecf6d0a3010f
@bind i PlutoUI.Slider(1:length(lastrices))

# ╔═╡ d967bfe9-da66-47be-9ecc-469ce324d633
wns = let X = lastrices[i]
	cs = map(sliding(X, 10)) do w
		std(w) / mean(X)
	end
	# map(2:length(cs)) do n
	# 	cov(cs[1:n])
	# end
end

# ╔═╡ daef4589-5ccf-4639-9f0d-4ecc15a1cb90
lines((wns / maximum(wns)))

# ╔═╡ b23c33cf-dffc-4bbd-b088-d51d67ec37df
function nmin_coeff(X, g; alpha=0.05, beta=0.95)
	V = length(X) - 1
	t = TDist(V)
	s = std(X)
	2 * s^2 / g * (quantile(t, 1 - alpha/2) + quantile(t, 1-beta/2))^2
end

# ╔═╡ 301520f9-1658-4e8a-b39b-7854ee6795d4
map(n -> nmin_coeff(lastrices[1][1:n], 5000), 2:1000)

# ╔═╡ 16375ff7-e018-4948-95c2-5c152afb3d0a
kmeans(lastrices[2] |> transpose, 2)

# ╔═╡ 02eae0af-570f-4a1d-9035-ba9be3efd7ef
const testresultdir = joinpath("..", "kiem-dinh")

# ╔═╡ 4068822f-5bd2-4b13-8919-18d5202ba54c
PlutoUI.TableOfContents()

# ╔═╡ 1b672c30-d567-4d1f-97e4-e462a314e434
md"# Sinh bảng"

# ╔═╡ 740dc4b4-4240-48fa-a164-4097353c48e2
testcsvs = joinpath.(testresultdir, readdir(testresultdir))

# ╔═╡ ee016a15-a598-451f-94c9-2b19f9c82673
@bind testcsv PlutoUI.Select(testcsvs)

# ╔═╡ 2c11fd74-711e-4781-b127-7fa01cfd5ed4
df = CSV.read(testcsv, DataFrame)

# ╔═╡ eea51e97-9f4e-41b5-b8f4-764e43bcb04f
function formaltable(df)
	s = latextabular(df; latex=false)
	lines = split(s, "\n")
	tbegin = replace(lines[1], "" => "")
	# header
	thead = strip(replace(lines[2], "\\\\" => ""))
	# thead = join("\\textbf{" .* split(thead, r" *& *") .* "}", " & ")
	tbody = join(lines[3:end-2], "\n\\midrule\n")
	s = """
	$(tbegin)
	\\toprule
	$(thead) \\\\
	\\midrule
	\\midrule
	$(tbody)
	\\bottomrule
	\\end{tabular}
	"""
	s = replace(s, r"([^a-zA-Z_.])([0-9]\.?[0-9]*)([^a-zA-Z_.])" => s"\1$\2$\3")
	# s = replace(s, r"([0-9]+\.[0-9]+)" => s"$\1$")
	s = replace(s, "1x2" => raw"\OneByTwo")
	s = replace(s, "3x3" => raw"\ThreeByThree")
	s = replace(s, "false" => raw"\Reject")
	s = replace(s, "true" => raw"\Accept")
	s = replace(s, "corner" => raw"\Corner")
	s = replace(s, "border" => raw"\Border")
	clipboard(s)
	return s
end

# ╔═╡ 41b5e3d7-5f39-4b55-b6be-3520c76c0924
function clean(df; ditch, group, key, val, newnames = nothing, sortcols=nothing)
	df = select(df, Not(ditch))
	function parseenvmap(m)
		ma = match(r"(\d+)-(\dx\d).csv", m)
		parse(Int, ma[1]), ma[2]
	end
	df = transform(df, :envmap => ByRow(parseenvmap) => [:flower, :layout])
	df = select(df, Not(:envmap))
	df = [unstack(d, key, val) for d in groupby(df, group)]
	df = vcat(df...)

	if !isnothing(sortcols)
		df = sort(df, sortcols)
		sortcols = String.(sortcols)
		othercols = setdiff(names(df), sortcols)
		df = select(df, [sortcols; othercols])
	end
	if !isnothing(newnames)
		df = rename(df, newnames)
	end
	return df
end

# ╔═╡ ead4bf05-6319-49aa-bcf9-6a5bb50469ea
sortcols = [:pr_killed, :layout, :init_position, :nb_bph_init]

# ╔═╡ ec65b6c9-6851-470a-a8cf-1db53e4f143b
newnames=[
	:pr_killed => raw"$\Pr_{\text{kill}}$",
	:init_position => raw"$P_0$",
	:layout => raw"$L$",
	:nb_bph_init => raw"$N$",
]

# ╔═╡ fa49c376-e47c-4089-9d40-f1237c74f749
function clean_for_h0(df)
	clean(df;
		ditch=[:k, :n],
		group=[:pr_killed, :layout],
		key=:flower,
		val=:accept,
		newnames=newnames,
		sortcols=sortcols,
	)
end

# ╔═╡ 3df31e23-605d-4833-ae3b-53f1770a487b
function clean_for_k(df)
	clean(df;
		ditch=[:accept, :n],
		group=[:pr_killed, :layout],
		key=:flower,
		val=:k,
		newnames=newnames,
		sortcols = sortcols
	)
end

# ╔═╡ c0780832-f855-447a-a03b-6327e08261c8
names(df)

# ╔═╡ 9f89b857-9dc2-45b2-a6e4-a9de4ec24370
nicedf = let df = df#@where(df, :pr_killed .== 0.12)
	s = clean_for_k(df) |> formaltable
end

# ╔═╡ 3fa30047-885d-4710-a291-c969ead989c9
with_terminal() do
	print(nicedf)
end

# ╔═╡ 8e0ba929-10fd-4586-96e1-28f59490aa69
dfs = let csvs = map(testcsvs) do f
		last(splitpath(f))
		m = match(r"test-p0\+([0-9.]+)-alpha\+([0-9.]+).csv", f)
		(
			p0 = parse(Float32, m[1]),
			alpha =parse(Float32, m[2]),
			path = f
		)
	end
	DataFrame(csvs)
end

# ╔═╡ 99b16a23-3eb1-46d4-a4e5-21eab89f3f52
let texpath = joinpath(@__DIR__, "..", "tex")
	mkpath(texpath)
	map(eachrow(dfs)) do row
		df = CSV.read(row.path, DataFrame)
		map(unique(df.pr_killed)) do pk
			df_ = @where(df, :pr_killed .== pk)
			df_ = clean_for_h0(df_)
			s = formaltable(df_)
			write(joinpath(texpath, "p+$pk-p0+$(row.p0)-alpha-$(row.alpha).tex"), s)
		end
	end
	df = CSV.read(dfs[1, :path], DataFrame)
	map(unique(df.pr_killed)) do pk
		df_ = @where(df, :pr_killed .== pk)
		df_ = clean_for_k(df_)
		s = formaltable(df_)
		write(joinpath(texpath, "p-$pk.tex"), s)
	end
end

# ╔═╡ 3b71c9e7-bd76-4a3a-ba13-e0be6072f417
let n = 1000
	d = Binomial(n, 0.001)
	x = 0:n
	y = pdf.(d, x)
	lines(x, y)
end

# ╔═╡ b9d9c8d9-7ad2-4acb-b10f-b69ddd8db2d6
cdf(Binomial(1000, 0.01), 9)

# ╔═╡ aad68c8f-3575-4d9a-8335-4f1e5339b574
t1 = BinomialTest(13, 1000, 0.001)

# ╔═╡ b8e982f6-a7ea-4031-b971-ce9bf2b0fb88
@bind p0 PlutoUI.Select([0.05, 0.01, 0.005, 0.001, 0] .|> string)

# ╔═╡ 24e13415-cee9-4408-be4f-c00d3f17d01e
# H0, hoa ok p <= p0
# Ha, hoa ko ok, p >= p0
# k: so lan hoa ok
# n: tong so lan (1000)

# ╔═╡ 08f01fc4-f0c0-4a5d-a733-eab8296a1618

# ╔═╡ 8a18a998-b6b6-4d30-8132-8a3dfd5260f1
@bind n PlutoUI.Slider(0:1000)

# ╔═╡ 6f1630ff-16c9-4721-9707-97fb60479914
t2 = BinomialTest(n, 1000, parse(Float64, p0))

# ╔═╡ afd473a4-0aa6-4bed-9f4a-d36ad86d1fa5
function is_flower_effective(t, alpha)
	T = t.x / t.n
	a, b = confint(t2; level=1-alpha, tail=:both)
	(T, a, b)
end

# ╔═╡ 78d7cd1c-c05c-497e-a31d-c79efc2e7c2e
@bind α PlutoUI.Select(string.([0.1, 0.05, 0.01, 0.005]))

# ╔═╡ 9446695b-59bf-4d36-bb1b-b66456eb5ac7
pvalue(t2; tail=:left)

# ╔═╡ a4855f8d-9182-4870-ab98-746b89df2017
let alphas = [0.1, 0.05, 0.01, 0.005]
	cons = [confint(t2; level=1-alpha, tail=:left) for alpha in alphas]
	as = getindex.(cons, 1)
	bs = getindex.(cons, 2)
	f = Figure()
	axis = Axis(f[1, 1])
	lines!(axis, as)
	lines!(axis, bs)
	f
end

# ╔═╡ ca31150c-8ffc-4547-8212-c277cc08163b
is_flower_effective(t2, parse(Float64, α))

# ╔═╡ b4ee62b9-9c93-4129-8146-929ac4690d99
t2.x / t2.n

# ╔═╡ 36544591-66ef-4df4-b1b6-728303ed7e7a
confint(t2; level=0.95, tail=:right)

# ╔═╡ fc715704-10a7-45b4-ba01-e2daefd1d5b8
pvalue(t1, tail=:left)

# ╔═╡ 4057976d-a3b2-4ba1-813b-7830c0341a71
pvalue(t2, tail=:left)

# ╔═╡ 03687c67-d14d-480f-b721-3b4598544560
confint(t2, tail=:left)

# ╔═╡ ff70e6d1-df98-4444-ad1c-40ba830ff348
isgood(n, 1000, 0.001, α)

# ╔═╡ e1d0fad7-9657-4085-b8eb-9fa522a70e20
OneSampleZTest(0.004, 0.001)

# ╔═╡ Cell order:
# ╠═5934e838-c887-11eb-3d34-5364828175f3
# ╠═6490f3c9-828c-4a61-a21c-201ca2134772
# ╠═c9739713-12b6-43a7-afe2-6e4fea25b732
# ╠═0ff69784-d617-476d-8793-7d82d8f1ef38
# ╠═b784697f-0644-40ef-97b7-e6fe4af4188e
# ╠═66571363-4357-4de9-b932-2731aa3a6b42
# ╠═909914e4-a7c1-4479-801d-9f92ccefafc4
# ╠═20b53174-e6ab-4035-a1aa-a1ad6a8b100b
# ╠═a0d379a9-75be-4c18-87a3-7410d2bc912c
# ╠═ff5d4c71-b9dd-4028-8c92-3951a4ebcf3c
# ╠═011014c1-ac43-44a8-a358-c9d86ddeb6bf
# ╠═26862094-b78f-4984-89ef-aa25177c0492
# ╠═c09ebe2c-e31a-49fd-9b36-2b653f682db5
# ╠═6d8c724a-cf10-4938-982c-7bff67480260
# ╠═da80c530-f5be-4a52-937b-abbf366fff67
# ╠═5f56093d-3c6f-46a0-85ce-e0f5a6fc792f
# ╠═7c96d734-9876-4215-810a-854afd73d90c
# ╠═9477c96b-8a6b-49c0-8a22-ed91cebf3e30
# ╠═8f252bb5-ae7e-4de3-a956-94ebfc056e1e
# ╠═eb8fbe5e-edf9-452b-a83c-37b8d5258344
# ╠═5445a823-f610-47e0-81da-ee8a302271e1
# ╠═38fb200c-b683-4b02-9341-ecf6d0a3010f
# ╠═d967bfe9-da66-47be-9ecc-469ce324d633
# ╠═daef4589-5ccf-4639-9f0d-4ecc15a1cb90
# ╠═5dfb64de-c28c-4aa7-abca-80933117ce09
# ╠═23956119-4eca-490d-8f4e-bb001ecbdd7d
# ╠═b23c33cf-dffc-4bbd-b088-d51d67ec37df
# ╠═301520f9-1658-4e8a-b39b-7854ee6795d4
# ╠═fd498b55-6d68-4418-8b50-de5cef64931a
# ╠═16375ff7-e018-4948-95c2-5c152afb3d0a
# ╠═b1b25bed-79ac-442a-8c10-985d8d40fdea
# ╠═02eae0af-570f-4a1d-9035-ba9be3efd7ef
# ╠═4068822f-5bd2-4b13-8919-18d5202ba54c
# ╠═1b672c30-d567-4d1f-97e4-e462a314e434
# ╠═740dc4b4-4240-48fa-a164-4097353c48e2
# ╠═ee016a15-a598-451f-94c9-2b19f9c82673
# ╠═2c11fd74-711e-4781-b127-7fa01cfd5ed4
# ╠═eea51e97-9f4e-41b5-b8f4-764e43bcb04f
# ╠═41b5e3d7-5f39-4b55-b6be-3520c76c0924
# ╠═fa49c376-e47c-4089-9d40-f1237c74f749
# ╠═ead4bf05-6319-49aa-bcf9-6a5bb50469ea
# ╠═ec65b6c9-6851-470a-a8cf-1db53e4f143b
# ╠═3df31e23-605d-4833-ae3b-53f1770a487b
# ╠═c0780832-f855-447a-a03b-6327e08261c8
# ╠═9f89b857-9dc2-45b2-a6e4-a9de4ec24370
# ╠═3fa30047-885d-4710-a291-c969ead989c9
# ╠═99b16a23-3eb1-46d4-a4e5-21eab89f3f52
# ╠═8e0ba929-10fd-4586-96e1-28f59490aa69
# ╠═46c9dbc9-6a9e-4393-a251-4b2bd1d7e638
# ╠═3b71c9e7-bd76-4a3a-ba13-e0be6072f417
# ╠═b9d9c8d9-7ad2-4acb-b10f-b69ddd8db2d6
# ╠═7338aa8b-63c8-440f-9e6b-172965b27c23
# ╠═aad68c8f-3575-4d9a-8335-4f1e5339b574
# ╠═6f1630ff-16c9-4721-9707-97fb60479914
# ╠═afd473a4-0aa6-4bed-9f4a-d36ad86d1fa5
# ╠═b8e982f6-a7ea-4031-b971-ce9bf2b0fb88
# ╠═24e13415-cee9-4408-be4f-c00d3f17d01e
# ╠═08f01fc4-f0c0-4a5d-a733-eab8296a1618
# ╠═8a18a998-b6b6-4d30-8132-8a3dfd5260f1
# ╠═78d7cd1c-c05c-497e-a31d-c79efc2e7c2e
# ╠═9446695b-59bf-4d36-bb1b-b66456eb5ac7
# ╠═a4855f8d-9182-4870-ab98-746b89df2017
# ╠═ca31150c-8ffc-4547-8212-c277cc08163b
# ╠═b4ee62b9-9c93-4129-8146-929ac4690d99
# ╠═36544591-66ef-4df4-b1b6-728303ed7e7a
# ╠═fc715704-10a7-45b4-ba01-e2daefd1d5b8
# ╠═4057976d-a3b2-4ba1-813b-7830c0341a71
# ╠═03687c67-d14d-480f-b721-3b4598544560
# ╠═ff70e6d1-df98-4444-ad1c-40ba830ff348
# ╠═e1d0fad7-9657-4085-b8eb-9fa522a70e20

max_stars_repo_name: ndgnuh/RiceBPH

max_issues_repo_name: ndgnuh/RiceBPH

max_forks_repo_name: ndgnuh/RiceBPH


text:
##############################################################################
##
## Sorting
##
##############################################################################

#########################################
## Permutation & Ordering types/functions
#########################################
# Sorting in Julia works through Orderings, where each ordering is a type
# which defines a comparison function lt(o::Ord, a, b).

# UserColOrdering: user ordering of a column; this is just a convenience container
#                  which allows a user to specify column specific orderings
#                  with "order(column, rev=true,...)"

type UserColOrdering{T<:ColumnIndex}
    col::T
    kwargs
end

# This is exported, and lets a user define orderings for a particular column
order{T<:ColumnIndex}(col::T; kwargs...) = UserColOrdering{T}(col, kwargs)

# Allow getting the column even if it is not wrapped in a UserColOrdering
_getcol(o::UserColOrdering) = o.col
_getcol(x) = x

###
# Get an Ordering for a single column
###
function ordering(col_ord::UserColOrdering, lt::Function, by::Function, rev::Bool, order::Ordering)
    for (k,v) in col_ord.kwargs
        if     k == :lt;    lt    = v
        elseif k == :by;    by    = v
        elseif k == :rev;   rev   = v
        elseif k == :order; order = v
        else
            error("Unknown keyword argument: ", string(k))
        end
    end

    Order.ord(lt,by,rev,order)
end

ordering(col::ColumnIndex, lt::Function, by::Function, rev::Bool, order::Ordering) =
             Order.ord(lt,by,rev,order)

# DFPerm: defines a permutation on a particular DataFrame, using
#         a single ordering (O<:Ordering) or a list of column orderings
#         (O<:AbstractVector{Ordering}), one per DataFrame column
#
#         If a user only specifies a few columns, the DataFrame
#         contained in the DFPerm only contains those columns, and
#         the permutation induced by this ordering is used to
#         sort the original (presumably larger) DataFrame

type DFPerm{O<:@compat(Union{Ordering, AbstractVector}), DF<:AbstractDataFrame} <: Ordering
    ord::O
    df::DF
end

function DFPerm{O<:Ordering}(ords::AbstractVector{O}, df::AbstractDataFrame)
    if length(ords) != ncol(df)
        error("DFPerm: number of column orderings does not equal the number of DataFrame columns")
    end
    DFPerm{AbstractVector{O}, typeof(df)}(ords, df)
end

DFPerm{O<:Ordering}(o::O, df::AbstractDataFrame) = DFPerm{O,typeof(df)}(o,df)

# For sorting, a and b are row indices (first two lt definitions)
# For issorted, the default row iterator returns DataFrameRows instead,
# so two more lt function is defined below
function Sort.lt{V<:AbstractVector}(o::DFPerm{V}, a, b)
    for i = 1:ncol(o.df)
        if lt(o.ord[i], o.df[a,i], o.df[b,i])
            return true
        end
        if lt(o.ord[i], o.df[b,i], o.df[a,i])
            return false
        end
    end
    false
end

function Sort.lt{O<:Ordering}(o::DFPerm{O}, a, b)
    for i = 1:ncol(o.df)
        if lt(o.ord, o.df[a,i], o.df[b,i])
            return true
        end
        if lt(o.ord, o.df[b,i], o.df[a,i])
            return false
        end
    end
    false
end

function Sort.lt{V<:AbstractVector}(o::DFPerm{V}, a::DataFrameRow, b::DataFrameRow)
    for i = 1:ncol(o.df)
        if lt(o.ord[i], a[i], b[i])
            return true
        end
        if lt(o.ord[i], b[i], a[i])
            return false
        end
    end
    false
end

function Sort.lt{O<:Ordering}(o::DFPerm{O}, a::DataFrameRow, b::DataFrameRow)
    for i = 1:ncol(o.df)
        if lt(o.ord, a[i], b[i])
            return true
        end
        if lt(o.ord, b[i], a[i])
            return false
        end
    end
    false
end

###
# Get an Ordering for a DataFrame
###

################
## Case 1: no columns requested, so sort by all columns using requested order
################
## Case 1a: single order
######
ordering(df::AbstractDataFrame, lt::Function, by::Function, rev::Bool, order::Ordering) =
    DFPerm(Order.ord(lt, by, rev, order), df)

######
## Case 1b: lt, by, rev, and order are Arrays
######
function ordering{S<:Function, T<:Function}(df::AbstractDataFrame,
                                            lt::AbstractVector{S}, by::AbstractVector{T},
                                            rev::AbstractVector{Bool}, order::AbstractVector)
    if !(length(lt) == length(by) == length(rev) == length(order) == size(df,2))
        throw(ArgumentError("Orderings must be specified for all DataFrame columns"))
    end
    DFPerm([Order.ord(_lt, _by, _rev, _order) for (_lt, _by, _rev, _order) in zip(lt, by, rev, order)], df)
end

################
## Case 2:  Return a regular permutation when there's only one column
################
## Case 2a: The column is given directly
######
ordering(df::AbstractDataFrame, col::ColumnIndex, lt::Function, by::Function, rev::Bool, order::Ordering) =
    Perm(Order.ord(lt, by, rev, order), df[col])

######
## Case 2b: The column is given as a UserColOrdering
######
ordering(df::AbstractDataFrame, col_ord::UserColOrdering, lt::Function, by::Function, rev::Bool, order::Ordering) =
    Perm(ordering(col_ord, lt, by, rev, order), df[col_ord.col])

################
## Case 3:  General case: cols is an iterable of a combination of ColumnIndexes and UserColOrderings
################
## Case 3a: None of lt, by, rev, or order is an Array
######
function ordering(df::AbstractDataFrame, cols::AbstractVector, lt::Function, by::Function, rev::Bool, order::Ordering)

    if length(cols) == 0
        return ordering(df, lt, by, rev, order)
    end

    if length(cols) == 1
        return ordering(df, cols[1], lt, by, rev, order)
    end

    # Collect per-column ordering info

    ords = Ordering[]
    newcols = Int[]

    for col in cols
        push!(ords, ordering(col, lt, by, rev, order))
        push!(newcols, index(df)[(_getcol(col))])
    end

    # Simplify ordering when all orderings are the same
    if all([ords[i] == ords[1] for i = 2:length(ords)])
        return DFPerm(ords[1], df[newcols])
    end

    return DFPerm(ords, df[newcols])
end

######
# Case 3b: cols, lt, by, rev, and order are all arrays
######
function ordering{S<:Function, T<:Function}(df::AbstractDataFrame, cols::AbstractVector,
                                            lt::AbstractVector{S}, by::AbstractVector{T},
                                            rev::AbstractVector{Bool}, order::AbstractVector)

    if !(length(lt) == length(by) == length(rev) == length(order))
        throw(ArgumentError("All ordering arguments must be 1 or the same length."))
    end

    if length(cols) == 0
        return ordering(df, lt, by, rev, order)
    end

    if length(lt) != length(cols)
        throw(ArgumentError("All ordering arguments must be 1 or the same length as the number of columns requested."))
    end

    if length(cols) == 1
        return ordering(df, cols[1], lt[1], by[1], rev[1], order[1])
    end

    # Collect per-column ordering info

    ords = Ordering[]
    newcols = Int[]

    for i in 1:length(cols)
        push!(ords, ordering(cols[i], lt[i], by[i], rev[i], order[i]))
        push!(newcols, index(df)[(_getcol(cols[i]))])
    end

    # Simplify ordering when all orderings are the same
    if all([ords[i] == ords[1] for i = 2:length(ords)])
        return DFPerm(ords[1], df[newcols])
    end

    return DFPerm(ords, df[newcols])
end

######
## At least one of lt, by, rev, or order is an array or tuple, so expand all to arrays
######
function ordering(df::AbstractDataFrame, cols::AbstractVector, lt, by, rev, order)
    to_array(src::AbstractVector, dims) = src
    to_array(src::Tuple, dims) = [src...]
    to_array(src, dims) = fill(src, dims)

    dims = length(cols) > 0 ? length(cols) : size(df,2)
    ordering(df, cols,
             to_array(lt, dims),
             to_array(by, dims),
             to_array(rev, dims),
             to_array(order, dims))
end

#### Convert cols from tuple to Array, if necessary
ordering(df::AbstractDataFrame, cols::Tuple, args...) = ordering(df, [cols...], args...)

###########################
# Default sorting algorithm
###########################

# TimSort is fast for data with structure, but only if the DataFrame is large enough
# TODO: 8192 is informed but somewhat arbitrary

Sort.defalg(df::AbstractDataFrame) = size(df, 1) < 8192 ? Sort.MergeSort : SortingAlgorithms.TimSort

# For DataFrames, we can choose the algorithm based on the column type and requested ordering
function Sort.defalg{T<:Real}(df::AbstractDataFrame, ::Type{T}, o::Ordering)
    # If we're sorting a single numerical column in forward or reverse,
    # RadixSort will generally be the fastest stable sort
    if isbits(T) && sizeof(T) <= 8 && (o==Order.Forward || o==Order.Reverse)
        SortingAlgorithms.RadixSort
    else
        Sort.defalg(df)
    end
end
Sort.defalg(df::AbstractDataFrame,        ::Type,            o::Ordering) = Sort.defalg(df)
Sort.defalg(df::AbstractDataFrame, col    ::ColumnIndex,     o::Ordering) = Sort.defalg(df, eltype(df[col]), o)
Sort.defalg(df::AbstractDataFrame, col_ord::UserColOrdering, o::Ordering) = Sort.defalg(df, col_ord.col, o)
Sort.defalg(df::AbstractDataFrame, cols,                     o::Ordering) = Sort.defalg(df)

function Sort.defalg(df::AbstractDataFrame, o::Ordering; alg=nothing, cols=[])
    alg != nothing && return alg
    Sort.defalg(df, cols, o)
end

########################
## Actual sort functions
########################

Base.issorted(df::AbstractDataFrame; cols=Any[], lt=isless, by=identity, rev=false, order=Forward) =
    issorted(eachrow(df), ordering(df, cols, lt, by, rev, order))

# sort and sortperm functions

for s in [:(Base.sort), :(Base.sortperm)]
    @eval begin
        function $s(df::AbstractDataFrame; cols=Any[], alg=nothing,
                    lt=isless, by=identity, rev=false, order=Forward)
            if !(isa(by, Function) || eltype(by) <: Function)
                msg = "'by' must be a Function or a vector of Functions. Perhaps you wanted 'cols'."
                throw(ArgumentError(msg))
            end
            ord = ordering(df, cols, lt, by, rev, order)
            _alg = Sort.defalg(df, ord; alg=alg, cols=cols)
            $s(df, _alg, ord)
        end
    end
end

Base.sort(df::AbstractDataFrame, a::Algorithm, o::Ordering) = df[sortperm(df, a, o),:]
Base.sortperm(df::AbstractDataFrame, a::Algorithm, o::@compat(Union{Perm,DFPerm})) = sort!([1:size(df, 1);], a, o)
Base.sortperm(df::AbstractDataFrame, a::Algorithm, o::Ordering) = sortperm(df, a, DFPerm(o,df))


text:
numbers = readlines("day3/numbers.txt")
numbers = permutedims(hcat(collect.(numbers)...))
numbers = parse.(Int, numbers)

function find_rates(numbers)
    γ_bits = sum(numbers, dims = 1) .> size(numbers, 1) / 2
    ε_bits = map(~, γ_bits)  # bitwise not
    γ = parse(Int, join(convert(Array{Int}, γ_bits)), base = 2)
    ε = parse(Int, join(convert(Array{Int}, ε_bits)), base = 2)
    γ, ε
end

γ, ε = find_rates(numbers)
p1 = γ * ε

@show p1

function find_rating(numbers, keep_bit)
    more_numbers = true
    bit_idx = 1
    while more_numbers
        bits = numbers[:, bit_idx]
        numbers = numbers[bits .== keep_bit(bits), :]
        more_numbers = size(numbers, 1) > 1
        bit_idx += 1
    end
    parse(Int, join(numbers), base = 2)
end

most_common_bit(bits) = sum(bits) >= length(bits) / 2 ? 1 : 0

least_common_bit(bits) = sum(bits) >= length(bits) / 2 ? 0 : 1

oxy = find_rating(numbers, most_common_bit)
co2 = find_rating(numbers, least_common_bit)
p2 = oxy * co2

@show p2


text:
using Test, YaoBlocks

@testset "test constant gates" begin
    include("const_gate.jl")
end

@testset "test phase gate" begin
    include("phase_gate.jl")
end

@testset "test shift gate" begin
    include("shift_gate.jl")
end

@testset "test rotation gate" begin
    include("rotation_gate.jl")
end

@testset "test time evolution" begin
    include("time_evolution.jl")
end

@testset "test general matrix gate" begin
    include("general_matrix_gate.jl")
end

@testset "test math gate" begin
    include("math_gate.jl")
end

@testset "test reflect gate" begin
    include("reflect_gate.jl")
end

@testset "test measure" begin
    include("measure.jl")
end

# it does nothing
@test chsubblocks(X, Y) === X


text:
function textDocument_signatureHelp_request(params::TextDocumentPositionParams, server::LanguageServerInstance, conn)
    doc = getdocument(server, URI2(params.textDocument.uri))
    sigs = SignatureInformation[]
    offset = get_offset(doc, params.position)
    x = get_expr(getcst(doc), offset)
    arg = 0
    if x isa EXPR && parentof(x) isa EXPR && CSTParser.iscall(parentof(x))
        if CSTParser.isidentifier(parentof(x).args[1])
            call_name = parentof(x).args[1]
        elseif CSTParser.iscurly(parentof(x).args[1]) && CSTParser.isidentifier(parentof(x).args[1].args[1])
            call_name = parentof(x).args[1].args[1]
        elseif CSTParser.is_getfield_w_quotenode(parentof(x).args[1])
            call_name = parentof(x).args[1].args[2].args[1]
        else
            call_name = nothing
        end
        if call_name !== nothing && (f_binding = refof(call_name)) !== nothing && (tls = StaticLint.retrieve_toplevel_scope(call_name)) !== nothing
            get_signatures(f_binding, tls, sigs, server)
        end
    end
    if (isempty(sigs) || (headof(x) === :RPAREN))
        return SignatureHelp(SignatureInformation[], 0, 0)
    end

    if headof(x) === :LPAREN
        arg = 0
    else
        arg = sum(headof(a) === :COMMA for a in parentof(x).trivia)
    end
    return SignatureHelp(filter(s -> length(s.parameters) > arg, sigs), 0, arg)
end

function get_signatures(b, tls::StaticLint.Scope, sigs::Vector{SignatureInformation}, server) end

function get_signatures(b::StaticLint.Binding, tls::StaticLint.Scope, sigs::Vector{SignatureInformation}, server)
    if b.val isa StaticLint.Binding
        get_signatures(b.val, tls, sigs, server)
    end
    if b.type == StaticLint.CoreTypes.Function || b.type == StaticLint.CoreTypes.DataType
        b.val isa SymbolServer.SymStore && get_signatures(b.val, tls, sigs, server)
        for ref in b.refs
            method = StaticLint.get_method(ref)
            if method !== nothing
                get_signatures(method, tls, sigs, server)
            end
        end
    end
end

function get_signatures(b::T, tls::StaticLint.Scope, sigs::Vector{SignatureInformation}, server) where T <: Union{SymbolServer.FunctionStore,SymbolServer.DataTypeStore}
    StaticLint.iterate_over_ss_methods(b, tls, server, function (m)
        push!(sigs, SignatureInformation(string(m), "", (a -> ParameterInformation(string(a[1]), string(a[2]))).(m.sig)))
        return false
    end)
end

function get_signatures(x::EXPR, tls::StaticLint.Scope, sigs::Vector{SignatureInformation}, server)
    if CSTParser.defines_function(x)
        sig = CSTParser.rem_where_decl(CSTParser.get_sig(x))
        params = ParameterInformation[]
        if sig isa EXPR && sig.args !== nothing
            for i = 2:length(sig.args)
                if (argbinding = bindingof(sig.args[i])) !== nothing
                    push!(params, ParameterInformation(valof(argbinding.name) isa String ? valof(argbinding.name) : "", missing))
                end
            end
            push!(sigs, SignatureInformation(string(Expr(sig)), "", params))
        end
    elseif CSTParser.defines_struct(x)
        args = x.args[3]
        if length(args) > 0
            if !any(CSTParser.defines_function, args.args)
                params = ParameterInformation[]
                for field in args.args
                    field_name = CSTParser.rem_decl(field)
                    push!(params, ParameterInformation(field_name isa EXPR && CSTParser.isidentifier(field_name) ? valof(field_name) : "", missing))
                end
                push!(sigs, SignatureInformation(string(Expr(x)), "", params))
            end
        end
    end
end


text:
@testset "games simulation" begin

MAX_ACTIONS_PER_GAME = 1000

SPIEL_GAMES_LIST = registered_games()

SPIEL_LOADABLE_GAMES_LIST = [
    g for g in SPIEL_GAMES_LIST if default_loadable(g)
]

@test length(SPIEL_LOADABLE_GAMES_LIST) >= 38

SPIEL_SIMULTANEOUS_GAMES_LIST = [
    g for g in SPIEL_LOADABLE_GAMES_LIST
    if dynamics(g) == SIMULTANEOUS
]

@test length(SPIEL_SIMULTANEOUS_GAMES_LIST) >= 14

SPIEL_MULTIPLAYER_GAMES_LIST = [
    (g, p)
    for g in SPIEL_LOADABLE_GAMES_LIST
    for p in max(min_num_players(g), 2) : min(max_num_players(g), 6)
    if (max_num_players(g) > 2) &&
        (max_num_players(g) > min_num_players(g)) &&
        (short_name(g) != "tiny_hanabi") &&  # default payoff only works for 2p
        (short_name(g) != "universal_poker")
]

@test length(SPIEL_MULTIPLAYER_GAMES_LIST) >= 35

function apply_action_test_clone(state, action)
    state_copy = copy(state)
    @test string(state) == string(state_copy)
    @test history(state) == history(state_copy)

    apply_action(state, action)
    apply_action(state_copy, action)

    @test string(state) == string(state_copy)
    @test history(state) == history(state_copy)
end

function serialize_deserialize(game, state)
    ser_str = serialize_game_and_state(game, state)
    new_game, new_state = deserialize_game_and_state(ser_str)
    @test string(game) == string(new_game)
    @test string(state) == string(new_state)
end

function simulate_game(game)
    @info "simulating game $(short_name(get_type(game)))"
    min_u, max_u = min_utility(game), max_utility(game)
    @test min_u < max_u

    state = new_initial_state(game)
    total_actions = 0

    next_serialize_check = 1

    while !is_terminal(state) && (total_actions <= MAX_ACTIONS_PER_GAME)
        total_actions += 1

        # Serialize/Deserialize is costly. Only do it every power of 2 actions.
        if total_actions >= next_serialize_check
            serialize_deserialize(game, state)
            next_serialize_check *= 2
        end

        # The state can be three different types: chance node,
        # simultaneous node, or decision node
        if is_chance_node(state)
            # Chance node: sample an outcome
            outcomes = chance_outcomes(state)
            @test length(outcomes) > 0
            action_list, prob_list = zip(outcomes...)
            action = action_list[sample(weights(collect(prob_list)))]
            apply_action(state, action)
        elseif is_simultaneous_node(state)
            chosen_actions = [
                rand(legal_actions(state, pid-1))
                for pid in 1:num_players(game)
            ]  # in julia, index starts with 1
            # Apply the joint action and test cloning states.
            apply_action_test_clone(state, chosen_actions)
        else
            # Decision node: sample action for the single current player
            action = rand(legal_actions(state, current_player(state)))
            # Apply action and test state cloning.
            apply_action_test_clone(state, action)
        end
    end

    @test total_actions > 0

    if is_terminal(state)
        # Check there are no legal actions.
        @test length(legal_actions(state)) == 0
        for player in 0:(num_players(game)-1)
            @test length(legal_actions(state, player)) == 0
        end

        utilities = returns(state)

        for u in utilities
            @test u >= min_utility(game)
            @test u <= max_utility(game)
        end

        @info "Simulation of game $game" total_actions utilities
    else
        @info "Simulation of game $game terminated after maximum number of actions $MAX_ACTIONS_PER_GAME"
    end
end

for game_info in SPIEL_LOADABLE_GAMES_LIST
    game = load_game(short_name(game_info))
    @test num_players(game) >= min_num_players(game_info)
    @test num_players(game) <= max_num_players(game_info)
    simulate_game(game)
end

for game_info in SPIEL_SIMULTANEOUS_GAMES_LIST
    converted_game = load_game_as_turn_based(short_name(game_info))
    simulate_game(converted_game)
end

for (game_info, n) in SPIEL_MULTIPLAYER_GAMES_LIST
    game = load_game(short_name(game_info); players=GameParameter(n))
    simulate_game(game)
end

simulate_game(load_game("breakthrough(rows=6,columns=6)"))
simulate_game(load_game("pig(players=2,winscore=15)"))

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 76e86802-19f0-11ec-0a0b-bfa6661be4f4
begin
    using Pkg
    Pkg.activate("..")
    using Javis, Random
    import Luxor
end

# ╔═╡ 630d8562-c784-4fc4-8180-222d913e5075
function ground(args...)
    background("blue")
    sethue("white")
end

# ╔═╡ a4477fbe-bdac-4e1c-8750-29fa62931968
mandel(z, c) = z^2 + c

# ╔═╡ 5be12534-6203-49a2-b64b-37076945bfaf
function mandelcheck(c, n, thr)
    out = mandel(0, c)
    k = 1
    while (k < n) | (abs(out) > thr)
        out = mandel(out, c)
        k += 1
    end
    out
end

# ╔═╡ a1b464f7-86c4-4eb4-933a-71a476c45c59
function draw_mandelbrot(width, height)
    thr = 100
    points = [
        Complex(i, j) for i in range(-2, 1 / 2, length = width),
        j in range(-1.5, 1.5, length = height)
    ]
    positions = CartesianIndices(points)
    origin(Point(0, 0))
    map(positions) do position
        if abs(mandelcheck(points[position], 20, thr)) < thr
            sethue(rand(["blue", "red", "orange"]))
            circle(Point(position[1], position[2]), 1, :fill)
        end
    end
end

# ╔═╡ 4eb94c31-dce3-4d09-8515-5a54ad077878
begin
    width = 500
    height = 500
    image = Luxor.@draw begin
        background("blue")
        draw_mandelbrot(500, 500)
    end width height
end

# ╔═╡ fc8bc0a2-602f-4700-9529-165692d044e2
function make_mandelbrot(width, height, n_frames, n_points; shape_lasting = 100)
    thr = 100
    points = [
        Complex(i, j) for i in range(-2, 1 / 2, length = n_points),
        j in range(-1.3, 1.3, length = n_points)
    ]
    positions = [
        (i, j) for i in range(1, height, length = size(points, 1)),
        j in range(1, width, length = size(points, 2))
    ]
    mask = abs.(mandelcheck.(points, 20, thr)) .< thr
    points = points[mask]
    positions = positions[mask]

    final_appearing_frame = n_frames - shape_lasting
    indices = shuffle(1:length(positions))
    map(zip(indices, positions, points)) do (idx, position, point)
        Object(
            floor(Int, rescale(idx, 2, length(positions), 1, final_appearing_frame)):n_frames,
            @JShape begin
                origin(O)
                sethue("orange")
                circle(Point(position...), 2, :fill)
            end
        )
    end
end

# ╔═╡ 15355ca2-64ca-41e9-b1ad-a3c4eaa14e01
let
    height, width = 500, 500
    n_frames = 500
    n_points = 200
    mandelvideo = Video(width, height)
    Background(1:n_frames, ground)
    make_mandelbrot(width, height, n_frames, n_points, shape_lasting = 100)
    render(mandelvideo, pathname = "../output/mandelbrot.gif")
end

# ╔═╡ Cell order:
# ╠═76e86802-19f0-11ec-0a0b-bfa6661be4f4
# ╠═630d8562-c784-4fc4-8180-222d913e5075
# ╠═a4477fbe-bdac-4e1c-8750-29fa62931968
# ╠═5be12534-6203-49a2-b64b-37076945bfaf
# ╠═4eb94c31-dce3-4d09-8515-5a54ad077878
# ╠═a1b464f7-86c4-4eb4-933a-71a476c45c59
# ╠═fc8bc0a2-602f-4700-9529-165692d044e2
# ╠═15355ca2-64ca-41e9-b1ad-a3c4eaa14e01


text:
# First-passage time distributions for diffusion models

# TODO:
# - add version that includes leak
# - add version that includes time-varying variance

# constant bound, constant drift
#
# algorithms adapted from
#    DJ Navarro & IG Fuss (2009). Fast and accurate calculations for
#    first-passage times in Wiener diffusion models. Journal of Mathematical
#    Psychology, 53(4), 222-230.

# Navarro & Fuss (2009), Eq. (13)
useshorttseries(t::Real, tol::Real) =
  2 + sqrt(-2t * log(2tol * sqrt(twoπ * t))) <
  sqrt(-2log(π * t * tol) / (t * abs2(π)))

abstract type PDFConstCacheBase end

struct AsymPDFConstCache <: PDFConstCacheBase
    c1::Float64
    c2::Float64
    c3::Float64
    c4::Float64
    w::Float64

    function AsymPDFConstCache(d::ConstDrift, b::ConstAsymBounds)
        mu, bu, bl = getmu(d), getubound(b), getlbound(b)
        return new(abs2(bu - bl), abs2(mu) / 2, mu * bu, mu * bl, -bl / (bu - bl))
    end
end

pdful(c::PDFConstCacheBase, t::Real, tol::Real=1.e-29) =
    pdfu(c, t, tol), pdfl(c, t, tol)
pdfu(c::AsymPDFConstCache, t::Real, tol::Real=1.e-29) =
    exp(c.c3 - c.c2 * t) / c.c1 * pdf_asymfastseries(t / c.c1, 1 - c.w, tol)
pdfl(c::AsymPDFConstCache, t::Real, tol::Real=1.e-29) =
    exp(c.c4 - c.c2 * t) / c.c1 * pdf_asymfastseries(t / c.c1, c.w, tol)

# fpt density for mu=0, constant bounds at 0 and 1, and starting point at w
pdf_asymfastseries(t::Real, w::Real, tol::Real) =
    t == 0.0 ? 0.0 :
        useshorttseries(t, tol) ? pdf_asymshortt(t, w, tol) :
                                  pdf_asymlongt(t, w, tol)

# Navarro & Fuss (2009), Eq. (6)
function pdf_asymshortt(t::Real, w::Real, tol::Real)
    b = t^-1.5 / sqrt2π
    tol *= b
    f = w * exp(-w * w / 2t)
    k = 1
    while true
        c = w + 2k
        incr = c * exp(-c * c / 2t)
        f += incr
        tol < abs(incr) || break
        c = w - 2k
        incr = c * exp(-c * c / 2t)
        f += incr
        tol < abs(incr) || break
        k += 1
    end
    return f * b
end

# Navarro & Fuss (2009), Eq. (5)
function pdf_asymlongt(t::Real, w::Real, tol::Real)
    tol *= π
    f, k = 0.0, 1
    while true
        incr = k * exp(-abs2(k * π) * t / 2) * sinpi(k * w)
        f += incr
        tol < abs(incr) || break
        k += 1
    end
    return f * π
end

struct SymPDFConstCache <: PDFConstCacheBase
    c1::Float64
    c2::Float64
    c3::Float64

    function SymPDFConstCache(d::ConstDrift, b::ConstSymBounds)
        mu, bound = getmu(d), getbound(b)
        new(4 * abs2(bound), abs2(mu) / 2, mu * bound)
    end
end

function pdful(c::SymPDFConstCache, t::Real, tol::Real=1.e-29)
    g = pdf_symfastseries(t / c.c1, tol) / c.c1
    return exp(c.c3 - c.c2 * t) * g, exp(- c.c3 - c.c2 * t) * g
end
pdfu(c::SymPDFConstCache, t::Real, tol::Real=1.e-29) =
    exp(c.c3 - c.c2 * t) / c.c1 * pdf_symfastseries(t / c.c1, tol)
pdfl(c::SymPDFConstCache, t::Real, tol::Real=1.e-29) =
    exp(- c.c3 - c.c2 * t) / c.c1 * pdf_symfastseries(t / c.c1, tol)

# fpt density for mu=0, constant bounds at 0 and 1, and starting point at 0.5,
pdf_symfastseries(t::Real, tol::Real) =
    t == 0.0 ? 0.0 :
        useshorttseries(t, tol) ?
            pdf_symseries(t, 1 / 8t, 1 / sqrt(8 * π * t^3), tol) :
            pdf_symseries(t, t * abs2(π) / 2, float(π), tol)

function pdf_symseries(t::Real, a::Real, b::Real, tol::Real)
    tol *= b
    f, twok = exp(-a), 3
    while true
        incr = twok * exp(- abs2(twok) * a)
        f -= incr
        tol < incr || break
        twok += 2
        incr = twok * exp(- abs2(twok) * a)
        f += incr
        tol < incr || break
        twok += 2
    end
    return b * f
end

PDFConstCache(d::ConstDrift, b::ConstAsymBounds) = AsymPDFConstCache(d, b)
PDFConstCache(d::ConstDrift, b::ConstSymBounds) = SymPDFConstCache(d, b)

pdfu(d::ConstDrift, b::ConstBounds, t::Real, tol::Real=1.e-29) = pdfu(PDFConstCache(d, b), t, tol)
pdfl(d::ConstDrift, b::ConstBounds, t::Real, tol::Real=1.e-29) = pdfl(PDFConstCache(d, b), t, tol)
pdful(d::ConstDrift, b::ConstBounds, t::Real, tol::Real=1.e-29) = pdful(PDFConstCache(d, b), t, tol)

function pdf(d::ConstDrift, b::ConstBounds, tmax::Real, tol::Real=1.e-29)
    tmax >= zero(tmax) || error("tmax needs to be non-negative")
    c = PDFConstCache(d, b)
    dt = getdt(d)
    @assert getdt(b) == dt
    maxn = length(0:dt:tmax)

    gu, gl, t = Array{Float64}(undef, maxn), Array{Float64}(undef, maxn), 0.0
    for n = 1:maxn
        gu[n], gl[n] = pdful(c, t, tol)
        t += dt
    end
    return gu, gl
end

# time-varying bound, time-varying drift
#
# using Volterra series expansion from
#     PL Smith (2000). Stochastic dynamic models of response time and
#     accuracy: a foundational primer. Journal of Mathematical
#     psychology 44, 408-463.

function pdf(d::AbstractDrift, b::AbstractBounds, tmax::Real)
    tmax >= zero(tmax) || error("tmax needs to be non-negative")

    dt = getdt(d)
    @assert getdt(b) == dt
    maxn = length(0:dt:tmax)
    @assert getmaxn(d) >= maxn && getmaxn(b) >= maxn

    g1, g2 = Array{Float64}(undef, maxn), Array{Float64}(undef, maxn)
    g1[1], g2[1] = 0.0, 0.0
    if maxn == 1
        return g1, g2
    end

    cΨ = 1.0 / √(twoπ)
    # computes diffusion kernel 2Ψ(b(t), t | y, tau) for parameters
    # mut = drift(t), bgradt = d b(t) / dt, dbm = b(t) - y - m(t) + m(tau),
    # dv = v(t) - v(tau), sqrtdv = √(dv)
    _Ψ(mut, bgradt, dbm, dv, sqrtdv) =
        cΨ * exp(-dbm * dbm / (2dv)) / sqrtdv * (bgradt - mut - dbm / dv)

    g1[2] = - _Ψ(getmu(d, 2), getuboundgrad(b, 2), getubound(b, 2) - getm(d, 2), dt, √(dt))
    g2[2] = _Ψ(getmu(d, 2), getlboundgrad(b, 2), getlbound(b, 2) - getm(d, 2), dt, √(dt))

    for n = 3:maxn
        g1n, g2n = 0.0, 0.0
        mun, mn = getmu(d, n), getm(d, n)
        bupn, blon = getubound(b, n), getlbound(b, n)
        bupgradn, blogradn = getuboundgrad(b, n), getlboundgrad(b, n)
        for j = 2:(n-1)
            dv = dt * (n-j)
            sqrtdv = √(dv)
            dbmu = getm(d, j) - mn - getubound(b, j)
            dbml = getm(d, j) - mn - getlbound(b, j)
            g1n += g1[j] * _Ψ(mun, bupgradn, bupn + dbmu, dv, sqrtdv) +
                g2[j] * _Ψ(mun, bupgradn, bupn + dbml, dv, sqrtdv)
            g2n += g1[j] * _Ψ(mun, blogradn, blon + dbmu, dv, sqrtdv) +
                g2[j] * _Ψ(mun, blogradn, blon + dbml, dv, sqrtdv)
        end
        g1[n] = - _Ψ(mun, bupgradn, bupn - mn, dt * (n-1), √(dt * (n-1))) + dt * g1n
        g2[n] = _Ψ(mun, blogradn, blon - mn, dt * (n-1), √(dt * (n-1))) - dt * g2n
    end
    return g1, g2
end


text:
mutable struct State
  ptype :: Int
  class :: Int
end

# outer constructors

State() = State(1, 1)
State(ptype) = State(ptype, 1)

get_ptype(x :: State) = x.ptype
get_neighbor_class(x :: State) = x.class

transform!(x :: State, l) = (x.ptype = l)
change_neighbor_class!(x :: State, k) = (x.class = k)

function Base.show(io :: IO, x :: State)
  println(io, "  type: ", x.ptype)
  print(io, "  nbclass: ", x.class)
end

struct Site{D,T}
  label :: Int
  state :: State
  coord :: MVector{D,T}

  function Site{D,T}(id :: Integer, state :: State, coord :: NTuple{D,T}) where {D,T}
      new{D,T}(id, state, MVector(coord...))
  end
end

Site(id, state, coord::NTuple{N,T}) where {N,T} = Site{N,T}(id, state, coord)

# accessors #

label(x :: Site) = x.label
state(x :: Site) = x.state
coordinates(x :: Site) = x.coord.data

# for search / sort compatibility
label(x::Int) = x
coordinates(x::NTuple) = x
coordinates(x::AbstractVector) = x

get_ptype(x :: Site) = get_ptype(state(x))
get_neighbor_class(x :: Site) = get_neighbor_class(state(x))

# state #

transform!(x :: Site, l) = transform!(state(x), l)
change_neighbor_class!(x :: Site, k) = change_neighbor_class!(state(x), k)

# iteration #
Base.iterate(x::Site) = iterate(coordinates(x))
Base.iterate(x::Site, state) = iterate(coordinates(x), state)

# indexing #

Base.getindex(x :: Site, i) = getindex(coordinates(x), i)
Base.setindex!(x :: Site, val, i) = setindex!(coordinates(x), val, i)
Base.eachindex(x :: Site) = eachindex(coordinates(x))

# IO #

Base.summary(io::IO, x::Site) = print(io, "Site")

function Base.show(io::IO, x::Site)
  # summary(io, x)
  print(io, coordinates(x))
end

function Base.show(io::IO, m::MIME"text/plain", x::Site)
  lk = string(get_ptype(x), "|", get_neighbor_class(x))
  print(io, lk)
  print(io, " @ ")
  print(io, coordinates(x))
end

##### recipes #####

@recipe function f(::Type{Site{1,T}}, x::Site{1,T}) where T
  [(0, x[1])]
end

@recipe function f(::Type{Site{2,T}}, x::Site{2,T}) where T
  [(x[1], x[2])]
end

@recipe function f(::Type{Site{3,T}}, x::Site{3,T}) where T
  [(x[1], x[2], x[3])]
end

function build_composition!(comp, lattice, x::Site{D}) where D
  fill!(comp, zero(eltype(comp)))
  comp[1] = capacity(topology(lattice), D)

  for id in neighborhood(lattice, x)
    y = get_site(lattice, id)
    l = get_ptype(y)

    comp[1] -= 1
    comp[l] += 1
  end

  return comp
end

##### sorting #####

const XD1 = (1,)

const XD2 = (1, 2)
const YD2 = (2, 1)

const XD3 = (1, 2, 3)
const YD3 = (1, 3, 2)
const ZD3 = (2, 3, 1)

Base.isless(x::Site, y::Site) = isless(label(x), label(y))

Base.:(<)(::Site{D}, ::Tuple{}) where D = false
Base.:(<)(x::Site{D}, y::NTuple{D}) where D = x.coord.data < y

# does not work for missing but should not be an issue
function lex_order(x, y, idx)
  for i in idx
    a, b = x[i], y[i]
    if !isequal(a, b)
      return a < b
    end
  end
  return false
end

lex_order_x1(x, y) = lex_order(x, y, XD1)

lex_order_x2(x, y) = lex_order(x, y, XD2)
lex_order_y2(x, y) = lex_order(x, y, YD2)

lex_order_x3(x, y) = lex_order(x, y, XD3)
lex_order_y3(x, y) = lex_order(x, y, YD3)
lex_order_z3(x, y) = lex_order(x, y, ZD3)

sort_by_x_1D!(sites) = sort!(sites,
  alg = Base.Sort.DEFAULT_STABLE,
  lt = lex_order_x1,
  by = coordinates)

sort_by_x_2D!(sites) = sort!(sites,
  alg = Base.Sort.DEFAULT_STABLE,
  lt = lex_order_x2,
  by = coordinates)

sort_by_y_2D!(sites) = sort!(sites,
  alg = Base.Sort.DEFAULT_STABLE,
  lt = lex_order_y2,
  by = coordinates)

sort_by_x_3D!(sites) = sort!(sites,
  alg = Base.Sort.DEFAULT_STABLE,
  lt = lex_order_x3,
  by = coordinates)

sort_by_y_3D!(sites) = sort!(sites,
  alg = Base.Sort.DEFAULT_STABLE,
  lt = lex_order_y3,
  by = coordinates)

sort_by_z_3D!(sites) = sort!(sites,
  alg = Base.Sort.DEFAULT_STABLE,
  lt = lex_order_z3,
  by = coordinates)


text:
""" Limits and colimits in a category.
"""
module Limits
export AbstractLimit, AbstractColimit, Limit, Colimit,
  ob, cone, cocone, apex, legs, limit, colimit, universal,
  Terminal, Initial, terminal, initial, delete, create, factorize,
  BinaryProduct, Product, product, proj1, proj2, pair,
  BinaryPullback, Pullback, BinaryEqualizer, Equalizer, pullback, incl,
  BinaryCoproduct, Coproduct, coproduct, coproj1, coproj2, copair,
  BinaryPushout, Pushout, BinaryCoequalizer, Coequalizer, pushout, proj,
  composite_pullback, composite_pushout

using Compat: only

using AutoHashEquals

using ...Theories
import ...Theories: ob, terminal, product, proj1, proj2, equalizer, incl,
  initial, coproduct, coproj1, coproj2, coequalizer, proj,
  delete, create, pair, copair, factorize
using ..FreeDiagrams
import ..FreeDiagrams: apex, legs

# Data types for limits
#######################

""" Abstract type for limit in a category.

The standard concrete subtype is [`Limit`](@ref), although for computational
reasons certain categories may use different subtypes to include extra data.
"""
abstract type AbstractLimit{Ob,Diagram} end

ob(lim::AbstractLimit) = apex(lim)
cone(lim::AbstractLimit) = lim.cone
apex(lim::AbstractLimit) = apex(cone(lim))
legs(lim::AbstractLimit) = legs(cone(lim))

Base.iterate(lim::AbstractLimit, args...) = iterate(cone(lim), args...)
Base.eltype(lim::AbstractLimit) = eltype(cone(lim))
Base.length(lim::AbstractLimit) = length(cone(lim))

""" Limit in a category.
"""
@auto_hash_equals struct Limit{Ob,Diagram,Cone<:Multispan{Ob}} <:
    AbstractLimit{Ob,Diagram}
  diagram::Diagram
  cone::Cone
end

const Terminal{Ob} = AbstractLimit{Ob,<:EmptyDiagram}
const BinaryProduct{Ob} = AbstractLimit{Ob,<:ObjectPair}
const Product{Ob} = AbstractLimit{Ob,<:DiscreteDiagram}
const BinaryPullback{Ob} = AbstractLimit{Ob,<:Cospan}
const Pullback{Ob} = AbstractLimit{Ob,<:Multicospan}
const BinaryEqualizer{Ob} = AbstractLimit{Ob,<:ParallelPair}
const Equalizer{Ob} = AbstractLimit{Ob,<:ParallelMorphisms}

proj1(lim::Union{BinaryProduct,BinaryPullback}) = first(legs(lim))
proj2(lim::Union{BinaryProduct,BinaryPullback}) = last(legs(lim))
incl(eq::Equalizer) = only(legs(eq))

# Data types for colimits
#########################

""" Abstract type for colimit in a category.

The standard concrete subtype is [`Colimit`](@ref), although for computational
reasons certain categories may use different subtypes to include extra data.
"""
abstract type AbstractColimit{Ob,Diagram} end

ob(colim::AbstractColimit) = apex(colim)
cocone(colim::AbstractColimit) = colim.cocone
apex(colim::AbstractColimit) = apex(cocone(colim))
legs(colim::AbstractColimit) = legs(cocone(colim))

Base.iterate(colim::AbstractColimit, args...) = iterate(cocone(colim), args...)
Base.eltype(colim::AbstractColimit) = eltype(cocone(colim))
Base.length(colim::AbstractColimit) = length(cocone(colim))

""" Colimit in a category.
"""
@auto_hash_equals struct Colimit{Ob,Diagram,Cocone<:Multicospan{Ob}} <:
    AbstractColimit{Ob,Diagram}
  diagram::Diagram
  cocone::Cocone
end

const Initial{Ob} = AbstractColimit{Ob,<:EmptyDiagram}
const BinaryCoproduct{Ob} = AbstractColimit{Ob,<:ObjectPair}
const Coproduct{Ob} = AbstractColimit{Ob,<:DiscreteDiagram}
const BinaryPushout{Ob} = AbstractColimit{Ob,<:Span}
const Pushout{Ob} = AbstractColimit{Ob,<:Multispan}
const BinaryCoequalizer{Ob} = AbstractColimit{Ob,<:ParallelPair}
const Coequalizer{Ob} = AbstractColimit{Ob,<:ParallelMorphisms}

coproj1(colim::Union{BinaryCoproduct,BinaryPushout}) = first(legs(colim))
coproj2(colim::Union{BinaryCoproduct,BinaryPushout}) = last(legs(colim))
proj(coeq::Coequalizer) = only(legs(coeq))

# Generic (co)limits
####################

""" Limit of a diagram.

To define limits in a category with objects `Ob`, override the method
`limit(::FreeDiagram{Ob})` for general limits or `limit(::D)` with suitable type
`D <: FixedShapeFreeDiagram{Ob}` for limits of specific shape, such as products
or equalizers.

See also: [`colimit`](@ref)
"""
function limit end

""" Colimit of a diagram.

To define colimits in a category with objects `Ob`, override the method
`colimit(::FreeDiagram{Ob})` for general colimits or `colimit(::D)` with
suitable type `D <: FixedShapeFreeDiagram{Ob}` for colimits of specific shape,
such as coproducts or coequalizers.

See also: [`limit`](@ref)
"""
function colimit end

""" Universal property of (co)limits.

Compute the morphism whose existence and uniqueness is guaranteed by the
universal property of (co)limits.

See also: [`limit`](@ref), [`colimit`](@ref).
"""
function universal end

# Specific (co)limits
#####################

""" Terminal object.

To implement for a type `T`, define the method `limit(::EmptyDiagram{T})`.
"""
terminal(T::Type) = limit(EmptyDiagram{T}())

""" Initial object.

To implement for a type `T`, define the method `colimit(::EmptyDiagram{T})`.
"""
initial(T::Type) = colimit(EmptyDiagram{T}())

""" Unique morphism into a terminal object.

To implement for a type `T`, define the method
`universal(::Terminal{T}, ::SMultispan{0,T})`.
"""
delete(A::T) where T = delete(terminal(T), A)
delete(lim::Terminal, A) = universal(lim, SMultispan{0}(A))

""" Unique morphism out of an initial object.

To implement for a type `T`, define the method
`universal(::Initial{T}, ::SMulticospan{0,T})`.
"""
create(A::T) where T = create(initial(T), A)
create(colim::Initial, A) = universal(colim, SMulticospan{0}(A))

""" Product of objects.

To implement for a type `T`, define the method `limit(::ObjectPair{T})` and/or
`limit(::DiscreteDiagram{T})`.
"""
product(A, B) = limit(ObjectPair(A, B))
product(As::AbstractVector) = limit(DiscreteDiagram(As))

""" Coproduct of objects.

To implement for a type `T`, define the method `colimit(::ObjectPair{T})` and/or
`colimit(::DiscreteDiagram{T})`.
"""
coproduct(A, B) = colimit(ObjectPair(A, B))
coproduct(As::AbstractVector) = colimit(DiscreteDiagram(As))

""" Equalizer of morphisms with common domain and codomain.

To implement for a type `T`, define the method `limit(::ParallelPair{T})` and/or
`limit(::ParallelMorphisms{T})`.
"""
equalizer(f, g) = limit(ParallelPair(f, g))
equalizer(fs::AbstractVector) = limit(ParallelMorphisms(fs))

""" Coequalizer of morphisms with common domain and codomain.

To implement for a type `T`, define the method `colimit(::ParallelPair{T})` or
`colimit(::ParallelMorphisms{T})`.
"""
coequalizer(f, g) = colimit(ParallelPair(f, g))
coequalizer(fs::AbstractVector) = colimit(ParallelMorphisms(fs))

""" Pullback of a pair of morphisms with common codomain.

To implement for a type `T`, define the method `limit(::Cospan{T})` and/or
`limit(::Multicospan{T})` or, if you have already implemented products and
equalizers, rely on the default implementation.
"""
pullback(f, g) = limit(Cospan(f, g))
pullback(fs::AbstractVector) = limit(Multicospan(fs))

""" Pushout of a pair of morphisms with common domain.

To implement for a type `T`, define the method `colimit(::Span{T})` and/or
`colimit(::Multispan{T})` or, if you have already implemented coproducts and
coequalizers, rely on the default implementation.
"""
pushout(f, g) = colimit(Span(f, g))
pushout(fs::AbstractVector) = colimit(Multispan(fs))

""" Pairing of morphisms: universal property of products/pullbacks.

To implement for products of type `T`, define the method
`universal(::BinaryProduct{T}, ::Span{T})` and/or
`universal(::Product{T}, ::Multispan{T})` and similarly for pullbacks.
"""
pair(f, g) = pair(product(codom(f), codom(g)), f, g)
pair(fs::AbstractVector) = pair(product(map(codom, fs)), fs)
pair(lim::Union{BinaryProduct,BinaryPullback}, f, g) =
  universal(lim, Span(f, g))
pair(lim::Union{Product,Pullback}, fs::AbstractVector) =
  universal(lim, Multispan(fs))

""" Copairing of morphisms: universal property of coproducts/pushouts.

To implement for coproducts of type `T`, define the method
`universal(::BinaryCoproduct{T}, ::Cospan{T})` and/or
`universal(::Coproduct{T}, ::Multicospan{T})` and similarly for pushouts.
"""
copair(f, g) = copair(coproduct(dom(f), dom(g)), f, g)
copair(fs::AbstractVector) = copair(coproduct(map(dom, fs)), fs)
copair(colim::Union{BinaryCoproduct,BinaryPushout}, f, g) =
  universal(colim, Cospan(f, g))
copair(colim::Union{Coproduct,Pushout}, fs::AbstractVector) =
  universal(colim, Multicospan(fs))

""" Factor morphism through (co)equalizer, via the universal property.

To implement for equalizers of type `T`, define the method
`universal(::Equalizer{T}, ::SMultispan{1,T})`. For coequalizers of type `T`,
define the method `universal(::Coequalizer{T}, ::SMulticospan{1,T})`.
"""
factorize(lim::Equalizer, h) = universal(lim, SMultispan{1}(h))
factorize(colim::Coequalizer, h) = universal(colim, SMulticospan{1}(h))

# Composite (co)limits
######################

""" Pullback formed as composite of product and equalizer.

The fields of this struct are an implementation detail; accessing them directly
violates the abstraction. Everything that you can do with a pushout, including
invoking its universal property, should be done through the generic interface
for limits.

See also: [`CompositePushout`](@ref).
"""
struct CompositePullback{Ob, Diagram<:Multicospan{Ob}, Cone<:Multispan{Ob},
    Prod<:Product{Ob}, Eq<:Equalizer{Ob}} <: AbstractLimit{Ob,Diagram}
  diagram::Diagram
  cone::Cone
  prod::Prod
  eq::Eq
end

""" Compute pullback as composite of product and equalizer.
"""
function composite_pullback(cospan::Cospan)
  f, g = cospan
  (π1, π2) = prod = product(dom(f), dom(g))
  (ι,) = eq = equalizer(π1⋅f, π2⋅g)
  CompositePullback(cospan, Span(ι⋅π1, ι⋅π2), prod, eq)
end
composite_pullback(f, g) = composite_pullback(Cospan(f, g))

function universal(lim::CompositePullback, cone::Multispan)
  factorize(lim.eq, universal(lim.prod, cone))
end

""" Pushout formed as composite of coproduct and equalizer.

See also: [`CompositePullback`](@ref).
"""
struct CompositePushout{Ob, Diagram<:Multispan{Ob}, Cocone<:Multicospan{Ob},
    Coprod<:Coproduct{Ob}, Coeq<:Coequalizer{Ob}} <: AbstractColimit{Ob,Diagram}
  diagram::Diagram
  cocone::Cocone
  coprod::Coprod
  coeq::Coeq
end

""" Compute pushout as composite of coproduct and coequalizer.
"""
function composite_pushout(span::Span)
  f, g = span
  (ι1, ι2) = coprod = coproduct(codom(f), codom(g))
  (π,) = coeq = coequalizer(f⋅ι1, g⋅ι2)
  CompositePushout(span, Cospan(ι1⋅π, ι2⋅π), coprod, coeq)
end
composite_pushout(f, g) = composite_pushout(Span(f, g))

function universal(lim::CompositePushout, cone::Multicospan)
  factorize(lim.coeq, universal(lim.coprod, cone))
end

end


text:
using SEAL
using Test

# Include files with example-specific tests
include("test_1_bfv_basics.jl")
include("test_2_encoders.jl")
include("test_3_levels.jl")
include("test_4_ckks_basics.jl")
include("test_5_rotation.jl")
include("test_6_serialization.jl")
include("test_extra.jl")


text:
struct QGramIterator{S <: Union{AbstractString, AbstractVector}}
	s::S   # Collection
	q::Int # Length of Qgram
end
Base.length(qgram::QGramIterator) = max(length(qgram.s) - qgram.q + 1, 0)

# q-grams of AbstractString
function Base.iterate(qgram::QGramIterator{<: AbstractString},
	state = (1, nextind(qgram.s, 0, qgram.q)))
	istart, iend = state
	iend > ncodeunits(qgram.s) && return nothing
	element = SubString(qgram.s, istart, iend)
	nextstate = nextind(qgram.s, istart), nextind(qgram.s, iend)
	element, nextstate
end
Base.eltype(qgram::QGramIterator{SubString{S}}) where {S} = SubString{S}
Base.eltype(qgram::QGramIterator{S}) where {S <: AbstractString} = SubString{S}

#q-grams of AbstractVector
# Alternatively, I could also use partition in IterTools but it creates a vector for each iteration
# so it does not seem to be worth it.
function Base.iterate(qgram::QGramIterator{<: AbstractVector}, state = firstindex(qgram.s))
	state + qgram.q - 1 > lastindex(qgram.s) && return nothing
	view(qgram.s, state:(state + qgram.q - 1)), state + 1
end
Base.eltype(qgram::QGramIterator{<: AbstractVector}) = typeof(first(qgram))

"""
Return an iterator on the q-gram of a string

### Arguments
* `s` iterator
* `q::Integer`: length of q-gram

## Examples
```julia
for x in qgrams("hello", 2)
	println(x)
end
```
"""
qgrams(s::Union{AbstractString, AbstractVector}, q::Integer) = QGramIterator(s, q)
qgrams(s, q::Integer) = QGramIterator(collect(s), q)

abstract type QGramDistance <: SemiMetric end

# For two iterators x1 and x2, that define a length and eltype method,
# this returns a dictionary which, for each element in x1 or x2,
# returns a tuple with the numbers of times it appears in x1 and x2
function count_map(s1, s2)
	K = promote_type(eltype(s1), eltype(s2))
	d = Dict{K, Tuple{Int, Int}}()
	sizehint!(d, length(s1) + length(s2))
	# I use a faster way to change a dictionary key
	# see setindex! in https://github.com/JuliaLang/julia/blob/master/base/dict.jl#L380
	for x1 in s1
		index = Base.ht_keyindex2!(d, x1)
		if index > 0
			d.age += 1
			@inbounds d.keys[index] = x1
			@inbounds d.vals[index] = (d.vals[index][1] + 1, 0)
		else
			@inbounds Base._setindex!(d, (1, 0), x1, -index)
		end
	end
	for x2 in s2
		index = Base.ht_keyindex2!(d, x2)
		if index > 0
			d.age += 1
			@inbounds d.keys[index] = x2
			@inbounds d.vals[index] = (d.vals[index][1], d.vals[index][2] + 1)
		else
			@inbounds Base._setindex!(d, (0, 1), x2, -index)
		end
	end
	return d
end

"""
	QGram(q::Int)

Creates a QGram metric.

The distance corresponds to

``||v(s1, q) - v(s2, q)||``

where ``v(s, q)`` denotes the vector on the space of q-grams of length q,
that contains the number of times a q-gram appears for the string s
"""
struct QGram <: QGramDistance
	q::Int
end

function evaluate(dist::QGram, s1, s2)
	(ismissing(s1) | ismissing(s2)) && return missing
	itr = values(count_map(qgrams(s1, dist.q), qgrams(s2, dist.q)))
	n = 0
	itr =
	for (n1, n2) in itr
		n += abs(n1 - n2)
	end
	n
end

"""
	Cosine(q::Int)

Creates a Cosine metric.

The distance corresponds to

`` 1 - v(s1, q).v(s2, q)  / ||v(s1, q)|| * ||v(s2, q)||``

where ``v(s, q)`` denotes the vector on the space of q-grams of length q,
that contains the  number of times a q-gram appears for the string s
"""
struct Cosine <: QGramDistance
	q::Int
end

function evaluate(dist::Cosine, s1, s2, max_dist = nothing)
	(ismissing(s1) | ismissing(s2)) && return missing
	itr = values(count_map(qgrams(s1, dist.q), qgrams(s2, dist.q)))
	norm1, norm2, prodnorm = 0, 0, 0
	for (n1, n2) in itr
		norm1 += n1^2
		norm2 += n2^2
		prodnorm += n1 * n2
	end
	1.0 - prodnorm / (sqrt(norm1) * sqrt(norm2))
end

"""
	Jaccard(q::Int)

Creates a Jaccard metric.

The distance corresponds to

``1 - |Q(s1, q) ∩ Q(s2, q)| / |Q(s1, q) ∪ Q(s2, q))|``

where ``Q(s, q)``  denotes the set of q-grams of length n for the string s
"""
struct Jaccard <: QGramDistance
	q::Int
end

function evaluate(dist::Jaccard, s1, s2, max_dist = nothing)
	(ismissing(s1) | ismissing(s2)) && return missing
	itr = values(count_map(qgrams(s1, dist.q), qgrams(s2, dist.q)))
	ndistinct1, ndistinct2, nintersect = 0, 0, 0
	for (n1, n2) in itr
		ndistinct1 += n1 > 0
		ndistinct2 += n2 > 0
		nintersect += (n1 > 0) & (n2 > 0)
	end
	1.0 - nintersect / (ndistinct1 + ndistinct2 - nintersect)
end

"""
	SorensenDice(q::Int)

Creates a SorensenDice metric

The distance corresponds to

``1 - 2 * |Q(s1, q) ∩ Q(s2, q)|  / (|Q(s1, q)| + |Q(s2, q))|)``

where ``Q(s, q)``  denotes the set of q-grams of length n for the string s
"""
struct SorensenDice <: QGramDistance
	q::Int
end

function evaluate(dist::SorensenDice, s1, s2, max_dist = nothing)
	(ismissing(s1) | ismissing(s2)) && return missing
	itr = values(count_map(qgrams(s1, dist.q), qgrams(s2, dist.q)))
	ndistinct1, ndistinct2, nintersect = 0, 0, 0
	for (n1, n2) in itr
		ndistinct1 += n1 > 0
		ndistinct2 += n2 > 0
		nintersect += (n1 > 0) & (n2 > 0)
	end
	1.0 - 2.0 * nintersect / (ndistinct1 + ndistinct2)
end

"""
	Overlap(q::Int)

Creates a Overlap metric

The distance corresponds to

``1 - |Q(s1, q) ∩ Q(s2, q)|  / min(|Q(s1, q)|, |Q(s2, q)|)``

where ``Q(s, q)``  denotes the set of q-grams of length n for the string s
"""
struct Overlap <: QGramDistance
	q::Int
end

function evaluate(dist::Overlap, s1, s2, max_dist = nothing)
	(ismissing(s1) | ismissing(s2)) && return missing
	itr = values(count_map(qgrams(s1, dist.q), qgrams(s2, dist.q)))
	ndistinct1, ndistinct2, nintersect = 0, 0, 0
	for (n1, n2) in itr
		ndistinct1 += n1 > 0
		ndistinct2 += n2 > 0
		nintersect += (n1 > 0) & (n2 > 0)
	end
	1.0 - nintersect / min(ndistinct1, ndistinct2)
end


text:
using Documenter, SplittablesBase, SplittablesTesting

makedocs(;
    modules=[SplittablesBase, SplittablesTesting],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/JuliaFolds/SplittablesBase.jl/blob/{commit}{path}#L{line}",
    sitename="SplittablesBase.jl",
    authors="Takafumi Arakaki <aka.tkf@gmail.com>",
)

deploydocs(;
    repo="github.com/JuliaFolds/SplittablesBase.jl",
    push_preview=true,
)


text:
using Clustering
using LinearAlgebra
using Plots
using Plots.PlotMeasures

"""
`plot_pmdi_data(dataFile,
                psm::Posterior_similarity_matrix;
                k::Union{Int64, Nothing} = nothing,
                h::Union{Float64, Nothing} = nothing,
                orderby::Int64 = 0,
                featureSelectProbs = nothing,
                z_Score::Bool = false)`

Generates a plot of a single raw data file, with the inferred clusters (as derived from the psm). If featureSelectProbs is supplied, features will be ordered according to the frequency they are selected in the model and the density of feature selection probabilities will be displayed along the margins.
## Input
- `psm` a posterior similarity matrix struct as output from `generate_psm`
- `k` optional - the number of clusters to highlight in the data
- `h` optional - the distance at which to cut the dendrogram. One of `k` or `h` must be specified.
- `orderby` observations in all plots are ordered in a way to separate clusters as well as possible. `orderby` specifies which dataset should be used to inform this ordering. `orderby = 0` will let the overall consensus dictate the ordering.
- `featureSelect` optional - a Vector of feature selection probabilities as output from `get_feature_select_probs()`.
- `z_score` - a Bool indicating whether the data should be standardised and discretised. Can help elucidate patterns.

## Output
- A heatmap of the data with clusters indicated.
"""
function plot_pmdi_data(data,
                        psm::Posterior_similarity_matrix;
                        k::Union{Int64, Nothing} = nothing,
                        h::Union{Float64, Nothing} = nothing,
                        orderby::Int64 = 0,
                        featureSelectProbs = nothing,
                        z_score::Bool = false,
                        linkage = :ward)
  dataFile = deepcopy(data)
  if featureSelectProbs != nothing
    @assert length(featureSelectProbs) == size(dataFile, 2) "Feature selection vector is not the same length as the number of features"
  end
  if (z_score) & (typeof(dataFile) == Array{Float64, 2})
    sd = std(dataFile, dims = 1)
    dataFile .-= mean(dataFile, dims = 1)
    dataFile ./= sd
    dataFile .= floor.(Int64, dataFile)
    dataFile[dataFile .> 2] .= 2
    dataFile[dataFile .< - 2] .= - 2
  end

  # Extract clusters
  if orderby == 0
      orderby = size(psm.psm, 1)
  end
  hc = hclust(1 .- Symmetric(psm.psm[orderby], :L), linkage = linkage)
  if h == nothing
    cuts = cutree(hc, k = k)[hc.order]
  elseif k == nothing
    cuts = cutree(hc, h = h)[hc.order]
  end

  nclust = length(unique(cuts))
  ticks = indexin(1:nclust, cuts) .- 0.5
  sort!(ticks)
  ticks = ticks[2:end]
  order_rows = hc.order

  if featureSelectProbs != nothing
    order_cols = sortperm(featureSelectProbs, rev = true)
    fs_plot = Plots.bar(featureSelectProbs[order_cols] .* size(data, 1),
                        legend = false,
                        fill = "#000000",
                        linealpha = 0,
                        yflip = true,
                        tickfontsize = 1,
                        xticks = false,
                        grid = false,
                        widen = false,
                        top_margin = - 7.5px,
                        left_margin = 0px,
                        right_margin = - 7.5px,
                        xlims = (0.5, size(dataFile, 2) + 0.5),
                        xlabel = "Features",
                        ylabel = "P(select)")
    l = @layout [a{0.97w, 0.8h} b{0.03w, 0.8h}; c{0.97w, 0.2h} d{0.03w, 0.2h}]
  else
    l = @layout [a{0.97w, 0.97h} b{0.03w, 0.97h};]
    order_cols = 1:size(dataFile, 2)
  end

  # clusters_matrix = Matrix{Int64}(undef, size(dataFile, 1), 1)
  # clusters_matrix .= cuts

  clust_mat = Matrix{Int64}(undef, nclust, 2)
  clust_mat[:, 1] =  unique(cuts)
  tmp = cumsum(map(x -> count(cuts .== x), clust_mat[:, 1]))
  clust_mat[:, 2] = tmp[end:-1:1]
  labs_y = vcat(tmp[1] / 2,
              tmp[1:end - 1] + (tmp[2:end] - tmp[1:end-1]) ./ 2)
  # C(n) = [RGB(cgrad(:viridis)[floor(Int, z)]) for z in collect(LinRange(1, 30, n))]
  # clusters_plot = plot(palette = C(nclust))
  clusters_plot = plot()
  bar!(clust_mat[:, 2]',
          seriescolor = clust_mat[end:-1:1, 1]',
          color = "#C8C8C8",
          legend = false,
          yflip = false,
          widen = false,
          yaxis = false,
          xticks = false,
          tickfontsize = 1,
          bottom_margin = -7.5px,
          left_margin = - 7.5px,
          annotation = (1, labs_y, clust_mat[:, 1]))

  #clusters_plot = Plots.heatmap(clusters_matrix,
  #tickfontsize = 1,
  #legend = false,
  #yaxis = false,
  #xticks = false,
  #palette = :viridis,
  #bottom_margin = -7.5px,
  #left_margin = - 7.5px)
  #hline!(ticks, palette = "#000000")

  clusters_plot2 = Plots.bar(clust_mat[:, 2]',
  tickfontsize = 1,
  legend = false,
  widen = false,
  α = 0,
  top_margin = - 7.5px,
  left_margin = - 7.5px
  )

  data_plot = Plots.heatmap(dataFile[order_rows, order_cols],
  tickfontsize = 1,
  xticks = false,
  yflip = false,
  legend = false,
  widen = false,
  left_margin = 0px,
  bottom_margin = -7.5px,
  right_margin = -7.5px,
  ylabel = "Observations",
  c = :viridis,
  xlims = (0.5, size(dataFile, 2) + 0.5)
  )
  hline!(ticks, linestyle = :dash, c = "#FFFFFF", linewidth = 1)

  if featureSelectProbs != nothing
    out = [data_plot, clusters_plot, fs_plot, clusters_plot2]
  else
    xlabel!("Features")
    out = [data_plot, clusters_plot]
  end

    Plots.plot(out...,
    layout = l,
    link = :both,
    framestyle = :none,
    right_margin = -5px,
    tickfontsize = 1)
end

"""
`get_feature_select_probs(featureSelect::String, dataName::String)`

Generates a vector of mean feature selection probabilities from the feature selection output of `pmdi()`.
## Input
- `featureSelect::String` - a file specifying a .csv file containing feature selection flags as returned from `pmdi()`
- `burnin::Int` - an integer specifying the number of iterations to discard from the beginning of the output
- `thin::Int` - an integer specifying the rate of thinning of the output, e.g. `thin = 2` will keep only every second iteration.
## Output
- A Vector of length `K`, where `K` is the number of datasets analysed, each element containing a feature selection probability for each feature.
"""
function get_feature_select_probs(featureSelect::String, burnin::Int = 0, thin::Int = 1)
  outputNames = split(readline(featureSelect), ',')
  dataNames = unique([replace(outputName, r"([A-Za-z0-9])(_d.+)" => s"\1") for outputName in outputNames])

  out = Vector{Any}(undef, length(dataNames))
  for i in 1:length(out)
      dataInd = map(outputNames) do str
          occursin(dataNames[i], str)
      end
      out[i] = [mean(readdlm(featureSelect, ',')[(2 + burnin):thin:end, dataInd], dims = 1)[:]][1]
  end
  return out
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module StateSelection

DAE-model to test manual state selection of simple multibody system.

This model is mainly used to test initialization
"""
module StateSelection

import ..ModiaMath
using  ..LinearAlgebra

mutable struct Model <: ModiaMath.AbstractSimulationModel
    simulationState::ModiaMath.SimulationState

    # Parameters
    n::Vector{Float64}
    m::Float64
    g::Float64
    s0::Float64
    v0::Float64

    function Model(;n=[0.0, 1.0, 0.0], m=1.0, g=9.81, s0=0.0, v0=0.0)
        @assert(length(n) == 3)
        @assert(m > 0.0)
        simulationState = ModiaMath.SimulationState("StateSelection", getModelResidues!, [s0,0.0,0.0,0.0,v0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], getVariableName;
                                x_fixed=[true, false,false,false, true, false,false,false, false, false,false,false],
								structureOfDAE = ModiaMath.DAE_NoSpecialStructure,
                                hev=1e-4)
        new(simulationState, n, m, g, s0, v0)
    end
end

getVariableName(model, vcat, vindex) = ModiaMath.getVariableName(model, vcat, vindex;
                                                xNames=["s", "f[1]", "f[2]", "f[3]", "sd", "der_der_r[1]","der_der_r[2]","der_der_r[3]",
                                                             "der_der_s", "der_v[1]", "der_v[2]", "der_v[3]"])

function getModelResidues!(m::Model, t::Float64, _x::Vector{Float64}, _derx::Vector{Float64}, _r::Vector{Float64}, _w::Vector{Float64})
    # map x
    s  = _x[1]
    f  = _x[2:4]
	g  = [0.0, 0.0, m.g]
    sd = _x[5]
    der_der_r = _x[6:8]
    der_der_s = _x[9]
    der_v     = _x[10:12]

    # map derx
    ders  = _derx[1]
    dersd = _derx[5]

    # compute residues
    _r[1]    = sd - ders
    _r[2:4]  = m.n * der_der_s - der_der_r
    _r[5:7]  = der_der_r - der_v
    _r[8:10] = f - m.m*g - m.m * der_v
    _r[11]   = dot(m.n, f)
    _r[12]   = der_der_s - dersd
    return nothing
end

end


text:
###
  Name     : c6_19.jl
  Book     : Hands-on Data Science with Anaconda )
  Publisher: Packt Publishing Ltd.
  Author   : Yuxing Yan and James Yan
  Date     : 1/25/2018
  email    : yany@canisius.edu
             paulyxy@hotmail.com
###

Pkg.add("Gadfly")
using Gadfly
draw(SVG("output.svg", 6inch, 3inch), plot([sin, cos], 0, 25))

function module_functions(modname)
    list = Symbol[]
    for nm in names(modname)
        typeof(eval(nm)) == Function && push!(list,nm)
    end
    return list
end

using PyPlot
module_functions(PyPlot)

using MTH229
f(x) = exp(-2*x^2)
plot(f, -3, 3)

using QuantEcon: meshgrid
n = 50
x = linspace(-3, 3, n)
y = x
z = Array{Float64}(n, n)
f(x, y) = cos(x^2 + y^2) / (1 + x^2 + y^2)
for i in 1:n
    for j in 1:n
        z[j, i] = f(x[i], y[j])
    end
end
xgrid, ygrid = meshgrid(x, y)
surf(xgrid, ygrid, z', cmap=ColorMap("jet"), alpha=0.7)
zlim(-0.5, 1.0)