text:
#=== elastic medium 2d ===#
function init2DEL{T1,T2,T3<:Real}(
    pvel::Union{String,<:Real,Array{<:Real,1},StepRange{Int64,Int64},UnitRange{Int64},StepRangeLen{Float32,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}},
    svel::Union{Void,String,Array{<:Real,1},StepRange{Int64,Int64},UnitRange{Int64},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}},
    rho::Union{String,<:Real},
    DZ::Union{<:Real,Array{<:Real,1},StepRange{Int64,Int64},UnitRange{Int64},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}},
    HX::T1,
    pkf::T2,
    T::T3;
    ext = 10,
    iflag = 2,
    dz = nothing,
    dx = nothing,
    dt= nothing,
    nT= nothing)

    medium = model(pvel, svel, rho, DZ, HX, pkf, T, ext, iflag, dx, dz,dt,nT)
    nwf = initnelwf(medium.nDZ, medium.nHX, ext, iflag)
    FDC = fdc(4)
# = unsplit PML boundary = #
    wf = initelwf(medium.nDZ, medium.nHX, ext, iflag)
    pml = BD(nwf, medium)
    # calpara = wfpara(medium)
    return elmod2d(medium, wf, nwf, FDC, pml)#, calpara)
end

function reset!(Model::elmod2d)
    Model.wf = initelwf(Model.medium.nDZ, Model.medium.nHX, Model.medium.ext, Model.medium.iflag)
    Model.pml = BD(Model.nwf, Model.medium)
end

#=== acoustic medium 2d ===#
function init2DAC{T1,T2,T3<:Real}(
    pvel::Union{String,<:Real,Array{<:Real,1},StepRange{Int64,Int64},UnitRange{Int64},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}},
    rho::Union{String,<:Real},
    DZ::Union{<:Real,Array{<:Real,1},StepRange{Int64,Int64},UnitRange{Int64},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}},
    HX::T1,
    pkf::T2,
    T::T3;
    ext = 10,
    iflag = 2,
    dz = nothing,
    dx = nothing,
    dt= nothing,
    nT= nothing)

    medium = model(pvel,rho, DZ, HX, pkf, T, ext, iflag, dx, dz,dt,nT)
    nwf = initnacwf(medium.nDZ, medium.nHX, ext, iflag)
    FDC = fdc(4)
# = unsplit PML boundary = #
    wf = initacwf(medium.nDZ, medium.nHX, ext, iflag)
    pml = BD(nwf, medium)
    # calpara = wfpara(medium)
    return acmod2d(medium, wf, nwf, FDC, pml)#, calpara)
end

function reset!(Model::acmod2d)
    Model.wf = initacwf(Model.medium.nDZ, Model.medium.nHX, Model.medium.ext, Model.medium.iflag)
    Model.pml = BD(Model.nwf, Model.medium)
end

############ elastic medium 3d ############
function init3DEL{T1,T2,T3,T4<:Real}(
    pvel::Union{String,<:Real,Array{<:Real,1},StepRange{Int64,Int64},UnitRange{Int64},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}},
    svel::Union{Void,String,<:Real,Array{<:Real,1},StepRange{Int64,Int64},UnitRange{Int64},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}},
    rho::Union{String,<:Real},
    DZ::Union{<:Real,Array{<:Real,1},StepRange{Int64,Int64},UnitRange{Int64},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}},
    HX::T1,
    HY::T2,
    pkf::T3,
    T::T4;
    ext = 10,
    iflag = 2,
    dz = nothing,
    dx = nothing,
    dy = nothing,
    dt= nothing,
    nT= nothing)

    medium = model(pvel, svel, rho, DZ, HX, HY, pkf, T, ext, iflag, dx, dy, dz, dt,nT)
    nwf = initnelwf(medium.nDZ, medium.nHX, medium.nHY, ext, iflag)
    FDC = fdc(4)
# = unsplit PML boundary = #
    wf = initelwf(medium.nDZ, medium.nHX, medium.nHY, ext, iflag)
    pml = BD(nwf, medium)
    # calpara = wfpara(medium)
    return elmod3d(medium, wf, nwf, FDC, pml)#, calpara)
end

function reset!(Model::elmod3d)
    Model.wf = initelwf(Model.medium.nDZ, Model.medium.nHX, Model.medium.nHY, Model.medium.ext, Model.medium.iflag)
    Model.pml = BD(Model.nwf, Model.medium)
end

############ acoustic medium 3d ############
function init3DAC{T1,T2,T3,T4<:Real}(
    pvel::Union{String,<:Real,Array{<:Real,1},StepRange{Int64,Int64},UnitRange{Int64},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}},
    rho::Union{String,<:Real},
    DZ::Union{<:Real,Array{<:Real,1},StepRange{Int64,Int64},UnitRange{Int64},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}},
    HX::T1,
    HY::T2,
    pkf::T3,
    T::T4;
    ext = 10,
    iflag = 2,
    dz = nothing,
    dx = nothing,
    dy = nothing,
    dt= nothing,
    nT= nothing)

    medium = model(pvel, rho, DZ, HX, HY, pkf, T, ext, iflag, dx, dy, dz, dt,nT)
    nwf = initnacwf(medium.nDZ, medium.nHX, medium.nHY, ext, iflag)
    FDC = fdc(4)
# = unsplit PML boundary = #
    wf = initacwf(medium.nDZ, medium.nHX, medium.nHY, ext, iflag)
    pml = BD(nwf, medium)
    # calpara = wfpara(medium)
    return acmod3d(medium, wf, nwf, FDC, pml)#, calpara)
end

function reset!(Model::acmod3d)
    Model.wf = initacwf(Model.medium.nDZ, Model.medium.nHX, Model.medium.nHY, Model.medium.ext, Model.medium.iflag)
    Model.pml = BD(Model.nwf, Model.medium)
end


text:
# ******************************************************************************************
# Notices:
#
# Copyright © 2022 United States Government as represented by the Administrator of the
# National Aeronautics and Space Administration.  All Rights Reserved.
#
# Disclaimers
#
# No Warranty: THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY KIND,
# EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY
# THAT THE SUBJECT SOFTWARE WILL CONFORM TO SPECIFICATIONS, ANY IMPLIED WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY
# WARRANTY THAT THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT DOCUMENTATION,
# IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE. THIS AGREEMENT DOES NOT, IN ANY MANNER,
# CONSTITUTE AN ENDORSEMENT BY GOVERNMENT AGENCY OR ANY PRIOR RECIPIENT OF ANY RESULTS,
# RESULTING DESIGNS, HARDWARE, SOFTWARE PRODUCTS OR ANY OTHER APPLICATIONS RESULTING FROM
# USE OF THE SUBJECT SOFTWARE.  FURTHER, GOVERNMENT AGENCY DISCLAIMS ALL WARRANTIES AND
# LIABILITIES REGARDING THIRD-PARTY SOFTWARE, IF PRESENT IN THE ORIGINAL SOFTWARE, AND
# DISTRIBUTES IT "AS IS."
#
# Waiver and Indemnity:  RECIPIENT AGREES TO WAIVE ANY AND ALL CLAIMS AGAINST THE UNITED
# STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT.
# IF RECIPIENT'S USE OF THE SUBJECT SOFTWARE RESULTS IN ANY LIABILITIES, DEMANDS, DAMAGES,
# EXPENSES OR LOSSES ARISING FROM SUCH USE, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON,
# OR RESULTING FROM, RECIPIENT'S USE OF THE SUBJECT SOFTWARE, RECIPIENT SHALL INDEMNIFY AND
# HOLD HARMLESS THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL
# AS ANY PRIOR RECIPIENT, TO THE EXTENT PERMITTED BY LAW.  RECIPIENT'S SOLE REMEDY FOR ANY
# SUCH MATTER SHALL BE THE IMMEDIATE, UNILATERAL TERMINATION OF THIS AGREEMENT.
# ******************************************************************************************

"""
Distance heuristic abstract type. Custom heuristics inherit this type and implement `reset!`
and `apply!` functions. Application can return scalar value or partially-applied function
transforming f(s,s') to f(s)(s').
"""
abstract type AbstractDistanceHeuristic end
reset!(::AbstractDistanceHeuristic) = nothing
apply!(::AbstractDistanceHeuristic, ::ASTMDP) = unimplemented()
(adh::AbstractDistanceHeuristic)(mdp::ASTMDP) = Functoid(apply!(adh, mdp))

"""
Gradient of conservative potential. Default and recommended.
"""
struct GradientHeuristic <: AbstractDistanceHeuristic end

function apply!(::GradientHeuristic, mdp::ASTMDP)
    d = distance(mdp.sim)
    return mdp′::ASTMDP -> d - distance(mdp′.sim)
end

"""
Minimum distance across episode. Warning: non-Markovian.
"""
Base.@kwdef mutable struct MinimumHeuristic <: AbstractDistanceHeuristic
    d_min::Float64 = Inf
end

function apply!(h::MinimumHeuristic, mdp::ASTMDP)
    h.d_min = min(h.d_min, distance(mdp.sim))
    return terminated(mdp) ? -h.d_min : 0.0
end

function reset!(h::MinimumHeuristic)
    h.d_min = Inf
end

"""
Final distance encountered. If MDP is episodic, this is equal to the minimum distance
(determined retroactively instead of accumulated as in MinimumHeuristic).
"""
struct FinalHeuristic <: AbstractDistanceHeuristic end
apply!(::FinalHeuristic, mdp::ASTMDP) = isterminal(mdp.sim) ? -distance(mdp.sim) : 0.0

"""
Null heuristic. Returns zero.
"""
struct NullHeuristic <: AbstractDistanceHeuristic end
apply!(::NullHeuristic, ::ASTMDP) = 0.0

ext: jl


text:
using ArgParseLite

function main()
    my_args = Arguments()

    push!(my_args, Argument("arg1"))
    push!(my_args, Argument("--opt1"))
    push!(my_args, Argument("--opt2", "-o"))
    push!(my_args, Argument("--flag1", action=:store_true))

    println("Parsed args:")
    for (arg,val) in ArgParseLite.parse_args(my_args)
        println("  $arg  =>  $val")
    end
end

main()


text:
#################### Slice Sampler ####################

#################### Types and Constructors ####################

typealias SliceForm Union{Univariate, Multivariate}

type SliceTune{F<:SliceForm} <: SamplerTune
  logf::Nullable{Function}
  width::Union{Float64, Vector{Float64}}

  SliceTune() = new()

  SliceTune{T<:Real}(x::Vector, width::ElementOrVector{T}) =
    SliceTune{F}(x, width, Nullable{Function}())

  SliceTune{T<:Real}(x::Vector, width::ElementOrVector{T}, logf::Function) =
    SliceTune{F}(x, width, Nullable{Function}(logf))

  SliceTune(x::Vector, width::Real, logf::Nullable{Function}) =
    new(logf, Float64(width))

  function SliceTune{T<:Real}(x::Vector, width::Vector{T},
                              logf::Nullable{Function})
    new(logf, convert(Vector{Float64}, width))
  end
end

typealias SliceUnivariate SamplerVariate{SliceTune{Univariate}}
typealias SliceMultivariate SamplerVariate{SliceTune{Multivariate}}

validate{F<:SliceForm}(v::SamplerVariate{SliceTune{F}}) =
  validate(v, v.tune.width)

validate{F<:SliceForm}(v::SamplerVariate{SliceTune{F}}, width::Float64) = v

function validate{F<:SliceForm}(v::SamplerVariate{SliceTune{F}}, width::Vector)
  n = length(v)
  length(width) == n ||
    throw(ArgumentError("length(width) differs from variate length $n"))
  v
end

#################### Sampler Constructor ####################

function Slice{T<:Real, F<:SliceForm}(params::ElementOrVector{Symbol},
                                      width::ElementOrVector{T},
                                      ::Type{F}=Multivariate;
                                      transform::Bool=false)
  samplerfx = function(model::Model, block::Integer)
    block = SamplingBlock(model, block, transform)
    v = SamplerVariate(block, width)
    sample!(v, x -> logpdf!(block, x))
    relist(block, v)
  end
  Sampler(params, samplerfx, SliceTune{F}())
end

#################### Sampling Functions ####################

sample!(v::Union{SliceUnivariate, SliceMultivariate}) = sample!(v, v.tune.logf)

function sample!(v::SliceUnivariate, logf::Function)
  logf0 = logf(v.value)

  n = length(v)
  lower = v - v.tune.width .* rand(n)
  upper = lower + v.tune.width

  for i in 1:n
    p0 = logf0 + log(rand())

    x = v[i]
    v[i] = rand(Uniform(lower[i], upper[i]))
    while true
      logf0 = logf(v.value)
      logf0 < p0 || break
      value = v[i]
      if value < x
        lower[i] = value
      else
        upper[i] = value
      end
      v[i] = rand(Uniform(lower[i], upper[i]))
    end
  end

  v
end

function sample!(v::SliceMultivariate, logf::Function)
  p0 = logf(v.value) + log(rand())

  n = length(v)
  lower = v - v.tune.width .* rand(n)
  upper = lower + v.tune.width

  x = v.tune.width .* rand(n) + lower
  while logf(x) < p0
    for i in 1:n
      value = x[i]
      if value < v[i]
        lower[i] = value
      else
        upper[i] = value
      end
      x[i] = rand(Uniform(lower[i], upper[i]))
    end
  end
  v[:] = x

  v
