text:
# Only for internal use.

# Indicates how `routing` should be understood:
# either as unitary flows (multiply by a demand matrix to write the capacity constraint)
# or directly flows (simple sum to write the capacity constraint)
# TODO: keep it as an enumaration or go to an abstract type?
@enum(FlowScaling, UnitaryFlows, UnscaledFlows)

mutable struct RoutingModel # TODO: use Unicode characters for Greek letters?
    data::RoutingData
    model::Model
    scaled_flows::FlowScaling

    ## Variables.
    mu # Either a maximum load or an oblivious ratio.
    tau # MMF objective (i.e. a threshold).
    routing # Main variable when optimising a routing.
    demand # When optimising for a traffic matrix (like oblivious subproblems).

    # For dual reformulations (robust and oblivious). Meaning heavily depends
    # on the actual model. These are not the dual values associated to the
    # constraints of the problem, but dual values coming from a reformulation.
    dual
    dual_alpha
    dual_beta
    dual_gamma
    dual_delta

    ## Constraints.
    constraints_capacity::Dict{Edge{Int}, Any} # Edge -> constraint ref
    constraints_mmf # Edge -> constraint ref (AbstractMatrix or Dict)

    # Flow-based.
    constraints_source_in::Dict{Edge{Int}, Any} # Demand -> constraint ref
    constraints_source_out::Dict{Edge{Int}, Any} # Demand -> constraint ref
    constraints_target_in::Dict{Edge{Int}, Any} # Demand -> constraint ref
    constraints_target_out::Dict{Edge{Int}, Any} # Demand -> constraint ref
    constraints_balance::Dict{Edge{Int}, Dict{Int, Any}} # Demand -> node -> constraint ref

    # Path-based.
    constraints_convexity # Demand -> constraint ref
    constraints_matrices::Dict{Dict{Edge{Int}, Float64}, Dict{Edge{Int}, Any}} # Matrix::Dict{Edge{Int}, Float64} -> edge -> constraint ref
    constraints_uncertainty_set::Dict{Edge{Int}, Dict{Edge{Int}, ConstraintRef}} # Edge -> demand -> constraint ref

    function RoutingModel(data::RoutingData, model::Model, scaled_flows::FlowScaling, routing;
        demand=nothing, mu=nothing, tau=nothing,
        dual=nothing, dual_alpha=nothing, dual_beta=nothing, dual_gamma=nothing, dual_delta=nothing,
        constraints_capacity::Dict{Edge{Int}, Any}=Dict{Edge{Int}, Any}(),
        constraints_mmf=nothing,
        constraints_source_in::Dict{Edge{Int}, Any}=Dict{Edge{Int}, Any}(), constraints_source_out::Dict{Edge{Int}, Any}=Dict{Edge{Int}, Any}(),
        constraints_target_in::Dict{Edge{Int}, Any}=Dict{Edge{Int}, Any}(), constraints_target_out::Dict{Edge{Int}, Any}=Dict{Edge{Int}, Any}(),
        constraints_balance::Dict{Edge{Int}, Dict{Int, Any}}=Dict{Edge{Int}, Dict{Int, Any}}(),
        constraints_convexity=nothing, constraints_matrices::Dict{Dict{Edge{Int}, Float64}, Dict{Edge{Int}, Any}}=Dict{Dict{Edge{Int}, Float64}, Dict{Edge{Int}, Any}}(),
        constraints_uncertainty_set::Dict{Edge{Int}, Dict{Edge{Int}, ConstraintRef}}=Dict{Edge{Int}, Dict{Edge{Int}, ConstraintRef}}()
    )
        @assert routing !== nothing

        return new(data, model, scaled_flows, mu, tau, routing, demand,
            dual, dual_alpha, dual_beta, dual_gamma, dual_delta,
            constraints_capacity, constraints_mmf,
            constraints_source_in, constraints_source_out, constraints_target_in, constraints_target_out, constraints_balance,
            constraints_convexity, constraints_matrices, constraints_uncertainty_set
        )
    end
end

function add_path!(rm::RoutingModel, demand::Edge, path::Vector{Edge})
    # TODO: propose to directly call `add_routing_var!`?
    return add_path!(rm.data, demand, path)
end

function add_routing_var!(rm::RoutingModel, demand::Edge, path_id::Int;
                          constraint_capacity::Bool=false)
    # Preconditions: the path already exists, this is a unitary-flow model,
    # this is a path-based formulation, this model has a convexity constraint.
    rd = rm.data
    @assert path_id in rd.demand_to_path_ids[demand]

    if ! (path_id in rd.demand_to_path_ids[demand])
        error("The path ID $(path_id) has not been added to the model.")
    end
    if rd.model_type.type != PathFormulation()
        error("This model is not a path-based formulation, but rather a $(rd.model_type). Impossible to add a new path in it.")
    end
    if length(rm.constraints_convexity) != n_demands(rd)
        error("This model does not look like it has a convexity constraint.")
    end

    # TODO: check for integer flows.
    # TODO: if constraint_capacity, ensure no ambiguity on traffic/capacity matrix.

    # Create the variable.
    if rm.scaled_flows == UnitaryFlows
        rm.routing[demand, path_id] = @variable(rm.model, lower_bound=0, upper_bound=1)
    else
        rm.routing[demand, path_id] = @variable(rm.model, lower_bound=0)
    end
    set_name(rm.routing[demand, path_id], "routing[$demand, $path_id]")

    # Add it into the convexity constraint.
    set_normalized_coefficient(rm.constraints_convexity[demand], rm.routing[demand, path_id], 1)
    # TODO: is it always 1? So far, yes...

    # Add it into the capacity constraint, if asked for.
    if constraint_capacity
        for e in rd.paths_edges[path_id]
            if rm.scaled_flows == UnitaryFlows
                coeff = rd.traffic_matrix[demand]
            elseif rm.scaled_flows == UnscaledFlows
                coeff = 1
            else
                error("Not implemented: $(rm.scaled_flows)")
            end

            set_normalized_coefficient(rm.constraints_capacity[e], rm.routing[demand, path_id], coeff)
        end
    end

    # TODO: other generic constraints (i.e. could be reused somewhere else)?
end

# function copy_model(rm::RoutingModel)
#     m2, refmap = copy_model(rm.model)
#     # TODO: translate more things? This is the bare minimum, and probably not even working in all needed cases...
#     return RoutingModel(rm.data, m2, refmap[rm.mu], refmap[rm.routing])
# end

# Channel the accessors through.
graph(rm::RoutingModel) = graph(rm.data)
edges(rm::RoutingModel) = edges(rm.data)
vertices(rm::RoutingModel) = vertices(rm.data)
demands(rm::RoutingModel) = edges(rm.data)

n_nodes(rm::RoutingModel) = n_nodes(rm.data)
n_edges(rm::RoutingModel) = n_edges(rm.data)
n_demands(rm::RoutingModel) = n_demands(rm.data)
n_paths(rm::RoutingModel) = n_paths(rm.data)

capacity(rm::RoutingModel, e::Edge{Int}) = capacity(rm.data, e)


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

insts = OrderedDict()
insts["minimal"] = [
    ((3, 2, true, 0),),
    ((3, 2, false, 0),),
    ((3, 2, true, 10),),
    ((3, 2, false, 10),),
    ]
insts["fast"] = [
    ((5, 3, true, 0),),
    ((5, 3, false, 0),),
    ((5, 3, true, 10),),
    ((5, 3, false, 10),),
    ((30, 10, true, 0),),
    ((30, 10, false, 0), (default_tol_relax = 1000,)),
    ((30, 10, true, 10),),
    ((30, 10, false, 10),),
    ]
insts["various"] = [
    ((10, 3, true, 10),),
    ((10, 3, false, 10),),
    ((20, 6, true, 10),),
    ((20, 6, false, 10),),
    ((40, 12, true, 10),),
    ((40, 12, false, 10),),
    ]
return (SparsePCANative, insts)

ext: jl


text:
module SIMDPoly

using SIMD

export evalpolydiv2x2c,
    packpoly2x2c,
    evalpoly1x2r,
    evalpoly1x4r,
    packpoly1x2r,
    packpoly1x4r

include("realpoly.jl")
include("complexpoly.jl")

end


text:
@testset "test_linalg" begin
a = [-3 2;1 1]
b = [0  1]'
c = [1 -1]
r = 3

@test care(a, b, c'*c, r) ≈ [0.5895174372762604 1.8215747248860816; 1.8215747248860823 8.818839806923107]
@test dare(a, b, c'*c, r) ≈ [240.73344504496302 -131.09928700089387; -131.0992870008943 75.93413176750603]

## Test dare method for non invertible matrices
A = [0 1; 0 0];
B0 = [0; 1];
Q = Matrix{Float64}(I, 2, 2);
R0 = 1
X = dare(A, B0, Q, R0);
# Reshape for matrix expression
B = reshape(B0, 2, 1)
R = fill(R0, 1, 1)
@test norm(A'X*A - X - (A'X*B)*((B'X*B + R)\(B'X*A)) + Q) < 1e-14
## Test dare for scalars
A = 1.0;
B = 1.0;
Q = 1.0;
R = 1.0;
X0 = dare(A, B, Q, R);
X = X0[1]
@test norm(A'X*A - X - (A'X*B)*((B'X*B + R)\(B'X*A)) + Q) < 1e-14
# Test for complex matrices
I2 = Matrix{Float64}(I, 2, 2)
@test dare([1.0 im; im 1.0], I2, I2, I2) ≈ [1 + sqrt(2) 0; 0 1 + sqrt(2)]
# And complex scalars
@test dare(1.0, 1, 1, 1) ≈ fill((1 + sqrt(5))/2, 1, 1)
@test dare(1.0im, 1, 1, 1) ≈ fill((1 + sqrt(5))/2, 1, 1)
@test dare(1.0, 1im, 1, 1) ≈ fill((1 + sqrt(5))/2, 1, 1)

## Test gram, ctrb, obsv
a_2 = [-5 -3; 2 -9]
C_212 = ss(a_2, [1; 2], [1 0; 0 1], [0; 0])
C_222 = ss(a_2, [1 0; 0 2], [1 0; 0 1], zeros(2,2))

a_3 = [-0.81  0.47 -0.43  1.58  0.25 -0.40  0.92;
       -0.61 -1.89  0.80 -1.59  2.01  0.98 -0.90;
        0.50 -1.18 -2.12 -1.56 -1.14  0.14 -0.87;
       -1.28  2.12  0.47 -1.16  3.65 -1.23 -1.26;
       -0.24 -0.08  1.58 -3.58 -1.26  1.72 -2.59;
        1.28 -0.96 -1.28 -0.57 -2.43 -2.43 -0.36;
       -0.16  1.53 -0.99  1.47  0.61 -2.22 -3.27]
b_3 =  [ 0.49 -0.80  0.10;
        -0.18  0.70  0.72;
         0.00  0.83  2.58;
         1.42  0.00  0.00;
         0.00  0.22  0.19;
         0.00  0.00 -0.08;
         0.00 -1.15 -1.93]
c_3 = [ 0.00  0.84  0.10  0.30  0.49  1.71  0.00
       -1.79 -0.89 -0.54 -0.60  0.74 -0.19  0.00]
d_3 = [ 1.36  0.96  1.44
       -1.07  0.00 -1.96]
C_732 = ss(a_3,b_3,c_3,d_3);

s = tf("s")
f_C_211 = (s+2)*(s+3)/((s+4)*(s+5))
# biquad passband
omega0 = 52.0; Q = 10
f_C_211_bis = (s/(Q*omega0)) / ((s/omega0)^2 + s/(Q*omega0) + 1)    # Works again in julia 1.0

C_22tf = [0 tf([3,0],[1,1,10]);tf([1,1],[1,5]) tf(2,[1,6])]

da_2 = [0.2 -0.8; -0.8 0.07]
D_221 = ss(da_2, [1 0; 0 2], [1 0], [0 0], 0.005)
D_222 = ss(da_2, [1 0; 0 2], [1 0; 0 1], zeros(2,2), 0.005)

a_4 = [-0.08 0.83 0.0 0.0;-0.83 -0.08 0.0 0.0;0.0 0.0 -0.7 9.0;0.0 0.0 -9.0 -0.7]
b_4 = [1 1;0 0;1 -1;0 0]
c_4 = [0.4 0.0 0.4 0.0;0.6 0.0 1.0 0.0]
d_4 = [0.3 0.0;0.0 -0.15]
D_422 = ss(a_4,b_4,c_4,d_4,0.2)
z = tf("z", 0.05)
D_311 = z^3/((z+0.5)*(z^2-1.4z+0.72801))

C_static = ss([-0.78593 -1.2707;  0.2638 -1.13143])
D_static = ss([0.704473 1.56483; -1.6661 -2.16041], 0.07)

@test gram(C_212,:c) ≈ [0.042016806722689065 0.09663865546218485
    0.09663865546218488 0.24369747899159663]
@test gram(C_212,:o) ≈ [0.09523809523809523 -0.0119047619047619
    -0.011904761904761911 0.059523809523809534]
@test gram(C_222,:c) ≈ [0.11764705882352941 -0.029411764705882346
    -0.02941176470588234 0.21568627450980393]
@test gram(C_222,:o) ≈ [0.09523809523809523 -0.0119047619047619
    -0.011904761904761911 0.059523809523809534]
@test gram(D_221,:c) ≈ [11.895658388619264 -7.526846649259635
    -7.526846649259634 12.517564258299078]
@test gram(D_221,:o) ≈ [3.1039794541519283 -1.7910988952221152
    -1.7910988952221152 2.197919733616834]
@test gram(D_222,:c) ≈ [11.895658388619264 -7.526846649259635
    -7.526846649259634 12.517564258299078]
@test gram(D_222,:o) ≈ [5.301899187768763 -3.2250358337314955
    -3.2250358337314955 4.777830864787395]

