text:
#
# formats.jl --
#
# Define OI-FITS data-blocks.
#
#------------------------------------------------------------------------------

# OI_TARGET definition (1st revision):
define("OI_TARGET", 1,
       ["OI_REVN    I      revision number of the table definition",
        "---------------------------------------------------------",
        "TARGET_ID  I(1)   index number",
        "TARGET     A(16)  target name",
        "RAEP0      D(1)   RA at mean equinox [deg]",
        "DECEP0     D(1)   DEC at mean equinox [deg]",
        "EQUINOX    E(1)   equinox [yr]",
        "RA_ERR     D(1)   error in RA at mean equinox [deg]",
        "DEC_ERR    D(1)   error in DEC at mean equino [deg]",
        "SYSVEL     D(1)   systemic radial velocity [m/s]",
        "VELTYP     A(8)   reference for radial velocity",
        "VELDEF     A(8)   definition of radial velocity",
        "PMRA       D(1)   proper motion in RA [deg/yr]",
        "PMDEC      D(1)   proper motion in DEC [deg/yr]",
        "PMRA_ERR   D(1)   error of proper motion in RA [deg/yr]",
        "PMDEC_ERR  D(1)   error of proper motion in DEC [deg/yr]",
        "PARALLAX   E(1)   parallax [deg]",
        "PARA_ERR   E(1)   error in parallax [deg]",
        "SPECTYP    A(16)  spectral type"])

# OI_TARGET definition (2nd revision):
define("OI_TARGET", 2,
       ["OI_REVN    I      revision number of the table definition",
        "---------------------------------------------------------",
        "TARGET_ID  I(1)   index number",
        "TARGET     A(16)  target name",
        "RAEP0      D(1)   RA at mean equinox [deg]",
        "DECEP0     D(1)   DEC at mean equinox [deg]",
        "EQUINOX    E(1)   equinox [yr]",
        "RA_ERR     D(1)   error in RA at mean equinox [deg]",
        "DEC_ERR    D(1)   error in DEC at mean equino [deg]",
        "SYSVEL     D(1)   systemic radial velocity [m/s]",
        "VELTYP     A(8)   reference for radial velocity",
        "VELDEF     A(8)   definition of radial velocity",
        "PMRA       D(1)   proper motion in RA [deg/yr]",
        "PMDEC      D(1)   proper motion in DEC [deg/yr]",
        "PMRA_ERR   D(1)   error of proper motion in RA [deg/yr]",
        "PMDEC_ERR  D(1)   error of proper motion in DEC [deg/yr]",
        "PARALLAX   E(1)   parallax [deg]",
        "PARA_ERR   E(1)   error in parallax [deg]",
        "SPECTYP    A(16)  spectral type",
        "CATEGORY  ?A(3)   CALibrator or SCIence target"])

# OI_ARRAY definition (1st revision):
define("OI_ARRAY", 1,
       ["OI_REVN     I     revision number of the table definition",
        "ARRNAME     A     array name for cross-referencing",
        "FRAME       A     coordinate frame",
        "ARRAYX      D     array center X-coordinate [m]",
        "ARRAYY      D     array center Y-coordinate [m]",
        "ARRAYZ      D     array center Z-coordinate [m]",
        "------------------------------------------------------------",
        "TEL_NAME   A(16)  telescope name",
        "STA_NAME   A(16)  station name",
        "STA_INDEX  I(1)   station index",
        "DIAMETER   E(1)   element diameter [m]",
        "STAXYZ     D(3)   station coordinates relative to array center [m]"])

# OI_ARRAY definition (2nd revision):
define("OI_ARRAY", 2,
       ["OI_REVN     I     revision number of the table definition",
        "ARRNAME     A     array name for cross-referencing",
        "FRAME       A     coordinate frame",
        "ARRAYX      D     array center X-coordinate [m]",
        "ARRAYY      D     array center Y-coordinate [m]",
        "ARRAYZ      D     array center Z-coordinate [m]",
        "------------------------------------------------------------",
        "TEL_NAME   A(16)  telescope name",
        "STA_NAME   A(16)  station name",
        "STA_INDEX  I(1)   station index",
        "DIAMETER   E(1)   element diameter [m]",
        "STAXYZ     D(3)   station coordinates relative to array center [m]",
        "FOV        D(1)   photometric field of view [arcsec]",
        "FOVTYPE    A(6)   model for FOV: 'FWHM' or 'RADIUS'"])

# OI_WAVELENGTH definition (1st and 2nd revisions):
for revn in (1, 2)
    define("OI_WAVELENGTH", revn,
           ["OI_REVN     I     revision number of the table definition",
            "INSNAME     A     name of detector for cross-referencing",
            "------------------------------------------------------------",
            "EFF_WAVE   E(1)   effective wavelength of channel [m]",
            "EFF_BAND   E(1)   effective bandpass of channel [m]"])
end

# OI_VIS definition (1st revision):
define("OI_VIS", 1,
       ["OI_REVN     I    revision number of the table definition",
        "DATE-OBS    A    UTC start date of observations",
        "ARRNAME    ?A    name of corresponding array",
        "INSNAME     A    name of corresponding detector",
        "------------------------------------------------------------",
        "TARGET_ID  I(1)  target number as index into OI_TARGET table",
        "TIME       D(1)  UTC time of observation [s]",
        "MJD        D(1)  modified Julian Day [day]",
        "INT_TIME   D(1)  integration time [s]",
        "VISAMP     D(W)  visibility amplitude",
        "VISAMPERR  D(W)  error in visibility amplitude",
        "VISPHI     D(W)  visibility phase [deg]",
        "VISPHIERR  D(W)  error in visibility phase [deg]",
        "UCOORD     D(1)  U coordinate of the data [m]",
        "VCOORD     D(1)  V coordinate of the data [m]",
        "STA_INDEX  I(2)  station numbers contributing to the data",
        "FLAG       L(W)  flag"])

# FIXME: In OI-FITS Rev. 2 only MJD and DATE-OBS must be used to express
#        time. The TIME column is retained only for backwards compatibility
#        and must contain zeros.

# OI_VIS definition (2nd revision):
define("OI_VIS", 2,
       ["OI_REVN          I       revision number of the table definition",
        "DATE-OBS         A       UTC start date of observations",
        "ARRNAME          A       name of corresponding array",
        "INSNAME          A       name of corresponding detector",
        "CORRNAME        ?A       name of corresponding correlation table",
        "AMPTYP          ?A       'ABSOLUTE', 'DIFFERENTIAL' or 'CORRELATED FLUX'",
        "PHITYP          ?A       'ABSOLUTE', or 'DIFFERENTIAL'",
        "AMPORDER        ?I       polynomial fit order for differential amplitudes",
        "PHIORDER        ?I       polynomial fit order for differential phases",
        "------------------------------------------------------------------------",
        "TARGET_ID        I(1)    target number as index into OI_TARGET table",
        "TIME             D(1)    UTC time of observation [s]",
        "MJD              D(1)    modified Julian Day [day]",
        "INT_TIME         D(1)    integration time [s]",
        "VISAMP           D(W)    visibility amplitude",
        "VISAMPERR        D(W)    error in visibility amplitude",
        "CORRINDX_VISAMP ?J(1)    index into correlation matrix for 1st VISAMP element",
        "VISPHI           D(W)    visibility phase [deg]",
        "VISPHIERR        D(W)    error in visibility phase [deg]",
        "CORRINDX_VISPHI ?J(1)    index into correlation matrix for 1st VISPHI element",
        "VISREFMAP       ?L(W,W)  true where spectral channels were taken as reference for differential visibility computation",
        "RVIS            ?D(W)    real part of complex coherent flux",
        "RVISERR         ?D(W)    error on RVIS",
        "CORRINDX_RVIS   ?J(1)    index into correlation matrix for 1st RVIS element",
        "IVIS            ?D(W)    imaginary part of complex coherent flux",
        "IVISERR         ?D(W)    error on IVIS",
        "CORRINDX_IVIS   ?J(1)    index into correlation matrix for 1st IVIS element",
        "UCOORD           D(1)    U coordinate of the data [m]",
        "VCOORD           D(1)    V coordinate of the data [m]",
        "STA_INDEX        I(2)    station numbers contributing to the data",
        "FLAG             L(W)    flag"])

# OI_VIS2 definition (1st revision):
define("OI_VIS2", 1,
       ["OI_REVN     I    revision number of the table definition",
        "DATE-OBS    A    UTC start date of observations",
        "ARRNAME    ?A    name of corresponding array",
        "INSNAME     A    name of corresponding detector",
        "------------------------------------------------------------",
        "TARGET_ID  I(1)  target number as index into OI_TARGET table",
        "TIME       D(1)  UTC time of observation [s]",
        "MJD        D(1)  modified Julian Day [day]",
        "INT_TIME   D(1)  integration time [s]",
        "VIS2DATA   D(W)  squared visibility",
        "VIS2ERR    D(W)  error in squared visibility",
        "UCOORD     D(1)  U coordinate of the data [m]",
        "VCOORD     D(1)  V coordinate of the data [m]",
        "STA_INDEX  I(2)  station numbers contributing to the data",
        "FLAG       L(W)  flag"])

# OI_VIS2 definition (2nd revision):
define("OI_VIS2", 2,
       ["OI_REVN            I    revision number of the table definition",
        "DATE-OBS           A    UTC start date of observations",
        "ARRNAME            A    name of corresponding array",
        "INSNAME            A    name of corresponding detector",
        "CORRNAME          ?A    name of corresponding correlation table",
        "------------------------------------------------------------",
        "TARGET_ID          I(1)  target number as index into OI_TARGET table",
        "TIME               D(1)  UTC time of observation [s]",
        "MJD                D(1)  modified Julian Day [day]",
        "INT_TIME           D(1)  integration time [s]",
        "VIS2DATA           D(W)  squared visibility",
        "VIS2ERR            D(W)  error in squared visibility",
        "CORRINDX_VIS2DATA ?J(1)  index into correlation matrix for 1st VIS2DATA element",
        "UCOORD             D(1)  U coordinate of the data [m]",
        "VCOORD             D(1)  V coordinate of the data [m]",
        "STA_INDEX          I(2)  station numbers contributing to the data",
        "FLAG               L(W)  flag"])

# OI_T3 definition (1st revision):
define("OI_T3", 1,
       ["OI_REVN     I    revision number of the table definition",
        "DATE-OBS    A    UTC start date of observations",
        "ARRNAME    ?A    name of corresponding array",
        "INSNAME     A    name of corresponding detector",
        "------------------------------------------------------------",
        "TARGET_ID  I(1)  target number as index into OI_TARGET table",
        "TIME       D(1)  UTC time of observation [s]",
        "MJD        D(1)  modified Julian Day [day]",
        "INT_TIME   D(1)  integration time [s]",
        "T3AMP      D(W)  triple product amplitude",
        "T3AMPERR   D(W)  error in triple product amplitude",
        "T3PHI      D(W)  triple product phase [deg]",
        "T3PHIERR   D(W)  error in triple product phase [deg]",
        "U1COORD    D(1)  U coordinate of baseline AB of the triangle [m]",
        "V1COORD    D(1)  V coordinate of baseline AB of the triangle [m]",
        "U2COORD    D(1)  U coordinate of baseline BC of the triangle [m]",
        "V2COORD    D(1)  V coordinate of baseline BC of the triangle [m]",
        "STA_INDEX  I(3)  station numbers contributing to the data",
        "FLAG       L(W)  flag"])

# OI_T3 definition (2nd revision):
define("OI_T3", 2,
       ["OI_REVN         I     revision number of the table definition",
        "DATE-OBS        A     UTC start date of observations",
        "ARRNAME         A     name of corresponding array",
        "INSNAME         A     name of corresponding detector",
        "CORRNAME       ?A     name of corresponding correlation table",
        "------------------------------------------------------------",
        "TARGET_ID       I(1)  target number as index into OI_TARGET table",
        "TIME            D(1)  UTC time of observation [s]",
        "MJD             D(1)  modified Julian Day [day]",
        "INT_TIME        D(1)  integration time [s]",
        "T3AMP           D(W)  triple product amplitude",
        "T3AMPERR        D(W)  error in triple product amplitude",
        "CORRINDX_T3AMP ?J(1)  index into correlation matrix for 1st T3AMP element",
        "T3PHI           D(W)  triple product phase [deg]",
        "T3PHIERR        D(W)  error in triple product phase [deg]",
        "CORRINDX_T3PHI ?J(1)  index into correlation matrix for 1st T3PHI element",
        "U1COORD         D(1)  U coordinate of baseline AB of the triangle [m]",
        "V1COORD         D(1)  V coordinate of baseline AB of the triangle [m]",
        "U2COORD         D(1)  U coordinate of baseline BC of the triangle [m]",
        "V2COORD         D(1)  V coordinate of baseline BC of the triangle [m]",
        "STA_INDEX       I(3)  station numbers contributing to the data",
        "FLAG            L(W)  flag"])

# OI_FLUX definition (1st revision):
define("OI_FLUX", 1,
       ["OI_REVN            I     revision number of the table definition",
        "DATE-OBS           A     UTC start date of observations",
        "INSNAME            A     name of corresponding detector",
        "ARRNAME           ?A     name of corresponding array",
        "CORRNAME          ?A     name of corresponding correlation table",
        "FOV               ?D     area of sky over which flux is integrated [arcsec]",
        "FOVTYPE           ?A     model for FOV: 'FWHM' or 'RADIUS'",
        "CALSTAT            A     'C': spectrum is calibrated, 'U': uncalibrated",
        "------------------------------------------------------------",
        "TARGET_ID          I(1)  target number as index into OI_TARGET table",
        "MJD                D(1)  modified Julian Day [day]",
        "INT_TIME           D(1)  integration time [s]",
        "FLUXDATA           D(W)  flux",
        "FLUXERR            D(W)  flux error",
        "CORRINDX_FLUXDATA ?J(1)  index into correlation matrix for 1st FLUXDATA element",
        "STA_INDEX         ?I(1)  station number contributing to the data",
        "FLAG               L(W)  flag"])