end


text:
type USDr <: StepData

	A # Current sum of gradients. Start with a zeros array
	f # Current array of gradient functions. Start with zero.
	  # Has all gradient function in it, no chunking here
 	y # Current array of compactly stored gradients
 	rf # Current array of restoring fucntions
	I # Current number of total chunks

	s::Function # Number of chunks to resample
	u::Function # Number of new chunks to sample
	beta::Function
	getStep::Function
	stepString::AbstractString
	numVars::Int64
	sagFsagaT::Bool

	function USDr(s::Function, u::Function, beta::Function, numVars::Int64, stepString::AbstractString,sagFsagaT::Bool)
		new(zeros(numVars,1),Function[], Matrix{Float64}[], Function[], 0, s, u, beta, urComputation, stepString, numVars, sagFsagaT)
	end
end

function onlyrAddBeta1(u::Function, uString::AbstractString, numVars::Int64,sagFsagaT::Bool)

	s = (k,I)-> 0
	beta = (k)->1

	stepString  = "urAb1."*uString*"."*string(sagFsagaT)

	USDr(s, u , beta, numVars, stepString,sagFsagaT)

end

function onlyrUpdateBeta1(s::Function, sString::AbstractString, numVars::Int64, ntp::Int64,sagFsagaT::Bool)

	u = (k,I)-> int(floor(k==0 ? int(floor(ntp)) : 0))
	beta = (k)->1

	stepString  = "urUb1."*sString*"."*string(sagFsagaT)

	USDr(s, u , beta, numVars, stepString,sagFsagaT)
end

function urLinBeta1(c::Float64,numVars::Int64, sagFsagaT::Bool)

	if c<1
		error("C must be >=1")
	end
	s = (k,I)-> int(floor(k==0 ? 0 : c))
	u = (k,I)-> int(c)
	beta = (k)->1

	stepString  = "urLb1"*".c="*string(c)*"."*string(sagFsagaT)

	USDr(s, u , beta, numVars, stepString, sagFsagaT)
end

function urQuadBeta1(c::Float64,numVars::Int64,sagFsagaT::Bool)

	s = (k,I)-> int(ceil(k==0 ? 0 : c*k))
	u = (k,I)-> int(ceil(c*(k+1)))
	beta = (k)->1

	stepString  = "urQb1"*".c="*string(c)*"."*string(sagFsagaT)

	USDr(s, u , beta, numVars, stepString,sagFsagaT)
end

function urExpBeta1(r::Float64, numVars::Int64,sagFsagaT::Bool)

	c = 1.0 / (r-1)
	s = (k,I)-> int(ceil(k==0 ? 0 : c*(r/(r-1))^(k-1)))
	u = (k,I)-> int(ceil(k==0 ? c*(r-1) : c*(r/(r-1))^(k-1)))
	beta = (k)->1

	stepString  = "urEb1"*".r="*string(r)*"."*string(sagFsagaT)

	USDr(s, u , beta, numVars, stepString,sagFsagaT)
end

function urExp2Beta1(r::Float64, numVars::Int64,sagFsagaT::Bool)

	c = 1.0
	s = (k,I)-> int(ceil(k==0 ? 0 : c*(r/(r-1))^(k-1)))
	u = (k,I)-> int(ceil(k==0 ? c*(r-1) : c*(r/(r-1))^(k-1)))
	beta = (k)->1

	stepString  = "urE2b1"*".r="*string(r)*"."*string(sagFsagaT)

	USDr(s, u , beta, numVars, stepString,sagFsagaT)
end

function urExpBeta0(r::Float64, numVars::Int64,sagFsagaT::Bool)

	c = 1.0 / (r-1)
	s = (k,I)-> int(ceil(k==0 ? 0 : c*(r/(r-1))^(k-1)))
	u = (k,I)-> int(ceil(k==0 ? c*(r-1) : c*(r/(r-1))^(k-1)))
	beta = (k)->0

	stepString  = "urEb0"*".r="*string(r)*"."*string(sagFsagaT)

	USDr(s, u , beta, numVars, stepString,sagFsagaT)
end

function urComputation(x, k, gnum, sd::USDr, problem::Problem; outputLevel = 0)

	outputLevel > 0 && println("Starting urComputation")

	U = (sd.I+1):(sd.I+sd.u(k,sd.I))

	outputLevel > 0  && println("sd.I=$(sd.I) sd.u(k,sd.I)=$(sd.u(k,sd.I)) sd.s(k,sd.I)=$(sd.s(k,sd.I))")

	for i in U
		try
			# println("before consume")
			# println(problem.getNextSampleFunction)
			kkkkkk=consume(problem.getNextSampleFunction)
			# println("after consume")
			push!(sd.f, kkkkkk[1])
			push!(sd.rf,kkkkkk[2])
		catch
			error("Not producing anymore! Probably out of datapoints.")
		end
	end

	# S is composed of sample indices
	S = StatsBase.sample(1:sd.I,sd.s(k,sd.I);replace=false)

	# B is the current sum of old gradients (for S)
	B = zeros(sd.numVars,1)
	for i in [S]
		B += sd.rf[i](sd.y[i]) # RF is the restoring function
	end

	# sumy is the new sum of gradients (for both U and S)
	sumy = zeros(sd.numVars,1)

	for i in [U]
		(f,sampleG,cs) = (sd.f[i])(x)
		push!(sd.y,cs)
		sumy += sampleG
	end

	for i in [S]
		(f,sampleG,cs) = (sd.f[i])(x)
		push!(sd.y,cs)
		sumy += sampleG
	end

	gnum += sd.s(k,sd.I)+sd.u(k,sd.I)

	if sd.sagFsagaT == false
		g = (sd.A - B + sumy )/(sd.I+sd.u(k,sd.I))
	else
		g = ((sd.s(k,sd.I) > 0  ? sd.s(k,sd.I)*((sd.beta(k)/sd.I)*sd.A) - sd.beta(k)*B : 0) + sumy )/(sd.s(k,sd.I)+sd.u(k,sd.I))
	end

	sd.I = sd.I + sd.u(k,sd.I)

	sd.A=sd.A-B+sumy

	(g, gnum)
end


text:
###############################################################################
#
# Update g_sw using different stomatal models
#
###############################################################################
"""
    prognostic_gsw!(
                clayer::CanopyLayer{FT},
                envir::AirLayer{FT},
                sm::EmpiricalStomatalModel{FT},
                β::FT,
                Δt::FT
    ) where {FT<:AbstractFloat}
    prognostic_gsw!(
                photo_set::AbstractPhotoModelParaSet{FT},
                clayer::CanopyLayer{FT},
                hs::LeafHydraulics{FT},
                envir::AirLayer{FT},
                sm::OSMWang{FT},
                Δt::FT
    ) where {FT<:AbstractFloat}

Update g_sw prognostically, given
- `clayer` A `CanopyLayer` type struct
- `envir` `AirLayer` type environmental conditions
- `sm` `EmpiricalStomatalModel` or `OSMWang` type stomatal model
- `β` Tune factor to stomatal g1. 1 for AbstractBetaV mode
- `Δt` Time interval for prognostic stomatal conductance
- `photo_set` `AbstractPhotoModelParaSet` type photosynthesis model, currently
    supports [`OSMWang`](@ref) model only
- `hs` Leaf hydraulic system
"""
function prognostic_gsw!(
            clayer::CanopyLayer{FT},
            envir::AirLayer{FT},
            sm::EmpiricalStomatalModel{FT},
            β::FT,
            Δt::FT
) where {FT<:AbstractFloat}
    # unpack values
    @unpack g_sw, n_leaf = clayer;

    # update g_sw
    for iLF in 1:n_leaf
        gsw_ss = max(sm.g0,
                     stomatal_conductance(sm, clayer, envir, β, iLF));
        g_sw[iLF] += (gsw_ss - g_sw[iLF]) / clayer.τ_esm * Δt;
    end

    return nothing
end

function prognostic_gsw!(
            photo_set::AbstractPhotoModelParaSet{FT},
            clayer::CanopyLayer{FT},
            hs::LeafHydraulics{FT},
            envir::AirLayer{FT},
            sm::OSMWang{FT},
            Δt::FT
) where {FT<:AbstractFloat}
    # unpack values
    @unpack APAR, g_sw, n_leaf = clayer;

    # update g_sw
    for iLF in 1:n_leaf
        if APAR[iLF] > 1
            diff = solution_diff!(g_sw[iLF], photo_set, clayer, hs, envir, sm,
                                  GswDrive(), iLF);
            g_sw[iLF] += diff * clayer.τ_osm * Δt;
        else
            diff = nocturnal_diff!(g_sw[iLF], photo_set, clayer, envir, sm);
            g_sw[iLF] += diff * clayer.τ_noc * Δt;
        end
    end

    return nothing
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
module PowerAnalyses

using Distributions:
    NoncentralChisq,
    NoncentralF,
    NoncentralT,
    TDist,
    UnivariateDistribution,
    ccdf,
    cdf,
    pdf,
    quantile
using Roots: find_zero

const PKGDIR = string(pkgdir(PowerAnalyses))::String

let
    path = joinpath(PKGDIR, "README.md")
    text = read(path, String)
    @doc text PowerAnalyses
end

include("types.jl")
include("power.jl")

export Tail, one_tail, two_tails

export IndependentSamplesTTest
export DependentSamplesTTest
export OneSampleTTest

export DeviationFromZeroMultipleRegression
export IncreaseMultipleRegression
export OneWayANOVA
export MultifactorFixedEffectsANOVA
export ConstantVectorHotellingTsqTest
export TwoVectorsHotellingTsqTest

export GoodnessOfFitChisqTest
export ConstantVarianceChisqTest

export get_alpha, get_power, get_es, get_n

end # module


text:
using Distributions, SpecialFunctions
import Distributions: mean

""" Exact Bayes (Adams and MacKay, 2007) for the Categorical task (section 2.3.2)
"""
struct TBayesFilter{D, A} <: BayesFilter
    ns::Int
    posteriors::Array{D, 1}
    prior::D
    probabilities::Array{Float64, 1}
    stochasticity::Float64
    changeprobability::Float64
    cutoff::Float64
    approximation::A
    Ps1::Array{Float64, 1}
end
function TBayesFilter(;ns=5, stochasticity = .01,
                    changeprobability = 0.01, cutoff = 1e-16)

    prior = MyDirichlet(ones(ns) * stochasticity)
    approximation = CutoffApprox(cutoff)
    TBayesFilter(ns, [deepcopy(prior)], prior, [1.], stochasticity,
            changeprobability, cutoff, approximation,
            Array{Float64, 1}(undef, ns))
end
meanparameter(f::Union{TBayesFilter, TMPN, TSOR}) = sum(mean.(f.posteriors) .* f.probabilities)
function update_posterior!(f::BayesFilter, x)
    πt= dataprobability.(f.posteriors, x)       # P(x_t|νʳ, χʳ)
    f.probabilities .*= πt
    f.probabilities ./= sum(f.probabilities)
    update!(f.posteriors, x)
end
function update_prior!(f::BayesFilter)
    p_switch = sum(f.probabilities) * f.changeprobability      # ∝ P(r_t = 0, x_{1:t})
    # @show sum(f.probabilities)
    f.probabilities .*= (1 - f.changeprobability)              # ∝ P(r_t = 1, x_{1:t})
    push!(f.probabilities, p_switch)
    push!(f.posteriors, deepcopy(f.prior))
end
function update!(f::Union{TBayesFilter, TMPN, TSOR}, x)
    update_posterior!(f, x)
    approximate!(f.approximation, f)
    parameter_estimate = meanparameter(f)       # E[θ | x_{1:t}]
    update_prior!(f)
    @. f.Ps1 .= copy(parameter_estimate)
end
struct MyDirichlet
    alpha::Vector{Float64}
end
function dataprobability(d::MyDirichlet, x)
    alpha0 = sum(d.alpha)
    tmp = lgamma(d.alpha[x] + 1) - lgamma(d.alpha[x]) +
            lgamma(alpha0) - lgamma(alpha0 + 1)
    exp(tmp)
end
mean(d::MyDirichlet) = d.alpha / sum(d.alpha)
function update!(d::Vector{<:MyDirichlet}, x)
    @inbounds for i in 1:length(d)
        d[i].alpha[x] += 1
    end
    d
end


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

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