@test obsv(C_222) == [1 0; 0 1; -5 -3; 2 -9]
@test ctrb(C_222) == [1 0 -5 -6; 0 2 2 -18]
@test obsv(C_212) == [1 0; 0 1; -5 -3; 2 -9]
@test ctrb(C_212) == [1 -11; 2 -16]
@test norm(C_222) ≈ 0.5773502691896258
@test norm(C_212) ≈ 0.5345224838248488
@test norm(D_222) ≈ 4.940973856125768
@test norm(D_221) ≈ 3.4490083195926426
@test norm(ss([1],[2],[3],[4])) == Inf

#
## Test Hinfinity norm computations
#
# Continuous
Ninf, ω_peak = hinfnorm(ss(-1, 1, 1, 0))
@test Ninf ≈ 1 rtol=1e-6
@test ω_peak ≈ 0 rtol=1e-6

Ninf, ω_peak = hinfnorm(ss(-1, 1, 1, 1))
@test Ninf ≈ 2 rtol=1e-6
@test ω_peak ≈ 0 rtol=1e-6

Ninf, ω_peak = hinfnorm(ss(1.5))
@test Ninf ≈ 1.5 rtol=1e-6
@test ω_peak ≈ 0 rtol=1e-6

tolHinf = 1e-12
@test norm(C_212, Inf, tol=tolHinf) ≈ 0.242535625036333 atol=tolHinf
@test norm(C_222, Inf, tol=tolHinf) ≈ 0.242535625036333 atol=tolHinf
Ninf, ω_peak = hinfnorm(C_732, tol=tolHinf)
@test Ninf ≈ 4.899135403568278 atol=(10*tolHinf)
@test ω_peak ≈ 6.112977387441163 atol=1e-6
@test norm(f_C_211, Inf, tol=tolHinf) ≈ 1.0 atol=(2*tolHinf)
Ninf, ω_peak =  hinfnorm(f_C_211_bis, tol=tolHinf)
@test Ninf ≈ 1.0
@test ω_peak ≈ 52.0

# unstable system
sys = 1/(s-1)
@inferred hinfnorm(sys)
@test norm(sys, Inf, tol=tolHinf) ≈ Inf
Ninf, ω_peak = hinfnorm(sys; tol=tolHinf)
@test Ninf == Inf
@test isnan(ω_peak)
Ninf, ω_peak = linfnorm(sys; tol=tolHinf)
@test Ninf ≈ 1
@test ω_peak ≈ 0

Ninf, ω_peak = hinfnorm(C_static, tol=tolHinf)
@test Ninf ≈ 1.760164138376307 atol=1e-12
@test ω_peak ≈ 0 atol=1e-12
@inferred hinfnorm(C_static, tol=tolHinf)

Ninf, ω_peak = hinfnorm(C_22tf, tol=tolHinf)
@test Ninf ≈ 3.014974550173459 atol=(10*tolHinf)
@test ω_peak ≈ 3.162123338668049 atol=1e-8

# Poles on the imaginary axis
# Only integrator
Ninf, ω_peak = hinfnorm(ss(0.0, 1, 1, 0))
@test Ninf == Inf
@test ω_peak ≈ 0

# Integrator and lag
Ninf, ω_peak = hinfnorm(ss([-2.0 1; 0 0], [0; 1], [1 0], 0))
@test Ninf == Inf
@test ω_peak ≈ 0

# 1/(s^2 + 1)
Ninf, ω_peak = hinfnorm(ss([0 1.0; -1.0 0], [1; 0], [0 1], 0))
@test Ninf == Inf
@test ω_peak ≈ 1

# 1/(s - im), complex coefficients
Ninf, ω_peak = hinfnorm(ss(im, 1, 1, 0))
@test Ninf == Inf
@test ω_peak ≈ 1

# Complex coefficient systems
Ninf, ω_peak = hinfnorm(ss(-1+im, 1, 1, 0))
@test Ninf ≈ 1 rtol=1e-6
@test ω_peak ≈ 1 rtol=1e-3

Ninf, ω_peak = hinfnorm(ss(-1-im, 1, 1, 0))
@test Ninf ≈ 1 rtol=1e-6
@test ω_peak ≈ -1 rtol=1e-3

# 1/(s + 1)/(s - 1)
@test hinfnorm(ss([-1 1; 0 1], [0; 1], [1 0], 0))[1] == Inf
@test linfnorm(ss([-1 1; 0 1], [0; 1], [1 0], 0))[1] ≈ 1 rtol=1e-6

# Second order resonant system with peaks at ±sqrt(2)
G = ss([0 1; -1 -1], [0; 1], [1 0], 0)
Ninf, ω_peak = hinfnorm(G)
@test Ninf ≈ 2/sqrt(3) rtol=1e-6
@test ω_peak ≈ 1/sqrt(2) rtol=1e-3
@test hinfnorm(G, tol=1e-8)[1] ≈ 2/sqrt(3) rtol=1e-8
@test hinfnorm(G, tol=1e-10)[1] ≈ 2/sqrt(3) rtol=1e-10

# Same version as above but frequency shifted by -1, i.e., two peaks at -1 ± sqrt(2)
Ninf, ω_peak = hinfnorm(ss([0 1; -im -(1+2im)], [0; 1], [1 0], 0))
@test Ninf ≈ 2/sqrt(3) rtol=1e-6
@test abs(ω_peak+1) ≈ 1/sqrt(2) rtol=1e-3

# Same version as above but frequency shifted by +1, i.e., two peaks at 1 ± sqrt(2)
Ninf, ω_peak = hinfnorm( ss([0 1; im -1+2im], [0; 1], [1 0], 0))
@test Ninf ≈ 2/sqrt(3) rtol=1e-6
@test abs(ω_peak-1) ≈ 1/sqrt(2) rtol=1e-3

# System with mutiple resonance peaks
A = [-0.1*I        -Diagonal(1:10)
     Diagonal(1:10) zeros(10,10)]

Ninf, ω_peak = hinfnorm(ss(A, 0.1*[ones(10); zeros(10)], [zeros(1,10) ones(1,10)], 0))
@test Ninf ≈ 1.05871265 rtol=1e-5
@test ω_peak ≈ 0.9874 rtol=1e-4

Ninf, ω_peak = hinfnorm(ss(A, 0.1*[ones(4); 2; ones(5); zeros(10)], [zeros(1,10) ones(1,10)], 0))
@test Ninf ≈ 2.00918 rtol=1e-5
@test ω_peak ≈ 4.9981 rtol=1e-4

# Discrete Time
Ninf, ω_peak =  hinfnorm(ss(0.8, 1, 1, 0, 1))
@test Ninf ≈ 5 rtol=1e-6
@test ω_peak ≈ 0

Ninf, ω_peak =  hinfnorm(ss(-0.8, 1, 1, 0, 1))
@test Ninf ≈ 5 rtol=1e-6
@test ω_peak ≈ pi

Ninf, ω_peak = hinfnorm(D_221, tol=tolHinf)
@test Ninf ≈ 17.794697451669421 atol=(20*tolHinf)
@test ω_peak ≈ 0 atol=1e-8

Ninf, ω_peak = linfnorm(D_422, tol=tolHinf)
@test Ninf ≈ 3.360351099392252 atol=(10*tolHinf)
@test ω_peak ≈ 8.320643111730551 atol=1e-6
Ninf, ω_peak = hinfnorm(D_422, tol=tolHinf)
@test Ninf == Inf
@test isnan(ω_peak)

Ninf, ω_peak = hinfnorm(D_311, tol=tolHinf)
@test Ninf ≈ 4.458729529942810 atol=(10*tolHinf)
@test ω_peak ≈ 11.878021287349698 atol=1e-6

Ninf, ω_peak = hinfnorm(D_static, tol=tolHinf)
@test Ninf ≈ 3.205246234285972 atol=1e-12
@test ω_peak ≈ 0 atol=1e-12

# Pole on unit circle
Ninf, ω_peak = hinfnorm(ss(1, 1, 1, 0, 1))
@test Ninf == Inf
@test ω_peak ≈ 0

# Complex coefficients
Ninf, ω_peak =  hinfnorm(ss(0.8*exp(1im), 1, 1, 0, 1))
@test Ninf ≈ 5 rtol=1e-6
@test ω_peak ≈ 1

Ninf, ω_peak =  hinfnorm(ss(0.8*exp(-1im), 1, 1, 0, 1))
@test Ninf ≈ 5 rtol=1e-6
@test ω_peak ≈ -1

sys = hinfnorm(ss(0.8*exp(2.5im), 1, 1, 0, 1))
Ninf, ω_peak =  hinfnorm(ss(0.8*exp(2.5im), 1, 1, 0, 1))
@test Ninf ≈ 5 rtol=1e-6
@test ω_peak ≈ 2.5
Ninf, ω_peak = linfnorm(ss(0.8*exp(2.5im), 1, 1, 0, 1))
@test Ninf ≈ 5 rtol=1e-6
@test ω_peak ≈ 2.5

# Unstable discrete-time system
sys = ss(2, 3, 1, 1, 1)
Ninf, ω_peak = hinfnorm(sys)
@test Ninf == Inf
@test isnan(ω_peak)
Ninf, ω_peak = linfnorm(sys)
@test Ninf ≈ 2.0 rtol=1e-6
@test ω_peak ≈ 0

A = [1  100  10000; .01  1  100; .0001  .01  1]
T, P, B = balance(A)
# The scaling is BLAS dependent. However, the ratio should be the same on all
# machines. We just need to check that T == res * constant
res_diag = [512, 8, 0.0625]
constant = T[1]/res_diag[1]
@test all(diag(T) == res_diag * constant)
@test P == [1 0 0; 0 1 0; 0 0 1]
@test B ≈ [1.0 1.5625 1.220703125; 0.64 1.0 0.78125; 0.8192 1.28 1.0]
end


text:
RGBAU8 = AlphaColorValue{RGB{Ufixed8}, Ufixed8}
Color.rgba(r::Real, g::Real, b::Real, a::Real)    = AlphaColorValue(RGB{Float32}(r,g,b), float32(a))
rgbaU8(r::Real, g::Real, b::Real, a::Real)  = AlphaColorValue(RGB{Ufixed8}(r,g,b), ufixed8(a))

#GLPlot.toopengl{T <: AbstractRGB}(colorinput::Input{T}) = toopengl(lift(x->AlphaColorValue(x, one(T)), RGBA{T}, colorinput))
tohsva(rgba)     = AlphaColorValue(convert(HSV, rgba.c), rgba.alpha)
torgba(hsva)     = AlphaColorValue(convert(RGB, hsva.c), hsva.alpha)
tohsva(h,s,v,a)  = AlphaColorValue(HSV(float32(h), float32(s), float32(v)), float32(a))

Base.convert{T <: AbstractAlphaColorValue}(typ::Type{T}, x::AbstractAlphaColorValue) = AlphaColorValue(convert(RGB{eltype(typ)}, x.c), convert(eltype(typ), x.alpha))


text:
# Obtained from https://ssd.jpl.nasa.gov/?constants
"""
    JD
One Julian Day is equal to exactly 86400 seconds
"""
const JD = 86400            # Julian Day {s}

"""
    JY
One Julian Year is equal to exactly 365.25 Julian Days or 31,557,600 seconds
"""
const JY = 365.25*JD           # Julian Year {s}

"""
    AU
One Astronomical Unit is equal to exactly 1.495978707e8 km
"""
const AU = 1.495978707e8     # [km] Astronomical Unit

"""
    G
The universal gravitational constant is approximately equal to 6.7430e-20 km³/(kg*s²)
"""
const G = 6.67430e-20        # [km³/(kg*s²)]


text:
immutable Vector3f
    x::Cfloat
    y::Cfloat
    z::Cfloat
end

abstract type Vector2 end
immutable Vector2i <: Vector2
    x::Cint
    y::Cint
end
immutable Vector2u <: Vector2
    x::UInt32
    y::UInt32
end
immutable Vector2f <: Vector2
    x::Cfloat
    y::Cfloat
end

function Vector2(x::Real, y::Real)
    typex = typeof(x); typey = typeof(y)
    if typex <: AbstractFloat || typey <: AbstractFloat
        Vector2f(x, y)
    elseif typex <: Unsigned || typey <: Unsigned
        Vector2u(x, y)
    else
        Vector2i(x, y)
    end
end

function to_vec2u(vec::Vector2)
    if typeof(vec) != Vector2u
        return Vector2u(UInt32(vec.x), UInt32(vec.y))
    else
        vec
    end
end

function to_vec2f(vec::Vector2)
    if typeof(vec) != Vector2f
        Vector2f(Float32(vec.x), Float32(vec.y))
    else
        vec
    end
end

function to_vec2i(vec::Vector2)
    if typeof(vec) != Vector2i
        Vector2i(Int32(vec.x), Int32(vec.y))
    else
        vec
    end
end

function distance_squared(vec1::Vector2, vec2::Vector2)
    (vec2.x - vec1.x)^2 + (vec2.y - vec1.y)^2
end

function distance(vec1::Vector2, vec2::Vector2)
    sqrt(distance_squared(vec1, vec2))
end

function *(vec1::Vector2, vec2::Vector2)
    Vector2(vec1.x * vec2.x, vec1.y * vec2.y)
end

function +(vec1::Vector2, vec2::Vector2)
    Vector2(vec1.x + vec2.x, vec1.y + vec2.y)
end

function -(vec1::Vector2, vec2::Vector2)
    Vector2(vec1.x - vec2.x, vec1.y - vec2.y)
end

function /(vec1::Vector2, vec2::Vector2)
    Vector2(vec1.x / vec2.x, vec1.y / vec2.y)
end

Base.contains(obj, v::Vector2) = contains(obj, v.x, v.y)

import Base.eltype
function eltype{T <: Vector2}(::Type{T})
    @assert isleaftype(T)
    typeof(T(0,0).x)
end
eltype(v::Vector2) = eltype(typeof(v))

import Base.*
*(s::Number, v::Vector2) = Vector2(s*v.x, s*v.y)
*(v::Vector2, s::Number) = s*v

import Base.==
==(v::Vector2, w::Vector2) = (v.x == w.x) & (v.y == w.y)

import Base.convert
convert{V <: Vector2}(::Type{V}, v::Vector2) = V(v.x, v.y)

