text:
module Netlist

using QSpice.State, QSpice.Gates, QSpice.Parser

export QuantumCircuit, readnetlist, simulate

# Any subtype of this may be a node in the DAG representing the
# quantum circuit
abstract QuantumCircuitNode

type GateNode <: QuantumCircuitNode
    fn::Function
    args::Vector{Any}

    qinputs::Vector{QuantumState}
    binputs::Vector{Vector{Int}}

    quantumedges::Vector{Int}
    bitedges::Vector{Int}
end

type QuantumNode <: QuantumCircuitNode
    state::QuantumState
    edges::Vector{Int}
end

type BitNode <: QuantumCircuitNode
    bits::Vector{Int}
    edges::Vector{Int}
end

type QuantumCircuit
    gatenodes::Dict{Int, GateNode}
    qnodes::Dict{Int, QuantumNode}
    bnodes::Dict{Int, BitNode}
end

# The list of nodes that don't need to wait for additional inputs
# in case of a quantum state change
const QUANTUM_PASSTHROUGH = [Gates.identity, hadamard, not, cnot,
                             ccnot, swap, cswap, sqrtswap, qft,
                             paulix, pauliy, pauliz, phaseshift,
                             measure, partialmeasure, probe]

# The list of nodes that don't need to wait for additional inputs
# in case of classical bit state change
const BIT_PASSTHROUGH = [probe]

function readnetlist(s)
    function addnotfound(gatenodes, index, gate)
        if !haskey(gatenodes, index)
            gatenodes[index] = GateNode(gate.fn, gate.args, [], [], [], [])
        end
    end

    gates, qstates, bstates = parsenetlist(s)

    gatenodes = Dict{Int, GateNode}()

    # First initialize quantum and classical states
    qnodes = [i[1] => QuantumNode(i[2], []) for i in qstates]
    bnodes = [i[1] => BitNode(i[2], []) for i in bstates]

    # Construct the DAG by iterating the gates and adding all necessary
    # edges (representing a dependency on either a quantum or classical state change)
    # to all of the nodes
    for (index, gate) in gates
        addnotfound(gatenodes, index, gate)

        for qi in gate.quantuminputs
            if haskey(qnodes, qi)
                push!(qnodes[qi].edges, index)
                continue
            end
            addnotfound(gatenodes, qi, gates[qi])
            push!(gatenodes[qi].quantumedges, index)
        end

        for bi in gate.bitinputs
            if haskey(bnodes, bi)
                push!(bnodes[bi].edges, index)
                continue
            end
            addnotfound(gatenodes, bi, gates[bi])
            push!(gatenodes[bi].bitedges, index)
        end
    end

    return QuantumCircuit(gatenodes, qnodes, bnodes)
end

# Broadcasts a quantum state change to all nodes connected by a quantum edge
function broadcast(netlist::QuantumCircuit, gatenode::GateNode, qstate::QuantumState)
    for qedge in gatenode.quantumedges
        accept(netlist, netlist.gatenodes[qedge], qstate)
    end
end

# Broadcasts a classical bit state change to all nodes connected by a bit edge
function broadcast(netlist::QuantumCircuit, gatenode::GateNode, bstate::Vector{Int})
    for bedge in gatenode.bitedges
        accept(netlist, netlist.gatenodes[bedge], bstate)
    end
end

# Broadcasts both a quantum and classical state change
function broadcast(netlist::QuantumCircuit, gatenode::GateNode,
                   states::Tuple{QuantumState, Vector{Int}})
    broadcast(netlist, gatenode, states[1])
    broadcast(netlist, gatenode, states[2])
end

function broadcast(netlist::QuantumCircuit, gatenode::GateNode, v::Void)
    # Ordinarily only probes should call this function
end

# Used by gates that have multiple inputs. It checks if the gate is ready to execute,
# and if it is, performs the operation.
# TODO(gustorn): this only supports unordered arguments right now, which is a pretty big limiation
# (and a wonder that it even works in case of superposition). Maybe make edges between inputs/outputs
# and not gates?
function executeready(netlist::QuantumCircuit, gatenode::GateNode)
    if (gatenode.fn == superposition || gatenode.fn == join) && length(gatenode.qinputs) == 2
        broadcast(netlist, gatenode, gatenode.fn(gatenode.qinputs..., gatenode.args...))
        empty!(gatenode.qinputs)
    elseif gatenode.fn == choose1 && !isempty(gatenode.qinputs) && !isempty(gatenode.binputs)
        broadcast(netlist, gatenode, gatenode.fn(gatenode.qinputs...,
                                                 gatenode.binputs...,
                                                 convert(Vector{Vector{Tuple{Function, Vector{Any}}}}, gatenode.args)))
        empty!(gatenode.qinputs)
        empty!(gatenode.binputs)
    end
end

function accept(netlist::QuantumCircuit, gatenode::GateNode, qstate::QuantumState)
    if gatenode.fn in QUANTUM_PASSTHROUGH
        broadcast(netlist, gatenode, gatenode.fn(qstate, gatenode.args...))
        return
    end
    push!(gatenode.qinputs, qstate)
    executeready(netlist, gatenode)
end

function accept(netlist::QuantumCircuit, gatenode::GateNode, bstate::Vector{Int})
    if gatenode.fn in BIT_PASSTHROUGH
        broadcast(netlist, gatenode, gatenode.fn(bstate, gatenode.args...))
        return
    end
    push!(gatenode.binputs, bstate)
    executeready(netlist, gatenode)
end

# Simulates the netlist by flushing all initial qubits and classical bitstates.
# If the DAG is correct, this will execute every reachable gate in the circuit
function simulate(netlist::QuantumCircuit)
    for qn in values(netlist.qnodes)
        for edge in qn.edges
            accept(netlist, netlist.gatenodes[edge], qn.state)
        end
    end
    for bn in values(netlist.bnodes)
        for edge in bn.edges
            accept(netlist, netlist.gatenodes[edge], bn.bits)
        end
    end
end

end

max_stars_repo_name: leios/quantum

max_issues_repo_name: leios/quantum

max_forks_repo_name: leios/quantum


text:
@naiveVariationalRule(:node_type     => Logit,
                      :outbound_type => Message{Bernoulli},
                      :inbound_types => (Nothing, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VBLogitOut)

@naiveVariationalRule(:node_type     => Logit,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (ProbabilityDistribution, Nothing, ProbabilityDistribution),
                      :name          => VBLogitIn1)

@naiveVariationalRule(:node_type     => Logit,
                      :outbound_type => Message{Function},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, Nothing),
                      :name          => VBLogitXi)

ext: jl


text:
module MathleticsFiles

using  DataFrames
using  Downloads: Downloads
using  Logging
using  Pkg.Artifacts: artifact_exists, artifact_hash, @artifact_str, bind_artifact!, create_artifact
using  Taro: Taro
using  ZipFile: ZipFile

export dataset

const URL = "https://excelwithwayne.com/wp-content/uploads/2019/05/MathleticsFiles.zip"
const ARTIFACT_TOML = joinpath(@__DIR__, "..", "Artifacts.toml")

function __init__()
    @info "Initializing Taro"
    Taro.init()
    try
        download_artifact()
    catch e
        @warn "Unable to download requisite artifacts; functionality will be limited until you run MathleticsFiles.download_artifact()"
    end
end

function extractall(zippath::AbstractString, dirpath::AbstractString)
    @debug "Extracting all contents of $zippath to $dirpath."
    open(zippath) do io
        reader = ZipFile.Reader(io)
        for file in reader.files
            outpath = joinpath(dirpath, file.name)
            @debug "Extracting $(file.name) to $outpath."
            open(outpath, write=true) do out
                write(out, read(file))
            end;
        end
    end;
end

function download_artifact(; sourceurl::AbstractString=URL)
    artifacthash = artifact_hash("mathletics_files", ARTIFACT_TOML)

    if isnothing(artifacthash) || !artifact_exists(artifacthash)
        @info "Generating MathleticsFiles artifact..."
        artifacthash = create_artifact() do artifactpath
            extractall(Downloads.download(sourceurl), artifactpath)
        end
        bind_artifact!(ARTIFACT_TOML, "mathletics_files", artifacthash, force=true)
        @info "MathleticsFiles artifact generated."
    end
end

function astemp(f::Function, path::AbstractString, mode::Integer)
    temppath = tempname()
    return try
        cp(path, temppath)
        chmod(temppath, mode)
        f(temppath)
    finally
        rm(temppath)
    end
end

function openartifact(f::Function, filename::AbstractString)
    # Taro.readxl requres write access; copy the path to a temporary location with write access.
    astemp(f, joinpath(artifact"mathletics_files", filename), 0o664)
end

"""
Retrieves a dataset by file name, sheetname, and range
"""
function dataset(filename::AbstractString, sheetname::AbstractString, range::AbstractString)::DataFrame
    openartifact(filename) do temppath
        return DataFrame(Taro.readxl(temppath, sheetname, range))
    end
end

"""
Retrieves a dataset by file name and range
"""
function dataset(filename::AbstractString, range::AbstractString)::DataFrame
    openartifact(filename) do temppath
        return DataFrame(Taro.readxl(temppath, range))
    end
end

"""
Retrieves a dataset (as a DataFrame) by alias.
"""
function dataset(alias::AbstractString)::DataFrame
    lookup = Dict{String,Tuple{Vararg{String}}}(
        "nfl_team_totals" => ("nflregression.xls", "data", "A5:M133"),
        "rushing" => ("firstdown.xls", "B4:E22"),
        "passing" => ("firstdown.xls", "B24:E42"),
    )
    args = lookup[alias]
    return dataset(args...)
end

end # module


text:
module BMath
    include("BPoint.jl")
    include("BPlot.jl")
end

using Main.BMath
using Main.BMath.BPlot

sets = []
xvals = []

for x in 1.0:1000.0
    push!(xvals, x)
    push!(sets, BPoint(x, [0.0]) + BPoint(x, [pi/2]))
end

plotb(xvals, sets)


text:
module PkgUnboundArgs

# Putting it in a submodule to test that `recursive=true` is used.
module M25341
_totuple(::Type{Tuple{Vararg{E}}}, itr, s...) where {E} = E
end
# `_totuple` is taken from
# https://github.com/JuliaLang/julia/blob/48634f9f8669e1dc1be0a1589cd5be880c04055a/test/ambiguous.jl#L257-L259

end  # module


text:
# ---
# title: 567. Permutation in String
# id: problem567
# author: zhwang
# date: 2022-02-17
# difficulty: Medium
# categories: Two Pointers, Sliding Window
# link: <https://leetcode.com/problems/permutation-in-string/description/>
# hidden: true
# ---
#
# Given two strings **s1** and **s2** , write a function to return true if
# **s2** contains the permutation of **s1**. In other words, one of the first
# string's permutations is the **substring** of the second string.
#
#
#
# **Example 1:**
#
#
#
#     Input: s1 = "ab" s2 = "eidbaooo"
#     Output: True
#     Explanation: s2 contains one permutation of s1 ("ba").
#
#
# **Example 2:**
#
#
#
#     Input: s1= "ab" s2 = "eidboaoo"
#     Output: False
#
#
#
#
# **Constraints:**
#
#   * The input strings only contain lower case letters.
#   * The length of both given strings is in range [1, 10,000].
#
#
## @lc code=start
using LeetCode

using DataStructures
function check_inclusion(s1::String, s2::String)::Bool
    equal2c1(c::Accumulator) = all(c[i]==c1[i] for i in keys(c1))
    c1, n1, n2 = counter(s1), length(s1), length(s2)
    n2 < n1 && return false
    c2 = counter(@view(s2[1:n1]))
    equal2c1(c2) && return true
    for (i, j) in zip(@view(s2[1:(n2 - n1)]), @view(s2[n1 + 1: n2]))
        c2[i] -= 1
        c2[j] += 1
        equal2c1(c2) && return true
    end
    false
end
## @lc code=end

## @lc test=start
@testset "567.permutation-in-string.jl" begin
    @test check_inclusion("abc", "ccccbbbbaaaa") == false
    @test check_inclusion("ab", "eidbaooo") == true
    @test check_inclusion("ab", "eidboaoo") == false
    @test check_inclusion("hello", "ooolleoooleh") == false
    @test check_inclusion("adc", "dcda") == true
end
## @lc test=end

## @lc info=start
# link: [solution to 567](https://leetcode-cn.com/problems/permutation-in-string/solution/567-zi-fu-chuan-de-pai-lie-by-rex-rs-ac0v/)
## @lc info=end


text:
"""
every model interface must implement a `fit` method of the form
`fit(model, verb::Integer, training_args...) -> fitresult, cache, report`
"""
function fit end

# fallback for static transformations
fit(::Static, ::Integer, a...) = (nothing, nothing, nothing)

# fallbacks for supervised models that don't support sample weights:
fit(m::Supervised, verb::Integer, X, y, w) = fit(m, verb, X, y)

# this operation can be optionally overloaded to provide access to
# fitted parameters (eg, coeficients of linear model):
fitted_params(::Model, fitres) = (fitresult=fitres,)

"""
each model interface may overload the `update` refitting method
"""
update(m::Model, verb::Integer, fitres, cache, a...) = fit(m, verb, a...)

"""
each model interface may overload the `update_data` refitting method for online learning
"""
function update_data end

"""
supervised methods must implement the `predict` operation
"""
function predict end

"""
probabilistic supervised models may overload `predict_mean`
"""
function predict_mean end

"""
probabilistic supervised models may overload `predict_mode`
"""
function predict_mode end

"""
probabilistic supervised models may overload `predict_median`
"""
function predict_median end

"""
`JointProbabilistic` supervised models MUST overload `predict_joint`.

`Probabilistic` supervised models MAY overload `predict_joint`.
"""
function predict_joint end

"""
unsupervised methods must implement the `transform` operation
"""
function transform end

"""
unsupervised methods may implement the `inverse_transform` operation
"""
function inverse_transform end