Author: Ariel Kellison, 12/2021
"""

fpath = @__DIR__

function plot_errors(in_file::String, err_file::String, hist::Bool, nsteps::Int, xlim::Float64, pbnd_x::Float64, pbnd_v::Float64)
    """
    Args: (1) in_file contains the random inputs for integration
          (2) err_file contains the ouput errors

    Outputs: two plot files to local dir.
          (1) "histo_errors.png", a histogram of errors
          (2) "scatter_errors.png", errors vs input val

    resource: https://docs.julialang.org/en/v1/stdlib/DelimitedFiles/
    """
    # read data as matrices from files

    in_rands = readdlm(in_file, Float64)
    out_xv   = readdlm(err_file, Float64)

    # data to vec

    vec_in       = [x::Float64 for x in in_rands]
    vec_xv       = [x::Float64 for x in out_xv]
    len          = convert(Int64,length(vec_xv)/2)
    in_dat_x     = [x for x in vec_in[1:len]]
    in_dat_v     = [v for v in vec_in[len+1:len*2]]
    out_dat_x    = [x for x in out_xv[1:len]]
    out_dat_v    = [v for v in out_xv[len+1:len*2]]

    nsteps > 1 ? fstr2 = "multi-step (n = $nsteps)" : fstr2 = "single step"

    # plot histograms

    if hist
    	plot([pbnd_x for i in 1:len],seriestype="vline",label="VCFloat bound=$pbnd_x",
		left_margin = 15mm, right_margin = 15mm
	 	)
    	plotx = histogram!(
        	out_dat_x, xlabel = "floating-point error",ylabel = "frequency",
        	label="empirical error",size = (900, 500),bottom_margin = 10mm,left_margin = 10mm,
        	legendfontsize=12, yguidefontsize=12, xguidefontsize=12, xtickfontsize = 12,
        	ytickfontsize = 12,
        	title = "Leapfrog integration: $fstr2 absolute position error"
    	)
    	savefig(plotx,fpath*"/histo_errors_x.png")

   	plot([pbnd_v for i in 1:len],seriestype="vline",label="VCFloat bound=$pbnd_v",
		left_margin = 15mm, right_margin = 15mm
		 )
    	plotv = histogram!(
        	out_dat_v, xlabel = "floating-point error",ylabel = "frequency",
        	label="empirical error",size = (900, 500),bottom_margin = 10mm,left_margin = 10mm,
        	legendfontsize=12, yguidefontsize=12, xguidefontsize=12, xtickfontsize = 12,
        	ytickfontsize = 12,
        	title = "Leapfrog integration: $fstr2 absolute momentum error"
    	)
    	savefig(plotv,fpath*"/histo_errors_v_$nsteps.png")
    end

    plot([pbnd_x for i in -xlim:xlim],seriestype="hline",label="VCFloat bound=$pbnd_x",
        left_margin = 15mm, right_margin = 15mm)

    plotx = scatter!(
        in_dat_x, out_dat_x, ylabel = "floating-point error",xlabel = "initial position",
        size = (900, 500),bottom_margin = 10mm,left_margin = 10mm,
        yguidefontsize=12, xguidefontsize=12, xtickfontsize = 12, ytickfontsize = 12,
        label="empirical error",
        xlims = (-xlim,xlim), title = "Leapfrog integration: $fstr2 absolute position error"
	)

    savefig(plotx,fpath*"/scatter_errors_x_$nsteps.png")

    plot([pbnd_v for i in -xlim:xlim],seriestype="hline",label="VCFloat bound=$pbnd_v",
        left_margin = 15mm, right_margin = 15mm)

    plotv = scatter!(
        in_dat_v, out_dat_v, ylabel = "floating-point error",xlabel = "initial momentum",
        size = (900, 500),bottom_margin = 10mm,left_margin = 10mm,
        yguidefontsize=12, xguidefontsize=12, xtickfontsize = 12, ytickfontsize = 12,
        label="empirical error",
        xlims = (-xlim,xlim), title = "Leapfrog integration: $fstr2 absolute momentum error"
	)
    savefig(plotv,fpath*"/scatter_errors_v_$nsteps.png")

end


text:
# In the 20×20 grid below, four numbers along a diagonal line have been marked
# in red.
#
# 08 02 22 97 38 15 00 40 00 75 04 05 07 78 52 12 50 77 91 08
# 49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48 04 56 62 00
# 81 49 31 73 55 79 14 29 93 71 40 67 53 88 30 03 49 13 36 65
# 52 70 95 23 04 60 11 42 69 24 68 56 01 32 56 71 37 02 36 91
# 22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80
# 24 47 32 60 99 03 45 02 44 75 33 53 78 36 84 20 35 17 12 50
# 32 98 81 28 64 23 67 10 26 38 40 67 59 54 70 66 18 38 64 70
# 67 26 20 68 02 62 12 20 95 63 94 39 63 08 40 91 66 49 94 21
# 24 55 58 05 66 73 99 26 97 17 78 78 96 83 14 88 34 89 63 72
# 21 36 23 09 75 00 76 44 20 45 35 14 00 61 33 97 34 31 33 95
# 78 17 53 28 22 75 31 67 15 94 03 80 04 62 16 14 09 53 56 92
# 16 39 05 42 96 35 31 47 55 58 88 24 00 17 54 24 36 29 85 57
# 86 56 00 48 35 71 89 07 05 44 44 37 44 60 21 58 51 54 17 58
# 19 80 81 68 05 94 47 69 28 73 92 13 86 52 17 77 04 89 55 40
# 04 52 08 83 97 35 99 16 07 97 57 32 16 26 26 79 33 27 98 66
# 88 36 68 87 57 62 20 72 03 46 33 67 46 55 12 32 63 93 53 69
# 04 42 16 73 38 25 39 11 24 94 72 18 08 46 29 32 40 62 76 36
# 20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74 04 36 16
# 20 73 35 29 78 31 90 01 74 31 49 71 48 86 81 16 23 57 05 54
# 01 70 54 71 83 51 54 69 16 92 33 48 61 43 52 01 89 19 67 48
#
# The product of these numbers is 26 × 63 × 78 × 14 = 1788696.
#
# What is the greatest product of four adjacent numbers in the same direction
# (up, down, left, right, or diagonally) in the 20×20 grid?

using ProjectEulerSolutions

# Simply cycle through rows, cols, and diagonals to get max product for adjacent
# numbers.  Not hard, just book-keeping.  Integer values don't promote to
# BigInt, so force BigInt.
function p011solution(matrix::Array{BigInt, 2}=ones(BigInt, 2, 2), adjacent::Integer=2)::BigInt
    nrows, ncols = size(matrix)
    maxprod = 0

    # Check rows
    for row = 1:nrows
        maxrowval = mapreduce(x -> prod(matrix[row, x:x+adjacent-1]), max, 1:(ncols-adjacent+1))
        maxprod = max(maxprod, maxrowval)
    end

    # Check cols
    for col = 1:ncols
        maxcolval = mapreduce(x -> prod(matrix[x:x+adjacent-1, col]), max, 1:(nrows-adjacent+1))
        maxprod = max(maxprod, maxcolval)
    end

    # Check \ diagonals
    for row = 1:nrows-adjacent+1
        for col = 1:ncols-adjacent+1
            localval = matrix[row, col]
            for i = 1:adjacent-1
                localval *= matrix[row + i, col + i]
            end
            # mapreduce is actually a little slower here
            #localval = mapreduce(i->matrix[row+i, col+i], *, 0:adjacent-1)
            maxprod = max(maxprod, localval)
        end
    end

    # Check / diagonals
    for row = 1:nrows-adjacent+1
        for col = adjacent:ncols
            localval = matrix[row, col]
            for i = 1:adjacent-1
               localval *= matrix[row + i, col - i]
            end
            # mapreduce is actually a little slower here
            #localval = mapreduce(i->matrix[row+i, col-i], *, 0:adjacent-1)
            maxprod = max(maxprod, localval)
        end
    end

    return maxprod
end

p011 = Problems.Problem(p011solution)

input = [08 02 22 97 38 15 00 40 00 75 04 05 07 78 52 12 50 77 91 08
         49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48 04 56 62 00
         81 49 31 73 55 79 14 29 93 71 40 67 53 88 30 03 49 13 36 65
         52 70 95 23 04 60 11 42 69 24 68 56 01 32 56 71 37 02 36 91
         22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80
         24 47 32 60 99 03 45 02 44 75 33 53 78 36 84 20 35 17 12 50
         32 98 81 28 64 23 67 10 26 38 40 67 59 54 70 66 18 38 64 70
         67 26 20 68 02 62 12 20 95 63 94 39 63 08 40 91 66 49 94 21
         24 55 58 05 66 73 99 26 97 17 78 78 96 83 14 88 34 89 63 72
         21 36 23 09 75 00 76 44 20 45 35 14 00 61 33 97 34 31 33 95
         78 17 53 28 22 75 31 67 15 94 03 80 04 62 16 14 09 53 56 92
         16 39 05 42 96 35 31 47 55 58 88 24 00 17 54 24 36 29 85 57
         86 56 00 48 35 71 89 07 05 44 44 37 44 60 21 58 51 54 17 58
         19 80 81 68 05 94 47 69 28 73 92 13 86 52 17 77 04 89 55 40
         04 52 08 83 97 35 99 16 07 97 57 32 16 26 26 79 33 27 98 66
         88 36 68 87 57 62 20 72 03 46 33 67 46 55 12 32 63 93 53 69
         04 42 16 73 38 25 39 11 24 94 72 18 08 46 29 32 40 62 76 36
         20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74 04 36 16
         20 73 35 29 78 31 90 01 74 31 49 71 48 86 81 16 23 57 05 54
         01 70 54 71 83 51 54 69 16 92 33 48 61 43 52 01 89 19 67 48]

Problems.benchmark(p011, convert(Array{BigInt}, input), 4)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
test = """
0,9 -> 5,9
8,0 -> 0,8
9,4 -> 3,4
2,2 -> 2,1
7,0 -> 7,4
6,4 -> 2,0
0,9 -> 2,9
3,4 -> 1,4
0,0 -> 8,8
5,5 -> 8,2
"""

function read_input(io)
    lines = readlines(io)
    numbers = hcat(split.(lines, r" -> |,")...) |> permutedims
    return parse.(Int, numbers) .+ 1
end

function board(values::Matrix; diagonal::Bool = false)
    m = maximum(values)
    M = fill(0, m, m)
    for (x1, y1, x2, y2) in eachrow(values)
        if x1 == x2 || y1 == y2
            (x1, x2) = extrema((x1, x2))
            (y1, y2) = extrema((y1, y2))
            M[y1:y2, x1:x2] .+= 1
        elseif diagonal
            for (x, y) in zip(_range(x1, x2), _range(y1, y2))
                M[y, x] += 1
            end
        end
    end
    return M
end

_range(a, b) = range(a, b; step = sign(b - a))

function count_overlaps(values::Matrix; diagonal::Bool = false)
    return count(>(1), board(values; diagonal))
end

let T = read_input(IOBuffer(test))
    @assert count_overlaps(T) == 5
    @assert count_overlaps(T; diagonal = true) == 12
end

M = read_input("day05.in")
println("Q1: ", count_overlaps(M))
println("Q2: ", count_overlaps(M; diagonal = true))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
struct Point{T}
    x::T
    y::T
    i::Int
    new_point::Bool
    function Point(x::T,y::T;i::Int=0,newpoint=false) where T
        @assert length(x) == length(y)
        new{T}(x,y,i,newpoint)
    end
end
function show(io::IO,p::Point{T}) where T
    print(io,"Point of type $T:\n")
    print(io,"      x = $(p.x)\n")
    print(io,"      y = $(p.y)\n")
    print(io,"    idx = $(p.i)\n")
    print(io,"    new = $(p.new_point)\n")
end

"""
# Line

A `Line` is composed of two vectors, `x` and `y`, representing abszissa and ordinate axis in a standard cartesian coordinate system.

## Fields

* `x`: x axis points
* `y`: y axis points
* `n`: number of points in line
* `dom`: the domain of the line, i.e. the range of `x`.
"""
mutable struct Line{T<:Number} <: AbstractArray{T<:Number,1}
    x::Vector{T}
    y::Vector{T}
    n::Int
    dom::Tuple
    function Line(x::Vector{T},y::Vector{T}) where {T<:Number}
        this = new{T}()
        this.x = copy(x)
        this.n = length(x)
        this.y = copy(y)
        this.dom = length(x) >0 ? extrema(x) : (0,0)
        @assert length(y)==this.n
        return this
    end
end
function show(io::IO,L::Line{T}) where {T<:Number}
    print(io,"$T Line\n")
    print(io,"number of grid points: $(L.n)\n")
    print(io,"domain: $(L.dom)\n")
    print(io,"range: $(extrema(L.y))\n")
    if L.n>5
        print(io,"first 5 points:\n")
        print(io,"x = $(L.x[1:5])...\n")
        print(io,"y = $(L.y[1:5])...\n")
    else
        print(io,"x = $(L.x)\n")
        print(io,"y = $(L.y)\n")
    end
end

function reconfigure!(m::Line)
    # after having updated some objects, need to recompute n
    m.n = length(m.x)
    m.dom = extrema(m.x)
    @assert m.n == length(m.y)
end

eltype(l::Line) = eltype(l.x)
size(l::Line) = (l.n,)
length(l::Line) = l.n
function getindex(l::Line,i::Int)
    (l.x[i],l.y[i])
end
function getindex(l::Line,i::UnitRange{Int})
    Line(l.x[i],l.y[i])
end
function setindex!(l::Line{T},x::T,y::T,i::Int) where {T<:Number}
    l.x[i] = x
    l.y[i] = y
end
function setindex!(l::Line{T},x::Vector{T},y::Vector{T},i::UnitRange{Int}) where {T<:Number}
    l.x[i] = x
    l.y[i] = y
end

# interpolating a line

"""
   interp(l::Line{T},ix::Vector{T},extrap::Bool=true) where {T<:Number}

