text:
function Base.string(a::ArbUInt; base::Int=16, pad::Int=1)
    bufSize = max(pad, stringsize(a, base))
    buf = Base.StringVector(bufSize)
    stringify!(buf, a, base)
    return String(buf)
end

"""
    stringsize(num::ArbUInt, base)

Returns the number of characters required to represent `num` in the given `base`.
"""
function stringsize(a::ArbUInt, base)
    base == 1 && throw(ArgumentError("Base has to be >= 2"))

    # trivial case
    iszero(a) && return 1

    if base == 2
        return length(a.data)*sizeof(ArbDigit)*8 - leading_zeros(a)
    elseif base == 16
        retsize = (length(a.data)-1)*16 # remaining data
        nleading_bits = sizeof(ArbDigit)*8 - leading_zeros(a)
        nnibblesLB = div(nleading_bits, 4)
        retsize += nnibblesLB
        retsize += nnibblesLB*4 < nleading_bits
        return retsize
    else
        throw(ArgumentError("stringsize for base $base is not implemented yet"))
    end
end

"""
    stringify!(buffer, num::ArbUInt, base)

Writes a string representation of `num` in base `base` to `buffer`. If the string representation of `num` is shorter than `buffer`, leading entries in `buffer` are filled with `'0'`.
"""
function stringify!(buf, a::ArbUInt, base)
    requiredSize = stringsize(a, base)
    length(buf) < requiredSize && throw(ArgumentError("buffer not large enough to hold string"))
    if base == 16
        bufIdx = length(buf) - requiredSize + 1
        buf[begin:bufIdx] .= 0x30
        iszero(a) && return

        numIdx = lastindex(a.data)
        digit = a.data[numIdx]
        # write the leading digit specially - we may not want to pad here
        shiftFactor = div(leading_zeros(digit), 4)
        sh = sizeof(ArbDigit)*8 - 4*shiftFactor
        @inbounds while !iszero(sh)
            sh -= 4
            lo = ((digit >> sh) % UInt8) & 0x0f
            # 'a' - '0' - 0xa = 0x27
            buf[bufIdx] = lo + 0x30 + (lo >= 0xa)*0x27
            bufIdx += 1
        end
        numIdx -= 1

        # now for the remaining data
        nwrites = 2*sizeof(ArbDigit)
        @inbounds while numIdx > 0
            sh = sizeof(ArbDigit)*8 - 4
            digit = a.data[numIdx]
            for _ in 0:nwrites-1
                lo = ((digit >> sh) % UInt8) & 0x0f
                buf[bufIdx] = lo + 0x30 + (lo >= 0xa)*0x27
                sh -= 4
                bufIdx += 1
            end
            numIdx -= 1
        end
    else
        throw(ArgumentError("stringification for base $base is not implemented yet"))
    end
end

function Base.parse(::Type{ArbUInt}, str::AbstractString; base=16)
    !isascii(str) && throw(ArgumentError("Given string is not ASCII"))
    if base == 16
        !all(isxdigit, str) && throw(ArgumentError("Non-hex digit detected"))
    end
    parse(ArbUInt, codeunits(str); base)
end

function Base.parse(::Type{ArbUInt}, units::AbstractVector{UInt8}; base::Int=16)
    isempty(units) && return throw(ArgumentError("input string is empty"))
    # TODO: parsing assumes lower case for now

    if isone(length(units))
        lo = units[1]
        return ArbUInt(parse(ArbDigit, lo; base))
    end

    # trim leading zeros
    if units[1] == 0x30 && units[2] == 0x30 # zero?
        nonZeroIdx = findfirst(!=(0x30), units)
        nonZeroIdx === nothing && return zero(ArbUInt)
        units = @view units[nonZeroIdx-1:end]
    end

    parsed = if base == 16
        # one hex digit fills 4 bits => 2 hex chars per byte
        # we have sizeof(ArbDigit) bytes, so we have sizeof(ArbDigit)*2
        unitsPerDigit = sizeof(ArbDigit)*2
        if units[1] == 0x30 && units[2] == 0x78
            # there may be a `0x` prefix
            units = @view units[3:end]
        end

        len_data = length(units)
        nWholeDigits = div(len_data, unitsPerDigit)
        shortDigits = len_data - nWholeDigits*unitsPerDigit

        data = Vector{ArbDigit}(undef, nWholeDigits + !iszero(shortDigits))
        dataIdx = 1

        @inbounds for unitIdx in unitsPerDigit:unitsPerDigit:len_data
            slice = @view units[end-unitIdx+1:end-unitIdx+unitsPerDigit]
            digit = zero(ArbDigit)
            for i in 0:2:unitsPerDigit-1
                lo = slice[end-i]
                hi = slice[end-i-1]
                digit |= ((lo - 0x30 - (lo >= 0x61)*0x27) % ArbDigit) << (i*4)
                digit |= ((hi - 0x30 - (hi >= 0x61)*0x27) % ArbDigit) << (i*4+4)
            end
            data[dataIdx] = digit
            dataIdx += 1
        end

        @inbounds if !iszero(shortDigits)
            restIdxStart = 1
            leftover = @view units[restIdxStart:restIdxStart+shortDigits-1]
            digit = zero(ArbDigit)
            for i in 0:lastindex(leftover)-1
                nibble = leftover[end-i]
                digit |= ((nibble - 0x30 - (nibble >= 0x61)*0x27) % ArbDigit) << (i*4)
            end
            data[dataIdx] = digit
        end

        data
    else
        throw(ArgumentError("parsing for base $base is not implemented yet"))
    end

    return ArbUInt(parsed)
end


text:
const K_List = K_Vector{TI0.number,K_,K}
Base.setindex!(x::K_List, el, i::Integer) = (r0(x.a[i]); x.a[i] = K_new(el))
function Base.push!(x::K_List, el)
    n′ = length(x) + 1
    a = K_new(el)
    p = K_(pointer(x)-16)
    p′ = jk(Ref{K_}(p), a)
    resetvector(x.a, n′, p′+16)
    x
end

ext: jl