# models can optionally overload these for enable serialization in a
# custom format:
function save end
function restore end

"""
some meta-models may choose to implement the `evaluate` operations
"""
function evaluate end


text:
# ------------------------------------------------------------
# Intervals

type Interval{T<:Real}
	left::T
	right::T

	function Interval(left, right)
		if left < right
			new(left, right)
		else
			throw(ArgumentError("Not an interval"))
		end
	end
end

function Interval(left, right)
	left, right = promote(left, right)
	T = typeof(left)
	Interval{T}(left, right)
end

typealias DaubSupport Interval{Int64}

@inline left(I::Interval) = I.left
@inline right(I::Interval) = I.right
Base.length(I::Interval) = right(I) - left(I)
@compat Base.:+(S::DaubSupport, k::Integer) = DaubSupport(left(S)+k, right(S)+k)

function Base.intersect(I1::Interval, I2::Interval)
	if left(I1) >= right(I2) || left(I2) >= right(I1)
		return Void
	end

	L = max( left(I1), left(I2) )
	R = min( right(I1), right(I2) )

	DaubSupport(L, R)
end

# ------------------------------------------------------------
# Vectors on DaubSupport

@doc """
	dyadic_rationals(I::DaubSupport, R) -> Vector

The dyadic rationals of resolution `R` in the interval `I`.
"""->
function dyadic_rationals(I::DaubSupport, R::Integer)
	R >= 0 || throw(DomainError())
	return left(I):2.0^(-R):right(I)
end

@doc """
	dyadic_rationals(I::DaubSupport, R, level) -> Vector

In a vector of dyadic rationals up to resolution `R` in `I`,
return the **indices** of those at exactly `level`.
"""->
function dyadic_rationals(I::DaubSupport, R::Integer, level::Integer)
	0 <= level <= R || throw(DomainError())

	N = length(I)*2^R + 1

	if level == 0
		step = 2^R
		return 1:step:N
	else
		step = 2^(R-level)
		return 1+step:2*step:N
	end
end

@doc """
	isuniform(x::AbstractVector) -> Bool

Test if the elements in `x` increments with the same amount (with sign).
"""->
function isuniform(x::AbstractVector)
	(Nx = length(x)) > 1 || throw(DomainError())
	Nx == 2 && return true

	diff = x[1] - x[2]
	@inbounds for n in 3:Nx
		if !isapprox(diff, x[n-1] - x[n])
			return false
		end
	end

	return true
end


text:
module Datasets

using Random: AbstractRNG, GLOBAL_RNG, MersenneTwister

abstract type AbstractDataset{D} end

ratio2num(n_data, ratio) = floor(Int, n_data * ratio)

### Synthetic datasets

include("ring.jl")
export RingDataset, TwoDimRingDataset, ThreeDimRingDataset
include("gaussian.jl")
export GaussianDataset

using StatsFuns: logit, logistic

function link_bin(x)
    a, b = zero(eltype(x)), one(eltype(x))
    return logit.((x .- a) ./ (b - a))
end
function invlink_bin(y)
    a, b = zero(eltype(y)), one(eltype(y))
    return (b - a) .* logistic.(y) .+ a
end

### Image datasets

using MLDataUtils: convertlabel, LabelEnc

abstract type ImageDataset{T, D} <: AbstractDataset{D} end

flatten(X) = reshape(X, :, last(size(X)))

function dequantize(rng, x, alpha::T) where {T}
    y = x + rand(rng, T) / 256
    y = clamp(y, 0, 1)
    y = alpha + (1 - 2alpha) * y
    return y
end

link(x) = link_bin(x)
invlink(d::ImageDataset{Val{:true}}, y) = invlink_bin(y)
invlink(d::ImageDataset{Val{:false}}, y) = y

function preprocess(rng, X, is_flatten, alpha, is_link::Bool=false)
    if is_flatten
        X = flatten(X)
    elseif ndims(X) == 3    # e.g. MNIST from MLDatasets has shape of (28, 28, 1, ?)
        X = reshape(X, (Base.front(size(X))..., 1, last(size(X))))
    end
    if !iszero(alpha)
        X = dequantize.(Ref(rng), X, alpha)
    end
    if is_link
        X = link(X)
    end
    return X
end

# TODO: implement flip augmentation
# Ref: https://github.com/gpapamak/maf/blob/master/datasets/cifar10.py#L52-L62
function get_image_data(
    IMAGE,
    n_data::Int,
    n_test::Int,
    is_flatten::Bool,
    alpha::T,
    is_link::Bool,
    K::Int,
    perm::Union{NTuple{3,Int}, NTuple{4,Int}};
    rng::AbstractRNG=GLOBAL_RNG,
) where {T}
    onehot_enc = LabelEnc.OneOfK(K)

    X = permutedims(IMAGE.traintensor(T, 1:n_data), perm)
    y = IMAGE.trainlabels(1:n_data) .+ 1
    X = preprocess(rng, X, is_flatten, alpha, is_link)
    Y = convertlabel(onehot_enc, y)

    Xt = permutedims(IMAGE.testtensor(T, 1:n_test), perm)
    yt = IMAGE.testlabels(1:n_test) .+ 1
    Xt = preprocess(rng, Xt, is_flatten, alpha, is_link)
    Yt = convertlabel(onehot_enc, yt)

    return X, y, Y, Xt, yt, Yt
end

include("mnist.jl")
export MNISTDataset
include("cifar10.jl")
export CIFAR10Dataset
include("features.jl")
export FeatureDataset, get_features_griffiths2011, get_features_xu2019

include("vis.jl")
export n_display, vis
include("utilites.jl")
export Dataset

export Datasets

end # module


text:
using SMTPClient

isdefined(@__MODULE__,:password) || (password = readline())
url = "smtps://smtp.gmail.com:465"

opt = SendOptions(
    isSSL = true,
    username = "akels14@gmail.com",
    passwd = password,verbose=true)

body = IOBuffer(
    """
    From: akels14@gmail.com
    To: graphitewriter@gmail.com
    Subject: Julia Test 2

    Test Message 23
    """
)

from = "akels14@gmail.com"

rcpt = ["graphitewriter@gmail.com"]

resp = send(url, rcpt, from, body, opt)

max_stars_repo_name: akels/DemocracyCredit

max_issues_repo_name: akels/DemocracyCredit

max_forks_repo_name: akels/DemocracyCredit


text:
using Test
using TightBindingLattice

@testset "irrepdatabase" begin
    # D10 group
    @testset "find" begin
        group1 = FiniteGroup([
            1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20;
            2  3  4  5  6  7  8  9  10  1  20  11  12  13  14  15  16  17  18  19;
            3  4  5  6  7  8  9  10  1  2  19  20  11  12  13  14  15  16  17  18;
            4  5  6  7  8  9  10  1  2  3  18  19  20  11  12  13  14  15  16  17;
            5  6  7  8  9  10  1  2  3  4  17  18  19  20  11  12  13  14  15  16;
            6  7  8  9  10  1  2  3  4  5  16  17  18  19  20  11  12  13  14  15;
            7  8  9  10  1  2  3  4  5  6  15  16  17  18  19  20  11  12  13  14;
            8  9  10  1  2  3  4  5  6  7  14  15  16  17  18  19  20  11  12  13;
            9  10  1  2  3  4  5  6  7  8  13  14  15  16  17  18  19  20  11  12;
            10  1  2  3  4  5  6  7  8  9  12  13  14  15  16  17  18  19  20  11;
            11  12  13  14  15  16  17  18  19  20  1  2  3  4  5  6  7  8  9  10;
            12  13  14  15  16  17  18  19  20  11  10  1  2  3  4  5  6  7  8  9;
            13  14  15  16  17  18  19  20  11  12  9  10  1  2  3  4  5  6  7  8;
            14  15  16  17  18  19  20  11  12  13  8  9  10  1  2  3  4  5  6  7;
            15  16  17  18  19  20  11  12  13  14  7  8  9  10  1  2  3  4  5  6;
            16  17  18  19  20  11  12  13  14  15  6  7  8  9  10  1  2  3  4  5;
            17  18  19  20  11  12  13  14  15  16  5  6  7  8  9  10  1  2  3  4;
            18  19  20  11  12  13  14  15  16  17  4  5  6  7  8  9  10  1  2  3;
            19  20  11  12  13  14  15  16  17  18  3  4  5  6  7  8  9  10  1  2;
            20  11  12  13  14  15  16  17  18  19  2  3  4  5  6  7  8  9  10  1;
        ])
        @test isnothing(IrrepDatabase.find(group1))
        group2 = FiniteGroup([1 2; 2 1])
        (irrep_data, elmap) = IrrepDatabase.find(group2)
        @test irrep_data.group == group2
        @test irrep_data.irreps == [[ones(ComplexF64, 1, 1), ones(ComplexF64, 1, 1)],
                                    [ones(ComplexF64, 1, 1), -ones(ComplexF64, 1, 1)]]
        @test irrep_data.conjugacy_classes == [[1], [2]]
        @test irrep_data.character_table == [1 1; 1 -1]
        @test elmap == [1,2]
    end

    @testset "find_irrep" begin
        (irrep, matrep_new, ϕ) = IrrepDatabase.find_irrep([ [1 0; 0 1], [1 0; 0 -1]])
        @test irrep.group == FiniteGroup([1 2; 2 1])
    end
end


text:
#================
Predictive
    Learning
        Models
================#
@doc """
| Model Type | Model Name |\n
|:---------- | ---------- |\n
| Baseline    | MeanBaseline |\n
| Baseline    | ClassBaseline |\n
| Continuous    | LinearRegression |\n
| Continuous    | LinearLeastSquare  |\n
| Categorical    | RandomForestClassifier  |\n
| Categorical    | DecisionTree  |\n
| Tools    | Pipeline  |\n
| Tools    | PowerLog  |
       """
module models
# [deps]
using Random
using DataFrames
using Lathe.stats
# [deps]
# Type heirarchy
abstract type Model <: LatheObject end
abstract type SupervisedModel <: Model end
abstract type ContinuousModel <: SupervisedModel end
abstract type Regressor <: ContinuousModel end
abstract type LinearModel <: Regressor end
abstract type CategoricalModel <: Model end
abstract type BaselineModel <: Model end
abstract type Tool <: Model end
abstract type UnsupervisedModel <: Model end
abstract type Classifier <: CategoricalModel end
#==
Continuous Models
==#
include("models/regressors.jl")
export LinearLeastSquare
export LinearRegression
export LassoRegression
#==
Categorical Models
==#
# Decision Tree Classifier/ RF
include("models/treeclass.jl")
export DecisionTreeClassifier, RandomForestClassifier
# Logistic Regression
include("models/logistic.jl")
export LogisticRegression
#==
Tools
==#
include("models/toolbox.jl")
export Pipeline, PowerLog, MeanBaseline, majClassBaseline
#----------------------------------------------
end


text:
using EchoJulia

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

x = transform(SimradRaw.load(EK60_SAMPLE))

@test length(x["t38"]) == 572

x = transform(SimradRaw.load(EK60_SAMPLE), starttime = 129053535177200000)

@test length(x["t38"]) == 391


text:
using Printf

for n in (0, 5, 50, 9000)
    @printf("%6i → %s\n", n, string(n, base=2))
end

# with pad
println("\nwith pad")
for n in (0, 5, 50, 9000)
    @printf("%6i → %s\n", n, string(n, base=2, pad=20))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    ModelParams

Support for calibrating models.

A `ModelObject` contains a `ParamVector` and the corresponding parameters. It is key to always keep those in sync.

Change:
    Reporting parameters
        make a "tree" of the model structure +++
            each node has name, type, pvec, children
            then one can re-engineer everything else from this (pvectors, child objects, etc)
"""
module ModelParams

import Base.show, Base.isempty, Base.isequal
import Base.append!, Base.length, Base.getindex, Base.values
import Random: AbstractRNG
using ArgCheck, DataStructures, DocStringExtensions, Formatting, Infiltrator, Lazy, Parameters, PrettyTables, Random
using ModelObjectsLH
using EconometricsLH # : RegressionTable, get_all_coeff_se, get_coeff_se_multiple, have_same_regressors

# Transformations
export LinearTransformation, transform_bounds, transform_param, untransform_param

# Parameters
export AbstractParam, Param
export calibrate!, fix!, fix_values!, set_bounds!, set_value!, set_default_value!, update!, validate, default_value, value, is_calibrated

# ParamVector
export ParamVector
export param_exists, make_dict, n_calibrated_params
export param_table, param_value, param_default_value, report_params, retrieve

# PVectorCollection
export PVectorCollection, set_calibration_status_all_params!, set_default_values_all_params!

# Model objects
export check_fixed_params, check_calibrated_params, check_own_fixed_params, check_own_calibrated_params, validate_all_params
export collect_pvectors, compare_params, check_params_match
export find_pvector, find_param, make_guess, perturb_guess_vector, perturb_params, params_equal, validate
export has_pvector, get_pvector, get_switches, param_tables, latex_param_table
export set_values_from_dicts!, sync_own_values!, sync_values!
export BoundedVector, IncreasingVector, values, set_pvector!

# ParamTable
export ParamTable, get_symbol, get_description, get_values, set_row!, latex_param_table

const ValueType = Float64;
# Bounds for parameter transformations into guess vectors.
const TransformationLb = 1;
const TransformationUb = 2;

include("param_table.jl");
include("types.jl");
# General purpose code copied from `CommonLH`
include("helpers.jl");
include("value_vector.jl");
include("guess.jl");

# Parameters
include("bounded_vector.jl");
include("bvector.jl");
# include("calibrated_array.jl");
include("cal_array_param.jl");
include("transformations.jl");
include("parameters.jl");
include("param_vector.jl");

include("m_objects.jl");
include("obj_pvectors.jl");

## ------------  Main user interface functions

"""
	$(SIGNATURES)