Interpolate a `Line` on a vector of valuels `x`
"""
function interp(l::Line{T},ix::Vector{T},extrap::Bool=false) where {T<:Number}
    # whenever
    xex = extrema(ix)
    # @debug(logger,"interpolating $ix ")
    if xex[1] < l.dom[1] || xex[2] > l.dom[2]
        if extrap
            itp = extrapolate(interpolate((l.x,),l.y,Gridded(Linear())),Linear())
        else
            itp = extrapolate(interpolate((l.x,),l.y,Gridded(Linear())),-Inf)
        end
    else
        itp = interpolate((l.x,),l.y,Gridded(Linear()))
    end
    return itp[ix]
end
function interp(e::Array{Line{T}},ix::Vector{T};extrap::Bool=true) where {T<:Number}
    y = zeros(T,length(e),length(ix))
    for i in eachindex(e)
        y[i,:] = interp(e[i],ix,extrap)
    end
    return y
end

# appending, prepending , deleting and splitting at

"prepend points to an Mline"
function prepend!(m::Line,x,y)
    prepend!(m.x,x)
    prepend!(m.y,y)
    reconfigure!(m)
end

"delete an index"
function delete!(m::Line,idx::Int)
    deleteat!(m.x,idx)
    deleteat!(m.y,idx)
    reconfigure!(m)
end

"append points to an Mline"
function append!(m::Line,x,y)
    append!(m.x,x)
    append!(m.y,y)
    reconfigure!(m)
end

"insert a single value at an interior index"
function insert!(m::Line,vx,vy,idx::Int )
    insert!(m.x,idx,vx)
    insert!(m.y,idx,vy)
    reconfigure!(m)
end

"""
    splitat(m::Line,j::Int,repeat_boundary::Bool=true)

Splits a `Line` object after given index and returns 2 new `Line`s as a tuple. If `repeat_boundary` is true, then the separating index is the first point of the second new `Line`.
"""
function splitat(m::Line,j::Int,repeat_boundary::Bool=true)
    m1 = Line(m.x[1:j],m.y[1:j])
    if repeat_boundary
        m2 = Line(m.x[j:end],m.y[j:end])
    else
        m2 = Line(m.x[j+1:end],m.y[j+1:end])
    end
    return (m1,m2)
end

"sort a `Line` along x-grid"
function sort!(m::Line)
    ix = sortperm(m.x)
    m.x = m.x[ix]
    m.y = m.y[ix]
end

max_stars_repo_name: floswald/MultiLine

max_issues_repo_name: floswald/MultiLine

max_forks_repo_name: floswald/MultiLine


text:
# module Ecto

module Changeset

import ..Schema
import ..Schema: Assoc

type t
    valid::Bool
    repo
    data::Union{Void,Schema.t}
    params::Dict{String,Any}
    changes::Assoc
    required::Vector{Symbol}
    prepare::Vector
    errors::Vector
    constraints::Vector
    validations::Assoc
    filters::Dict{Symbol,Any}
    action
    types::Vector
    t(data) = new(true, nothing, data, Dict{String,Any}(), Assoc([]), Vector{Symbol}(), [], [], [], Assoc([]), Dict{Symbol,Any}(), nothing, [])
end

function change(schema::Schema.t; kw...)::Changeset.t
    changeset = Changeset.t(schema)
    changeset.changes = setdiff(Assoc(kw), schema.struct)
    changeset
end

function change(changeset::Changeset.t; kw...)::Changeset.t
    changeset
end

function change(typ::Type; kw...)::Changeset.t
    mirror = plain_to_mirrorstruct(typ)
    schema = mirrorstruct_to_schema(mirror)
    Changeset.t(schema)
end

function cast(schema::Schema.t, params, allowed)::Changeset.t
    Changeset.t(schema)
end

function cast(assoc::Assoc, params)::Function
    (schema) -> Changeset.t(schema)
end

include("changeset/relation.jl")

end # module Ecto.Changeset

import Base: rem
import .MirrorTypes: plain_to_mirrorstruct, mirrorstruct_to_schema

function rem(schema::Ecto.Schema.t; kw...)::Changeset.t
    changes = Assoc([])
    changeset = Changeset.t(schema)
    changeset.changes = changes
    changeset
end

function rem(::Type{Changeset.t}; kw...)::Changeset.t
    opts = Dict(kw)
    data = opts[:data]
    changeset = Changeset.t(data)
    pop!(opts, :data)
    for (k,v) in kw
        setfield!(changeset, k, v)
    end
    changeset
end

function rem(changeset::Changeset.t; kw...)::Changeset.t
    opts = Dict(kw)
    for (k,v) in opts
        setfield!(changeset, k, v)
    end
    changeset
end

## Errors
type InvalidChangesetError
    messsage
    function InvalidChangesetError(action::Symbol, changeset::Changeset.t)
        new("could not perform $action because changeset is invalid")
    end
end


text:
"""
    NonlinearContact{T,N} <: Contact{T,N}

    contact object for impact and friction with a nonlinear friction cone

    friction_coefficient: value of friction coefficient
    surface_projector: mapping from world frame to surface tangent frame
    surface_normal_projector: inverse/complement of surface_projector
    contact_point: position of contact on Body relative to center of mass
    offset: position of contact relative to contact_point
"""
mutable struct NonlinearContact{T,N} <: Contact{T,N}
    friction_coefficient::T
    surface_projector::SMatrix{2,3,T,6}
    surface_normal_projector::Adjoint{T,SVector{3,T}} # inverse matrix
    contact_point::SVector{3,T}
    offset::SVector{3,T}

    function NonlinearContact(body::Body{T}, normal::AbstractVector, friction_coefficient;
        contact_point=szeros(T, 3),
        offset::AbstractVector=szeros(T, 3)) where T
        V1, V2, V3 = orthogonal_columns(normal)
        A = [V1 V2 V3]
        Ainv = inv(A)
        surface_normal_projector = Ainv[3,SA[1; 2; 3]]'
        surface_projector = SA{T}[
            1 0 0
            0 1 0
        ]
        new{T,8}(friction_coefficient, surface_projector, surface_normal_projector, contact_point, offset)
    end
end

function constraint(mechanism, contact::ContactConstraint{T,N,Nc,Cs}) where {T,N,Nc,Cs<:NonlinearContact{T,N}}
    model = contact.model
    body = get_body(mechanism, contact.parent_id)
    x2, v25, q2, ϕ25 = current_configuration_velocity(body.state)
    x3, q3 = next_configuration(body.state, mechanism.timestep)

    constraint(model, contact.impulses_dual[2], contact.impulses[2], x3, q3, v25, ϕ25)
end

function constraint(model::NonlinearContact, s::AbstractVector{T}, γ::AbstractVector{T},
        x3::AbstractVector{T}, q3::Quaternion{T},
        v25::AbstractVector{T}, ϕ25::AbstractVector{T}) where T

    # transforms the velocities of the origin of the link into velocities
    vp = v25 + skew(vector_rotate(ϕ25, q3)) * (vector_rotate(model.contact_point, q3) - model.offset)
    SVector{4,T}(
        model.surface_normal_projector * (x3 + vector_rotate(model.contact_point, q3) - model.offset) - s[1],
        model.friction_coefficient * γ[1] - γ[2],
        (model.surface_projector * vp - s[@SVector [3,4]])...)
end

function constraint_jacobian_velocity(model::NonlinearContact{T},
    x3::AbstractVector{T}, q3::Quaternion{T},
    x2::AbstractVector{T}, v25::AbstractVector{T}, q2::Quaternion{T}, ϕ25::AbstractVector{T},
    λ, timestep::T) where T
    V = [model.surface_normal_projector * timestep;
        szeros(1,3);
        model.surface_projector]

    ∂v∂q3 = skew(vector_rotate(ϕ25, q3)) * ∂vector_rotate∂q(model.contact_point, q3)
    ∂v∂q3 += skew(model.offset - vector_rotate(model.contact_point, q3)) * ∂vector_rotate∂q(ϕ25, q3)
    ∂v∂ϕ25 = skew(model.offset - vector_rotate(model.contact_point, q3)) * ∂vector_rotate∂p(ϕ25, q3)
    Ω = [model.surface_normal_projector * ∂vector_rotate∂q(model.contact_point, q3) * rotational_integrator_jacobian_velocity(q2, ϕ25, timestep);
        szeros(1,3);
        model.surface_projector * (∂v∂ϕ25 + ∂v∂q3 * rotational_integrator_jacobian_velocity(q2, ϕ25, timestep))]
    return [V Ω]
end

function constraint_jacobian_configuration(model::NonlinearContact{T},
    x3::AbstractVector{T}, q3::Quaternion{T},
    x2::AbstractVector{T}, v25::AbstractVector{T}, q2::Quaternion{T}, ϕ25::AbstractVector{T},
    λ, timestep::T) where T
    X = [model.surface_normal_projector;
        szeros(1,3);
        szeros(2,3)]

    ∂v∂q3 = skew(vector_rotate(ϕ25, q3)) * ∂vector_rotate∂q(model.contact_point, q3)
    ∂v∂q3 += skew(model.offset - vector_rotate(model.contact_point, q3)) * ∂vector_rotate∂q(ϕ25, q3)
    Q = [model.surface_normal_projector * ∂vector_rotate∂q(model.contact_point, q3);
        szeros(1,4);
        model.surface_projector * ∂v∂q3]
    return [X Q]
end

function set_matrix_vector_entries!(mechanism::Mechanism, matrix_entry::Entry, vector_entry::Entry,
    contact::ContactConstraint{T,N,Nc,Cs,N½}) where {T,N,Nc,Cs<:NonlinearContact{T,N},N½}
    # ∇impulses[impulses .* impulses - μ; g - s] = [diag(impulses); -diag(0,1,1)]
    # ∇impulses[impulses .* impulses - μ; g - s] = [diag(impulses); -diag(1,0,0)]
    # (friction_coefficient γ - ψ) dependent of ψ = impulses[2][1:1]
    # B(z) * zdot - sβ dependent of sβ = impulses[2][2:end]
    friction_coefficient = contact.model.friction_coefficient
    γ = contact.impulses[2] + REG * neutral_vector(contact.model) # TODO need to check this is legit
    s = contact.impulses_dual[2] + REG * neutral_vector(contact.model) # TODO need to check this is legit

    # ∇s = [contact.impulses[2][1] szeros(1,3); szeros(3,1) cone_product_jacobian(contact.impulses[2][2:4]); Diagonal([-1, 0, -1, -1])]
    ∇s1 = [γ[SA[1]]; szeros(T,3)]'
    ∇s2 = [szeros(T,3,1) cone_product_jacobian(γ[SA[2,3,4]])]
    ∇s3 = Diagonal(SVector{4,T}(-1, 0, -1, -1))
    ∇s = [∇s1; ∇s2; ∇s3]

    # ∇γ = [contact.impulses_dual[2][1] szeros(1,3); szeros(3,1) cone_product_jacobian(contact.impulses_dual[2][2:4]); szeros(1,4); friction_coefficient -1 0 0; szeros(2,4)]
    ∇γ1 = [s[SA[1]]; szeros(T,3)]'
    ∇γ2 = [szeros(T,3,1) cone_product_jacobian(s[SA[2,3,4]])]
    ∇γ3 = SA[0   0 0 0;
             friction_coefficient -1 0 0;
             0   0 0 0;
             0   0 0 0;]
    ∇γ = [∇γ1; ∇γ2; ∇γ3]

    # matrix_entry.value = [[contact.impulses[2][1] szeros(1,3); szeros(3,1) cone_product_jacobian(contact.impulses[2][2:4]); Diagonal([-1, 0, -1, -1])] [contact.impulses_dual[2][1] szeros(1,3); szeros(3,1) cone_product_jacobian(contact.impulses_dual[2][2:4]); szeros(1,4); friction_coefficient -1 0 0; szeros(2,4)]]
    matrix_entry.value = [∇s ∇γ]

    # [-impulses .* impulses + μ; -g + s]
    vector_entry.value = vcat(-complementarityμ(mechanism, contact), -constraint(mechanism, contact))
    return
end

neutral_vector(model::NonlinearContact{T,N}) where {T,N} = [sones(T, 2); szeros(T, Int(N/2) -2)]

cone_degree(model::NonlinearContact) = 2


text:
using Gadfly

set_default_plot_size(14cm, 8cm)

# Issue #1344

datas = Gadfly.Data(color= ["D", "A", "C", "D", "A", "C", "D", "D", "A", "B"],
        shape=["E", "F", "F", "E", "G", "E", "E", "G", "H"])
scales = [Scale.color_discrete(levels=["A","B","C","D"]), Scale.shape_discrete(levels=["E","F","G","H"]) ]
aes = Scale.apply_scales(scales, datas)
theme1 = Theme(key_max_columns=4)
guides = [render(g, theme1, aes[1])[1].ctxs  for g in (Guide.colorkey(), Guide.shapekey())]

gridstack([guides[i][j] for i in 1:2, j in 1:4])


text:
# simple program to test affinity propagation

using Distance
using Clustering

x = rand(100, 500)
S = -pairwise(Euclidean(), x, x)

# set diagonal value to median value
S = S - diagm(diag(S)) + median(S)*eye(size(S,1))

result = affinity_propagation(S; max_iter=200, display=:iter)


text:
module MANA
# Packge Main for MANA GWAS package
using Distributed, DelimitedFiles, Statistics, LinearAlgebra, Random, SparseArrays, BEDFiles, Distributions, JuliaDB, IndexedTables, OnlineStats, JWAS, DataFrames
#
include("types.jl")
include("read/read.jl")
include("model/model.jl")
include("mutate/mutate.jl")
include("run/run.jl")
include("summarise/summarise.jl")
include("ycreator/ycreator.jl")
include("pedigree/pedigree.jl")

# Functions for basic model
export build_gwas
export add_snps
export setup_phenotypes
export initiate_model
export run_model
export concatenate_summaries
# Functions for adjusting marker effects
export make_alphas
export filter_alphas
export print_alphas
# Functions for adjusting phenotypes
export build_ycreator
export add_phenotypes
export add_alphas
export add_markers_gensel
export sort_ycreator
export initiate_ycreator
export create_yk_loso
# Functions for pedigree models
export run_pedigree_model
export run_spliceQTL
export run_pedigree_model_multirecord
export make_genotype_table
export generate_genotypeclass_columns
export center_phenotype

end

max_stars_repo_name: egmreynolds/MANA

max_issues_repo_name: egmreynolds/MANA

max_forks_repo_name: egmreynolds/MANA


text:
################################################################################
#
#   random.jl : random functions
#
################################################################################

# Random suite over finite fields in Nemo

export randomPolynomial

"""
    randomElem(finField::Nemo.FinField)