Base.start(v::Vector2) = 1
Base.next(v::Vector2, i) = (v[i], i+1)
Base.done(v::Vector2, i) = i > 2
function Base.getindex(v::Vector2, i)
    if i == 1
        return v.x
    elseif i == 2
        return v.y
    else
        throw(BoundsError(v,i))
    end
end
Base.length(v::Vector2) = 2


text:
# not run during tests, but reproduce teh effective sample size simulations from
# Fearnhead (2004).
using Distributed
addprocs()

@everywhere begin
    using ParallelDataTransfer, SharedArrays
    using Particles, Distributions, StatsBase, Random
    using Particles: ncomponents
    using Distributions: MixtureModel
    using ConjugatePriors: NormalInverseGamma, NormalInverseChisq

    function m_m2(ps::ParticleFilter)
        ns = ncomponents.(particles(ps))
        ws = Weights(weight.(particles(ps)))
        (mean(ns, ws), var(ns, ws))
    end

    function m_m2(gc::GibbsCRP, n; burnin=500)
        for _ in 1:burnin
            sample!(gc)
        end
        ns = Int[]
        for _ in 1:n
            sample!(gc)
            ncomp = length(gc.components) - length(gc.empties)
            push!(ns, ncomp)
        end
        (mean(ns), var(ns))
    end

    n_samp = 5000
    n_dat = 200
    α = 0.5
    prior = convert(NormalInverseChisq, NormalInverseGamma(0., 25., 1., 1.))
    stateprior = ChineseRestaurantProcess(α)

end

function ess(f::Function; n::Int=100)
    sendto(workers(), f=f)
    ms = SharedArray{Float64}(n)
    m2s = SharedArray{Float64}(n)
    @sync @distributed for i in 1:n
        Random.seed!(i)
        ms[i], m2s[i] = f()
    end
    return mean(m2s) / var(ms)
end

n_ess = 100
results = []
for σ in [0.5, 2.5]
    for μ in [0.5, 1, 2, 5]
        mm = MixtureModel([Normal(0, 0.5), Normal(μ, 0.5), Normal(2μ, σ)],
                          [1/2, 1/6, 1/3])
        x = rand(mm, n_dat)
        fc = ess(() -> m_m2(filter!(FearnheadParticles(n_samp, prior, stateprior), x, false)), n=n_ess)
        cl = ess(() -> m_m2(filter!(ChenLiuParticles(n_samp, prior, stateprior, rejuv=50.), x, false)), n=n_ess)
        gs = ess(() -> m_m2(GibbsCRP(prior, α, x), n_samp, burnin=500), n=n_ess)
        @show μ, σ, fc, cl, gs
        push!(results, (μ, σ, fc, cl, gs))
    end
end

using JLD2
@save "ess_results_$(DateTime(now())).jld2" results

# run using ESS formula from Fearnhead(2004)
# (μ, σ, fc, gs) = (0.5, 0.5, 707.4197999272996, 250.4401078254187)
# (μ, σ, fc, gs) = (1.0, 0.5, 4.986123628551047, 231.97429021977032)
# (μ, σ, fc, gs) = (2.0, 0.5, 24.706860608800266, 712.9339627024202)
# (μ, σ, fc, gs) = (5.0, 0.5, 2294.7314794665976, 517.5777000665694)
# (μ, σ, fc, gs) = (0.5, 2.5, 3.6015506679644673, 192.154516933825)
# (μ, σ, fc, gs) = (1.0, 2.5, 4.093719484784503, 317.6933662619361)
# (μ, σ, fc, gs) = (2.0, 2.5, 3.3444126682400066, 310.6637042075222)
# (μ, σ, fc, gs) = (5.0, 2.5, 9.241825755336542, 482.340345642232)

# run using mean(vars) / var(means)
# (μ, σ, fc, gs) = (0.5, 0.5, 117.49800834618993, 363.32269048299077)
# (μ, σ, fc, gs) = (1.0, 0.5, 6.220696850386848, 413.1589709680883)
# (μ, σ, fc, gs) = (2.0, 0.5, 19.496984108864652, 468.6555274028235)
# (μ, σ, fc, gs) = (5.0, 0.5, 4553.633313173106, 1009.8919191035734)
# (μ, σ, fc, gs) = (0.5, 2.5, 2.3728181956466163, 245.8546338400684)
# (μ, σ, fc, gs) = (1.0, 2.5, 1.751769156292295, 238.46679386421778)
# (μ, σ, fc, gs) = (2.0, 2.5, 3.9716813025135163, 373.4610563047813)
# (μ, σ, fc, gs) = (5.0, 2.5, 23.66274369175917, 321.8129868916317)


text:
using PackageCompiler

target_dir = get(ENV, "OUTDIR", "$(@__DIR__)/PLextCompiled")
target_dir = replace(target_dir, "\\"=>"/")       # Change Windows paths to use "/"

package_dir = get(ENV, "PKGDIR", "$(@__DIR__)")

println("Creating library in $target_dir")
PackageCompiler.create_library(package_dir, target_dir;
                                lib_name="plext",
                                precompile_execution_file=["$(@__DIR__)/generate_precompile.jl"],
                                precompile_statements_file=["$(@__DIR__)/additional_precompile.jl"],
                                incremental=false,
                                filter_stdlibs=true,
                                header_files = ["$(@__DIR__)/plext.h", "$(@__DIR__)/jl_init.h"],
                                force=true
                            )

max_stars_repo_name: mlxd/PLext

max_issues_repo_name: mlxd/PLext

max_forks_repo_name: mlxd/PLext


text:
abstract type InterpretableModelOptimizer end

abstract type HeuristicOptimizer <: InterpretableModelOptimizer end

fit!(m::InterpretableModel, o::InterpretableModelOptimizer) = notdef((m, o), :fit!)


text:
"""
    binary2OneOrNegativeOne(val)

Converts an integer that is either 0 or 1 to -1 or 1. Values that are 0 are
returned as -1 and values that are 1 are returned as 1.

Required Arguments:

- `val`: integer value that is either 0 or 1

Returns:

- a value that is either -1 or 1
"""
function binary2OneOrNegativeOne(val)
    if val == 1
        return 1
    elseif val == 0
        return -1
    else
        error("Value $(val) must be either 0 or 1.")
    end
end

"""
    calc_code_val(signal::ReplicaSignal, t)

Calculates the value of the generic code with parameters defined by `signal` at
a time specified by `t` in seconds. Returns a Complex{Int} value that is either
-1 or 1 and a phase in `rad`. Use the `code_ϕ`, or carrier phase` result to
perform QPSK correctly.

Required Arguments:

- `signal::ReplicaSignal`: struct containing signal data and parameters
- `t::Float64`: current time in seconds

Returns:

- `(code_val, code_ϕ)::Tuple` where
    * `code_val::Complex{Int}`: unnormalized complex number with the values of
                                the I and Q channel codes being the real and
                                imaginary parts of this value, respectivly
    * `code_ϕ::Float64`: the carrier phase in `rad` corresponding to the values
                         of the I and Q channel codes

The carrier phase value for a given set of I and Q channel code values at a
given time `t` is given by:

```math
    code_\\phi = \\tan\\left(\\frac{Q}{I}\\right)^{-1}
```

where for codes on **both** the I and Q channel the phase is:

| I Code Value  | Q Code Value  | Phase (degrees)  |
|--------------:|--------------:|-----------------:|
| 1             | 1             | 45               |
| -1            | 1             | 135              |
| -1            | -1            | 225              |
| 1             | -1            | 315              |

and for codes on **only** the I channel:

| I Code Value  | Q Code Value  | Phase (degrees)  |
|--------------:|--------------:|-----------------:|
| 1             | 0             | 0                |
| -1            | 0             | 180              |

and for codes on **only** the Q channel:

| I Code Value  | Q Code Value  | Phase (degrees) |
|--------------:|--------------:|----------------:|
| 0             | 1             | 90              |
| 0             | -1            | 270             |
"""
function calc_code_val(signal::ReplicaSignal, t)
    prn = signal.prn
    # I channel
    # If I channel codes exist, Multiply the value of each layer of code at a given
    # `t` with the value of `I_val`
    if signal.signal_type.include_I
        # I_val = 0
        I_val = 1
        I_codes = signal.signal_type.I_codes
        # Determine whether databits are included in I channel and whether
        # `include_databits_I` flag is true. If so, databits are multiplied by
        # the value of `I_val`. If this flag is set to false or if
        # `signal.isreplica` is set to true, no databits are included.
        if I_codes.databits
            if ~signal.isreplica && signal.include_databits_I
                code_num = I_codes.code_num
            else
                code_num = I_codes.code_num - 1
            end
        else
            code_num = I_codes.code_num
        end
        # Loop through I channel codes and multiply current value at time `t` to
        # value of `I_val`
        for i in 1:code_num
            code_idx = calccodeidx(signal.init_code_phases_I[i],
                                   signal.f_code_d_I[i], signal.f_code_dd_I[i],
                                   t, I_codes.code_lengths[i])
            # I_val = xor(I_val, I_codes.codes[i][prn][code_idx])
            I_val *= binary2OneOrNegativeOne(I_codes.codes[i][prn][code_idx])
        end
        # Change from 0:1 state to -1:1 state. This prevensts divide over zero
        # error in `atan` function at the end of this function.
        # I_val = binary2OneOrNegativeOne(I_val)
    else
        I_val = 0
    end
    # Q channel
    # If Q channel codes exist, Multiply the value of each layer of code at a given
    # `t` with the value of `Q_val`
    if signal.signal_type.include_Q
        # Q_val = 0
        Q_val = 1
        Q_codes = signal.signal_type.Q_codes
        # Determine whether databits are included in Q channel and whether
        # `include_databits_Q` flag is true. If so, databits are multiplied by
        # the value of `Q_val`. If this flag is set to false or if
        # `signal.isreplica` is set to true, no databits are included.
        if Q_codes.databits
            if ~signal.isreplica && signal.include_databits_Q
                code_num = Q_codes.code_num
            else
                code_num = Q_codes.code_num - 1
            end
        else
            code_num = Q_codes.code_num
        end
        # Loop through Q channel codes and multiply current value at time `t` to
        # value of `Q_val`
        for i in 1:code_num
            code_idx = calccodeidx(signal.init_code_phases_Q[i],
                                   signal.f_code_d_Q[i], signal.f_code_dd_Q[i],
                                   t, Q_codes.code_lengths[i])
            # Q_val = xor(Q_val, Q_codes.codes[i][prn][code_idx])
            Q_val *= binary2OneOrNegativeOne(Q_codes.codes[i][prn][code_idx])
        end
        # Change from 0:1 state to -1:1 state. This prevensts divide over zero
        # error in `atan` function at the end of this function.
        # Q_val = binary2OneOrNegativeOne(Q_val)
    else
        Q_val = 0
    end
    # Calculate the phase of the carrier based off the values of `Q_val` and
    # `I_val`. This simulates a QPSK signal.
    code_ϕ = atan(Q_val, I_val)
    # Construct a `Complex{Int}` that has the value of the I channel codes in
    # `real(code_val)` and the value of the Q channel codes in `iamg(code_val)`.
    code_val = complex(I_val, Q_val)
    return (code_val, code_ϕ)
end

"""
    calc_code_val(signal::ReplicaSignal, t, code_chips_I, code_chips_Q)

Calculates the value of the generic code with parameters defined by `signal` at
a time specified by `t` in seconds. Returns a Complex{Int} value that is either
-1 or 1 and a phase in `rad`. Use the `code_ϕ`, or carrier phase` result to
perform QPSK correctly.

This method accepts arrays `code_chips_I` and `code_chips_Q`, which are arrays
containing interpolated functions that calculate the current code chip as a
function of `t` in seconds. This method is used when the user supplies
`generatesignal!` with `doppler_curve` and `doppler_t` arrays, where
`doppler_curve` is assumed to be a non-linear Doppler frequency curve.

Required Arguments:

- `signal::ReplicaSignal`: struct containing signal data and parameters
- `t::Float64`: current time in seconds
- `code_chips_I`: N element array containing interpolated functions
                  that compute the current code chip for a given `t`
                  for each code on the I channel
- `code_chips_Q`: M element array containing interpolated functions
                  that compute the current code chip for a given `t`
                  for each code on the Q channel

Returns:

- `(code_val, code_ϕ)::Tuple` where
    * `code_val::Complex{Int}`: unnormalized complex number with the values of
                                the I and Q channel codes being the real and
                                imaginary parts of this value, respectivly
    * `code_ϕ::Float64`: the carrier phase in `rad` corresponding to the values
                         of the I and Q channel codes

The carrier phase value for a given set of I and Q channel code values at a
given time `t` is given by:

```math
    code_\\phi = \\tan\\left(\\frac{Q}{I}\\right)^{-1}
```

where for codes on **both** the I and Q channel the phase is:

| I Code Value  | Q Code Value  | Phase (degrees)  |
|--------------:|--------------:|-----------------:|
| 1             | 1             | 45               |
| -1            | 1             | 135              |
| -1            | -1            | 225              |
| 1             | -1            | 315              |

and for codes on **only** the I channel:

| I Code Value  | Q Code Value  | Phase (degrees)  |
|--------------:|--------------:|-----------------:|
| 1             | 0             | 0                |
| -1            | 0             | 180              |

and for codes on **only** the Q channel:

| I Code Value  | Q Code Value  | Phase (degrees) |
|--------------:|--------------:|----------------:|
| 0             | 1             | 90              |
| 0             | -1            | 270             |
"""
function calc_code_val(signal::ReplicaSignal, t, code_chips_I, code_chips_Q)
    prn = signal.prn
    # I channel
    # If I channel codes exist, Multiply the value of each layer of code at a
    # given `t` with the value of `I_val`
    if signal.signal_type.include_I
        # I_val = 0
        I_val = 1
        I_codes = signal.signal_type.I_codes
        # Determine whether databits are included in I channel and whether
        # `include_databits_I` flag is true. If so, databits are multiplied by
        # the value of `I_val`. If this flag is set to false or if
        # `signal.isreplica` is set to true, no databits are included.
        if I_codes.databits
            if ~signal.isreplica && signal.include_databits_I
                code_num = I_codes.code_num
            else
                code_num = I_codes.code_num - 1
            end
        else
            code_num = I_codes.code_num
        end
        # Loop through I channel codes and multiply current value at time `t` to
        # value of `I_val`
        for i in 1:code_num
            code_chip = code_chips_I[i](t)
            code_idx = Int(floor(code_chip%signal.signal_type.I_codes.code_lengths[i])) + 1
            # I_val = xor(I_val, I_codes.codes[i][prn][code_idx])
            I_val *= binary2OneOrNegativeOne(I_codes.codes[i][prn][code_idx])
        end
        # Change from 0:1 state to -1:1 state. This prevensts divide over zero
        # error in `atan` function at the end of this function.
        # I_val = binary2OneOrNegativeOne(I_val)
    else
        I_val = 0
    end
    # Q channel
    # If Q channel codes exist, Multiply the value of each layer of code at a
    # given `t` with the value of `Q_val`
    if signal.signal_type.include_Q
        # Q_val = 0
        Q_val = 1
        Q_codes = signal.signal_type.Q_codes
        # Determine whether databits are included in Q channel and whether
        # `include_databits_Q` flag is true. If so, databits are multiplied by
        # the value of `Q_val`. If this flag is set to false or if
        # `signal.isreplica` is set to true, no databits are included.
        if Q_codes.databits
            if ~signal.isreplica && signal.include_databits_Q
                code_num = Q_codes.code_num
            else
                code_num = Q_codes.code_num - 1
            end
        else
            code_num = Q_codes.code_num
        end
        # Loop through Q channel codes and multiply current value at time `t` to
        # value of `Q_val`
        for i in 1:code_num
            code_chip = code_chips_Q[i](t)
            code_idx = Int(floor(code_chip%signal.signal_type.Q_codes.code_lengths[i])) + 1
            # Q_val = xor(Q_val, Q_codes.codes[i][prn][code_idx])
            Q_val *= binary2OneOrNegativeOne(Q_codes.codes[i][prn][code_idx])
        end
        # Change from 0:1 state to -1:1 state. This prevensts divide over zero
        # Q_val = binary2OneOrNegativeOne(Q_val)
    else
        Q_val = 0
    end
    # Calculate the phase of the carrier based off the values of `Q_val` and
    # `I_val`. This simulates a QPSK signal.
    code_ϕ = atan(Q_val, I_val)
    # Construct a `Complex{Int}` that has the value of the I channel codes in
    # `real(code_val)` and the value of the Q channel codes in `iamg(code_val)`.
    code_val = complex(I_val, Q_val)
    return (code_val, code_ϕ)