Perturb guesses at indices `dIdx` by amount `dGuess`.
Ensure that guesses stay in bounds.
The user typically calls the method that accepts a `ValueVector` instead.
"""
function perturb_guess_vector(m :: ModelObject, g :: Guess{F1},
    guessV :: AbstractVector{F1}, dGuess; dIdx = nothing) where F1

    # vv = make_guess(m);
    if isnothing(dIdx)
        guess2V = guessV .+ dGuess;
    else
        guess2V = copy(guessV);
        guess2V[dIdx] = guessV[dIdx] .+ dGuess;
    end
    enforce_bounds!(m, g, guess2V);
    return guess2V
end

function enforce_bounds!(m :: ModelObject, g :: Guess{F1},
    guessV :: AbstractVector{F1}) where F1
    pvecV = collect_pvectors(m);

    for (_, pvec) in pvecV
        enforce_bounds!(pvec, g, guessV);
    end
end

function enforce_bounds!(pvec :: ParamVector, g :: Guess{F1}, guessV :: AbstractVector{F1}) where F1
    clamp!(guessV, lb(pvec.pTransform), ub(pvec.pTransform));
    # for (j, v) in enumerate(guessV)
    #     if v < lb(pvec.pTransform)
    #         guessV[j] = lb(pvec.pTransform);
    #     end
    #     if v > ub(pvec.pTransform)
    #         guessV[j] = ub(pvec.pTransform);
    #     end
    # end
end

# """
# 	$(SIGNATURES)

# Perturb a guess, provided as a `ValueVector`. Return a new `ValueVector`.
# """
# function perturb_guess(m :: ModelObject, guess :: ValueVector, dGuess;
#     dIdx = nothing)

#     guessV = perturb_guess(m, get_values(guess), dGuess; dIdx = dIdx);
#     return ValueVector(guessV, lb(guess), ub(guess), pnames(guess))
# end

"""
	$(SIGNATURES)

Perturb calibrated model parameters. Including child objects.
"""
function perturb_params(m :: ModelObject, g :: Guess{F1}, dGuess; dIdx = nothing) where F1
    guessV = get_values(m, g);
    guess2V = perturb_guess_vector(m, g, guessV, dGuess; dIdx = dIdx);
    set_params_from_guess!(m, g, guess2V);
    return nothing
end

"""
    $(SIGNATURES)

Report all parameters by calibration status. For all ModelObjects contained in `o`.

Intended for reporting at the end (or during) a calibration run. Not formatted for inclusion in papers.

Each table row looks like:
"Description (name): value"
"""
function report_params(o :: ModelObject, isCalibrated :: Bool;
    io :: IO = stdout,
    closeToBounds :: Bool = false)

    # One model object + children at a time. So it can be skipped if no entries.
    mV = collect_model_objects(o; flatten = false);
    for m in mV
        lineV = report_params_one(m, isCalibrated;  closeToBounds);
        if length(lineV) > 1
            for line in lineV
                println(io, line);
            end
        end
    end
end

# This is for a Vector, such as [o2, o2.child1, o2.child2, o2.child2.gc1]
function report_params_one(mV :: Vector{T}, isCalibrated; closeToBounds) where T
    lineV = Vector{String}();
    for (j, m) in enumerate(mV)
        mLineV = report_params_one(m, isCalibrated; closeToBounds);
        # First object Id is always reported, so structure is visible.
        if !isempty(mLineV)  ||  (j == 1)
            push!(lineV, report_obj_id(get_object_id(m)));
        end
        if !isempty(mLineV)
            append!(lineV, mLineV);
        end
    end
    return lineV
end

function report_params_one(m, isCalibrated; closeToBounds)
    if has_pvector(m)
        pvec = get_pvector(m);
        lineV = report_pvec_params(pvec, isCalibrated; closeToBounds);
    else
        lineV = Vector{String}();
    end
    return lineV
end

# Return empty if no matching parameters found.
function report_pvec_params(pvec :: ParamVector, isCalibrated; closeToBounds)
    oId = get_object_id(pvec);
    nParents = n_parents(oId);
    lineV = Vector{String}();
    tbM = param_table(pvec, isCalibrated;  closeToBounds);
    if !isnothing(tbM)
        for ir = 1 : length(tbM)
            push!(lineV,
                indent_string(nParents + 1) *
                get_description(tbM, ir) *
                " (" * get_name(tbM, ir) * "):  " *
                get_value(tbM, ir));
        end
    end
    return lineV
end

indent_string(n) = "   " ^ n;

# Report ObjectId with indent depending on no of parents
function report_obj_id(oId :: ObjectId; nParents = n_parents(oId))
    if isempty(description(oId))
        descrStr = "";
    else
        descrStr = ":  " * description(oId);
    end
    return indent_string(nParents) * string(own_name(oId)) * descrStr;
end

"""
	$(SIGNATURES)

Generate a set of `ParamTable`s. One table per model object.
These are stored in a Dict with ObjectId's as keys.

The purpose is to make it easy to generate nicely formatted parameter tables that are grouped in a sensible way.
"""
function param_tables(o :: ModelObject, isCalibrated :: Bool)
    d = Dict{ObjectId, ParamTable}();

    # Iterating over ModelObjects rather than ParamVectors allows specific objects to override how their parameters are reported.
    objV = collect_model_objects(o);
    for obj in objV
        pt = param_table(obj, isCalibrated);
        if !isnothing(pt)
            d[get_object_id(obj)] = pt;
        end
    end
    return d
end

function param_table(o :: ModelObject, isCalibrated :: Bool)
    if has_pvector(o)
        pt = param_table(get_pvector(o), isCalibrated);
    else
        pt = nothing;
    end
    return pt
end

"""
	$(SIGNATURES)

Make a Latex parameter table for a set of model objects identified by their ObjectIds. Not all objects have ParamVectors.

Returns the table body only as a vector of String. Each element is a line in a Latex table. This can be embedded into a 3 column latex table with headers "Symbol & Description & Value".

# Arguments
- warnWhenMissing: Warn when an `ObjectId` does not have a parameter table (no own calibrated parameters). Default: `false`.
"""
function latex_param_table(o :: ModelObject, isCalibrated :: Bool,
    objIdV :: AbstractVector{ObjectId}, descrV :: AbstractVector{String};
    warnWhenMissing :: Bool = false)

    @assert size(objIdV) == size(descrV)

    # Param tables for all model objects. Not efficient, but convenient.
    d = param_tables(o, isCalibrated);

    lineV = Vector{String}();
    for (j, objId) ∈ enumerate(objIdV)
        if haskey(d, objId)
            newLineV = latex_param_table(d[objId], descrV[j]);
            if !isnothing(newLineV)
                append!(lineV, newLineV);
            end
        elseif warnWhenMissing
            @warn "$objId not found in $o"
        end
    end
    return lineV
end

# Latex parameter table for one `ParamVector`.
# `descr` is an optional description that is made into a multicolumn header.
function latex_param_table(pvec :: ParamVector, isCalibrated :: Bool,
    descr)

    tbM = param_table(pvec, isCalibrated);
    if isnothing(tbM)
        return nothing
    end
    return latex_param_table(tbM, descr)
end

"""
    $(SIGNATURES)

Number of calibrated parameters in object and its children. Also returns the number of elements (scalar values) in these parameters.
"""
function n_calibrated_params(o :: ModelObject; isCalibrated :: Bool = true)
    pvecV = collect_pvectors(o);
    nParam = 0;
    nElem = 0;
    if !isempty(pvecV)
        for (_, pvec) in pvecV
            nParam2, nElem2 = n_calibrated_params(pvec; isCalibrated);
            nParam += nParam2;
            nElem += nElem2;
        end
    end
    return nParam, nElem
end

end # module


text:
# NOTE: this model file is run from `/deps/generate_models.jl` or in the build phase.

# Model as in https://www.sas.upenn.edu/~jesusfv/ARE_Estimation.pdf

function FVGQ20()
    # Parameter values
    # δ = 0.025
    # ε = 10
    # ϕ = 0
    # γ2 = 0.001

    # Estimated parameters, taken from FV(2010), Table 3, p. 38, median estimate parameters
    # β = 0.998
    # h = 0.97
    # ϑ = 1.17
    # κ = 9.51
    # α = 0.21

    # θp = 0.82
    # χ = 0.63
    # γR = 0.77
    # γy = 0.19

    # γΠ = 1.29
    # Πbar = 1.01
    # ρd  = 0.12
    # ρφ = 0.93
    # σ_A = -3.97
    # σ_d = -1.51

    # σ_φ = -2.36
    # σ_μ = -5.43
    # σ_m = -5.85
    # Λμ = 3.4e-3
    # ΛA = 2.8e-3

    # 7 variable relabelings
    # 7 observations
    # 38 equations

    # 28 parameters. 23 of them to estimate, 5 parameters as constants
    @parameters δ ε ϕ γ2 β h ϑ κ α θp χ γR γy γΠ Πbar ρd ρφ ρg g_bar σ_A σ_d σ_φ σ_μ σ_m σ_g Λμ ΛA Ω_ii
    # States, 14
    @variables c_m x_m Π_m w_m R_m yd_m vp_m k d φ μ_I μ_A mshock g
    # Controls, 24
    @variables c x Π w R yd vp Πstar u q λ μ_z r l mc g1 g2 ob_μz ob_Π ob_R ob_dw ob_dy ob_l ob_μ
    # Markovize everything
    @make_markov c_m x_m Π_m w_m R_m yd_m vp_m k d φ μ_I μ_A mshock g c x Π w R yd vp Πstar u q λ μ_z r l mc g1 g2 ob_μz ob_Π ob_R ob_dw ob_dy ob_l ob_μ

    # Precompute some parameters, and steady state values
    ΛYd = (ΛA + α * Λμ) / (1 - α)
    Λx = exp(ΛYd)
    μ_z_val = exp(ΛYd)
    μ_I_val = exp(Λμ)
    γ1 = μ_z_val * μ_I_val / β - (1 - δ)
    Rbar = Πbar * μ_z_val / β
    Πstar_val = ((1 - θp * Πbar^(-(1 - ε) * (1 - χ))) / (1 - θp))^(1 / (1 - ε))
    vp_val = (1 - θp) / (1 - θp * Πbar^((1 - χ) * ε)) * Πstar_val^(-ε)
    mc_val = (1 - β * θp * Πbar^((1 - χ) * ε)) / (1 - β * θp * Πbar^(-(1 - ε) * (1 - χ))) * (ε - 1) / ε * Πstar_val
    w_val = (1 - α) * (mc_val * (α / γ1)^α)^(1 / (1 - α))
    k_val = α / (1 - α) * w_val / γ1 * μ_z_val * μ_I_val
    yd_val = (exp(ΛA) / μ_z_val * k_val^α - ϕ) / vp_val
    x_val = (μ_z_val * μ_I_val - (1 - δ)) / (μ_z_val * μ_I_val) * k_val
    c_val = (1 - g_bar) * yd_val - x_val
    λ_val = (1 - h * β / μ_z_val) / (1 - h / μ_z_val) / c_val
    ψ = w_val * λ_val
    g1_val = λ_val * yd_val / (1 - β * θp * Πbar^(-(1 - ε) * (1 - χ))) * (ε - 1) / ε * Πstar_val
    g2_val = λ_val * yd_val / (1 - β * θp * Πbar^(-(1 - ε) * (1 - χ))) * Πstar_val

    x_sym = [c_m, x_m, Π_m, w_m, R_m, yd_m, vp_m, k, d, φ, μ_I, μ_A, mshock, g]
    y_sym = [c, x, Π, w, R, yd, vp, Πstar, u, q, λ, μ_z, r, l, mc, g1, g2, ob_μz, ob_Π, ob_R, ob_dw, ob_dy, ob_l, ob_μ]
    p = [β, h, ϑ, κ, α, θp, χ, γR, γy, γΠ, Πbar, ρd, ρφ, ρg, g_bar, σ_A, σ_d, σ_φ, σ_μ, σ_m, σ_g, Λμ, ΛA]
    p_f = [δ, ε, ϕ, γ2, Ω_ii]

    H = [
         # Household Problem
         d * (c - h * c_m * μ_z^(-1))^(-1) - h * β * d_p * (c_p * μ_z_p - h * c)^(-1) - λ,
         ψ * φ * l^ϑ - w * λ, λ - β * λ_p * μ_z_p^(-1) * R / Π_p, γ1 + γ2 * (u - 1) - r,
         q -
         β * λ_p / λ / μ_z_p / μ_I_p *
         ((1 - δ) * q_p + r_p * u_p - (γ1 * (u_p - 1) + γ2 / 2 * (u_p - 1)^2)),
         q *
         (1 - (κ / 2 * (x / x_m * μ_z - Λx)^2) - (κ * (x / x_m * μ_z - Λx) * x / x_m * μ_z)) +
         β * q_p * λ_p / λ / μ_z_p * κ * (x_p / x * μ_z_p - Λx) * (x_p / x * μ_z_p)^2 - 1,

         # Firm Problem
         λ * mc * yd + β * θp * (Π^χ / Π_p)^(-ε) * g1_p - g1,
         λ * Πstar * yd + β * θp * (Π^χ / Π_p)^(1 - ε) * Πstar / Πstar_p * g2_p - g2,
         ε * g1 - (ε - 1) * g2, u * k / l - α / (1 - α) * w / r * μ_z * μ_I,
         (1 / (1 - α))^(1 - α) * (1 / α)^α * w^(1 - α) * r^α - mc,

         # Price Evolution
         θp * (Π_m^χ / Π)^(1 - ε) + (1 - θp) * Πstar^(1 - ε) - 1,

         # Taylor Rule
         R / Rbar -
         (R_m / Rbar)^γR *
         ((Π / Πbar)^γΠ * ((yd / yd_m * μ_z) / exp(ΛYd))^γy)^(1 - γR) *
         mshock,

         # Market Clear
         c + x + g + μ_z^(-1) * μ_I^(-1) * (γ1 * (u - 1) + γ2 / 2 * (u - 1)^2) * k - yd,
         (μ_A / μ_z * (u * k)^α * l^(1 - α) - ϕ) / vp - yd,
         θp * (Π_m^χ / Π)^(-ε) * vp_m + (1 - θp) * Πstar^(-ε) - vp,
         k_p * μ_z * μ_I - (1 - δ) * k -
         μ_z * μ_I * (1 - κ / 2 * (x / x_m * μ_z - Λx)^2) * x,

         # Variable Mapping
         c - c_m_p, x - x_m_p, Π - Π_m_p, w - w_m_p, R - R_m_p, yd - yd_m_p, vp - vp_m_p,

         # Shock Evolution
         μ_z - μ_A^(1 / (1 - α)) * μ_I^(α / (1 - α)), log(d_p) - ρd * log(d),
         log(φ_p) - ρφ * log(φ), log(μ_I_p) - Λμ, log(μ_A_p) - ΛA, log(mshock_p) - 0,
         log(g_p) - ρg * log(g) - (1 - ρg) * log(g_bar * yd_val),

         # Observation equations Π, R, dw, dy, l, μ_I(-1)
         ob_μz - log(μ_z), ob_Π - log(Π), ob_R - log(R), ob_dw - (log(w) - log(w_m)),
         ob_dy - (log(yd) - log(yd_m)), ob_l - log(l), ob_μ - log(μ_I)]

    x̄ = [c_m_ss ~ c_val, x_m_ss ~ x_val, Π_m_ss ~ Πbar, w_m_ss ~ w_val, R_m_ss ~ Rbar,
          yd_m_ss ~ yd_val, vp_m_ss ~ vp_val, k_ss ~ k_val, d_ss ~ 1.0, φ_ss ~ 1.0,
          μ_I_ss ~ exp(Λμ), μ_A_ss ~ exp(ΛA), mshock_ss ~ 1.0, g_ss ~ g_bar * yd_val]
    ȳ = [c_ss ~ c_val, x_ss ~ x_val, Π_ss ~ Πbar, w_ss ~ w_val, R_ss ~ Rbar,
          yd_ss ~ yd_val, vp_ss ~ vp_val, Πstar_ss ~ Πstar_val, u_ss ~ 1.0, q_ss ~ 1.0,
          λ_ss ~ λ_val, μ_z_ss ~ exp(ΛYd), r_ss ~ γ1, l_ss ~ 1.0, mc_ss ~ mc_val,
          g1_ss ~ g1_val, g2_ss ~ g2_val, ob_μz_ss ~ ΛYd, ob_Π_ss ~ log(Πbar),
          ob_R_ss ~ log(Rbar), ob_dw_ss ~ 0.0, ob_dy_ss ~ 0.0, ob_l_ss ~ 0.0, ob_μ_ss ~ Λμ]

    # x̄_iv = deepcopy(x̄)
    # ȳ_iv = deepcopy(ȳ)

    n_ϵ = 6
    n_x = length(x_sym)
    n_y = length(y_sym)
    n_z = 6
    Γ = zeros(Operation, n_ϵ, n_ϵ) # Make sure it is not a float64 matrix
    Γ[1, 1] = σ_d
    Γ[2, 2] = σ_φ
    Γ[3, 3] = σ_μ
    Γ[4, 4] = σ_A
    Γ[5, 5] = σ_m
    Γ[6, 6] = σ_g
    η_mat = zeros(n_x, n_ϵ)
    # The 9-14 states are: d φ μ_I μ_A mshock g
    η_mat[9, 1] = 1
    η_mat[10, 2] = 1
    η_mat[11, 3] = 1
    η_mat[12, 4] = 1
    η_mat[13, 5] = 1
    η_mat[14, 6] = 1
    # Define Q
    Q = zeros(n_z, n_x + n_y)
    Q[:, (n_y - 5):n_y] = [1.0 0.0 0.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.0 0.0 0.0; 0.0 0.0 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0]
    # Define Ω
    Ω = fill(Ω_ii, n_z)

    return H,
           (x = x_sym, y = y_sym, x̄ = x̄, ȳ = ȳ, Γ = Γ, η = η_mat, p = p, p_f = p_f,
            x̄_iv = nothing, ȳ_iv = nothing, Q = Q, Ω = Ω), "FVGQ20"