Return a random element in `finField`.
"""
function randomElem(finField::Nemo.FinField)
    x = gen(finField)
    c::Int = characteristic(finField) - 1
    n = finField.mod_length - 2
    res = finField(rand(0:c))
    for j in 1:n
        res += rand(0:c)*x^j
    end
    return res
end

"""
    randomList(finField::Nemo.FinField, len::Integer)

Return an `Array` of length `len` with random elements in `finField`.
"""
function randomList(finField::Nemo.FinField, len::Integer)
    A = Array(finField, len)
    for i in 1:len
        A[i] = randomElem(finField)
    end
    return A
end

"""
    randomPolynomial(polyRing::Nemo.PolyRing, degree::Integer)

Return a random polynomial of degree `degree` in the ring `polyRing`.
"""
function randomPolynomial(polyRing::Nemo.PolyRing, degree::Integer)
    L = randomList(base_ring(polyRing), degree + 1)
    while L[degree + 1] == 0
        L = randomList(base_ring(polyRing), degree + 1)
    end
    return polyRing(L)
end

"""
    randomIrrPolynomial(polyRing::Nemo.PolyRing, degree::Integer)

Return a randon irreducible polynomial of degree `degree` in the ring
`polyRing`.
"""
function randomIrrPolynomial(polyRing::Nemo.PolyRing, degree::Integer)
    P = randomPolynomial(polyRing, degree)
    while length(factor(P)) != 1
        P = randomPolynomial(polyRing, degree)
    end
    return P
end

# Random B such that X^(q+1) - BX + B factors

"""
    randomSplitElem(polyRing::Nemo.Ring, q::Int)

Return a random element B in the base ring of `polyRing` such that the
polynomial X^(q+1) - BX + B splits into linear factors.

See Reference **[4]** for an explanation about the formula. The references can
be found in the module documentation.
"""
function randomSplitElem(polyRing::Nemo.Ring, q::Int)

    # We pick a random element in F_q^k\F_q^2, where F_q^k is the base ring of
    # `polyRing`
    F = base_ring(polyRing)
    randEl = randomElem(F)
    while randEl == randEl^(q^2)
        randEl = randomElem(F)
    end

    # Then we apply the formula
    res = (randEl-randEl^(q^2))^(q+1)*(randEl-randEl^q)^(-q^2-1)

    return res
end


text:
module SlaterKosterTests

using Test
using SlaterKoster

const DIR = joinpath(@__DIR__, "datasets")

@testset "Load" begin
    @test SlaterKoster.load(joinpath(DIR, "C-C.skf")) == SlaterKoster.load(Float64, joinpath(DIR, "C-C.skf"))

    for T in (Float32, Float64)
        skt = SlaterKoster.loaddir(T, DIR)
        cc = SlaterKoster.load(T, joinpath(DIR, "C-C.skf"))
        ch = SlaterKoster.load(T, joinpath(DIR, "C-H.skf"), element_symbols=(:C, :H))

        @test skt[(:C, :C)] == cc
        @test skt[(:C, :H)] == ch
    end
end

end


text:
"""
    abstract type AbstractImpurityDensity{T <: SSDFloat} end

Struct defining an impurity density inside a [`Semiconductor`](@ref).

For each impurity density, there should be a method
[`get_impurity_density`](@ref) which returns the impurity density in SI units (1/m³) at a given point.

## Examples
* [`ConstantImpurityDensity`](@ref)
* [`LinearImpurityDensity`](@ref)
* [`CylindricalImpurityDensity`](@ref)

!!! note
    The sign of the impurity density is important. It is taken into account in the conversion to a charge density
    and, thus, defines where the semiconductor is n-type or p-type.
"""
abstract type AbstractImpurityDensity{T <: SSDFloat} end

@inline function ImpurityDensity(T::DataType, dict::Union{Dict{String, Any}, Dict{Any, Any}}, input_units::NamedTuple)
    return ImpurityDensity(T, Val{Symbol(dict["name"])}(), dict, input_units)
end

"""
    get_impurity_density(id::AbstractImpurityDensity, pt::AbstractCoordinatePoint)

Returns the impurity density at a given point, `pt`, based on the impurity density model `id`.

## Arguments
* `id::AbstractImpurityDensity`: The [`AbstractImpurityDensity`](@ref) defining the impurity density inside a [`Semiconductor`](@ref).
* `pt::AbstractCoordinatePoint`: The point at which `id` is to be evaluated.

!!! note
    The value returned by `get_impurity_density` is in units of 1/m³ (SI Units).
"""
function get_impurity_density end

include("ConstantImpurityDensity.jl")
include("LinearImpurityDensity.jl")
include("CylindricalImpurityDensity.jl")


text:
using ReplMaker

function apl_reader(str)
    parse_apl(str) |> eval_apl
end

function init_repl(; prompt_text="APL.jl> ", prompt_color = :red, start_key = "/", sticky=true)
    ReplMaker.initrepl(str -> eval_apl(parse_apl(str)),
	               repl = Base.active_repl,
	               prompt_text = prompt_text,
	               prompt_color = prompt_color,
	               start_key = start_key,
                       sticky_mode=sticky,
	               mode_name = "APL Mode")
end


text:
using DrinkGuiden
using Test

@testset "DrinkGuiden.jl" begin

@testset "Datalayer" begin
include("testdata.jl")
end

@testset "Logiclayer" begin
include("testlogic.jl")
end

@testset "Presentationlayer" begin
include("testpresentation.jl")
end

end


text:
struct CustomEvent{K, T} <: AbstractEvent
    state::Dict{K, T}
    extras::Dict
    old_state::Function # f(e, args...) compute the old state
    new_state::Function # f(e, args...) compute the new state
    update!::Function    # f(e, args...) sync old and new states
    trigger::Function   # f(e, args...) definition that the event happens
end

function CustomEvent{K, T}(;
        state = Dict{K, T}(),
        extras = Dict(),
        old_state, new_state, update!, trigger
    ) where {K, T}
    CustomEvent{K, T}(state, extras, old_state, new_state, update!, trigger)
end

states(e::CustomEvent) = e.state
extras(e::CustomEvent) = e.extras
old_state(e::CustomEvent, args...) = e.old_state(e, args...)
new_state(e::CustomEvent, args...) = e.new_state(e, args...)
trigger(e::CustomEvent, args...) = e.trigger(e, args...)
update!(e::CustomEvent, arg, args...) = e.update!(e, arg, args...)
update!() = update!(e, keys(states(e)))
reset!(e::CustomEvent) = (empty!(extras(e)); empty!(states(e)))


text:
using MetabolicFluxAnalysis
using CoordinateTransformations

    # current tests..
    fluxi_1_f = open("testdata/test_01.fluxi")
    fluxi_1_s = readstring(fluxi_1_f)
    fluxi1 = MetabolicFluxAnalysis.parseFluxi( fluxi_1_s )

    fluxi1.net.fluxes_db[fluxi1.net.fluxes[1]].transition
    emudecomp_1 = MetabolicFluxAnalysis.compute_emu_decomp( fluxi1.net , [ fluxi1.measurements["MassDistribution_01"][1].sass ] )
    emu_graph_1 = MetabolicFluxAnalysis.create_emu_graph(emudecomp_1)

    configured_sas = MetabolicFluxAnalysis.get_configured_reactands(fluxi1,"SC_010")
    #fluxi_ast_1 = MetabolicFluxAnalysis.create_ast_from_emu_graph( emu_graph_1 , configured_sas )
    fluxi_ast_1 = MetabolicFluxAnalysis.create_ast_from_emu_graph( emu_graph_1 )

    # get the target semus from the fluxi measurements..
    target_smrs = MetabolicFluxAnalysis.construct_all_mass_ratios_for_measurement(fluxi1,"MassDistribution_01")
    julia_fw_sim_safe = MetabolicFluxAnalysis.ast_to_juliacode_01(fluxi1.net,fluxi_ast_1,target_smrs)
    julia_fw_sim_fast = MetabolicFluxAnalysis.ast_to_juliacode_01(fluxi1.net,fluxi_ast_1,target_smrs;fast_lu=true)

    ocl_fw_sim_fast = MetabolicFluxAnalysis.ast_to_opencl_01(fluxi1.net,fluxi_ast_1,target_smrs)

    julia_fw_sim_1 = julia_fw_sim_safe

    println("\n\nJULIA CODE:\n\n")
    print(julia_fw_sim_safe.code)

    jfw1 = eval( parse( julia_fw_sim_safe.code ) )
    jfw1(rand(11),rand(9))

    mds_a =  MetabolicFluxAnalysis.construct_mds(fluxi1.substrate_configs["SC_010"],julia_fw_sim_1.semu_mass_ratios)
    print("\nInput Config:\n")
    for zi=1:length(mds_a)
        print("$(julia_fw_sim_1.semu_mass_ratios[zi]) --> $(mds_a[zi])\n")
    end

    # check if everythings right:
    #                            v1  v2    v3   v4   v5   v6      Ain  Eout  Fout
    results_test = jfw1(mds_a,[ 100.; 110 ; 50 ; 20 ; 20 ; 80   ;  100 ; 60 ; 80 ])

    #m_ref_01 = fluxi1.measurements["MassDistribution_01"]
    m1_ref = MetabolicFluxAnalysis.get_mass_distribution_measurement_value(fluxi1,"MassDistribution_01",target_smrs[1])
    m2_ref = MetabolicFluxAnalysis.get_mass_distribution_measurement_value(fluxi1,"MassDistribution_01",target_smrs[2])
    m3_ref = MetabolicFluxAnalysis.get_mass_distribution_measurement_value(fluxi1,"MassDistribution_01",target_smrs[3])
    m4_ref = MetabolicFluxAnalysis.get_mass_distribution_measurement_value(fluxi1,"MassDistribution_01",target_smrs[4])

    @test abs( results_test[1] - m1_ref[1,1] ) < 0.0001
    @test abs( results_test[2] - m2_ref[1,1] ) < 0.0001
    @test abs( results_test[3] - m3_ref[1,1] ) < 0.0001
    @test abs( results_test[4] - m4_ref[1,1] ) < 0.0001

    #reload("MetabolicFluxAnalysis")
    # now get and then parametrize the flux polytope
    fp = MetabolicFluxAnalysis.get_flux_polyhedron( fluxi1.net , fluxi1.linear_constraints["ConstraintsExchange"] , max_val=200 )

    fp_min = MetabolicFluxAnalysis.find_minimal_affine_subspace(fp)

    if(false)
        # benchmark: # currently this does not work nicely.. we need steady state flux distr. to get ok results..
        n_benchmark_runs = 200_000
        fluxes_i  = rand(9,n_benchmark_runs) + 0.02;
        #smrs_i    = rand(11,n_benchmark_runs) + 0.005;
        results_i = NaN * ones(4,n_benchmark_runs);
        tic()
            for zi=1:size(fluxes_i,2)
                results_i[:,zi] = jfw1(mds_a,fluxes_i[:,zi])
            end
        time_1 = toc()
        println("Eval Speed: $(n_benchmark_runs/time_1) / sec")
    end

# run this benchmark..
MetabolicFluxAnalysis.benchmark_matrix_solving_01()
# and run this..
MetabolicFluxAnalysis.test_lu_solving_01()
##

#wccs = weakly_connected_components( gg1.g )


text:
module ZWarpedArrays

using Images, AxisArrays, ImageTransformations, Interpolations, CoordinateTransformations
const axes = Base.axes #for name conflict with AxisArrays

using Unitful #just for match_axisspacing

using CachedArrays
import CachedArrays: AbstractCachedArray,
                        update_cache!,
                        parent,
                        cache,
                        current_I,
                        set_I!,
                        cached_axes,
                        noncached_axes,
                        axisspacing,
                        match_axisspacing
import Base: show

export ZWarpedArray, warp_and_resample

include("util.jl")
include("zwarpedarray.jl")

end # module


text:
module Condor

#using Falcons
using Healpix
using PyCall
using Base.Threads
using StaticArrays
using ReferenceFrameRotations
using LinearAlgebra
using ProgressMeter
using HDF5
using Interpolations
const hp = pyimport("healpy")
np = pyimport("numpy")
pd = pyimport("pandas")
using WignerD

function __init__()
    copy!(hp, pyimport_conda("healpy", "healpy"))
end

include("./function/beam_func.jl")
include("./function/convolution_func.jl")
include("./function/adding_YN_func.jl")

export get_2Ddata, symmetrizer, dBi, gen_beammap, beam_pointor, AlmPair, gen_Blm, gen_GaussBeammap
export EffectiveBeamConvolution, GaussBeamConvolution, convolvor

end


text:
using Distributions
import Sobol

#=
References
----------
    [1] Sobol, I. M. (2001).  "Global sensitivity indices for nonlinear
        mathematical models and their Monte Carlo estimates."  Mathematics
        and Computers in Simulation, 55(1-3):271-280,
        doi:10.1016/S0378-4754(00)00270-6.
    [2] Saltelli, A. (2002).  "Making best use of model evaluations to
        compute sensitivity indices."  Computer Physics Communications,
        145(2):280-297, doi:10.1016/S0010-4655(02)00280-1.
    [3] Saltelli, A., P. Annoni, I. Azzini, F. Campolongo, M. Ratto, and
        S. Tarantola (2010).  "Variance based sensitivity analysis of model
        output.  Design and estimator for the total sensitivity index."
        Computer Physics Communications, 181(2):259-270,
        doi:10.1016/j.cpc.2009.09.018.
=#

"""
    sample(data::SobolData)