end


text:
# Main GaussianProcess type

abstract type GPBase end

"""
    The abstract CovarianceStrategy type is for types that control how
    the covariance matrices and their positive definite representation
    are obtained or approximated. See SparseStrategy for examples.
"""
abstract type CovarianceStrategy end
struct FullCovariance <: CovarianceStrategy end
KernelData(k::Kernel, X1::AbstractMatrix, X2::AbstractMatrix, covstrat::CovarianceStrategy) = EmptyData()
KernelData(k::Kernel, X1::AbstractMatrix, X2::AbstractMatrix, covstrat::FullCovariance) = KernelData(k, X1, X2)
function alloc_cK(covstrat::CovarianceStrategy, nobs)
    # create placeholder PDMat
    m = Matrix{Float64}(undef, nobs, nobs)
    chol = Matrix{Float64}(undef, nobs, nobs)
    cK = PDMats.PDMat(m, Cholesky(chol, 'U', 0))
    return cK
end

#===============================
  Predictions
================================#
function predictMVN!(Kxx, Kff, Kfx, mx, αf)
    mu = mx + Kfx' * αf
    Lck = whiten!(Kff, Kfx)
    subtract_Lck!(Kxx, Lck)
    return mu, Kxx
end

"""
        predictMVN(xpred::AbstractMatrix, xtrain::AbstractMatrix, ytrain::AbstractVector,
                   kernel::Kernel, meanf::Mean, alpha::AbstractVector,
                   covstrat::CovarianceStrategy, Ktrain::AbstractPDMat)

Compute predictions using the standard multivariate normal conditional distribution formulae.
"""
function predictMVN(xpred::AbstractMatrix, xtrain::AbstractMatrix, ytrain::AbstractVector,
                   kernel::Kernel, meanf::Mean, alpha::AbstractVector,
                   covstrat::CovarianceStrategy, Ktrain::AbstractPDMat)
    crossdata = KernelData(kernel, xtrain, xpred)
    priordata = KernelData(kernel, xpred, xpred)
    Kcross = cov(kernel, xtrain, xpred, crossdata)
    Kpred = cov(kernel, xpred, xpred, priordata)
    mx = mean(meanf, xpred)
    mu, Sigma_raw = predictMVN!(Kpred, Ktrain, Kcross, mx, alpha)
    return mu, Sigma_raw
end

@inline function subtract_Lck!(Sigma_raw::AbstractArray{<:AbstractFloat}, Lck::AbstractArray{<:AbstractFloat})
    LinearAlgebra.BLAS.syrk!('U', 'T', -1.0, Lck, 1.0, Sigma_raw)
    LinearAlgebra.copytri!(Sigma_raw, 'U')
end
@inline subtract_Lck!(Sigma_raw, Lck) = Sigma_raw .-= Lck'Lck

"""
    predict_f(gp::GPBase, X::Matrix{Float64}[]; full_cov::Bool = false)

Return posterior mean and variance of the Gaussian Process `gp` at specfic points which are
given as columns of matrix `X`. If `full_cov` is `true`, the full covariance matrix is
returned instead of only variances.
"""
function predict_f(gp::GPBase, x::AbstractMatrix; full_cov::Bool=false)
    size(x,1) == gp.dim || throw(ArgumentError("Gaussian Process object and input observations do not have consistent dimensions"))
    if full_cov
        return predict_full(gp, x)
    else
        ## Calculate prediction for each point independently
        μ = Array{eltype(x)}(undef, size(x,2))
        σ2 = similar(μ)
        for k in 1:size(x,2)
            m, sig = predict_full(gp, x[:,k:k])
            μ[k] = m[1]
            σ2[k] = max(diag(sig)[1], 0.0)
        end
        return μ, σ2
    end
end