end


text:
include("Blocks.jl")
include("Electric.jl")
include("Rotational.jl")


text:
using Documenter, Finch

makedocs(
    sitename="Finch",
    pages = [
        "index.md",
        "start.md",
        "Usage" => [
            "configuration.md",
            "mesh.md",
            "entities.md",
            "conditions.md",
            "equation.md",
            "solution.md",
            "reorder.md",
            "misc.md"
        ],
        "examples.md"
    ]
)

deploydocs(
    repo = "github.com/paralab/Finch.git",
)

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
# kspace.jl
# kspace sampling patterns
# 2018-06-13 Jeff Fessler, University of Michigan

using Plots: plot, plot!, scatter, gui
using Test

"""
`kspace = ir_mri_kspace_ga_radial(; Nro=?, Nspoke=?, ...)`

Generate k-space sampling pattern for "golden angle" radial sampling.

option
* `Nro:Int`		number of samples in each readout/spoke, default 256
* `Nspoke::Int`	number of spokes, default 1
* `start::Real`	first angle in series [radians], default pi/2
* `angle::Real`	angular spacing [radians], default GA
* `delta_ro::Real`	readout spacing, default 1/Nro
* `shift::Real`		shift due to gradient delays, default 0
		radial sample locations are `ir * delta_ro`
		where `ir = [-(Nro/2 - 1):1:Nro/2] + shift`
* `show::Bool`		plot k-space locations, default `false`

out
* `kspace`	`[Nro Nspoke 2]` (Float32)
			kx and ky k-space locations for Nspoke*Nro samples
			in interval (-0.5 0.5] for default shift, delta_ro
			so default units are "cycles / sample"

2015-07 Mai Le, original Matlab version
2015-07-04 Jeff Fessler, minor changes to Matlab version
"""
function ir_mri_kspace_ga_radial(;
		Nro::Int = 256,
		Nspoke::Int = 1,
		delta_ro::Real = 1/Nro,
		shift::Real = 0,
		angle::Real = pi*(sqrt(5)-1)/2, # golden angle [radians]
		start::Real = pi/2,
		show::Bool = false,
	)

	isodd(Nro) && throw(ArgumentError("Number along read out not even!"))

	rho = (-(Nro/2 - 1):1:Nro/2) * delta_ro # radial samples
	rho = rho .+ shift * delta_ro # apply gradient delay shift

	phi = (0:Nspoke-1) * angle .+ start # angles

	rho = Float32.(rho)
	phi = Float32.(phi)

	kx = rho * cos.(phi)' # [Nro Nspoke]
	ky = rho * sin.(-phi)'
	kspace = cat(dims=3, kx, ky) # [Nro Nspoke 2]

	if show
		scatter(kx, ky, label="", aspect_ratio=1)
		kmax = delta_ro * Nro / 2
		plot!(xlim = [-1,1]*kmax)
		plot!(ylim = [-1,1]*kmax)
		plot!(xtick = (-1:1) * kmax, xlabel = "k_x")
		plot!(ytick = (-1:1) * kmax, ylabel = "k_y")
		gui()
	end

	return kspace
end

"""
`ir_mri_kspace_ga_radial`
self test
"""
function ir_mri_kspace_ga_radial(test::Symbol)
	test != :test && throw(":test")
	k = ir_mri_kspace_ga_radial(Nro=30, Nspoke=13, start=pi,
		show=true, shift=-0.5, delta_ro=1)
	@test size(k) == (30, 13, 2)
	true
end

# ir_mri_kspace_ga_radial(:test)


text:
# Events for timing

export CuEvent, record, synchronize, elapsed, @elapsed

const CuEvent_t = Ptr{Cvoid}

@enum(CUevent_flags, EVENT_DEFAULT = Cint(0),
                     EVENT_BLOCKING_SYNC = Cint(1),
                     EVENT_DISABLE_TIMING = Cint(2),
                     EVENT_INTERPROCESS = Cint(4))

# FIXME: EnumSet from JuliaLang/julia#19470
Base.:|(x::CUevent_flags, y::CUevent_flags) =
    reinterpret(CUevent_flags, Base.cconvert(Unsigned, x) | Base.cconvert(Unsigned, y))

"""
    CuEvent()

Create a new CUDA event.
"""
mutable struct CuEvent
    handle::CuEvent_t
    ctx::CuContext

    function CuEvent(flags=EVENT_DEFAULT)
        handle_ref = Ref{CuEvent_t}()
        @apicall(:cuEventCreate, (Ptr{CuEvent_t}, Cuint), handle_ref, flags)

        ctx = CuCurrentContext()
        obj = new(handle_ref[], ctx)
        finalizer(unsafe_destroy!, obj)
        return obj
    end
end

function unsafe_destroy!(e::CuEvent)
    if isvalid(e.ctx)
        @apicall(:cuEventDestroy, (CuEvent_t,), e)
    end
end

Base.unsafe_convert(::Type{CuEvent_t}, e::CuEvent) = e.handle

Base.:(==)(a::CuEvent, b::CuEvent) = a.handle == b.handle
Base.hash(e::CuEvent, h::UInt) = hash(e.handle, h)

"""
    record(e::CuEvent, stream=CuDefaultStream())

Record an event on a stream.
"""
record(e::CuEvent, stream::CuStream=CuDefaultStream()) =
    @apicall(:cuEventRecord, (CuEvent_t, CuStream_t), e, stream)

"""
    synchronize(e::CuEvent)

Waits for an event to complete.
"""
synchronize(e::CuEvent) = @apicall(:cuEventSynchronize, (CuEvent_t,), e)

"""
    wait(e::CuEvent, stream=CuDefaultStream())

Make a stream wait on a event. This only makes the stream wait, and not the host; use
[`synchronize(::CuEvent)`](@ref) for that.
"""
wait(e::CuEvent, stream::CuStream=CuDefaultStream()) =
    @apicall(:cuStreamWaitEvent, (CuStream_t, CuEvent_t, Cuint), stream, e, 0)

"""
    elapsed(start::CuEvent, stop::CuEvent)

Computes the elapsed time between two events (in seconds).
"""
function elapsed(start::CuEvent, stop::CuEvent)
    time_ref = Ref{Cfloat}()
    @apicall(:cuEventElapsedTime, (Ptr{Cfloat}, CuEvent_t, CuEvent_t),
                                  time_ref, start, stop)
    return time_ref[]/1000
end

"""
    @elapsed stream ex
    @elapsed ex

A macro to evaluate an expression, discarding the resulting value, instead returning the
number of seconds it took to execute on the GPU, as a floating-point number.
"""
macro elapsed(stream, ex)
    quote
        t0, t1 = CuEvent(), CuEvent()
        record(t0, $(esc(stream)))
        $(esc(ex))
        record(t1, $(esc(stream)))
        synchronize(t1)
        elapsed(t0, t1)
    end
end
macro elapsed(ex)
    quote
        @elapsed(CuDefaultStream(), $(esc(ex)))
    end
end


text:
# ---
# title: "`erode`"
# id: op_erode
# cover: assets/erode.png
# ---

# The dilation operator `erode` is essentially a min filter. This is a basic term in
# mathematical morphology -- many operations are built on top of `erode` and the
# [`dilate`](@ref) operator.

using ImageMorphology
using TestImages
using ImageBase
using ImageShow

img = restrict(testimage_dip3e("fig1005")) # wirebond mask

# For each pixel, dilation is the minimum of the pixels in the neighborhood. In mathematics,
# `erode(A, Ω)` is defined as ``\delta_A[p] = \inf\{A[p+o] | o \in \Omega\}`` where `Ω` is
# the [structuring element](@ref concept_se).

# The simplest usage is `erode(img; [dims], [r])`, where `dims` and `r` controls the
# neighborhood shape.

out1 = erode(img) # default: all spatial dimensions, r=1, a box-shape SE
out2 = erode(img; dims=(2,)) # only apply to the second dimension
out3 = erode(img; r=5) # half-size r=5
mosaic(out1, out2, out3; nrow=1)

# It uses the [`strel_box`](@ref) to create a box-shaped structuring element. You can also
# provide a custom SE via the `erode(img, se)` interface.

out1 = erode(img, strel_box((3, 3))) # default se for`erode(img)`
se = centered(Bool[1 1 1; 1 1 0; 0 0 0]) # match top-left region
out2 = erode(img, se)
mosaic(out1, out2; nrow=1)

# An in-place version `erode!` is also provided, for instance

out1 = similar(img)
erode!(out1, img)

out2 = similar(img)
erode!(out2, img, strel_diamond((3, 3)))
#md nothing #hide

# ## See also

# [`erode`](@ref op_erode) is the dual operator of `dilate` in the following sense:

complement.(erode(img)) == erode(complement.(img))

# For bool arrays and [symmetric SEs](@ref concept_symmetric), erosion becomes equivalent
# to the minkowski difference on sets: ``A \ominus B = \{ a-b | a \in A, b \in B \}``.

# For a comprehensive and more accurate documentation, please check the [`erode`](@ref)
# reference page.

# save cover image #src
using FileIO #src
mkpath("assets") #src
img = Gray.(restrict(testimage("blobs"))) #src
save("assets/erode.png", erode(img; r=1)) #src