Generate a matrix containing the model inputs for Sobol sensitivity analysis with
the information in the `data`. In this function we apply Saltelli's
extension of the Sobol  sequence. Saltelli's scheme extends the Sobol sequence in
a way to reduce the error rates in the resulting sensitivity index calculations.
The resulting matrix has `N` * (`D` + 2) rows, where `D` is the number of parameters
and `N` is the number of samples.
"""
function sample(data::SobolData)

    # check for parameters
    if data.params === nothing
        error("Cannot sample to generate model inputs because there are no parameters")
    end

    # set number of values to skip from the initial sequence
    numskip = 1000

    # constants
    calc_second_order = data.calc_second_order
    D = length(data.params) # number of uncertain parameters in problem
    N = data.N # number of samples

    seq = Sobol.SobolSeq(2 * D)
    base_seq = hcat([Sobol.next!(seq) for i = 1:N + numskip - 1]...)'
    base_seq = base_seq[numskip:end, :] # SALIb includes first row of zeros, so skip one less
    scale_sobol_seq!(base_seq, [values(data.params)...]) # scale

    # create the Saltelli sequence
    if calc_second_order
        saltelli_seq = Array{Float64}(undef, N * (2 * D + 2), D)
    else
        saltelli_seq = Array{Float64}(undef, N * (D + 2), D)
    end

    # The Saltelli sequence is made up of N blocks of (D + 2) rows (if don't
    # calculate second order), or (2 * D +2) rows (if do calculate second order)
    #  where each block j contains a first row from A_j, a last row B_j and D
    # middle rows that form AB_j. [Saltelli et al., 2010 see Radial Sampling].

    index = 1

    for i in 1:N

        # copy matrix "A" (first row of each block)
        for j in 1:D
            saltelli_seq[index, j] = base_seq[i, j]
        end
        index += 1

        # for each parameter, place elements of "B" into "A" and insert those D rows
        # of "AB" (top set of middle rows of each block)
        for k in 1:D
            for j in 1:D
                if j == k
                    saltelli_seq[index, j] = base_seq[i, j + D] # from B
                else
                    saltelli_seq[index, j] = base_seq[i, j] # from A
                end
            end
            index += 1
        end

        # for each parameter, place elements of "A" into "B" and insert those D
        # rows of "BA" (bottom set of middle rows of each block)
        # Only needed if you're doing second-order indices
        if calc_second_order
            for k in 1:D
                for j in 1:D
                    if j == k
                        saltelli_seq[index, j] = base_seq[i, j]
                    else
                        saltelli_seq[index, j] = base_seq[i, j + D]
                    end
                end
                index += 1
            end
        end

        # copy matrix "B" (last row of each block)
        for j in 1:D
            saltelli_seq[index, j] = base_seq[i, j + D]
        end
        index += 1
    end

    return saltelli_seq
end


text:
## COV_EXCL_START

# TODO
# - serial version for lower latency
# - block-stride loop to delay need for second kernel launch

# Reduce a value across a warp
@inline function reduce_warp(op, val)
    # offset = CUDAnative.warpsize() ÷ 2
    # while offset > 0
    #     val = op(val, shfl_down_sync(0xffffffff, val, offset))
    #     offset ÷= 2
    # end

    # Loop unrolling for warpsize = 32
    val = op(val, shfl_down_sync(0xffffffff, val, 16, 32))
    val = op(val, shfl_down_sync(0xffffffff, val, 8, 32))
    val = op(val, shfl_down_sync(0xffffffff, val, 4, 32))
    val = op(val, shfl_down_sync(0xffffffff, val, 2, 32))
    val = op(val, shfl_down_sync(0xffffffff, val, 1, 32))

    return val
end

# Reduce a value across a block, using shared memory for communication
@inline function reduce_block(op, val::T, neutral, shuffle::Val{true}) where T
    # shared mem for 32 partial sums
    shared = @cuStaticSharedMem(T, 32)  # NOTE: this is an upper bound; better detect it

    wid, lane = fldmod1(threadIdx().x, CUDAnative.warpsize())

    # each warp performs partial reduction
    val = reduce_warp(op, val)

    # write reduced value to shared memory
    if lane == 1
        @inbounds shared[wid] = val
    end

    # wait for all partial reductions
    sync_threads()

    # read from shared memory only if that warp existed
    val = if threadIdx().x <= fld1(blockDim().x, CUDAnative.warpsize())
         @inbounds shared[lane]
    else
        neutral
    end

    # final reduce within first warp
    if wid == 1
        val = reduce_warp(op, val)
    end

    return val
end
@inline function reduce_block(op, val::T, neutral, shuffle::Val{false}) where T
    threads = blockDim().x
    thread = threadIdx().x

    # shared mem for a complete reduction
    shared = @cuDynamicSharedMem(T, (2*threads,))
    @inbounds shared[thread] = val

    # perform a reduction
    d = threads>>1
    while d > 0
        sync_threads()
        if thread <= d
            shared[thread] = op(shared[thread], shared[thread+d])
        end
        d >>= 1
    end

    # load the final value on the first thread
    if thread == 1
        val = @inbounds shared[thread]
    end

    return val
end

# Partially reduce an array across the grid. The reduction is partial, with multiple
# blocks `gridDim_reduce` working on reducing data from `A` and writing it to multiple
# outputs in `R`. All elements to be processed can be addressed by the product of the
# two iterators `Rreduce` and `Rother`, where the latter iterator will have singleton
# entries for the dimensions that should be reduced (and vice versa). The output array
# is expected to have an additional dimension with as size the number of reduced values
# for every reduction (i.e. more than one if there's multiple blocks participating).
function partial_mapreduce_grid(f, op, A, R, neutral, Rreduce, Rother, gridDim_reduce, shuffle)
    # decompose the 1D hardware indices into separate ones for reduction (across threads
    # and possibly blocks if it doesn't fit) and other elements (remaining blocks)
    threadIdx_reduce = threadIdx().x
    blockDim_reduce = blockDim().x
    blockIdx_other, blockIdx_reduce = fldmod1(blockIdx().x, gridDim_reduce)

    # block-based indexing into the values outside of the reduction dimension
    # (that means we can safely synchronize threads within this block)
    iother = blockIdx_other
    @inbounds if iother <= length(Rother)
        Iother = Rother[iother]

        # load the neutral value
        Iout = CartesianIndex(Tuple(Iother)..., blockIdx_reduce)
        neutral = if neutral === nothing
            R[Iout]
        else
            neutral
        end

        # get a value that should be reduced
        ireduce = threadIdx_reduce + (blockIdx_reduce - 1) * blockDim_reduce
        val = if ireduce <= length(Rreduce)
            Ireduce = Rreduce[ireduce]
            J = max(Iother, Ireduce)
            f(A[J])
        else
            neutral
        end
        val = op(val, neutral)

        val = reduce_block(op, val, neutral, shuffle)

        # write back to memory
        if threadIdx_reduce == 1
            R[Iout] = val
        end
    end

    return
end

## COV_EXCL_STOP

NVTX.@range function GPUArrays.mapreducedim!(f, op, R::CuArray{T}, A::AbstractArray, init=nothing) where T
    Base.check_reducedims(R, A)
    isempty(A) && return R

    f = cufunc(f)
    op = cufunc(op)

    # be conservative about using shuffle instructions
    shuffle = true
    shuffle &= capability(device()) >= v"3.0"
    shuffle &= T in (Int32, Int64, Float32, Float64, ComplexF32, ComplexF64)
    # TODO: add support for Bool (CUDAnative.jl#420)

    # iteration domain, split in two: one part covers the dimensions that should
    # be reduced, and the other covers the rest. combining both covers all values.
    Rall = CartesianIndices(A)
    Rother = CartesianIndices(R)
    Rreduce = CartesianIndices(ifelse.(axes(A) .== axes(R), Ref(Base.OneTo(1)), axes(A)))
    # NOTE: we hard-code `OneTo` (`first.(axes(A))` would work too) or we get a
    #       CartesianIndices object with UnitRanges that behave badly on the GPU.
    @assert length(Rall) == length(Rother) * length(Rreduce)

    # allocate an additional, empty dimension to write the reduced value to.
    # this does not affect the actual location in memory of the final values,
    # but allows us to write a generalized kernel supporting partial reductions.
    R′ = reshape(R, (size(R)..., 1))

    # determine how many threads we can launch
    args = (f, op, A, R′, init, Rreduce, Rother, 1, Val(shuffle))
    kernel_args = cudaconvert.(args)
    kernel_tt = Tuple{Core.Typeof.(kernel_args)...}
    kernel = cufunction(partial_mapreduce_grid, kernel_tt)
    kernel_config =
        launch_configuration(kernel.fun; shmem = shuffle ? 0 : threads->2*threads*sizeof(T))

    # determine the launch configuration
    dev = device()
    reduce_threads = shuffle ? nextwarp(dev, length(Rreduce)) : nextpow(2, length(Rreduce))
    if reduce_threads > kernel_config.threads
        reduce_threads = shuffle ? prevwarp(dev, kernel_config.threads) : prevpow(2, kernel_config.threads)
    end
    reduce_blocks = cld(length(Rreduce), reduce_threads)
    other_blocks = length(Rother)
    threads, blocks = reduce_threads, reduce_blocks*other_blocks
    shmem = shuffle ? 0 : 2*threads*sizeof(T)

    # perform the actual reduction
    if reduce_blocks == 1
        # we can cover the dimensions to reduce using a single block
        @cuda threads=threads blocks=blocks shmem=shmem partial_mapreduce_grid(
            f, op, A, R′, init, Rreduce, Rother, 1, Val(shuffle))
    else
        # we need multiple steps to cover all values to reduce
        partial = similar(R, (size(R)..., reduce_blocks))
        if init === nothing
            # without an explicit initializer we need to copy from the output container
            sz = prod(size(R))
            for i in 1:reduce_blocks
                # TODO: async copies (or async fill!, but then we'd need to load first)
                #       or maybe just broadcast since that extends singleton dimensions
                copyto!(partial, (i-1)*sz+1, R, 1, sz)
            end
        end
        @cuda threads=threads blocks=blocks shmem=shmem partial_mapreduce_grid(
            f, op, A, partial, init, Rreduce, Rother, reduce_blocks, Val(shuffle))

        GPUArrays.mapreducedim!(identity, op, R′, partial, init)
    end

    return R
end


text:
using AlignmentStatistics
using Base.Test
using FastaIO

# write your own tests here
aaix = read_AAindex1() #test func.
@test size(aaix) == (534,22) #change the 534 accordingly if more indexes are contained

daaix = AAindex1_to_Dict(aaix,1) #test func.
@test daaix["A"] == 4.35
@test daaix["F"] == 4.66

fasta = readfasta("interleukin-4-aln.fa")
labs = get_sequence_labels(fasta) #test func.
@test labs[1] == "P07750"
@test labs[length(labs)] == "O77762"

run(`touch test.fa`)
run(`rm -f test.fa`)
rfasta = rectangularize_alignment(fasta) #test func.
@test rfasta[2,1] == "M"
@test rfasta[2,155] == "S"
@test rfasta[3,155] == "-"
export_fasta("test.fa",labs,rfasta) #test func.
fasta2 = readfasta("test.fa")
@test fasta2 == fasta

sfasta = split_fasta_sequences(fasta)
x = sequence_composition(sfasta[1][1:10]) #test func.
@test x["G"] == 1
@test x["V"] == 2
@test x["L"] == 2

@test PDB_AA_to_1_letter_code["ALA"] == "A" #test Dict
@test PDB_AA_to_1_letter_code["GLN"] == "Q"

seqs = readfasta("interleukin-4.fa")
seqs_only = split_fasta_sequences(seqs)
nseqs = sequences_to_numerical_properties(daaix, seqs_only) #test func.
@test nseqs[1][1] == 4.52

x = sequence_lengths(seqs) #test func.
@test length(x) == 5
@test x[135] == 1

sis = symbols_in_sequences(sfasta[1], rfasta) #test func.
@test sis[2][1] == 5
@test sis[2][length(sis[2])] == 4
@test sis[3][length(sis[2])] == 0.8

mc = majority_consensus(rfasta) #test func.
@test mc[2][1] == 5
@test mc[2][length(mc[2])] == 4
@test mc[3][length(mc[2])] == 0.8

lfasta = readfasta("interleukin-4-labs.fa")
labs2 = extract_label_element(get_sequence_labels(lfasta),".",2) #test func.
slabs2 = sorted_label_frequencies(labs2) #test func.
@test slabs2[1][1] == "alpha"
@test slabs2[2][2] == 2


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

module Directions
export Constant, RectGrid, UniformCone, HexapolarCone

using ....OpticSim
using ...Emitters
using ...Geometry
using LinearAlgebra

abstract type AbstractDirectionDistribution{T<:Real} end

Base.iterate(a::AbstractDirectionDistribution, state = 1) = state > length(a) ? nothing : (generate(a, state - 1), state + 1)
Base.getindex(a::AbstractDirectionDistribution, index) = generate(a, index)
Base.firstindex(a::AbstractDirectionDistribution) = 0
Base.lastindex(a::AbstractDirectionDistribution) = length(a) - 1
Base.copy(a::AbstractDirectionDistribution) = a # most don't have any heap allocated stuff so don't really need copying

"""
    Constant{T} <: AbstractDirectionDistribution{T}

Encapsulates a single ray direction, where the default direction is unitZ3 [0, 0, 1].

```julia
Constant(direction::Vec3{T}) where {T<:Real}
Constant(::Type{T} = Float64) where {T<:Real}
```
"""
struct Constant{T} <: AbstractDirectionDistribution{T}
    direction::Vec3{T}

    function Constant(dirx::T,diry::T,dirz::T) where{T<:Real}
        return new{T}(Vec3(dirx,diry,dirz))
    end

    function Constant(direction::Vec3{T}) where {T<:Real}
        return new{T}(direction)
    end

    function Constant(::Type{T} = Float64) where {T<:Real}
        direction = unitZ3(T)
        return new{T}(direction)
    end
end

Base.length(d::Constant) = 1
Emitters.generate(d::Constant, ::Integer) = d.direction

"""
    RectGrid{T} <: AbstractDirectionDistribution{T}

Encapsulates a single ray direction, where the default direction is unitZ3 [0, 0, 1].

```julia
Constant(direction::Vec3{T}) where {T<:Real}
Constant(::Type{T} = Float64) where {T<:Real}
```
"""
struct RectGrid{T} <: AbstractDirectionDistribution{T}
    direction::Vec3{T}
    halfangleu::T
    halfanglev::T
    nraysu::Integer
    nraysv::Integer
    uvec::Vec3{T}
    vvec::Vec3{T}

    function RectGrid(direction::Vec3{T}, halfangleu::T, halfanglev::T, numraysu::Integer, numraysv::Integer) where {T<:Real}
        (uvec, vvec) = get_orthogonal_vectors(direction)
        return new{T}(direction, halfangleu, halfanglev, numraysu, numraysv, uvec, vvec)
    end

    function RectGrid(halfangleu::T, halfanglev::T, numraysu::Integer, numraysv::Integer) where {T<:Real}
        direction, uvec, vvec = (unitZ3(T), unitX3(T), unitY3(T))
        return new{T}(direction, halfangleu, halfanglev, numraysu, numraysv, uvec, vvec)
    end
end

Base.length(d::RectGrid) = d.nraysu * d.nraysv
function Emitters.generate(d::RectGrid{T}, n::Integer) where {T<:Real}
    direction = d.direction
    uvec = d.uvec
    vvec = d.vvec

    # distributing evenly across the area of the rectangle which subtends the given angle (*not* evenly across the angle range)
    dindex = mod(n, d.nraysu * d.nraysv)
    v = d.nraysv == 1 ? zero(T) : 2 * Integer(floor(dindex / d.nraysu)) / (d.nraysv - 1) - 1.0
    u = d.nraysu == 1 ? zero(T) : 2 * mod(dindex, d.nraysu) / (d.nraysu - 1) - 1.0
    θu = atan(u * tan(d.halfangleu) / 2) * d.halfangleu
    θv = atan(v * tan(d.halfanglev) / 2) * d.halfanglev
    dir = cos(θv) * (cos(θu) * direction + sin(θu) * uvec) + sin(θv) * vvec
    return dir
end

"""
    UniformCone{T} <: AbstractDirectionDistribution{T}

Encapsulates `numsamples` rays sampled uniformly from a cone with max angle θmax.

```julia
UniformCone(direction::Vec3{T}, θmax::T, numsamples::Integer) where {T<:Real}
UniformCone(θmax::T, numsamples::Integer) where {T<:Real}
```
"""
struct UniformCone{T} <: AbstractDirectionDistribution{T}
    direction::Vec3{T}
    θmax::T
    numsamples::Integer
    uvec::Vec3{T}
    vvec::Vec3{T}

    function UniformCone(direction::Vec3{T}, θmax::T, numsamples::Integer) where {T<:Real}
        (uvec, vvec) = get_orthogonal_vectors(direction)
        return new{T}(direction, θmax, numsamples, uvec, vvec)
    end

    function UniformCone(θmax::T, numsamples::Integer) where {T<:Real}
        direction, uvec, vvec = (unitZ3(T), unitX3(T), unitY3(T))
        return new{T}(direction, θmax, numsamples, uvec, vvec)
    end
end

Base.length(d::UniformCone) = d.numsamples
function Emitters.generate(d::UniformCone{T}, ::Integer) where {T<:Real}
    direction = d.direction
    θmax = d.θmax
    uvec = d.uvec
    vvec = d.vvec

    ϕ = rand(T) * 2π
    θ = acos(clamp(one(T) + rand(T) * (cos(θmax) - 1), -one(T), one(T)))
    return normalize(sin(θ) * (cos(ϕ) * uvec + sin(ϕ) * vvec) + cos(θ) * direction)
end

"""
    HexapolarCone{T} <: AbstractDirectionDistribution{T}

Rays are generated by sampling a cone with θmax angle in an hexapolar fashion. The number of rays depends on the requested rings and is computed using the following formula:
`1 + round(Integer, (nrings * (nrings + 1) / 2) * 6)`

```julia
HexapolarCone(direction::Vec3{T}, θmax::T, nrings::Integer) where {T<:Real}
HexapolarCone(θmax::T, nrings::Integer = 3) where {T<:Real}
```
"""
struct HexapolarCone{T} <: AbstractDirectionDistribution{T}
    direction::Vec3{T}
    θmax::T
    nrings::Integer
    uvec::Vec3{T}
    vvec::Vec3{T}

    function HexapolarCone(direction::Vec3{T}, θmax::T, nrings::Integer) where {T<:Real}
        (uvec, vvec) = get_orthogonal_vectors(direction)
        return new{T}(direction, θmax, nrings, uvec, vvec)
    end

    # assume canonical directions
    function HexapolarCone(θmax::T, nrings::Integer = 3) where {T<:Real}
        direction, uvec, vvec = (unitZ3(T), unitX3(T), unitY3(T))
        return new{T}(direction, θmax, nrings, uvec, vvec)
    end
end

Base.length(d::HexapolarCone) = 1 + round(Integer, (d.nrings * (d.nrings + 1) / 2) * 6)
function Emitters.generate(d::HexapolarCone{T}, n::Integer) where {T<:Real}
    dir = d.direction
    θmax = d.θmax
    uvec = d.uvec
    vvec = d.vvec

    n = mod(n, length(d))
    if n == 0
        return normalize(dir)
    else
        t = 1
        ringi = 1
        for i in 1:(d.nrings)
            t += 6 * i
            if n < t
                ringi = i
                break
            end
        end
        ρ = ringi / d.nrings
        pind = n - (t - 6 * ringi)

        ϕ = (pind / (6 * ringi)) * 2π
        # elevation calculated as ring fraction multipled by max angle
        θ = acos(clamp(one(T) +  (cos(ρ * θmax) - 1), -one(T), one(T)))
        return normalize(sin(θ) * (cos(ϕ) * uvec + sin(ϕ) * vvec) + cos(θ) * dir)
    end
end

end # module Directions


text:
module Datasets

using Reexport

@reexport using GraphMLDatasets

end


text:
using JSON
using Spark

function start(worker::Worker)
    server = listen(IPv4(0), worker.port)
    println("Starting worker")
    while worker.active
        sock = accept(server)
        @async while worker.active
            try
                line = readline(sock)
                result = handle(worker, line)
                println(sock, json({"result" => result}))
            catch e
                showerror(STDERR, e)
                break
            end
        end
    end
end

function handle(worker::Worker, line::ASCIIString)
    # General format of a message:
    # {:call => "funcname", :args => {anything}}
    # this is dispatched to any function call - fine since we're assuming
    # a private non-adversarial network.
    msg = JSON.parse(strip(line))
    r = false
    if "call" in keys(msg)
        r = eval(Expr(:call, symbol(msg["call"]), worker, msg["args"]))
    end
    return r
end


text:
using BibInternal
using Test

# TODO: make error testing for new version
@test true


text:
# This file is a part of InverseFunctions.jl, licensed under the MIT License (MIT).

using Test
using InverseFunctions

foo(x) = inv(exp(-x) + 1)
inv_foo(y) = log(y / (1 - y))

InverseFunctions.inverse(::typeof(foo)) = inv_foo
InverseFunctions.inverse(::typeof(inv_foo)) = foo

struct Bar{MT<:AbstractMatrix}
    A::MT
end

(f::Bar)(x) = f.A * x
InverseFunctions.inverse(f::Bar) = Bar(inv(f.A))

@testset "inverse" begin
    f_without_inverse(x) = 1
    @test inverse(f_without_inverse) isa NoInverse
    @test_throws ErrorException inverse(f_without_inverse)(42)
    @test inverse(inverse(f_without_inverse)) === f_without_inverse

    for f in (f_without_inverse ∘ exp, exp ∘ f_without_inverse, Base.Fix1(broadcast, f_without_inverse), Base.Fix1(map, f_without_inverse))
        @test inverse(f) == NoInverse(f)
        @test inverse(inverse(f)) == f
    end

    InverseFunctions.test_inverse(inverse, log, compare = ===)

    InverseFunctions.test_inverse(!, false)

    x = rand()
    for f in (
            foo, inv_foo, log, log2, log10, log1p, sqrt,
            Base.Fix2(^, rand()), Base.Fix2(^, rand([-10:-1; 1:10])), Base.Fix1(^, rand()), Base.Fix1(log, rand()), Base.Fix2(log, rand()),
        )
        InverseFunctions.test_inverse(f, x)
    end
    for f in (
            +, -, exp, exp2, exp10, expm1, cbrt, deg2rad, rad2deg, conj,
            Base.Fix1(+, rand()), Base.Fix2(+, rand()), Base.Fix1(-, rand()), Base.Fix2(-, rand()),
            Base.Fix1(*, rand()), Base.Fix2(*, rand()), Base.Fix1(/, rand()), Base.Fix2(/, rand()), Base.Fix1(\, rand()), Base.Fix2(\, rand()),
            Base.Fix2(^, rand(-11:2:11)),
        )
        InverseFunctions.test_inverse(f, x)
        InverseFunctions.test_inverse(f, -x)
    end
    InverseFunctions.test_inverse(conj, 2 - 3im)

    # ensure that inverses have domains compatible with original functions
    @test_throws DomainError inverse(Base.Fix1(*, 0))
    @test_throws DomainError inverse(Base.Fix2(^, 0))
    @test_throws DomainError inverse(Base.Fix1(log, 2))(-5)
    InverseFunctions.test_inverse(Base.Fix1(log, 2), -5 + 0im)
    @test_throws DomainError inverse(Base.Fix2(^, 0.5))(-5)
    InverseFunctions.test_inverse(Base.Fix2(^, 0.5), -5 + 0im)

    A = rand(5, 5)
    for f in (
            identity, inv, adjoint, transpose,
            log, sqrt, +, -, exp,
            Base.Fix1(+, rand(5, 5)), Base.Fix2(+, rand(5, 5)), Base.Fix1(-, rand(5, 5)), Base.Fix2(-, rand(5, 5)),
            Base.Fix1(*, rand()), Base.Fix2(*, rand()), Base.Fix1(*, rand(5, 5)), Base.Fix2(*, rand(5, 5)),
            Base.Fix2(/, rand()), Base.Fix1(/, rand(5, 5)), Base.Fix2(/, rand(5, 5)),
            Base.Fix1(\, rand()), Base.Fix1(\, rand(5, 5)), Base.Fix2(\, rand(5, 5)),
        )
        if f != log || VERSION >= v"1.6"
            # exp(log(A::AbstractMatrix)) ≈ A is broken on at least Julia v1.0
            InverseFunctions.test_inverse(f, A)
        end
    end

    X = rand(5)
    for f in (Base.Fix1(broadcast, foo), Base.Fix1(map, foo))
        for x in (x, fill(x, 3), X)
            InverseFunctions.test_inverse(f, x)
        end
    end

    InverseFunctions.test_inverse(Bar(rand(3,3)), rand(3))

    @static if VERSION >= v"1.6"
        InverseFunctions.test_inverse(log ∘ foo, x)
    end
end


text:
module FullRegisterGate

export expandGateToFullRegister

# expandGateToFullRegister expands the given gate with optional control qubits to entire quantum register with the given size.
function expandGateToFullRegister(register_size::Integer,
	gate::AbstractMatrix{Complex{Float64}},
	gate_lowest_index::Integer,
	control_bit_indexes::AbstractArray{Int64,1} = Array{Int64,1}([])
	)::AbstractMatrix{Complex{Float64}}

	small_gate_size = size(gate)[1]
	small_gate_qubit_count = Int(log2(small_gate_size))
	gate_bitmask = gate_full_register_bitmask(small_gate_qubit_count, gate_lowest_index)
	if gate_bitmask == 0
		error("gate_bitmask cannot be 0")
	end
	n = 2 ^ register_size
	big_gate = zeros(Complex{Float64}, n, n)
	control_bitmask = control_full_register_bitmask(control_bit_indexes)
	for big_gate_column_index ∈ 0:n-1
		if !all_control_bits_set(big_gate_column_index, control_bitmask)
			big_gate[big_gate_column_index+1, big_gate_column_index+1] = Complex(1)
			continue
		end
		target_bits = (big_gate_column_index & gate_bitmask) >>> gate_lowest_index
		output_state = gate[:,target_bits+1] # Selecting a column here yields the result of matrix-vector multiplication.
		for state_index ∈ 0:small_gate_size-1
			big_gate_row_index = (big_gate_column_index & ~gate_bitmask) | (state_index << gate_lowest_index)
			big_gate[big_gate_row_index+1, big_gate_column_index+1] = Complex(output_state[state_index+1])
		end
	end
	return big_gate
end

function gate_full_register_bitmask(gate_qubit_count::Integer, gate_lowest_index::Integer)::UInt64
	bitmask = zero(UInt64)
	for _ ∈ 1:gate_qubit_count
		bitmask <<= 1
		bitmask |= 1
	end
	return bitmask << gate_lowest_index
end

function control_full_register_bitmask(control_bit_indexes::AbstractArray{Int64,1})::UInt64
	bitmask = zero(UInt64)
	for i ∈ control_bit_indexes
		bitmask |= one(UInt64) << i
	end
	return bitmask
end

@inline function all_control_bits_set(i::Integer, control_bitmask::Integer)
	return i & control_bitmask == control_bitmask
end

end # module


text:
# ---
# title: Row/Column-wise Distances
# author: Johnny Chen
# date: 2021-5-17
# ---

# This example shows how `LazyDistances.PairwiseDistance` can be used to get a better
# `Distances.pairwise` performance.
#
# !!! warning
#     The benchmark result `@btime` is not generated automatically so it may be outdated.
#     If so please file an issue/PR to update it.
#

using Distances, LazyDistances
using LinearAlgebra
using Random
Random.seed!(0)

# ## pairwise distance

X = rand(1:10, 100, 2);
Y = rand(1:10, 100, 3);

dist = pairwise(Euclidean(), X, Y)

# `Distances.pairwise` explicitly generates a large array, which does not perform well for massive
# scale data. [`PairwiseDistance`](@ref) mimics the pairwise operation in a lazy and optionally
# cached manner, so that you don't need to worry about memory allocation and the overhead due to GC.
# Unlike `pairwise`, `PairwiseDistance` does not assume the data to be row or column shaped; instead
# you need to manually specify how data pair is constructed:

lazy_dist = let X=X, Y=Y
    ## defines:
    ##   ĩ := σ(i)
    ##   j̃ := σ(j)
    ## where `i` and `j` are `CartesianIndex`.
    ## We don't need to operate on index `i` and `j` here, so use `identity`.
    index_map = (
        identity,
        identity
    )
    ## defines:
    ##   p = A[ĩ] := @view X[:, ĩ]
    ##   q = B[j̃] := @view Y[:, j̃]
    ## Here the input A is axes(X, 2) and is not used.
    getindex_op = (
        (ax, i) -> view(X, :, i),
        (ax, j) -> view(Y, :, j)
    )
    ## defines:
    ##   out[i, j] := d(p, q)
    d = Euclidean()
    PairwiseDistance(index_map, getindex_op, d, (axes(X, 2), axes(Y, 2)))
end

#-
dist == lazy_dist

# The actual computation is not done until we retrieve the result:
#
# ```julia
# julia> @btime getindex($dist, 2, 3);
#  1.706 ns (0 allocations: 0 bytes)
#
# julia> @btime getindex($lazy_dist, 2, 3)
#  27.086 ns (0 allocations: 0 bytes)
# ```
#
# To properly benchmark the performance, let's do some simple statistics, e.g, `sum`, over the
# result.
#
# ```julia
# X = rand(1:10, 100, 200);
# Y = rand(1:10, 100, 300);
#
# function naive_pairwise_sum(d, X, Y)
#     return sum(pairwise(d, X, Y))
# end
#
# function lazy_pairwise_sum(d, X, Y)
#     getindex_op = (
#         (ax, i) -> view(X, :, i),
#         (ax, j) -> view(Y, :, j)
#     )
#     dist = PairwiseDistance(identity, getindex_op, d, map(A->axes(A, 2), (X, Y)))
#     sum(dist)
# end
#
# naive_pairwise_sum(Euclidean(), X, Y) ≈ lazy_pairwise_sum(Euclidean(), X, Y) # true
#
# @btime naive_pairwise_sum(Euclidean(), $X, $Y); # 4.513 ms (5 allocations: 473.23 KiB)
# @btime lazy_pairwise_sum(Euclidean(), $X, $Y); # 2.265 ms (13 allocations: 352 bytes)
# ```
#
# Great that we have successfully removed the memory allocation overhead and it does improve the
# performance. Of course, if you need to reuse the data repeatedly, you still need to `collect` it
# so that you don't suffer from the repeat calculation.

# ## columnwise distance
#
# Strictly speaking, `Distances.columnwise` distance is just a simple map function. Again, in
# contrast, `PairwiseDistance` does not assume the data layout, and you have to specify it

X = rand(1:10, 100, 2);
Y = rand(1:10, 100, 2);
dist = colwise(Euclidean(), X, Y)

lazy_dist = let X=X, Y=Y
    getindex_op = (
        (ax, i) -> view(X, :, i),
        (ax, j) -> view(Y, :, j)
    )
    PairwiseDistance(identity, getindex_op, Euclidean(), (axes(X, 2), axes(Y, 2)))
end

diag(lazy_dist) == dist

# Let's also benchmark the results:
#
# ```julia
# X = rand(1:10, 100, 200);
# Y = rand(1:10, 100, 200);
#
# function naive_colwise_sum(d, X, Y)
#     return sum(colwise(d, X, Y))
# end
#
# function lazy_colwise_sum(d, X, Y)
#     getindex_op = (
#         (ax, i) -> view(X, :, i),
#         (ax, j) -> view(Y, :, j)
#     )
#     dist = PairwiseDistance(identity, getindex_op, d, map(A->axes(A, 2), (X, Y)))
#     sum(diag(dist))
# end
#
# naive_colwise_sum(Euclidean(), X, Y) ≈ lazy_colwise_sum(Euclidean(), X, Y) # true
#
# @btime naive_colwise_sum(Euclidean(), $X, $Y); # 6.002 μs (1 allocation: 1.77 KiB)
# @btime lazy_colwise_sum(Euclidean(), $X, $Y); # 8.166 μs (14 allocations: 2.11 KiB)
# ```
#
# As you can see, this time our lazy version becomes slower than the colwise version. This is
# because the `diag` operation itself gives some overhead. `PairwiseDistance` itself also
# contributes some of the overhead.


text:
using BIDSTools
using Test
using DataStructures
using DataFrames

@testset "BIDSTools.jl" begin

    include("File.jl")
    include("Session.jl")
    include("Subject.jl")
    include("Layout.jl")

end


text:
# Overlap Iterator
# ================

struct OverlapIterator{Sa,Sb,F,G}
    intervals_a::Sa
    intervals_b::Sb
    isless::F
    filter::G
end

function Base.eltype(::Type{OverlapIterator{Sa,Sb,F,G}}) where {Sa,Sb,F,G}
    return Tuple{Interval{metadatatype(Sa)},Interval{metadatatype(Sb)}}
end

function Base.IteratorSize(::Type{OverlapIterator{Sa,Sb,F,G}}) where {Sa,Sb,F,G}
    return Base.SizeUnknown()
end

"""
    eachoverlap(intervals_a, intervals_b, [seqname_isless=Base.isless])

Create an iterator of overlapping intervals between `intervals_a` and `intervals_b`.

This function assumes elements of `intervals_a` and `intervals_b` are sorted by
its sequence name and left position.  If the element type is not a subtype of
`GenomicFeatures.Interval`, elements are converted to `Interval` objects.

The third optional argument is a function that defines the order of sequence
names. The default function is `Base.isless`, which is the lexicographical
order.
"""
function eachoverlap(intervals_a, intervals_b, seqname_isless=Base.isless; filter=true_cmp)
    return OverlapIterator(intervals_a, intervals_b, seqname_isless, filter)
end

struct OverlapIteratorState{Sa,Sb,Ta,Tb}
    next_a::Sa
    next_b::Sb
    queue::Queue{Interval{Tb}}
    queue_index::Int
end

function OverlapIteratorState(
        Ta::Type, Tb::Type, next_a::Sa, next_b::Sb, queue::Queue, queue_index::Int) where {Sa, Sb}
    return OverlapIteratorState{Sa,Sb,Ta,Tb}(next_a, next_b, queue, queue_index)
end

function OverlapIteratorState(Ta::Type, Tb::Type, next_a::Sa, next_b::Sb) where {Sa, Sb}
    queue = Queue{Interval{Tb}}()
    return OverlapIteratorState{Sa,Sb,Ta,Tb}(next_a, next_b, queue, 1)
end

function Base.iterate(iter::OverlapIterator)
    next_a = iterate(iter.intervals_a)
    next_b = iterate(iter.intervals_b)

    Ta = metadatatype(iter.intervals_a)
    Tb = metadatatype(iter.intervals_b)
    state = OverlapIteratorState(Ta, Tb, next_a, next_b)

    return iterate(iter, state)
end

# check i1 and i2 are ordered
function check_ordered(i1, i2, compare_func)
    if !isordered(i1, i2, compare_func)
        error("intervals are not sorted")
    end
    return nothing
end

function Base.iterate(iter::OverlapIterator, state::OverlapIteratorState{Sa,Sb,Ta,Tb}) where {Sa,Sb,Ta,Tb}
    next_a      = state.next_a
    next_b      = state.next_b
    queue       = state.queue
    queue_index = state.queue_index

    if next_a === nothing
        return nothing
    end

    entry_a, state_a = next_a
    interval_a = convert(Interval{Ta}, entry_a)

    while true
        if queue_index > lastindex(state.queue)
            # end of queue: add more to queue, or advance a
            if next_b === nothing
                next_a = iterate(iter.intervals_a, state_a)
                if next_a === nothing
                    return break
                end

                entry_a, state_a = next_a
                next_interval_a = convert(Interval{Ta}, entry_a)
                check_ordered(interval_a, next_interval_a, iter.isless)
                interval_a = next_interval_a
                queue_index = firstindex(state.queue)
            else
                entry_b, state_b = next_b
                interval_b = convert(Interval{Tb}, entry_b)
                if !isempty(queue)
                    check_ordered(queue[end], interval_b, iter.isless)
                end
                push!(queue, interval_b)
                next_b = iterate(iter.intervals_b, state_b)
            end
        else
            entry_a, state_a = next_a
            interval_a = convert(Interval{Ta}, entry_a)
            interval_b = queue[queue_index]
            c = compare_overlap(interval_a, interval_b, iter.isless)
            queue_index += 1

            if c < 0
                # No more possible intersections with interval_a, advance
                next_a = iterate(iter.intervals_a, state_a)
                if next_a === nothing
                    break
                end
                entry_a, state_a = next_a
                next_interval_a = convert(Interval{Ta}, entry_a)

                check_ordered(interval_a, next_interval_a, iter.isless)
                interval_a = next_interval_a
                queue_index = firstindex(state.queue)
            elseif c == 0
                if iter.filter(interval_a, interval_b)
                    return ((interval_a, interval_b),
                        OverlapIteratorState(Ta, Tb, next_a, next_b, queue, queue_index))
                end
            else
                if queue_index == firstindex(queue) + 1
                    # noting else can intersect front of queue
                    popfirst!(queue)
                end
            end
        end
    end

    # no more intersections found
    return nothing
end

# Return:
#   -1 when `i1` precedes `i2`,
#   0 when `i1` overlaps with `i2`, and
#   +1 when `i1` follows `i2`.
function compare_overlap(i1::Interval, i2::Interval, isless::Function)
    if isless(i1.seqname, i2.seqname)::Bool
        return -1
    elseif isless(i2.seqname, i1.seqname)::Bool
        return +1
    else  # i1.seqname == i2.seqname
        if i1.last < i2.first
            return -1
        elseif i1.first > i2.last
            return +1
        else
            return 0
        end
    end
end

# Faster comparison for `Base.isless`.  Note that `Base.isless` must be
# consistent wtih `Base.cmp` to work correctly.
function compare_overlap(i1::Interval, i2::Interval, ::typeof(Base.isless))
    c = cmp(i1.seqname, i2.seqname)
    if c != 0
        return c
    end
    if i1.last < i2.first
        return -1
    elseif i1.first > i2.last
        return +1
    else
        return 0
    end
end