# OI_CORR definition (1st revision):
define("OI_CORR", 1,
       ["OI_REVN            I     revision number of the table definition",
        "CORRNAME           A     name of correlation data set",
        "NDATA              I     number of correlated data",
        "------------------------------------------------------------",
        "IINDX              J(1)  1st index of correlation matrix element",
        "JINDX              J(1)  2nd index of correlation matrix element",
        "CORR               D(1)  matrix element"])

# OI_INSPOL definition (1st revision):
define("OI_INSPOL", 1,
       ["OI_REVN    I     revision number of the table definition",
        "DATE-OBS   A     UTC start date of observations",
        "NPOL       I     number of polarization types in this table",
        "ARRNAME    A     identifies corresponding OI_ARRAY",
        "ORIENT     A     orientation of the Jones Matrix: 'NORTH' (for on-sky orientation), or 'LABORATORY'",
        "MODEL      A     describe the way the Jones matrix is estimated",
        "---------------------------------------------------------------------------------------",
        "TARGET_ID  I(1)  target number as index into OI_TARGET table",
        "INSNAME    A(16) INSNAME of this polarization",
        "MJD_OBS    D(1)  modified Julian day, start of time lapse",
        "MJD_END    D(1)  modified Julian day, end of time lapse",
        "JXX        C(W)  complex Jones Matrix component along X axis",
        "JYY        C(W)  complex Jones Matrix component along Y axis",
        "JXY        C(W)  complex Jones Matrix component between Y and X axis",
        "JYX        C(W)  complex Jones Matrix component between Y and X axis",
        "STA_INDEX ?I(1)  station number for the above matrices"])


text:
function init()
    envcp = get(ENV, "CLASSPATH", "")
    sparkjlassembly = joinpath(dirname(@__FILE__), "..", "jvm", "sparkjl", "target", "sparkjl-0.1-assembly.jar")
    classpath = @windows? "$envcp;$sparkjlassembly" : "$envcp:$sparkjlassembly"
    try
        # prevent exceptions in REPL on code reloading
        JavaCall.init(["-ea", "-Xmx1024M", "-Djava.class.path=$classpath"])
    end
end

init()

ext: jl


text:
using LinearAlgebra.BLAS

function BLAS.axpy!(a::Single32, X::T, Y::T) where {T<:AbstractArray}
    axpy!(reinterpret(Float64,a), X, Y)
    return reinterpret(Single32,Y)
end


text:
const TupType = Tuple{Float64, Float64, Float64, Float128, Float128}

struct bfgsObjFunc{pType<:AbstractParams}
	params::pType
	g_cash::G_Cash
    psamples::DArray{TupType, 1, Vector{TupType}}
end

struct bfgsObjFunc_repul{pType<:AbstractParams, GT<:AbstractRepulsion}
	params::pType
	g_cash::G_Cash
    psamples::DArray{TupType, 1, Vector{TupType}}
    grep::GT
end

function (fg!::bfgsObjFunc)(F,G, bs)
	if !(G==nothing)
		ΔF, ∂F = precompute_step!(fg!.g_cash, bs, fg!.params, fg!.psamples)
		G.=∂F
		if !(F==nothing)
			return ΔF
		end
	else
		if !(F==nothing)
			return precompute_step(bs, fg!.params, fg!.psamples)
		end
	end
end

function (fg!::bfgsObjFunc_repul)(F,G, bs)
	if !(G==nothing)
		ΔF, ∂F = precompute_step!(fg!.g_cash, bs, fg!.params, fg!.psamples)
        f, ∂f = compute_repulsion!(fg!.grep, bs)
        ΔF += f
        ∂F[1:length(∂f)] .+= ∂f
		G.=∂F
		if !(F==nothing)
			return ΔF
		end
	else
		if !(F==nothing)
            return precompute_step(bs, fg!.params, fg!.psamples) + compute_repulsion(fg!.grep, bs)
		end
	end
end

function construct_objective(params::AbstractParams, psamples::DArray{TupType, 1, Vector{TupType}}, bs0)
    @assert get_length(params)==length(bs0)
	fg! = bfgsObjFunc(params, G_Cash(params), psamples)
    return OnceDifferentiable(only_fg!(fg!), bs0)
end

function construct_objective(params::AbstractParams, psamples::DArray{TupType, 1, Vector{TupType}}, bs0, ω₀::Float64)
    @assert get_length(params)==length(bs0)
    fg! = bfgsObjFunc_repul(params, G_Cash(params), psamples, GaugeRepulsion2(params, ω₀))
    return OnceDifferentiable(only_fg!(fg!), bs0)
end
function construct_objective(params::AbstractParams, psamples::DArray{TupType, 1, Vector{TupType}}, bs0, reptyp::AbstractRepulsionType)
    @assert get_length(params)==length(bs0)
    fg! = bfgsObjFunc_repul(params, G_Cash(params), psamples, construct_repulsion(reptyp, params))
    return OnceDifferentiable(only_fg!(fg!), bs0)
end


text:
@diffusion_process FitzHughDiffusionConjugAuxSimple begin
    :dimensions
    process --> 2
    wiener --> 1

    :parameters
    (ϵ, s, γ, β, σ, t, T) --> Float64
    (u, v) --> ℝ{2}

    :additional
    constdiff --> true
    linear --> true
end

B(t, P::FitzHughDiffusionConjugAuxSimple) = @SMatrix [0.0  1.0; 0.0 0.0]

β(t, P::FitzHughDiffusionConjugAuxSimple) = ℝ{2}(0.0, 0.0)

σ(t, P::FitzHughDiffusionAux) = ℝ{2}(0.0, P.σ)


text:
function transfer_learn(model; opt=Descent(0.01), epochs=30_000,
                        train_batch_size=50, eval_batch_size=10, eval_interval=1000)

    weights = params(model)
    dist = Uniform(-5, 5)
    testx = Float32.(range(-5, stop=5, length=50))

    for i in 1:epochs
        task = SineWave()
        x = Float32.(rand(dist, train_batch_size))
        y = task(x)
        l = Flux.mse(model(x'), y')
        Flux.back!(l)
        Flux.Optimise._update_params!(opt, weights)

        if i % eval_interval == 0
            @printf("Iteration %d, evaluating model on random task...\n", i)
            evalx = Float32.(rand(dist, eval_batch_size))
            eval_model(model, evalx, testx, SineWave())
        end
    end
end

max_stars_repo_name: JuliaTagBot/MetaLearning

max_issues_repo_name: JuliaTagBot/MetaLearning

max_forks_repo_name: JuliaTagBot/MetaLearning


text:
using ArgParse
using JLD,MAT,JSON

include("util.jl")

function main(args)
    s = ArgParseSettings()
    s.description = "CNN-HMM hybrid ASR system in Knet by Ilker Kesen."

    @add_arg_table s begin
        ("--jsonfile"; required=true)
        ("--savefile"; required=true)
    end

    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true); println(o); flush(STDOUT)

    jsondata = JSON.parsefile(abspath(o[:jsonfile]))
    p2i, i2p = build_vocab(jsondata)

    posteriors, priors = get_probabilities(jsondata, p2i, i2p)
    save(o[:savefile], "posteriors", posteriors, "priors", priors)
end

!isinteractive() && !isdefined(Core.Main, :load_only) && main(ARGS)

max_stars_repo_name: ilkerkesen/phoneme-convnet

max_issues_repo_name: ilkerkesen/phoneme-convnet

max_forks_repo_name: ilkerkesen/phoneme-convnet


text:
struct CellMap
    yL::Any
    yR::Any
    jacobian::Any
    dim::Any
    function CellMap(yL, yR)
        dim = length(yL)
        @assert length(yR) == dim
        @assert all(yR .>= yL)
        jacobian = 0.5 * (yR - yL)
        new(yL, yR, jacobian, dim)
    end
end

function Base.show(io::IO, cellmap::CellMap)
    dim = cellmap.dim
    yL = cellmap.yL
    yR = cellmap.yR
    str = "CellMap\n\tLeft : $yL\n\tRight: $yR"
    print(io, str)
end

function dimension(C::CellMap)
    return C.dim
end

function jacobian(C::CellMap)
    return C.jacobian
end

function inverse_jacobian(C::CellMap)
    return 1.0 ./ jacobian(C)
end

function determinant_jacobian(C)
    return prod(jacobian(C))
end

function face_determinant_jacobian(C)
    jac = jacobian(C)
    return [jac[1], jac[2], jac[1], jac[2]]
end

function (C::CellMap)(x)
    dim = dimension(C)
    return C.yL .+ (jacobian(C) .* (x .+ ones(dim)))
end

function inverse(c::CellMap,x)
    jac = jacobian(c)
    return ((x .- c.yL) ./ jac) .- 1.0
end


text:
"""
    embed_system(m;name)

Embeds system `m` as the single component of a larger system.
This helps to match the naming of the states, parameters, observables to
the namespace of the system.

```jldocstest; output=false
using ModelingToolkit, DifferentialEquations
using MTKHelpers
function samplesystem(;name)
    @variables t
    D = Differential(t)
    sts = @variables x(t) RHS(t)  # RHS is observed
    ps = @parameters τ       # parameters
    ODESystem([ RHS  ~ (1 - x)/τ, D(x) ~ RHS ], t, sts, ps; name)
end
@named m = samplesystem()
@named sys = embed_system(m)

# note that keys are `m.x`,`m.τ` or `m.RHS`.
# Hence only m needs to be defined rather than all the states, parameters,
# and observables of the system.
prob = ODEProblem(sys, [m.x => 0.0], (0.0,10.0), [m.τ => 3.0])
sol = solve(prob);
dx1 = sol[m.RHS][1:3]
dx2 = sol[getproperty(m,:RHS)][1:3]  # access by symbols

# using Plots
# plot(sol, vars=[m.x,m.RHS])
# plot(sol, vars=getproperty.(Ref(m),[:x, :RHS])) # access by symbols
length(dx2) == 3
# output
true
```
"""
function embed_system(m;name, simplify=true)
    @named _sys_embed = ODESystem(Equation[], ModelingToolkit.get_iv(m))
    sys = compose(_sys_embed, [m]; name)
    if simplify
        sys = structural_simplify(sys)
    end
    sys
end

"""
    symbol(t)

Extract the inner symbol from a Term or Num object.
"""
function symbol(t::Term); symbol(t.f); end
symbol(num::Num) = symbol(num.val)
symbol(s) = Symbol(s)

"""
    strip_namespace(s)

Omit the part before the first dot.
"""
function strip_namespace(s::String); match(r"[^.₊]+$",s).match; end
function strip_namespace(s::Symbol); Symbol(strip_namespace(string(s))); end

"""
    symbols_state(sys::ODESystem)
    symbols_par(sys::ODESystem)

Extract the basic symbols without namespace of system states and system parameters.
"""
function symbols_state(sys::ODESystem); symbol.(states(sys)); end
symbols_par(sys::ODESystem) = symbol.(parameters(sys))

# "apply fun to x until fun(x) == x"
# function fixpoint(fun, x, nrecur_max=12; fmap=identity)
#     nrecur_max == 0 && error("cound not find fixpoint for $f and $x.")
#     px = fun(x)
#     #@show nrecur_max, x, fmap(x), fmap(px)
#     fmap(px) == fmap(x) && return(px)
#     fixpoint(fun, px, nrecur_max-1)
# end


text:
include("gemm.jl")

using UnicodePlots
function plot_data(ns, data)
  scaling = 1e-9 .* (2 * ns.^3)
  init = true

  minv = floatmax()
  maxv = floatmin()
  for (s, v) = pairs(data)
    if s === :ref
      continue
    end
    minv = min(minv, minimum(scaling ./ v.min))
    maxv = max(maxv, maximum(scaling ./ v.min))
  end

  plt = lineplot(ns, scaling ./ data.ref.min, xlim = (ns[1], ns[end]),
                 ylim = (minv, maxv), title="max GFLOPS", xlabel="n",
                 ylabel="time", name = "ref")

  for (s, v) = pairs(data)
    if s == :ref
      continue
    end
    lineplot!(plt, ns, scaling ./ v.min, name = string(s))
  end
  plt
end

using BenchmarkTools

function run(ns, nsamples = 10)

  ref = (min=zeros(length(ns)), avg = zeros(length(ns)))
  ijp = (min=zeros(length(ns)), avg = zeros(length(ns)))
  ipj = (min=zeros(length(ns)), avg = zeros(length(ns)))
  pij = (min=zeros(length(ns)), avg = zeros(length(ns)))
  pji = (min=zeros(length(ns)), avg = zeros(length(ns)))
  jpi = (min=zeros(length(ns)), avg = zeros(length(ns)))
  jip = (min=zeros(length(ns)), avg = zeros(length(ns)))

  for (i, n) = enumerate(ns)
    @show i,n
    m = k = n
    C = rand(m, n)
    A = rand(m, k)
    B = rand(k, n)
    Cref = C + A * B

    bref = @benchmark $C += $A * $B samples=nsamples
    ref.min[i] = minimum(bref.times)
    ref.avg[i] = sum(bref.times) / length(bref.times)

    bijp = @benchmark mygemm_ijp($C, $A, $B) samples=nsamples
    ijp.min[i] = minimum(bijp.times)
    ijp.avg[i] = sum(bijp.times) / length(bijp.times)

    bipj = @benchmark mygemm_ipj($C, $A, $B) samples=nsamples
    ipj.min[i] = minimum(bipj.times)
    ipj.avg[i] = sum(bipj.times) / length(bipj.times)

    bpij = @benchmark mygemm_pij($C, $A, $B) samples=nsamples
    pij.min[i] = minimum(bpij.times)
    pij.avg[i] = sum(bpij.times) / length(bpij.times)

    bpji = @benchmark mygemm_pji($C, $A, $B) samples=nsamples
    pji.min[i] = minimum(bpji.times)
    pji.avg[i] = sum(bpji.times) / length(bpji.times)

    bjpi = @benchmark mygemm_jpi($C, $A, $B) samples=nsamples
    jpi.min[i] = minimum(bjpi.times)
    jpi.avg[i] = sum(bjpi.times) / length(bjpi.times)

    bjip = @benchmark mygemm_jip($C, $A, $B) samples=nsamples
    jip.min[i] = minimum(bjip.times)
    jip.avg[i] = sum(bjip.times) / length(bjip.times)
  end

  data = (ref = ref, ijp = ijp, ipj = ipj, pij = pij, pji = pji, jpi = jpi,
          jip = jip)

  return data
end

ns = 50:50:500
# data = run(ns)
plt = plot_data(ns, data)

max_stars_repo_name: jkozdon/Julia-LAFF-On-PfHP

max_issues_repo_name: jkozdon/Julia-LAFF-On-PfHP

max_forks_repo_name: jkozdon/Julia-LAFF-On-PfHP


text:
using DynamicFactorModeling
using Test

@testset "DynamicFactorModeling.jl" begin
    @test sayhi() == "Hi!"
end


text:
defaultOutputFormat(plt::PlottingObject) = "png"

function png(plt::PlottingObject, fn::@compat(AbstractString))
  fn = addExtension(fn, "png")
  io = open(fn, "w")
  writemime(io, MIME("image/png"), plt)
  close(io)
end
png(fn::@compat(AbstractString)) = png(current(), fn)

function svg(plt::PlottingObject, fn::@compat(AbstractString))
  fn = addExtension(fn, "svg")
  io = open(fn, "w")
  writemime(io, MIME("image/svg+xml"), plt)
  close(io)
end
svg(fn::@compat(AbstractString)) = svg(current(), fn)

function pdf(plt::PlottingObject, fn::@compat(AbstractString))
  fn = addExtension(fn, "pdf")
  io = open(fn, "w")
  writemime(io, MIME("application/pdf"), plt)
  close(io)
end
pdf(fn::@compat(AbstractString)) = pdf(current(), fn)

function ps(plt::PlottingObject, fn::@compat(AbstractString))
  fn = addExtension(fn, "ps")
  io = open(fn, "w")
  writemime(io, MIME("application/postscript"), plt)
  close(io)
end
ps(fn::@compat(AbstractString)) = ps(current(), fn)

function tex(plt::PlottingObject, fn::@compat(AbstractString))
  fn = addExtension(fn, "tex")
  io = open(fn, "w")
  writemime(io, MIME("application/x-tex"), plt)
  close(io)
end
tex(fn::@compat(AbstractString)) = tex(current(), fn)

# ----------------------------------------------------------------

@compat const _savemap = Dict(
    "png" => png,
    "svg" => svg,
    "pdf" => pdf,
    "ps"  => ps,
    "tex" => tex,
  )

function getExtension(fn::@compat(AbstractString))
  pieces = split(fn, ".")
  length(pieces) > 1 || error("Can't extract file extension: ", fn)
  ext = pieces[end]
  haskey(_savemap, ext) || error("Invalid file extension: ", fn)
  ext
end

function addExtension(fn::@compat(AbstractString), ext::@compat(AbstractString))
  try
    oldext = getExtension(fn)
    if oldext == ext
      return fn
    else
      return "$fn.$ext"
    end
  catch
    return "$fn.$ext"
  end
end

function savefig(plt::PlottingObject, fn::@compat(AbstractString))

  # get the extension
  local ext
  try
    ext = getExtension(fn)
  catch
    # if we couldn't extract the extension, add the default
    ext = defaultOutputFormat(plt)
    fn = addExtension(fn, ext)
  end

  # save it
  func = get(_savemap, ext) do
    error("Unsupported extension $ext with filename ", fn)
  end
  func(plt, fn)
end
savefig(fn::@compat(AbstractString)) = savefig(current(), fn)

# savepng(args...; kw...) = savepng(current(), args...; kw...)
# savepng(plt::PlottingObject, fn::@compat(AbstractString); kw...) = (io = open(fn, "w"); writemime(io, MIME("image/png"), plt); close(io))

# ---------------------------------------------------------

gui(plt::PlottingObject = current()) = display(PlotsDisplay(), plt)

# override the REPL display to open a gui window
Base.display(::Base.REPL.REPLDisplay, ::MIME"text/plain", plt::PlottingObject) = gui(plt)


text:
struct ForwardLocal <: MainPathAlgorithm end

struct ForwardLocalResult{T<:Integer, U<:Real} <: MainPathResult
  edges::Vector{Edge{T}}
  weights::Vector{U}
end

function forward_local(
  g::AbstractGraph{T},
  s::Vector{R},
  weights::AbstractMatrix{U}
  ) where {T,R <: Integer} where U <: Real

  sink = ifelse.(outdegree(g) .== 0, true, false)
  visited = falses(nv(g))
  S = copy(s)
  mp_e = Vector{Edge{T}}()

  while !all(sink[S])
    next = Vector{T}()
    for u in S
      (sink[u] || visited[u]) && continue
      maxnb = max_outneighbors(g, u, weights)
      append!(next, maxnb)
      [push!(mp_e, Edge(u, nb)) for nb in maxnb]
      visited[u] = true
    end
    S = next
  end

  from = src.(mp_e)
  to = dst.(mp_e)
  w_e = [weights[s,t] for (s, t) in zip(from, to)]
  return ForwardLocalResult(mp_e, w_e)
end

function main_path(
  g::AbstractGraph{T},
  s::AbstractVector{S},
  weights::AbstractMatrix{U},
  ::ForwardLocal
  ) where {S, T <: Integer} where U <: Real

  return forward_local(g, s, weights)
end

function main_path(
  g::AbstractGraph{T},
  s::AbstractVector{S},
  ::ForwardLocal
  ) where {S, T <: Integer}

  return forward_local(g, s, weights(g))
end

function main_path(
  g::AbstractGraph{T},
  s::U,
  ::ForwardLocal
  ) where {T, U <: Integer}

   return forward_local(g, [s], weights(g))
 end

function main_path(
  g::AbstractGraph{T},
  ::ForwardLocal
  ) where T <: Integer

  w = weights(g)
  sources = get_sources(g)
  sweights = [max_outweight(g, v, w) for v in sources]
  s = sources[sweights .== maximum(sweights)]

  return forward_local(g, s, w)
end


text:
##############################################################################
##
## PooledDataArray type definition
##
## An AbstractDataArray with efficient storage when values are repeated. A
## PDA wraps an array of unsigned integers, which are used to index into a
## compressed pool of values. NA's are 0's in the refs array.  The unsigned
## integer type used for the refs array defaults to UInt32.  The `compact`
## function converts to a smallest integer size that will index the entire pool.
##
## TODO: Allow ordering of factor levels
## TODO: Add metadata for dummy conversion
##
##############################################################################

# This is used as a wrapper during PooledDataArray construction only, to distinguish
# arrays of pool indices from normal arrays
type RefArray{R<:Integer,N}
    a::Array{R,N}
end

type PooledDataArray{T, R<:Integer, N} <: AbstractDataArray{T, N}
    refs::Array{R, N}
    pool::Vector{T}

    function PooledDataArray(rs::RefArray{R, N},
                             p::Vector{T})
        # refs mustn't overflow pool
        if length(rs.a) > 0 && maximum(rs.a) > prod(size(p))
            throw(ArgumentError("Reference array points beyond the end of the pool"))
        end
        new(rs.a,p)
    end
end
typealias PooledDataVector{T,R} PooledDataArray{T,R,1}
typealias PooledDataMatrix{T,R} PooledDataArray{T,R,2}

##############################################################################
##
## PooledDataArray constructors
##
# How do you construct a PooledDataArray from an Array?
# From the same sigs as a DataArray!
#
# Algorithm:
# * Start with:
#   * A null pool
#   * A pre-allocated refs
#   * A hash from T to Int
# * Iterate over d
#   * If value of d in pool already, set the refs accordingly
#   * If value is new, add it to the pool, then set refs
##############################################################################

# Echo inner constructor as an outer constructor
function PooledDataArray{T,R<:Integer,N}(refs::RefArray{R, N},
                                         pool::Vector{T})
    PooledDataArray{T,R,N}(refs, pool)
end

# A no-op constructor
PooledDataArray(d::PooledDataArray) = d

# Constructor from array, w/ pool, missingness, and ref type
function PooledDataArray{T,R<:Integer,N}(d::AbstractArray{T, N},
                                         pool::Vector{T},
                                         m::AbstractArray{Bool, N},
                                         r::Type{R} = DEFAULT_POOLED_REF_TYPE)
    if length(pool) > typemax(R)
        throw(ArgumentError("Cannot construct a PooledDataVector with type $R with a pool of size $(length(pool))"))
    end

    newrefs = Array(R, size(d))
    poolref = Dict{T, R}()

    # loop through once to fill the poolref dict
    for i = 1:length(pool)
        poolref[pool[i]] = i
    end

    # fill in newrefs
    for i = 1:length(d)
        if m[i]
            newrefs[i] = 0
        else
            newrefs[i] = get(poolref, d[i], 0)
        end
    end
    return PooledDataArray(RefArray(newrefs), pool)
end

# Constructor from array, w/ missingness and ref type
function PooledDataArray{T,R<:Integer,N}(d::AbstractArray{T, N},
                                         m::AbstractArray{Bool, N},
                                         r::Type{R} = DEFAULT_POOLED_REF_TYPE)
    pool = convert(Array, unique(d[!m]))
    if method_exists(isless, (T, T))
        sort!(pool)
    end
    PooledDataArray(d, pool, m, r)
end

# Construct an all-NA PooledDataVector of a specific type
PooledDataArray(t::Type, dims::@compat Tuple{Vararg{Int}}) = PooledDataArray(Array(t, dims), trues(dims))
PooledDataArray(t::Type, dims::Int...) = PooledDataArray(Array(t, dims), trues(dims))
PooledDataArray{R<:Integer}(t::Type, r::Type{R}, dims::@compat Tuple{Vararg{Int}}) = PooledDataArray(Array(t, dims), trues(dims), r)
PooledDataArray{R<:Integer}(t::Type, r::Type{R}, dims::Int...) = PooledDataArray(Array(t, dims), trues(dims), r)

# Construct an empty PooledDataVector of a specific type
PooledDataArray(t::Type) = PooledDataArray(similar(Array(t,1),0), trues(0))
PooledDataArray{R<:Integer}(t::Type, r::Type{R}) = PooledDataArray(similar(Array(t,1),0), trues(0), r)

# Convert a BitArray to an Array{Bool} (m = missingness)
# For some reason an additional method is needed but even that doesn't work
# For a BitArray a refs type of UInt8 will always be sufficient as the size of the pool is 0, 1 or 2
PooledDataArray{N}(d::BitArray{N}) = PooledDataArray(Array(d), falses(size(d)), UInt8)
PooledDataArray{N}(d::BitArray{N}, m::AbstractArray{Bool, N}) = PooledDataArray(Array(d), m, UInt8)

# Convert a DataArray to a PooledDataArray
PooledDataArray{T,R<:Integer}(da::DataArray{T},
                              r::Type{R} = DEFAULT_POOLED_REF_TYPE) = PooledDataArray(da.data, da.na, r)
PooledDataArray{T,R<:Integer}(da::DataArray{T},
                              pool::Vector{T},
                              r::Type{R} = DEFAULT_POOLED_REF_TYPE) = PooledDataArray(da.data, pool, da.na, r)

# Convert a Array{T} to a PooledDataArray
PooledDataArray{T,R<:Integer}(d::Array{T},
                              r::Type{R} = DEFAULT_POOLED_REF_TYPE) = PooledDataArray(d, falses(size(d)), r)
PooledDataArray{T,R<:Integer}(d::Array{T},
                              pool::Vector{T},
                              r::Type{R} = DEFAULT_POOLED_REF_TYPE) = PooledDataArray(d, pool, falses(size(d)), r)

# Explicitly convert Ranges into a PooledDataVector
PooledDataArray{R<:Integer}(rs::Range,
                            r::Type{R} = DEFAULT_POOLED_REF_TYPE) = PooledDataArray(collect(rs), falses(length(rs)), r)

# Initialized constructors with 0's, 1's
for (f, basef) in ((:pdatazeros, :zeros), (:pdataones, :ones))
    @eval begin
        ($f)(dims::Int...) = PooledDataArray(($basef)(dims...), falses(dims...))
        ($f)(t::Type, dims::Int...) = PooledDataArray(($basef)(t, dims...), falses(dims...))
        ($f){R<:Integer}(t::Type, r::Type{R}, dims::Int...) = PooledDataArray(($basef)(t, dims...), falses(dims...), r)
    end
end

##############################################################################
##
## Basic size properties of all Data* objects
##
##############################################################################

Base.size(pda::PooledDataArray) = size(pda.refs)
Base.length(pda::PooledDataArray) = length(pda.refs)
Base.endof(pda::PooledDataArray) = endof(pda.refs)