text:
"""
    DataFrame(m::WSVarLmmModel)

Convert the data in `WSVarLmmModel` to a `DataFrame`.
"""
function DataFrame(m::WSVarLmmModel)
    p, q, l, n = m.p, m.q, m.l, m.nsum
    # preallocate arrays
    id = Vector{eltype(m.ids)}(undef, n)
    y  = Matrix{Float64}(undef, n, 1)
    X  = Matrix{Float64}(undef, n, p)
    Z  = Matrix{Float64}(undef, n, q)
    W  = Matrix{Float64}(undef, n, l)
    addweight = !isempty(m.obswts)
    if addweight
        weights = Vector{Float64}(undef, n)
    end

    #get names to reconstruct Dataframe properly
    meannames = map(x -> join(x[2:end], ": "),  (split.(m.meannames, ": ")))
    renames = map(x -> join(x[2:end], ": "),  (split.(m.renames, ": ")))
    wsvarnames = map(x -> join(x[2:end], ": "),  (split.(m.wsvarnames, ": ")))

    # gather data
    offset = 1
    for (i, vlmmobs) in enumerate(m.data)
        ni            = length(vlmmobs.y)
        rangei        = offset:(offset + ni - 1)
        id[rangei]   .= m.ids[i]
        y[rangei]    = vlmmobs.y
        X[rangei, :] = transpose(vlmmobs.Xt)
        Z[rangei, :] = transpose(vlmmobs.Zt)
        W[rangei, :] = transpose(vlmmobs.Wt)
        if addweight
            weights[rangei] .= m.obswts[i]
        end
        offset       += ni
    end

    # Skip repeated variables in X, Z, W and don't copy interaction terms
    Xcols = .![occursin("&", meanname) for meanname in meannames]
    Zcols = .![(zname in meannames || occursin("&", zname)) for zname in renames]
    Wcols = .![(wname in meannames || wname in renames ||
        occursin("&", wname)) for wname in wsvarnames]

    # Construct DataFrame
    df = hcat(DataFrame(id = id),
        DataFrame(y, [m.respname]),
        DataFrame(X[:, Xcols], meannames[Xcols]),
        DataFrame(Z[:, Zcols], renames[Zcols]),
        DataFrame(W[:, Wcols], wsvarnames[Wcols]))
    if addweight
        df[!, :obswts] = weights
    end

    allnames = union(meannames, renames, wsvarnames)
    filter!(x -> !occursin("&", x), allnames)

    # Create columns for factored variables (for reusing original formula and GWAS)
    isfactor = occursin.(": ", allnames)
    splitnames = map(x -> split(x, ": ")[1], allnames)

    for name in unique(splitnames)
        numfactors = findall(isequal.(name, splitnames))
        if (length(numfactors) == 1 && .!(any(isfactor[numfactors])))
            continue
        end
        factors = map(x -> join(split(x, ": ")[2:end], ": "),
            allnames[numfactors])
        valholder = repeat([""], n)
        for factor in factors
            dfname = join([name, factor], ": ")
            valholder[df[!, dfname] .== 1] .= factor
        end
        valholder[valholder .== ""] .= "Reference"
        df[!, name] = valholder
        df[!, name] = levels!(CategoricalArray(df[!, name]),
            ["Reference"; factors]);
    end
    df
end

"""
    WSVarLmmModel(meanformula::FormulaTerm, reformula::FormulaTerm,
    wsvarformula::FormulaTerm, idvar::Union{String, Symbol}, tbl)

Constructor of `WSVarLmmModel` from a Tables.jl compatible source.

# Positional arguments
- `meanformula`: formula for the mean fixed effects β (variables in X matrix).
- `reformula`: formula for the mean random effects γ (variables in Z matrix).
- `wsvarformula`: formula for the within-subject variance effects τ (variables in W matrix).
- `idvar`: id variable for groupings.
- `tbl:` data table holding all of the data for the model. It can be a
`DataFrame` or column-based table such as an `IndexedTable` from JuliaDB.

# Keyword arguments
- `wtvar`: variable name corresponding to the observation weights in the datatable.

# Example
```
vlmm3 = WSVarLmmModel(@formula(y ~ 1 + x2 + x3 + x4 + x5),
    @formula(y ~ 1 + z2 + z3), @formula(y ~ 1 + w2 + w3 + w4 + w5), "id", df)
```
"""
WSVarLmmModel(
    meanformula  :: FormulaTerm,
    reformula    :: FormulaTerm,
    wsvarformula :: FormulaTerm,
    idvar        :: Union{String, Symbol},
    tbl;
    wtvar        :: Union{String, Symbol} = ""
    ) = WSVarLmmModel(meanformula, reformula, wsvarformula,
     idvar, columntable(tbl); wtvar = wtvar)

function WSVarLmmModel(
    meanformula  :: FormulaTerm,
    reformula    :: FormulaTerm,
    wsvarformula :: FormulaTerm,
    idvar        :: Union{String, Symbol},
    tbl          :: T;
    wtvar        :: Union{String, Symbol} = ""
    ) where T <: Tables.ColumnTable

    idvar = Symbol(idvar)
    iswtvar = !isempty(string(wtvar))

    # keep original non-schema applied formulas for model object
    orig_meanformla = deepcopy(meanformula)
    orig_reformla = deepcopy(reformula)
    orig_wsvarformla = deepcopy(wsvarformula)

    function varlmmobs(tab)
        y, X = modelcols(meanformula, tab)
        Z    = modelmatrix(reformula, tab)
        W    = modelmatrix(wsvarformula, tab)
        return WSVarLmmObs(y, X, Z, W)
    end

    # collect all terms to perform dropping properly
    if iswtvar
        alltermformula = meanformula.lhs ~ sum(term.(union(terms(meanformula.rhs),
            terms(reformula.rhs), terms(wsvarformula.rhs)))) +
            term(idvar) + term(Symbol(wtvar))
    else
        alltermformula = meanformula.lhs ~ sum(term.(union(terms(meanformula.rhs),
          terms(reformula.rhs), terms(wsvarformula.rhs)))) + term(idvar)
    end

    tbl, _ = StatsModels.missing_omit(tbl, alltermformula)

    # apply df-wide schema
    meanformula  = apply_schema(meanformula, schema(meanformula, tbl))#, ydict))
    reformula    = apply_schema(reformula, schema(reformula, tbl))#, ydict))
    wsvarformula = apply_schema(wsvarformula, schema(wsvarformula, tbl))#, ydict))

    # variable names
    respname = string(meanformula.lhs)
    meannames = StatsModels.coefnames(meanformula.rhs)
    # either array{Names} or string of one variable
    meannames = typeof(meannames) <: Array ? ["β$i: " for i in 1:length(meannames)] .*
        meannames : ["β1: " * meannames]
    renames = StatsModels.coefnames(reformula.rhs)
    renames = typeof(renames) <: Array ? ["γ$i: " for i in 1:length(renames)] .*
        renames : ["γ1: " * renames]
    wsvarnames = StatsModels.coefnames(wsvarformula.rhs)
    wsvarnames = typeof(wsvarnames) <: Array ?  ["τ$i: " for i in 1:length(wsvarnames)] .*
        wsvarnames : ["τ1: " * wsvarnames]

    # build grouped dataframe and record ID order
    gdf = DataFrames.groupby(DataFrame(tbl, copycols = false), idvar)
    ids = map(x -> x[1], keys(gdf))

    if isempty(string(wtvar))
        wts = []
    else
        cnames = Tables.columnnames(tbl)
        wtvar  = Symbol(wtvar)
        wtvar in cnames ||
            error("weight variable $wtvar not in datatable $tbl")
        wts = combine(gdf, wtvar => first)[!, 2]
    end
    obsvec = combine(varlmmobs, gdf)[!, 2]
    varlmm = WSVarLmmModel(obsvec, respname = respname,
        meannames = meannames, renames = renames,
        wsvarnames = wsvarnames, ids = ids, obswts = wts,
        meanformula = orig_meanformla, reformula = orig_reformla,
        wsvarformula = orig_wsvarformla)
    return varlmm
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module TestHelpers

include("dimensionful.jl")
export Furlong

mutable struct FakeTerminal <: Base.Terminals.UnixTerminal
    in_stream::Base.IO
    out_stream::Base.IO
    err_stream::Base.IO
    hascolor::Bool
    raw::Bool
    FakeTerminal(stdin,stdout,stderr,hascolor=true) =
        new(stdin,stdout,stderr,hascolor,false)
end

Base.Terminals.hascolor(t::FakeTerminal) = t.hascolor
Base.Terminals.raw!(t::FakeTerminal, raw::Bool) = t.raw = raw
Base.Terminals.size(t::FakeTerminal) = (24, 80)

function open_fake_pty()
    @static if is_windows()
        error("Unable to create a fake PTY in Windows")
    end

    const O_RDWR = Base.Filesystem.JL_O_RDWR
    const O_NOCTTY = Base.Filesystem.JL_O_NOCTTY

    fdm = ccall(:posix_openpt, Cint, (Cint,), O_RDWR|O_NOCTTY)
    fdm == -1 && error("Failed to open PTY master")
    rc = ccall(:grantpt, Cint, (Cint,), fdm)
    rc != 0 && error("grantpt failed")
    rc = ccall(:unlockpt, Cint, (Cint,), fdm)
    rc != 0 && error("unlockpt")

    fds = ccall(:open, Cint, (Ptr{UInt8}, Cint),
        ccall(:ptsname, Ptr{UInt8}, (Cint,), fdm), O_RDWR|O_NOCTTY)

    # slave
    slave = RawFD(fds)
    master = Base.TTY(RawFD(fdm); readable = true)
    slave, master
end

function with_fake_pty(f)
    slave, master = open_fake_pty()
    try
        f(slave, master)
    finally
        ccall(:close,Cint,(Cint,),slave) # XXX: this causes the kernel to throw away all unread data on the pty
        close(master)
    end
end

function challenge_prompt(code::AbstractString, challenges; timeout::Integer=10, debug::Bool=true)
    output_file = tempname()
    wrapped_code = """
    result = let
        $code
    end
    open("$output_file", "w") do fp
        serialize(fp, result)
    end
    """
    cmd = `$(Base.julia_cmd()) --startup-file=no -e $wrapped_code`
    try
        challenge_prompt(cmd, challenges, timeout=timeout, debug=debug)
        return open(output_file, "r") do fp
            deserialize(fp)
        end
    finally
        isfile(output_file) && rm(output_file)
    end
    return nothing
end

function challenge_prompt(cmd::Cmd, challenges; timeout::Integer=10, debug::Bool=true)
    function format_output(output)
        debug ? "Process output found:\n\"\"\"\n$(readstring(seekstart(out)))\n\"\"\"" : ""
    end
    out = IOBuffer()
    with_fake_pty() do slave, master
        p = spawn(detach(cmd), slave, slave, slave)
        # Kill the process if it takes too long. Typically occurs when process is waiting for input
        @async begin
            sleep(timeout)
            kill(p)
            close(master)
        end
        try
            for (challenge, response) in challenges
                process_exited(p) && error("Too few prompts. $(format_output(out))")

                write(out, readuntil(master, challenge))
                if !isopen(master)
                    error("Could not locate challenge: \"$challenge\". $(format_output(out))")
                end
                write(master, response)
            end
            wait(p)
        finally
            kill(p)
        end
        # Determine if the process was explicitly killed
        killed = process_exited(p) && (p.exitcode != 0 || p.termsignal != 0)
        killed && error("Too many prompts. $(format_output(out))")
    end
    nothing
end

# OffsetArrays (arrays with indexing that doesn't start at 1)

# This test file is designed to exercise support for generic indexing,
# even though offset arrays aren't implemented in Base.

module OAs

using Base: Indices, IndexCartesian, IndexLinear, tail

export OffsetArray

struct OffsetArray{T,N,AA<:AbstractArray} <: AbstractArray{T,N}
    parent::AA
    offsets::NTuple{N,Int}
end
OffsetVector{T,AA<:AbstractArray} = OffsetArray{T,1,AA}

OffsetArray{T,N}(A::AbstractArray{T,N}, offsets::NTuple{N,Int}) = OffsetArray{T,N,typeof(A)}(A, offsets)
OffsetArray{T,N}(A::AbstractArray{T,N}, offsets::Vararg{Int,N}) = OffsetArray(A, offsets)

(::Type{OffsetArray{T,N}}){T,N}(inds::Indices{N}) = OffsetArray{T,N,Array{T,N}}(Array{T,N}(map(length, inds)), map(indsoffset, inds))
(::Type{OffsetArray{T}}){T,N}(inds::Indices{N}) = OffsetArray{T,N}(inds)

Base.IndexStyle{T<:OffsetArray}(::Type{T}) = Base.IndexStyle(parenttype(T))
parenttype{T,N,AA}(::Type{OffsetArray{T,N,AA}}) = AA
parenttype(A::OffsetArray) = parenttype(typeof(A))

Base.parent(A::OffsetArray) = A.parent

errmsg(A) = error("size not supported for arrays with indices $(indices(A)); see https://docs.julialang.org/en/latest/devdocs/offset-arrays/")
Base.size(A::OffsetArray) = errmsg(A)
Base.size(A::OffsetArray, d) = errmsg(A)
Base.eachindex(::IndexCartesian, A::OffsetArray) = CartesianRange(indices(A))
Base.eachindex(::IndexLinear, A::OffsetVector) = indices(A, 1)

# Implementations of indices and indices1. Since bounds-checking is
# performance-critical and relies on indices, these are usually worth
# optimizing thoroughly.
@inline Base.indices(A::OffsetArray, d) = 1 <= d <= length(A.offsets) ? indices(parent(A))[d] + A.offsets[d] : (1:1)
@inline Base.indices(A::OffsetArray) = _indices(indices(parent(A)), A.offsets)  # would rather use ntuple, but see #15276
@inline _indices(inds, offsets) = (inds[1]+offsets[1], _indices(tail(inds), tail(offsets))...)
_indices(::Tuple{}, ::Tuple{}) = ()
Base.indices1{T}(A::OffsetArray{T,0}) = 1:1  # we only need to specialize this one

function Base.similar(A::OffsetArray, T::Type, dims::Dims)
    B = similar(parent(A), T, dims)
end
function Base.similar(A::AbstractArray, T::Type, inds::Tuple{UnitRange,Vararg{UnitRange}})
    B = similar(A, T, map(length, inds))
    OffsetArray(B, map(indsoffset, inds))
end

Base.similar(f::Union{Function,Type}, shape::Tuple{UnitRange,Vararg{UnitRange}}) = OffsetArray(f(map(length, shape)), map(indsoffset, shape))