text:
@testset verbose=true "Projected CG                " begin
	n = 1000
	m = 10

	# construct some random data
	A = 0.01*randn(n, n)
	A = A*A' + 0.5*I

	b = randn(n)
	c = randn(m)

	F = qr(randn(n, m))
	U = Matrix(F.Q)

	x = zeros(n)
	λ = zeros(m)

	bigMat = [A U; U' zeros(m,m)]
	rhs = vcat(b, c)

	@testset "Accuracy" begin
		# do pcg with different tolerances and test for accuracy
		for tol in 10.0.^(-6:-1:-20)
			_, nr = LFPSQP.projcg!(x, λ, A, U, b, c, tol=tol)

			@test nr < tol
			@test norm(U'*x - c) < 1e-14
			@test norm(bigMat*vcat(x, λ) - rhs) < max(tol, 1e-13)
		end
	end

	@testset "Memory allocation" begin
		work = LFPSQP.ProjCGWork(n, m)
		mem = @allocated LFPSQP.projcg!(x, λ, A, U, b, c, tol=1e-20, work=work)

		@test mem == 0
	end

	@testset "Negative direction" begin
		# generate symmetric matrix with positive and negative eigenvalues
		F = qr(randn(n, n))
		S = Matrix(F.Q)
		Λ = vcat(rand(n-2*m) .+ 1, -1 .- rand(2*m))

		A = S*diagm(0 => Λ)*S'

		c = zeros(m)	# "negative direction" really only makes sense when c is 0

		i, nr = LFPSQP.projcg!(x, λ, A, U, b, c, tol=1e-20)

		@test isinf(nr)
		@test all(isnan.(λ))
		@test norm(U'*x - c) < 1e-14
		@test dot(x, A*x) <= 0.0
	end
end


text:
import MLJBase
import DataFrames
import CategoricalArrays

mutable struct XGBClassifier <: MLJBase.Probabilistic
    maxdepth::Integer
    nrounds::Integer
 end

function XGBClassifier(; maxdepth = 6, nrounds = 2)
    model = XGBClassifier(maxdepth, nrounds)
    return model
end

function MLJBase.scitypes(X::Vector{<:AbstractFactor})
    [MLJBase.Multiclass for x in X]
end

function MLJBase.scitype_union(y::AbstractFactor)
    MLJBase.Multiclass
end

function MLJBase.scitype_union(y::AbstractCovariate)
    MLJBase.Multiclass
end

function MLJBase.target_scitype_union(model::XGBClassifier)
    MLJBase.Multiclass
end

function MLJBase.selectrows(X::Vector{<:AbstractFactor}, r)
    (X, r)
end

function MLJBase.selectrows(X::JuML.DataFrame, r)
    (X, r)
end

function MLJBase.selectrows(y::AbstractFactor, r)
    (y, r)
end

function MLJBase.selectrows(y::AbstractCovariate, r)
    (y, r)
end

function MLJBase.fit(model::XGBClassifier, verbosity::Integer, X::DataFrames.DataFrame,
                     y::CategoricalArrays.CategoricalArray)

    X = JuML.DataFrame(X)
    levels = [string(level) for level in CategoricalArrays.index(y.pool)]
    T = eltype(y.refs)
    y = Factor{T}("y", levels, y.refs)
    label = covariate(y, (level -> level == levels[1] ? 1.0 : 0.0))
    trainsel = BoolVariate("", BitArray{1}(undef, 0))
    fitres = xgblogit(label, X.factors; trainselector = trainsel, maxdepth = model.maxdepth,
                      nrounds = model.nrounds)
    return fitres, nothing, nothing
end

function MLJBase.fit(model::XGBClassifier, verbosity::Integer, X::Tuple{Vector{<:AbstractFactor}, Any},
                     y::Tuple{AbstractCovariate, Any})

    trainsel = X[2] <: Colon ?
        begin
            BoolVariate("", BitArray{1}(undef, 0))
        end :
        begin
            n = length(y[1])
            b = falses(n)
            for i in X[2]
                b[i] = true
            end
            BoolVariate("", b)
        end

    factors = X[1]
    label = y[1]
    fitres = xgblogit(label, factors; trainselector = trainsel, maxdepth = model.maxdepth,
                      nrounds = model.nrounds)
    return fitres, nothing, nothing

end

function MLJBase.fit(model::XGBClassifier, verbosity::Integer, X::Tuple{JuML.DataFrame, Any},
                     y::Tuple{AbstractCovariate, Any})

    factors = X[1].factors
    MLJBase.fit(model, verbosity, (factors, X[2]), y)
end

function MLJBase.fit(model::XGBClassifier, verbosity::Integer, X::Tuple{JuML.DataFrame, Any},
                     y::Tuple{AbstractFactor, Any})

    factors = X[1].factors
    ylevels = getlevels(y[1])
    label = covariate(getname(y), (level -> level == ylevels[1] ? 1.0 : 0.0))
    MLJBase.fit(model, verbosity, (factors, X[2]), (label, y[2]))

end

function MLJBase.predict(model::XGBClassifier, fitresult, Xnew::JuML.DataFrame)

    predict(fitresult, Xnew)

end

function MLJBase.predict(model::XGBClassifier, fitresult, Xnew)

    predict(fitresult, JuML.DataFrame(Xnew))

end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
module NonLinearSolversTests

using Test
using Gridap
using ..NonLinearOperatorMocks

op = NonLinearOperatorMock()

x = [2.0, 3.0]
x = [1.0, 3.0]
@test residual(op,x) ≈ [0.0,0.0]
@test jacobian(op,x) ≈ [-1.0 0.0; 0.0 1.0]
@test create_in_domain(op) == zeros(2)

ls = LUSolver()
tol = 1.e-10
maxiters = 20
nls = NewtonRaphsonSolver(ls,tol,maxiters)

x = solve(nls,op)
@test x ≈ [1.0, 3.0]

x = [2.1,2.9]
cache = solve!(x,nls,op)
@test x ≈ [2.0, 3.0]

x = [2.1,2.9]
solve!(x,nls,op,cache)
@test x ≈ [2.0, 3.0]

end


text:
import Statistics
using FastGaussQuadrature

struct TIEstimate
    burn_in::Int
    integration_nodes::Int
    num_samples::Int
end

name(x::TIEstimate) = "TI"

struct ThermodynamicIntegrationResult{V} <: SimulationResult
    integration_weights::Vector{Float64}
    inv_temps::Vector{Float64}
    energies::Array{V,2}
    acceptance::Array{Bool,2}
end

function blocks(result::ThermodynamicIntegrationResult, block_size=2^10)
    block_averages(array) = map(mean, Iterators.partition(array, block_size))
    mapreduce(block_averages, hcat, eachcol(result.energies))
end

# perform the quadrature integral
log_marginal(result::ThermodynamicIntegrationResult{<:Real}) = dot(result.integration_weights, vec(mean(result.energies, dims=1)))
log_marginal(result::ThermodynamicIntegrationResult{<:Vector}) = sum(result.integration_weights .* vec(mean(result.energies, dims=1)))

function Statistics.var(result::ThermodynamicIntegrationResult, block_size=2^10)
    b = blocks(result, block_size)
    σ² = var(b, dims=1) ./ size(b, 1)
    dot(result.integration_weights.^2, σ²)
end

function summary(results::ThermodynamicIntegrationResult...)
    block_size = 2^10
    energy_blocks = [blocks(r, block_size) for r in results]
    acceptance = [mean(r.acceptance, dims=1) for r in results]
    inv_temps = [r.inv_temps for r in results]
    integration_weights = [r.integration_weights for r in results]

    DataFrame(EnergyBlocks=energy_blocks, Acceptance=acceptance, Theta=inv_temps, IntegrationWeights=integration_weights)
end

# Monte-Carlo computation of the marginal probability for the given configuration
function simulate(algorithm::TIEstimate, initial, system; kwargs...)
    # Generate the array of θ values for which we want to simulate the system.
    # We use Gauss-Legendre quadrature which predetermines the choice of θ.
    nodes, weights = gausslegendre(algorithm.integration_nodes)
    θrange = 0.5 .* nodes .+ 0.5
    # The factor 0.5 comes from rescaling the integration limits from [-1,1] to [0,1].
    weights = 0.5 .* weights

    initial_energy = energy_difference(initial, system)
    energies = fill(initial_energy, algorithm.num_samples, length(θrange))
    accept = Array{Bool}(undef, algorithm.num_samples, length(θrange))
    for i in eachindex(θrange)
        chn = chain(system; θ=θrange[i], kwargs...)
        sampler = MetropolisSampler(algorithm.burn_in, 0, energy(initial, chn), copy(initial), chn)
        for (j, was_accepted) in Iterators.enumerate(Iterators.take(sampler, algorithm.num_samples))
            energies[j, i] = -energy_difference(sampler.state, system)
            accept[j, i] = was_accepted != 0
        end
        initial = sampler.state
    end

    ThermodynamicIntegrationResult(weights, θrange, energies, accept)
end


text:
using Test

include("test_Samsara.jl")

max_stars_repo_name: leikanger/Samsara

max_issues_repo_name: leikanger/Samsara

max_forks_repo_name: leikanger/Samsara


text:
const VecI  = AbstractVector # Input  Vector
const VecO  = AbstractVector # Output Vector
const VecB  = AbstractVector # Buffer Vector
const VecIO = AbstractVector # In/Out Vector
const MatI  = AbstractMatrix # Input  Matrix
const MatO  = AbstractMatrix # Output Matrix
const MatB  = AbstractMatrix # Buffer Matrix
const MatIO = AbstractMatrix # In/Out Matrix

using Printf

import LinearAlgebra: BLAS
BLAS.set_num_threads(4)

import PyPlot as plt
import ApproxInference
import ApproxInference: @get, @objective, unsafe_copy!, diagmax, nrm2, dot, cubic,
                        lm_get_residual!, lm_get_jacobian!, vi_get_apprhess!,
                        lm_get_gradient!, lm_get_cholesky!, lm_get_blastrsv!,
                        lm_get_squaresum, vi_get_evidence!, vi_get_gradient!,
                        VariationalInferenceOptimizer

meshgrid(x::VecI{T}, y::VecI{T}) where T = Matrix{T}(undef, length(y), length(x))

function init_rc!()
    rcParams = plt.PyDict(plt.matplotlib."rcParams")
    pgfpreamble = string("\\usepackage{amsthm}",
                         "\\usepackage{mathtools}",
                         "\\usepackage[libertine]{newtxmath}",
                         "\\usepackage[tt=false]{libertine}",
                         "\\usepackage{siunitx}",
                         "\\usepackage{mhchem}")
    # Backend
    rcParams["backend"] = "qt5cairo"
    # Latex
    rcParams["pgf.rcfonts"] = false
    rcParams["pgf.preamble"] = pgfpreamble
    rcParams["text.usetex"] = true
    rcParams["text.latex.preamble"] = "\\usepackage{siunitx}\\usepackage{mhchem}"
    # Font
    rcParams["font.size"] = 14
    rcParams["font.family"] = "serif"
    # Axes
    rcParams["axes.linewidth"] = 1.5
    rcParams["xtick.major.width"] = 1.5
    rcParams["ytick.major.width"] = 1.5
    # Grid
    rcParams["grid.alpha"] = 0.5
    rcParams["grid.linewidth"] = 1.5
    # Lines
    rcParams["lines.linewidth"] = 2.0
    # Legend
    rcParams["legend.framealpha"] = 0.0
    return nothing
end

@objective struct BoxBOD
    x::Vector{Float64}; y::Vector{Float64}

    BoxBOD() = new([1., 2., 3., 5., 7., 10.], [109., 149., 149., 191., 213., 224.])
end

function ApproxInference.get_residual!(r::VecIO, θ::VecI, f::BoxBOD)
    x, y = f.x, f.y
    @inbounds θ1 = θ[1]
    @inbounds θ2 = θ[2]
    @simd for i in eachindex(x)
        @inbounds r[i] = y[i] - θ1 * (1.0 - exp(-θ2 * x[i]))
    end
    return nothing
end

function ApproxInference.get_jacobian!(J::MatIO, θ::VecI, f::BoxBOD)
    x = f.x
    @inbounds θ1 = θ[1]
    @inbounds θ2 = θ[2]
    @inbounds for i in eachindex(x)
        tmp = exp(-θ2 * x[i])
        J[i,1] = 1.0 - tmp
        J[i,2] = θ1 * x[i] * tmp
    end
    return nothing
end

function main(fname::String, μ0::VecI, Λ0::MatI, Λy::Union{VecI, MatI}, τ::Real, h::Real, itmax::Int)
    fn = BoxBOD()
    vi = VariationalInferenceOptimizer(2, 6)
    x1 = Float64[]; @inbounds push!(x1, μ0[1])
    x2 = Float64[]; @inbounds push!(x2, μ0[2])
    #### Allocations ####
    @get vi nd ny rs rt μs μt dμ δμ Δμ
    @get vi Js Jt ΛJ Λs Λt Λb Ft Fb
    #### Initialization
    ih1 = inv(h)
    ih2 = ih1 * ih1
    toN = eachindex(1:nd)
    unsafe_copy!(μs, μ0)
    unsafe_copy!(Λs, Λ0)
    #### Compute precision matrix
    lm_get_residual!(rs, μs, fn)
    lm_get_jacobian!(Js, μs, fn)
    vi_get_apprhess!(Λt, Λy, Js, Λ0, ΛJ) # ΛJ = Λy * Js + Λ0
    # Controlling parameters of trust region
    a = τ * diagmax(Λt, nd); b = 2.0
    #### Compute free energy
    fill!(Δμ, 0.0)
    unsafe_copy!(Ft, Λs)
    unsafe_copy!(Fb, Λt)
    Fnow = vi_get_evidence!(Λy, rs, ny, Λ0, Δμ, nd, Ft, Fb)
    #### Iteration
    it = 0
    while it < itmax
        it += 1
        # Damped Gauss-Newton approximation
        unsafe_copy!(Ft, Λt)
        λ = a * diagmax(Ft, nd)
        @simd for i in toN
            @inbounds Ft[i,i] += λ
        end
        cholesky_state = lm_get_cholesky!(Ft)
        # Levenberg-Marquardt step by 1st order linearization
        vi_get_gradient!(dμ, ΛJ, rs, Λ0, Δμ)
        lm_get_blastrsv!(dμ, Ft)
        # Finite-difference of 2nd order directional derivative
        @simd for i in toN
            @inbounds μt[i] = μs[i] + h * dμ[i]
        end
        lm_get_residual!(rt, μt, fn)
        lm_get_gradient!(rt, rs, Js, dμ, ih1, ih2)
        # Geodesic acceleration step by 2nd order linearization
        lm_get_gradient!(δμ, ΛJ, rt)
        lm_get_blastrsv!(δμ, Ft)
        gratio = 2.0 * nrm2(δμ) / nrm2(dμ)
        gratio > 0.75 && (a *= b; b *= 2.0; continue)
        # (Levenberg-Marquardt velocity) + (Geodesic acceleration)
        BLAS.axpy!(1.0, δμ, dμ) # dμ ← dμ + δμ
        # Compute predicted gain of the least square cost by linear approximation
        LinApprox = lm_get_squaresum(Λt, dμ, nd) + λ * dot(dμ, dμ, nd)
        # Compute gain ratio
        @simd for i in toN
            @inbounds μt[i] = μs[i] + dμ[i]
        end
        lm_get_residual!(rt, μt, fn)
        lm_get_jacobian!(Jt, μt, fn)
        vi_get_apprhess!(Λb, Λy, Jt, Λ0, ΛJ)
        @simd for i in toN
            @inbounds δμ[i] = μt[i] - μ0[i]
        end
        unsafe_copy!(Ft, Λt)
        unsafe_copy!(Fb, Λb)
        Fnew = vi_get_evidence!(Λy, rt, ny, Λ0, δμ, nd, Ft, Fb)

        Δ = Fnow - Fnew; ρ = Δ / LinApprox
        if ρ > 0.0
            Fnow = Fnew
            unsafe_copy!(μs, μt)
            ###
            @inbounds push!(x1, μs[1])
            @inbounds push!(x2, μs[2])
            ###
            unsafe_copy!(Λs, Λt)
            (maximum(δμ) < 1e-15 || Δ < 1e-10) && break # Check localized convergence
            unsafe_copy!(rs, rt)
            unsafe_copy!(Λt, Λb)
            unsafe_copy!(Js, Jt)
            unsafe_copy!(Δμ, δμ)
            a = ρ < 0.9367902323681495 ? a * (1.0 - cubic(2.0 * ρ - 1.0)) : a / 3.0
            b = 2.0
        else
            a > 1e10 && break # Check localized divergence
            a *= b; b *= 2.0
        end
    end

    contourX = collect(range(-50.0, 250.0; step=1.0)) # length(X) = N
    contourY = collect(range(-0.10, 14.0; step=0.05)) # length(Y) = M
    contourZ = meshgrid(contourX, contourY)           # size(Z) = (M, N)

    for j in eachindex(contourX)
        @inbounds μs[1] = contourX[j]
        @inbounds for i in eachindex(contourY)
            @inbounds μs[2] = contourY[i]
            lm_get_residual!(rs, μs, fn)
            contourZ[i,j] = log10(lm_get_squaresum(Λy, rs, ny))
        end
    end

    init_rc!()

    _cmap = plt.get_cmap("GnBu")
    fig   = plt.figure(num=1, figsize=(8.8, 6.6), clear=true, linewidth=6.)
    ax1   = fig.add_subplot(111)
    _fill = ax1.contourf(contourX, contourY, contourZ, cmap=_cmap, levels=50)
    _line = ax1.plot(x1, x2, "--*", c="k", lw=1)

    ax1.set_xlabel(raw"$\theta_{1}$")
    ax1.set_ylabel(raw"$\theta_{2}$")

    fig.colorbar(_fill, pad=0.02)
    ax1.patch.set_alpha(0.)
    fig.patch.set_alpha(0.)
    fig.set_tight_layout(true)
    fig.savefig(fname, format="pdf", dpi=600, backend="pgf", bbox_inches="tight", pad_inches=0.05)
    return nothing
end

function main()
    parμ = [10.0, 5.5]
    parΛ = zeros(Float64, 2, 2)
    βmat = zeros(Float64, 6, 6)

    @simd for i in eachindex(1:6)
        @inbounds βmat[i,i] = 1.0
    end

    for λ in (1e-7, 1e-5, 1e-3, 1e-1, 5e-1, 1e0)
        @simd for i in eachindex(parμ)
            @inbounds parΛ[i,i] = λ
        end
        fname = @sprintf "./BoxBOD by Variational-Inference Method λ=%.1e.pdf" λ
        main(fname, parμ, parΛ, βmat, 1e-3, 0.1, 100)
    end

    return nothing
end

main()


text:
# E. Dynamic Parallelism

## streams

export CuDeviceStream

struct CuDeviceStream
    handle::cudaStream_t

    function CuDeviceStream(flags=cudaStreamNonBlocking)
        handle_ref = Ref{cudaStream_t}()
        cudaStreamCreateWithFlags(handle_ref, flags)
        return new(handle_ref[])
    end
end

Base.unsafe_convert(::Type{cudaStream_t}, s::CuDeviceStream) = s.handle

function unsafe_destroy!(s::CuDeviceStream)
    cudaStreamDestroy(s)
    return
end

## execution

# device-side counterpart of launch
@inline function device_launch(f, blocks, threads, shmem, stream, args...)
    blocks = CuDim3(blocks)
    threads = CuDim3(threads)

    buf = parameter_buffer(f, blocks, threads, shmem, args...)
    cudaLaunchDeviceV2(buf, stream)

    return
end

@generated function parameter_buffer(f, blocks, threads, shmem, args...)
    # allocate a buffer
    ex = quote
        Base.@_inline_meta
        buf = cudaGetParameterBufferV2(f, blocks, threads, shmem)
    end

    # store the parameters
    #
    # > Each individual parameter placed in the parameter buffer is required to be aligned.
    # > That is, each parameter must be placed at the n-th byte in the parameter buffer,
    # > where n is the smallest multiple of the parameter size that is greater than the
    # > offset of the last byte taken by the preceding parameter. The maximum size of the
    # > parameter buffer is 4KB.
    offset = 0
    for i in 1:length(args)
        T = args[i]
        align = sizeof(T)
        buf_index = Base.ceil(Int, offset / align) + 1
        offset = buf_index * align
        ptr = :(Base.unsafe_convert(Ptr{$T}, buf))
        push!(ex.args, :(
            Base.pointerset($ptr, args[$i], $buf_index, $align)
        ))
    end

    push!(ex.args, :(return buf))

    return ex
end

## synchronization

"""
    device_synchronize()

Wait for the device to finish. This is the device side version,
and should not be called from the host.

`device_synchronize` acts as a synchronization point for
child grids in the context of dynamic parallelism.
"""
device_synchronize() = cudaDeviceSynchronize()


text:
using Juno
using Base.Test

@test Juno.isactive() == false

@test_throws Exception selector(["foo", "bar", "baz"])

let i = 0
  @progress for _ = 1:100
    i += 1
  end
  @test i == 100
end


text:
#--------------------- MODEL ----------------------------------------

function readFixedFile(dfx::DataFrame,mname::String)
    println("LOADING DFX Fixed")
    vout = dfx[(dfx[:modelType].=="GLM")&(dfx[:model].==mname),[:parameter,:coef,:stderr,:zval,:pval]]
    names!(vout,[:x,:estimate, :std_error, :z_value,:pr_z_])
    return vout
end
#readFixedFile(dfx,"occ")

function readRandFile(dfx::DataFrame,mname::String)
    println("LOADING DFX Rand")
    vout = dfx[(dfx[:modelType].=="GLMM")&(dfx[:model].==mname),[:ranef,:parameter,:coef,:stderr,:zval,:pval]]
    vout[:adj_coef] = 0.0
    vout[:adj_stderr]=  0.0
    vout[:row] = 0
    vout[:levelse] = 0
    vout[:exposed_orig]= 1
    vout[:exposed] = true
    vout[:key] = map( (c,l) ->  c*" ("*l*")"  , vout[:ranef],vout[:parameter] )
    for row in eachrow(vout)
        n=vout[(vout[:ranef].==row[:ranef])&(vout[:parameter].=="none"),:coef][1]
        s=vout[(vout[:ranef].==row[:ranef])&(vout[:parameter].=="none"),:stderr][1]
        println(row[:ranef]," ~~ ",n," ~~ ",s)
        row[:adj_coef] = row[:coef]-n
        row[:adj_stderr] = sqrt(row[:stderr]^2+s^2)
        row[:coef] = 0
    end
    vout=vout[[  :ranef,:row,:parameter,:coef,:adj_coef,:levelse,:stderr,:adj_stderr,:exposed_orig,:exposed,:key ]]
    names!(vout,[:class,:row,:level,    :B0  ,:B1      ,:levelse,:SE0   ,:SE1       ,:exposed_orig,:exposed,:key ])
    return vout[vout[:level].!="none",:]
end
#readRandFile(dfx,"occ")

"""
function readFixedFile(fname::AbstractString)
    if isfile(fname)
        idf = readtable(fname, header=true)
        dfnames = [:x,:Estimate,:Std_Error,:z_value,:Pr_z_]
        dfnames2 = [:x,:Estimate,:Std_Error,:t_value,:Pr_z_]
        dfnames3 = [:x,:Estimate,:Std_Error,:t_value,:Pr_t_]
        if (names(idf) == dfnames)|(names(idf) == dfnames2)|(names(idf) == dfnames3)
            idf[:x] = map(x-> lowercase(replace(replace(x,"(",""),")","")), idf[:x])
            names!(idf,[:x,:estimate,:std_error,:z_value,:pr_z_])
        else
            error("Invalid Headers in Fixed Src File :",fname)
        end
    else
        error("Can't find file: ",fname)
    end
    return idf
end

function readRandFile(fname::AbstractString, cfg::OrderedDict)
    defaultDF=DataFrame(class=Any[],row=Int64[],level=UTF8String[],B0=Float64[],B1=Float64[],levelse=Int64[],SE0=Float64[],SE1=Float64[],exposed_orig=Int64[],exposed=Any[],key=UTF8String[])

    if isfile(fname)
        idf = readtable(fname, header=true)
        dfnames = [:class,:row,:level,:group0,:group1,:levelSE,:SE_group0,:SE_group1,:exposed]
        dfnames2  = [:class,:row,:level,:group0,:group1,:levelSE,:SE_group0,:SE_group1,:exposed]   ## Same for now
        if (names(idf) == dfnames)|(names(idf) == dfnames2)
            idf[:class] = map(x-> lowercase(x), idf[:class])
            names!(idf,[:class,:row,:level,:B0,:B1,:levelse,:SE0,:SE1,:exposed_orig])
            idf[:exposed] = map(x-> x .==1 ? true:false, idf[:exposed_orig])
            # --- GROUP1 -----
#            idf[:Z1] = idf[:B1] ./ idf[:SE1]
#	        idf[:P1] = 2.0 * ccdf(Normal(), abs(idf[:Z1]))
#            idf[:SIG1] = map( x -> ((isnan(x))|(x > 0.2)) ?  false : true, idf[:P1] )
#            # --- GROUP0 -----
#            idf[:Z0] = idf[:B0] ./ idf[:SE0]
#            idf[:P0] = 2.0 * ccdf(Normal(), abs(idf[:Z0]))
#            idf[:SIG0] = map( x -> (isnan(x)|(x > 0.2)) ?  false : true, idf[:P0] )
            idf[:key] = map((c,l) ->  rkey(c,l), idf[:class],idf[:level])

            x=get(cfg,:excludedKeys, NA)
            if typeof(x) != NAtype
                for xv in x
                    idf=idf[idf[:key].!=xv,:]
                end
            end
            x=get(cfg,:excludedBreaks, NA)
            if typeof(x) != NAtype
                for xv in x
                    idf=idf[idf[:class].!=xv,:]
                end
            end
            x=get(cfg,:excludedLevels, NA)
            if typeof(x) != NAtype
                for xv in x
                    idf=idf[idf[:level].!=xv,:]
                end
            end

        else
            error("Invalid Headers in Random Src File",fname)
        end
    else
        #error("Can't find file: ",fname)
        println("Can't find file: ",fname)
        idf=defaultDF
    end
    return idf
end
"""

abstract MModel

function genScoreTotals(m::MModel, df_data::DataFrame, totalOnly::Bool)  # Still need to implement totalOnly -- cureently checks for reff_xxxx
    #println(methods(exp))
    if m.modelName == "occ"
        #df_data[:pre_occ_score0]=map(Float64,df_data[:pre_occ_score0])
        #df_data[:pre_occ_score1]=map(Float64,df_data[:pre_occ_score1])
        if :reff_occ in names(df_data)
            #df_data[:reff_occ] = map(Float64, df_data[:reff_occ])
            println("using reff : ",typeof(df_data[:pre_occ_score0]))
            df_data[:occ_score0] = exp(df_data[:pre_occ_score0] + df_data[:reff_occ]  )
            df_data[:occ_score1] = exp(df_data[:pre_occ_score1] + df_data[:reff_occ] )
        else
            df_data[:occ_score0] = exp(df_data[:pre_occ_score0])
            df_data[:occ_score1] = exp(df_data[:pre_occ_score1])
        end
    end
    #---------------------
        #DataArrays
    if m.modelName == "dolocc"
        if :reff_dolocc in names(df_data)
            println("using reff")
            df_data[:dolocc_score0] = exp(df_data[:pre_dolocc_score0] + df_data[:reff_dolocc]  )
            df_data[:dolocc_score1] = exp(df_data[:pre_dolocc_score1] + df_data[:reff_dolocc] )
        else
            df_data[:dolocc_score0] = exp(df_data[:pre_dolocc_score0])
            df_data[:dolocc_score1] = exp(df_data[:pre_dolocc_score1])
        end
    end
    #----------------------
    if m.modelName == "pen"
        if :reff_pen in names(df_data)
            println("using reff")
            df_data[:pen_score0] = exp((df_data[:pre_pen_score0] + df_data[:reff_pen])) ./ ( exp((df_data[:pre_pen_score0] + df_data[:reff_pen])) +1)
            df_data[:pen_score1] = exp((df_data[:pre_pen_score1] + df_data[:reff_pen])) ./ (exp((df_data[:pre_pen_score1] + df_data[:reff_pen])) +1)
        else
            df_data[:pen_score0] =  exp(df_data[:pre_pen_score0]) ./ ( exp(df_data[:pre_pen_score0]) +1)
            df_data[:pen_score1] = exp(df_data[:pre_pen_score1]) ./ ( exp(df_data[:pre_pen_score1]) +1)
        end
    end
end

type MOcc <: MModel
    hasBreaks::Bool
    fList::Vector{Symbol}
    rList::Vector{Symbol}
    dList::Vector{Symbol}
    feff::FOcc
    reff::ROcc
    fmula::AbstractString
    modelName::AbstractString
    logvar::AbstractString
    function MOcc(df_data::DataFrame,cfg::OrderedDict,dfx::DataFrame) this=new(); this.modelName="occ"; this.logvar="trps_pre_p1"; #this.sdf=pushSDFrow!(deepcopy(SDF),"Total Campaign",this.v_model);
        this.fList=Symbol[]; this.rList=Symbol[]; this.dList=Symbol[]
        #this.hasBreaks=false
        this.fmula=""
        idf = readFixedFile(dfx,"occ")
        #idf=readFixedFile("Occasion/Occ_fixed_effects.csv")
        this.feff=FOcc(df_data, idf)
        idf = readRandFile(dfx,"occ")
        #idf=readRandFile("Occasion/Occ_random_effects.csv",cfg)
        this.reff=ROcc(df_data, idf)
        this.hasBreaks = length(idf[1])==0 ? false : true
        genScoreTotals(this,df_data,cfg[:TotalModelsOnly])
        return this
    end

end

type MDolOcc <: MModel
    hasBreaks::Bool
    fList::Vector{Symbol}
    rList::Vector{Symbol}
    dList::Vector{Symbol}
    feff::FDolOcc
    reff::RDolOcc
    fmula::AbstractString
    modelName::AbstractString
    logvar::AbstractString
    function MDolOcc(df_data::DataFrame,cfg::OrderedDict,dfx::DataFrame) this=new(); this.modelName="dolocc"; this.logvar="dol_per_trip_pre_p1"; #this.sdf=pushSDFrow!(deepcopy(SDF),"Total Campaign",this.v_model);
        this.fList=Symbol[]; this.rList=Symbol[]; this.dList=Symbol[]
        this.fmula=""
        idf = readFixedFile(dfx,"dolocc")
        #idf=readFixedFile("DollarsOccasion/DolOcc_fixed_effects.csv")
        this.feff=FDolOcc(df_data, idf)
        idf = readRandFile(dfx,"dolocc")
        #idf=readRandFile("DollarsOccasion/DolOcc_random_effects.csv",cfg)
        this.reff=RDolOcc(df_data, idf)
        this.hasBreaks = length(idf[1])==0 ? false : true
        genScoreTotals(this,df_data,cfg[:TotalModelsOnly])
        return this
    end

end

type MPen <: MModel
    hasBreaks::Bool
    fList::Vector{Symbol}
    rList::Vector{Symbol}
    dList::Vector{Symbol}
    feff::FPen
    reff::RPen
    fmula::AbstractString
    modelName::AbstractString
    logvar::AbstractString
    function MPen(df_data::DataFrame,cfg::OrderedDict,dfx::DataFrame) this=new(); this.modelName="pen"; this.logvar="buyer_pre_p1"; #this.sdf=pushSDFrow!(deepcopy(SDF),"Total Campaign",this.v_model);
        this.fList=Symbol[]; this.rList=Symbol[]; this.dList=Symbol[]
        this.fmula=""
        idf = readFixedFile(dfx,"pen")
        #idf=readFixedFile("Penetration/Pen_fixed_effects.csv")
        this.feff=FPen(df_data, idf)
        idf = readRandFile(dfx,"pen")
        #idf=readRandFile("Penetration/Pen_random_effects.csv",cfg)
        this.reff=RPen(df_data, idf)
        this.hasBreaks = length(idf[1])==0 ? false : true
        genScoreTotals(this,df_data,cfg[:TotalModelsOnly])
        return this
    end

end

type MDolHH
    sdf::DataFrame
    mocc::MOcc
    mdolocc::MDolOcc
    mpen::MPen
    modelName::AbstractString
    function MDolHH(mocc::MOcc,mdolocc::MDolOcc,mpen::MPen) this=new(); this.modelName="dolhh";
        this.mocc=mocc
        this.mdolocc=mdolocc
        this.mpen=mpen
        this.sdf=pushSDFrow!(deepcopy(SDF),"Total Campaign",this.modelName);
        #by(vcat(mocc.reff.src[[:key, :class, :level,:exposed]], mdolocc.reff.src[[:key, :class, :level,:exposed]], mpen.reff.src[[:key, :class, :level,:exposed]]), [:key,:class,:level], df ->  maximum(df[:exposed]) )

        occlst = mocc.hasBreaks ? mocc.reff.src[:key] : DataArray[]
        dolocclst = mocc.hasBreaks ? mocc.reff.src[:key] : DataArray[]
        penlst = mocc.hasBreaks ? mocc.reff.src[:key] : DataArray[]

        for k in unique(vcat(occlst,dolocclst,penlst))
            pushSDFrow!(this.sdf,k,this.modelName);
        end
        #for k in unique(vcat(mocc.reff.src[[:key]], mdolocc.reff.src[[:key]], mpen.reff.src[[:key]]))[1]
        #    pushSDFrow!(this.sdf,k,this.modelName);
        #end

        delete!(this.sdf, :mean_score0)
        delete!(this.sdf, :mean_score1)
        delete!(this.sdf, :B)
        delete!(this.sdf, :SE)
        delete!(this.sdf, :P)
        delete!(this.sdf, :M)
        delete!(this.sdf, :Mt)
        delete!(this.sdf, :Mc)

        return this
    end
end

function Base.show(io::IO, mdolhh::MDolHH)
    show(io,mdolhh.sdf)
end

function Base.show(io::IO, mmodel::MModel)
    show(io,fieldnames(mmodel))
end

function genCnts(df_data::DataFrame, mocc::MOcc,mdolocc::MDolOcc,mpen::MPen,TotalModelsOnly::Bool)
    tot=DataFrame(key=["Total Campaign"],class=[""],level=[""],exposed=[true], M=[0.0],Mt=[0.0],Mc=[0.0],N=[0.0],Nt=[0.0],Nc=[0.0])
    tot[:M] =  length(df_data[ df_data[:buyer_pos_p1] .== 1 ,1])
    tot[:Mt] = length(df_data[ (df_data[:group] .== 1) & (df_data[:buyer_pos_p1] .== 1 )  ,1])
    tot[:Mc] = length(df_data[ (df_data[:group] .== 0) & (df_data[:buyer_pos_p1] .== 1 ) ,1])
    tot[:N] = length(df_data[1])
    tot[:Nc] = length(df_data[ df_data[:group] .== 0 ,1])
    tot[:Nt] = length(df_data[ df_data[:group] .== 1 ,1])

    if !TotalModelsOnly
        randlist = vcat(mocc.reff.src[[:key, :class, :level,:exposed]],mdolocc.reff.src[[:key, :class, :level,:exposed]],mpen.reff.src[[:key, :class, :level,:exposed]])
        if length(randlist[1]) > 0
            sdf = by(randlist, [:key,:class,:level], df ->  maximum(df[:exposed]) )
            #sdf = by(vcat(mocc.reff.src[[:key, :class, :level,:exposed]],mdolocc.reff.src[[:key, :class, :level,:exposed]],mpen.reff.src[[:key, :class, :level,:exposed]]), [:key,:class,:level], df ->  maximum(df[:exposed]) )
            rename!(sdf,:x1,:exposed)
            sdf[:M]=0.0
            sdf[:Mt]=0.0
            sdf[:Mc]=0.0
            sdf[:N]=0.0
            sdf[:Nt]=0.0
            sdf[:Nc]=0.0
            for row in eachrow(sdf)
                k = row[:key]
                ranfx = row[:class]
                v_level = row[:level]
                exposed = row[:exposed]
                sdf[sdf[:key].==k, :Mt] = length(df_data[ (df_data[:group] .== 1) & (df_data[:buyer_pos_p1] .== 1 ) & (df_data[Symbol(ranfx)] .== v_level) ,1])
                sdf[sdf[:key].==k, :M] = length(df_data[ (df_data[:buyer_pos_p1] .== 1) & (df_data[Symbol(ranfx)] .== v_level) ,1])
                sdf[sdf[:key].==k, :Nt] = length(df_data[ (df_data[:group] .== 1) & (df_data[Symbol(ranfx)] .== v_level) ,1])
                sdf[sdf[:key].==k, :N] = length(df_data[ (df_data[Symbol(ranfx)] .== v_level) ,1])
                if exposed    #  If Exposed - default Nc & Mc to total -- else count by rndfx
                    sdf[sdf[:key].==k, :Mc] = length(df_data[ (df_data[:group] .== 0) & (df_data[:buyer_pos_p1] .== 1 ) ,1])
                    sdf[sdf[:key].==k, :Nc] = length(df_data[ (df_data[:group] .== 0) ,1])
                    sdf[sdf[:key].==k, :M] = sdf[sdf[:key].==k, :Mt] + sdf[sdf[:key].==k, :Mc]
                    sdf[sdf[:key].==k, :N] = sdf[sdf[:key].==k, :Nt] + sdf[sdf[:key].==k, :Nc]
                else
                    sdf[sdf[:key].==k, :Mc] = length(df_data[ (df_data[:group] .== 0) & (df_data[:buyer_pos_p1] .== 1 ) & (df_data[Symbol(ranfx)] .== v_level) ,1])
                    sdf[sdf[:key].==k, :Nc] = length(df_data[ (df_data[:group] .== 0) & (df_data[Symbol(ranfx)] .== v_level) ,1])
                end
                println("-- ",ranfx,"-- ",v_level,"-- ", exposed )
            end
            return vcat(tot,sdf)

        else
            return tot
        end
    else
        return tot
    end
end
#genCnts(df_data, mocc,mdolocc,mpen)

type Cnts
    sdf::DataFrame
#    get::Function
    function Cnts(df_data::DataFrame, mocc::MOcc,mdolocc::MDolOcc,mpen::MPen,TotalModelsOnly::Bool)
        this = new()
        this.sdf = genCnts(df_data, mocc,mdolocc,mpen,TotalModelsOnly)
#        function getk(k::AbstractString)
#            return this.sdf[this.sdf[:key].==k,:]
#        end
#        function getk(k::AbstractString,s::Symbol)
#            return this.sdf[this.sdf[:key].==k,s][1]
#        end
#        this.get = getk
        return this
    end
end

function hhcounts(cnts::Cnts)
    sdf=cnts.sdf
    sdf[:hh]=0.0
    sdf[:impressions]=0.0
    sdf[:weight]=0.0
    df_hhc = readtable("hhcounts.csv",header=true); lowercase(df_hhc)
    df_hhc[:class] = map(x->lowercase(x)  ,df_hhc[:class])
    for row in eachrow(df_hhc[df_hhc[:class].!="total campaign",:])
        k=rkey(row)
        sdf[sdf[:key].==k,:hh] = row[:hh]
        sdf[sdf[:key].==k,:impressions] = row[:impressions]
    end
    #tot = sum(sdf[:hh])
    tot = by(sdf,:class, df -> sum(df[:hh]))[:x1][2]
    toti = sum(sdf[:impressions])
    sdf[sdf[:key].=="Total Campaign",:hh] = tot
    sdf[sdf[:key].=="Total Campaign",:impressions] = toti
    for row in eachrow(cnts.sdf[cnts.sdf[:key].!="Total Campaign",:])
        k=rkey(row)
        sdf[sdf[:key].==k,:weight] = row[:hh] / tot
        println(typeof(row[:hh]))
        println(row[:hh])
    end
end

function hhcounts_OLD(cnts::Cnts)
    sdf=cnts.sdf
    sdf[:hh]=0
    sdf[:impressions]=0
    sdf[:weight]=0.0
    df_hhc = readtable("hhcounts.csv",header=true); lowercase(df_hhc)
    df_hhc[:class] = map(x->lowercase(x)  ,df_hhc[:class])
    tot = df_hhc[df_hhc[:class].=="total campaign",:hh][1]
    toti = df_hhc[df_hhc[:class].=="total campaign",:impressions][1]
    sdf[sdf[:key].=="Total Campaign",:hh] = tot
    sdf[sdf[:key].=="Total Campaign",:impressions] = toti
    for row in eachrow(df_hhc[df_hhc[:class].!="total campaign",:])
        k=rkey(row)
        #println(k,"-",row[:level],"-",row[:hh])
        sdf[sdf[:key].==k,:hh] = row[:hh]
        sdf[sdf[:key].==k,:impressions] = row[:impressions]
        sdf[sdf[:key].==k,:weight] = row[:hh] / tot
        println(typeof(row[:hh]))
    end
    #sum(cnts.sdf[:hh])
    class=""
    for row in eachrow(cnts.sdf[cnts.sdf[:key].!="Total Campaign",:])
        if row[:class] == class

        else
            class=row[:class]
        end
        println(row[:hh])
    end
end

function getk(sdf::DataFrame,k::AbstractString)
    return sdf[sdf[:key].==k,:]
end
function getk(sdf::DataFrame,k::AbstractString,s::Symbol)
    return sdf[sdf[:key].==k,s][1]
end

#function getk(cnts::Cnts,k::AbstractString)
#    return cnts.sdf[cnts.sdf[:key].==k,:]
#end
#function getk(cnts::Cnts,k::AbstractString,s::Symbol)
#    return cnts.sdf[cnts.sdf[:key].==k,s][1]
#end

#function getk(mdf::MDF,k::AbstractString,cnts::Cnts)
#    mdf.sdf = genMDF(cnts,this.mocc,this.mdolocc, this.mpen)
#    return this.sdf[this.sdf[:key].==k,:]
#end

function genRndMDF(cnts::Cnts, mocc::MOcc, mdolocc::MDolOcc, mpen::MPen, withNA::Bool=true, includeRawData::Bool=false)
    mdfcols=[:key,:B1,:SE1,:P1,:o_mean_score0,:o_mean_score1,:inOcc,:isSig1,:B1_orig,:SE1_orig,:o_adj_mean_cntrl_grp, :o_adj_mean_expsd_grp,
        :B2,:SE2,:P2,:y_mean_score0,:y_mean_score1,:inDolOcc,:isSig2,:B2_orig,:SE2_orig,:y_adj_mean_cntrl_grp, :y_adj_mean_expsd_grp,
             :B3,:SE3,:P3,:p_mean_score0,:p_mean_score1,:inPen,:isSig3,:B3_orig,:SE3_orig,:p_adj_mean_cntrl_grp, :p_adj_mean_expsd_grp,
        :M,:Mt,:Mc,:N,:Nt,:Nc,:o_B0,:o_SE0,:y_B0,:y_SE0,:p_B0,:p_SE0
            ]
    if mocc.hasBreaks & mdolocc.hasBreaks & mpen.hasBreaks & isdefined(mocc, :reff) & isdefined(mdolocc, :reff) & isdefined(mpen, :reff)
        #or=deepcopy(mocc.reff.sdf[[:key,:B1_combo,:SE1,:SE1_combo,:P1_combo,:mean_score0,:mean_score1]])
        or=deepcopy(mocc.reff.sdf[[:key,:B1_combo,:SE1_combo,:P1_combo,:mean_score0,:mean_score1,:B0,:SE0,:B1,:SE1,:adj_mean_cntrl_grp, :adj_mean_expsd_grp]])  #:adj_mean_cntrl_grp, :adj_mean_expsd_grp
        rename!(or,:B1,:B1_orig)
        rename!(or,:SE1,:SE1_orig)

        rename!(or,:B1_combo,:B1)
        rename!(or,:P1_combo,:P1)
        #rename!(or,:SE1,:SE1)
        rename!(or,:SE1_combo,:SE1)
        rename!(or,:B0,:o_B0)
        rename!(or,:SE0,:o_SE0)

        rename!(or,:mean_score0,:o_mean_score0)
        rename!(or,:mean_score1,:o_mean_score1)

        rename!(or,:adj_mean_cntrl_grp,:o_adj_mean_cntrl_grp); rename!(or,:adj_mean_expsd_grp ,:o_adj_mean_expsd_grp)

        or[:inOcc] = true
        or[:isSig1] = true

        #yr=deepcopy(mdolocc.reff.sdf[[:key,:B1_combo,:SE1,:SE1_combo,:P1_combo,:mean_score0,:mean_score1]])
        yr=deepcopy(mdolocc.reff.sdf[[:key,:B1_combo,:SE1_combo,:P1_combo,:mean_score0,:mean_score1,:B0,:SE0,:B1,:SE1,:adj_mean_cntrl_grp, :adj_mean_expsd_grp]])
        rename!(yr,:B1,:B2_orig)
        rename!(yr,:SE1,:SE2_orig)

        rename!(yr,:B1_combo,:B2)
        rename!(yr,:P1_combo,:P2)
        #rename!(yr,:SE1,:SE2)
        rename!(yr,:SE1_combo,:SE2)
        rename!(yr,:B0,:y_B0)
        rename!(yr,:SE0,:y_SE0)

        rename!(yr,:mean_score0,:y_mean_score0)
        rename!(yr,:mean_score1,:y_mean_score1)

        rename!(yr,:adj_mean_cntrl_grp,:y_adj_mean_cntrl_grp); rename!(yr, :adj_mean_expsd_grp ,:y_adj_mean_expsd_grp)

        yr[:inDolOcc] = true
        yr[:isSig2] = true

        #pr=deepcopy(mpen.reff.sdf[[:key,:B1_combo,:SE1,:SE1_combo,:P1_combo,:mean_score0,:mean_score1]])
        pr=deepcopy(mpen.reff.sdf[[:key,:B1_combo,:SE1_combo,:P1_combo,:mean_score0,:mean_score1,:B0,:SE0,:B1,:SE1,:adj_mean_cntrl_grp, :adj_mean_expsd_grp]])
        rename!(pr,:B1,:B3_orig)
        rename!(pr,:SE1,:SE3_orig)

        rename!(pr,:B1_combo,:B3)
        rename!(pr,:P1_combo,:P3)
        #rename!(pr,:SE1,:SE3)
        rename!(pr,:SE1_combo,:SE3)

        rename!(pr,:B0,:p_B0)
        rename!(pr,:SE0,:p_SE0)

        rename!(pr,:mean_score0,:p_mean_score0)
        rename!(pr,:mean_score1,:p_mean_score1)

        rename!(pr,:adj_mean_cntrl_grp,:p_adj_mean_cntrl_grp); rename!(pr, :adj_mean_expsd_grp ,:p_adj_mean_expsd_grp)

        pr[:inPen] = true
        pr[:isSig3] = true

        ztmp= join(join(or,yr, on = :key, kind = :left),pr, on = :key, kind = :outer)   #:left
        z=join(ztmp,cnts.sdf[[:key,:M,:Mt,:Mc,:N,:Nt,:Nc,]], on = :key)[mdfcols]
        #z= join(join(or,yr, on = :key, kind = :left),pr, on = :key, kind = :left)[mdfcols]
        #showcols(z)
        z[find(isna(z[:inOcc])), :inOcc] = false
        z[find(isna(z[:isSig1])), :isSig1] = false
        z[find(isna(z[:inDolOcc])), :inDolOcc] = false
        z[find(isna(z[:isSig2])), :isSig2] = false
        z[find(isna(z[:inPen])), :inPen] = false
        z[find(isna(z[:isSig3])), :isSig3] = false
        #t=join(z,rcnts.sdf[[:key,:M,:Mt,:Mc,:N,:Nt,:Nc]], on = :key)
        #t[[:M,:M_1,:Mt,:Mt_1,:Mc,:Mc_1,:N,:N_1,:Nt,:Nt_1,:Nc,:Nc_1]]

        z[find(isna(z[:o_B0])), :o_B0] = 0.0
        z[find(isna(z[:o_SE0])), :o_SE0] = 0.0
        z[find(isna(z[:y_B0])), :y_B0] = 0.0
        z[find(isna(z[:y_SE0])), :y_SE0] = 0.0
        z[find(isna(z[:p_B0])), :p_B0] = 0.0
        z[find(isna(z[:p_SE0])), :p_SE0] = 0.0

        xdf=z

        #if !withNA
            #for c in 2:length(z)
            #    r=xdf[1,c]
	        #    xdf[find(isna(xdf[c])), c]=r
            #end
            #for i in 1:length(xdf)
            #    o_fval=mocc.reff.sdf[names(mocc.reff.sdf)[i]][1]
            #    y_fval=mdolocc.reff.sdf[names(mdolocc.reff.sdf)[i]][1]
            #    p_fval=mpen.reff.sdf[names(mpen.reff.sdf)[i]][1]
            #    println(i," ~~ ",o_fval, " - " ,y_fval, " - " ,p_fval )
            #end

            #for i in mdfcols
            #    r=
            #    xdf[find(isna(xdf[c])), c]=r
            #end

        #end
    end
    xdf = join(xdf, cnts.sdf[[:key,:class,:level,:exposed]], on = :key)
    return xdf
end
#genRndMDF(cnts, mocc, mdolocc, mpen, false)

function genFixedMDF(cnts::Cnts, mocc::MOcc, mdolocc::MDolOcc, mpen::MPen)
    o=deepcopy(mocc.feff.sdf[[:key,:B,:SE,:P,:unadj_mean_score0,:unadj_mean_score1,:M,:Mt,:Mc,:adj_mean_expsd_grp,:adj_mean_cntrl_grp]])
    rename!(o,:B,:B1)
    rename!(o,:P,:P1)
    rename!(o,:SE,:SE1)
    #rename!(o,:unadj_mean_score0,:o_mean_score0)
    #rename!(o,:unadj_mean_score1,:o_mean_score1)
    rename!(o,:adj_mean_cntrl_grp,:o_mean_score0)
    rename!(o,:adj_mean_expsd_grp,:o_mean_score1)
    o[:inOcc] = true
    o[:isSig1] = true
      o[:B1_orig]=o[:B1]
      o[:SE1_orig]=o[:SE1]

    y=deepcopy(mdolocc.feff.sdf[[:key,:B,:SE,:P,:unadj_mean_score0,:unadj_mean_score1,:adj_mean_expsd_grp,:adj_mean_cntrl_grp]])
    rename!(y,:B,:B2)
    rename!(y,:P,:P2)
    rename!(y,:SE,:SE2)
    #rename!(y,:unadj_mean_score0,:y_mean_score0)
    #rename!(y,:unadj_mean_score1,:y_mean_score1)
    rename!(y,:adj_mean_cntrl_grp,:y_mean_score0)
    rename!(y,:adj_mean_expsd_grp,:y_mean_score1)
    y[:inDolOcc] = true
    y[:isSig2] = true
       y[:B2_orig]=y[:B2]
       y[:SE2_orig]=y[:SE2]

    p=deepcopy(mpen.feff.sdf[[:key,:B,:SE,:P,:unadj_mean_score0,:unadj_mean_score1,:M,:Mt,:Mc,:adj_mean_expsd_grp,:adj_mean_cntrl_grp]])
    rename!(p,:B,:B3)
    rename!(p,:P,:P3)
    rename!(p,:SE,:SE3)
    #rename!(p,:unadj_mean_score0,:p_mean_score0)
    #rename!(p,:unadj_mean_score1,:p_mean_score1)
    rename!(p,:adj_mean_cntrl_grp,:p_mean_score0)
    rename!(p,:adj_mean_expsd_grp,:p_mean_score1)

    rename!(p,:M,:N)
    rename!(p,:Mt,:Nt)
    rename!(p,:Mc,:Nc)
    p[:inPen] = true
    p[:isSig3] = true
    p[:B3_orig]=p[:B3]
    p[:SE3_orig]=p[:SE3]

    mdfcols=[:key,:B1,:SE1,:P1,:o_mean_score0,:o_mean_score1,:inOcc,:isSig1,:B1_orig,:SE1_orig,
             :B2,:SE2,:P2,:y_mean_score0,:y_mean_score1,:inDolOcc,:isSig2,:B2_orig,:SE2_orig,
             :B3,:SE3,:P3,:p_mean_score0,:p_mean_score1,:inPen,:isSig3,:B3_orig,:SE3_orig,
             :M,:Mt,:Mc,:N,:Nt,:Nc
            ]
    xdf= join(join(o,y, on = :key),p, on = :key)[mdfcols]
    xdf = join(xdf, cnts.sdf[[:key,:class,:level]], on = :key)
    return xdf
end
#genFixedMDF(cnts, mocc, mdolocc, mpen)

function genMDF(cnts::Cnts, mocc::MOcc, mdolocc::MDolOcc, mpen::MPen)
    mdfcols=[:key,:B1,:SE1,:P1,:o_mean_score0,:o_mean_score1,:inOcc,:isSig1,:B1_orig,:SE1_orig,
             :B2,:SE2,:P2,:y_mean_score0,:y_mean_score1,:inDolOcc,:isSig2,:B2_orig,:SE2_orig,
             :B3,:SE3,:P3,:p_mean_score0,:p_mean_score1,:inPen,:isSig3,:B3_orig,:SE3_orig,
             :M,:Mt,:Mc,:N,:Nt,:Nc
            ]

    fmdf=genFixedMDF(cnts, mocc, mdolocc, mpen)
    if length(cnts.sdf[1]) > 1
        rmdf=genRndMDF(cnts, mocc, mdolocc, mpen)
        xdf = vcat(fmdf,rmdf)[mdfcols]
    else
        xdf = fmdf
    end

    #for i in [:B1_orig,:SE1_orig,:B2_orig,:SE2_orig,:B3_orig,:SE3_orig]
    #    xdf[find(isna(xdf[i])), i]=0.0
    #end

    for i in mdfcols   # default to total
        r=xdf[1,i]
        xdf[find(isna(xdf[i])), i]=r
    end

    return xdf
end
#genMDF(cnts, mocc, mdolocc, mpen)

type MDF
    sdf::DataFrame
    withNA::Bool
    mocc::MOcc
    mdolocc::MDolOcc
    mpen::MPen
    cnts::Cnts
#    get::Function
    function MDF(cnts::Cnts,mocc::MOcc,mdolocc::MDolOcc,mpen::MPen,withNA::Bool=true)
        this=new()
        this.cnts = cnts
        this.mocc = mocc
        this.mdolocc = mdolocc
        this.mpen = mpen
        this.withNA=withNA
#        function getk(k::AbstractString)
#            #if !isdefined(this, :sdf)
#            #    this.sdf = genMDF(this.cnts,this.mocc,this.mdolocc, this.mpen, false)
#            #end
#            this.sdf = genMDF(this.cnts,this.mocc,this.mdolocc, this.mpen)
#            return this.sdf[this.sdf[:key].==k,:]
#        end

#        this.get = getk
        return this
    end
end

function Base.show(io::IO, mdf::MDF)
    mdf.sdf = genMDF(mdf.cnts,mdf.mocc,mdf.mdolocc, mdf.mpen)
    show(io,mdf.sdf)
end

#function genRDF_OLD(mocc::MOcc,mdolocc::MDolOcc,mpen::MPen, mdolhh::MDolHH)
#    rdfFmt=[:key,:model,
#            :unadj_avg_expsd_hh_pre,:unadj_avg_cntrl_hh_pre,:unadj_avg_expsd_hh_pst,:unadj_avg_cntrl_hh_pst,:adj_mean_expsd_grp,:adj_mean_cntrl_grp,
#            :adj_dod_effct,:twotail_pval,:onetail_pval,
#            :onetail_80_pct_intrvl_lb,:onetail_80_pct_intrvl_ub,:twotail_80_pct_intrvl_lb,:twotail_80_pct_intrvl_ub,
#            :onetail_90_pct_intrvl_lb,:onetail_90_pct_intrvl_ub,:twotail_90_pct_intrvl_lb,:twotail_90_pct_intrvl_ub
#           ]
#    x=vcat(mocc.feff.sdf[rdfFmt], mdolocc.feff.sdf[rdfFmt],mpen.feff.sdf[rdfFmt], mdolhh.sdf[1,rdfFmt])
#    if isdefined(mocc, :reff)
#        #x=vcat(x, mocc.reff.sdf[rdfFmt], mdolocc.reff.sdf[rdfFmt], mpen.reff.sdf[rdfFmt], mdolhh.sdf[2:end,rdfFmt])
#        x=vcat(x,mocc.reff.sdf[map(x -> contains(x, "(none)") ? 1:0, mocc.reff.sdf[:key]) .== 0,:][rdfFmt],
#                 mdolocc.reff.sdf    [map(x -> contains(x, "(none)") ? 1:0, mdolocc.reff.sdf[:key]) .== 0,:]     [rdfFmt],
#                 mpen.reff.sdf       [map(x -> contains(x, "(none)") ? 1:0, mpen.reff.sdf[:key]) .== 0,:]     [rdfFmt],
#                 mdolhh.sdf          [map(x -> contains(x, "(none)") ? 1:0, mdolhh.sdf[:key]) .== 0,:]     [2:end,rdfFmt]
#              )
#    end
#    sdf=sort(x, cols = :key)
#end
##genRDF(mocc,mdolocc,mpen, mdolhh)

function genRDF(mocc::MOcc,mdolocc::MDolOcc,mpen::MPen, mdolhh::MDolHH,cfg::OrderedDict=OrderedDict(:TotalModelsOnly => false))
    rdfFmt=[:key,:model,
            :unadj_avg_expsd_hh_pre,:unadj_avg_cntrl_hh_pre,:unadj_avg_expsd_hh_pst,:unadj_avg_cntrl_hh_pst,:adj_mean_expsd_grp,:adj_mean_cntrl_grp,
            :adj_dod_effct,:twotail_pval,:onetail_pval,
            :onetail_80_pct_intrvl_lb,:onetail_80_pct_intrvl_ub,:twotail_80_pct_intrvl_lb,:twotail_80_pct_intrvl_ub,
            :onetail_90_pct_intrvl_lb,:onetail_90_pct_intrvl_ub,:twotail_90_pct_intrvl_lb,:twotail_90_pct_intrvl_ub
           ]
    x=vcat(mocc.feff.sdf[rdfFmt], mdolocc.feff.sdf[rdfFmt],mpen.feff.sdf[rdfFmt], mdolhh.sdf[1,rdfFmt])

    if !cfg[:TotalModelsOnly]
    #if isdefined(mocc, :reff)
        #x=vcat(x,mocc.reff.sdf[map(x -> contains(x, "(none)") ? 1:0, mocc.reff.sdf[:key]) .== 0,:][rdfFmt],
        #         mdolocc.reff.sdf[map(x -> contains(x, "(none)") ? 1:0, mdolocc.reff.sdf[:key]) .== 0,:][rdfFmt],
        #         mpen.reff.sdf[map(x -> contains(x, "(none)") ? 1:0, mpen.reff.sdf[:key]) .== 0,:][rdfFmt],
        #         mdolhh.sdf[map(x -> contains(x, "(none)") ? 1:0, mdolhh.sdf[:key]) .== 0,:][2:end,rdfFmt]
        #      )
        z=vcat(mocc.reff.sdf[rdfFmt], mdolocc.reff.sdf[rdfFmt], mpen.reff.sdf[rdfFmt],mdolhh.sdf[2:end,rdfFmt])
        if length(z[1]) > 0  x=vcat(x,z) end
    end
    x2 = similar(x,0)
    for k in unique(x[:key])
        for m in ["occ","dolocc","pen","dolhh"]
            len = length(x[ (x[:model].== m) & (x[:key].== k)  ,:key])
            if len == 1
                v = collect(values(df2dict(x[(x[:model].== m) & (x[:key].== k),:])))
                push!(x2, v)
            else
                push!(x2, [k;m;
                           x[(x[:key].=="Total Campaign") & (x[:model].==m),:unadj_avg_expsd_hh_pre];
                           x[(x[:key].=="Total Campaign") & (x[:model].==m),:unadj_avg_cntrl_hh_pre];
                           x[(x[:key].=="Total Campaign") & (x[:model].==m),:unadj_avg_expsd_hh_pst];
                           x[(x[:key].=="Total Campaign") & (x[:model].==m),:unadj_avg_cntrl_hh_pst];
                           x[(x[:key].=="Total Campaign") & (x[:model].==m),:adj_mean_expsd_grp];
                           x[(x[:key].=="Total Campaign") & (x[:model].==m),:adj_mean_cntrl_grp];
                           x[(x[:key].=="Total Campaign") & (x[:model].==m),:adj_dod_effct];
                           x[(x[:key].=="Total Campaign") & (x[:model].==m),:twotail_pval];
                           x[(x[:key].=="Total Campaign") & (x[:model].==m),:onetail_pval];
                           NaN;  #:onetail_80_pct_intrvl_lb
                           NaN;  #:onetail_80_pct_intrvl_ub
                           NaN;  #:twotail_80_pct_intrvl_lb
                           NaN;  #:twotail_80_pct_intrvl_ub
                           NaN;  #:onetail_90_pct_intrvl_lb
                           NaN;  #:onetail_90_pct_intrvl_ub
                           NaN;  #:twotail_90_pct_intrvl_lb
                           NaN   #:twotail_90_pct_intrvl_ub
                          ]
                    )
            end
        end
    end
    sdf=sort(x2, cols = :key)
end

type RDF
    sdf::DataFrame
    mocc::MOcc
    mdolocc::MDolOcc
    mpen::MPen
    mdolhh::MDolHH
    function RDF(mocc::MOcc,mdolocc::MDolOcc,mpen::MPen,mdolhh::MDolHH)
        this=new()
        this.mocc = mocc
        this.mdolocc = mdolocc
        this.mpen = mpen
        this.mdolhh = mdolhh
        return this
    end
end

function Base.show(io::IO, rdf::RDF)
    #x=vcat(rdf.mocc.feff.sdf[rdfFmt], rdf.mdolocc.feff.sdf[rdfFmt],rdf.mpen.feff.sdf[rdfFmt],rdf.mdolhh.sdf[1,rdfFmt])
    #if isdefined(rdf.mocc, :reff)
    #    x=vcat(x, rdf.mocc.reff.sdf[rdfFmt], rdf.mdolocc.reff.sdf[rdfFmt], rdf.mpen.reff.sdf[rdfFmt], rdf.mdolhh.sdf[2:end,rdfFmt])
    #end
    #rdf.sdf=sort(x, cols = :key)
    rdf.sdf=genRDF(rdf.mocc,rdf.mdolocc,rdf.mpen, rdf.mdolhh)
    show(io,rdf.sdf)
end

"""
# scoreme file
function aggregate!(df_data::DataFrame, focc::FOcc, cnts::Cnts)
    sdf = focc.sdf
    sdf[:M] = cnts.get("Total Campaign")[:M][1]
    sdf[:Mt] = cnts.get("Total Campaign")[:Mt][1]
    sdf[:Mc] = cnts.get("Total Campaign")[:Mc][1]
#    if :reff_occ in names(df_data)
#        println("using reff")
#        df_data[:occ_score0] = exp(df_data[:pre_occ_score0] + df_data[:reff_occ]  )
#        df_data[:occ_score1] = exp(df_data[:pre_occ_score1] + df_data[:reff_occ] )
#    else
#        df_data[:occ_score0] = exp(df_data[:pre_occ_score0])
#        df_data[:occ_score1] = exp(df_data[:pre_occ_score1])
#    end
    sdf[:mean_score0]=mean( df_data[ df_data[:buyer_pos_p1] .== 1 , :occ_score0] )
    sdf[:mean_score1]=mean( df_data[ df_data[:buyer_pos_p1] .== 1 , :occ_score1] )

    sdf[:unadj_mean_score0]=mean(df_data[(df_data[:buyer_pos_p1] .== 1) & (df_data[:group] .== 0),:occ_score0])
    sdf[:unadj_mean_score1]=mean(df_data[(df_data[:buyer_pos_p1] .== 1) & (df_data[:group] .== 1),:occ_score1])

    sdf[:adj_mean_cntrl_grp ] = mean( df_data[ df_data[:buyer_pos_p1] .== 1 , :occ_score0] )
    sdf[:adj_mean_expsd_grp ] = mean( df_data[ df_data[:buyer_pos_p1] .== 1 , :occ_score1] )
    # -------  Occ Score
    sdf[:unadj_avg_cntrl_hh_pre] = mean(df_data[ (df_data[:buyer_pre_p1] .== 1) & (df_data[:group] .== 0), :trps_pre_p1] )
    sdf[:unadj_avg_expsd_hh_pre] = mean(df_data[ (df_data[:buyer_pre_p1] .== 1) & (df_data[:group] .== 1), :trps_pre_p1] )
    sdf[:unadj_avg_cntrl_hh_pst] = mean(df_data[ (df_data[:buyer_pos_p1] .== 1) & (df_data[:group] .== 0), :trps_pos_p1] )
    sdf[:unadj_avg_expsd_hh_pst] = mean(df_data[ (df_data[:buyer_pos_p1] .== 1) & (df_data[:group] .== 1), :trps_pos_p1] )
end

function aggregate!(df_data::DataFrame, fdolocc::FDolOcc, cnts::Cnts)
    sdf = fdolocc.sdf
    sdf[:M] = cnts.get("Total Campaign")[:M][1]
    sdf[:Mt] = cnts.get("Total Campaign")[:Mt][1]
    sdf[:Mc] = cnts.get("Total Campaign")[:Mc][1]
#    if :reff_dolocc in names(df_data)
#        println("using reff")
#        df_data[:dolocc_score0] = exp(df_data[:pre_dolocc_score0] + df_data[:reff_dolocc]  )
#        df_data[:dolocc_score1] = exp(df_data[:pre_dolocc_score1] + df_data[:reff_dolocc] )
#    else
#        df_data[:dolocc_score0] = exp(df_data[:pre_dolocc_score0])
#        df_data[:dolocc_score1] = exp(df_data[:pre_dolocc_score1])
#    end
    sdf[:mean_score0]=mean( df_data[ df_data[:buyer_pos_p1] .== 1 , :dolocc_score0] )
    sdf[:mean_score1]=mean( df_data[ df_data[:buyer_pos_p1] .== 1 ,  :dolocc_score1] )

    sdf[:unadj_mean_score0]=mean(df_data[(df_data[:buyer_pos_p1] .== 1) & (df_data[:group] .== 0),:dolocc_score0])
    sdf[:unadj_mean_score1]=mean(df_data[(df_data[:buyer_pos_p1] .== 1) & (df_data[:group] .== 1),:dolocc_score1])

    sdf[:adj_mean_cntrl_grp ] =  mean( df_data[ df_data[:buyer_pos_p1] .== 1 , :dolocc_score0] )
    sdf[:adj_mean_expsd_grp ] =  mean( df_data[ df_data[:buyer_pos_p1] .== 1 ,  :dolocc_score1] )
    # -------  DolOcc Score
    sdf[:unadj_avg_cntrl_hh_pre] =  mean(df_data[ (df_data[:buyer_pre_p1] .== 1) & (df_data[:group] .== 0), :dol_per_trip_pre_p1] )
    sdf[:unadj_avg_expsd_hh_pre] =  mean(df_data[ (df_data[:buyer_pre_p1] .== 1) & (df_data[:group] .== 1), :dol_per_trip_pre_p1] )
    sdf[:unadj_avg_cntrl_hh_pst] =  mean(df_data[ (df_data[:buyer_pos_p1] .== 1) & (df_data[:group] .== 0), :dol_per_trip_pos_p1] )
    sdf[:unadj_avg_expsd_hh_pst] =  mean(df_data[ (df_data[:buyer_pos_p1] .== 1) & (df_data[:group] .== 1), :dol_per_trip_pos_p1] )
end

function aggregate!(df_data::DataFrame, fpen::FPen, cnts::Cnts)
    sdf = fpen.sdf
    sdf[:M] = cnts.get("Total Campaign")[:N][1]
    sdf[:Mt] = cnts.get("Total Campaign")[:Nt][1]
    sdf[:Mc] = cnts.get("Total Campaign")[:Nc][1]
#    if :reff_pen in names(df_data)
#        println("using reff")
#        df_data[:pen_score0] = exp((df_data[:pre_pen_score0] + df_data[:reff_pen])) ./ ( exp((df_data[:pre_pen_score0] + df_data[:reff_pen])) +1)
#        df_data[:pen_score1] = exp((df_data[:pre_pen_score1] + df_data[:reff_pen])) ./ (exp((df_data[:pre_pen_score1] + df_data[:reff_pen])) +1)
#    else
#        df_data[:pen_score0] =  exp(df_data[:pre_pen_score0]) ./ ( exp(df_data[:pre_pen_score0]) +1)
#        df_data[:pen_score1] = exp(df_data[:pre_pen_score1]) ./ ( exp(df_data[:pre_pen_score1]) +1)
#    end
    sdf[:mean_score0]=mean( df_data[ :pen_score0] )
    sdf[:mean_score1]=mean( df_data[ :pen_score1] )
    sdf[:unadj_mean_score0]=mean(df_data[(df_data[:group] .== 0),:pen_score0])
    sdf[:unadj_mean_score1]=mean(df_data[(df_data[:group] .== 1),:pen_score1])

    sdf[:adj_mean_cntrl_grp ] = mean( df_data[ :pen_score0] )
    sdf[:adj_mean_expsd_grp ] = mean( df_data[ :pen_score1] )
    # -------  Pen Score   - buyer_pre_p1
    sdf[:unadj_avg_cntrl_hh_pre ] = mean(df_data[ (df_data[:group] .== 0)  , :buyer_pre_p1] )
    sdf[:unadj_avg_expsd_hh_pre ] = mean(df_data[ (df_data[:group] .== 1)  , :buyer_pre_p1] )
    sdf[:unadj_avg_cntrl_hh_pst ] = mean(df_data[ (df_data[:group] .== 0)  , :buyer_pos_p1] )
    sdf[:unadj_avg_expsd_hh_pst ] = mean(df_data[ (df_data[:group] .== 1)  , :buyer_pos_p1] )
end

function aggregateCommon!(fixedeff::FixedEffect)
    sdf = fixedeff.sdf
    #src = fixedeff.src
    sdf[:adj_dod_effct] = ((sdf[:adj_mean_expsd_grp] - sdf[:adj_mean_cntrl_grp])  ./ sdf[:adj_mean_cntrl_grp] ) * 100
    sdf[:twotail_pval] = sdf[:P]   #twotail_pval
    sdf[:onetail_pval] = 1 - (sdf[:twotail_pval] ./ 2)   #onetail_pval
    sdf[:twotail_pval] = 1-sdf[:P]    #redo :twotail_pval
    #For NON-SIG Models, update - adj_mean_expsd_grp (to = adj_mean_cntrl_grp) & adj_dod_effct (= 0)
    sdf[(sdf[:twotail_pval] .< 0.80) & (sdf[:twotail_pval] .!= 0),:adj_mean_expsd_grp]      =  sdf[(sdf[:twotail_pval] .< 0.80) & (sdf[:twotail_pval] .!= 0),:adj_mean_cntrl_grp]
end

function ConfidenceIntervals(feff::FixedEffect)
    sdf = feff.sdf
    for row in eachrow(sdf)
            k = row[:key]
            md = df2dict(row)
            md[:mean_score0] = md[:unadj_mean_score0];
            md[:mean_score1] = md[:unadj_mean_score1];
            cis =  CIs(md, ZDict, 0)
            for (zscore_key,zscore) in  ZDict
                    ubk = Symbol(zscore_key*"_ub")
                    lbk = Symbol(zscore_key*"_lb")
                    sdf[sdf[:key].==k, lbk] = md[lbk]*100
                    sdf[sdf[:key].==k, ubk] = md[ubk]*100
            end
    end
end

# === Random Effects ====
function init!(cnts::Cnts, m::MModel)
        src=m.reff.src
        # --- GROUP1 -----
        src[:Z1] = src[:B1] ./ src[:SE1]
        src[:P1] = 2.0 * ccdf(Normal(), abs(src[:Z1]))
        src[:SIG1] = map( x -> ((isnan(x))|(x > 0.2)) ?  false : true, src[:P1] )
        # --- GROUP0 -----
        src[:Z0] = src[:B0] ./ src[:SE0]
        src[:P0] = 2.0 * ccdf(Normal(), abs(src[:Z0]))
        src[:SIG0] = map( x -> (isnan(x)|(x > 0.2)) ?  false : true, src[:P0] )

        m.reff.sdf = join(m.reff.sdf,m.reff.src[[:class,:level,:B0,:B1,:SE0,:SE1,:exposed,:P1,:key]], on = :key)
        sdf=m.reff.sdf
        sdf[:B_fixed] = m.feff.sdf[1,:B]
        sdf[:P_fixed] = m.feff.sdf[1,:P]
        sdf[:SE_fixed] = m.feff.sdf[1,:SE]
        #sdf[:B1_combo] = sdf[:B1] + sdf[:B_fixed] #- sdf[:B0]
        #sdf[:SE1_combo] = sqrt(sdf[:SE_fixed].^2+sdf[:SE1].^2)
        sdf[:B1_combo] = sdf[:B1] + sdf[:B_fixed] - sdf[:B0]
        sdf[:SE1_combo] = sqrt(sdf[:SE_fixed].^2+sdf[:SE1].^2+sdf[:SE0].^2)
        sdf[:P1_combo] = 2.0 * ccdf(Normal(), abs( (sdf[:B1_combo]) ./ sdf[:SE1_combo]))

        if m.reff.v_model == "pen"
            for k in sdf[:key]
                sdf[sdf[:key].==k ,:M] = cnts.get(k)[:N][1]
                sdf[sdf[:key].==k ,:Mt] = cnts.get(k)[:Nt][1]
                sdf[sdf[:key].==k ,:Mc] = cnts.get(k)[:Nc][1]
            end
        else
            for k in sdf[:key]
                sdf[sdf[:key].==k ,:M] = cnts.get(k)[:M][1]
                sdf[sdf[:key].==k ,:Mt] = cnts.get(k)[:Mt][1]
                sdf[sdf[:key].==k ,:Mc] = cnts.get(k)[:Mc][1]
            end
        end
end

function aggregate!(df_data::DataFrame, mocc::MOcc)
    sdf=mocc.reff.sdf
    for row = eachrow(mocc.reff.sdf)
        k = row[:key]
        ranfx = lowercase(row[:class])
        sranfx = Symbol(ranfx)
        v_level = row[:level]
        exposed = row[:exposed]

        dest_colname0=Symbol(k*"_occ_0")
        dest_colname1=Symbol(k*"_occ_1")
        xcol=Symbol("reff_"*ranfx*"_occ")   # Column to exclude
        fmula0="df_data[:pre_occ_score0]+"*getRandomFormula(mocc.reff.rcols, xcol)*"+"*string(row[:B0])
        fmula0=replace(fmula0,"++","+")    #if only one col
        fmula1="df_data[:pre_occ_score1]+"*getRandomFormula(mocc.reff.rcols, xcol)*"+"*string(row[:B1])
        fmula1=replace(fmula1,"++","+")

        df_data[dest_colname0] = exp(eval(parse(fmula0)))
        df_data[dest_colname1] = exp(eval(parse(fmula1)))
        println("\n\n\n",xcol," ~~ ",fmula1)
        println("occ ~ ",k," --> ",dest_colname0)

        mean_score0 = mean(df_data[df_data[:buyer_pos_p1] .== 1,  dest_colname0 ])
        mean_score1 = mean(df_data[ df_data[:buyer_pos_p1] .== 1 , dest_colname1 ])
        sdf[sdf[:key].==k ,:mean_score0] = mean_score0
        sdf[sdf[:key].==k ,:mean_score1] = mean_score1
        sdf[sdf[:key].==k ,:adj_mean_cntrl_grp] = mean_score0
        sdf[sdf[:key].==k ,:adj_mean_expsd_grp] = mean_score1

        #sdf[sdf[:key].==k ,:unadj_avg_expsd_hh_pre] = mean(df_data[ (df_data[:group] .== 1) & (df_data[:buyer_pre_p1] .== 1) & (df_data[sranfx] .== v_level) , :trps_pre_p1] )
        #sdf[sdf[:key].==k ,:unadj_avg_expsd_hh_pst] = mean(df_data[ (df_data[:group] .== 1) & (df_data[:buyer_pos_p1] .== 1) & (df_data[sranfx] .== v_level) , :trps_pos_p1] )
        sdf[sdf[:key].==k ,:unadj_avg_expsd_hh_pre] = mean(df_data[ (df_data[:group] .== 1) & (df_data[sranfx] .== v_level) , :trps_pre_p1] )
        sdf[sdf[:key].==k ,:unadj_avg_expsd_hh_pst] = mean(df_data[ (df_data[:group] .== 1) & (df_data[sranfx] .== v_level) , :trps_pos_p1] )

        if exposed
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pre] = mocc.feff.sdf[1,:unadj_avg_cntrl_hh_pre]
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pst] = mocc.feff.sdf[1,:unadj_avg_cntrl_hh_pst]
        else
            #sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pre] = mean(df_data[ (df_data[:group] .== 0) & (df_data[:buyer_pre_p1] .== 1) & (df_data[sranfx] .== v_level) , :trps_pre_p1] )
            #sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pst] = mean(df_data[ (df_data[:group] .== 0) & (df_data[:buyer_pos_p1] .== 1) & (df_data[sranfx] .== v_level) , :trps_pos_p1] )
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pre] = mean(df_data[ (df_data[:group] .== 0) & (df_data[sranfx] .== v_level) , :trps_pre_p1] )
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pst] = mean(df_data[ (df_data[:group] .== 0) & (df_data[sranfx] .== v_level) , :trps_pos_p1] )
        end
    end
end

function aggregate!(df_data::DataFrame, mdolocc::MDolOcc)
    sdf=mdolocc.reff.sdf
    for row = eachrow(mdolocc.reff.sdf)
        k = row[:key]
        ranfx = lowercase(row[:class])
        sranfx = Symbol(ranfx)
        v_level = row[:level]
        exposed = row[:exposed]
        dest_colname0=Symbol(k*"_dolocc_0")
        dest_colname1=Symbol(k*"_dolocc_1")
        xcol=Symbol("reff_"*ranfx*"_dolocc")
        fmula0="df_data[:pre_dolocc_score0]+"*getRandomFormula(mdolocc.reff.rcols, xcol)*"+"*string(row[:B0])
        fmula0=replace(fmula0,"++","+")    #if only one col
        fmula1="df_data[:pre_dolocc_score1]+"*getRandomFormula(mdolocc.reff.rcols, xcol)*"+"*string(row[:B1])
        fmula1=replace(fmula1,"++","+")

        df_data[dest_colname0] = exp(eval(parse(fmula0)))
        df_data[dest_colname1] = exp(eval(parse(fmula1)))
        println("\n\n\n",xcol," ~~ ",fmula1)
        println("dolocc ~ ",k," --> ",dest_colname0)

        mean_score0 = mean(df_data[df_data[:buyer_pos_p1] .== 1,  dest_colname0 ])
        mean_score1 = mean(df_data[ df_data[:buyer_pos_p1] .== 1 , dest_colname1 ])
        sdf[sdf[:key].==k ,:mean_score0] = mean_score0
        sdf[sdf[:key].==k ,:mean_score1] = mean_score1
        sdf[sdf[:key].==k ,:adj_mean_cntrl_grp] = mean_score0
        sdf[sdf[:key].==k ,:adj_mean_expsd_grp] = mean_score1

        #sdf[sdf[:key].==k ,:unadj_avg_expsd_hh_pre] = mean(df_data[ (df_data[:group] .== 1) & (df_data[:buyer_pre_p1] .== 1) & (df_data[sranfx] .== v_level) , :dol_per_trip_pre_p1] )
        #sdf[sdf[:key].==k ,:unadj_avg_expsd_hh_pst] = mean(df_data[ (df_data[:group] .== 1) & (df_data[:buyer_pos_p1] .== 1) & (df_data[sranfx] .== v_level) , :dol_per_trip_pos_p1] )
        sdf[sdf[:key].==k ,:unadj_avg_expsd_hh_pre] = mean(df_data[ (df_data[:group] .== 1) & (df_data[sranfx] .== v_level) , :dol_per_trip_pre_p1] )
        sdf[sdf[:key].==k ,:unadj_avg_expsd_hh_pst] = mean(df_data[ (df_data[:group] .== 1) & (df_data[sranfx] .== v_level) , :dol_per_trip_pos_p1] )

        if exposed
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pre] = mdolocc.feff.sdf[1,:unadj_avg_cntrl_hh_pre]
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pst] = mdolocc.feff.sdf[1,:unadj_avg_cntrl_hh_pst]
        else
            #sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pre] = mean(df_data[ (df_data[:group] .== 0) & (df_data[:buyer_pre_p1] .== 1) & (df_data[sranfx] .== v_level) , :dol_per_trip_pre_p1] )
            #sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pst] = mean(df_data[ (df_data[:group] .== 0) & (df_data[:buyer_pos_p1] .== 1) & (df_data[sranfx] .== v_level) , :dol_per_trip_pos_p1] )
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pre] = mean(df_data[ (df_data[:group] .== 0) & (df_data[sranfx] .== v_level) , :dol_per_trip_pre_p1] )
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pst] = mean(df_data[ (df_data[:group] .== 0) & (df_data[sranfx] .== v_level) , :dol_per_trip_pos_p1] )
        end
    end
end

function aggregate!(df_data::DataFrame,mpen::MPen)
    sdf=mpen.reff.sdf
    for row = eachrow(mpen.reff.sdf)
        k = row[:key]
        ranfx = lowercase(row[:class])
        sranfx = Symbol(ranfx)
        v_level = row[:level]
        exposed = row[:exposed]
        dest_colname0=Symbol(k*"_pen_0")
        dest_colname1=Symbol(k*"_pen_1")
        xcol=Symbol("reff_"*ranfx*"_pen")
        fmula0="df_data[:pre_pen_score0]+"*getRandomFormula(mpen.reff.rcols, xcol)*"+"*string(row[:B0])
        fmula0=replace(fmula0,"++","+")    #if only one col
        fmula1="df_data[:pre_pen_score1]+"*getRandomFormula(mpen.reff.rcols, xcol)*"+"*string(row[:B1])
        fmula1=replace(fmula1,"++","+")

        df_data[dest_colname0] = exp(eval(parse(fmula0))) ./ (1+exp(eval(parse(fmula0))))
        df_data[dest_colname1] = exp(eval(parse(fmula1))) ./ (1+exp(eval(parse(fmula1))))
        println("\n\n\n",xcol," ~~ ",fmula1)
        println("pen ~ ",k," --> ",dest_colname1)
        mean_score0 = mean(df_data[dest_colname0])
        mean_score1 = mean(df_data[dest_colname1])

        println("..... ",mean_score0," ~ ",mean_score1)

        sdf[sdf[:key].==k ,:mean_score0] = mean_score0
        sdf[sdf[:key].==k ,:mean_score1] = mean_score1
        sdf[sdf[:key].==k ,:adj_mean_cntrl_grp] = mean_score0
        sdf[sdf[:key].==k ,:adj_mean_expsd_grp] = mean_score1

        sdf[sdf[:key].==k ,:unadj_avg_expsd_hh_pre] = mean(df_data[ (df_data[:group] .== 1) & (df_data[Symbol(ranfx)] .== v_level) , :buyer_pre_p1] )
        sdf[sdf[:key].==k ,:unadj_avg_expsd_hh_pst] = mean(df_data[ (df_data[:group] .== 1) & (df_data[Symbol(ranfx)] .== v_level) , :buyer_pos_p1] )
        if exposed
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pre] = mpen.feff.sdf[1,:unadj_avg_cntrl_hh_pre]
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pst] = mpen.feff.sdf[1,:unadj_avg_cntrl_hh_pst]
        else
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pre] = mean(df_data[ (df_data[:group] .== 0) & (df_data[Symbol(ranfx)] .== v_level) , :buyer_pre_p1] )
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pst] = mean(df_data[ (df_data[:group] .== 0) & (df_data[Symbol(ranfx)] .== v_level) , :buyer_pos_p1] )
        end
    end
end
"""

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
using Random, Distributions
Random.seed!(0)

lambda = 0.1
dist = Exponential(1/lambda)
actualMedian = median(dist)

M = 10^3
N = 10^4
nRange = 2:2:10
alpha = 0.05

for n in nRange
    coverageCount = 0
    for _ in 1:M
        sampleData = rand(dist, n)
        bootstrapSampleMeans = [median(rand(sampleData, n)) for _ in 1:N]
        L = quantile(bootstrapSampleMeans, alpha/2)
        U = quantile(bootstrapSampleMeans, 1-alpha/2)
        coverageCount += L < actualMedian && actualMedian < U
    end
    println("n = ",n,"\t coverage = ", coverageCount/M)
end

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using YaoPlots
using Compose
using Test
using YaoBlocks

@testset "gate styles" begin
	@test YaoPlots.get_brush_texts(X)[][2] == "X"
	@test YaoPlots.get_brush_texts(Rx(0.5))[][2] == "Rx(0.5)"
	@test YaoPlots.get_brush_texts(shift(0.5))[][2] == "ϕ(0.5)"
	@test YaoPlots.get_brush_texts(YaoBlocks.phase(0.5))[][2] == "^0.5"
end

@testset "circuit canvas" begin
	c = CircuitGrid(5)
	@test YaoPlots.nline(c) == 5
	@test YaoPlots.frontier(c, 2, 3) == 0
	@test YaoPlots.depth(c) == 0
	circuit_canvas(5) do c
		YaoPlots.draw!(c, put(5, 3=>X), 1:5, [])
		@test YaoPlots.frontier(c, 1, 2) == 0
		@test YaoPlots.frontier(c, 3, 5) == 1
		@test YaoPlots.depth(c) == 1
	end

	gg = circuit_canvas(5) do c
		put(3=>X) >> c
		control(2, 3=>X) >> c
		control(2, 3=>shift(0.5)) >> c
		chain(5, control(2, 3=>X), put(1=>X)) >> c
		@test YaoPlots.depth(c) == 4
	end
	@test gg isa Context

	g = put(5, (3, 4)=>SWAP) |> vizcircuit(; scale=0.7, w_line=0.8, w_depth=0.9)
	@test g isa Context
	@test vizcircuit(put(5, (3,4)=>kron(X, Y)); scale=0.7, w_line=0.8, w_depth=0.9) isa Context

	@test vizcircuit(control(10, (2, -3), 6=>X)) isa Context
	@test vizcircuit(control(10, (2, -3), 6=>im*X)) isa Context
	@test plot(put(7, (2,3)=>matblock(randn(4,4)))) isa Context
end

@testset "fix #3" begin
	@test (control(4, 2, (1,3)=>kron(X, X)) |> vizcircuit) isa Context
end

@testset "pretty_angle" begin
	@test YaoPlots.pretty_angle(π) == "π"
	@test YaoPlots.pretty_angle(π*1.0) == "π"
	@test YaoPlots.pretty_angle(-π*1.0) == "-π"
	@test YaoPlots.pretty_angle(-π*0.0) == "0"
	@test YaoPlots.pretty_angle(-π*0.5) == "-π/2"
	@test YaoPlots.pretty_angle(1.411110) == "1.41"
end


text:
@testset "Generate Types" begin
    menu = """
    {
      "menu": {
        "header": "SVG Viewer",
        "items": [
          {
            "id": "Open"
          },
          {
            "id": "OpenNew",
            "label": "Open New"
          },
          {
            "id": "ZoomIn",
            "label": "Zoom In"
          },
          {
            "id": "ZoomOut",
            "label": "Zoom Out"
          },
          {
            "id": "OriginalView",
            "label": "Original View"
          },
          {
            "id": "Quality"
          },
          {
            "id": "Pause"
          },
          {
            "id": "Mute"
          },
          null,
          {
            "id": "Find",
            "label": "Find..."
          },
          {
            "id": "FindAgain",
            "label": "Find Again"
          },
          {
            "id": "Copy"
          },
          {
            "id": "CopyAgain",
            "label": "Copy Again"
          },
          {
            "id": "CopySVG",
            "label": "Copy SVG"
          },
          {
            "id": "ViewSVG",
            "label": "View SVG"
          },
          {
            "id": "ViewSource",
            "label": "View Source"
          },
          {
            "id": "SaveAs",
            "label": "Save As"
          },
          null,
          {
            "id": "Help"
          },
          {
            "id": "About",
            "label": "About Adobe SVG Viewer..."
          }
        ]
      }
    }
    """

    menu_array = """
    [
        {
          "menu": {
            "header": "SVG Viewer",
            "items": [
              {
                "id": "Open"
              },
              {
                "id": "OpenNew",
                "label": "Open New"
              },
              {
                "id": "ZoomIn",
                "label": "Zoom In"
              },
              {
                "id": "ZoomOut",
                "label": "Zoom Out"
              },
              {
                "id": "OriginalView",
                "label": "Original View"
              },
              {
                "id": "Quality"
              },
              {
                "id": "Pause"
              },
              {
                "id": "Mute"
              },
              null,
              {
                "id": "Find",
                "label": "Find..."
              },
              {
                "id": "FindAgain",
                "label": "Find Again"
              },
              {
                "id": "Copy"
              },
              {
                "id": "CopyAgain",
                "label": "Copy Again"
              },
              {
                "id": "CopySVG",
                "label": "Copy SVG"
              },
              {
                "id": "ViewSVG",
                "label": "View SVG"
              },
              {
                "id": "ViewSource",
                "label": "View Source"
              },
              {
                "id": "SaveAs",
                "label": "Save As"
              },
              null,
              {
                "id": "Help"
              },
              {
                "id": "About",
                "label": "About Adobe SVG Viewer..."
              }
            ]
          }
        },
        {
          "menu": {
            "header": "A second item",
            "items": [
              {
                "id": "Open"
              },
              {
                "id": "OpenNew",
                "label": "Open New"
              }
            ]
          }
        }
    ]
    """

    @testset "Object & Array" begin
        json = JSON3.read(menu) # JSON.Object

        # build a type for the JSON
        raw_json_type = JSON3.generate_type(json)
        @test raw_json_type <: NamedTuple

        # turn the type into struct expressions, including replacing sub types with references to a struct
        json_exprs = JSON3.generate_exprs(raw_json_type; root_name=:MyStruct)
        @test length(json_exprs) == 3

        # write the types to a file, then can be edited/included as needed
        path = mktempdir()
        file = joinpath(path, "test_type.jl")
        JSON3.write_exprs(json_exprs, file)
        lines = readlines(file)
        @test length(lines) > 3

        json_arr = JSON3.read(menu_array) # JSON.Array

        # build a type for the JSON
        raw_json_arr_type = JSON3.generate_type(json_arr)
        @test raw_json_arr_type <: Array

        # turn the type into struct expressions, including replacing sub types with references to a struct
        json_arr_exprs = JSON3.generate_exprs(raw_json_arr_type; root_name=:MyStruct)
        @test length(json_arr_exprs) == 3

        # write the types to a file, then can be edited/included as needed
        path = mktempdir()
        file_arr = joinpath(path, "test_arr.jl")
        JSON3.write_exprs(json_arr_exprs, file_arr)
        arr_lines = readlines(file_arr)
        @test length(arr_lines) > 3
        @test arr_lines == lines
    end

    @testset "Write & Read" begin
        path = mktempdir()
        file_path_default = joinpath(path, "default.jl")
        JSON3.writetypes(menu, file_path_default)

        include(file_path_default)
        json = JSON3.read(menu, JSONTypes.Root)

        @test json.menu.header == "SVG Viewer"
        @test length(json.menu.items) > 0
        @test json.menu.items[1].id == "Open"

        @test fieldnames(JSONTypes.Item) == (:id, :label)

        file_path_json = joinpath(path, "menu_array.json")
        Base.write(file_path_json, menu_array)
        file_path_mod = joinpath(path, "my_mod.jl")
        JSON3.writetypes(file_path_json, file_path_mod; module_name=:MyMod, root_name=:MenuArray)

        include(file_path_mod)
        json_arr = JSON3.read(menu_array, Vector{MyMod.MenuArray})
        @test length(json_arr) == 2
        @test json_arr[1].menu.header == "SVG Viewer"
    end

    @testset "Struct" begin
        json = """
            [
                {"a": 1, "b": 2, "c": {"d": 4}},
                {"a": "w", "b": 5, "c": {"d": 4}},
                {"a": 3, "b": 4, "c": {"d": 6}},
                {"a": 7, "b": 7, "c": {"d": 7}}
            ]
        """

        path = mktempdir()
        file_path = joinpath(path, "struct.jl")

        JSON3.writetypes(json, file_path; mutable=false)
        include(file_path)
        parsed = JSON3.read(json, Vector{JSONTypes.Root})

        @test !(JSONTypes.Root.mutable)
        @test parsed[1].c.d == 4
        @test fieldtype(JSONTypes.Root, 1) == Union{Int64, String}
    end

    @testset "Raw Types" begin
        json = """
            [
                {"a": 1, "b": 2, "c": 5},
                {"a": "w", "b": 5, "c": 9},
                {"a": 3, "b": 4, "c": 2},
                {"a": 7, "b": 7, "c": 0}
            ]
        """
        parsed = JSON3.read(json) # JSON.Array

        # build a type for the JSON
        raw_json_type = JSON3.generate_type(parsed)
        @test raw_json_type <: Array

        res = JSON3.read(json, raw_json_type)
        @test isa(res, raw_json_type)

        @test res[4].b == 7

        json = """{"a": 1, "b": 2, "c": 5}"""
        parsed = JSON3.read(json) # JSON.Object

        # build a type for the JSON
        raw_json_type = JSON3.generate_type(parsed)
        @test raw_json_type <: NamedTuple

        res = JSON3.read(json, raw_json_type)
        @test isa(res, raw_json_type)

        @test res.b == 2
    end

    @testset "Pascal Case" begin
        @test JSON3.pascalcase(:items) == :Item
        @test JSON3.pascalcase(:i) == :I
        @test JSON3.pascalcase(:is) == :I
        @test JSON3.pascalcase(:daisies) == :Daisie # don't let perfect be the enemy of good?
        @test JSON3.pascalcase(:ties) == :Tie
        @test JSON3.pascalcase(:my_types) == :MyType
        @test JSON3.pascalcase(:witness) == :Witness
    end

    @testset "Macro" begin
        JSON3.@generatetypes(menu)
        json = JSON3.read(menu, JSONTypes.Root)

        @test json.menu.header == "SVG Viewer"
        @test length(json.menu.items) > 0
        @test json.menu.items[1].id == "Open"

        @test fieldnames(JSONTypes.Item) == (:id, :label)

        JSON3.@generatetypes(menu_array, :MyMod)
        json_arr = JSON3.read(menu_array, Vector{MyMod.Root})
        @test length(json_arr) == 2
        @test json_arr[1].menu.header == "SVG Viewer"
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder
name = "libass"
version = v"0.15.1"

# Collection of sources required to build libass
sources = [
    ArchiveSource("https://github.com/libass/libass/releases/download/$(version)/libass-$(version).tar.xz",
                  "1cdd39c9d007b06e737e7738004d7f38cf9b1e92843f37307b24e7ff63ab8e53"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libass-*
apk add nasm
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --disable-require-system-font-provider --disable-static
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(p -> arch(p) != "armv6l", supported_platforms(; experimental=true))

# The products that we will ensure are always built
products = [
    LibraryProduct("libass", :libass),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("FreeType2_jll"),
    Dependency("FriBidi_jll"),
    Dependency("HarfBuzz_jll"; compat="2.8.1"),
    Dependency("Bzip2_jll"; compat="1.0.8"),
    Dependency("Zlib_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
include("posterior_discrete_examinee.jl")
include("posterior_discrete_examinee_matrix.jl")
#include("posterior_discrete_examinee_quick.jl")
include("posterior_pg_examinee.jl")
include("posterior_pg_item.jl")


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
]activate .

# %%
]st

# %%
]add CairoMakie LaTeXStrings

# %%
]st

# %%
using CairoMakie, LaTeXStrings

xs = 0.0:0.01:2.0
lines(xs, sin.(π.*xs), axis=(title=L"\sin{x}",))

# %%
using Pkg
println("Julia v", VERSION)
Pkg.status("CairoMakie")
Pkg.status("Makie"; mode=PKGMODE_MANIFEST)
Pkg.status("LaTeXStrings")

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
@inline function l2_norm(z::AbstractVector)
    return sqrt(sum(abs2, z))
end

@inline function l2_distance(x::AbstractVector, y::AbstractVector)
    z = x - y
    return l2_norm(z)
end


text:
using Images, FileIO
using CSV
using ImageTransformations
using Distributed

function read_images(path, resize_image_dim_x, resize_image_dim_y)
    image_paths = map(path -> "$image_folder_path/$path", filter(path -> contains(path, ".png"), readdir(image_folder_path)))
    loaded_images = map(path -> Float32.(channelview(load(path))), (image_paths))
    channel_size =  size(loaded_images[1])[begin]
    loaded_images = pmap(img -> imresize(img, channel_size, resize_image_dim_x, resize_image_dim_y), loaded_images)
    single_image_size = size(loaded_images[1])
    vcatted_images = vcat(loaded_images...)
    reshape(vcatted_images, (length(image_paths), single_image_size...))
end

function read_labels_as_mask(path, image_x_dim, image_y_dim, number_of_categories)
    labels = []
    for row in CSV.File(path)
        if !contains(row.scorer, ".png")
            continue
        end
        points = values(row)[2:end]
        points = map(e -> parse(Float32, e), points)
        reshaped_points = reshape(points, (number_of_categories, 2))
        label = ones(image_x_dim, image_y_dim)
        for (i, point) in  enumerate(eachrow(reshaped_points))
            rounded_x = convert(Int64, round(round(point[1]), digits=0))
            rounded_y = convert(Int64, round(round(point[2]), digits=0))
            label[rounded_x, rounded_y] = i + 1
        end
        push!(labels, label)
    end
    vcatted_labels = vcat(labels...)
    reshaped_labels = reshape(vcatted_labels, (length(labels), image_x_dim, image_y_dim))
    return reshaped_labels
end

# compress data between 0 - 1
function read_labels_as_keypoints(path, image_x_dim, image_y_dim, number_of_categories)
    labels = []
    for row in CSV.File(path)
        if !contains(row.scorer, ".png")
            continue
        end
        points = values(row)[2:end]
        points = map(e -> parse(Float32, e), points)

        for (i, point) in enumerate(points)
            points[i] = point / ( isodd(i) ?  image_x_dim : image_y_dim )
        end

        reshaped_points = reshape(points, (number_of_categories * 2, 1))
        push!(labels, reshaped_points)
    end
    vcatted_labels = vcat(labels...)
    reshaped_labels = reshape(vcatted_labels, (length(labels), number_of_categories * 2))
    return reshaped_labels
end

function read_all_data(path,
        image_x_dim,
        image_y_dim,
        number_of_categories;
        resized_image_dim_x,
        resized_image_dim_y
     )
    image_data = read_images(path, resized_image_dim_x, resized_image_dim_y)
    label_data = read_labels_as_keypoints(path, image_x_dim, image_y_dim, number_of_categories)
    return (image_data, label_data)
end


text:
rdata000,cdata000=simpleprover("data/data000.cnf",2,2)
rdata000a,cdata000a=simpleprover("data/data000a.cnf",2,2)
r20Q01,c20Q01=simpleprover("data/20Q01.cnf",5,2)
rdata001,cdata001=simpleprover("data/data001.cnf",3,5)
rdata002,cdata002=simpleprover("data/data002.cnf",3,5)
rdata003,cdata003=simpleprover("data/data003.cnf",3,5)
rdata004,cdata004=simpleprover("data/data004.cnf",5,2)
rdata005,cdata005=simpleprover("data/data005.cnf",5,5)
rdata006,cdata006=simpleprover("data/data006.cnf",5,2)
rdata007,cdata007=simpleprover("data/data007.cnf",15,5)
rdata008,cdata008=simpleprover("data/data008.cnf",3,5)
rdata010,cdata010=simpleprover("data/data010.cnf",5,2)
rdata011,cdata011=simpleprover("data/data011.cnf",5,2) #0
rdata012,cdata012=simpleprover("data/data012.cnf",5,2) #0
rdata013,cdata013=simpleprover("data/data013.cnf",5,2) #0
rdata014,cdata014=simpleprover("data/data014.cnf",5,2) #0
# if config.jl has evalon=true, then 15 make 1[], else 0[]
rdata015,cdata015=simpleprover("data/data015.cnf",5,2)
rdata016,cdata016=simpleprover("data/data016.cnf",5,2)
rdata017,cdata017=simpleprover("data/data017.cnf",5,2)
rdata018,cdata018=simpleprover("data/data018.cnf",10,4)
rdata019,cdata019=simpleprover("data/data019.cnf",5,5)
rdata020,cdata020=simpleprover("data/data020.cnf",5,1)
rdata021,cdata021=simpleprover("data/data021.cnf",5,2)
rdata022,cdata022=simpleprover("data/data022.cnf",3,1) #1
rdb1,cdb1=simpleprover("data/db1.cnf",5,2)
reveryonedie,ceveryonedie=simpleprover("data/everyonedie.cnf",5,2)
rfriend1,cfriend1=simpleprover("data/friend1.cnf",5,2)
rfriend2,cfriend2=simpleprover("data/friend2.cnf",5,2)
rfriend3,cfriend3=simpleprover("data/friend3.cnf",5,1)
rfriend4,cfriend4=simpleprover("data/friend4.cnf",5,2)
rfriend41,cfriend41=simpleprover("data/friend41.cnf",5,2)
rfriend5,cfriend5=simpleprover("data/friend5.cnf",5,2)
rfriend6,cfriend6=simpleprover("data/friend6.cnf",5,2)
rfriend7,cfriend7=simpleprover("data/friend7.cnf",5,5) #9 []?
rfriend8,cfriend8=simpleprover("data/friend8.cnf",5,5) #8 [] ?
rfriend9,cfriend9=simpleprover("data/friend9.cnf",15,2) #0
rhave001,chave001=simpleprover("data/have001.cnf",7,1)
rhave002,chave002=simpleprover("data/have002.cnf",7,1)
rid1,cid1=simpleprover("data/id1.cnf",5,2)
rid2,cid2=simpleprover("data/id2.cnf",5,2) # 0
rid3,cid3=simpleprover("data/id3.cnf",5,2) # 2
rid4,cid4=simpleprover("data/id4.cnf",5,2) # 0
rimplya,cimplya=simpleprover("data/implya.cnf",5,2) #2
rimplyn,cimplyn=simpleprover("data/implyn.cnf",5,2) #1
rkuukai,ckuukai=simpleprover("data/kuukai.cnf",5,2) #1
rkuukai2,ckuukai2=simpleprover("data/kuukai2.cnf",5,2)
###
rmagia,cmagia=simpleprover("data/magia.cnf",5,2)

#rml001,cml001=simpleprover("data/ml001.cnf",5,2) # inpro
### from here, cheaplogic and naivelogic are different
rml002,cml002=simpleprover("data/ml002.cnf",5,2)  # evalon
rml003,cml003=simpleprover("data/ml003.cnf",5,2)  # evalon
rml004,cml004=simpleprover("data/ml004.cnf",5,2)  # evalon
rml005,cml005=simpleprover("data/ml005.cnf",5,2)  # evalon
rml006,cml006=simpleprover("data/ml006.cnf",6,2)  # evalon #1[] determined img3 battle. but no others.
rml007,cml007=simpleprover("data/ml007.cnf",50,20)  # evalon #5[] can't determin img3. 50 steps no effect. contralimit=20 needed
rml007a,cml007a=simpleprover("data/ml007a.cnf",50,20)  # evalon #5[] can't determin img3. 50 steps no effect. contralimit=20 needed
rml007b,cml007b=simpleprover("data/ml007b.cnf",50,20)  # evalon #5[] can't determin img3. 50 steps no effect. contralimit=20 needed
rml007c,cml007c=simpleprover("data/ml007c.cnf",50,20)  # evalon #5[] can't determin img3. 50 steps no effect. contralimit=20 needed
rml008,cml008=simpleprover("data/ml008.cnf",15,2)  # evalon

#rproblogic0,cproblogic0=simpleprover("data/problogic0.cnf",5,2)

rprop001,cprop001=simpleprover("data/prop001.cnf",5,2)

rtime,ctime=simpleprover("data/time.cnf",5,2) # unless repeat check, infinitely loop

rtime2,ctime2=simpleprover("data/time2.cnf",50,2)
rtime3,ctime3=simpleprover("data/time3.cnf",50,2)
rto1,cto1=simpleprover("data/to1.cnf",5,2)
rto2,cto2=simpleprover("data/to2.cnf",5,2)
rto3,cto3=simpleprover("data/to3.cnf",5,2)
rto4,cto4=simpleprover("data/to4.cnf",5,2)

rser001,cser001=simpleprover("data/ser001.cnf",5,2)
rpar001,cpar001=simpleprover("data/par001.cnf",10,2) # too many proofs free for y

rhis001,chis001=simpleprover("data/his001.cnf",10,5)
rhis002,chis002=simpleprover("data/his002.cnf",10,5)

rjp1,cjp1=simpleprover("data/jp001.cnf",2,2)

#rpprop001,cpprop001=simpleprover("data/pprop001.cnf",5,2)
rprob00,cprob00=simpleprover("data/prob00.cnf",5,2) # 2[] 2ways found
rprob01,cprob01=simpleprover("data/prob01.cnf",5,4) # intrinsic 3[]in 7  2ways found
rpar002,cpar002=simpleprover("data/par002.cnf",10,5)

rdop001,cdop001=simpleprover("data/dop001.cnf",10,10) # contralimit must > 5
rident001,cident001=simpleprover("data/ident001.cnf",5,5) # identity
rident002,cident002=simpleprover("data/ident002.cnf",5,5) # identity with near time

#### end of dirrerent results

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
@testset "$(rpad("Deterministic equations",80))" begin

    ################################################################################
    # Test ODE: Ordinary Differential Equation
    ################################################################################

    function f_ode(t, x, f)
        f[1] = x[1]
    end

    ode  = ODE(eltype(q₀), 1, 1, 1, f_ode, t₀, q₀)
    ode1 = ODE(f_ode, t₀, q₀)
    ode2 = ODE(f_ode, q₀)

    @test ode == ode1
    @test ode == ode2

    @test hash(ode1) == hash(ode2)

    @test ode == similar(ode, t₀, q₀)
    @test ode == similar(ode, q₀)

    ################################################################################
    # Test PODE: Partitioned Ordinary Differential Equation
    ################################################################################

    function v_pode(t, q, p, v)
        v[1] = q[1]
    end

    function f_pode(t, q, p, f)
        f[1] = 2p[1]
    end

    pode  = PODE(eltype(q₀), 1, 1, 1, v_pode, f_pode, t₀, q₀, p₀)
    pode1 = PODE(v_pode, f_pode, t₀, q₀, p₀)
    pode2 = PODE(v_pode, f_pode, q₀, p₀)

    @test pode == pode1
    @test pode == pode2

    @test hash(pode1) == hash(pode2)

    @test pode == similar(pode, t₀, q₀, p₀)
    @test pode == similar(pode, q₀, p₀)

    ################################################################################
    # Test IODE: Implicit Ordinary Differential Equation
    ################################################################################

    function iode_ϑ(t, q, v, p)
        p[1] = v[1]
    end

    function iode_f(t, q, v, f)
        f[1] = sin(q[1])
    end

    function iode_g(t, q, λ, g)
        g[1] = λ[1]
    end

    function iode_v(t, q, p, v)
        v[1] = p[1]
    end

    iode  = IODE(eltype(q₀), 1, 1, 1, iode_ϑ, iode_f, iode_g, t₀, q₀, p₀, λ₀; v=iode_v)
    iode1 = IODE(iode_ϑ, iode_f, iode_g, t₀, q₀, p₀, λ₀; v=iode_v)
    iode2 = IODE(iode_ϑ, iode_f, iode_g, t₀, q₀, p₀; v=iode_v)
    iode3 = IODE(iode_ϑ, iode_f, iode_g, q₀, p₀; v=iode_v)

    @test iode == iode1
    @test iode == iode2
    @test iode == iode3

    @test hash(iode1) == hash(iode2)

    @test iode == similar(iode, t₀, q₀, p₀, λ₀)
    @test iode == similar(iode, t₀, q₀, p₀)
    @test iode == similar(iode, q₀, p₀)

    ################################################################################
    # Test VODE: Variational Ordinary Differential Equation
    ################################################################################

    vode  = VODE(eltype(q₀), 1, 1, 1, iode_ϑ, iode_f, iode_g, t₀, q₀, p₀, λ₀; v=iode_v)
    vode1 = VODE(iode_ϑ, iode_f, iode_g, t₀, q₀, p₀, λ₀; v=iode_v)
    vode2 = VODE(iode_ϑ, iode_f, iode_g, t₀, q₀, p₀; v=iode_v)
    vode3 = VODE(iode_ϑ, iode_f, iode_g, q₀, p₀; v=iode_v)

    @test vode == vode1
    @test vode == vode2
    @test vode == vode3

    @test hash(vode1) == hash(vode2)

    @test vode == similar(vode, t₀, q₀, p₀, λ₀)
    @test vode == similar(vode, t₀, q₀, p₀)
    @test vode == similar(vode, q₀, p₀)

    ################################################################################
    # Test DAE: Differential Algebraic Equation
    ################################################################################

    function v_dae(t, x, v)
        v[1] = x[1]
        v[2] = x[2]
    end

    function u_dae(t, x, λ, u)
        u[1] = +λ[1]
        u[2] = -λ[1]
    end

    function ϕ_dae(t, x, λ, ϕ)
        ϕ[1] = x[2] - x[1]
    end

    dae  = DAE(eltype(q₀), 1, 2, 1, 1, v_dae, u_dae, ϕ_dae, t₀, x₀, λ₀)
    dae1 = DAE(v_dae, u_dae, ϕ_dae, t₀, x₀, λ₀)
    dae2 = DAE(v_dae, u_dae, ϕ_dae, x₀, λ₀)

    @test dae == dae1
    @test dae == dae2

    @test hash(dae1) == hash(dae2)

    @test dae == similar(dae, t₀, x₀, λ₀)
    @test dae == similar(dae, x₀, λ₀)

    ################################################################################
    # Test PDAE: Partitioned Differential Algebraic Equation
    ################################################################################

    function pdae_v(t, q, p, v)
        v[1] = q[1]
    end

    function pdae_f(t, q, p, f)
        f[1] = p[1]
    end

    function pdae_p(t, q, v, p)
        p[1] = v[1]
    end

    function pdae_u(t, q, p, λ, u)
        u[1] = λ[1]
    end

    function pdae_g(t, q, p, λ, g)
        g[1] = λ[1]
    end

    function pdae_ϕ(t, q, p, λ, ϕ)
        ϕ[1] = p[1] - q[1]
    end

    pdae  = PDAE(eltype(q₀), 1, 1, 1, 1, pdae_v, pdae_f, pdae_u, pdae_g, pdae_ϕ, t₀, q₀, p₀, λ₀)
    pdae1 = PDAE(pdae_v, pdae_f, pdae_u, pdae_g, pdae_ϕ, t₀, q₀, p₀, λ₀)
    pdae2 = PDAE(pdae_v, pdae_f, pdae_u, pdae_g, pdae_ϕ, q₀, p₀, λ₀)

    @test pdae == pdae1
    @test pdae == pdae2

    @test hash(pdae1) == hash(pdae2)

    @test pdae == similar(pdae, t₀, q₀, p₀, λ₀)
    @test pdae == similar(pdae, t₀, q₀, p₀)
    @test pdae == similar(pdae, q₀, p₀)

    idae  = IDAE(eltype(q₀), 1, 1, 1, 1, pdae_p, pdae_f, pdae_u, pdae_g, pdae_ϕ, t₀, q₀, p₀, λ₀)
    idae1 = IDAE(pdae_p, pdae_f, pdae_u, pdae_g, pdae_ϕ, t₀, q₀, p₀, λ₀)
    idae2 = IDAE(pdae_p, pdae_f, pdae_u, pdae_g, pdae_ϕ, q₀, p₀, λ₀)

    @test idae == idae1
    @test idae == idae2

    @test hash(idae1) == hash(idae2)

    @test idae == similar(idae, t₀, q₀, p₀, λ₀)
    @test idae == similar(idae, t₀, q₀, p₀)
    @test idae == similar(idae, q₀, p₀)

    ################################################################################
    # Test VDAE: Variational Differential Algebraic Equation
    ################################################################################

    function vdae_ψ(t, q, p, λ, μ, ψ)
        ψ[1] = μ[1] - λ[1]
    end

    vdae  = VDAE(eltype(q₀), 1, 1, 1, 1, iode_ϑ, iode_f, iode_g, iode_g, pdae_ϕ, vdae_ψ, t₀, q₀, p₀, λ₀, λ₀; v=iode_v)
    vdae1 = VDAE(iode_ϑ, iode_f, iode_g, iode_g, pdae_ϕ, vdae_ψ, t₀, q₀, p₀, λ₀, λ₀; v=iode_v)
    vdae2 = VDAE(iode_ϑ, iode_f, iode_g, iode_g, pdae_ϕ, vdae_ψ, t₀, q₀, p₀, λ₀; v=iode_v)
    vdae3 = VDAE(iode_ϑ, iode_f, iode_g, iode_g, pdae_ϕ, vdae_ψ, t₀, q₀, p₀; v=iode_v)
    vdae4 = VDAE(iode_ϑ, iode_f, iode_g, iode_g, pdae_ϕ, vdae_ψ, q₀, p₀; v=iode_v)

    @test vdae == vdae1
    @test vdae == vdae2
    @test vdae == vdae3
    @test vdae == vdae4

    @test hash(vdae1) == hash(vdae2)
    @test hash(vdae3) == hash(vdae4)

    @test vdae == similar(vdae, t₀, q₀, p₀, λ₀, λ₀)
    @test vdae == similar(vdae, t₀, q₀, p₀, λ₀)
    @test vdae == similar(vdae, t₀, q₀, p₀)
    @test vdae == similar(vdae, q₀, p₀)

    ################################################################################
    # Test SODE: Split Ordinary Differential Equation
    ################################################################################

    function f_sode_1(t, x, f)
        f[1] = x[1]
    end

    function f_sode_2(t, x, f)
        f[1] = x[1]^2
    end

    f_sode = (f_sode_1, f_sode_2)

    sode  = SODE(eltype(q₀), 1, 1, 1, f_sode, t₀, q₀)
    sode1 = SODE(f_sode, t₀, q₀)
    sode2 = SODE(f_sode, q₀)

    @test sode == sode1
    @test sode == sode2

    @test hash(sode1) == hash(sode2)

    @test sode == similar(sode, t₀, q₀)
    @test sode == similar(sode, q₀)

end