##############################################################################
##
## Copying Data* objects
##
##############################################################################

Base.copy(pda::PooledDataArray) = PooledDataArray(RefArray(copy(pda.refs)),
                                                  copy(pda.pool))
# TODO: Implement copy_to()

function Base.resize!{T,R}(pda::PooledDataArray{T,R,1}, n::Int)
    oldn = length(pda.refs)
    resize!(pda.refs, n)
    pda.refs[oldn+1:n] = zero(R)
    pda
end

##############################################################################
##
## Predicates, including the new isna()
##
##############################################################################

function Base.isnan(pda::PooledDataArray)
    PooledDataArray(RefArray(copy(pda.refs)), isnan(pda.pool))
end

function Base.isfinite(pda::PooledDataArray)
    PooledDataArray(RefArray(copy(pda.refs)), isfinite(pda.pool))
end

isna(pda::PooledDataArray) = pda.refs .== 0

#' @description
#'
#' Safe and type-stable way to determine if element `i` of an
#' PooledDataArray is `NA`.
#'
#' @param pda::PooledDataArray The PooledDataArray whose missingness will
#'        be assessed.
#' @param i::Integer The index of the element to be checked for `NA`.
#'
#' @returns na::Bool Is the element `NA` or not?
#'
#' @examples
#'
#' a = @pdata([1, 2, 3])
#' isna(a, 1)
isna(pda::PooledDataArray, i::Real) = pda.refs[i] == 0 # -> Bool

#' @description
#'
#' Determine if any of the entries of an PooledDataArray are `NA`.
#'
#' @param pda::PooledDataArray The PooledDataArray whose elements will
#'        be assessed.
#'
#' @returns out::Bool Are any of the elements of `pda` an `NA` value?
#'
#' @examples
#'
#' pda = @pdata([1, 2, 3])
#' anyna(pda)
function anyna(pda::PooledDataArray)
    for ref in pda.refs
        ref == 0 && return true
    end
    return false
end

#' @description
#'
#' Determine if all of the entries of an PooledDataArray are `NA`.
#'
#' @param pda::PooledDataArray The PooledDataArray whose elements will
#'        be assessed.
#'
#' @returns out::Bool Are all of the elements of `pda` an `NA` value?
#'
#' @examples
#'
#' pda = @pdata([1, 2, 3])
#' allna(pda)
function allna(pda::PooledDataArray)
    for ref in pda.refs
        ref == 0 || return false
    end
    return true
end

##############################################################################
##
## PooledDataArray utilities
##
## TODO: Add methods with these names for DataArray's
##
##############################################################################

function compact{T,R<:Integer,N}(d::PooledDataArray{T,R,N})
    sz = length(d.pool)

    REFTYPE = sz <= typemax(UInt8)  ? UInt8 :
              sz <= typemax(UInt16) ? UInt16 :
              sz <= typemax(UInt32) ? UInt32 :
                                      UInt64

    if REFTYPE == R
        return d
    end

    newrefs = convert(Array{REFTYPE, N}, d.refs)
    PooledDataArray(RefArray(newrefs), d.pool)
end

#' @description
#'
#' Return a DataVector containing the unique values of a `PooledDataArray`,
#' in the order they appear in the data, including `NA` if any missing entries
#' are encountered. For `PooledDataArray`s, this function is much less efficient
#' than `levels`, which does not return the values in the same order.
#'
#' @param da::DataArray{T} `DataArray` whose unique values are desired.
#'
#' @returns dv::DataVector{T} `DataVector` containing the unique values
#'          from `pda`, in the order they appear, including `NA` if there are
#'          any missing entries in `pda`.
#'
#' @examples
#'
#' pdv = @pdata [1, -2, 1, NA, 4]
#' distinct_values = unique(pdv)
function Base.unique{T}(pda::PooledDataArray{T})
    n = length(pda)
    nlevels = length(pda.pool)
    unique_values = Array(T, 0)
    sizehint!(unique_values, nlevels)
    seen = Set{eltype(pda.refs)}()

    firstna = 0
    for i in 1:n
        if isna(pda, i)
            if firstna == 0
                firstna = length(unique_values) + 1
            end
        elseif !in(pda.refs[i], seen)
            push!(seen, pda.refs[i])
            push!(unique_values, pda.pool[pda.refs[i]])
        else
            continue
        end

        if firstna > 0 && length(unique_values) == nlevels
            break
        end
    end

    if firstna > 0
        res = DataArray(Array(T, nlevels + 1))
        i = 0
        for val in unique_values
            i += 1
            if i == firstna
                res.na[i] = true
                i += 1
            end
            res.data[i] = val
        end

        if firstna == nlevels + 1
            res.na[nlevels + 1] = true
        end

        return res
    else
        return DataArray(unique_values)
    end
end

#' @description
#'
#' Return a DataVector containing the levels of a `PooledDataArray`,
#' excluding `NA`. For `PooledDataArray`s, this function is much more
#' efficient than `unique`, and it returns the levels in the order used
#' when calling `setlevels!` rather than in the order of appearance.
#' Contrary to `unique`, it also returns levels which are not present in
#' the vector.
#'
#'
#' @param pda::PooledDataArray{T} PooledDataArray whose levels values are
#' desired.
#'
#' @returns dv::DataVector{T} DataVector containing the levels
#'          of `da`, excluding `NA`.
#'
#' @examples
#'
#' pdv = @pdata [1, -2, 1, NA, 4]
#' distinct_values = levels(pdv)
levels{T}(pda::PooledDataArray{T}) = copy(pda.pool)

function PooledDataArray{S,R,N}(x::PooledDataArray{S,R,N},
                                newpool::Vector{S})
    # QUESTION: should we have a ! version of this? If so, needs renaming?
    tidx::Array{R} = findat(newpool, x.pool)
    refs = zeros(R, length(x))
    for i in 1:length(refs)
        if x.refs[i] != 0
            refs[i] = tidx[x.refs[i]]
        end
    end
    PooledDataArray(RefArray(refs), newpool)
end

myunique(x::AbstractVector) = x[sort(unique(findat(x, x)))]  # gets the ordering right
myunique(x::AbstractDataVector) = myunique(dropna(x))   # gets the ordering right; removes NAs

function setlevels{T,R}(x::PooledDataArray{T,R}, newpool::AbstractVector)
    pool = myunique(newpool)
    refs = zeros(R, length(x))
    tidx::Array{R} = findat(pool, newpool)
    tidx[isna(newpool)] = 0
    for i in 1:length(refs)
        if x.refs[i] != 0
            refs[i] = tidx[x.refs[i]]
        end
    end
    return PooledDataArray(RefArray(refs), pool)
end

function setlevels!{T,R}(x::PooledDataArray{T,R}, newpool::AbstractVector{T})
    if newpool == myunique(newpool) # no NAs or duplicates
        x.pool = newpool
        return x
    else
        x.pool = myunique(newpool)
        tidx::Array{R} = findat(x.pool, newpool)
        tidx[isna(newpool)] = 0
        for i in 1:length(x.refs)
            if x.refs[i] != 0
                x.refs[i] = tidx[x.refs[i]]
            end
        end
        return x
    end
end

setlevels!{T, R}(x::PooledDataArray{T, R},
                 newpool::AbstractVector) = setlevels!(x, convert(Array{T}, newpool))

function setlevels(x::PooledDataArray, d::Dict)
    newpool = copy(DataArray(x.pool))
    # An NA in `v` is put in the pool; that will cause it to become NA
    for (k,v) in d
        idx = findin(newpool, [k])
        if length(idx) == 1
            newpool[idx[1]] = v
        end
    end
    setlevels(x, newpool)
end

function setlevels!{T,R}(x::PooledDataArray{T,R}, d::Dict{T,T})
    for (k,v) in d
        idx = findin(x.pool, [k])
        if length(idx) == 1
            x.pool[idx[1]] = v
        end
    end
    x
end

function setlevels!{T,R}(x::PooledDataArray{T,R}, d::Dict{T,Any}) # this version handles NAs in d's values
    newpool = copy(DataArray(x.pool))
    # An NA in `v` is put in the pool; that will cause it to become NA
    for (k,v) in d
        idx = findin(newpool, [k])
        if length(idx) == 1
            newpool[idx[1]] = v
        end
    end
    setlevels!(x, newpool)
end

reorder(x::PooledDataArray) = PooledDataArray(x, sort(levels(x)))  # just re-sort the pool

reorder(x::PooledDataArray, y::AbstractVector...) = reorder(mean, x, y...)

### FIXME: this can't work because we don't know about DataFrames
# reorder(fun::Function, x::PooledDataArray, y::AbstractVector...) =
#     reorder(fun, x, DataFrame({y...}))

Base.reverse(x::PooledDataArray) = PooledDataArray(RefArray(reverse(x.refs)), x.pool)

function Base.permute!!{T<:Integer}(x::PooledDataArray, p::AbstractVector{T})
    Base.permute!!(x.refs, p)
    x
end

function Base.ipermute!!{T<:Integer}(x::PooledDataArray, p::AbstractVector{T})
    Base.ipermute!!(x.refs, p)
    x
end

##############################################################################
##
## similar()
##
##############################################################################

function Base.similar{T,R}(pda::PooledDataArray{T,R}, S::Type, dims::Dims)
    PooledDataArray(RefArray(zeros(R, dims)), S[])
end

##############################################################################
##
## find()
##
##############################################################################

Base.find(pdv::PooledDataVector{Bool}) = find(convert(Vector{Bool}, pdv, false))

##############################################################################
##
## setindex!() definitions
##
##############################################################################

function getpoolidx{T,R}(pda::PooledDataArray{T,R}, val::Any)
    val::T = convert(T,val)
    pool_idx = findfirst(pda.pool, val)
    if pool_idx <= 0
        push!(pda.pool, val)
        pool_idx = length(pda.pool)
        if pool_idx > typemax(R)
            throw(ArgumentError(
                "You're using a PooledDataArray with ref type $R,\n" *
                "which can only hold $(Int(typemax(R))) values, and you just tried to add the $(Int(typemax(R))+1)th\n" *
                "reference. Please change the ref type to a larger int type, or use the\n" *
                "default ref type ($DEFAULT_POOLED_REF_TYPE)."
            ))
        end
    end
    return pool_idx
end

getpoolidx{T,R}(pda::PooledDataArray{T,R}, val::NAtype) = zero(R)

##############################################################################
##
## show() and similar methods
##
##############################################################################

function Base.string(x::PooledDataVector)
    tmp = join(x, ", ")
    return "[$tmp]"
end

# # Need setindex!()'s to make this work
# This is broken now because the inner show returns to the outer show.
# function show(io::IO, pda::PooledDataArray)
#     invoke(show, Tuple{Any,AbstractArray}, io, pda)
#     print(io, "\nlevels: ")
#     print(io, levels(pda))
# end

##############################################################################
##
## Replacement operations
##
##############################################################################

function replace!(x::PooledDataArray{NAtype}, fromval::NAtype, toval::NAtype)
    NA # no-op to deal with warning
end
function replace!(x::PooledDataArray, fromval::NAtype, toval::NAtype)
    NA # no-op to deal with warning
end
function replace!{S, T}(x::PooledDataArray{S}, fromval::T, toval::NAtype)
    fromidx = findfirst(x.pool, fromval)
    if fromidx == 0
        throw(ErrorException("can't replace a value not in the pool in a PooledDataVector!"))
    end

    x.refs[x.refs .== fromidx] = 0

    return NA
end
function replace!{S, T}(x::PooledDataArray{S}, fromval::NAtype, toval::T)
    toidx = findfirst(x.pool, toval)
    # if toval is in the pool, just do the assignment
    if toidx != 0
        x.refs[x.refs .== 0] = toidx
    else
        # otherwise, toval is new, add it to the pool
        push!(x.pool, toval)
        x.refs[x.refs .== 0] = length(x.pool)
    end

    return toval
end
function replace!{R, S, T}(x::PooledDataArray{R}, fromval::S, toval::T)
    # throw error if fromval isn't in the pool
    fromidx = findfirst(x.pool, fromval)
    if fromidx == 0
        throw(ErrorException("can't replace a value not in the pool in a PooledDataArray!"))
    end

    # if toval is in the pool too, use that and remove fromval from the pool
    toidx = findfirst(x.pool, toval)
    if toidx != 0
        x.refs[x.refs .== fromidx] = toidx
        #x.pool[fromidx] = None    TODO: what to do here??
    else
        # otherwise, toval is new, swap it in
        x.pool[fromidx] = toval
    end

    return toval
end

##############################################################################
##
## Sorting can use the pool to speed things up
##
##############################################################################

function Base.sortperm(pda::PooledDataArray; alg::Base.Sort.Algorithm=Base.Sort.DEFAULT_UNSTABLE,
                       lt::Function=isless, by::Function=identity,
                       rev::Bool=false, order=Base.Sort.Forward)
    order = Base.ord(lt, by, rev, order)

    # TODO handle custom ordering efficiently
    if !isa(order, Base.Order.ForwardOrdering) && !isa(order, Base.Order.ReverseOrdering)
        return sort!(collect(1:length(pda)), alg, Base.Order.Perm(order,pda))
    end

    # TODO handle non-sorted keys without copying
    perm = issorted(pda.pool) ? groupsort_indexer(pda, true)[1] : sortperm(reorder(pda))
    isa(order, Base.Order.ReverseOrdering) && reverse!(perm)
    perm
end

Base.sort(pda::PooledDataArray; kw...) = pda[sortperm(pda; kw...)]

type FastPerm{O<:Base.Sort.Ordering,V<:AbstractVector} <: Base.Sort.Ordering
    ord::O
    vec::V
end
Base.sortperm{V}(x::AbstractVector, a::Base.Sort.Algorithm, o::FastPerm{Base.Sort.ForwardOrdering,V}) = x[sortperm(o.vec)]
Base.sortperm{V}(x::AbstractVector, a::Base.Sort.Algorithm, o::FastPerm{Base.Sort.ReverseOrdering,V}) = x[reverse(sortperm(o.vec))]
Perm{O<:Base.Sort.Ordering}(o::O, v::PooledDataVector) = FastPerm(o, v)

##############################################################################
##
## PooledDataVecs: EXPLANATION SHOULD GO HERE
##
##############################################################################

function PooledDataVecs{S,Q<:Integer,R<:Integer,N}(v1::PooledDataArray{S,Q,N},
                                                   v2::PooledDataArray{S,R,N})
    pool = sort(unique([v1.pool; v2.pool]))
    sz = length(pool)

    REFTYPE = sz <= typemax(UInt8)  ? UInt8 :
              sz <= typemax(UInt16) ? UInt16 :
              sz <= typemax(UInt32) ? UInt32 :
                                      UInt64

    tidx1 = convert(Vector{REFTYPE}, findat(pool, v1.pool))
    tidx2 = convert(Vector{REFTYPE}, findat(pool, v2.pool))
    refs1 = zeros(REFTYPE, length(v1))
    refs2 = zeros(REFTYPE, length(v2))
    for i in 1:length(refs1)
        if v1.refs[i] != 0
            refs1[i] = tidx1[v1.refs[i]]
        end
    end
    for i in 1:length(refs2)
        if v2.refs[i] != 0
            refs2[i] = tidx2[v2.refs[i]]
        end
    end
    return (PooledDataArray(RefArray(refs1), pool),
            PooledDataArray(RefArray(refs2), pool))
end

function PooledDataVecs{S,R<:Integer,N}(v1::PooledDataArray{S,R,N},
                                        v2::AbstractArray{S,N})
    return PooledDataVecs(v1,
                          PooledDataArray(v2))
end

####
function PooledDataVecs{S,R<:Integer,N}(v1::AbstractArray{S,N},
                                        v2::PooledDataArray{S,R,N})
    return PooledDataVecs(PooledDataArray(v1),
                          v2)
end

function PooledDataVecs(v1::AbstractArray,
                        v2::AbstractArray)

    ## Return two PooledDataVecs that share the same pool.

    ## TODO: allow specification of REFTYPE
    refs1 = Array(DEFAULT_POOLED_REF_TYPE, size(v1))
    refs2 = Array(DEFAULT_POOLED_REF_TYPE, size(v2))
    poolref = Dict{promote_type(eltype(v1), eltype(v2)), DEFAULT_POOLED_REF_TYPE}()
    maxref = 0

    # loop through once to fill the poolref dict
    for i = 1:length(v1)
        if !isna(v1[i])
            poolref[v1[i]] = 0
        end
    end
    for i = 1:length(v2)
        if !isna(v2[i])
            poolref[v2[i]] = 0
        end
    end

    # fill positions in poolref
    pool = sort(collect(keys(poolref)))
    i = 1
    for p in pool
        poolref[p] = i
        i += 1
    end

    # fill in newrefs
    zeroval = zero(DEFAULT_POOLED_REF_TYPE)
    for i = 1:length(v1)
        if isna(v1[i])
            refs1[i] = zeroval
        else
            refs1[i] = poolref[v1[i]]
        end
    end
    for i = 1:length(v2)
        if isna(v2[i])
            refs2[i] = zeroval
        else
            refs2[i] = poolref[v2[i]]
        end
    end

    return (PooledDataArray(RefArray(refs1), pool),
            PooledDataArray(RefArray(refs2), pool))
end

Base.convert{S,T,R1<:Integer,R2<:Integer,N}(::Type{PooledDataArray{S,R1,N}}, pda::PooledDataArray{T,R2,N}) =
    PooledDataArray(RefArray(convert(Array{R1,N}, pda.refs)), convert(Vector{S}, pda.pool))
Base.convert{S,T,R<:Integer,N}(::Type{PooledDataArray{S,R,N}}, pda::PooledDataArray{T,R,N}) =
    PooledDataArray(RefArray(copy(pda.refs)), convert(Vector{S}, pda.pool))
Base.convert{T,R<:Integer,N}(::Type{PooledDataArray{T,R,N}}, pda::PooledDataArray{T,R,N}) = pda
Base.convert{S,T,R1<:Integer,R2<:Integer,N}(::Type{PooledDataArray{S,R1}}, pda::PooledDataArray{T,R2,N}) =
    convert(PooledDataArray{S,R1,N}, pda)
Base.convert{S,T,R<:Integer,N}(::Type{PooledDataArray{S}}, pda::PooledDataArray{T,R,N}) =
    convert(PooledDataArray{S,R,N}, pda)
Base.convert{T,R<:Integer,N}(::Type{PooledDataArray}, pda::PooledDataArray{T,R,N}) = pda

Base.convert{T,R<:Integer,N}(::Type{PooledDataArray{T,R,N}}, a::DataArray{T,N}) =
    PooledDataArray(a.data, a.na, R)
Base.convert{S,T,R<:Integer,N}(::Type{PooledDataArray{S,R,N}}, a::DataArray{T,N}) =
    convert(PooledDataArray{S,R,N}, PooledDataArray(a.data, a.na, R))
Base.convert{S,T,R<:Integer,N}(::Type{PooledDataArray{S,R}}, a::DataArray{T,N}) =
    convert(PooledDataArray{S,R,N}, PooledDataArray(a.data, a.na, R))
Base.convert{S,T,N}(::Type{PooledDataArray{S}}, a::DataArray{T,N}) =
    convert(PooledDataArray{S}, PooledDataArray(a.data, a.na))
Base.convert(::Type{PooledDataArray}, a::DataArray) =
    PooledDataArray(a.data, a.na)

Base.convert{S,T,R<:Integer,N}(::Type{PooledDataArray{S,R,N}}, a::AbstractArray{T,N}) =
    PooledDataArray(convert(Array{S,N}, a), falses(size(a)), R)
Base.convert{S,T,R<:Integer,N}(::Type{PooledDataArray{S,R}}, a::AbstractArray{T,N}) =
    PooledDataArray(convert(Array{S,N}, a), falses(size(a)), R)
Base.convert{S,T,N}(::Type{PooledDataArray{S}}, a::AbstractArray{T,N}) =
    PooledDataArray(convert(Array{S,N}, a), falses(size(a)))
Base.convert(::Type{PooledDataArray}, a::AbstractArray) =
    PooledDataArray(a, falses(size(a)))

function Base.convert{S,T,R<:Integer,N}(::Type{DataArray{S,N}},
                                        pda::PooledDataArray{T,R,N})
    res = DataArray(Array(S, size(pda)), BitArray(size(pda)))
    for i in 1:length(pda)
        r = pda.refs[i]
        if r == 0 # TODO: Use zero(R)
            res.na[i] = true
        else
            res.na[i] = false
            res.data[i] = pda.pool[r]
        end
    end
    return res
end
Base.convert{S,T,R<:Integer,N}(::Type{DataArray{S}}, pda::PooledDataArray{T,R,N}) =
    convert(DataArray{S,N}, pda)
Base.convert{T,R<:Integer,N}(::Type{DataArray}, pda::PooledDataArray{T,R,N}) =
    convert(DataArray{T,N}, pda)

pdata(a::AbstractArray) = convert(PooledDataArray, a)

function Base.convert{S, T, R, N}(
    ::Type{Array{S, N}},
    pda::PooledDataArray{T, R, N}
)
    res = Array(S, size(pda))
    for i in 1:length(pda)
        if pda.refs[i] == zero(R)
            throw(NAException())
        else
            res[i] = pda.pool[pda.refs[i]]
        end
    end
    return res
end

function Base.convert{T, R}(::Type{Vector}, pdv::PooledDataVector{T, R})
    return convert(Array{T, 1}, pdv)
end

function Base.convert{T, R}(::Type{Matrix}, pdm::PooledDataMatrix{T, R})
    return convert(Array{T, 2}, pdm)
end

function Base.convert{T, R, N}(::Type{Array}, pda::PooledDataArray{T, R, N})
    return convert(Array{T, N}, pda)
end

function Base.convert{S, T, R, N}(
    ::Type{Array{S, N}},
    pda::PooledDataArray{T, R, N},
    replacement::Any
)
    res = Array(S, size(pda))
    replacementS = convert(S, replacement)
    for i in 1:length(pda)
        if pda.refs[i] == zero(R)
            res[i] = replacementS
        else
            res[i] = pda.pool[pda.refs[i]]
        end
    end
    return res
end

function Base.convert{T, R}(::Type{Vector}, pdv::PooledDataVector{T, R}, replacement::Any)
    return convert(Array{T, 1}, pdv, replacement)
end

function Base.convert{T, R}(::Type{Matrix}, pdm::PooledDataMatrix{T, R}, replacement::Any)
    return convert(Array{T, 2}, pdm, replacement)
end

function Base.convert{T, R, N}(::Type{Array}, pda::PooledDataArray{T, R, N}, replacement::Any)
    return convert(Array{T, N}, pda, replacement)
end

function dropna{T}(pdv::PooledDataVector{T})
    n = length(pdv)
    res = Array(T, n)
    total = 0
    for i in 1:n
        if pdv.refs[i] > 0
            total += 1
            res[total] = pdv.pool[pdv.refs[i]]
        end
    end
    resize!(res, total)
    return res
end


text:
"""
Function to find the labels which seem to be the most relevant and should be called next.
It finds the x labels with the highest score in the scorestatistic.
Arguments:
unibeta::Unibeta -> Contains the scores to find the new labels
numberofcovs::Int -> Number of covariances which should be called additionally
usedlabels::Array{ASCIIString, 1} -> Array with the labels which already been used (no need to ask for them again)
"""
function selectionofcovs(unibeta::Unibeta, numberofcovs::Int, usedlabels::Array{ASCIIString,1} = Array{ASCIIString,1}())
	dummy = deepcopy(unibeta)
	dummy.unibeta = dummy.unibeta.^2
	if(!isempty(usedlabels))
		for i = 1 : length(usedlabels)
			select = findfirst(dummy.labels, usedlabels[i])
			deleteat!(dummy.unibeta, select)
			deleteat!(dummy.labels, select)
		end
	end
	wantedlabels = Array{ASCIIString,1}(numberofcovs)
	for i = 1 : length(wantedlabels)
		select = indmax(dummy.unibeta)
		deleteat!(dummy.unibeta, select)
		wantedlabels[i] = dummy.labels[select]
		deleteat!(dummy.labels, select)
	end
	return wantedlabels
end
"""
Function to find the labels which seem to be the most relevant and should be called next.
It finds the most relevant labels with respect to the heuristic approach and calls the function selectionofcovs() for the x
labels which should be called additionally.
Arguments:
myscratch::Boostscratch -> Container with data which are relevant to perform the boostinsteps
newselval::Float64 -> value which is used as a criteria for choosing the next labels
"""
function getselections(myscratch::Boostscratch, newselval::Float64)
	wantedlabels = Array{ASCIIString,1}()
	for i = 1 : length(myscratch.actualnom)
		if(myscratch.actualnom[i]^2 >= newselval && findfirst(myscratch.pooledcovarmat.labels, myscratch.pooledunibeta.labels[i]) == 0)
			wantedlabels = vcat(wantedlabels, myscratch.pooledunibeta.labels[i])
		end
	end
	if(myscratch.x > 0 && length(wantedlabels) > 0)
		wantedlabels = vcat(wantedlabels, selectionofcovs(Unibeta(myscratch.actualnom, ones(1), myscratch.pooledunibeta.labels), myscratch.x, vcat(myscratch.pooledcovarmat.labels, wantedlabels)))
	end
	return wantedlabels
end

function getcovarfromtriangular(cov::Array{Float64,2}, index::Int)
	return reshape(hcat(reshape(cov[1:index,index],1,index),cov[index,(index+1):size(cov)[2]]),size(cov)[2])
end

function merge_covarmats!(myscratch::Boostscratch, expansion::Array{Float64, 2})
	myscratch.pooledcovarmat.covarmat = hcat(vcat(myscratch.pooledcovarmat.covarmat,zeros(size(expansion)[1] - size(myscratch.pooledcovarmat.covarmat)[1], size(myscratch.pooledcovarmat.covarmat)[2])),expansion)
	myscratch.pooledcovarmat.labels = vcat(myscratch.pooledcovarmat.labels, myscratch.wantedlabels)
end
"""
Function to provide the covarmats and univariate estimators of a cohort
Arguments:
path_data::ASCIIString -> path which leads to the individual data
path_data_out::ASCIIString -> path which defines where the calculated data should be stores
path_labels::ASCIIString -> path to the file with the wanted labels
"""
function provide(path_data::ASCIIString, path_data_out::ASCIIString, path_labels::ASCIIString = "")
	mycohort = read_data(path_data)
	mycohort = standcohort(mycohort)
	if(length(path_labels) > 0)
		labels = read_labels(path_labels)
		wantedlabels = labels[1]
		usedlabels = labels[2]
		if(length(usedlabels) > 0)
			covarmat = calc_covarmat(mycohort, wantedlabels, usedlabels)
		else
			covarmat = calc_covarmat(mycohort, wantedlabels)
		end
		write_covarmat(covarmat, path_data_out)
	else
		unibeta = calc_unibeta(mycohort)
		write_unibeta(unibeta, path_data_out)
	end
end

function labelmap(unibeta::Unibeta,labels)
	lindex = zeros(Int,length(labels))
	cdict = Dict{ASCIIString,Int}([unibeta.labels[i] => i for i=1:length(unibeta.labels)])
	[get(cdict,arg,0) for arg in labels]
end

function labelmap(cohort::Discohort,labels)
	lindex = zeros(Int,length(labels))
	cdict = Dict{ASCIIString,Int}([cohort.labels[i] => i for i=1:length(cohort.labels)])
	[get(cdict,arg,0) for arg in labels]
end
"""
Function which performs standardization to achieve comparability across different cohorts
cohort::Discohort -> contains all information of one cohort
"""
function standcohort(cohort::Discohort)
	cohort_n = length(cohort.y)
	p = size(cohort.x,2)
	xmean = mean(cohort.x,1)
	xsd = std(cohort.x,1)
	xstd = Float64[(cohort.x[i,j] - xmean[j])/xsd[j] for i in 1:cohort_n, j in 1:p]
	ystd = deepcopy(cohort.y)
	ystd -= mean(ystd)
	ystd /= std(ystd)

	return Discohort(xstd,ystd,cohort.labels)
end

max_stars_repo_name: danielazoeller/DistributedBoosting

max_issues_repo_name: danielazoeller/DistributedBoosting

max_forks_repo_name: danielazoeller/DistributedBoosting


text:
"""
    iaaft(ts::AbstractArray{T, 1} where T;
            n_maxiter = 200, tol = 1e-6, n_windows = 50)

Generate an iteratively adjusted amplitude adjusted Fourier transform (IAAFT)
([Schreiber & Schmitz, 1996](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.77.635))
surrogate realization of `ts`.

## Arguments

- **`ts`**: the time series for which to generate an AAFT surrogate realization.

- **`n_maxiter`**: sets the maximum number of iterations to allow before ending
    the algorithm (if convergence is slow).

- **`tol`**: the relative tolerance for deciding if convergence is achieved.

- **`n_window`**: the number is windows used when binning the periodogram (used
    for determining convergence).

## References

T. Schreiber; A. Schmitz (1996). "Improved Surrogate Data for Nonlinearity
Tests". Phys. Rev. Lett. 77 (4): 635–638. doi:10.1103/PhysRevLett.77.635.
PMID 10062864. [https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.77.635](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.77.635)

"""
function iaaft(ts::AbstractVector{T} where T;
                n_maxiter = 200, tol = 1e-6, n_windows = 50)
    any(isnan.(ts)) && throw(DomainError(NaN,"The input must not contain NaN values."))

    # Sorted version of the original time series
    original_sorted = sort(ts)

    # Fourier transform and its amplitudes
    original_fft = fft(ts)
    original_fft_amplitudes = abs.(original_fft)

    # RANK ORDERING.
    # Create some Gaussian noise, and find the indices that sorts it with
    # increasing amplitude. Then sort the original time series according to
    # the indices rendering the Gaussian noise sorted.
    n = length(ts)
    g = rand(Normal(), n)
    inds = sortperm(g)
    ts_sorted = ts[inds]

    iter = 1
    success = false

    spectrum = Vector{Complex{Float64}}(undef, n)
    surrogate = Vector{Float64}(undef, n)

    diffs = zeros(Float64, 2)
    while iter <= n_maxiter
        # Take the Fourier transform of `ts_sorted` and get the phase angles of
        # the resulting complex numbers.
        FT = fft(ts_sorted)
        phase_angles = angle.(FT)

        # The new spectrum preserves the amplitudes of the Fourier transform of
        # the original time series, but randomises the phases (because the
        # phases are derived from the *randomly sorted* version of the original
        # time series).
        spectrum[:] = original_fft_amplitudes .* exp.(phase_angles .* 1im)

        # Now, let the surrogate time series be the values of the original time
        # series, but sorted according to the new spectrum. The shuffled series
        # is generated by taking the inverse Fourier transform of the spectrum
        # consisting of the original amplitudes, but having randomised phases.
        surrogate[:] = real(ifft(spectrum)) # ifft normalises by default

        # map original values onto shuffle
        surrogate[sortperm(surrogate)] = original_sorted

        ts_sorted[:] = surrogate

        # Convergence check
        periodogram = DSP.mt_pgram(ts)
        periodogram_surr = DSP.mt_pgram(surrogate)

        power_binned = interp([x for x in periodogram.freq],
                            periodogram.power,
                            n_windows)

        power_binned_surr = interp([x for x in periodogram_surr.freq],
                            periodogram_surr.power,
                            n_windows)

        if iter == 1
            diffs[1] = sum((power_binned[2] .- power_binned_surr[2]).^2) /
                        sum(power_binned[2].^2)
        else
            diffs[2] = sum((power_binned[2] .- power_binned_surr[2]).^2) /
                        sum(power_binned[2].^2)
            if abs(diffs[1] - diffs[2]) < tol
                break
            end
            diffs[1] = copy(diffs[2])
        end

        iter += 1
    end
    surrogate
end

"""
    iaaft_iters(ts::AbstractArray{T, 1} where T;
                n_maxiter = 100, tol = 1e-5, n_windows = 50)

Generate an iteratively adjusted amplitude adjusted Fourier transform (IAAFT) [1]
surrogate series for `ts` and return a vector containing the surrogate series
from each iteration. The last vector contains the final surrogate.

# Literature references
1. T. Schreiber; A. Schmitz (1996). "Improved Surrogate Data for Nonlinearity
Tests". Phys. Rev. Lett. 77 (4): 635–638. doi:10.1103/PhysRevLett.77.635. PMID
10062864.
"""
function iaaft_iters(ts::AbstractArray{T, 1} where T;
                        n_maxiter = 100, tol = 1e-5, n_windows = 50)

    # Sorted version of the original time series
    original_sorted = sort(ts)

    # Fourier transform and its amplitudes
    original_fft = fft(ts)
    original_fft_amplitudes = abs.(original_fft)

    # Create some Gaussian noise, and find the indices that sorts it with
    # increasing amplitude. Then sort the original time series according to the
    # indices rendering the Gaussian noise sorted.
    n = length(ts)
    g = rand(Normal(), n)
    inds = sortperm(g)
    ts_sorted = ts[inds]

    iter = 1
    success = false

    spectrum = Vector{Complex{Float64}}(undef, n)
    surrogate = Vector{Float64}(undef, n)
    surrogates = Vector{Vector{Float64}}(undef, 0)

    diffs = zeros(Float64, 2)
    while iter <= n_maxiter
        # Take the Fourier transform of `ts_sorted` and get the phase angles
        # of the resulting complex numbers.
        ts_sorted_fft = fft(ts_sorted)
        phase_angles = angle.(ts_sorted_fft)

        # The new spectrum preserves the amplitudes of the Fourier transform
        # of the original time series, but randomises the phases (because the
        # phases are derived from the *randomly sorted* version of the original
        # time series).
        spectrum .= original_fft_amplitudes .* exp.(phase_angles .* 1im)

        # Now, let the surrogate time series be the values of the original time
        # series, but sorted according to the new spectrum. The shuffled series
        # is generated by taking the inverse Fourier transform of the spectrum
        # consisting of the original amplitudes, but having randomised phases.
        surrogate[:] = real(ifft(spectrum)) # ifft normalises by default
        surrogate[sortperm(surrogate)] = original_sorted # map original values onto shuffle
        push!(surrogates, surrogate[:])

        ts_sorted .= surrogate

        # Convergence check
        periodogram = DSP.mt_pgram(ts)
        periodogram_surr = DSP.mt_pgram(surrogate)

        power_binned = interp([x for x in periodogram.freq],
                                periodogram.power,
                                n_windows)

        power_binned_surr = interp([x for x in periodogram_surr.freq],
                                    periodogram_surr.power,
                                    n_windows)

        if iter == 1
            diffs[1] = sum((power_binned[2] .- power_binned_surr[2]).^2) /
                        sum(power_binned[2].^2)
        else
            diffs[2] = sum((power_binned[2] .- power_binned_surr[2]).^2) /
                        sum(power_binned[2].^2)

            if abs(diffs[1] - diffs[2]) < tol
                break
            end
            diffs[1] = copy(diffs[2])
        end

        iter += 1
    end
    surrogates
end

export iaaft, iaaft_iters


text:
using DataFrames
using CSV

function save_simulations(D,TP_local,state,file_date,onset_dates,rng)
    """
    Saves a CSV in the same format as required for the ensemble forecast.
    Also saves a CSV for the TP paths used per simulation.
    """
    df_observed = DataFrame()
    df_TP_local = DataFrame()
    df_observed[!,"state"] = [state for _ in 1:length(onset_dates)]
    df_observed[!,"onset date"] = onset_dates
    df_TP_local[!,"state"] = [state for _ in 1:length(onset_dates)]
    df_TP_local[!,"onset date"] = onset_dates
    # total observed cases
    D_observed = D[:,1,:] + D[:,2,:]

    # indexer for naming the columns
    i = 0
    for (d,tp) in zip(eachcol(D_observed),eachcol(TP_local))
        df_observed[!,"sim"*string(i)] = d[(end-length(onset_dates)+1):end]
        df_TP_local[!,"sim"*string(i)] = tp[(end-length(onset_dates)+1):end]
        i += 1
        if i == 2000
            break
        end
    end

    if i == 0
        for i in 0:1999
            # if nothing sampled, set observations to be missing
            df_observed[!,"sim"*string(i)] = [missing for _ in 1:length(onset_dates)]
            df_TP_local[!,"sim"*string(i)] = [missing for _ in 1:length(onset_dates)]
        end
    elseif i < 2000
        # if we are under the required number of sims, sample randomly from the good sims
        for j in i:1999
            # sample a random column
            ind = rand(rng, 1:(i-1))
            d = D_observed[:,ind]
            tp = TP_local[:,ind]
            df_observed[!,"sim"*string(j)] = d[(end-length(onset_dates)+1):end]
            df_TP_local[!,"sim"*string(j)] = tp[(end-length(onset_dates)+1):end]
        end
    end

    # this adds a column for the date we are supplied the NINDSS file
    df_observed[!,"data date"] .= file_date
    df_TP_local[!,"data date"] .= file_date

    # directory path, and check to see whether it's good
    dir_name = joinpath("results", "UoA_forecast_output", file_date)
    if !ispath(dir_name)
        mkpath(dir_name)
    end

    # file name is just the state and file date
    sim_file_name = state*"_"*file_date*"_sim.csv"
    TP_file_name = state*"_"*file_date*"_TP.csv"
    CSV.write(dir_name*"/"*sim_file_name, df_observed)
    CSV.write(dir_name*"/"*TP_file_name, df_TP_local)

    return nothing

end

function merge_simulation_files(file_date)
    """
    Merge the simulation files into a single file.
    """
    # set the dir name and read in the state files for merging
    dir_name = joinpath("results", "UoA_forecast_output", file_date)
    state_file_names = readdir(dir_name)
    # need to remove the non-output files
    states = []
    for f in state_file_names
        ind = findfirst("_",f)[1]
        if (f[1:ind-1] ∉ states) && (f[1:ind-1] != "UoA")
            push!(states, f[1:ind-1])
        end
    end

    # now we make sure the filenames are as we want
    for (i,f) in enumerate(states)
        g = f*"_"*file_date*"_sim.csv"
        states[i] = g
    end

    df_merged = DataFrame()
    df_tmp = DataFrame()

    # loop over the states and read in and add to the merged df
    for s in states
        df_tmp = CSV.read(dir_name*"/"*s, DataFrame)
        df_merged = [df_merged; df_tmp]
    end

    CSV.write(dir_name*"/"*"UoA_samples_"*file_date*".csv", df_merged)

    return nothing

end

function merge_TP_files(file_date)
    """
    Merge the TP files in a similar format to the overall simulation file.
    """
    # set the dir name and read in the state files for merging
    dir_name = joinpath("results", "UoA_forecast_output", file_date)
    state_file_names = readdir(dir_name)
    # need to remove the non-output files
    states = []
    for f in state_file_names
        ind = findfirst("_",f)[1]
        if (f[1:ind-1] ∉ states) && (f[1:ind-1] != "UoA")
            push!(states, f[1:ind-1])
        end
    end

    # now we make sure the filenames are as we want
    for (i,f) in enumerate(states)
        g = f*"_"*file_date*"_TP.csv"
        states[i] = g
    end

    df_merged = DataFrame()
    df_tmp = DataFrame()

    # loop over the states and read in and add to the merged df
    for s in states
        df_tmp = CSV.read(dir_name*"/"*s, DataFrame)
        df_merged = [df_merged; df_tmp]
    end

    CSV.write(dir_name*"/"*"UoA_TP_"*file_date*".csv", df_merged)

    return nothing

end

function summarise_forecast_for_plotting(D_observed)
    """
    Summarises the simulations into a dataframe.
    """

    df_observed_summary = DataFrame(
        "median" => median.(eachrow(D_observed)),
        "bottom" => quantile.(eachrow(D_observed), 0.025),
        "lower" => quantile.(eachrow(D_observed), 0.25),
        "upper" => quantile.(eachrow(D_observed), 0.75),
        "top" => quantile.(eachrow(D_observed), 0.975),
    )

    return df_observed_summary

end


text:
"""
    @filter_api api_name [filter_type=AbstractImageFilter]

This macro it generate two methods:

* `api_name([::Type,] img, f::filter_type, args...; kwargs...)`
* `api_name!([out,] img, f::filter_type, args...; kwargs...)`

For in-place method `api_name!`, `out` will be changed after calling the method.
When `out` is not explicitly passed, `img` will be changed after calling the method.

!!! info

    * Any api implementation needs to support a `f(out, in, args...)` method.
    * This macro is designed to be used in ImagesAPIs, and not downstream packages

## Example:

### 1. register the api in ImagesAPI.jl

```julia
abstract type AbstractImageNoise <: AbstractImageFilter end
@filter_api apply_noise AbstractImageNoise
```

### 2. implement the api in ImageNoise.jl
```julia
import Main.ImagesAPI: AbstractImageNoise, AbstractImageFilter, apply_noise, apply_noise!

export
    apply_noise, apply_noise!,
    AbstractImageNoise,
    AdditiveWhiteGaussianNoise

struct AdditiveWhiteGaussianNoise{T<:AbstractFloat} <: AbstractImageNoise
    mean::T
    std::T
end

function (noise::AdditiveWhiteGaussianNoise)(out, in::AbstractArray)
    @. out = in + noise.std * randn(eltype(out), size(in)) + noise.mean
end
```

### 3. user call the API in a consistent way
```julia
using ImageNoise
noise = AdditiveWhiteGaussianNoise(0.0, 0.1)

# simple usage
apply_noise(ones(3,3), noise)

# inplace changing
img = ones(3,3)
apply_noise!(img, noise)

# preallocation output
img = ones(3,3)
out = zeros(3, 3)
apply_noise!(out, img, noise)
```
"""
macro filter_api(func_name, filter_type = AbstractImageFilter)
    inplace_func_name = Symbol(String(func_name) * "!")
    # TODO: there's little performance improvement on inplce operation over normal one
    @eval begin
        function $(inplace_func_name)(out, img, f::$(filter_type), args...; kwargs...)
            f(out, img, args...; kwargs...)
            out
        end

        function $(inplace_func_name)(img, f::$(filter_type), args...; kwargs...)
            tmp = copy(img)
            f(img, tmp, args...; kwargs...)
            img
        end

        function $(func_name)(::Type{T}, img, f::$(filter_type), args...; kwargs...) where T
            out = Array{T}(undef, size(img))
            $(inplace_func_name)(out, img, f, args...; kwargs...)
            return out
        end

        function $(func_name)(img, f::$(filter_type), args...; kwargs...)
            $(func_name)(eltype(img), img, f, args...; kwargs...)
        end
    end
end


text:
#!/usr/bin/env julia

# Pandoc filter to convert all regular text to uppercase.
# Code, link URLs, etc. are not affected.

using PandocFiltersLiveJuliaCode

function caps(t, v, format, meta)
   (t == "Str") || return nothing
   return Str(uppercase(v))
end

if abspath(PROGRAM_FILE) == @__FILE__
   PandocFiltersLiveJuliaCode.filter(caps)
end


text:
# twobinomial.jl
# Two sample hypothesis tests for proportions
#
# Copyright (C) 2020 Vladimir Arnautov
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

export TwoSampleBinomialTest

struct TwoSampleBinomialTask
    x1::Int
    n1::Int
    x2::Int
    n2::Int
    δ::Float64
    ptype::Symbol
    htype::Symbol
    alpha::Float64
    method::Symbol
end

struct TwoSampleBinomialTest <: HypothesisTest
    test::TwoSampleBinomialTask
    est
    se
    ci
    result
end

function cihalpha(htype, alpha)
    if  htype == :equivalence || htype == :superiority
        cialpha = alpha * 2
    elseif htype == :equality
        cialpha = alpha
    else
        throw(ArgumentError("unknown hypothesis=$(htype)"))
    end
    cialpha
end

"""
    TwoSampleBinomialTest(x1::Real, n1::Real, x2::Real, n2::Real, δ::Real; ptype::Symbol, htype::Symbol, alpha::Real = 0.05, method::Symbol = :default)

Perform a two sample binomial hypothesis check via computing confidence intrval
with coverage `1 - alpha` for `equality` test or `1 - 2 * alpha` for `equivalence`
or `superiority` test.

Hypothesis type `htype` are:
- `:equality`
- `:equivalence`
- `:superiority`

Parameters type `ptype` are:
    - `:difference` : risk difference;
    - `:riskratio` : risk ratio;
    - `:oddratio` : odd ratio;

For each parameters can be used  following methods.

Possible values for `method` for `difference` are:

  - `:mn`
  - `:wald`
  - `:waldcc`
  - `:nhs`
  - `:nhscc`
  - `:ac`
  - `:ha`
  - `:mnmee`
  - `:mover`

Possible values for `method` for `riskratio` are:

  - `:mn`
  - `:cli`
  - `:li`
  - `:mover`

Possible values for `method` for `oddratio` are:

  - `:mn`
  - `:woolf`
  - `:awoolf`
  - `:mover`

"""
function TwoSampleBinomialTest(x1::Real, n1::Real, x2::Real, n2::Real, δ::Real; ptype::Symbol, htype::Symbol, alpha::Real = 0.05, method::Symbol = :default)
    test = TwoSampleBinomialTask(x1, n1, x2, n2, δ, ptype, htype, alpha, method)
    cialpha = cihalpha(htype, alpha)
    if ptype == :difference
        est, se, lci, uci = ci_diff(test; level = 1 - cialpha, method = method)
    elseif ptype == :oddratio
        est, se, lci, uci = ci_or(test; level = 1 - cialpha, method = method)
    elseif ptype == :riskratio
        est, se, lci, uci = ci_rr(test; level = 1 - cialpha, method = method)
    else
        throw(ArgumentError("unknown parameter=$(ptype)"))
    end
    result = false
    if htype == :equivalence
        if lci > -abs(δ) && uci < abs(δ) result = true end
    elseif htype == :equality
        if lci > δ || uci < δ result = true end
    elseif htype == :superiority
        if lci > δ result = true end
    end
    TwoSampleBinomialTest(test, est, se, (lci, uci), result)
end
################################################################################
function ci_prop_wilson(x::Int, n::Int, alpha::Real)
    z   = abs(quantile(Normal(), 1 - alpha / 2))
    p   = x / n
    d   = 1 + (z^2) / n
    se  = sqrt((p * (1 - p) + (z^2) / (4 * n)) / n) / d
    est = (p + (z^2) / (2 * n)) / d
    return est, se, est - z * se, est + z * se
end
function ci_prop_wilson_cc(x::Int, n::Int, alpha::Real)
    z = abs(quantile(Normal(), 1 - alpha / 2))
    p = x / n
    l = (2*n*p+z*z-1-z*sqrt(z*z-2-1/n+4*p*(n*(1-p)+1)))/2/(n+z*z)
    u = (2*n*p+z*z+1+z*sqrt(z*z+2-1/n+4*p*(n*(1-p)-1)))/2/(n+z*z)
    return p, NaN, min(p, l), max(p, u)
end
################################################################################
function ci_diff(test; level, method)
    if method == :mn || method == :default
        ci_diff_mn(test, level)
    elseif method == :wald
        ci_diff_wald(test, level)
    elseif method == :waldcc
        ci_diff_wald_cc(test, level)
    elseif method == :nhs
        ci_diff_nhs(test, level)
    elseif method == :nhscc
        ci_diff_nhs_cc(test, level)
    elseif method == :ac
        ci_diff_ac(test, level)
    elseif method == :ha
        ci_diff_ha(test, level)
    elseif method == :mnmee
        ci_diff_mnmee(test, level)
    elseif method == :mover
        ci_diff_mover(test, level)
    else
        throw(ArgumentError("unknown ci method=$(method)"))
    end
end
function ci_rr(test; level, method)
    if method == :mn || method == :default
        ci_rr_mn(test, level)
    elseif method == :cli
        ci_rr_cli(test, level)
    elseif method == :li
        ci_rr_li(test, level)
    elseif method == :mover
        ci_rr_mover(test, level)
    else
        throw(ArgumentError("unknown ci method=$(method)"))
    end
end
function ci_or(test; level, method)
    if method == :mn || method == :default
        ci_or_mn(test, level)
    elseif method == :woolf
        ci_or_woolf(test, level)
    elseif method == :awoolf
        ci_or_awoolf(test, level)
    elseif method == :mover
        ci_or_mover(test, level)
    else
        throw(ArgumentError("unknown ci method=$(method)"))
    end
end
################################################################################
#Not implementsd yet
function pvalue(tr::TwoSampleBinomialTest)
    if tr.test.method ∈ [:wald, :waldcc, :ac, :ha]
        #return ccdf(Normal(), z)*2.
        throw(ArgumentError("not implemented"))
    elseif tr.test.method ∈ [:mn, :nhs, :nhscc, :mnmee]
        #return NaN
        throw(ArgumentError("not implemented"))
    else
        throw(ArgumentError("unknown ci method=$(method)"))
    end
end
function confint(tr::TwoSampleBinomialTest)
    throw(ArgumentError("not implemented"))
end
################################################################################
#Method of Mee 1984 with Miettinen and Nurminen modification n / (n - 1) Newcombe 1998
#Score intervals for the difference of two binomial proportions
@inline function mle_diff(p1::Float64, n1::Int, p2::Float64, n2::Int, δ::Float64)
    if p1 - p2 - δ == 0 return 0.0 end
    θ = n2 / n1
    a = 1 + θ
    b = -(1 + θ + p1 + θ * p2 + δ * (θ + 2))
    c = δ^2 + δ * (2 * p1 + θ + 1) + p1 + θ * p2
    d = -p1 * δ * (1 + δ)
    v = (b / a / 3)^3 - b * c/(6 * a * a) + d / 2 / a
    u = sign(v) * sqrt((b / 3 / a)^2 - c / 3 / a)
    w = (pi + acos(v / u^3)) / 3
    p1n = 2 * u * cos(w) - b / 3 /a
    p2n = p1n - δ
    return p1n, p2n
end
@inline function mn_diff_z_val(p1::Real, n1::Int, p2::Real, n2::Int, est::Real, δ::Real)
    p1n, p2n = mle_diff(p1, n1, p2, n2, δ)
    return (est - δ)^2 / ((n1 + n2) / (n1 + n2 - 1) * (p1n * (1 - p1n) / n1 + p2n * (1 - p2n) / n2))
end
function ci_diff_mn(test, level; atol::Float64 = 1E-8)
    ests, ses, lcis, ucis = ci_diff_nhs_cc(test, level)
    p1       = test.x1 / test.n1
    p2       = test.x2 / test.n2
    est      = p1 - p2
    z        = quantile(Chisq(1), level)
    fmnd(x)  = mn_diff_z_val(p1, test.n1, p2, test.n2, est, x) - z

    if fmnd(lcis) * fmnd(est - eps()) < 0.0
        ll = lcis
        lu = est - eps()
    else
        ll = -1.0 + eps()
        lu = lcis
    end
    if fmnd(ucis) * fmnd(est + eps()) < 0.0
        ul = est + eps()
        uu = ucis
    else
        ul = ucis
        uu = 1.0 - eps()
    end
    lci = find_zero(fmnd, (ll, lu), atol = atol)
    uci = find_zero(fmnd, (ul, uu), atol = atol)
    return est, NaN, lci, uci
end
#Wald
function ci_diff_wald(test, level)
    alpha    = 1 - level
    p1       = test.x1 / test.n1
    p2       = test.x2 / test.n2
    est      = p1 - p2
    z        = quantile(Normal(), 1 - alpha / 2)
    se       = sqrt(p1 * (1 - p1) / test.n1 + p2 * (1 - p2) / test.n2)
    return est, se, est - z * se, est + z * se
end
#Wald continuity correction
function ci_diff_wald_cc(test, level)
    alpha    = 1 - level
    p1       = test.x1 / test.n1
    p2       = test.x2 / test.n2
    est      = p1 - p2
    z        = quantile(Normal(), 1 - alpha / 2)
    se       = sqrt(p1 * (1 - p1) / test.n1 + p2 * (1 - p2) / test.n2)
    cc       = 0.5*(1 / test.n1 + 1 / test.n2)
    return est, se, est - z * se - cc, est + z * se + cc
end
#Newcombes Hybrid (wilson) Score interval for the difference of proportions
#Newcombe 1998
function ci_diff_nhs(test, level)
    alpha    = 1 - level
    p1       = test.x1 / test.n1
    p2       = test.x2 / test.n2
    est      = p1 - p2
    z        = quantile(Normal(), 1 - alpha / 2)
    est1, se1, lci1, uci1 = ci_prop_wilson(test.x1, test.n1, alpha)
    est2, se2, lci2, uci2 = ci_prop_wilson(test.x2, test.n2, alpha)
    return est, NaN, est - z * sqrt(lci1 * (1 - lci1)/test.n1 + uci2 * (1 - uci2) / test.n2), est + z * sqrt(uci1 * (1 - uci1) / test.n1 + lci2 * (1 - lci2) / test.n2)
end
#Newcombes Hybrid Score continuity correction
function ci_diff_nhs_cc(test, level)
    alpha    = 1 - level
    p1       = test.x1 / test.n1
    p2       = test.x2 / test.n2
    est      = p1 - p2
    z        = quantile(Normal(), 1 - alpha / 2)
    est1, se1, lci1, uci1 = ci_prop_wilson_cc(test.x1, test.n1, alpha)
    est2, se2, lci2, uci2 = ci_prop_wilson_cc(test.x2, test.n2, alpha)
    return est, NaN, est - sqrt((p1 - lci1)^2 + (uci2 - p2)^2), est + sqrt((uci1 - p1)^2 + (p2 - lci2)^2)
end
#Agresti-Caffo interval for the difference of proportions
#Agresti A, Caffo B., “Simple and effective confidence intervals for proportions and differences of proportions result from adding two successes and two failures”, American Statistician 54: 280–288 (2000)
function ci_diff_ac(test, level)
    alpha    = 1 - level
    z        = quantile(Normal(), 1 - alpha / 2)
    p1I      = (test.x1 + 1) / (test.n1 + 2)
    p2I      = (test.x2 + 1) / (test.n2 + 2)
    n1I      = test.n1 + 2
    n2I      = test.n2 + 2
    est      = p1I - p2I
    se       = sqrt(p1I * (1 - p1I) / n1I + p2I * (1 - p2I) / n2I)
    return est, se, est - z * se, est + z * se
end
function ci_diff_ha(test, level)
    alpha    = 1 - level
    p1       = test.x1 / test.n1
    p2       = test.x2 / test.n2
    est      = p1 - p2
    z        = quantile(Normal(), 1 - alpha / 2)
    se       = sqrt(p1 * (1 - p1) / (test.n1 - 1) + p2 * (1 - p2) / (test.n2 - 1))
    cc       = 1 / min(test.n1, test.n2)
    return est, se, est - z * se - cc, est + z * se + cc
end
@inline function mnmee_z_val(p1::Real, n1::Int, p2::Real, n2::Int, est::Real, δ::Real)
    p1n, p2n = mle_diff(p1, n1, p2, n2, δ)
    return (est - δ)^2 / (p1n * (1 - p1n) / n1 + p2n * (1 - p2n) / n2)
end
#Mee RW (1984) Confidence bounds for the difference between two probabilities,Biometrics40:1175-1176
#MN - no correction
function ci_diff_mnmee(test, level; atol::Float64 = 1E-8)
    ests, ses, lcis, ucis = ci_diff_nhs_cc(test, level)
    alpha    = 1 - level
    p1       = test.x1 / test.n1
    p2       = test.x2 / test.n2
    est      = p1 - p2
    z        = quantile(Chisq(1), 1 - alpha)
    fmnd(x)  = mnmee_z_val(p1, test.n1, p2, test.n2, est, x) - z
    if fmnd(lcis) * fmnd(est - eps()) < 0.0
        ll = lcis
        lu = est - eps()
    else
        ll = -1.0 + eps()
        lu = lcis
    end
    if fmnd(ucis) * fmnd(est + eps()) < 0.0
        ul = est + eps()
        uu = ucis
    else
        ul = ucis
        uu = 1.0 - eps()
    end
    lci = find_zero(fmnd, (ll, lu), atol = atol)
    uci = find_zero(fmnd, (ul, uu), atol = atol)
    return est, NaN, lci, uci
end
#Brown, Li's Jeffreys
function ci_diff_jeffreys(test, level)
    p1   = (test.x1 + 0.5) / (test.n1 + 1)
    p2   = (test.x2 + 0.5) / (test.n2+1)
    se   = sqrt(p1*(1 - p1) / test.n1 + p2 * (1 - p2) / test.n2)
    z    = quantile(Normal(), 1 - alpha / 2)
    est  = p1 - p2
    return test.x1/test.n1-test.x2/test.n2, se, max(-1.0, est - z * se), min(1.0, est + z * se)
end
#Method of variance estimates recovery
function ci_diff_mover(test, level)
    alpha    = 1 - level
    p1       = test.x1 / test.n1
    p2       = test.x2 / test.n2
    est      = p1 - p2
    Z        = quantile(Normal(), 1 - alpha / 2)
    est1, se1, lci1, uci1   = ci_prop_wilson(test.x1, test.n1, alpha)
    est2, se2, lci2, uci2   = ci_prop_wilson(test.x2, test.n2, alpha)
    lci      = est - sqrt((p1 - lci1)^2 + (uci2 - p2)^2)
    uci      = est + sqrt((uci1 - p1)^2 + (p2 - lci2)^2)
    return est, NaN, lci, uci
end
################################################################################
@inline function mle_or(φ::Real, x1::Int, n1::Int, x2::Int, n2::Int)::Tuple{Float64,Float64}
    a  = n2 * (φ-1)
    b  = φ * n1 + n2 - (x1 + x2) * (φ - 1)
    c  = -(x1 + x2)
    p2 = (-b + sqrt(b * b - 4 * a * c)) / a / 2
    p1 = p2 * φ/(1 + p2 * (φ - 1))
    return p1, p2
end
@inline function mle_or_z_val(φ::Real, x1::Int, n1::Int, x2::Int, n2::Int)::Float64
    p1 = x1 / n1
    pmle1, pmle2 = mle_or(φ, x1, n1, x2, n2)
    return (n1 * (p1 - pmle1))^2 * (1 / (n1 * pmle1 * (1 - pmle1)) + 1/(n2 * pmle2 * (1 - pmle2))) / ((n1 + n2)/(n1 + n2 - 1))
end
################################################################################
#Miettinen O. S., Nurminen M. (1985) Comparative analysis of two rates.Statistics in Medicine4,213–226
#MN Score
function ci_or_mn(test, level; atol::Float64 = 1E-8)
    z        = quantile(Chisq(1), level)
    fmnor(x) = mle_or_z_val(x, test.x1, test.n1, test.x2, test.n2) - z
    if (test.x1==0 && test.x2==0) || (test.x1==test.n1 && test.x2==test.n2)
        return NaN, NaN, 0.0, Inf
    elseif test.x1==0 || test.x2==test.n2
        return 0.0, NaN, 0.0, find_zero(fmnor, 1e-6, atol = atol)
    elseif test.x1==test.n1 || test.x2 == 0
        return Inf, NaN, find_zero(fmnor, 1e-6, atol = atol), Inf
    else
        est  = (test.x1 / (test.n1 - test.x1)) / (test.x2 / (test.n2 - test.x2))
        return est, NaN, find_zero(fmnor, 1e-6, atol = atol), find_zero(fmnor, est+1e-6, atol = atol)
    end
end
#Woolf logit
#Woolf, B. (1955). On estimating the relation between blood group and disease. Annals of human genetics, 19(4):251-253.
function ci_or_woolf(test, level)
    alpha     = 1 - level
    xa        = test.x1
    xb        = test.n1 - test.x1
    xc        = test.x2
    xd        = test.n2 - test.x2
    est       = xa*xd/xc/xb
    estI      = log(est)
    se        = sqrt(1/xa + 1/xb + 1/xc + 1/xd)
    z         = quantile(Normal(), 1 - alpha / 2)
    return est, se, exp(estI - z * se), exp(estI + z * se)
end
#Adjusted Woolf interval (Gart adjusted logit) Lawson, R (2005):Smallsample confidence intervals for the odds ratio.  Communication in Statistics Simulation andComputation, 33, 1095-1113.
#Gart, J. J. (1966). Alternative analyses of contingency tables. Journal of the Royal Statistical Society. Series B (Methodological), 28:164-179.
function ci_or_awoolf(test, level)
    alpha     = 1 - level
    xa        = test.x1 + 0.5
    xb        = test.n1 - test.x1 + 0.5
    xc        = test.x2 + 0.5
    xd        = test.n2 - test.x2 + 0.5
    est       = xa*xd/xc/xb
    estI      = log(est)
    se        = sqrt(1/xa + 1/xb + 1/xc + 1/xd)
    z         = quantile(Normal(), 1 - alpha/2)
    return est, se, exp(estI - z * se), exp(estI + z * se)
end
#Method of variance estimates recovery
#Donner, A. and Zou, G. (2012). Closed-form confidence intervals for functions of the normal mean and standard deviation. Statistical Methods in Medical Research, 21(4):347-359.
function ci_or_mover(test, level)
    alpha    = 1 - level
    p1       = (test.x1/(test.n1-test.x1))
    p2       = (test.x2/(test.n2-test.x2))
    est      = p1/p2
    z        = quantile(Normal(), 1-alpha/2)
    est1, se1, lci1, uci1   = ci_prop_wilson(test.x1, test.n1, alpha)
    est2, se2, lci2, uci2   = ci_prop_wilson(test.x2, test.n2, alpha)
    vl1      = lci1/(1 - lci1)
    vu1      = uci1/(1 - uci1)
    vl2      = lci2/(1 - lci2)
    vu2      = uci2/(1 - uci2)
    lci      = (p1*p2-sqrt((p1*p2)^2 - vl1*vu2*(2*p1-vl1)*(2*p2-vu2)))/(vu2*(2*p2 - vu2))
    uci      = (p1*p2+sqrt((p1*p2)^2 - vu1*vl2*(2*p1-vu1)*(2*p2-vl2)))/(vl2*(2*p2 - vl2))
    return est, NaN, lci, uci
end
################################################################################
@inline function mle_rr(φ, x1::Int, n1::Int, x2::Int, n2::Int)
    a = (n1 + n2) * φ
    b = -(φ * (x1 + n2) + x2 + n1)
    c = x1 + x2
    p2 = (-b - sqrt(b * b - 4 * a * c)) / 2 / a
    p1 = p2 * φ
    return p1, p2
end
@inline function mle_rr_z_val(φ, x1::Int, n1::Int, x2::Int, n2::Int)
    p1 = x1 / n1
    p2 = x2 / n2
    pmle1, pmle2 = mle_rr(φ, x1, n1, x2, n2)
    return (p1 - φ * p2)^2 / ((pmle1 * (1 - pmle1) / n1 + φ * φ * pmle2 * (1 - pmle2) / n2) * ((n1 + n2 - 1) / (n1 + n2)))
end
################################################################################
#Miettinen-Nurminen Score interval
#Miettinen, O. and Nurminen, M. (1985), Comparative analysis of two rates. Statist. Med., 4: 213-226. doi:10.1002/sim.4780040211
function ci_rr_mn(test, level; atol::Float64 = 1E-8)
    z        = quantile(Chisq(1), level)
    fmnrr(x) = mle_rr_z_val(x, test.x1, test.n1, test.x2, test.n2) - z
    if (test.x1==0 && test.x2==0) || (test.x1==test.n1 && test.x2==test.n2)
        return NaN, NaN, 0.0, Inf
    elseif test.x1==0 || test.x2==test.n2
        return 0.0, NaN, 0.0, find_zero(fmnrr, 1e-8, atol=atol)
    elseif test.x1==test.n1 || test.x2 == 0
        return Inf, NaN, find_zero(fmnrr, 1e-8, atol=atol), Inf
    else
        est = (test.x1 / test.n1) / (test.x2 / test.n2)
        return est, NaN, find_zero(fmnrr, 1e-8, atol=atol), find_zero(fmnrr, est+1e-6, atol=atol)
    end
end
#Crude log interval
#Gart, JJand Nam, J (1988): Approximate interval estimation of the ratio of binomial parameters: Areview and corrections for skewness. Biometrics 44, 323-338.
function ci_rr_cli(test, level)
    alpha     = 1 - level
    x1I       = test.x1 + 0.5
    x2I       = test.x2 + 0.5
    n1I       = test.n1 + 0.5
    n2I       = test.n2 + 0.5
    estI      = log((x1I / n1I) / (x2I / n2I))
    se        = sqrt(1 / x2I + 1 / x1I - 1 / n2I - 1 / n1I)
    est       = (test.x1 / test.n1) / (test.x2 / test.n2)
    z         =  quantile(Normal(), 1 - alpha / 2)
    return est, se, exp(estI - z * se), exp(estI + z * se)
end
#Katz D, Baptista J, Azen SP and Pike MC. Obtaining confidence intervals for the risk ratio in cohort studies. Biometrics 1978; 34: 469–474
function ci_rr_li(test, level)
    alpha     = 1 - level
    est       = (test.x1 / test.n1) / ( test.x2 / test.n2)
    estI      = log(est)
    se        = sqrt(1 / test.x2 + 1 / test.x1 - 1 / test.n2 - 1 / test.n1)
    z         = quantile(Normal(), 1 - alpha / 2)
    return est, se, exp(estI - z * se), exp(estI + z * se)
end
#Method of variance estimates recovery (Donner, Zou, 2012)
function ci_rr_mover(test, level)
    alpha    = 1 - level
    p1       = test.x1 / test.n1
    p2       = test.x2 / test.n2
    est      = p1 / p2
    Z        = quantile(Normal(), 1 - alpha / 2)
    est1, se1, lci1, uci1   = ci_prop_wilson(test.x1, test.n1, alpha)
    est2, se2, lci2, uci2   = ci_prop_wilson(test.x2, test.n2, alpha)
    lci      = (p1 * p2 - sqrt((p1 * p2)^2 - lci1 * uci2 * (2 * p1 - lci1) * (2 * p2 - uci2))) / (uci2 * (2 * p2 - uci2))
    uci      = (p1 * p2 + sqrt((p1 * p2)^2 - uci1 * lci2 * (2 * p1 - uci1) * (2 * p2 - lci2))) / (lci2 * (2 * p2 - lci2))
    return est, NaN, lci, uci
end
################################################################################

function Base.show(_io::IO, res::TwoSampleBinomialTest)
    println(_io, "Two sample binomial hypothesis test")
    println(_io, "Parameter type: ", res.test.ptype)
    println(_io, "Hypothesis type: ", res.test.htype)
    println(_io, "A: ", res.test.x1, " / ", res.test.n1)
    println(_io, "B: ", res.test.x2, " / ", res.test.n2)
    println(_io, "Estimate: ", res.est)
    println(_io, "SE: ", res.se)

    cialpha = cihalpha(res.test.htype, res.test.alpha)

    println(_io, "$((1 - cialpha)*100)% CI: ", res.ci)
    println(_io, "H₀ rejected: ", res.result)
end

################################################################################

#function StatsBase.confint(x::TwoSampleBinomialTest; level::Float64=0.95, tail=:both, method=:default)
#end