Base.reshape(A::AbstractArray, inds::Tuple{UnitRange,Vararg{UnitRange}}) = OffsetArray(reshape(A, map(length, inds)), map(indsoffset, inds))

@inline function Base.getindex{T,N}(A::OffsetArray{T,N}, I::Vararg{Int,N})
    checkbounds(A, I...)
    @inbounds ret = parent(A)[offset(A.offsets, I)...]
    ret
end
# Vectors don't support one-based linear indexing; they always use the offsets
@inline function Base.getindex(A::OffsetVector, i::Int)
    checkbounds(A, i)
    @inbounds ret = parent(A)[offset(A.offsets, (i,))[1]]
    ret
end
# But multidimensional arrays allow one-based linear indexing
@inline function Base.getindex(A::OffsetArray, i::Int)
    checkbounds(A, i)
    @inbounds ret = parent(A)[i]
    ret
end
@inline function Base.setindex!{T,N}(A::OffsetArray{T,N}, val, I::Vararg{Int,N})
    checkbounds(A, I...)
    @inbounds parent(A)[offset(A.offsets, I)...] = val
    val
end
@inline function Base.setindex!(A::OffsetVector, val, i::Int)
    checkbounds(A, i)
    @inbounds parent(A)[offset(A.offsets, (i,))[1]] = val
    val
end
@inline function Base.setindex!(A::OffsetArray, val, i::Int)
    checkbounds(A, i)
    @inbounds parent(A)[i] = val
    val
end

@inline function Base.deleteat!(A::OffsetArray, i::Int)
    checkbounds(A, i)
    @inbounds deleteat!(parent(A), offset(A.offsets, (i,))[1])
end

@inline function Base.deleteat!{T,N}(A::OffsetArray{T,N}, I::Vararg{Int, N})
    checkbounds(A, I...)
    @inbounds deleteat!(parent(A), offset(A.offsets, I)...)
end

@inline function Base.deleteat!(A::OffsetArray, i::UnitRange{Int})
    checkbounds(A, first(i))
    checkbounds(A, last(i))
    first_idx = offset(A.offsets, (first(i),))[1]
    last_idx = offset(A.offsets, (last(i),))[1]
    @inbounds deleteat!(parent(A), first_idx:last_idx)
end

# Computing a shifted index (subtracting the offset)
offset{N}(offsets::NTuple{N,Int}, inds::NTuple{N,Int}) = _offset((), offsets, inds)
_offset(out, ::Tuple{}, ::Tuple{}) = out
@inline _offset(out, offsets, inds) = _offset((out..., inds[1]-offsets[1]), Base.tail(offsets), Base.tail(inds))

indsoffset(r::Range) = first(r) - 1
indsoffset(i::Integer) = 0

Base.resize!(A::OffsetVector, nl::Integer) = (resize!(A.parent, nl); A)

end

end

max_stars_repo_name: omus/julia

max_issues_repo_name: omus/julia

max_forks_repo_name: omus/julia


text:
module IntegralsCuba

using Integrals, Cuba
import Integrals: transformation_if_inf, scale_x, scale_x!

abstract type AbstractCubaAlgorithm <: SciMLBase.AbstractIntegralAlgorithm end
struct CubaVegas <: AbstractCubaAlgorithm end
struct CubaSUAVE <: AbstractCubaAlgorithm end
struct CubaDivonne <: AbstractCubaAlgorithm end
struct CubaCuhre <: AbstractCubaAlgorithm end

function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgorithm, sensealg,
    lb, ub, p, args...;
    reltol=1e-8, abstol=1e-8,
    maxiters=alg isa CubaSUAVE ? 1000000 : typemax(Int),
    kwargs...)
    prob = transformation_if_inf(prob) #intercept for infinite transformation
    p = p
    if lb isa Number && prob.batch == 0
        _x = Float64[lb]
    elseif lb isa Number
        _x = zeros(length(lb), prob.batch)
    elseif prob.batch == 0
        _x = zeros(length(lb))
    else
        _x = zeros(length(lb), prob.batch)
    end
    ub = ub
    lb = lb

    if prob.batch == 0
        if isinplace(prob)
            f = function (x, dx)
                prob.f(dx, scale_x!(_x, ub, lb, x), p)
                dx .*= prod((y) -> y[1] - y[2], zip(ub, lb))
            end
        else
            f = function (x, dx)
                dx .= prob.f(scale_x!(_x, ub, lb, x), p) .* prod((y) -> y[1] - y[2], zip(ub, lb))
            end
        end
    else
        if lb isa Number
            if isinplace(prob)
                f = function (x, dx)
                    #todo check scale_x!
                    prob.f(dx', scale_x!(view(_x, 1:length(x)), ub, lb, x), p)
                    dx .*= prod((y) -> y[1] - y[2], zip(ub, lb))
                end
            else
                if prob.f([lb ub], p) isa Vector
                    f = function (x, dx)
                        dx .= prob.f(scale_x(ub, lb, x), p)' .* prod((y) -> y[1] - y[2], zip(ub, lb))
                    end
                else
                    f = function (x, dx)
                        dx .= prob.f(scale_x(ub, lb, x), p) .* prod((y) -> y[1] - y[2], zip(ub, lb))
                    end
                end
            end
        else
            if isinplace(prob)
                f = function (x, dx)
                    prob.f(dx, scale_x(ub, lb, x), p)
                    dx .*= prod((y) -> y[1] - y[2], zip(ub, lb))
                end
            else
                if prob.f([lb ub], p) isa Vector
                    f = function (x, dx)
                        dx .= prob.f(scale_x(ub, lb, x), p)' .* prod((y) -> y[1] - y[2], zip(ub, lb))
                    end
                else
                    f = function (x, dx)
                        dx .= prob.f(scale_x(ub, lb, x), p) .* prod((y) -> y[1] - y[2], zip(ub, lb))
                    end
                end
            end
        end
    end

    ndim = length(lb)

    nvec = prob.batch == 0 ? 1 : prob.batch

    if alg isa CubaVegas
        out = Cuba.vegas(f, ndim, prob.nout; rtol=reltol,
            atol=abstol, nvec=nvec,
            maxevals=maxiters, kwargs...)
    elseif alg isa CubaSUAVE
        out = Cuba.suave(f, ndim, prob.nout; rtol=reltol,
            atol=abstol, nvec=nvec,
            maxevals=maxiters, kwargs...)
    elseif alg isa CubaDivonne
        out = Cuba.divonne(f, ndim, prob.nout; rtol=reltol,
            atol=abstol, nvec=nvec,
            maxevals=maxiters, kwargs...)
    elseif alg isa CubaCuhre
        out = Cuba.cuhre(f, ndim, prob.nout; rtol=reltol,
            atol=abstol, nvec=nvec,
            maxevals=maxiters, kwargs...)
    end

    if isinplace(prob) || prob.batch != 0
        val = out.integral
    else
        if prob.nout == 1 && prob.f(lb, p) isa Number
            val = out.integral[1]
        else
            val = out.integral
        end
    end

    SciMLBase.build_solution(prob, alg, val, out.error,
        chi=out.probability, retcode=:Success)
end

export CubaVegas, CubaSUAVE, CubaDivonne, CubaCuhre

end


text:
@testset "asymptotic_block.jl" begin
    @testset "estimate, stderr, and z" begin
        kernel = (ExponentialKernel() ∘ ScaleTransform(0.1)) ⊗ WhiteKernel()
        for nclasses in (2, 10, 100), nsamples in (10, 50, 100)
            for blocksize in (2, 5, 10, 50)
                # blocksize may no be greater than number of samples
                blocksize < nsamples || continue

                # sample predictions and targets
                dist = Dirichlet(nclasses, 1)
                predictions = [rand(dist) for _ in 1:nsamples]
                targets_consistent = [
                    rand(Categorical(prediction)) for prediction in predictions
                ]
                targets_onlyone = ones(Int, length(predictions))

                skce = BlockUnbiasedSKCE(kernel, blocksize)

                # for both sets of targets
                for targets in (targets_consistent, targets_onlyone)
                    test = AsymptoticBlockSKCETest(kernel, blocksize, predictions, targets)

                    @test test.blocksize == blocksize
                    @test test.nblocks == nsamples ÷ blocksize
                    @test test.estimate ≈ skce(predictions, targets)
                    @test test.z == test.estimate / test.stderr

                    @test pvalue(test) ==
                          pvalue(test; tail=:right) ==
                          ccdf(Normal(), test.z)
                    @test_throws ArgumentError pvalue(test; tail=:left)
                    @test_throws ArgumentError pvalue(test; tail=:both)

                    for α in 0.55:0.05:0.95
                        q = quantile(Normal(), α)
                        @test confint(test; level=α) ==
                              confint(test; level=α, tail=:right) ==
                              (max(0, test.estimate - q * test.stderr), Inf)
                        @test_throws ArgumentError confint(test; level=α, tail=:left)
                        @test_throws ArgumentError confint(test; level=α, tail=:both)
                    end
                end
            end
        end
    end

    @testset "consistency" begin
        kernel1 = ExponentialKernel() ∘ ScaleTransform(0.1)
        kernel2 = WhiteKernel()
        kernel = kernel1 ⊗ kernel2
        αs = 0.05:0.1:0.95
        nsamples = 100

        pvalues_consistent = Vector{Float64}(undef, 100)
        pvalues_onlyone = similar(pvalues_consistent)

        for blocksize in (2, 5, 10)
            for nclasses in (2, 10)
                dist = Dirichlet(nclasses, 1)
                predictions = [Vector{Float64}(undef, nclasses) for _ in 1:nsamples]
                targets_consistent = Vector{Int}(undef, nsamples)
                targets_onlyone = ones(Int, nsamples)

                for i in eachindex(pvalues_consistent)
                    # sample predictions and targets
                    for j in 1:nsamples
                        rand!(dist, predictions[j])
                        targets_consistent[j] = rand(Categorical(predictions[j]))
                    end

                    # define test
                    test_consistent = AsymptoticBlockSKCETest(
                        kernel, blocksize, predictions, targets_consistent
                    )
                    test_onlyone = AsymptoticBlockSKCETest(
                        kernel, blocksize, predictions, targets_onlyone
                    )

                    # estimate pvalues
                    pvalues_consistent[i] = pvalue(test_consistent)
                    pvalues_onlyone[i] = pvalue(test_onlyone)
                end

                # compute empirical test errors
                ecdf_consistent = ecdf(pvalues_consistent)
                @test all(ecdf_consistent(α) ≤ α + 0.1 for α in αs)
                @test all(p < 0.05 for p in pvalues_onlyone)
            end
        end
    end
end


text:
module Individual

include("Sampling.jl")
include("SchemaBase.jl")
include("SchemaEvents.jl")

end # module


text:
# f(x)=vcat(x...)
Fvcat(x...)=(vcat(x...),nothing)
Fvcat_inplace(value::Array,auxvalue,x...)=copy!(value,vcat(x...))

function Dvcat(derivativeIDX,f_c,faux_c,grad_c,grad_n,x...)
    startind=1
    for i=1:length(x)
        endind=startind+length(x[i])-1
        if pointer(x[i])==pointer(x[derivativeIDX]) # use pointers so we can deal with vcat([A B A])
            axpy!(1.0,grad_c[startind:endind],grad_n)
        end
        startind=endind+1
    end
end

if PROC=="GPU"

    function Fvcat_inplace(value::CudaArray,auxvalue,x::CudaArray...) # inplace
        totallength=1
        for i in 1:length(x)
            copyinto!(value,x[i],totallength)
            totallength+=length(x[i])
        end
    end

    function Dvcat(derivativeIDX,f_c,faux_c,grad_c,grad_n::CudaArray,x...)
        startind=1
        for i=1:length(x)
            endind=startind+length(x[i])-1
            if pointer(x[i])==pointer(x[derivativeIDX]) # use pointers so we can deal with vcat([A B A])
                copyfrom_update!(grad_n,grad_c,startind)
            end
            startind=endind+1
        end
    end

end

Derivative[Fvcat]=Dvcat # Define dictionary lookup
Inplace[Fvcat]=Fvcat_inplace

import Base.vcat
vcat(n::ADnode)=ADnode(Fvcat,n)

vcat(n::Array{ADnode,1})=n
vcat(n::ArrayADnode)=ADnode(Fvcat,n)

export vcat

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
using PMDPs
using PMDPs.LP
using DrWatson

using POMDPs
using BSON, CSV
using DrWatson
using RandomNumbers.Xorshifts
using Random
using DataFrames
using StaticArrays, Distributions # load

using MCTS, DiscreteValueIteration
using POMDPSimulators # load histories
using POMDPPolicies

using LightGraphs
using GraphPlot
using Cairo, Compose

using Debugger

using TableView

include(srcdir("MDPPricing.jl"))

function prepare_traces(pp::PMDPs.PMDPProblem, pp_params::Dict, vi::Bool, name::String, N::Int64; trace_folder = "test_traces", seed = 1, verbose=false)
    mg = PMDPs.PMDPg(pp)
    rnd = Xorshift128Plus(seed)
    fname = savename("$(name)_N=$(N)", pp_params,  "bson")
    fpath = datadir(trace_folder, fname)

    if isfile(fpath)
        data = PMDPs.load_traces(fpath)
        verbose ? println("Loading $fpath") : nothing
    else
        traces = [PMDPs.simulate_trace(mg, rnd) for i in 1:N]
        data =  @dict(name, pp, pp_params, traces, vi)
        @tagsave(fpath, data)
        verbose ? println("Saving $fpath") : nothing
    end
    return data
end

"""
Get input data
"""

problems = get_fast_benchmarks()
inputs = [prepare_traces(pp, params, vi, name, 100; verbose=true) for (pp, params, vi, name) in problems[1:end] ]

"""
Evaluate
"""
N_sim = 20

dpw_solver_params = (;depth=50,
    exploration_constant=40.0, max_time=1.,
    enable_state_pw = false,
    keep_tree=true, show_progress=false, rng=Xorshift128Plus())

mcts_solver_params = (;depth=50,
    exploration_constant=40.0, max_time=1.,
    rng=Xorshift128Plus())

# pp_params = Dict(pairs((nᵣ=3, c=3, T=10, expected_res=3., res_budget_μ=5., objective=objective)))
# name = "linear_problem"
out_folder="test"

for (i, data) in enumerate(inputs)
    print("\t Data $i - Evaluating $(data[:name]) with $(data[:pp_params]): ")
    print("flatrate..."); PMDPs.process_data(data, PMDPs.flatrate; folder=out_folder, N=N_sim)
    print("hindsight..."); PMDPs.process_data(data, PMDPs.hindsight; folder=out_folder, N=N_sim)
    print("vi..."); data[:vi] && PMDPs.process_data(data, PMDPs.vi; folder=out_folder, N=N_sim)
    # print("vi..."); data[:vi] && PMDPs.process_data(data, PMDPs.fhvi; folder=out_folder, N=N_sim)

    print("dpw..."); PMDPs.process_data(data, PMDPs.mcts; folder=out_folder, N=N_sim,
            method_info="dpw_$(savename(dpw_solver_params))", solver=DPWSolver(;dpw_solver_params...))
    println("mcts..."); PMDPs.process_data(data, PMDPs.mcts; folder=out_folder, N=N_sim,
            method_info="vanilla_$(savename(mcts_solver_params))", solver=MCTSSolver(;mcts_solver_params...))
end
println("LP Done.")

"""
Collect results
"""

results = folder_report(datadir("results", "test", "linear_problem"); raw_result_array=false);
agg_res = format_result_table(results.results; N=N_sim)

# using WebIO
vscodedisplay(agg_res)


text:
"""
The advent of quantum technology requires new figures of merit to compare the performance
of quantum and classical systems.
The signaling dimension quantifies the classical
simulation cost of black-box devices.
That is, it specifies the minimum amount of noiseless, classical communication
needed to exactly simulate the input-output statistics of the device in question.
The signaling dimension is discussed in greater detail in the Background section
of this documentation.

## Features:
* Tools to compute signaling correlations and polytopes (see [Signaling Correlations](@ref) and [Signaling Polytopes](@ref)).
* A catalog of Bell inequalities for witnessing signaling dimension (see [Bell Inequalities](@ref)).
* Methods to certify the signaling dimension of signaling systems (see [Certifying Signaling Dimension](@ref)).
"""
module SignalingDimension

using LinearAlgebra
using Combinatorics
using BellScenario

include("./bell_inequalities.jl")
include("./affinely_independent_enumerations.jl")
include("./certifying_signaling_dimension.jl")
include("./verifications.jl")

end


text:
@testset "Rating Kernel Tests" begin
    preference = Persa.Preference([1, 2, 3, 4, 5])

    @test Persa.value(Persa.MissingRating()) === missing
    @test Persa.value(Persa.MissingRating{Int}()) === missing

    @test Persa.value(Persa.Rating(4, preference)) == 4
    @test Persa.value(Persa.Rating(6, preference)) == Persa.maximum(preference)
    @test Persa.value(Persa.Rating(-1, preference)) == Persa.minimum(preference)

    @test Persa.value(Persa.PredictRating(4, preference)) == 4
    @test Persa.value(Persa.PredictRating(6, preference)) == Persa.maximum(preference)
    @test Persa.value(Persa.PredictRating(-1, preference)) == Persa.minimum(preference)
    @test Persa.value(Persa.PredictRating(5.5, preference)) == Persa.maximum(preference)
    @test Persa.value(Persa.PredictRating(-1.1, preference)) == Persa.minimum(preference)
    @test Persa.value(Persa.PredictRating(4.1, preference)) == 4.1

    predict = Persa.PredictRating(4.5, preference)
    label = Persa.Rating(4, preference)
end


text:
"""
    two_sum(a, b)

Computes `hi = fl(a+b)` and `lo = err(a+b)`, using the algorithm by D. E. Knuth,
"The Art of Computer Programming: Seminumerical Algorithms", 1969.

This algorithm does not use any branch. See also `fast_two_sum` for an
alternative algorithm which branches but does fewer arithmetic operations.
"""
@inline function two_sum(a::T, b::T) where {T}
    @fastmath  hi = a + b
    v  = hi - a
    lo = (a - (hi - v)) + (b - v)
    return hi, lo
end

"""
    fast_two_sum(a, b)

Computes `hi = fl(a+b)` and `lo = err(a+b)`, using the algorithm by
T. J. Dekker, "A Floating-Point Technique for Extending the Available
Precision", 1971

Despite its name, this algorithm may not be as fast as expected on all hardware
architectures because it branches. See also `two_sum` for an alternative
algorithm which performs more arithmetic operations, but does not branch.
"""
@inline function fast_two_sum(a::T, b::T) where T <: Real
    x = a + b

    if abs(a) < abs(b)
        (a_,b_) = (b,a)
    else
        (a_,b_) = (a,b)
    end

    z = x - a_
    e = b_ - z

    x, e
end

@inline function fast_two_sum(a::T, b::T) where T <: AbstractSIMD
    x = a + b

    t = abs(a) < abs(b)
    a_ = ifelse(t, b, a)
    b_ = ifelse(t, a, b)

    z = x - a_
    e = b_ - z

    x, e
end

"""
   three_sum(a, b, c)

Computes `hi = fl(a+b+c)` and `md = err(a+b+c), lo = err(md)`.
"""
function three_sum(a::T,b::T,c::T) where {T}
    s, t   = two_sum(b, c)
    hi, u  = two_sum(a, s)
    md, lo = two_sum(u, t)
    hi, md = two_hilo_sum(hi, md)
    return hi, md, lo
end

"""
    two_sum(a, b, c)

Computes `hi = fl(a+b+c)` and `lo = err(a+b+c)`.
"""
function two_sum(a::T,b::T,c::T) where {T}
    s, t   = two_sum(b, c)
    hi, u  = two_sum(a, s)
    lo     = u + t
    hi, lo = two_hilo_sum(hi, lo)
    return hi, lo
end

"""
    four_sum(a, b, c, d)

Computes `s1 = fl(a+b+c+d)` and `s2 = err(a+b+c+d),  s3 = err(himd), s4 = err(lomd)`.
- Unchecked Precondition: !(isinf(a) | isinf(b) | isinf(c) | isinf(d))
"""
function four_sum(a::T,b::T,c::T,d::T) where {T}
    t0, t1 = two_sum(a,  b)
    t2, t3 = two_sum(c,  d)
    hi, t4 = two_sum(t0, t2)
    t5, lo = two_sum(t1, t3)
    hm, ml = two_sum(t4, t5)
    ml, lo = two_hilo_sum(ml, lo)
    hm, ml = two_hilo_sum(hm, ml)
    hi, hm = two_hilo_sum(hi,hm)
    return hi, hm, ml, lo
end

"""
    three_sum(a, b, c, d)

Computes `hi = fl(a+b+c+d)` and `md = err(a+b+c+d), lo = err(md)`.
"""
function three_sum(a::T,b::T,c::T,d::T) where {T}
    t0, t1 = two_sum(a ,  b)
    t0, t2 = two_sum(t0,  c)
    hi, t3 = two_sum(t0,  d)
    t0, t1 = two_sum(t1, t2)
    hm, t2 = two_sum(t0, t3) # here, t0 >= t3
    ml     = t1 + t2
    return hi, hm, ml
end

"""
    two_sum(a, b, c, d)

Computes `hi = fl(a+b+c+d)` and `lo = err(a+b+c+d)`.
"""
function two_sum(a::T,b::T,c::T,d::T) where {T}
    t0, t1 = two_sum(a ,  b)
    t0, t2 = two_sum(t0,  c)
    hi, t3 = two_sum(t0,  d)
    t0, t1 = two_sum(t1, t2)
    lo     = t0 + t3
    return hi, lo
end

"""
    five_sum(a, b, c, d, e)

Computes `s = fl(a+b+c+d+e)` and
    `e1 = err(a+b+c+d), e2 = err(e1), e3 = err(e2), e4 = err(e3)`.
"""
function five_sum(v::T, w::T, x::T, y::T, z::T) where {T}
    t0, t4 = two_sum(y, z)
    t0, t3 = two_sum(x, t0)
    t0, t2 = two_sum(w, t0)
    a, t1  = two_sum(v, t0)
    t0, t3 = two_sum(t3, t4)
    t0, t2 = two_sum(t2, t0)
    b, t1  = two_sum(t1, t0)
    t0, t2 = two_sum(t2, t3)
    c, t1  = two_sum(t1, t0)
    d, e   = two_sum(t1, t2)
    return a, b, c, d, e
end

"""
    two_diff(a, b)
Computes `s = fl(a-b)` and `e = err(a-b)`.
"""
@inline function two_diff(a::T, b::T) where {T}
    hi = a - b
    a1 = hi + b
    b1 = hi - a1
    lo = (a - a1) - (b + b1)
    return hi, lo
end
"""
    three_diff(a, b, c)

Computes `s = fl(a-b-c)` and `e1 = err(a-b-c), e2 = err(e1)`.
"""
function three_diff(a::T,b::T,c::T) where {T}
    s, t = two_diff(-b, c)
    x, u = two_sum(a, s)
    y, z = two_sum(u, t)
    x, y = two_hilo_sum(x, y)
    return x, y, z
end

"""
    four_diff(a, b, c, d)

Computes `hi = fl(a-b-c-d)` and `hm = err(a-b-c-d), ml = err(hm), lo = err(ml)`.
"""
function four_diff(a::T,b::T,c::T,d::T) where {T}
    t0, t1 = two_diff(a ,  b)
    t0, t2 = two_diff(t0,  c)
    hi, t3 = two_diff(t0,  d)
    t0, t1 = two_sum(t1, t2)
    hm, t2 = two_sum(t0, t3) # here, t0 >= t3
    ml, lo = two_sum(t1, t2)
    return hi, hm, ml, lo
end

"""
    two_square(a)
Computes `s = fl(a*a)` and `e = err(a*a)`.
"""
@inline function two_square(a::T) where {T}
    p = a * a
    e = fma(a, a, -p)
    p, e
end

"""
    two_prod(a, b)
Computes `s = fl(a*b)` and `e = err(a*b)`.
"""
@inline function two_prod(a::T, b::T) where {T}
    p = a * b
    e = fma(a, b, -p)
    p, e
end

"""
    three_prod(a, b, c)

Computes `hi = fl(a*b*c)` and `md = err(a*b*c), lo = err(md)`.
"""
function three_prod(a::T, b::T, c::T) where {T}
    abhi, ablo = two_prod(a, b)
    hi, abhiclo = two_prod(abhi, c)
    ablochi, abloclo = two_prod(ablo, c)
    md, lo, tmp  = three_sum(ablochi, abhiclo, abloclo)
    return hi, md, lo
end

#=
   three_fma algorithm from
   Sylvie Boldo and Jean-Michel Muller
   Some Functions Computable with a Fused-mac
=#

"""
   three_fma(a, b, c)
Computes `s = fl(fma(a,b,c))` and `e1 = err(fma(a,b,c)), e2 = err(e1)`.
"""
function three_fma(a::T, b::T, c::T) where {T}
     x = fma(a, b, c)
     y, z = two_prod(a, b)
     t, z = two_sum(c, z)
     t, u = two_sum(y, t)
     y = ((t - x) + u)
     y, z = two_hilo_sum(y, z)
     return x, y, z
end

# with arguments sorted by magnitude

"""
    two_hilo_sum(a, b)
*unchecked* requirement `|a| ≥ |b|`
Computes `hi = fl(a+b)` and `lo = err(a+b)`.
"""
@inline function two_hilo_sum(a::T, b::T) where {T}
    hi = a + b
    lo = b - (hi - a)
    return hi, lo
end

"""
    two_lohi_sum(a, b)
*unchecked* requirement `|b| ≥ |a|`
Computes `hi = fl(a+b)` and `lo = err(a+b)`.
"""
@inline function two_lohi_sum(a::T, b::T) where {T}
    hi = b + a
    lo = a - (hi - b)
    return hi, lo
end

"""
    two_hilo_diff(a, b)

*unchecked* requirement `|a| ≥ |b|`
Computes `hi = fl(a-b)` and `lo = err(a-b)`.
"""
@inline function two_hilo_diff(a::T, b::T) where {T}
    hi = a - b
    lo = (a - hi) - b
    hi, lo
end

"""
    two_lohi_diff(a, b)

*unchecked* requirement `|b| ≥ |a|`
Computes `hi = fl(a-b)` and `lo = err(a-b)`.
"""
@inline function two_lohi_diff(a::T, b::T) where {T}
    hi = b - a
    lo = (b - hi) - a
    hi, lo
end

"""
    three_hilo_sum(a, b, c)

*unchecked* requirement `|a| ≥ |b| ≥ |c|`
Computes `x = fl(a+b+c)` and `y = err(a+b+c), z = err(y)`.
"""
function three_hilo_sum(a::T,b::T,c::T) where {T}
    s, t = two_hilo_sum(b, c)
    x, u = two_hilo_sum(a, s)
    y, z = two_hilo_sum(u, t)
    x, y = two_hilo_sum(x, y)
    return x, y, z
end

"""
    three_lohi_sum(a, b, c)

*unchecked* requirement `|c| ≥ |b| ≥ |a|`
Computes `x = fl(a+b+c)` and `y = err(a+b+c), z = err(y)`.
"""
function three_lohi_sum(a::T,b::T,c::T) where {T}
    s, t = two_hilo_sum(b, a)
    x, u = two_hilo_sum(c, s)
    y, z = two_hilo_sum(u, t)
    x, y = two_hilo_sum(x, y)
    return x, y, z
end

"""
    three_hilo_diff(a, b, c)

*unchecked* requirement `|a| ≥ |b| ≥ |c|`
Computes `x = fl(a-b-c)` and `y = err(a-b-c), z = err(y)`.
"""
function three_hilo_diff(a::T,b::T,c::T) where {T}
    s, t = two_hilo_diff(b, -c)
    x, u = two_hilo_sum(a, s)
    y, z = two_hilo_sum(u, t)
    x, y = two_hilo_sum(x, y)
    return x, y, z
end

"""
    three_lohi_diff(a, b, c)

*unchecked* requirement `|c| ≥ |b| ≥ |a|`
Computes `x = fl(a-b-c)` and `y = err(a-b-c), z = err(y)`.
"""
function three_lohi_diff(c::T,b::T,a::T) where {T}
    s, t = two_hilo_diff(b, -c)
    x, u = two_hilo_sum(a, s)
    y, z = two_hilo_sum(u, t)
    x, y = two_hilo_sum(x, y)
    return x, y, z
end

"""
    four_hilo_sum(a, b, c, d)

*unchecked* requirement `|a| ≥ |b| ≥ |c| ≥ |d|`
Computes `hi = fl(a+b+c+d)` and `hm = err(a+b+c+d), ml = err(hm), lo = err(ml)`.
"""
function four_hilo_sum(a::T,b::T,c::T,d::T) where {T}
    t0, t1 = two_hilo_sum(a ,  b)
    t0, t2 = two_hilo_sum(t0,  c)
    hi, t3 = two_hilo_sum(t0,  d)
    t0, t1 = two_hilo_sum(t1, t2)
    hm, t2 = two_hilo_sum(t0, t3) # here, t0 >= t3
    ml, lo = two_hilo_sum(t1, t2)
    return hi, hm, ml, lo
end

"""
    four_lohi_sum(a, b, c, d)

*unchecked* requirement `|d| ≥ |c| ≥ |b| ≥ |a|`
Computes `hi = fl(a+b+c+d)` and `hm = err(a+b+c+d), ml = err(hm), lo = err(ml)`.
"""
function four_lohi_sum(d::T,c::T,b::T,a::T) where {T}
    t0, t1 = two_hilo_sum(a ,  b)
    t0, t2 = two_hilo_sum(t0,  c)
    hi, t3 = two_hilo_sum(t0,  d)
    t0, t1 = two_hilo_sum(t1, t2)
    hm, t2 = two_hilo_sum(t0, t3)
    ml, lo = two_hilo_sum(t1, t2)
    return hi, hm, ml, lo
end

"""
    four_hilo_diff(a, b, c, d)

*unchecked* requirement `|a| ≥ |b| ≥ |c| ≥ |d|`
Computes `hi = fl(a-b-c-d)` and `hm = err(a-b-c-d), ml = err(hm), lo = err(ml)`.
"""
function four_hilo_diff(a::T,b::T,c::T,d::T) where {T}
    t0, t1 = two_hilo_diff(a,  b)
    t0, t2 = two_hilo_diff(t0,  c)
    hi, t3 = two_hilo_diff(t0,  d)
    t0, t1 = two_hilo_sum(t1, t2)
    hm, t2 = two_hilo_sum(t0, t3) # here, t0 >= t3
    ml, lo = two_hilo_sum(t1, t2)
    return hi, hm, ml, lo
end

"""
    four_hilo_diff(a, b, c, d)

*unchecked* requirement `|d| ≥ |c| ≥ |b| ≥ |a|`
Computes `hi = fl(a-b-c-d)` and `hm = err(a-b-c-d), ml = err(hm), lo = err(ml)`.
"""
function four_lohi_diff(d::T,c::T,b::T,a::T) where {T}
    t0, t1 = two_hilo_diff(a,  b)
    t0, t2 = two_hilo_diff(t0,  c)
    hi, t3 = two_hilo_diff(t0,  d)
    t0, t1 = two_hilo_sum(t1, t2)
    hm, t2 = two_hilo_sum(t0, t3) # here, t0 >= t3
    ml, lo = two_hilo_sum(t1, t2)
    return hi, hm, ml, lo
end


text:
using Iris
using Test

@testset "VectorFields" begin
    @testset "VectorField" begin
        N = 10 # dimension
        M = 2 # number of components
        n = 15 # number of sites
        nmodes = 3
        r = rand(ComplexF64,M*n,nmodes)
        pos = Vector{Point{N}}(undef,n)
        for i ∈ 1:n pos[i] = Point(randn(Float64,N)...) end
        e = VectorField{M}(pos,r)
        @test e==conj(conj(e))
        @test e + e == 2e
        @test e - e == 0e
        @test 4e/2==2e
        @test firstindex(e)==1
        @test lastindex(e)==M*n*nmodes
        @test length(e) == M*n*nmodes
        @test size(e) == (M*n,nmodes)
        @test eltype(e) == ComplexF64
        @test typeof(e(2))<:VectorField
        @test_nowarn VectorField(e,rand(ComplexF64,size(e)...))
        @test_nowarn update!(e,rand(ComplexF64,size(e)...))
        @test_nowarn component(1,e)
        @test_throws BoundsError component(3,e; verbose=false)
    end
    @testset "ScalarField" begin
        N = 1 # dimension
        M = 1 # number of components
        n = 25 # number of sites
        nmodes = 4
        r = rand(ComplexF64,M*n,nmodes)
        pos = Vector{Point{N}}(undef,n)
        for i ∈ 1:n pos[i] = Point(randn(Float64,N)...) end
        e = ScalarField(pos,r)
        @test e==conj(conj(e))
        @test e + e == 2e
        @test e - e == 0e
        @test 4e/2==2e
        @test firstindex(e)==1
        @test lastindex(e)==M*n*nmodes
        @test length(e) == M*n*nmodes
        @test size(e) == (M*n,nmodes)
        @test eltype(e) == ComplexF64
        @test typeof(e(2))<:ScalarField
        @test_nowarn ScalarField(e,rand(ComplexF64,size(e)...))
        @test_nowarn update!(e,rand(ComplexF64,size(e)...))
        @test_nowarn component(1,e)
        @test_throws BoundsError component(3,e; verbose=false)
    end
    @testset "ElectricField" begin
        N = 1 # dimension
        M = 3 # number of components
        n = 25 # number of sites
        nmodes = 4
        r = rand(ComplexF64,M*n,nmodes)
        pos = Vector{Point{N}}(undef,n)
        for i ∈ 1:n pos[i] = Point(randn(Float64,N)...) end
        e = ElectricField(pos,r)
        @test e==conj(conj(e))
        @test e + e == 2e
        @test e - e == 0e
        @test 4e/2==2e
        @test firstindex(e)==1
        @test lastindex(e)==M*n*nmodes
        @test length(e) == M*n*nmodes
        @test size(e) == (M*n,nmodes)
        @test eltype(e) == ComplexF64
        @test typeof(e(2))<:ElectricField
        @test_nowarn ElectricField(e,rand(ComplexF64,size(e)...))
        @test_nowarn update!(e,rand(ComplexF64,size(e)...))
        @test_nowarn component(1,e)
        @test_throws BoundsError component(4,e; verbose=false)
    end
end


text:
# Noncentral hypergeometric distribution
# TODO: this distribution needs clean-up and testing

@compat abstract type NoncentralHypergeometric{T<:Real} <: DiscreteUnivariateDistribution end

### handling support

function insupport(d::NoncentralHypergeometric, x::Real)
    isinteger(x) && minimum(d) <= x <= maximum(d)
end

@distr_support NoncentralHypergeometric max(d.n - d.nf, 0) min(d.ns, d.n)

# Functions

function quantile{T<:Real}(d::NoncentralHypergeometric{T}, q::Real)
    if !(zero(q) <= q <= one(q))
        T(NaN)
    else
        range = support(d)
        if q > 1/2
            q = 1 - q
            range = reverse(range)
        end

        qsum, i = zero(T), 0
        while qsum < q
            i += 1
            qsum += pdf(d, range[i])
        end
        range[i]
    end
end

params(d::NoncentralHypergeometric) = (d.ns, d.nf, d.n, d.ω)
@inline partype{T<:Real}(d::NoncentralHypergeometric{T}) = T

## Fisher's noncentral hypergeometric distribution

immutable FisherNoncentralHypergeometric{T<:Real} <: NoncentralHypergeometric{T}
    ns::Int    # number of successes in population
    nf::Int    # number of failures in population
    n::Int     # sample size
    ω::T # odds ratio

    function (::Type{FisherNoncentralHypergeometric{T}}){T}(ns::Real, nf::Real, n::Real, ω::T)
        @check_args(FisherNoncentralHypergeometric, ns >= zero(ns) && nf >= zero(nf))
        @check_args(FisherNoncentralHypergeometric, zero(n) < n < ns + nf)
        @check_args(FisherNoncentralHypergeometric, ω > zero(ω))
        new{T}(ns, nf, n, ω)
    end
end

FisherNoncentralHypergeometric{T<:Real}(ns::Integer, nf::Integer, n::Integer, ω::T) = FisherNoncentralHypergeometric{T}(ns, nf, n, ω)

FisherNoncentralHypergeometric(ns::Integer, nf::Integer, n::Integer, ω::Integer) = FisherNoncentralHypergeometric(ns, nf, n, Float64(ω))

# Conversions
convert{T<:Real}(::Type{FisherNoncentralHypergeometric{T}}, ns::Real, nf::Real, n::Real, ω::Real) = FisherNoncentralHypergeometric(ns, nf, n, T(ω))
convert{T<:Real, S<:Real}(::Type{FisherNoncentralHypergeometric{T}}, d::FisherNoncentralHypergeometric{S}) = FisherNoncentralHypergeometric(d.ns, d.nf, d.n, T(d.ω))

# Properties
@compat function _P(d::FisherNoncentralHypergeometric, k::Int)
    y = support(d)
    p = -log.(d.ns + 1) - lbeta.(d.ns - y + 1, y + 1) -
            log.(d.nf + 1) - lbeta.(d.nf - d.n + y + 1, d.n - y + 1) +
            xlogy.(y, d.ω) + xlogy.(k, y)
    logsumexp(p)
end

function _mode(d::FisherNoncentralHypergeometric)
    A = d.ω - 1
    B = d.n - d.nf - (d.ns + d.n + 2)*d.ω
    C = (d.ns + 1)*(d.n + 1)*d.ω
    -2C / (B - sqrt(B^2-4A*C))
end

mean(d::FisherNoncentralHypergeometric) = @compat(exp.(_P(d,1) - _P(d,0)))
@compat var(d::FisherNoncentralHypergeometric) = exp.(_P(d,2) - _P(d,0)) - exp.(2*(_P(d,1) - _P(d,0)))
mode(d::FisherNoncentralHypergeometric) = floor(Int, _mode(d))

testfd(d::FisherNoncentralHypergeometric) = d.ω^3

logpdf(d::FisherNoncentralHypergeometric, k::Int) =
    -log(d.ns + 1) - lbeta(d.ns - k + 1, k + 1) -
    log(d.nf + 1) - lbeta(d.nf - d.n + k + 1, d.n - k + 1) +
    xlogy(k, d.ω) - _P(d, 0)

pdf(d::FisherNoncentralHypergeometric, k::Int) = exp(logpdf(d, k))

## Wallenius' noncentral hypergeometric distribution

immutable WalleniusNoncentralHypergeometric{T<:Real} <: NoncentralHypergeometric{T}
    ns::Int    # number of successes in population
    nf::Int    # number of failures in population
    n::Int     # sample size
    ω::T # odds ratio

    function (::Type{WalleniusNoncentralHypergeometric{T}}){T}(ns::Real, nf::Real, n::Real, ω::T)
        @check_args(WalleniusNoncentralHypergeometric, ns >= zero(ns) && nf >= zero(nf))
        @check_args(WalleniusNoncentralHypergeometric, zero(n) < n < ns + nf)
        @check_args(WalleniusNoncentralHypergeometric, ω > zero(ω))
        new{T}(ns, nf, n, ω)
    end
end

WalleniusNoncentralHypergeometric{T<:Real}(ns::Integer, nf::Integer, n::Integer, ω::T) = WalleniusNoncentralHypergeometric{T}(ns, nf, n, ω)

WalleniusNoncentralHypergeometric(ns::Integer, nf::Integer, n::Integer, ω::Integer) = WalleniusNoncentralHypergeometric(ns, nf, n, Float64(ω))

# Conversions
convert{T<:Real}(::Type{WalleniusNoncentralHypergeometric{T}}, ns::Real, nf::Real, n::Real, ω::Real) = WalleniusNoncentralHypergeometric(ns, nf, n, T(ω))
convert{T<:Real, S<:Real}(::Type{WalleniusNoncentralHypergeometric{T}}, d::WalleniusNoncentralHypergeometric{S}) = WalleniusNoncentralHypergeometric(d.ns, d.nf, d.n, T(d.ω))

# Properties
mean(d::WalleniusNoncentralHypergeometric) = sum(support(d) .* pdf(d, support(d)))
var(d::WalleniusNoncentralHypergeometric)  = sum((support(d) - mean(d)).^2 .* pdf(d, support(d)))
mode(d::WalleniusNoncentralHypergeometric) = support(d)[indmax(pdf(d, support(d)))]

entropy(d::WalleniusNoncentralHypergeometric) = 1

testfd(d::WalleniusNoncentralHypergeometric) = d.ω^3

function logpdf(d::WalleniusNoncentralHypergeometric, k::Int)
    D = d.ω * (d.ns - k) + (d.nf - d.n + k)
    f(t) = (1 - t^(d.ω / D))^k * (1 - t^(1 / D))^(d.n - k)
    I, _ = quadgk(f, 0, 1)
    return -log(d.ns + 1) - lbeta(d.ns - k + 1, k + 1) -
    log(d.nf + 1) - lbeta(d.nf - d.n + k + 1, d.n - k + 1) + log(I)
end

pdf(d::WalleniusNoncentralHypergeometric, k::Int) = exp(logpdf(d, k))