# 1D Case for prediction of process
predict_f(gp::GPBase, x::AbstractVector, args...; kwargs...) = predict_f(gp, x', args...; kwargs...)
# 1D Case for prediction of observations
predict_y(gp::GPBase, x::AbstractVector, args...; kwargs...) = predict_y(gp, x', args...; kwargs...)
@deprecate predict predict_y

wrap_cK(cK::PDMat, Σbuffer, chol::Cholesky) = PDMat(Σbuffer, chol)
mat(cK::PDMat) = cK.mat
cholfactors(cK::PDMat) = cK.chol.factors
"""
    make_posdef!(m::Matrix{Float64}, chol_factors::Matrix{Float64})

Try to encode covariance matrix `m` as a positive definite matrix.
The `chol_factors` matrix is recycled to store the cholesky decomposition,
so as to reduce the number of memory allocations.

Sometimes covariance matrices of Gaussian processes are positive definite mathematically
but have negative eigenvalues numerically. To resolve this issue, small weights are added
to the diagonal (and hereby all eigenvalues are raised by that amount mathematically)
until all eigenvalues are positive numerically.
"""
function make_posdef!(m::AbstractMatrix, chol_factors::AbstractMatrix)
    n = size(m, 1)
    size(m, 2) == n || throw(ArgumentError("Covariance matrix must be square"))
    for _ in 1:10 # 10 chances
        try
            # return m, cholesky(m)
            copyto!(chol_factors, m)
            chol = cholesky!(Symmetric(chol_factors, :U))
            return m, chol
        catch err
            if typeof(err)!=LinearAlgebra.PosDefException
                throw(err)
            end
            # that wasn't (numerically) positive definite,
            # so let's add some weight to the diagonal
            ϵ = 1e-6 * tr(m) / n
            @inbounds for i in 1:n
                m[i, i] += ϵ
            end
        end
    end
    copyto!(chol_factors, m)
    chol = cholesky!(Symmetric(chol_factors, :U))
    return m, chol
end
function make_posdef!(m::AbstractMatrix)
    chol_buffer = similar(m)
    return make_posdef!(m, chol_buffer)
end

#———————————————————————————————————————————————————————————
# Sample random draws from the GP
function Random.rand!(gp::GPBase, x::AbstractMatrix, A::DenseMatrix)
    nobs = size(x,2)
    n_sample = size(A,2)

    if gp.nobs == 0
        # Prior mean and covariance
        μ = mean(gp.mean, x);
        Σraw = cov(gp.kernel, x, x);
        Σraw, chol = make_posdef!(Σraw)
        Σ = PDMat(Σraw, chol)
    else
        # Posterior mean and covariance
        μ, Σraw = predict_f(gp, x; full_cov=true)
        Σraw, chol = make_posdef!(Σraw)
        Σ = PDMat(Σraw, chol)
    end
    return broadcast!(+, A, μ, unwhiten!(Σ,randn(nobs, n_sample)))
end

Random.rand(gp::GPBase, x::AbstractMatrix, n::Int) = rand!(gp, x, Array{Float64}(undef, size(x, 2), n))

# Sample from 1D GPBase
Random.rand(gp::GPBase, x::AbstractVector, n::Int) = rand(gp, x', n)

# Generate only one sample from the GPBase and returns a vector
Random.rand(gp::GPBase, x::AbstractMatrix) = vec(rand(gp,x,1))
Random.rand(gp::GPBase, x::AbstractVector) = vec(rand(gp,x',1))

# Sample random draws from the GP with a defined random number generator
function Random.rand!(rng::AbstractRNG, gp::GPBase, x::AbstractMatrix, A::DenseMatrix)
    nobs = size(x,2)
    n_sample = size(A,2)

    if gp.nobs == 0
        # Prior mean and covariance
        μ = mean(gp.mean, x);
        Σraw = cov(gp.kernel, x, x);
        Σraw, chol = make_posdef!(Σraw)
        Σ = PDMat(Σraw, chol)
    else
        # Posterior mean and covariance
        μ, Σraw = predict_f(gp, x; full_cov=true)
        Σraw, chol = make_posdef!(Σraw)
        Σ = PDMat(Σraw, chol)
    end
    return broadcast!(+, A, μ, unwhiten!(Σ,randn(rng, nobs, n_sample)))
end

Random.rand(rng::AbstractRNG, gp::GPBase, x::AbstractMatrix, n::Int) = rand!(rng, gp, x, Array{Float64}(undef, size(x, 2), n))

# Sample from 1D GPBase
Random.rand(rng::AbstractRNG, gp::GPBase, x::AbstractVector, n::Int) = rand(rng, gp, x', n)

# Generate only one sample from the GPBase and returns a vector
Random.rand(rng::AbstractRNG, gp::GPBase, x::AbstractMatrix) = vec(rand(rng, gp,x,1))
Random.rand(rng::AbstractRNG, gp::GPBase, x::AbstractVector) = vec(rand(rng, gp,x',1))


text:
module TestGLM

using Test

using MLJBase
using RDatasets
using LinearAlgebra

import MLJModels
import Distributions
import GLM
using MLJModels.GLM_
using Random: seed!
using RDatasets

###
### OLSREGRESSOR
###

boston = dataset("MASS", "Boston")
X = MLJBase.selectcols(boston, [:Crim, :Zn, :Indus, :NOx, :Rm, :Age,
                                :Dis, :Rad, :Tax, :PTRatio, :Black,
                                :LStat])
y = MLJBase.selectcols(boston, :MedV)

train, test = partition(eachindex(y), 0.7)

atom_ols = LinearRegressor()

Xtrain = selectrows(X, train)
ytrain = selectrows(y, train)
Xtest  = selectrows(X, test)

fitresult, _, report = fit(atom_ols, 1, Xtrain, ytrain)
θ = MLJBase.fitted_params(atom_ols, fitresult)

p = predict_mean(atom_ols, fitresult, Xtest)

# hand made regression to compare

Xa = MLJBase.matrix(X) # convert(Matrix{Float64}, X)
Xa1 = hcat(Xa, ones(size(Xa, 1)))
coefs = Xa1[train, :] \ y[train]
p2 = Xa1[test, :] * coefs

@test p ≈ p2

infos = info_dict(atom_ols)

@test infos[:name] == "LinearRegressor"
@test infos[:package_name] == "GLM"
@test infos[:is_pure_julia]
@test infos[:is_supervised]
@test infos[:package_license] == "MIT"
@test infos[:prediction_type] == :probabilistic
@test infos[:hyperparameters] == (:fit_intercept, :allowrankdeficient)
@test infos[:hyperparameter_types] == ("Bool", "Bool")

p_distr = predict(atom_ols, fitresult, selectrows(X, test))

@test p_distr[1] == Distributions.Normal(p[1], GLM.dispersion(fitresult))

###
### Logistic regression
###

seed!(0)

# data drawn from https://stats.idre.ucla.edu/r/dae/poisson-regression/
data = dataset("MASS", "Melanoma")

X = data[:, [:Status, :Sex, :Age, :Year, :Thickness]]
y_plain = data.Ulcer
y = categorical(y_plain)

n = length(y)

baseline_y = convert.(Int, rand(n) .> 0.5)
baseline_mse = sum((baseline_y - y_plain).^2)/n

@test baseline_mse ≤ 0.55

lr = LinearBinaryClassifier()
fitresult, _, report = fit(lr, 1, X, y)

p_mean  = predict_mean(lr, fitresult, X)
p_mode1 = convert.(Int, p_mean .> 0.5)
@test sum((p_mode1 - y_plain).^2)/n < 0.26

p_mode = predict_mode(lr, fitresult, X)

@test p_mode1 == convert.(Int, p_mode)

pr = LinearBinaryClassifier(link=GLM.ProbitLink())
fitresult, _, report = fit(pr, 1, X, y)
p_mode = convert.(Int, predict_mode(pr, fitresult, X))
@test sum((p_mode - y_plain).^2)/n < 0.26

# info_dict(atom_glmcount)

###
### Count regression
###

seed!(1512)

X = randn(500, 5)
θ = randn(5)
y = rand.(Distributions.Poisson.(exp.(X*θ)))

XTable = MLJBase.table(X)

lcr = LinearCountRegressor(fit_intercept=false)
fitresult, _, _ = fit(lcr, 1, XTable, y)

θ̂ = fitted_params(lcr, fitresult).coef

@test norm(θ̂ .- θ)/norm(θ) ≤ 5e-3

end # module
true


text:
struct SRKFTUIntermediate{T}
    x_apri::Vector{T}
    p_apri::Matrix{T}
    qr_zeros::Vector{T}
    qr_space::Vector{T}
    puft_vcat_q::Matrix{T}
end

function SRKFTUIntermediate(T::Type, num_x::Number)
    qr_zeros = zeros(T, 2 * num_x)
    puft_vcat_q = Matrix{T}(undef, 2 * num_x, num_x)
    qr_space_length = calc_gels_working_size(puft_vcat_q, qr_zeros)
    SRKFTUIntermediate(
        Vector{T}(undef, num_x),
        Matrix{T}(undef, num_x, num_x),
        qr_zeros,
        Vector{T}(undef, qr_space_length),
        puft_vcat_q
    )
end

SRKFTUIntermediate(num_x::Number) = SRKFTUIntermediate(Float64, num_x)

function calc_upper_triangular_of_qr!(A)
    R, = LAPACK.gels!('N', A, zeros(size(A, 1), 1)) # Q, R = qr(A)
    R
end

function calc_upper_triangular_of_qr_inplace!(res, A, B, space)
    mygels!(res, A, B, space)
end

# This implementation is based on
# https://github.com/rlabbe/filterpy/blob/master/filterpy/kalman/square_root.py
@inline function calc_apriori_covariance(
    P::Cholesky,
    F::Union{Number, AbstractMatrix},
    Q::Cholesky
)
    A = vcat(P.U * F', Q.U)
    R = calc_upper_triangular_of_qr!(A)
    Cholesky(R, 'U', 0)
end

function time_update!(tu::SRKFTUIntermediate, x, P::Cholesky, F::Union{Number, AbstractMatrix}, Q::Cholesky)
    x_apri = calc_apriori_state!(tu.x_apri, x, F)
    tu.puft_vcat_q[1:size(F, 1),:] .= @~ P.U * F'
    tu.puft_vcat_q[size(F, 1) + 1:end,:] .= Q.U
    R = calc_upper_triangular_of_qr_inplace!(tu.p_apri, tu.puft_vcat_q, tu.qr_zeros, tu.qr_space)
    P_apri = Cholesky(R, 'U', 0)
    KFTimeUpdate(x_apri, P_apri)
end

function measurement_update(
    x,
    P::Cholesky,
    y,
    H::Union{Number, AbstractVector, AbstractMatrix},
    R::Cholesky;
    consider = nothing
)
    ỹ = calc_innovation(H, x, y)
    PHᵀ, S, P_post = calc_cross_cov_innovation_posterior(P, H, R, consider)
    K = calc_kalman_gain(PHᵀ, S.L, consider)
    x_post = calc_posterior_state(x, K, ỹ, consider)
    KFMeasurementUpdate(x_post, P_post, ỹ, S, K)
end

function create_matrix_for_qr(P::Cholesky, H, R::Cholesky)
    num_y = size(R, 1)
    num_x = size(P, 1)
    M = zeros(num_y + num_x, num_y + num_x)
    M[1:num_y, 1:num_y] .= R.U
    M[num_y + 1:end, 1:num_y] .= P.U * H'
    M[num_y + 1:end, num_y + 1:end] .= P.U
    M
end

function calc_cross_cov_innovation_posterior(P::Cholesky, H, R::Cholesky, consider)
    num_y = size(R, 1)
    M = create_matrix_for_qr(P, H, R)
    RU = calc_upper_triangular_of_qr!(M)
    PHᵀ = extract_cross_covariance(RU, num_y)
    S = extract_innovation_covariance(RU, num_y)
    P_post = extract_posterior_covariance(RU, num_y, P, consider)
    PHᵀ, S, P_post
end

@inline extract_cross_covariance(RU, num_y) = transpose(RU[1:num_y, num_y + 1:end])
@inline extract_innovation_covariance(RU, num_y) = Cholesky(RU[1:num_y, 1:num_y], 'U', 0)
@inline extract_posterior_covariance(RU, num_y, P, consider::Nothing) =
    Cholesky(RU[num_y + 1:end, num_y + 1:end], 'U', 0)

struct SRKFMUIntermediate{T,K<:Union{<:AbstractVector{T},<:AbstractMatrix{T}}}
    innovation::Vector{T}
    kalman_gain::K
    m::Matrix{T}
    qr_zeros::Vector{T}
    qr_space::Vector{T}
    R::Matrix{T}
    x_posterior::Vector{T}
end

function SRKFMUIntermediate(T::Type, num_x::Number, num_y::Number)
    qr_zeros = zeros(T, num_x + num_y)
    M = Matrix{T}(undef, num_x + num_y, num_x + num_y)
    qr_space_length = calc_gels_working_size(M, qr_zeros)
    SRKFMUIntermediate(
        Vector{T}(undef, num_y),
        Matrix{T}(undef, num_x, num_y),
        M,
        qr_zeros,
        Vector{T}(undef, qr_space_length),
        Matrix{T}(undef, num_x + num_y, num_x + num_y),
        Vector{T}(undef, num_x)
    )
end

SRKFMUIntermediate(num_x::Number, num_y::Number) = SRKFMUIntermediate(Float64, num_x, num_y)

function measurement_update!(
    mu::SRKFMUIntermediate,
    x,
    P::Cholesky,
    y,
    H::Union{Number, AbstractVector, AbstractMatrix},
    R::Cholesky
)
    ỹ = calc_innovation!(mu.innovation, H, x, y)
    dim_y = length(y)
    M = mu.m
    M[1:dim_y, 1:dim_y] .= R.U
    M[dim_y + 1:end, 1:dim_y] .= @~ P.U * H'
    M[dim_y + 1:end, dim_y + 1:end] .= P.U
    RU = calc_upper_triangular_of_qr_inplace!(mu.R, M, mu.qr_zeros, mu.qr_space)
    PHᵀ = transpose(@view(RU[1:dim_y, dim_y + 1:end]))
    S = Cholesky(@view(RU[1:dim_y, 1:dim_y]), 'U', 0)
    K = calc_kalman_gain!(mu.kalman_gain, PHᵀ, S.L)
    x_post = calc_posterior_state!(mu.x_posterior, x, K, ỹ)
    P_post = Cholesky(@view(RU[dim_y + 1:end, dim_y + 1:end]), 'U', 0)
    KFMeasurementUpdate(x_post, P_post, ỹ, S, K)
end

function calc_kalman_gain!(K, PHᵀ, SL::LowerTriangular)
    K .= PHᵀ
    rdiv!(K, SL)
end


text:
"""
    TraitSubstitutionModel

Abstract type for discrete trait substitution models,
using a continous time Markov model on a phylogeny.
Adapted from the substitutionModels module in BioJulia.
The same Q and P function names are used for the transition rates and probabilities.
"""

abstract type TraitSubstitutionModel end
const SM = TraitSubstitutionModel
const Bmatrix = SMatrix{2, 2, Float64}

"""
    Q(model)

Substitution rate matrix for a given substitution model:
Q[i,j] is the rate of transitioning from state i to state j.
"""
Q(mod::TraitSubstitutionModel) = error("rate matrix Q not defined for $(typeof(mod)).")

"""
    showQ(model)

Print the Q matrix to the screen, with trait states as labels on rows and columns.
adapted from prettyprint function by mcreel, found 2017/10 at
https://discourse.julialang.org/t/display-of-arrays-with-row-and-column-names/1961/6
"""
function showQ(object::TraitSubstitutionModel)
    M = Q(object)
    pad = max(8,maximum(length.(object.label))+1)
    for i = 1:size(M,2) # print the header
        print(lpad(object.label[i],(i==1? 2*pad : pad), " "))
    end
    print("\n")
    for i = 1:size(M,1) # print one row per state
        if object.label != ""
            print(lpad(object.label[i],pad," "))
        end
        for j = 1:size(M,2)
            if j == i
                print(lpad("*",pad," "))
            else
                fmt = "%$(pad).4f"
                @eval(@printf($fmt,$(M[i,j])))
            end
        end
        print("\n")
    end
end

"""
    P(mod, t)

Probability transition matrix for a [`TraitSubstitutionModel`](@ref), of the form

    P[1,1] ... P[1,k]
       .          .
       .          .
    P[k,1] ... P[k,k]

where P[i,j] is the probability of ending in state j after time t,
given that the process started in state i.
"""
@inline function P(mod::SM, t::Float64)
    t >= 0.0 || error("t must be positive")
    return expm(Q(mod) * t)
end

"""
    P(mod, t::Array{Float64})

When applied to a general substitution model, matrix exponentiation is used.
The time argument `t` can be an array.
"""
function P(mod::SM, t::Array{Float64})
    all(t .>= 0.0) || error("t's must all be positive")
    try
        eig_vals, eig_vecs = eig(Q(mod)) # Only hermitian matrices are diagonalizable by
        # *StaticArrays*. Non-Hermitian matrices should be converted to `Array`first.
        return [eig_vecs * expm(diagm(eig_vals)*i) * eig_vecs' for i in t]
    catch
        eig_vals, eig_vecs = eig(Array(Q(mod)))
        k = nStates(mod)
        return [SMatrix{k,k}(eig_vecs * expm(diagm(eig_vals)*i) * inv(eig_vecs)) for i in t]
    end
end

"""
    BinaryTraitSubstitutionModel(α, β [, label])

[`TraitSubstitutionModel`](@ref) for binary traits (with 2 states).
Default labels are "0" and "1".
α is the rate of transition from "0" to "1", and β from "1" to "0".
"""
mutable struct BinaryTraitSubstitutionModel <: TraitSubstitutionModel
    rate::Vector{Float64}
    label::SVector{2, String}
    function BinaryTraitSubstitutionModel(α::Float64, β::Float64, label::SVector{2, String})
	    α >= 0. || error("parameter α must be non-negative")
	    β >= 0. || error("parameter β must be non-negative")
	    ab = α+β
        ab > 0. || error("α+β must be positive")
        new([α,β], label)
    end
end
const BTSM = BinaryTraitSubstitutionModel
BinaryTraitSubstitutionModel(α, β) = BinaryTraitSubstitutionModel(α, β, SVector("0", "1"))

"""
For a BinaryTraitSubstitutionModel, the rate matrix Q is of the form:

    -α  α
     β -β
"""
@inline function Q(mod::BTSM)
    return Bmatrix(-mod.rate[1], mod.rate[2], mod.rate[1], -mod.rate[2])
end

function show(io::IO, object::BinaryTraitSubstitutionModel)
    str = "Binary Trait Substitution Model:\n"
    str *= "rate $(object.label[1])→$(object.label[2]) α=$(object.rate[1])\n"
    str *= "rate $(object.label[2])→$(object.label[1]) β=$(object.rate[2])\n"
    print(io, str)
end

@inline function P(mod::BTSM, t::Float64)
    t >= 0.0 || error("t must be positive")
    ab = mod.rate[1] + mod.rate[2]
    e1 = exp(-ab*t)
    p0 = mod.rate[2]/ab # asymptotic frequency of state "0"
    p1 = mod.rate[1]/ab # asymptotic frequency of state "1"
    a0= p0 *e1
    a1= p1*e1
    return Bmatrix(p0+a1, p0-a0, p1-a1, p1+a0) # by columns
end

"""
    nStates(model)

Number of character states for a given trait evolution model.
"""
nStates(mod::TraitSubstitutionModel) = error("nStates not defined for $(typeof(mod)).")

"""
# Examples

```julia-repl
julia> m1 = BinaryTraitSubstitutionModel(1.0, 2.0)
julia> nStates(m1)
2
```
"""
function nStates(mod::BTSM)
    return 2::Int
end

"""
    TwoBinaryTraitSubstitutionModel(rate [, label])

[`TraitSubstitutionModel`](@ref) for two binary traits, possibly correlated.
Default labels are "x0", "x1" for trait 1, and "y0", "y1" for trait 2.
If provided, `label` should be a vector of size 4, listing labels for
trait 1 first then labels for trait 2.
`rate` should be a vector of substitution rates of size 8.
rate[1],...,rate[4] describe rates of changes in trait 1.
rate[5],...,rate[8] describe rates of changes in trait 2.

In the transition matrix, trait combinations are listed in the following order:
x0-y0, x0-y1, x1-y0, x1-y1.

# example

```julia
model = TwoBinaryTraitSubstitutionModel([2.0,1.2,1.1,2.2,1.0,3.1,2.0,1.1],
        ["carnivory", "noncarnivory", "wet", "dry"]);
model
using PhyloPlots
plot(model) # to visualize states and rates
```
"""
mutable struct TwoBinaryTraitSubstitutionModel <: TraitSubstitutionModel
    rate::Vector{Float64}
    label::Vector{String}
    function TwoBinaryTraitSubstitutionModel(α::AbstractVector{Float64}, label::AbstractVector{String})
	    all( x -> x >= 0., α) || error("rates must be non-negative")
        new(α, [string(label[1], "-", label[3]),
                string(label[1], "-", label[4]),
                string(label[2], "-", label[3]),
                string(label[2], "-", label[4])])
    end
end
const TBTSM = TwoBinaryTraitSubstitutionModel
TwoBinaryTraitSubstitutionModel(α) = TwoBinaryTraitSubstitutionModel(α, ["x0", "x1", "y0", "y1"])

function Q(mod::TwoBinaryTraitSubstitutionModel)
    M = fill(0.0,(4,4))
    a = mod.rate
    M[1,3] = a[1]
    M[3,1] = a[2]
    M[2,4] = a[3]
    M[4,2] = a[4]
    M[1,2] = a[5]
    M[2,1] = a[6]
    M[3,4] = a[7]
    M[4,3] = a[8]
    M[1,1] = -M[1,2] - M[1,3]
    M[2,2] = -M[2,1] - M[2,4]
    M[3,3] = -M[3,4] - M[3,1]
    M[4,4] = -M[4,3] - M[4,2]
    return M
end

function show(io::IO, object::TwoBinaryTraitSubstitutionModel)
    print("Substitution model for 2 binary traits, with rate matrix:\n")
    showQ(object)
end

"""
    EqualRatesSubstitutionModel(numberStates, α, labels)

[`TraitSubstitutionModel`](@ref) for traits with any number of states
and equal substitution rates α between all states.
Default labels are "1","2",...
"""
mutable struct EqualRatesSubstitutionModel <: TraitSubstitutionModel
    k::Int
    α::Float64
    label::Vector{String}
    function EqualRatesSubstitutionModel(k::Int, α::Float64, label::Vector{String})
        k >= 2 || error("parameter k must be greater than or equal to 2")
        α > 0 || error("parameter α must be positive")
        length(label)==k || error("label vector of incorrect length")
        new(k, α, label)
    end
end
EqualRatesSubstitutionModel(k, α) = EqualRatesSubstitutionModel(k, α, string.(1:k))

function show(io::IO, object::EqualRatesSubstitutionModel)
    str = "Equal Rates Substitution Model with k=$(object.k),\n"
    str *= "all rates equal to α=$(object.α).\n"
    str *= "rate matrix Q:\n"
    print(io, str)
    showQ(object)
end

function nStates(mod::EqualRatesSubstitutionModel)
    return mod.k
end

function Q(mod::EqualRatesSubstitutionModel)
    M = fill(mod.α,(mod.k,mod.k))
    d = -(mod.k-1)*mod.α
    for i in 1:mod.k
        M[i,i] = d
    end
    return M
end

"""
    randomTrait(model, t, start)
    randomTrait!(end, model, t, start)

Simulate traits along one edge of length t.
`start` must be a vector of integers, each representing the starting value of one trait.
The bang version (ending with !) uses the vector `end` to store the simulated values.

# Examples
```julia-repl
julia> m1 = PhyloNetworks.BinaryTraitSubstitutionModel(1.0, 2.0)
julia> srand(12345);
julia> randomTrait(m1, 0.2, [1,2,1,2,2])
 5-element Array{Int64,1}:
 1
 2
 1
 2
 2
```
"""
function randomTrait(mod::SM, t::Float64, start::AbstractVector{Int})
    res = Vector{Int}(length(start))
    randomTrait!(res, mod, t, start)
end

function randomTrait!(endTrait::AbstractVector{Int}, mod::SM, t::Float64, start::AbstractVector{Int})
    Pt = P(mod, t)
    k = size(Pt, 1) # number of states
    w = [aweights(Pt[i,:]) for i in 1:k]
    for i in 1:length(start)
        endTrait[i] =sample(1:k, w[start[i]])
    end
    return endTrait
end

"""
    randomTrait(mod, net; ntraits=1, keepInternal=true, checkPreorder=true)

Simulate evolution of discrete traits on a rooted evolutionary network based on
the supplied evolutionary model. Trait sampling is uniform at the root.

# Arguments
- ntraits: number of traits to be simulated (default: 1 trait).
- keepInternal: if true, export character states at all nodes, including
  internal nodes. if false, export character states at tips only.

# Output
- matrix of character states with one row per trait, one column per node
- array of node labels (for tips) or node numbers (for internal nodes)
  in the same order as columns in the character state matrix

# Examples
```julia-repl
julia> m1 = PhyloNetworks.BinaryTraitSubstitutionModel(1.0, 2.0)
julia> net = readTopology("(A:1.0,(B:1.0,(C:1.0,D:1.0):1.0):1.0);")
julia> srand(12345);
julia> a,b = randomTrait(m1, net)
 ([1 2 … 1 2], String["-2", "-3", "-4", "D", "C", "B", "A"])
julia> a
 1×7 Array{Int64,2}:
 1  2  1  1  1  1  2
julia> b
 7-element Array{String,1}:
 "-2"
 "-3"
 "-4"
 "D"
 "C"
 "B"
 "A"
```
"""

function randomTrait(mod::SM, net::HybridNetwork;
    ntraits=1::Int, keepInternal=true::Bool, checkPreorder=true::Bool)
    net.isRooted || error("net needs to be rooted for preorder recursion")
    if(checkPreorder)
        preorder!(net)
    end
    nnodes = net.numNodes
    M = Matrix{Int}(ntraits, nnodes) # M[i,j]= trait i for node j
    randomTrait!(M,mod,net)
    if !keepInternal
        M = getTipSubmatrix(M, net, indexation=:cols) # subset columns only. rows=traits
        nodeLabels = [n.name for n in net.nodes_changed if n.leaf]
    else
        nodeLabels = [n.name == "" ? string(n.number) : n.name for n in net.nodes_changed]
    end
    return M, nodeLabels
end

function randomTrait!(M::Matrix{Int}, mod::SM, net::HybridNetwork)
    recursionPreOrder!(net.nodes_changed, M, # updates M in place
            updateRootRandomTrait!,
            updateTreeRandomTrait!,
            updateHybridRandomTrait!,
            mod)
end

function updateRootRandomTrait!(V::AbstractArray, i::Int, mod)
    sample!(1:nStates(mod), view(V, :, i)) # uniform at the root
    return
end

function updateTreeRandomTrait!(V::Matrix,
    i::Int,parentIndex::Int,edge::Edge,
    mod)
    randomTrait!(view(V, :, i), mod, edge.length, view(V, :, parentIndex))
end

function updateHybridRandomTrait!(V::Matrix,
        i::Int, parentIndex1::Int, parentIndex2::Int,
        edge1::Edge, edge2::Edge, mod)
    randomTrait!(view(V, :, i), mod, edge1.length, view(V, :, parentIndex1))
    tmp = randomTrait(mod, edge2.length, view(V, :, parentIndex2))
    for j in 1:size(V,1) # loop over traits
        if V[j,i] == tmp[j] # both parents of the hybrid node have the same trait
            continue # skip the rest: go to next trait
        end
        if rand() > edge1.gamma
            V[j,i] = tmp[j] # switch to inherit trait of parent 2
        end
    end
end


text:
function nlp3(oct::Bool = false)
    m = JuMP.Model()
    @variable(m, 0 <= x[1:10])
    for i = 1:8
        JuMP.set_lower_bound(x[[1, 2, 4, 6, 7, 8, 9, 10][i]],
                             [1,1,1,85,90,3,1.2,145][i])
    end
    for i = 1:10
        JuMP.set_upper_bound(x[i], [2000, 16000, 120, 5000, 2000,
                                    93, 95, 12, 4, 162][i])
        JuMP.set_start_value(x[i], [1724.90452208, 16000,
                                    98.0900813608, 3049.1211364,
                                    1995.02326433, 90.718089075,
                                    94.2274481766, 10.432474977,
                                    2.59051951438, 149.682344530][i])
    end
    @constraint(m, e1, x[1] - 1.22*x[4] + x[5] == 0)
    @constraint(m, e2, x[9] + 0.222*x[10] == 35.82)
    @constraint(m, e3, 3*x[7] - x[10] == 133)
    if !oct
        @NLconstraint(m, e4, x[7] - 1.098*x[8] + 0.038*(x[8]^2) - 0.325*(x[6] - 89) == 86.35)
        @NLconstraint(m, e5, x[4]*x[9]*x[6] + 1000*x[3]*x[6] - 98000*x[3] == 0)
        @NLconstraint(m, e6, x[2] + x[5] - x[1]*x[8] == 0)
        @NLconstraint(m, e7, 1.12*x[1] + 0.13167*x[8]*x[1] - 0.00667*(x[8]^2)*x[1] - x[4] >= 0)
        @NLobjective(m, Min, 5.04*x[1] + 0.035*x[2] + 10*x[3] + 3.36*x[5] - 0.063*x[4]*x[7])
        return m
    else
        @variable(m, obj)
        @objective(m, Min, obj)
        gm = GlobalModel(model = m, name = "nlp3")
        add_nonlinear_constraint(gm, :(x -> x[7] - 1.098*x[8] + 0.038*(x[8]^2) - 0.325*(x[6] - 89) - 86.35),
                                vars = [x[6], x[7], x[8]], name = "e4", equality=true)
        add_nonlinear_constraint(gm, :(x -> x[4]*x[9]*x[6] + 1000*x[3]*x[6] - 98000*x[3]),
                                vars = [x[3], x[4], x[6], x[9]], name = "e5", equality = true)
        add_nonlinear_constraint(gm, :(x -> x[2] + x[5] - x[1]*x[8]),
                                vars = [x[1], x[2], x[5], x[8]], name = "e6", equality = true)
        add_nonlinear_constraint(gm, :(x -> 1.12*x[1] + 0.13167*x[8]*x[1] - 0.00667*(x[8]^2)*x[1] - x[4]),
                                vars = [x[1], x[4], x[8]], name = "e7")
        add_nonlinear_constraint(gm, :(x -> 5.04*x[1] + 0.035*x[2] + 10*x[3] + 3.36*x[5] - 0.063*x[4]*x[7]),
                                vars = [x[1], x[2], x[3], x[4], x[5], x[7]],
                                dependent_var = obj, name = "obj")
        return gm
    end
end


text:
module YAXArrayBase
using Requires: @require
using DataStructures: OrderedDict

include("datasets/datasetinterface.jl")
include("axisarrays/axisinterface.jl")
include("axisarrays/namedtuple.jl")

defaultfillval(T::Type{<:AbstractFloat}) = convert(T,1e32)
defaultfillval(::Type{Float16}) = Float16(3.2e4)
defaultfillval(T::Type{<:Integer}) = typemax(T)
defaultfillval(T::Type{<:AbstractString}) = ""

function __init__()
  """
  YAXArrayBase.backendlist

  List of symbols defining dataset backends.
  """
  backendlist = OrderedDict{Symbol, Any}(
    :array => Array,
  )

  backendregex = Pair[]

  @require NamedDims="356022a1-0364-5f58-8944-0da4b18d706f" include("axisarrays/nameddims.jl")

  @require DimensionalData="0703355e-b756-11e9-17c0-8b28908087d0" include("axisarrays/dimensionaldata.jl")

  @require AxisArrays="39de3d68-74b9-583c-8d2d-e117c070f3a9" include("axisarrays/axisarrays.jl")

  @require AxisIndices="f52c9ee2-1b1c-4fd8-8546-6350938c7f11" include("axisarrays/axisindices.jl")

  @require ArchGDAL="c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3" begin
    include("axisarrays/archgdal.jl")
    include("datasets/archgdal.jl")
  end

  @require Zarr="0a941bbe-ad1d-11e8-39d9-ab76183a1d99" include("datasets/zarr.jl")

  @require NetCDF="30363a11-5582-574a-97bb-aa9a979735b9" include("datasets/netcdf.jl")

end

export dimvals, dimname, dimnames, iscontdim, iscompressed,
getattributes, getdata, yaxconvert
export get_var_handle, get_varnames, get_var_dims, get_var_attrs,
create_empty, add_var, allow_parallel_write,
to_dataset, allow_missings
end # module


text:
###### n = 2 ##########
puzz21 = [
    3 0 4 0
    0 1 0 2
    0 4 0 3
    2 0 1 0
];

puzz22 = [
    0 1 3 0
    0 0 0 0
    0 0 0 0
    3 4 2 1
];

puzz23 = [
    0 3 0 1
    0 0 2 0
    0 0 3 0
    0 4 0 2
];

###### n = 3 ###########
puzz1 = [
    9 3 0 0 0 0 0 4 0
    0 0 0 0 4 2 0 9 0
    8 0 0 1 9 6 7 0 0
    0 0 0 4 7 0 0 0 0
    0 2 0 0 0 0 0 6 0
    0 0 0 0 2 3 0 0 0
    0 0 8 5 3 1 0 0 2
    0 9 0 2 8 0 0 0 0
    0 7 0 0 0 0 0 5 3
];

puzz2 = [
    0 0 9 0 1 0 0 6 0
    0 0 0 0 6 7 0 0 3
    0 3 5 0 0 0 0 0 0
    3 0 0 0 0 2 0 0 0
    5 1 7 0 0 0 2 3 4
    0 0 0 3 0 0 0 0 7
    0 0 0 0 0 0 9 1 0
    1 0 0 6 8 0 0 0 0
    0 7 0 0 3 0 5 0 0
];

puzz3 = [
    0 3 7 8 6 0 0 4 0
    0 0 6 0 0 7 0 0 0
    0 2 0 0 3 0 0 0 6
    0 8 0 2 0 0 0 0 0
    9 0 0 0 0 0 0 0 1
    0 0 0 0 0 6 0 9 0
    5 0 0 0 7 0 0 3 0
    0 0 0 3 0 0 8 0 0
    0 1 0 0 8 4 2 6 0
];

#### n = 4 #####
puzz41 = [
    0 16 0 0 0 0 0 0 0 0 0 0 0 0 0 1
    10 0 0 14 2 0 5 7 15 0 0 11 0 0 13 9
    0 0 0 11 0 0 0 0 13 0 0 12 0 0 0 0
    0 0 3 0 0 10 14 0 0 0 0 0 0 0 5 0
    0 0 0 7 0 5 0 15 0 6 13 10 14 0 0 3
    9 11 0 0 3 0 12 1 0 14 0 16 7 6 0 8
    4 0 0 0 0 6 0 0 12 0 1 0 0 0 0 15
    0 0 12 0 0 0 0 11 0 0 0 0 0 0 0 16
    3 0 0 9 1 0 11 0 0 0 0 0 16 10 0 0
    0 0 0 0 7 4 0 9 0 2 0 5 0 0 0 11
    0 1 0 0 0 16 2 0 0 0 8 0 15 0 7 0
    0 0 0 0 0 0 0 8 3 0 9 0 4 0 14 0
    0 0 0 0 15 0 0 0 2 0 0 13 0 0 1 0
    1 0 0 5 0 8 0 0 0 10 0 14 11 0 0 12
    16 15 0 8 0 0 0 14 0 1 12 4 10 0 3 0
    0 0 11 0 5 0 0 0 0 0 0 15 0 8 9 2
];


text:
using BinaryBuilder

cuda_version = v"10.1.243"   # NOTE: could be less specific

sources_linux_x64 = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.0.2.39/10.1_20200626/cudnn-10.1-linux-x64-v8.0.2.39.tgz",
                  "82148a68bd6bdaab93af5e05bb1842b8ccb3ab7de7bed41f609a7616c102213d")
]
sources_linux_ppc64le = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.0.2.39/10.1_20200626/cudnn-10.1-linux-ppc64le-v8.0.2.39.tgz",
                  "8196ec4f031356317baeccefbc4f61c8fccb2cf0bdef0a6431438918ddf68fb9")
]
sources_windows = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.0.2.39/10.1_20200626/cudnn-10.1-windows10-x64-v8.0.2.39.zip",
                  "4922e878681502c7135e1b9c7702d6d9b05b5c0d5a175bb91d878bf076657948")
]

include("../common.jl")

max_stars_repo_name: jpsamaroo/Yggdrasil

max_issues_repo_name: jpsamaroo/Yggdrasil

max_forks_repo_name: jpsamaroo/Yggdrasil


text:
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = page_model()
    include("../src/components/CO2forcing.jl")

    add_comp!(m, co2forcing)

    set_param!(m, :co2forcing, :c_CO2concentration, readpagedata(m,"test/validationdata/$valdir/c_co2concentration.csv"))

    ##running Model
    run(m)

    forcing=m[:co2forcing,:f_CO2forcing]
    forcing_compare=readpagedata(m,"test/validationdata/$valdir/f_co2forcing.csv")

    @test forcing ≈ forcing_compare rtol=1e-3
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
module Backend

using LibPQ
using Tables

export conn, init_db, get_conn, update_contact_list, get_key, set_key, get_state, set_state

conn = -1

function update_contact_list(message)
    id = message["message"]["from"]["id"]
    res = execute(conn, "SELECT * from pizza_bot WHERE id=$id;")
    if length(columntable(res)[1]) == 0
       res = execute(conn, "INSERT INTO pizza_bot VALUES ($id,0);")
    end
end

function send_message(params)
    req = HTTP.request("POST",string(url,key,"/sendMessage"),["Content-Type" => "application/json"],JSON.json(params))
end

function get_key(id, key)
    res = execute(conn, "SELECT * from pizza_bot WHERE id=$id;")
    data = columntable(res)
    keys = Dict("id"=>data[1][1], "state"=>data[2][1], "program"=>data[3][1])
    return keys[key]
end

function set_key(id, key, value)
    res = execute(conn, "UPDATE pizza_bot SET $(key)='$value' WHERE id=$id;")
end

function get_state(id)
    return get_key(id,"state")
end

function set_state(id, state)
    set_key(id,"state", state)
end

function init_db()
    conn = LibPQ.Connection("dbname=pizza_bot user=bot password=pizza host=database")

    result = execute(conn, """
        CREATE TABLE IF NOT EXISTS pizza_bot (
            id    integer PRIMARY KEY,
            state   integer,
            program varchar(10)
        );
    """)
    println(columntable(result))

    close(conn)
end

function get_conn()
    global conn
    if conn == -1
        conn = LibPQ.Connection("dbname=pizza_bot user=bot password=pizza host=database")
    end
    return conn
end

end


text:
@testset "contrasts" begin

    d = DataFrame(y = rand(6),
                  x = [:b, :a, :c, :a, :a, :b])

    mf = ModelFrame(@formula(y ~ x), d)

    ## testing equality of ContrastsMatrix
    should_equal = [ContrastsMatrix(DummyCoding(), [:a, :b, :c]),
                    ContrastsMatrix(DummyCoding(base=:a), [:a, :b, :c]),
                    ContrastsMatrix(DummyCoding(base=:a, levels=[:a, :b, :c]), [:a, :b, :c]),
                    ContrastsMatrix(DummyCoding(levels=[:a, :b, :c]), [:a, :b, :c])]

    for c in should_equal
        @test mf.schema[term(:x)].contrasts == c
        @test hash(mf.schema[term(:x)].contrasts) == hash(c)
    end

    should_not_equal = [ContrastsMatrix(EffectsCoding(), [:a, :b, :c]),
                        ContrastsMatrix(DummyCoding(), [:b, :c]),
                        ContrastsMatrix(DummyCoding(), [:b, :c, :a]),
                        ContrastsMatrix(DummyCoding(), [:b, :c, :a]),
                        ContrastsMatrix(DummyCoding(), [:b, :c, :a]),
                        ContrastsMatrix(DummyCoding(base=:b), [:a, :b, :c]),
                        ContrastsMatrix(DummyCoding(base=:a, levels=[:b, :a, :c]), [:a, :b, :c])]

    for c in should_not_equal
        @test mf.schema[term(:x)].contrasts != c
    end

    # Dummy coded contrasts by default:
    @test ModelMatrix(mf).m == [1  1  0
                                1  0  0
                                1  0  1
                                1  0  0
                                1  0  0
                                1  1  0]
    @test coefnames(mf) == ["(Intercept)"; "x: b"; "x: c"]

    mmm = ModelMatrix(mf).m
    setcontrasts!(mf, x = DummyCoding())
    @test ModelMatrix(mf).m == mmm

    setcontrasts!(mf, x = EffectsCoding())
    @test ModelMatrix(mf).m == [1  1  0
                                1 -1 -1
                                1  0  1
                                1 -1 -1
                                1 -1 -1
                                1  1  0]
    @test coefnames(mf) == ["(Intercept)"; "x: b"; "x: c"]

    # change base level of contrast
    setcontrasts!(mf, x = EffectsCoding(base = :b))
    @test ModelMatrix(mf).m == [1 -1 -1
                                1  1  0
                                1  0  1
                                1  1  0
                                1  1  0
                                1 -1 -1]
    @test coefnames(mf) == ["(Intercept)"; "x: a"; "x: c"]

    # change levels of contrast
    setcontrasts!(mf, x = EffectsCoding(levels = [:c, :b, :a]))
    @test ModelMatrix(mf).m == [1  1  0
                                1  0  1
                                1 -1 -1
                                1  0  1
                                1  0  1
                                1  1  0]
    @test coefnames(mf) == ["(Intercept)"; "x: b"; "x: a"]

    # change levels and base level of contrast
    setcontrasts!(mf, x = EffectsCoding(levels = [:c, :b, :a], base = :a))
    @test ModelMatrix(mf).m == [1  0  1
                                1 -1 -1
                                1  1  0
                                1 -1 -1
                                1 -1 -1
                                1  0  1]
    @test coefnames(mf) == ["(Intercept)"; "x: c"; "x: b"]

    # respect order of levels

    data = DataFrame(x = levels!(categorical(['A', 'B', 'C', 'C', 'D']), ['C', 'B', 'A', 'D']))
    f = apply_schema(@formula(x ~ 1), schema(data))
    @test modelcols(f.lhs, data) == [0 1 0; 1 0 0; 0 0 0; 0 0 0; 0 0 1]

    # Helmert coded contrasts
    setcontrasts!(mf, x = HelmertCoding())
    @test ModelMatrix(mf).m == [1  1 -1
                                1 -1 -1
                                1  0  2
                                1 -1 -1
                                1 -1 -1
                                1  1 -1]
    @test coefnames(mf) == ["(Intercept)"; "x: b"; "x: c"]

    # Mismatching types of data and contrasts levels throws an error:
    @test_throws ArgumentError setcontrasts!(mf, x = EffectsCoding(levels = ["a", "b", "c"]))

    # Missing data is handled gracefully, dropping columns when a level is lost
    allowmissing!(d,:x)
    d[3, :x] = missing
    mf_missing = ModelFrame(@formula(y ~ x), d,
                            contrasts = Dict(:x => EffectsCoding()))
    @test ModelMatrix(mf_missing).m == [1  1
                                        1 -1
                                        1 -1
                                        1 -1
                                        1  1]
    @test coefnames(mf_missing) == ["(Intercept)"; "x: b"]

    # Things that are bad to do:
    # Applying contrasts that only have a subset of data levels:
    @test_throws ArgumentError setcontrasts!(mf, x = EffectsCoding(levels = [:a, :b]))
    # Applying contrasts that expect levels not found in data:
    @test_throws ArgumentError setcontrasts!(mf, x = EffectsCoding(levels = [:a, :b, :c, :d]))
    # Asking for base level that's not found in data
    @test_throws ArgumentError setcontrasts!(mf, x = EffectsCoding(base = :e))

    # Manually specified contrasts
    contrasts = [0  1
                 -1 -.5
                 1  -.5]
    setcontrasts!(mf, x = ContrastsCoding(contrasts))
    @test ModelMatrix(mf).m == [1 -1 -.5
                                1  0  1
                                1  1 -.5
                                1  0  1
                                1  0  1
                                1 -1 -.5]

    # throw argument error if number of levels mismatches
    @test_throws ArgumentError setcontrasts!(mf, x = ContrastsCoding(contrasts[1:2, :]))
    @test_throws ArgumentError setcontrasts!(mf, x = ContrastsCoding(hcat(contrasts, contrasts)))

    # contrasts types must be instantiated (should throw ArgumentError, currently
    # MethodError on apply_schema)
    @test_broken setcontrasts!(mf, x = DummyCoding)

end


text:
module Rayly

using LinearAlgebra
using StaticArrays
using Statistics
using ColorTypes
using ProgressMeter
using FileIO
import GeometryTypes: HomogenousMesh
using Images
using Random

for (n,f) in [(:compmin, :min), (:compmax, :max)]
    @eval $n(a::SVector{3,T}, b::SVector{3,T}) where T =
        SVector{3,T}($f(a[1], b[1]), $f(a[2], b[2]), $f(a[3], b[3]))
end

struct Ray{T<:AbstractFloat}
    pos::SVector{3,T}
    dir::SVector{3,T}
    Ray(pos::SVector{3,T}, dir::SVector{3,T}) where T =
        new{T}(pos, normalize(dir))
end

abstract type Intersectable{T<:AbstractFloat} end
Base.eltype(::Intersectable{T}) where T = T

abstract type Accelerator{T<:AbstractFloat,O<:Intersectable{T}} end

include("intersection.jl")

include("aabb.jl")
include("sphere.jl")
include("triangle.jl")
include("scene.jl")

include("camera.jl")
include("list_accelerator.jl")
include("bvh.jl")
include("bvh_simple_tree.jl")
include("bvh_int_tree.jl")
include("bvh_simple_builder.jl")
include("bvh_accelerator.jl")
include("bvh_afra.jl")
include("bvh_debug.jl")
include("sampling.jl")
include("render.jl")
include("utils.jl")

# Default BVH builder is bvh_simple_build
Base.convert(::Type{Tree}, scene::Scene{O}) where {T,Tree<:AbstractTree{T},O<:Intersectable{T}} =
    bvh_simple_build(Tree, scene.objs)

export Ray, Intersectable, Accelerator

function __init__()
    add_format(format"RSC", "RSC", ".rsc", [:Rayly])
end

end # module


text:
"""
    ranks(x; groups = 10, rank = tiedrank, rev = false)

Return a variable which bins `x` into `groups` number of bins.
The `rank` keyword allows different ranking method;
use `rev = true` to reverse sort so that small bin number is large value of `x`.
Missing values are assigned to group `missing`.

Default values of `rank = tiedrank` and `rev = false` results
in similar grouping as SAS PROC RANK groups=n tied=mean.
"""
function ranks(x::Vector{T} where {T<:Real}; groups::Int = 10, rank = tiedrank, rev = false)::Vector{Int32}
    nx = length(x)
    nx >= 2 || throw(ArgumentError("x length should be >= 2"))
    groups >= 2 || throw(ArgumentError("groups should be >= 2"))

    floor.(Int32, rank(x, rev = rev) * groups / (nx+1))
end

function ranks(x::Vector{T} where {T<:Union{Missing,Real}}; groups::Int = 10, rank = tiedrank, rev = false)::Vector{Union{Missing,Int32}}
    nx = length(x)
    nx >= 2 || throw(ArgumentError("x length should be >= 2"))
    groups >= 2 || throw(ArgumentError("groups should be >= 2"))

    floor.(Union{Missing,Int32}, rank(x, rev = rev) * groups / (nx+1))
end

"""
    strstd(s; empty = "")

Standardizes string `s` by stripping leading and trailing blanks,
embedded blanks are replaced with underscore _.
If `s` is missing or all blanks, the result is value of `empty`.

Default of `empty` is "", missing values and all blank strings become "".
"""
function strstd(s::AbstractString; empty = "")
    str = strip(s)                          # strip leading and trailing blanks
    if isempty(str)                         # replace empty string
        str = empty                         # empty is not subject to blank to underscore
    else
        str = replace(str, ' ' => '_')        # embedded blanks become underscore
    end

    return str
end

function strstd(s::Missing; empty = "")
    return empty
end

"""
    sumxm(x...)     sum treating missings as 0

Returns the sum of `x` treating missing values as 0's.
`x` can be varying number of scalars, in this case, their sum is returned.
If `x` is varying number of vectors, they are summed element-wise across the vectors.
"""
sumxm(x...) = _sumxm.(x...)                         # . over "rows" across varying vectors
_sumxm(x::T...) where T<:Real = sum(x)              # no missing
_sumxm(x...) = sum(z -> ismissing(z) ? 0 : z, x)    # missing is 0 in sum

"""
    kstest(class, var; rev = true)

Calculate empirical 2 sample Kolmogorov-Smirnov statistic and its location.
`class` is a 2 level categorical variable, `var` is the distribution to analyze.

Returns:
- n, total number of observations
- n1, number of observations of class 1
- n0, number of observations of class 0
- baserate, incidence rate of class 1
- ks, the maximum separation between the two cumulative distributions
- ksarg, the value of `var` at which maximum separation is achieved
- ksdep, depth of ksarg in the sorted values of `var`;
`rev` = true counts depth from high value towards low value.
"""
function kstest(class::BitArray{1}, var::Vector; rev = true)
    n = length(class)
    n == length(var) || error(ArgumentError("class and var should have the same length"))
    n1 = sum(class)
    n1 == 0 && error(ArgumentError("there are no class 1"))
    n0 = n - n1
    n0 == 0 && error(ArgumentError("there are no class 0"))

    baserate = n1 / n
    idx = sortperm(var, rev = rev)
    tgt = class[idx]
    cdf1 = cumsum(tgt) ./ n1            # TPR: True Positive Rate
    cdf0 = cumsum(.!tgt) ./ n0          # FPR: False Positive Rate
    sep = cdf1 .- cdf0
    ks, ksidx = findmax(sep)
    ksarg = var[idx[ksidx]]             # var @ max sep
    ksdep = ksidx / n                   # depth @ max sep

    (n = n, n1 = n1, n0 = n0, baserate = baserate, ks = ks, ksarg = ksarg, ksdep = ksdep)
end

function kstest(class::Vector, var::Vector; rev = true)
    uc = unique(class)
    length(uc) == 2 || error(ArgumentError("class should have 2 levels"))

    kstest(class .== uc[2], var; rev = rev)
end

"""
    auroc(class, var; tie = 1e-6)

Calculate area under Receiver Operating Characteristics (ROC) curve,
`class` is a 2 level categorical variable, `var` is the distribution to analyze.
Pair-wise comparison between class 1 values with class 0 values are made as follows:
- class 1 value > class 0 value is Concordant
- class 1 value ≈ class 0 value (within `tie`) is Tied
- class 1 value < class 0 value is Discordant

Returns:
- concordant, number of concordant comparisons
- tied, number of tied comparisons
- discordant, number of discordant comparisons
- auc, or C, is (Concordant + 0.5Tied) / Total comparisons; same as numeric integration of ROC curve
- gini, 2C-1, also known as Somer's D, is (Concordant - Discordant) / Total comparisons

Note there are other rank order measures:
- Goodman-Kruskal Gamma is (Concordant - Discordant) / (Concordant + Discordant), no penalty for ties
- Kendall's Tau is (Concordant - Discordant) / (0.5 N(N-1))
"""
function auroc(class::BitArray{1}, var::Vector; tie = 1e-6)
    n = length(class)
    n == length(var) || error(ArgumentError("class and var should have the same length"))
    n1 = sum(class)
    n1 == 0 && error(ArgumentError("there are no class 1"))
    n0 = n - n1
    n0 == 0 && error(ArgumentError("there are no class 0"))

    c1 = sort!(var[class])              # sorted var of class 1
    c0 = sort!(var[.!class])            # sorted var of class 0

    conc = 0
    tied = 0
    l_ix = 1
    u_ix = 1
    @inbounds for v in c1               # loop over this usually smaller array
        l_v = v - tie
        u_v = v + tie

        while l_ix <= n0                # first index within [l_v, u_v] window
            l_v <= c0[l_ix] && break
            l_ix += 1                   # could exit as n0 + 1
        end

        while u_ix <= n0                # first index beyond [l_v, u_v] window
            u_v < c0[u_ix] && break
            u_ix += 1                   # could exit as n0 + 1
        end

        conc += (l_ix - 1)              # lower than l_ix are Concordant
        tied += (u_ix - l_ix)           # within window are tied
    end

    pairs = n1 * n0
    disc = pairs - conc - tied
    auc = (conc + 0.5tied) / pairs

    (conc = conc, tied = tied, disc = disc, auc = auc, gini = 2auc - 1)
end

function auroc(class::Vector, var::Vector; tie = 1e-6)
    uc = unique(class)
    length(uc) == 2 || error(ArgumentError("class should have 2 levels"))

    auroc(class .== uc[2], var; tie = tie)
end

"""
    concordance(class, var, tie)

Concordance calculation with flexible tied region (`auroc` uses fixed width region).
`class` is a 2 level categorical variable, `var` is the distribution to analyze,
`tie`(x) returns the lower and upper bound of tied region of x.

Pair-wise comparison between class 1 values with class 0 values are made as follows:
class 1 value > class 0 value is Concordant; class 1 value ≈ class 0 value (within `tie`) is Tied;
class 1 value < class 0 value is Discordant.

Returns:
- concordant, number of concordant comparisons
- tied, number of tied comparisons
- discordant, number of discordant comparisons
- auroc, or C, is (Concordant + 0.5Tied) / Total comparisons; same as numeric integration of ROC curve
- gini, 2C-1, also known as Somer's D, is (Concordant - Discordant) / Total comparisons

Note Goodman-Kruskal Gamma is (Concordant - Discordant) / (Concordant + Discordant);
and Kendall's Tau is (Concordant - Discordant) / (0.5 x Total count x (Total count - 1))
"""
function concordance(class::BitArray{1}, var::Vector, tie)
    n = length(class)
    n == length(var) || error(ArgumentError("class and var should have the same length"))
    n1 = sum(class)
    n1 == 0 && error(ArgumentError("there are no class 1"))
    n0 = n - n1
    n0 == 0 && error(ArgumentError("there are no class 0"))

    c1 = sort!(var[class])              # sorted var of class 1
    c0 = sort!(var[.!class])            # sorted var of class 0

    conc = 0
    tied = 0
    l_ix = 1
    u_ix = 1
    @inbounds for v in c1               # loop over this usually smaller array
        l_v, u_v = tie(v)               # <== tie function

        while l_ix <= n0                # first index within [l_v, u_v] window
            l_v <= c0[l_ix] && break
            l_ix += 1                   # could exit as n0 + 1
        end

        while u_ix <= n0                # first index beyond [l_v, u_v] window
            u_v < c0[u_ix] && break
            u_ix += 1                   # could exit as n0 + 1
        end

        conc += (l_ix - 1)              # lower than l_ix are Concordant
        tied += (u_ix - l_ix)           # within window are tied
    end

    pairs = n1 * n0
    disc = pairs - conc - tied
    auc = (conc + 0.5tied) / pairs

    (conc = conc, tied = tied, disc = disc, auc = auc, gini = 2auc - 1)
end

function concordance(class::Vector, var::Vector, tie)
    uc = unique(class)
    length(uc) == 2 || error(ArgumentError("class should have 2 levels"))

    concordance(class .== uc[2], var, tie)
end

"""
    infovalue(g::Vector{Integer}, b::Vector{Integer}))

Information value calculation of `g`, `b` vector of binned frequency counts

- weight of evidence = log(density g / density b), 0 adjusted
- infovalue = sum (density g - density b) * weight of evidence

Industry rule of thumb:
- iv <= 0.1         no significant change
- 0.1 < iv <= 0.25  minor change
- 0.25 < iv         major change
"""
function infovalue(g::Vector{TG} where TG <: Integer,
                   b::Vector{TB} where TB <: Integer; woeadjust = 0.5)
    length(g) >= 2 || error("frequency vectors should have at least 2 elements")
    length(g) == length(b) || error("frequency vectors should be the same length")

    # infovalue is class invariant - swapping g/b doesn't matter
    # infovalue is order invariant - re-ordering rows doesn't matter
    # implementation modelled after SAS PROC BINNING where
    # weight of evidence is adjusted to avoid divide by 0

    pb = b ./ sum(b)
    pg = g ./ sum(g)
    pba = (b + (b .== 0) * woeadjust) ./ sum(b)     # adjust only pure bins
    pga = (g + (g .== 0) * woeadjust) ./ sum(g)
    woe = log.(pga ./ pba)
    sum((pg .- pb) .* woe)
end


text:
using Test
import Weave: Markdown2HTML
import Markdown

# Test markdown2html writer

html = Markdown2HTML.html(Markdown.parse("""

# H1

## H2

## H3

Some **text** with different [^note] *formatting* and \$math\$ and text.

`some code` with [link](http://github.com)

[^note]: test note

---

\$more math\$

* List one
* List two

1. List one
2. List two

```julia
x = 3
```

!!! note "Something"

    Test admonition with ![Image](link/to/image.png)

> Some important quote

"""))

ref_html = """
H1

H2

H3

Some text with different [note] formatting and \$math\$ and text.

some code with link

note

test note

\\[
\\begin{align}
more math
\\end{align}
\\]

List one

List two

List one

List two

x = 3

Something

Test admonition with

Some important quote

"""

@test html == ref_html


text:
module CommunalHelperDreamMoveBlock

using ..Ahorn, Maple

@mapdef Entity "CommunalHelper/DreamMoveBlock" DreamMoveBlock(x::Integer, y::Integer,
	width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight,
	direction::String="Right", fast::Bool=false, noCollide::Bool=false,
	featherMode::Bool=false, oneUse::Bool=false, doubleRefill::Bool=false)

const placements = Ahorn.PlacementDict(
    "Dream Move Block ($direction) (Communal Helper)" => Ahorn.EntityPlacement(
        DreamMoveBlock,
        "rectangle",
        Dict{String, Any}(
            "direction" => direction
        )
    ) for direction in Maple.move_block_directions
)

Ahorn.editingOptions(entity::DreamMoveBlock) = Dict{String, Any}(
    "direction" => Maple.move_block_directions
)
Ahorn.minimumSize(entity::DreamMoveBlock) = 16, 16
Ahorn.resizable(entity::DreamMoveBlock) = true, true

Ahorn.selection(entity::DreamMoveBlock) = Ahorn.getEntityRectangle(entity)

const arrows = Dict{String, String}(
    "up" => "objects/CommunalHelper/dreamMoveBlock/arrow02",
    "left" => "objects/CommunalHelper/dreamMoveBlock/arrow04",
    "right" => "objects/CommunalHelper/dreamMoveBlock/arrow00",
    "down" => "objects/CommunalHelper/dreamMoveBlock/arrow06",
)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::DreamMoveBlock, room::Maple.Room)
    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

	 Ahorn.Cairo.save(ctx)

    Ahorn.set_antialias(ctx, 1)
    Ahorn.set_line_width(ctx, 1)

    fillColor = get(entity.data, "featherMode", false) ? (0.31, 0.69, 1.0, 0.4) : (0.0, 0.0, 0.0, 0.4)
	 lineColor = get(entity.data, "oneUse", false) ? (1.0, 0.0, 0.0, 1.0) : (1.0, 1.0, 1.0, 1.0)
    Ahorn.drawRectangle(ctx, 0, 0, width, height, fillColor, lineColor)

    Ahorn.restore(ctx)

    direction = lowercase(get(entity.data, "direction", "up"))
    arrowSprite = Ahorn.getSprite(arrows[lowercase(direction)], "Gameplay")
    Ahorn.drawImage(ctx, arrowSprite, div(width - arrowSprite.width, 2), div(height - arrowSprite.height, 2))
end

end

max_stars_repo_name: FlynnSC/CommunalHelper

max_issues_repo_name: FlynnSC/CommunalHelper

max_forks_repo_name: FlynnSC/CommunalHelper


text:
using Colors
using Gadfly
using Random

grad = [RGB24(r/255, 128/255, b/255) for b=255:-1:0, r=0:255]

cm = colormap("Blues",11)
xs = 0:0.01:1

Random.seed!(1234)
n = 8
kf = 18 .^ rand(Float64, n) * 1.5
kf[1] = 0.0
ka1 = rand(Float64, n) * 2
kp1 = rand(Float64, n) * 6
ka2 = rand(Float64, n) * 2
kp2 = rand(Float64, n) * 6
f(x, t) = sum(i->(ka2[i] * t + ka1[i] * (1-t)) * sin(kf[i] * x + (kp2[i] * t + kp1[i] * (1-t))), 1:n)

p = plot(
    Coord.cartesian(ymin=-10, ymax=10),
    Guide.xticks(ticks=nothing), Guide.yticks(ticks=nothing),
    layer(y=f.(xs,0.0),x=xs, Geom.line, Theme(default_color=cm[1])),
    layer(y=f.(xs,0.1),x=xs, Geom.line, Theme(default_color=cm[2])),
    layer(y=f.(xs,0.2),x=xs, Geom.line, Theme(default_color=cm[3])),
    layer(y=f.(xs,0.3),x=xs, Geom.line, Theme(default_color=cm[4])),
    layer(y=f.(xs,0.4),x=xs, Geom.line, Theme(default_color=cm[5])),
    layer(y=f.(xs,0.5),x=xs, Geom.line, Theme(default_color=cm[6])),
    layer(y=f.(xs,0.6),x=xs, Geom.line, Theme(default_color=cm[7])),
    layer(y=f.(xs,0.7),x=xs, Geom.line, Theme(default_color=cm[8])),
    layer(y=f.(xs,0.8),x=xs, Geom.line, Theme(default_color=cm[9])),
    layer(y=f.(xs,0.9),x=xs, Geom.line, Theme(default_color=cm[10])),
    layer(y=f.(xs,1.0),x=xs, Geom.line, Theme(default_color=cm[11])),
)

p |> SVG("waves.svg", 100mm, 100mm)


text:
module CalleeError

inner(A, i) = A[i]
function outer(A)
    s = zero(eltype(A))
    for i = 1:length(A)+1
        s += inner(A, i)
    end
    return s
end

s = outer([1,2,3])

foo(x::Float32) = 1   # this line is tested for being on line number 14

end


text:
tup = (1,2,3)
@show typeof(tup)

@show tup[1]
for t in tup
    println(t)
end

#=
tup[1]=9
ERROR: LoadError: MethodError: no method matching setindex!(::Tuple{Int64,Int64,Int64}, ::Int64, ::Int64)
=#

namedtup = (integer = 3, realnumber=3.14,complexnumber=2+3im)
@show namedtup.integer
@show namedtup.realnumber
@show namedtup.complexnumber
@show namedtup[:integer]
@show namedtup[:realnumber]
@show namedtup[:complexnumber]

for k in keys(namedtup)
    println(k)
end

for value in values(namedtup)
    println(value)
end

for value in namedtup
    println(value)
end

for (k,v) in pairs(namedtup)
    println(k,v)
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer