text:
function build_jlModel()
    fields = Tuple{Symbol,Type}[(:m, mjModel), (:cptr, Ptr{mjModel})]
    ctorargs = []

    getpropbody = :(getfield(value, name))
    setpropbody = :(setfield!(value, name, x))
    propnames = Symbol[]
    skippedfields = [:buffer, :stack]

    #forwardedfields = []
    sinfo = STRUCT_INFO[:mjModel]
    for (name, T) in zip(fieldnames(mjModel), fieldtypes(mjModel))
        if name in skippedfields
            continue
        end
        push!(propnames, name)
        if T <: Ptr
            if haskey(sinfo, name)
                finfo = sinfo[name]
                N = length(finfo.mjsize)
                T = Array{eltype(T),N}
                push!(fields, (name, T))
                push!(ctorargs, quote
                    ptr = getfield(m, $(QuoteNode(name)))
                    size = lookupsize(:mjModel, m, $(QuoteNode(name)))
                    Base.unsafe_wrap(Array, ptr, size)
                end)
                if !isnothing(finfo.dsize)
                    getpropbody = Expr(
                        :elseif,
                        :(name === $(QuoteNode(name))),
                        quote
                            A = getfield(value, $(QuoteNode(name)))
                            m = getfield(value, :m)
                            sz = lookup_dsize(
                                Val{mjModel}(),
                                m,
                                Val{$(QuoteNode(name))}(),
                            )
                            UnsafeArray(pointer(A), sz)
                        end,
                        getpropbody,
                    )
                end
            else
                @warn "Found Ptr $name but no size, skipping"
            end
        elseif T <: StaticArray
            N = ndims(T)
            elT = eltype(T)
            sz = size(T)
            Tnew = Array{eltype(T),N}
            foffset = fieldoffset(mjModel, Base.fieldindex(mjModel, name))

            push!(fields, (name, Tnew))
            push!(ctorargs, quote
                ptr = Ptr{$elT}(pm + $foffset)
                Base.unsafe_wrap(Array, ptr, $sz)
            end)
        else
            # forward
            getpropbody = Expr(
                :elseif,
                :(name === $(QuoteNode(name))),
                :(getfield(getfield(value, :m), $(QuoteNode(name)))),
                getpropbody,
            )
            setpropbody = Expr(
                :elseif,
                :(name === $(QuoteNode(name))),
                :(setfield!(getfield(value, :m), $(QuoteNode(name)), convert($T, x))),
                setpropbody,
            )
        end
    end
    structdef = build_structdef(
        :jlModel,
        fields,
        mutable = true,
        supertype = AbstractMJModel,
    )
    structexpr = combinestructdef(structdef)

    getpropbody = Expr(:if, getpropbody.args...)
    setpropbody = Expr(:if, setpropbody.args...)

    quote
        $structexpr
        function Base.getproperty(value::jlModel, name::Symbol)
            $getpropbody
        end
        function Base.setproperty!(value::jlModel, name::Symbol, x)
            $setpropbody
        end
        function Base.propertynames(::Union{jlModel,Type{jlModel}})
            $propnames
        end
        function jlModel(pm::Ptr{mjModel}; own::Bool = false)
            m = unsafe_load(pm)
            jm = jlModel(m, pm, $(ctorargs...))
            # must free C-allocated Ptr, not pointer to m!
            own && Base.finalizer(jm -> MJCore.mj_deleteModel(jm.cptr), jm)
            jm
        end
        function jlModel(path::String)
            kind = MJCore.check_modelpath(path)
            pm = kind === :MJCF ? mj_loadXML(path) : mj_loadModel(path)
            pm == C_NULL && @mjerror "$path not a valid MJCF or MJB file"
            jlModel(pm, own = true)
        end
        Base.cconvert(::Type{Ptr{mjModel}}, jm::jlModel) = Base.cconvert(Ptr{mjModel}, jm.m)
        jlModel(jm::jlModel) = Base.copy(jm)
        Base.deepcopy_internal(jm::jlModel, ::IdDict) = Base.copy(jm)
        Base.copy(jm::jlModel) = jlModel(mj_copyModel(jm.cptr), own = true)
    end
end

@eval $(build_jlModel())


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Author: Maur좩cio Carvalho Maciel, Paula B. Staudt, Rafael P. Soares
#* $Id$
#*--------------------------------------------------------------------
type splitter3
	splitter3()=begin
		new(
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet stream",
				:PosX=>0,
				:PosY=>0.5001,
				:Symbol=>"_{in}"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Outlet stream 1",
				:PosX=>1,
				:PosY=>0.25,
				:Symbol=>"_{Out1}"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Outlet stream 2",
				:PosX=>1,
				:PosY=>0.5059,
				:Symbol=>"_{Out2}"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Outlet stream 3",
				:PosX=>1,
				:PosY=>0.75,
				:Symbol=>"_{Out3}"
			)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Distribution of Outlets",
				:Default=>0.33,
				:Symbol=>"\\phi"
			)),(3)),
			[
				:(sum(FlowRatios) = 1),
				:(Outlet1.F = Inlet.F*FlowRatios(1)),
				:(Outlet2.F = Inlet.F*FlowRatios(2)),
				:(Outlet3.F = Inlet.F*FlowRatios(3)),
				:(Outlet1.z = Inlet.z),
				:(Outlet2.z = Inlet.z),
				:(Outlet3.z = Inlet.z),
				:(Outlet1.P = Inlet.P),
				:(Outlet2.P = Inlet.P),
				:(Outlet3.P = Inlet.P),
				:(Outlet1.h = Inlet.h),
				:(Outlet2.h = Inlet.h),
				:(Outlet3.h = Inlet.h),
				:(Outlet1.T = Inlet.T),
				:(Outlet2.T = Inlet.T),
				:(Outlet3.T = Inlet.T),
				:(Outlet1.v = Inlet.v),
				:(Outlet2.v = Inlet.v),
				:(Outlet3.v = Inlet.v),
			],
			[
				"Normalize Flow Ratios","Outlet1 Flow","Outlet2 Flow","Outlet3 Flow","Outlet1 Composition","Outlet2 Composition","Outlet3 Composition","Outlet1 Pressure","Outlet2 Pressure","Outlet3 Pressure","Outlet1 Enthalpy","Outlet2 Enthalpy","Outlet3 Enthalpy","Outlet1 Temperature","Outlet2 Temperature","Outlet3 Temperature","Outlet1 Vapourisation Fraction","Outlet2 Vapourisation Fraction","Outlet3 Vapourisation Fraction",
			],
			[:Inlet,:Outlet1,:Outlet2,:Outlet3,:FlowRatios,]
		)
	end
	Inlet::stream
	Outlet1::stream
	Outlet2::stream
	Outlet3::stream
	FlowRatios::Array{fraction}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export splitter3
function setEquationFlow(in::splitter3)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	addEquation(11)
	addEquation(12)
	addEquation(13)
	addEquation(14)
	addEquation(15)
	addEquation(16)
	addEquation(17)
	addEquation(18)
	addEquation(19)
end
function atributes(in::splitter3,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/splitter"
	fields[:Brief]="Model of a splitter with 3 outlet streams"
	fields[:Info]="== Assumptions ==
* Thermodynamics equilibrium
* Adiabatic

== Specify ==
*The inlet stream
*Two FlowRatios of split of the outlet streams:

	FlowRatios(i) = (Mole Flow of the outlet stream i /
				Mole Flow of the inlet stream)
				where i = 1, 2, 3
"
	drive!(fields,_)
	return fields
end
splitter3(_::Dict{Symbol,Any})=begin
	newModel=splitter3()
	newModel.attributes=atributes(newModel,_)
	newModel
end

ext: jl


text:
function solve!(dvars, kvars, auxvars, Δτ, n₀, n, rtol)
    @unpack F, Fˢ, ζ = dvars
    @unpack svars, Λ = kvars
    @unpack S, Sˢ, B, Bˢ, w, υ = svars
    @unpack A, Aˢ, Z, Fᵢ, Fˢᵢ, ΔF₁, ΔFˢ₁, Δζ, ζoζ = auxvars

    fill_aux_svars!(A, Aˢ, ΔF₁, ΔFˢ₁, S, Sˢ, Bˢ, Λ, F, Fˢ, ζ[1], Δτ)
    fill_aux_tvars!(Δζ, ζoζ, ζ, n₀ - 1)

    for i = n₀:n
        perform_linstep!(Fᵢ, Fˢᵢ, F, Fˢ, ΔF₁, ΔFˢ₁, B, Bˢ, ζ, Δζ, ζoζ, Λ, Δτ, i)
        perform_nlstep!(F, Fˢ, Fᵢ, Fˢᵢ, ΔF₁, ΔFˢ₁, A, Aˢ, B, Bˢ, Z, ζ, Λ, w, υ, Δτ, i, rtol)
        diff!(Δζ, ζ, i)
        conv!(ζoζ, ζ, i)
    end

    return dvars
end

function fill_aux_svars!(A, Aˢ, ΔF₁, ΔFˢ₁, S, Sˢ, Bˢ, Λ, F, Fˢ, ζ₁, Δτ)
    Δτ⁻¹ = inv(Δτ)

    @inbounds for j in eachindex(S)
        ΔF₁[j]  = S[j] - F[1, j]
        ΔFˢ₁[j] = Sˢ[j] - Fˢ[1, j]

        m = Λ[j] * ζ₁
        a = alpha_a(Δτ, Δτ⁻¹, m)
        b = alpha_b(Bˢ[j], Δτ, m)
        A[j]  = inv(a + b * inv(S[j]))
        Aˢ[j] = inv(a + b)
    end

    return A, Aˢ, ΔF₁, ΔFˢ₁
end

alpha_a(Δτ, Δτ⁻¹, m) = (Δτ⁻¹ + m)
function alpha_a(Δτ, Δτ⁻¹, m::TR)
    t = Δτ⁻¹ + m.t
    r = t + m.r * (1 + Δτ * m.t)
    return TR(t, r)
end

alpha_b(k²D₀, Δτ, m) = k²D₀
function alpha_b(k²D₀::LDProjections{1}, Δτ, m::TR)
    t = k²D₀.t
    r = t * (1 + Δτ * m.r) + (1 + Δτ * m.t) * k²D₀.r[1]
    return LDProjections(t, SVector(r))
end
function alpha_b(k²D₀, Δτ, m::TR)
    t = k²D₀.t
    r = t * (1 + Δτ * m.r) .+ (1 + Δτ * m.t) .* k²D₀.r
    return LDProjections(t, r)
end

function fill_aux_tvars!(Δζ, ζoζ::Nothing, ζ, n)
    @inbounds for i = 2:n
        diff!(Δζ, ζ, i)
    end
    return Δζ
end

function fill_aux_tvars!(Δζ, ζoζ, ζ, n)
    conv!(ζoζ, ζ, 1)
    @inbounds for i = 2:n
        diff!(Δζ, ζ, i)
        conv!(ζoζ, ζ, i)
    end
    return Δζ
end

### Linear contributions of the SCGLE equations
function perform_linstep!(Fᵢ, Fˢᵢ, F, Fˢ, ΔF₁, ΔFˢ₁, B, Bˢ, ζ, Δζ, ζoζ, Λ, Δτ, n)
    n₋ = n - 1
    ζn₋ = ζ[n₋]

    @inbounds for j in eachindex(Λ)
        λ, β, βˢ, ΔF₁ⱼ, ΔFˢ₁ⱼ, Δτλλ = index_lin_svars(Λ, B, Bˢ, ΔF₁, ΔFˢ₁, Δτ, j)

        Fₐ  = initialize_linterm(F[1, j], ζn₋, ζoζ, n₋)
        Fˢₐ = initialize_linterm(Fˢ[1, j], ζn₋, ζoζ, n₋)

        for i = 2:n₋
            l = n - i + 1
            Δζₗ, Δτζₗ, ζᵢζₗ, Δζoζₗ = index_lin_tvars(Δζ, ζ, ζoζ, Δτ, l, i)

            Fₐ  = accumulate_linterm(Fₐ, F[i, j], Δζₗ, Δτζₗ, ζᵢζₗ, Δζoζₗ, β, ΔF₁ⱼ)
            Fˢₐ = accumulate_linterm(Fˢₐ, Fˢ[i, j], Δζₗ, Δτζₗ, ζᵢζₗ, Δζoζₗ, βˢ, ΔFˢ₁ⱼ)
        end

        Fᵢ[j]  = reduce_linterm(Fₐ, λ, F[n₋, j], Δτ, Δτλλ)
        Fˢᵢ[j] = reduce_linterm(Fˢₐ, λ, Fˢ[n₋, j], Δτ, Δτλλ)
    end

    return Fᵢ, Fˢᵢ
end

index_lin_svars(Λ, B, Bˢ, ΔF₁, ΔFˢ₁, Δτ, j) =
    (Λ[j], nothing, nothing, nothing, nothing, nothing)
#
index_lin_svars(Λ::Vector{<:TR}, B, Bˢ, ΔF₁, ΔFˢ₁, Δτ, j) =
    (Λ[j], B[j], Bˢ[j], ΔF₁[j], ΔFˢ₁[j], Δτ * Λ[j].t * Λ[j].r)

index_lin_tvars(Δζ, ζ, ζoζ::Nothing, Δτ, l, i) =
    (Δζ[l], nothing, nothing, nothing)
#
index_lin_tvars(Δζ, ζ, ζoζ, Δτ, l, i) =
    (Δζ[l], -Δτ * ζ[l], ζ[i].t * ζ[l].r, ζoζ[l - 1] - ζoζ[l])

initialize_linterm(F₁ⱼ, ζₙ, ζoζ::Nothing, n) = F₁ⱼ * ζₙ
#
function initialize_linterm(Fᵢⱼ, ζₙ, ζoζ, n)
    Fₐᵀ   = Fᵢⱼ.t * ζₙ.t  # Center-of-mass projection with coefficient λᵀ
    Fₐᴿᵀ  = Fᵢⱼ.r * ζₙ.t  # Projections with coefficient λᵀ
    Fₐᴿᴿ  = Fᵢⱼ.r * ζₙ.r  # Projections with coefficient λᴿ
    Fₐᴿᵀᴿ = Fᵢⱼ.r * ζoζ[n] # Projections with coefficient λᵀλᴿ

    return Fₐᵀ, Fₐᴿᵀ, Fₐᴿᴿ, Fₐᴿᵀᴿ
end

accumulate_linterm(Fₐ, Fᵢⱼ, Δζₗ, Δτζₗ::Nothing, ζᵢζₗ, Δζoζₗ, β, ΔF₁ⱼ) = muladd(Fᵢⱼ, Δζₗ, Fₐ)
#
function accumulate_linterm(Fₐ, Fᵢⱼ, Δζₗ, Δτζₗ, ζᵢζₗ, Δζoζₗ, β, ΔF₁ⱼ)
    Fₐᵀ, Fₐᴿᵀ, Fₐᴿᴿ, Fₐᴿᵀᴿ = Fₐ

    Fₐᵀ = muladd(Δζₗ.t , Fᵢⱼ.t, Fₐᵀ)
    Fₐᴿᵀ  = muladd.(muladd.(Δτζₗ.t, getr(β.r), Δζₗ.t), Fᵢⱼ.r, Fₐᴿᵀ)
    Fₐᴿᴿ  = muladd.(muladd.(Δτζₗ.r, getr(β.t), Δζₗ.r), Fᵢⱼ.r, Fₐᴿᴿ)
    Fₐᴿᵀᴿ = muladd.(Δζoζₗ, Fᵢⱼ.r, muladd.(ζᵢζₗ, ΔF₁ⱼ.r, Fₐᴿᵀᴿ))

    return Fₐᵀ, Fₐᴿᵀ, Fₐᴿᴿ, Fₐᴿᵀᴿ
end

reduce_linterm(Fₐ, λ, Fᵢⱼ, Δτ, Δτλλ::Nothing) = muladd(λ, Fₐ, Fᵢⱼ / Δτ)
#
function reduce_linterm(Fₐ, λ, Fᵢⱼ, Δτ, Δτλλ)
    Fₐᵀ, Fₐᴿᵀ, Fₐᴿᴿ, Fₐᴿᵀᴿ = Fₐ

    F̃ = Fᵢⱼ / Δτ
    t = muladd(λ.t, Fₐᵀ, F̃.t)
    r = muladd.(λ.t, Fₐᴿᵀ, muladd.(λ.r, Fₐᴿᴿ, muladd.(Δτλλ, Fₐᴿᵀᴿ, F̃.r)))

    return Projections.constructorname(F̃)(t, r)
end

### Nonlinear contributions of the SCGLE equations
function perform_nlstep!(F, Fˢ, Fᵢ, Fˢᵢ, ΔF₁, ΔFˢ₁, A, Aˢ, B, Bˢ, Z, ζ, Λ, w, υ, Δτ, i, rtol)
    ζ₁ = ζ[1]
    ζ[i] = init_guess(ζ, i)

    F₁ = view_firststep(F, 1)
    Fˢ₁ = view_firststep(Fˢ, 1)

    @inbounds while true
        ζᵢ = ζ[i]

        for j in eachindex(Λ)
            m, βm, βˢm = compute_nlmemory(Λ, ζᵢ, ζ₁, B, Bˢ, Δτ, j)
            F[i, j]  = reduce_nlterm(A[j], m, ΔF₁[j], Fᵢ[j], F₁, βm, j)
            Fˢ[i, j] = reduce_nlterm(Aˢ[j], m, ΔFˢ₁[j], Fˢᵢ[j], Fˢ₁, βˢm, j)
            Z[j] = product(w[j], F[i, j], Fˢ[i, j])
        end

        ζ[i] = υ * sum(Z)

        if isapprox(ζᵢ, ζ[i]; rtol = rtol, nans = true)
            break
        end
    end
end

init_guess(ζ, i) = ζ[i - 1]

view_firststep(F, i) = nothing
view_firststep(F::Array{<:Projections.AbstractProjections}, i) = view(F, i, :)

compute_nlmemory(Λ, ζᵢ, ζ₁, B, Bˢ, Δτ, j) = (Λ[j] * ζᵢ, nothing, nothing)
#
function compute_nlmemory(Λ, ζᵢ::TR, ζ₁, B, Bˢ, Δτ, j)
    λ = Λ[j]

    m = λ * ζᵢ
    moᵢm = m.t * (1 + Δτ * λ.r * ζ₁.r) + m.r * (1 + Δτ * λ.t * ζ₁.t)
    βm  = -Δτ .* (getr(B[j].t) * m.r .+ getr(B[j].r) .* m.t)
    βˢm = -Δτ .* (Bˢ[j].t * m.r .+ Bˢ[j].r .* m.t)

    return (m.t, moᵢm), βm, βˢm
end

reduce_nlterm(α, m, ΔF₁ⱼ, Fᵢⱼ, F₁ⱼ, βm::Nothing, j) = α * muladd(m, ΔF₁ⱼ, Fᵢⱼ)
#
function reduce_nlterm(α, m, ΔF₁ⱼ, Fᵢⱼ, F₁, βm, j)
    mᵀ, moᵢm = m
    F₁ⱼ = F₁[j]

    t = α.t * muladd(mᵀ, ΔF₁ⱼ.t, Fᵢⱼ.t)
    r = α.r .* muladd.(moᵢm, ΔF₁ⱼ.r, muladd.(βm, F₁ⱼ.r, Fᵢⱼ.r))

    return Projections.constructorname(Fᵢⱼ)(t, r)
end


text:
module Configuration

using DocStringExtensions

export Config

"""
$(TYPEDEF)
$(TYPEDFIELDS)

- `threshold` : A value between `0` and `1`, which controls the probability threshold at
    which expected improvement criteria is modeled.
- `linear_forgetting` : A positive value which controls the number of historic points which
    are used for probabilistic modelling, and older points beyond this are linearly
    de-weighted.
- `draws` : A positive value which controls the number of samples to draw when making a
    recommendation for next optimisation canidate.
- `random_trials` : A positive value which controls the number of trials of randomly
    generated candidate points before TreeParzen optimisation is used.
- `prior_weight` : A value between `0` and `1`, which controls importance of user specified
    probabilistic parameters vs the history of trials.
"""
struct Config
    threshold::Float64
    linear_forgetting::Int
    draws::Int
    random_trials::Int
    prior_weight::Float64

    function Config(
        threshold::Float64, linear_forgetting::Int, draws::Int, random_trials::Int,
        prior_weight::Float64
    )
        if !(0 <= threshold <= 1)
            throw(ArgumentError("threshold must be between 0 and 1"))
        end
        if linear_forgetting < 0
            throw(ArgumentError("linear_forgetting must be positive"))
        end
        if draws < 0
            throw(ArgumentError("draws must be positive"))
        end
        if random_trials < 0
            throw(ArgumentError("random_trials must be positive"))
        end
        if !(0 <= prior_weight <= 1)
            throw(ArgumentError("prior_weight must be between 0 and 1"))
        end

        return new(threshold, linear_forgetting, draws, random_trials, prior_weight)
    end
end

"""
$(TYPEDSIGNATURES)

Keyword argument outer constructor
- `threshold::Float64`; default 0.25.
- `linear_forgetting::Int`; default 25.`
- `draws::Int`; default 24. `
- `random_trials::Int`; default 20.
- `prior_weight::Float64`; default 1.0.
"""
function Config(;
    threshold::Float64 = 0.25,
    linear_forgetting::Int = 25,
    draws::Int = 24,
    random_trials::Int = 20,
    prior_weight::Float64 = 1.0,
    )
    return Config(threshold, linear_forgetting, draws, random_trials, prior_weight)
end

end # module Configuration


text:
immutable ScatteringLayer <: Mocha.Layer
    bank::AbstractBank
    bottoms::Vector{Symbol}
    name::AbstractString
    previous_bank::AbstractBank
    sibling_mask_factor::Float64
    tops::Vector{Symbol}
end

Mocha.can_do_bp(::ScatteringLayer) = false
Mocha.has_neuron(::ScatteringLayer) = false
Mocha.has_param(::ScatteringLayer) = false

function ScatteringLayer( ;
        bank::AbstractBank = NullBank(),
        bottoms::Vector{Symbol} = Symbol[],
        name::AbstractString = "scattering",
        previous_bank::AbstractBank = NullBank(),
        sibling_mask_factor = 1.0,
        tops::Vector{Symbol} = Symbol[])
    @assert length(bottoms) == 1
    @assert length(tops) == 1
    ScatteringLayer(bank, bottoms, name, previous_bank, sibling_mask_factor,
        tops)
end


text:
using PLEXOSUtils

# TODO: Actually test things
function runtest(zippath::String)
    PLEXOSSolutionDatasetSummary(zippath)
    PLEXOSSolutionDataset(zippath)
end

testfolder = dirname(@__FILE__) * "/"

zipfiles = ["Model Base_8200 Solution.zip",
            "Model DAY_AHEAD_NO_TX Solution.zip",
            "Model DAY_AHEAD_ALL_TX Solution.zip",
            "Model DA_h2hybrid_SCUC_select_lines_Test_1day Solution.zip",
            "Model DAY_AHEAD_PRAS Solution.zip"]

for zipfile in zipfiles
    println(zipfile)
    runtest(testfolder * zipfile)
end


text:
abstract type AbstractConvectiveAdjustmentParameters <: AbstractParameters end

const ModelWithConvectiveAdjustment = Model{L, K, <:AbstractConvectiveAdjustmentParameters} where {L, K}

@inline Cᴬu(m) = Cᴬu(m.convective_adjustment)
@inline Cᴬc(m) = Cᴬc(m.convective_adjustment)
@inline Cᴬe(m) = Cᴬe(m.convective_adjustment)

Base.@kwdef struct UnityPrandtlConvectiveAdjustment{T} <: AbstractConvectiveAdjustmentParameters
    Cᴬ :: T = 1.0
end

@inline Cᴬu(p::UnityPrandtlConvectiveAdjustment) = p.Cᴬ
@inline Cᴬc(p::UnityPrandtlConvectiveAdjustment) = p.Cᴬ
@inline Cᴬe(p::UnityPrandtlConvectiveAdjustment) = p.Cᴬ

Base.@kwdef struct VariablePrandtlConvectiveAdjustment{T} <: AbstractConvectiveAdjustmentParameters
    Cᴬu :: T = 0.5
    Cᴬc :: T = 1.0
    Cᴬe :: T = 0.5
end

@inline Cᴬu(p::VariablePrandtlConvectiveAdjustment) = p.Cᴬu
@inline Cᴬc(p::VariablePrandtlConvectiveAdjustment) = p.Cᴬc
@inline Cᴬe(p::VariablePrandtlConvectiveAdjustment) = p.Cᴬe

@inline KU(m::ModelWithConvectiveAdjustment, i) = KU₀(m) + K_adjustment(Cᴷu(m, i), Cᴬu(m), m, i)
@inline KV(m::ModelWithConvectiveAdjustment, i) = KV₀(m) + K_adjustment(Cᴷv(m, i), Cᴬu(m), m, i)
@inline KT(m::ModelWithConvectiveAdjustment, i) = KT₀(m) + K_adjustment(CᴷT(m, i), Cᴬc(m), m, i)
@inline KS(m::ModelWithConvectiveAdjustment, i) = KS₀(m) + K_adjustment(CᴷS(m, i), Cᴬc(m), m, i)
@inline Ke(m::ModelWithConvectiveAdjustment, i) = Ke₀(m) + K_adjustment(Cᴷe(m, i), Cᴬe(m), m, i)

@inline K_adjustment(Cᴷϕ, Cᴬϕ, m, i) =
    @inbounds ifelse(∂B∂z(m, i) < 0 && m.state.Qb > 0,
                     Cᴬϕ * onface(m.solution.e, i)^2 / m.state.Qb,
                     Cᴷϕ * onface(m.state.K, i))


text:
function node2elem(nodes,re,ce)

# # node-2-elemsent connectivity
# [re,ce] = mkv2t(t,ne);

# elemsent-to-node connectivity in nonoverlapping subdomain Omega_i
# nodes = t(:,intelem);

# remove duplications
nodes = sort(unique(nodes[:]));

# remove zero or negative
nodes = nodes[nodes.>0];

# find all elemsents connected to nodes
nt = sum(re[nodes.+1]-re[nodes]);
elems = zeros(Int,nt,1);
j = 1;
for i = 1:length(nodes)
    ni = nodes[i];
    ri = (re[ni]+1):re[ni+1];
    elems[j:(j+length(ri)-1)] = ce[ri];
    j = j+length(ri);
end
elems = sort(unique(elems));

return elems

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
using CompScienceMeshes
using SauterSchwabQuadrature

Accuracy = 12
pd = PositiveDistance(Accuracy)
pI = point(0,0,-1)
pII = point(0,0,1)

function integrand(x,y)
			return(((x-pI)'*(y-pII))*exp(-im*1*norm(x-y))/(4pi*norm(x-y)))
end

function Sc(û)
	u = [(pi/2)*û[1]/(1-û[2]) , û[2]]	#mapping from referencetriangle to rectangle
	return(u)
end

function Tc(v̂)
	v = [v̂[1]/(1-v̂[2]), v̂[2]]			#mapping from referencetriangle to square
	return(v)
end

function INTEGRAND(û,v̂)

	ϴ = Sc(û)[1] + pi/2
	ϕ = Sc(û)[2]
	ϴ1= Tc(v̂)[1]
	ϕ1= Tc(v̂)[2]

	x = [sin(ϴ)*cos(ϕ),   sin(ϴ)*sin(ϕ),   cos(ϴ) ]		#spherical coordinates
	y = [sin(ϴ1)*cos(ϕ1), sin(ϴ1)*sin(ϕ1), cos(ϴ1)]		#spherical coordiantes

	output = integrand(x,y) * sin(ϴ)*sin(ϴ1) * (pi/2)*(1/(1-û[2]))*(1/(1-v̂[2]))
	#sin(ϴ)*sin(ϴ1) = surface element of spherical coordinates
	#(pi/2)*(1/(1-û[2]))*(1/(1-v̂[2])) = surface element of first two mappings

return(output)

end

result = sauterschwab_parameterized(INTEGRAND, pd)

println(result)


text:
lm = LakeModel()
e_0 = 0.92     # Initial employment rate
u_0 = 1 - e_0  # Initial unemployment rate
T = 50         # Simulation length

xbar = rate_steady_state(lm)
x_0 = [e_0; u_0]
x_path = simulate_rate_path(lm, x_0, T)

titles = ["Employment rate" "Unemployment rate"]
dates = collect(1:T)

plot(dates, x_path', layout=(2, 1), title=titles, legend=:none)
hline!(xbar', layout=(2, 1), color=:red, linestyle=:dash,
       ylims=[0.9999*minimum(x_path[1,:]) Inf;
       -Inf 1.001*maximum(x_path[2,:])])


text:
module tst_cov
    using Test
    using CloudCovErr
    using FITSIO
    using PyCall

    test_dir = dirname(@__FILE__)

    @testset "DECam" begin
        ref = test_dir*"/data/decapsi/c4d_170119_085651_ood_r_v1.I.fits.fz"
        @test ref == CloudCovErr.inject_rename(test_dir*"/data/decaps/c4d_170119_085651_ood_r_v1.fits.fz")

        ref_im, d_im = CloudCovErr.read_decam(test_dir*"/data/decaps/c4d_","170119_085651","r","v1","S6";corrects7=true)
        @test d_im[1,1] == 1
        @test abs(ref_im[1]-416.4855651855) < 1e-7

        decam_dir = dirname(pwd())*"/src/decam_dir"

        py"""
        import os
        os.environ['DECAM_DIR'] = $decam_dir
        """
        ref_im, d_im = CloudCovErr.read_decam(test_dir*"/data/decaps/c4d_","170119_085651","r","v1","S7";corrects7=true)
        @test d_im[1,1] == 1
        @test abs(ref_im[1]-375.392761230468) < 1e-7

        out = CloudCovErr.read_crowdsource(test_dir*"/data/decaps/","170119_085651","r","v1","S6")
        @test length(out) == 9
        wcol = out[end-1]
        w = out[end]

        out = CloudCovErr.load_psfmodel_cs(test_dir*"/data/decaps/","170119_085651","r","v1","S6")
        @test size(out(1,1,11)) == (11,11)

        if isfile(test_dir*"/data/decaps/cer/c4d_170119_085651_ooi_r_v1.cat.cer.fits")
            rm(test_dir*"/data/decaps/cer/c4d_170119_085651_ooi_r_v1.cat.cer.fits")
        end
        CloudCovErr.save_fxn(wcol,w,test_dir*"/data/decaps/","170119_085651","r","v1","S6")
        f = FITS(test_dir*"/data/decaps/cer/c4d_170119_085651_ooi_r_v1.cat.cer.fits")
        @test length(read(f["S6_CAT"],"x")) == 14572
        close(f)
        if isfile(test_dir*"/data/decaps/cer/c4d_170119_085651_ooi_r_v1.cat.cer.fits")
            rm(test_dir*"/data/decaps/cer/c4d_170119_085651_ooi_r_v1.cat.cer.fits")
        end

        f = FITS(test_dir*"/data/decaps/cat/c4d_170119_085651_ooi_r_v1.cat.fits")
        @test CloudCovErr.get_catnames(f) == ["S6","S7"]

        # run one real example end to end
        @test_nowarn CloudCovErr.proc_all(test_dir*"/data/decaps/c4d_","170119_085651","r","v1",test_dir*"/data/decaps/",ccdlist=["S6"],resume=true,corrects7=true,thr=20,outthr=20000,Np=33,widx=129,tilex=8)
        @test_nowarn CloudCovErr.proc_all(test_dir*"/data/decaps/c4d_","170119_085651","r","v1",test_dir*"/data/decaps/",ccdlist=["S6"],resume=true,corrects7=true,thr=20,outthr=20000,Np=33,widx=129,tilex=8)
        #based = test_dir*"/data/decaps/c4d_"
        #catd = test_dir*"/data/decaps/"
        #decap2f = test_dir*"/decaps2.jl"
        #@test_nowarn run(`julia $decap2f --ccdlist S6 --cS7 -r --tilex8 $based 170119_085651 r v1 $catd`)
    end
end


text:
struct EMConstantCache <: StochasticDiffEqConstantCache end
@cache struct EMCache{uType,rateType,rateNoiseType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  rtmp1::rateType
  rtmp2::rateNoiseType
end

alg_cache(alg::EM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} = EMConstantCache()

function alg_cache(alg::EM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  tmp = zero(u); rtmp1 = zero(rate_prototype);
  if noise_rate_prototype !== nothing
    rtmp2 = zero(noise_rate_prototype)
  else
    rtmp2 = nothing
  end
  EMCache(u,uprev,tmp,rtmp1,rtmp2)
end

struct SplitEMConstantCache <: StochasticDiffEqConstantCache end
@cache struct SplitEMCache{uType,rateType,rateNoiseType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  rtmp1::rateType
  rtmp2::rateNoiseType
end

alg_cache(alg::SplitEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} = SplitEMConstantCache()

function alg_cache(alg::SplitEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  tmp = zero(u); rtmp1 = zero(rate_prototype);
  rtmp2 = zero(noise_rate_prototype)
  SplitEMCache(u,uprev,tmp,rtmp1,rtmp2)
end

struct EulerHeunConstantCache <: StochasticDiffEqConstantCache end
@cache struct EulerHeunCache{uType,rateType,rateNoiseType,rateNoiseCollectionType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  ftmp1::rateType
  ftmp2::rateType
  nrtmp::rateNoiseCollectionType
  gtmp1::rateNoiseType
  gtmp2::rateNoiseType
end

alg_cache(alg::EulerHeun,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} = EulerHeunConstantCache()

function alg_cache(alg::EulerHeun,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  tmp = zero(u); ftmp1 = zero(rate_prototype); ftmp2 = zero(rate_prototype)
  nrtmp = zero(rate_prototype)
  gtmp1 = zero(noise_rate_prototype); gtmp2 = zero(noise_rate_prototype)
  EulerHeunCache(u,uprev,tmp,ftmp1,ftmp2,nrtmp,gtmp1,gtmp2)
end

struct RandomEMConstantCache <: StochasticDiffEqConstantCache end
@cache struct RandomEMCache{uType,rateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  rtmp::rateType
end

alg_cache(alg::RandomEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} = RandomEMConstantCache()

function alg_cache(alg::RandomEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  tmp = zero(u); rtmp = zero(rate_prototype)
  RandomEMCache(u,uprev,tmp,rtmp)
end

struct SimplifiedEMConstantCache <: StochasticDiffEqConstantCache end
@cache struct SimplifiedEMCache{randType,uType,rateType,rateNoiseType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  _dW::randType
  rtmp1::rateType
  rtmp2::rateNoiseType
end

alg_cache(alg::SimplifiedEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} = SimplifiedEMConstantCache()

function alg_cache(alg::SimplifiedEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
  else
    _dW = zero(ΔW)
  end

  rtmp1 = zero(rate_prototype)
  rtmp2 = zero(noise_rate_prototype)
  SimplifiedEMCache(u,uprev,_dW, rtmp1,rtmp2)
end

struct RKMilConstantCache <: StochasticDiffEqConstantCache end
@cache struct RKMilCache{uType,rateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  du1::rateType
  du2::rateType
  K::rateType
  tmp::uType
  L::rateType
end

alg_cache(alg::RKMil,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} = RKMilConstantCache()

function alg_cache(alg::RKMil,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  du1 = zero(rate_prototype); du2 = zero(rate_prototype)
  K = zero(rate_prototype); tmp = zero(u); L = zero(rate_prototype)
  RKMilCache(u,uprev,du1,du2,K,tmp,L)
end

struct RKMilCommuteConstantCache{WikType} <: StochasticDiffEqConstantCache
  WikJ::WikType
end
@cache struct RKMilCommuteCache{uType,rateType,rateNoiseType,WikType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  du1::rateType
  du2::rateType
  K::rateType
  gtmp::rateNoiseType
  L::rateNoiseType
  WikJ::WikType
  mil_correction::rateType
  Kj::uType
  Dgj::rateNoiseType
  tmp::uType
end

function alg_cache(alg::RKMilCommute,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} WikJ = WikJ = get_WikJ(ΔW,prob,alg)
  RKMilCommuteConstantCache{typeof(WikJ)}(WikJ)
end

function alg_cache(alg::RKMilCommute,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  du1 = zero(rate_prototype); du2 = zero(rate_prototype)
  K = zero(rate_prototype); gtmp = zero(noise_rate_prototype);
  L = zero(noise_rate_prototype); tmp = zero(rate_prototype)
  WikJ = get_WikJ(ΔW,prob,alg)
  mil_correction = zero(rate_prototype)
  Kj = zero(u); Dgj = zero(noise_rate_prototype)
  RKMilCommuteCache(u,uprev,du1,du2,K,gtmp,L,WikJ,mil_correction,Kj,Dgj,tmp)
end

struct RKMilGeneralConstantCache{WikType} <: StochasticDiffEqConstantCache
  WikJ::WikType
end

@cache struct RKMilGeneralCache{uType, rateType, rateNoiseType, WikType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  du₁::rateType
  du₂::rateType
  K::uType
  L::rateNoiseType
  mil_correction::uType
  ggprime::rateNoiseType
  WikJ::WikType
end

function alg_cache(alg::RKMilGeneral,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  WikJ = get_WikJ(ΔW,prob,alg)
  RKMilGeneralConstantCache{typeof(WikJ)}(WikJ)
end

function alg_cache(alg::RKMilGeneral,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  tmp = zero(u)
  du₁ = zero(rate_prototype)
  du₂ = zero(rate_prototype)
  K = zero(u)
  L = zero(noise_rate_prototype)
  mil_correction = zero(u)
  ggprime = zero(noise_rate_prototype)
  WikJ = get_WikJ(ΔW,prob,alg)
  RKMilGeneralCache{typeof(u), typeof(rate_prototype), typeof(noise_rate_prototype), typeof(WikJ)}(u, uprev, tmp, du₁, du₂, K, L, mil_correction, ggprime, WikJ)
end


text:
export getA
function computeFunc_ellipsoids!(m::Vector{Float64},k::Int64,Xc,h::Vector{Float64},n::Vector{Int64},
		u::Vector{Float64},xmin::Vector{Float64},beta::Float64)
numParamOfRBF = 10;
(x1,x2,x3) = getXt(k,m,numParamOfRBF);
i = round.(Int64,([x1;x2;x3] - xmin)./h .+ 0.5);
# ## here we take a reasonable box around the bf.
# ## The constant depends on the type of RBF that we take.
(a1,a2,a3,a4,a5,a6) = getA(k,m);
boxL = ceil.(Int64,0.9./(h*beta));
imax = min.(i + boxL,n);
imin = max.(i - boxL,[1;1;1]);
alpha = getAlpha(k,m,numParamOfRBF);

iishift2 = 0;
for l = imin[3]:imax[3]
	@inbounds iishift3 = (l-1)*n[1]*n[2];
	for j = imin[2]:imax[2]
		@inbounds iishift2 = iishift3 + (j-1)*n[1];
		for q = imin[1]:imax[1]
			ii = iishift2 + q;
			nX = 0.0;
			@inbounds y1 = Xc[ii,1]-x1;
			@inbounds y2 = Xc[ii,2]-x2;
			@inbounds y3 = Xc[ii,3]-x3;
			Ay1,Ay2,Ay3 = mulA(a1,a2,a3,a4,a5,a6,y1,y2,y3);
			y1*=Ay1;y2*=Ay2;y3*=Ay3;
			nX+=y1;nX+=y2;nX+=y3;
			if (nX <= 0.81) # 0.77~0.875^2
				argii = radiust(nX);
				argii = psi1(argii);
				argii *= alpha;
				@inbounds u[ii] += argii;
			end
		end
	end
end
end

function updateJacobian_ellipsoids!(k::Int64,m::Vector{Float64},dsu::Vector{Float64},h::Vector{Float64},n::Vector{Int64},Jbuilder::SpMatBuilder,u::Vector{Float64},iifunc::Function,Xc,xmin::Vector{Float64},beta::Float64)
numParamOfRBF = 10;
(x1,x2,x3) = getXt(k,m,numParamOfRBF);
alpha = getAlpha(k,m,numParamOfRBF);
(a1,a2,a3,a4,a5,a6) = getA(k,m);

i = round.(Int64,([x1;x2;x3] - xmin)./h .+ 0.5);
# ## here we take a reasonable box around the bf.
# ## The constant depends on the type of RBF that we take.
boxL = ceil.(Int64,0.9./(h*beta));
imax = min.(i + boxL,n);
imin = max.(i - boxL,[1;1;1]);

offset = convert(Int64,(k-1)*numParamOfRBF + 1);
md = 1e-3*maximum(abs.(dsu));
temp = 0.0;

iishift2 = 0;
for l = imin[3]:imax[3]
	@inbounds iishift3 = (l-1)*n[1]*n[2];
	for j = imin[2]:imax[2]
		@inbounds iishift2 = iishift3 + (j-1)*n[1];
		for q = imin[1]:imax[1]
			ii = iishift2 + q;
			temp = dsu[ii];
			if temp >= md
				nX = 0.0;
				@inbounds z1 = x1 - Xc[ii,1];
				@inbounds z2 = x2 - Xc[ii,2];
				@inbounds z3 = x3 - Xc[ii,3];
				(Az1,Az2,Az3) = mulA(a1,a2,a3,a4,a5,a6,z1,z2,z3);
				nX = z1*Az1 + z2*Az2 + z3*Az3;
				if (nX <= 0.81) # 0.77~0.875^2
					radii = radiust(nX);
					psi,dpsi = dpsi1_t(radii);
					psi*=temp;
					ii = iifunc(ii);
					setNext!(Jbuilder,ii,offset,psi);
					temp *= alpha;
					temp /= radii;
					temp *= dpsi;
					psi = z1; psi*=temp;
					setNext!(Jbuilder,ii,offset+1,0.5*z1*psi);
					setNext!(Jbuilder,ii,offset+2,z2*psi);
					setNext!(Jbuilder,ii,offset+3,z3*psi);
					psi = z2; psi*=temp;z2*=0.5;z2*=psi;
					setNext!(Jbuilder,ii,offset+4,z2);
					setNext!(Jbuilder,ii,offset+5,z3*psi);
					setNext!(Jbuilder,ii,offset+6,0.5*z3*z3*temp);
					Az1*=temp;Az2 *= temp; Az3 *= temp;
					setNext!(Jbuilder,ii,offset+7,Az1);
					setNext!(Jbuilder,ii,offset+8,Az2);
					setNext!(Jbuilder,ii,offset+9,Az3);
				end
			end
		end
	end
end
end

function getA(k::Int64,m::Vector{Float64})
	return (m[10*(k-1) + 2],m[10*(k-1) + 3],m[10*(k-1) + 4],m[10*(k-1) + 5],m[10*(k-1) + 6],m[10*(k-1) + 7]);
end
export idx_trilA
function idx_trilA()
	return [1;2;3;5;6;9];
end

## A:
# a1 a2 a3
# a2 a4 a5
# a3 a5 a6

function mulA(a1::Float64,a2::Float64,a3::Float64,a4::Float64,a5::Float64,a6::Float64,x1::Float64,x2::Float64,x3::Float64)
t = 0.0;
a1 *= x1; t = a2; t*=x2; a1+=t; t=a3; t*=x3; a1+=t;
a4 *= x2; a2*=x1; a4+=a2;  t = a5; t*=x3; a4+=t; ## a2 is now run over.
a6 *= x3; a3*=x1; a6 += a3; a5*=x2; a6 += a5;
return (a1,a4,a6);
end

# A = randn(3,3); A = A'*A + eye(3); #L = 0.5*(L+L');
# x = randn(3);
# dA = 0.005*randn(3,3); dA = 0.5*(dA+dA');

# function forward(x,A)
# return dot(x,A*x);
# end

# function der(x,A)
# C = (x*x');
# return C;
# end

# u0 = forward(x,A);
# for k=1:5
	# dA.*=0.5;
	# u = forward(x,A+dA);
	# ut = u0 + vec(der(x,A))'*vec(dA);
	# println("norm(ut-u0): ",abs(u-u0),", norm(ut - u0 - J0*dm): ",norm(u - ut));
# end

#########################################################################################################
#########################################################################################################
#########################################################################################################
# In the old code below we use a representation using a Cholesky factor.
# This code is still with Ihuge
# If no rotations are applied, using L frees us from the SPD regularization.

# function computeFuncAndUpdateIhuge_ellipsoids!(m::Vector{Float64},k::Int64,Xc::Array{Float64},h::Vector{Float64},n::Vector{Int64},
		# Ihuge::Vector{Int32},Istarts::Vector{Int64},u::Vector{Float64},xmin::Vector{Float64},ii_count_huge::Int64)
# numParamOfRBF = 10;
# (x1,x2,x3) = getXt(k,m,numParamOfRBF);
# i = round.(Int64,(getX(k,m,numParamOfRBF) - xmin)./h + 0.5);
# # ## here we take a reasonable box around the bf.
# # ## The constant depends on the type of RBF that we take.
# (l1,l2,l3,l4,l5,l6) = getL(k,m);
# L = [l1 0.0 0.0 ; l2 l4 0.0 ; l3 l5 l6];
# beta = 0.5*sqrt(maximum(eig(L*L')[1]));
# boxL = ceil.(Int32,0.85./(h*beta));
# imax = min.(i + boxL,n);
# imin = max.(i - boxL,[1;1;1]);
# alpha = getAlpha(k,m,numParamOfRBF);

# for l = imin[3]:imax[3]
	# iishift3 = (l-1)*n[1]*n[2];
	# for j = imin[2]:imax[2]
		# iishift2 = iishift3 + (j-1)*n[1];
		# for q = imin[1]:imax[1]
			# ii = iishift2 + q;
			# nX = 0.0;
			# @inbounds y1 = Xc[ii,1]-x1;
			# @inbounds y2 = Xc[ii,2]-x2;
			# @inbounds y3 = Xc[ii,3]-x3;
			# y1,y2,y3 = mulLT(l1,l2,l3,l4,l5,l6,y1,y2,y3);
			# y1*=y1;y2*=y2;y3*=y3;
			# nX+=y1;nX+=y2;nX+=y3;

			# if (nX <= 0.77) # 0.77~0.875^2
				# argii = radiust(nX);
				# argii = psi1(argii);
				# argii *= alpha;
				# @inbounds u[ii] += argii;
				# Ihuge[ii_count_huge] = ii;
				# ii_count_huge += 1;
			# end
		# end
	# end
# end
# Istarts[k+1] = ii_count_huge;
# return ii_count_huge;
# end

# function updateJacobianArrays_ellipsoids!(k::Int64,m::Vector{Float64},curr::Int64,dsu::Vector{Float64},Is::Vector{Int32},Js::Vector{Int32},Vs::Vector{Float64},u::Vector{Float64},Ihuge::Vector{Int32},Istarts::Vector{Int64},iifunc::Function,Xc::Array{Float64})
# numParamOfRBF = 10;
# (x1,x2,x3) = getXt(k,m,numParamOfRBF);
# alpha = getAlpha(k,m,numParamOfRBF);
# (l1,l2,l3,l4,l5,l6) = getL(k,m);

# offset = convert(Int32,(k-1)*numParamOfRBF + 1);
# temp = 0.0;
# for iih = (Istarts[k]):(Istarts[k+1]-1)
	# @inbounds ii = Ihuge[iih];
	# @inbounds temp = dsu[ii];
	# @inbounds z1 = x1 - Xc[ii,1];
	# @inbounds z2 = x2 - Xc[ii,2];
	# @inbounds z3 = x3 - Xc[ii,3];
	# (Ltz1,Ltz2,Ltz3) = mulLT(l1,l2,l3,l4,l5,l6,z1,z2,z3);
	# radii = Ltz1*Ltz1 + Ltz2*Ltz2 + Ltz3*Ltz3;
	# # radii = z1*z1 + z2*z2 + z3*z3;
	# radii = radiust(radii);
	# psi,dpsi = dpsi1_t(radii);
	# psi*=temp;
	# @inbounds Vs[curr] = psi;
	# temp *= alpha;
	# temp /= radii;
	# temp *= dpsi;
	# Ltz1*=temp; Ltz2*=temp; Ltz3*=temp;
	# @inbounds Vs[curr+1] = Ltz1*z1;
	# @inbounds Vs[curr+2] = Ltz1*z2;
	# @inbounds Vs[curr+3] = Ltz1*z3;
	# @inbounds Vs[curr+4] = Ltz2*z2;
	# @inbounds Vs[curr+5] = Ltz2*z3;
	# @inbounds Vs[curr+6] = Ltz3*z3;
	# (Az1, Az2, Az3) = mulL(l1,l2,l3,l4,l5,l6,Ltz1,Ltz2,Ltz3);
	# # Az1 = z1*temp;Az2 = z2*temp;Az3 = z3*temp;
	# @inbounds Vs[curr+7] = Az1;
	# @inbounds Vs[curr+8] = Az2;
	# @inbounds Vs[curr+9] = Az3;
	# ii = iifunc(ii);
	# for p = 0:(numParamOfRBF-1)
		# @inbounds Is[curr+p] = ii;
		# @inbounds Js[curr+p] = offset + Int32(p);
	# end
	# curr += numParamOfRBF;
# end
# return curr;
# end

# m is ordered as alpha, l1, l2, l3, l4, l5, l6, x1,x2,x3

## L:
# L1 0 0
# L2 L4 0
# L3 L5 L6

# function getL(k::Int64,m::Vector{Float64})
	# return (m[10*(k-1) + 2],m[10*(k-1) + 3],m[10*(k-1) + 4],m[10*(k-1) + 5],m[10*(k-1) + 6],m[10*(k-1) + 7]);
# end

# function mulL(L1::Float64,L2::Float64,L3::Float64,L4::Float64,L5::Float64,L6::Float64,x1::Float64,x2::Float64,x3::Float64)
# L1 *= x1;
# L4 *= x2; L2*=x1; L4+=L2;
# L3 *= x1; L5*=x2; L6*=x3; L6+=L5; L6+=L3;
# return (L1,L4,L6);
# end

# function mulLT(L1::Float64,L2::Float64,L3::Float64,L4::Float64,L5::Float64,L6::Float64,x1::Float64,x2::Float64,x3::Float64)
# L1 *= x1; L2*=x2; L3*=x3; L1 += L2; L1+=L3;
# L4 *= x2; L5*=x3; L4+=L5;
# L6*=x3;
# return (L1,L4,L6);
# end

# L = randn(3,3); #L = 0.5*(L+L');
# x = randn(3);
# dL = 0.5*randn(3,3); #dL = 0.5*(dL+dL');

# function forward(x,L)
# return dot(L'*x,L'*x);
# end

# function der(x,L)
# C = 2.0*(L'*x)*x';
# return C';
# end

# u0 = forward(x,L);
# for k=1:5
	# dL.*=0.5;
	# u = forward(x,L+dL);
	# ut = u0 + vec(der(x,L))'*vec(dL);
	# # ut = u0 + 2.0*dot(L'*x,dL'*x);
	# # ut = u0 + 2.0*trace(dL'*x*x'*L);
	# # ut = u0 + 2.0*vecdot(dL,x*x'*L);
	# println("norm(ut-u0): ",abs(u-u0),", norm(ut - u0 - J0*dm): ",norm(u - ut));
# end


text:
"""
	contribution(s,o)

Return the contribution of the given node or link in the Stream.
The contribution of a node v in a stream is defined as:
```math
n_v=\\frac{\\left| T_v \\right|}{\\left| T \\right|}
```
The contribution of a link uv in a stream is defined as:
```math
m_{uv}=\\frac{\\left| T_{uv} \\right|}{\\left| T \\right|}
```
Note: If the stream has no duration, this function returns 0.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function contribution(s::AbstractStream, o::StreamObject)
	duration(s)!=0 ? duration(o) / duration(s) : 0.0
end

"""
	contribution(ls,"node_name")

Return the contribution of the given node in the link stream.
By definition, this is always 1 since nodes have no dynamics in
link streams.
"""
contribution(ls::Union{LinkStream,DirectedLinkStream},node_name::AbstractString)=1.0

"""
	contribution(s,"node_name")

Return the contribution of the given node in the stream graph.
"""
function contribution(s::Union{StreamGraph,DirectedStreamGraph},node_name::AbstractString)
	contribution(s,s.W[node_name])
end

"""
	contribution(ls,"node_name",t)

Return the contribution of the given node at the given time in the link stream.
By definition, this is always 1 since nodes have no dynamics in
link streams.
"""
contribution(ls::Union{LinkStream,DirectedLinkStream},node_name::AbstractString,t::Float64)=1.0

"""
	contribution(s,"from","to")

Return the contribution of the directed link in the directed stream.
"""
function contribution(s::AbstractDirectedStream, from::AbstractString, to::AbstractString)
	((haskey(s.E,from))&&(haskey(s.E[from],to))) ? contribution(s,s.E[from][to]) : 0
end

"""
	contribution(s,"from","to")

Return the contribution of the undirected link in the undirected stream.
This means that contribution(s,"from","to")==contribution(s,"to","from")
"""
function contribution(s::AbstractUndirectedStream, from::AbstractString, to::AbstractString)
	if haskey(s.E,from) && haskey(s.E[from],to)
		contribution(s,s.E[from][to])
	elseif haskey(s.E,to) && haskey(s.E[to],from)
		contribution(s,s.E[to][from])
	else
		throw("Unknown link ($from,$to).")
	end
end

"""
	node_contribution(s,t)

Return the node contribution of the given time instant in the stream.
The node contribution of a time instant is defined as:
```math
k_t = \\frac{\\left| V_t \\right|}{\\left| V \\right|}
```
Note: If the stream has no node, this function returns 0.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function node_contribution(s::AbstractStream,t::Float64)
	length(s.V) != 0 ? length(nodes(s,t))/length(s.V) : 0.0
end

"""
	node_contribution(s,tc,t)

Return the node contribution of the given time instant in the stream
using a time cursor.
"""
function node_contribution(s::AbstractStream,tc::TimeCursor,t::Float64)
	length(s.V) != 0 ? length(nodes(tc,t))/length(s.V) : 0.0
end

"""
	link_contribution(s,t)

Return the link contribution of the given time instant in the stream.
The link contribution of a time instant is defined as:
```math
l_t = \\frac{\\left| E_t \\right|}{\\left| V \\otimes V \\right|}
```
Note: If the stream has no node, this function returns 0.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function link_contribution(s::AbstractStream,t::Float64)
	length(s.V) != 0 ? length(links(s,t))/length(s.V ⊗ s.V) : 0.0
end

"""
	link_contribution(s,tc,t)

Return the link contribution of the given time instant in the stream
using a time cursor.
"""
function link_contribution(s::AbstractStream,tc::TimeCursor,t::Float64)
	length(s.V) != 0 ? length(links(tc,t))/length(s.V ⊗ s.V) : 0.0
end


text:
export logistic_regression

import ..LogisticRegression

using mlpack._Internal.io

import mlpack_jll
const logistic_regressionLibrary = mlpack_jll.libmlpack_julia_logistic_regression

# Call the C binding of the mlpack logistic_regression binding.
function logistic_regression_mlpackMain()
  success = ccall((:logistic_regression, logistic_regressionLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module logistic_regression_internal
  import ..logistic_regressionLibrary

import ...LogisticRegression

# Get the value of a model pointer parameter of type LogisticRegression.
function IOGetParamLogisticRegression(paramName::String)::LogisticRegression
  LogisticRegression(ccall((:IO_GetParamLogisticRegressionPtr, logistic_regressionLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type LogisticRegression.
function IOSetParamLogisticRegression(paramName::String, model::LogisticRegression)
  ccall((:IO_SetParamLogisticRegressionPtr, logistic_regressionLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializeLogisticRegression(stream::IO, model::LogisticRegression)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializeLogisticRegressionPtr, logistic_regressionLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializeLogisticRegression(stream::IO)::LogisticRegression
  buffer = read(stream)
  LogisticRegression(ccall((:DeserializeLogisticRegressionPtr, logistic_regressionLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    logistic_regression(; [batch_size, decision_boundary, input_model, labels, lambda, max_iterations, optimizer, step_size, test, tolerance, training, verbose])

An implementation of L2-regularized logistic regression using either the L-BFGS
optimizer or SGD (stochastic gradient descent).  This solves the regression
problem

  y = (1 / 1 + e^-(X * b))

where y takes values 0 or 1.

This program allows loading a logistic regression model (via the `input_model`
parameter) or training a logistic regression model given training data
(specified with the `training` parameter), or both those things at once.  In
addition, this program allows classification on a test dataset (specified with
the `test` parameter) and the classification results may be saved with the
`predictions` output parameter. The trained logistic regression model may be
saved using the `output_model` output parameter.

The training data, if specified, may have class labels as its last dimension.
Alternately, the `labels` parameter may be used to specify a separate matrix of
labels.

When a model is being trained, there are many options.  L2 regularization (to
prevent overfitting) can be specified with the `lambda` option, and the
optimizer used to train the model can be specified with the `optimizer`
parameter.  Available options are 'sgd' (stochastic gradient descent) and
'lbfgs' (the L-BFGS optimizer).  There are also various parameters for the
optimizer; the `max_iterations` parameter specifies the maximum number of
allowed iterations, and the `tolerance` parameter specifies the tolerance for
convergence.  For the SGD optimizer, the `step_size` parameter controls the step
size taken at each iteration by the optimizer.  The batch size for SGD is
controlled with the `batch_size` parameter. If the objective function for your
data is oscillating between Inf and 0, the step size is probably too large.
There are more parameters for the optimizers, but the C++ interface must be used
to access these.

For SGD, an iteration refers to a single point. So to take a single pass over
the dataset with SGD, `max_iterations` should be set to the number of points in
the dataset.

Optionally, the model can be used to predict the responses for another matrix of
data points, if `test` is specified.  The `test` parameter can be specified
without the `training` parameter, so long as an existing logistic regression
model is given with the `input_model` parameter.  The output predictions from
the logistic regression model may be saved with the `predictions` parameter.

Note : The following parameters are deprecated and will be removed in mlpack 4:
`output`, `output_probabilities`
Use `predictions` instead of `output`
Use `probabilities` instead of `output_probabilities`

This implementation of logistic regression does not support the general
multi-class case but instead only the two-class case.  Any labels must be either
0 or 1.  For more classes, see the softmax_regression program.

As an example, to train a logistic regression model on the data '`data`' with
labels '`labels`' with L2 regularization of 0.1, saving the model to
'`lr_model`', the following command may be used:

```julia
julia> using CSV
julia> data = CSV.read("data.csv")
julia> labels = CSV.read("labels.csv"; type=Int)
julia> _, lr_model, _, _, _ = logistic_regression(labels=labels,
            lambda=0.1, training=data)
```

Then, to use that model to predict classes for the dataset '`test`', storing the
output predictions in '`predictions`', the following command may be used:

```julia
julia> using CSV
julia> test = CSV.read("test.csv")
julia> predictions, _, _, _, _ =
            logistic_regression(input_model=lr_model, test=test)
```

# Arguments

 - `batch_size::Int`: Batch size for SGD.  Default value `64`.

 - `decision_boundary::Float64`: Decision boundary for prediction; if the
      logistic function for a point is less than the boundary, the class is
      taken to be 0; otherwise, the class is 1.  Default value `0.5`.

 - `input_model::LogisticRegression`: Existing model (parameters).
 - `labels::Array{Int, 1}`: A matrix containing labels (0 or 1) for the
      points in the training set (y).
 - `lambda::Float64`: L2-regularization parameter for training.  Default
      value `0`.

 - `max_iterations::Int`: Maximum iterations for optimizer (0 indicates no
      limit).  Default value `10000`.

 - `optimizer::String`: Optimizer to use for training ('lbfgs' or 'sgd').
      Default value `lbfgs`.

 - `step_size::Float64`: Step size for SGD optimizer.  Default value
      `0.01`.

 - `test::Array{Float64, 2}`: Matrix containing test dataset.
 - `tolerance::Float64`: Convergence tolerance for optimizer.  Default
      value `1e-10`.

 - `training::Array{Float64, 2}`: A matrix containing the training set
      (the matrix of predictors, X).
 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `output::Array{Int, 1}`: If test data is specified, this matrix is
      where the predictions for the test set will be saved.
 - `output_model::LogisticRegression`: Output for trained logistic
      regression model.
 - `output_probabilities::Array{Float64, 2}`: If test data is specified,
      this matrix is where the class probabilities for the test set will be
      saved.
 - `predictions::Array{Int, 1}`: If test data is specified, this matrix is
      where the predictions for the test set will be saved.
 - `probabilities::Array{Float64, 2}`: If test data is specified, this
      matrix is where the class probabilities for the test set will be saved.

"""
function logistic_regression(;
                             batch_size::Union{Int, Missing} = missing,
                             decision_boundary::Union{Float64, Missing} = missing,
                             input_model::Union{LogisticRegression, Missing} = missing,
                             labels = missing,
                             lambda::Union{Float64, Missing} = missing,
                             max_iterations::Union{Int, Missing} = missing,
                             optimizer::Union{String, Missing} = missing,
                             step_size::Union{Float64, Missing} = missing,
                             test = missing,
                             tolerance::Union{Float64, Missing} = missing,
                             training = missing,
                             verbose::Union{Bool, Missing} = missing,
                             points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, logistic_regressionLibrary), Nothing, ());

  IORestoreSettings("L2-regularized Logistic Regression and Prediction")

  # Process each input argument before calling mlpackMain().
  if !ismissing(batch_size)
    IOSetParam("batch_size", convert(Int, batch_size))
  end
  if !ismissing(decision_boundary)
    IOSetParam("decision_boundary", convert(Float64, decision_boundary))
  end
  if !ismissing(input_model)
    logistic_regression_internal.IOSetParamLogisticRegression("input_model", convert(LogisticRegression, input_model))
  end
  if !ismissing(labels)
    IOSetParamURow("labels", labels)
  end
  if !ismissing(lambda)
    IOSetParam("lambda", convert(Float64, lambda))
  end
  if !ismissing(max_iterations)
    IOSetParam("max_iterations", convert(Int, max_iterations))
  end
  if !ismissing(optimizer)
    IOSetParam("optimizer", convert(String, optimizer))
  end
  if !ismissing(step_size)
    IOSetParam("step_size", convert(Float64, step_size))
  end
  if !ismissing(test)
    IOSetParamMat("test", test, points_are_rows)
  end
  if !ismissing(tolerance)
    IOSetParam("tolerance", convert(Float64, tolerance))
  end
  if !ismissing(training)
    IOSetParamMat("training", training, points_are_rows)
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("output")
  IOSetPassed("output_model")
  IOSetPassed("output_probabilities")
  IOSetPassed("predictions")
  IOSetPassed("probabilities")
  # Call the program.
  logistic_regression_mlpackMain()

  return IOGetParamURow("output"),
         logistic_regression_internal.IOGetParamLogisticRegression("output_model"),
         IOGetParamMat("output_probabilities", points_are_rows),
         IOGetParamURow("predictions"),
         IOGetParamMat("probabilities", points_are_rows)
end


text:
using WiltonInts84

using Base.Test
using FixedSizeArrays

include("num_quad.jl")

if !isdefined(:record)
    record = false
end

p1 = Vec(1.0, 0.0, 0.0)
p2 = Vec(0.0, 1.0, 0.0)
p3 = Vec(0.0, 0.0, 0.0)
z  = Vec(0.0,0.0,1.0)

N = 2
T = eltype(Float64)
P = Vec{3,T}
J = Vector{Vec{N+3,T}}()
L = Vector{Vec{N+3,P}}()
I = similar(J)
K = similar(L)

h = 1.2

# code 01 and 00
c = p3
r, R = 0.1, 1.1
c += Vec(0,0,h)
r, R = sqrt(r^2+h^2), sqrt(R^2+h^2)
q = contour(p1,p2,p3,c,r,R)
@test length(q.segments) == 3
@test length(q.arcs) == 1
@test length(q.circles) == 0
ctr = contour(p1,p2,p3,c,r,R)
A,B = wiltonints(ctr,c, Val{N})
push!(I,A); push!(K,B)
if record
    P,Q = dblquadints1(p1,p2,p3,c,Val{N},r,R)
    push!(J,P); push!(L,Q)
end
#error("stop")

# code 02
c = (p1+p2+p3)/3
r, R = 0.35, 200.0
c += Vec(0,0,h)
r, R = sqrt(r^2+h^2), sqrt(R^2+h^2)
q = contour(p1,p2,p3,c,r,R)
@test length(q.segments) == 6
@test length(q.arcs) == 3
@test length(q.circles) == 0
ctr = contour(p1,p2,p3,c,r,R)
A,B = wiltonints(ctr,c, Val{N})
push!(I,A); push!(K,B)
if record
    P,Q = dblquadints1(p1,p2,p3,c,Val{N},r,R)
    push!(J,P); push!(L,Q)
end

# code 10
c = (p1+p2+p3)/3
r, R = 0.1, 0.49
c += Vec(0,0,h)
r, R = sqrt(r^2+h^2), sqrt(R^2+h^2)
q = contour(p1,p2,p3,c,r,R)
@test length(q.segments) == 3
@test length(q.arcs) == 2
@test length(q.circles) == 1
ctr = contour(p1,p2,p3,c,r,R)
A,B = wiltonints(ctr,c, Val{N})
push!(I,A); push!(K,B)
if record
    P,Q = dblquadints1(p1,p2,p3,c,Val{N},r,R)
    push!(J,P); push!(L,Q)
end

# code 11
c = p3
r, R = 0.1, 0.3
c += Vec(0,0,h)
r, R = sqrt(r^2+h^2), sqrt(R^2+h^2)
q = contour(p1,p2,p3,c,r,R)
@test length(q.segments) == 2
@test length(q.arcs) == 2
@test length(q.circles) == 0
ctr = contour(p1,p2,p3,c,r,R)
A,B = wiltonints(ctr,c, Val{N})
push!(I,A); push!(K,B)
if record
    P,Q = dblquadints1(p1,p2,p3,c,Val{N},r,R)
    push!(J,P); push!(L,Q)
end

# code 12
c = Vec(0.25,0.1,0.0)
r, R = 0.2, 0.3
c += Vec(0,0,h)
r, R = sqrt(r^2+h^2), sqrt(R^2+h^2)
q = contour(p1,p2,p3,c,r,R)
@test length(q.segments) == 3
@test length(q.arcs) == 2
@test length(q.circles) == 0
ctr = contour(p1,p2,p3,c,r,R)
A,B = wiltonints(ctr,c, Val{N})
push!(I,A); push!(K,B)
if record
    P,Q = dblquadints1(p1,p2,p3,c,Val{N},r,R)
    push!(J,P); push!(L,Q)
end

c = Vec(0.1,0.25,0.0)
r, R = 0.2, 0.3
c += Vec(0,0,h)
r, R = sqrt(r^2+h^2), sqrt(R^2+h^2)
q = contour(p1,p2,p3,c,r,R)
@test length(q.segments) == 3
@test length(q.arcs) == 2
@test length(q.circles) == 0
ctr = contour(p1,p2,p3,c,r,R)
A,B = wiltonints(ctr,c, Val{N})
push!(I,A); push!(K,B)
if record
    P,Q = dblquadints1(p1,p2,p3,c,Val{N},r,R)
    push!(J,P); push!(L,Q)
end

c = Vec(0.15,0.15,0.0)
r, R = 0.2, 0.3
c += Vec(0,0,h)
r, R = sqrt(r^2+h^2), sqrt(R^2+h^2)
q = contour(p1,p2,p3,c,r,R)
@test length(q.segments) == 4
@test length(q.arcs) == 3
@test length(q.circles) == 0
ctr = contour(p1,p2,p3,c,r,R)
A,B = wiltonints(ctr,c, Val{N})
push!(I,A); push!(K,B)
if record
    P,Q = dblquadints1(p1,p2,p3,c,Val{N},r,R)
    push!(J,P); push!(L,Q)
end

# code 20
c = p3
r, R = 0.1, 0.9
c += Vec(0,0,h)
r, R = sqrt(r^2+h^2), sqrt(R^2+h^2)
q = contour(p1,p2,p3,c,r,R)
@test length(q.segments) == 3
@test length(q.arcs) == 3
@test length(q.circles) == 0
ctr = contour(p1,p2,p3,c,r,R)
A,B = wiltonints(ctr,c, Val{N})
push!(I,A); push!(K,B)
if record
    P,Q = dblquadints1(p1,p2,p3,c,Val{N},r,R)
    push!(J,P); push!(L,Q)
end

c = (p1+p2+p3)/3
r, R = 0.1, 0.45
c += Vec(0,0,h)
r, R = sqrt(r^2+h^2), sqrt(R^2+h^2)
q = contour(p1,p2,p3,c,r,R)
@test length(q.segments) == 3
@test length(q.arcs) == 3
@test length(q.circles) == 1
ctr = contour(p1,p2,p3,c,r,R)
A,B = wiltonints(ctr,c, Val{N})
push!(I,A); push!(K,B)
if record
    P,Q = dblquadints1(p1,p2,p3,c,Val{N},r,R)
    push!(J,P); push!(L,Q)
end

# code 21 is impossible

# code 22
c = (p1+p2+p3)/3
r, R = 0.35, 0.45
c += Vec(0,0,h)
r, R = sqrt(r^2+h^2), sqrt(R^2+h^2)
q = contour(p1,p2,p3,c,r,R)
@test length(q.segments) == 6
@test length(q.arcs) == 6
@test length(q.circles) == 0
ctr = contour(p1,p2,p3,c,r,R)
A,B = wiltonints(ctr,c, Val{N})
push!(I,A); push!(K,B)
if record
    P,Q = dblquadints1(p1,p2,p3,c,Val{N},r,R)
    push!(J,P); push!(L,Q)
end

# cover the remaining arc configurations
c = (p1+p2+p3)/3
r, R = 0.5, 0.6
c += Vec(0,0,h)
r, R = sqrt(r^2+h^2), sqrt(R^2+h^2)
q = contour(p3,p1,p2,c,r,R)
@test length(q.segments) == 4
@test length(q.arcs) == 4
@test length(q.circles) == 0
ctr = contour(p3,p1,p2,c,r,R)
A,B = wiltonints(ctr,c, Val{N})
push!(I,A); push!(K,B)
if record
    P,Q = dblquadints1(p3,p1,p2,c,Val{N},r,R)
    push!(J,P); push!(L,Q)
end

c = (p1+p2+p3)/3
r, R = 0.5, 0.6
c += Vec(0,0,h)
r, R = sqrt(r^2+h^2), sqrt(R^2+h^2)
q = contour(p2,p3,p1,c,r,R)
@test length(q.segments) == 4
@test length(q.arcs) == 4
@test length(q.circles) == 0
ctr = contour(p2,p3,p1,c,r,R)
A,B = wiltonints(ctr,c, Val{N})
push!(I,A); push!(K,B)
if record
    P,Q = dblquadints1(p2,p3,p1,c,Val{N},r,R)
    push!(J,P); push!(L,Q)
end

c = (p1+p2+p3)/3
r, R = 0.5, 0.6
c += Vec(0,0,h)
r, R = sqrt(r^2+h^2), sqrt(R^2+h^2)
q = contour(p3,p1,p2,c,r,R)
@test length(q.segments) == 4
@test length(q.arcs) == 4
@test length(q.circles) == 0
ctr = contour(p3,p1,p2,c,r,R)
A,B = wiltonints(ctr,c, Val{N})
push!(I,A); push!(K,B)
if record
    P,Q = dblquadints1(p3,p1,p2,c,Val{N},r,R)
    push!(J,P); push!(L,Q)
end

c = (p1+p2+p3)/3
r, R = 0.5, 0.6
c += Vec(0,0,h)
r, R = sqrt(r^2+h^2), sqrt(R^2+h^2)
q = contour(p1,p3,p2,c,r,R)
@test length(q.segments) == 4
@test length(q.arcs) == 4
@test length(q.circles) == 0
ctr = contour(p1,p3,p2,c,r,R)
A,B = wiltonints(ctr,c, Val{N})
push!(I,A); push!(K,B)
if record
    P,Q = dblquadints1(p1,p3,p2,c,Val{N},r,R)
    push!(J,P); push!(L,Q)
end

# remaining circle case
c = (p1+p2+p3)/3
r, R = 0.1, 0.2
c += Vec(0,0,h)
r, R = sqrt(r^2+h^2), sqrt(R^2+h^2)
q = contour(p1,p2,p3,c,r,R)
@test length(q.segments) == 0
@test length(q.arcs) == 0
@test length(q.circles) == 2
ctr = contour(p1,p3,p2,c,r,R)
A,B = wiltonints(ctr,c, Val{N})
push!(I,A); push!(K,B)
if record
    P,Q = dblquadints1(p1,p3,p2,c,Val{N},r,R)
    push!(J,P); push!(L,Q)
end

# singularity outside, h = 0
c = Vec(0.5,-0.1,0.0)
r, R = 0.2, 0.25
c += Vec(0,0,h)
r, R = sqrt(r^2+h^2), sqrt(R^2+h^2)
q = contour(p1,p2,p3,c,r,R)
@test length(q.segments) == 2
@test length(q.arcs) == 2
@test length(q.circles) == 0
ctr = contour(p1,p3,p2,c,r,R)
A,B = wiltonints(ctr,c, Val{N})
push!(I,A); push!(K,B)
if record
    P,Q = dblquadints1(p1,p3,p2,c,Val{N},r,R)
    push!(J,P); push!(L,Q)
end

# sphere tangent to triangle at vertex
c = Vec(0.0,1.0,0.1)
r, R = 0.0, 0.1
q = contour(p1,p2,p3,c,r,R)
@test length(q.segments) == 0
@test length(q.arcs) == 0
@test length(q.circles) == 0
A,B = wiltonints(q,c, Val{N})
push!(I,A); push!(K,B)
if record
    P,Q = dblquadints1(p1,p2,p3,c,Val{N},r,R)
    push!(J,P); push!(L,Q)
end

# circle through triangle vertex
c = Vec(-0.11,0.98,0.1)
r, R = 0.0, 0.15
q = contour(p1,p2,p3,c,r,R)
@test length(q.segments) == 1
@test length(q.arcs) == 1
@test length(q.circles) == 0
A,B = wiltonints(q,c, Val{N})
push!(I,A); push!(K,B)
if record
    P,Q = dblquadints1(p1,p2,p3,c,Val{N},r,R)
    push!(J,P); push!(L,Q)
end

# circle through triangle vertex
c = Vec(-0.11,-0.02,0.1)
r, R = 0.0, 0.15
q = contour(p1,p2,p3,c,r,R)
@test length(q.segments) == 0
@test length(q.arcs) == 0
@test length(q.circles) == 0
A,B = wiltonints(q,c, Val{N})
push!(I,A); push!(K,B)
if record
    P,Q = dblquadints1(p1,p2,p3,c,Val{N},r,R)
    push!(J,P); push!(L,Q)
end

using JLD
fn = joinpath(dirname(@__FILE__),"dblquad1.jld")

# convert data to Julia data type
I = T[I[m][n] for m in eachindex(I), n in 1:length(eltype(I))]
K = T[K[m][n][p] for m in eachindex(K), n in 1:length(eltype(K)), p in 1:3]

if record == true
    J = T[J[m][n] for m in eachindex(J), n in 1:length(eltype(J))]
    L = T[L[m][n][p] for m in eachindex(L), n in 1:length(eltype(L)), p in 1:3]
    jldopen(fn,"w") do file
        write(file, "J", J)
        write(file, "L", L)
    end
else
    J,L = jldopen(fn,"r") do file
        read(file, "J"), read(file, "L")
    end
end

@test maximum(I-J) < 1.0e-5
@test maximum(K-L) < 1.0e-5


text:
# https://www.hackerrank.com/challenges/30-binary-numbers/problem

function count_bin(n::Int)
    sn = string(n, base=2)
    consecutive_1 = Int[]
    count_1 = 0
    for i in 1:length(sn)
        if sn[i] == '1'
            count_1 += 1
            append!(consecutive_1, count_1)
        else
            count_1 = 0
            append!(consecutive_1, 0)
        end
    end
    return maximum(consecutive_1)
end

n = parse(Int, readline())
println(count_bin(n))

max_stars_repo_path: Practice/Tutorials/30 Days of Code/Day 10 Binary Numbers.jl

max_stars_repo_name: mtaufanr/hackerrank-code-julia

max_issues_repo_name: mtaufanr/hackerrank-code-julia

max_forks_repo_name: mtaufanr/hackerrank-code-julia


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -4.847398418942631 <= q <= 3.5461678827819525)

                     add_NL_constraint(m, :(1/(1 + exp(-(-0.19745953189550702 + -0.698583212677272*1/(1 + exp(-(-0.2024495957518173 + -0.1501654471281335*1/(1 + exp(-(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2])))) + 0.7965619538518913*1/(1 + exp(-(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2])))) + 0.8656395810315112*1/(1 + exp(-(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))))))) + 0.07925279857906231*1/(1 + exp(-(-0.4774715210476401 + -0.18966514690787006*1/(1 + exp(-(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2])))) + 0.1687755600264067*1/(1 + exp(-(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2])))) + 0.9882336248352352*1/(1 + exp(-(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))))))) + -0.8658353105986021*1/(1 + exp(-(-0.5967554189497073 + 0.7632066186026765*1/(1 + exp(-(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2])))) + -0.767873838532454*1/(1 + exp(-(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2])))) + 0.02410968604346264*1/(1 + exp(-(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2])))))))))) + 1/(1 + exp(-(0.9047173820312868 + 0.09375520961098927*1/(1 + exp(-(-0.2024495957518173 + -0.1501654471281335*1/(1 + exp(-(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2])))) + 0.7965619538518913*1/(1 + exp(-(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2])))) + 0.8656395810315112*1/(1 + exp(-(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))))))) + -0.9282166334997988*1/(1 + exp(-(-0.4774715210476401 + -0.18966514690787006*1/(1 + exp(-(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2])))) + 0.1687755600264067*1/(1 + exp(-(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2])))) + 0.9882336248352352*1/(1 + exp(-(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))))))) + 0.3928411969784107*1/(1 + exp(-(-0.5967554189497073 + 0.7632066186026765*1/(1 + exp(-(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2])))) + -0.767873838532454*1/(1 + exp(-(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2])))) + 0.02410968604346264*1/(1 + exp(-(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2])))))))))) + 1/(1 + exp(-(0.16193188153327842 + 0.4964335766578394*1/(1 + exp(-(-0.2024495957518173 + -0.1501654471281335*1/(1 + exp(-(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2])))) + 0.7965619538518913*1/(1 + exp(-(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2])))) + 0.8656395810315112*1/(1 + exp(-(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))))))) + -0.6165233713312808*1/(1 + exp(-(-0.4774715210476401 + -0.18966514690787006*1/(1 + exp(-(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2])))) + 0.1687755600264067*1/(1 + exp(-(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2])))) + 0.9882336248352352*1/(1 + exp(-(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))))))) + -0.4911670927879337*1/(1 + exp(-(-0.5967554189497073 + 0.7632066186026765*1/(1 + exp(-(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2])))) + -0.767873838532454*1/(1 + exp(-(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2])))) + 0.02410968604346264*1/(1 + exp(-(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2])))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module Mesh2d

export nverts, nedges

abstract Element
abstract Element2d <: Element

type Quad <: Element2d
    id::Int
    nv::Int
    ne::Int
    verts::Array{Float64,2}
    neigh::Array{Int,2}
    curve::Array{Int,1}
    curved::Bool
    hasbc::Bool
    orient::Array{Bool,1}
    Quad() = new(0, 4, 4, zeros(Float64,4,2), zeros(Int,4,2), zeros(Int,4),
                 false, false, [true, true, false, false])
end

type Tri <: Element2d
    id::Int
    nv::Int
    ne::Int
    verts::Array{Float64,2}
    neigh::Array{Int,2}
    curve::Array{Int,1}
    curved::Bool
    hasbc::Bool
    Tri() = new(0, 3, 3, zeros(Float64,3,2), zeros(Int,3,2), zeros(Int,3), false, false)
end

nverts(e::Quad) = 4
nedges(e::Quad) = 4
nverts(e::Tri) = 3
nedges(e::Tri) = 3

xverts(e::Element2d) = sub(e.verts, :, 1)
yverts(e::Element2d) = sub(e.verts, :, 2)

vertex_in_edge_quad = reshape([1,2,2,3,3,4,4,1], 2,4)
vertex_in_edge_tri = reshape([1,2,2,3,3,1], 2,3)
vinedge(el::Quad, edge) = sub(vertex_in_edge_quad, :, edge)
vinedge(edge) = sub(vertex_in_edge_quad, :, edge)

vinedge(el::Tri, edge) = sub(vertex_in_edge_tri, :, edge)

edge_in_vertex_quad = reshape([4,1,1,2,2,3,3,4], 2,4)
edge_in_vertex_tri = reshape([3,1,1,2,2,3], 2,3)
einvert(el::Quad, vert) = sub(edge_in_vertex_quad, :, vert)
einvert(vert) = sub(edge_in_vertex_quad, :, vert)

einvert(el::Tri, vert) = sub(edge_in_vertex_tri, :, vert)

abstract Bndry

abstract Curve
abstract Curve2d <: Curve

type Segment <: Curve2d
    seg::Int
    Segment() = new(0)
end

type Circle <: Curve2d
    xc::Float64
    yc::Float64
    rad::Float64
    inside::Bool
    Circle(x, y, r, s=false) = new(x, y, r, s)
end

type Mesh
    elems::Array{Element2d,1}
    curves::Array{Curve2d,1}
    bcs::Array{Bndry,1}
    other::Any
end

include("curve.jl")

end


text:
module LightNLP

using JLD2, FileIO

function save(filename, obj)
    #@save filename obj
end

function load(filename)
    #@load filename "obj"
end

include("intdict.jl")
#include("depparser/DepParser.jl")
include("ner/NER.jl")
include("lm/LM.jl")
#include("tokenization/Tokenization.jl")

export namedtuple

function namedtuple(config::Dict)
    k = map(Symbol, collect(keys(config)))
    NamedTuple{tuple(k...)}(tuple(values(config)...))
end

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule LERC_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("LERC")
JLLWrappers.@generate_main_file("LERC", UUID("88015f11-f218-50d7-93a8-a6af411a945d"))
end  # module LERC_jll


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

"""
    add_fano_noise(E_dep::RealQuantity, E_ionisation::RealQuantity, f_fano::Real)::RealQuantity

Add Fano noise to an energy deposition `E_dep`, assuming a detector material
ionisation energy `E_ionisation` and a Fano factor `f_fano`.
"""
function add_fano_noise end
export add_fano_noise

function add_fano_noise(rng::AbstractRNG, E_dep::RealQuantity, E_ionisation::RealQuantity, f_fano::Real)
    target_unit = unit(E_dep)
    n_expected = uconvert(Unitful.NoUnits, E_dep/E_ionisation)

    # dist = Distributions.Poisson(n_expected)
    # n_observed = n_expected + f_fano * (rand(dist) - n_expected)

    # Use normal distribution instead of Poisson, makes no differenc for
    # typical numbers of charges:
    σ = sqrt(n_expected)
    n_observed = n_expected + f_fano * σ * randn(rng)
    n_observed_pos = max(0, n_observed)

    uconvert(target_unit, E_ionisation * n_observed_pos)
end

function add_fano_noise(rng::AbstractRNG, events::DetectorHitEvents, E_ionisation::RealQuantity, f_fano::Real)
    noisy_edep = deepmap(x -> add_fano_noise(x, E_ionisation, f_fano), events.edep)

    TypedTables.Table(merge(
        Tables.columns(events),
        (edep = noisy_edep,)
    ))
end

add_fano_noise(x, E_ionisation::RealQuantity, f_fano::Real) =
    add_fano_noise(Random.GLOBAL_RNG, x, E_ionisation, f_fano)


text:
using StatsBase
using LinearAlgebra, Random, Test

@testset "StatsBase.Covariance" begin
weight_funcs = (weights, aweights, fweights, pweights)

@testset "$f" for f in weight_funcs
    X = randn(3, 8)

    Z1 = X .- mean(X, dims = 1)
    Z2 = X .- mean(X, dims = 2)

    w1 = rand(3)
    w2 = rand(8)

    # varcorrection is negative if sum of weights is smaller than 1
    if f === fweights
        w1[1] += 1
        w2[1] += 1
    end

    wv1 = f(w1)
    wv2 = f(w2)

    Z1w = X .- mean(X, wv1, 1)
    Z2w = X .- mean(X, wv2, 2)

    ## reference results

    S1 = Z1'Z1
    S2 = Z2 * Z2'

    Sz1 = X'X
    Sz2 = X * X'

    S1w = Z1w' * Matrix(Diagonal(w1)) * Z1w
    S2w = Z2w * Matrix(Diagonal(w2)) * Z2w'

    Sz1w = X' * Matrix(Diagonal(w1)) * X
    Sz2w = X * Matrix(Diagonal(w2)) * X'

    @testset "Scattermat" begin
        @test scattermat(X)    ≈ S1
        @test scattermat(X, 2) ≈ S2

        @test StatsBase.scattermatm(X, 0)    ≈ Sz1
        @test StatsBase.scattermatm(X, 0, 2) ≈ Sz2

        @test StatsBase.scattermatm(X, mean(X, dims = 1))    ≈ S1
        @test StatsBase.scattermatm(X, mean(X, dims = 2), 2) ≈ S2

        @test StatsBase.scattermatm(X, zeros(1,8))  ≈ Sz1
        @test StatsBase.scattermatm(X, zeros(3), 2) ≈ Sz2

        @testset "Weighted" begin
            @test scattermat(X, wv1)    ≈ S1w
            @test scattermat(X, wv2, 2) ≈ S2w

            @test StatsBase.scattermatm(X, 0, wv1)    ≈ Sz1w
            @test StatsBase.scattermatm(X, 0, wv2, 2) ≈ Sz2w

            @test StatsBase.scattermatm(X, mean(X, wv1, 1), wv1)    ≈ S1w
            @test StatsBase.scattermatm(X, mean(X, wv2, 2), wv2, 2) ≈ S2w

            @test StatsBase.scattermatm(X, zeros(1,8), wv1)  ≈ Sz1w
            @test StatsBase.scattermatm(X, zeros(3), wv2, 2) ≈ Sz2w
        end
    end

    @testset "Uncorrected" begin
        @testset "Weighted Covariance" begin
            @test cov(X, wv1; corrected=false)    ≈ S1w ./ sum(wv1)
            @test cov(X, wv2, 2; corrected=false) ≈ S2w ./ sum(wv2)

            @test StatsBase.covm(X, 0, wv1, 1; corrected=false) ≈ Sz1w ./ sum(wv1)
            @test StatsBase.covm(X, 0, wv2, 2; corrected=false) ≈ Sz2w ./ sum(wv2)

            @test StatsBase.covm(X, mean(X, wv1, 1), wv1, 1; corrected=false) ≈ S1w ./ sum(wv1)
            @test StatsBase.covm(X, mean(X, wv2, 2), wv2, 2; corrected=false) ≈ S2w ./ sum(wv2)

            @test StatsBase.covm(X, zeros(1,8), wv1, 1; corrected=false) ≈ Sz1w ./ sum(wv1)
            @test StatsBase.covm(X, zeros(3), wv2, 2; corrected=false)   ≈ Sz2w ./ sum(wv2)
        end

        @testset "Mean and covariance" begin
            (m, C) = mean_and_cov(X; corrected=false)
            @test m == mean(X, dims = 1)
            @test C == cov(X, dims = 1, corrected=false)

            (m, C) = mean_and_cov(X, 1; corrected=false)
            @test m == mean(X, dims = 1)
            @test C == cov(X, dims = 1, corrected = false)

            (m, C) = mean_and_cov(X, 2; corrected=false)
            @test m == mean(X, dims = 2)
            @test C == cov(X, dims = 2, corrected = false)

            (m, C) = mean_and_cov(X, wv1; corrected=false)
            @test m == mean(X, wv1, 1)
            @test C == cov(X, wv1, 1, corrected=false)

            (m, C) = mean_and_cov(X, wv1, 1; corrected=false)
            @test m == mean(X, wv1, 1)
            @test C == cov(X, wv1, 1, corrected=false)

            (m, C) = mean_and_cov(X, wv2, 2; corrected=false)
            @test m == mean(X, wv2, 2)
            @test C == cov(X, wv2, 2, corrected=false)
        end
        @testset "Conversions" begin
            std1 = std(X, wv1, 1; corrected=false)
            std2 = std(X, wv2, 2; corrected=false)

            cov1 = cov(X, wv1, 1; corrected=false)
            cov2 = cov(X, wv2, 2; corrected=false)

            cor1 = cor(X, wv1, 1)
            cor2 = cor(X, wv2, 2)

            @testset "cov2cor" begin
                @test cov2cor(cov(X, dims = 1), std(X, dims = 1)) ≈ cor(X, dims = 1)
                @test cov2cor(cov(X, dims = 2), std(X, dims = 2)) ≈ cor(X, dims = 2)
                @test cov2cor(cov1, std1) ≈ cor1
                @test cov2cor(cov2, std2) ≈ cor2
            end
            @testset "cor2cov" begin
                @test cor2cov(cor(X, dims = 1), std(X, dims = 1)) ≈ cov(X, dims = 1)
                @test cor2cov(cor(X, dims = 2), std(X, dims = 2)) ≈ cov(X, dims = 2)
                @test cor2cov(cor1, std1) ≈ cov1
                @test cor2cov(cor2, std2) ≈ cov2
            end
        end
    end

    @testset "Corrected" begin
        @testset "Weighted Covariance" begin
            if isa(wv1, Weights)
                @test_throws ArgumentError cov(X, wv1; corrected=true)
            else
                var_corr1 = StatsBase.varcorrection(wv1, true)
                var_corr2 = StatsBase.varcorrection(wv2, true)

                @test cov(X, wv1; corrected=true)    ≈ S1w .* var_corr1
                @test cov(X, wv2, 2; corrected=true) ≈ S2w .* var_corr2

                @test StatsBase.covm(X, 0, wv1, 1; corrected=true) ≈ Sz1w .* var_corr1
                @test StatsBase.covm(X, 0, wv2, 2; corrected=true) ≈ Sz2w .* var_corr2

                @test StatsBase.covm(X, mean(X, wv1, 1), wv1, 1; corrected=true) ≈ S1w .* var_corr1
                @test StatsBase.covm(X, mean(X, wv2, 2), wv2, 2; corrected=true) ≈ S2w .* var_corr2

                @test StatsBase.covm(X, zeros(1,8), wv1, 1; corrected=true) ≈ Sz1w .* var_corr1
                @test StatsBase.covm(X, zeros(3), wv2, 2; corrected=true)   ≈ Sz2w .* var_corr2
            end
        end
        @testset "Mean and covariance" begin
            (m, C) = mean_and_cov(X; corrected=true)
            @test m == mean(X, dims =1)
            @test C == cov(X, dims = 1, corrected = true)

            (m, C) = mean_and_cov(X, 1; corrected=true)
            @test m == mean(X, dims = 1)
            @test C == cov(X, dims = 1, corrected = true)

            (m, C) = mean_and_cov(X, 2; corrected=true)
            @test m == mean(X, dims = 2)
            @test C == cov(X, dims = 2, corrected = true)

            if isa(wv1, Weights)
                @test_throws ArgumentError mean_and_cov(X, wv1; corrected=true)
            else
                (m, C) = mean_and_cov(X, wv1; corrected=true)
                @test m == mean(X, wv1, 1)
                @test C == cov(X, wv1, 1; corrected=true)

                (m, C) = mean_and_cov(X, wv1, 1; corrected=true)
                @test m == mean(X, wv1, 1)
                @test C == cov(X, wv1, 1; corrected=true)

                (m, C) = mean_and_cov(X, wv2, 2; corrected=true)
                @test m == mean(X, wv2, 2)
                @test C == cov(X, wv2, 2; corrected=true)
            end
        end
        @testset "Conversions" begin
            if !isa(wv1, Weights)
                std1 = std(X, wv1, 1; corrected=true)
                std2 = std(X, wv2, 2; corrected=true)

                cov1 = cov(X, wv1, 1; corrected=true)
                cov2 = cov(X, wv2, 2; corrected=true)

                cor1 = cor(X, wv1, 1)
                cor2 = cor(X, wv2, 2)

                @testset "cov2cor" begin
                    @test cov2cor(cov(X, dims = 1), std(X, dims = 1)) ≈ cor(X, dims = 1)
                    @test cov2cor(cov(X, dims = 2), std(X, dims = 2)) ≈ cor(X, dims = 2)
                    @test cov2cor(cov1, std1) ≈ cor1
                    @test cov2cor(cov2, std2) ≈ cor2
                end

                @testset "cov2cor!" begin
                    tmp_cov1 = copy(cov1)
                    @test !(tmp_cov1 ≈ cor1)
                    StatsBase.cov2cor!(tmp_cov1, std1)
                    @test tmp_cov1 ≈ cor1

                    tmp_cov2 = copy(cov2)
                    @test !(tmp_cov2 ≈ cor2)
                    StatsBase.cov2cor!(tmp_cov2, std2)
                    @test tmp_cov2 ≈ cor2
                end

                @testset "cor2cov" begin
                    @test cor2cov(cor(X, dims = 1), std(X, dims = 1)) ≈ cov(X, dims = 1)
                    @test cor2cov(cor(X, dims = 2), std(X, dims = 2)) ≈ cov(X, dims = 2)
                    @test cor2cov(cor1, std1) ≈ cov1
                    @test cor2cov(cor2, std2) ≈ cov2
                end

                @testset "cor2cov!" begin
                    tmp_cor1 = copy(cor1)
                    @test !(tmp_cor1 ≈ cov1)
                    StatsBase.cor2cov!(tmp_cor1, std1)
                    @test tmp_cor1 ≈ cov1

                    tmp_cor2 = copy(cor2)
                    @test !(tmp_cor2 ≈ cov2)
                    StatsBase.cor2cov!(tmp_cor2, std2)
                    @test tmp_cor2 ≈ cov2
                end
            end
        end
    end

    @testset "Correlation" begin
        @test cor(X, f(ones(3)), 1) ≈ cor(X, dims = 1)
        @test cor(X, f(ones(8)), 2) ≈ cor(X, dims = 2)

        cov1 = cov(X, wv1, 1; corrected=false)
        std1 = std(X, wv1, 1; corrected=false)
        cov2 = cov(X, wv2, 2; corrected=false)
        std2 = std(X, wv2, 2; corrected=false)
        expected_cor1 = StatsBase.cov2cor!(cov1, std1)
        expected_cor2 = StatsBase.cov2cor!(cov2, std2)

        @test cor(X, wv1, 1) ≈ expected_cor1
        @test cor(X, wv2, 2) ≈ expected_cor2
    end
end
end # @testset "StatsBase.Covariance"


text:
# setup project paths -
const _PATH_TO_SRC = dirname(pathof(@__MODULE__))
const _PATH_TO_BASE = joinpath(_PATH_TO_SRC, "base")
const _PATH_TO_CONFIG = joinpath(_PATH_TO_SRC, "config")

# packages that I'm going to use ...
using Test
using TOML
using DataFrames
using CSV
using JSON
using BioSequences
using BioSymbols
using FASTX
using DelimitedFiles
using Logging
using BSON

# load my codes ...
include(joinpath(_PATH_TO_BASE, "VLTypes.jl"))
include(joinpath(_PATH_TO_BASE, "VLBase.jl"))
include(joinpath(_PATH_TO_BASE, "VLSequenceUtilities.jl"))
include(joinpath(_PATH_TO_BASE, "VLMetabolicUtilities.jl"))
include(joinpath(_PATH_TO_BASE, "VLFileUtilities.jl"))


text:
using Statistics
using Flux: onehotbatch, onecold, crossentropy, throttle
using Base.Iterators: repeated
#using CuArrays if you want to use GPU

imgs = Flux.Data.MNIST.images()
labels = Flux.Data.MNIST.labels();

X = hcat(float.(reshape.(imgs, :))...) #stack all the images
Y = onehotbatch(labels, 0:9) # just a common way to encode categorical variables

m = Chain(
  Dense(28^2, 32, relu),
  Dense(32, 10),
  softmax)

loss(x, y) = crossentropy(m(x), y)
opt = ADAM(); # popular stochastic gradient descent variant

accuracy(x, y) = mean(onecold(m(x)) .== onecold(y)) # cute way to find average of correct guesses

dataset = repeated((X,Y),200) # repeat the data set 200 times, as opposed to @epochs 200 ...
evalcb = () -> @show(loss(X, Y)) # callback to show loss

Flux.train!(loss, params(m), dataset, opt, cb = throttle(evalcb, 10)); #took me ~5 minutes to train on CPU

# Same preprocessing
test_X = hcat(float.(reshape.(Flux.Data.MNIST.images(:test), :))...)
test_Y = onehotbatch(Flux.Data.MNIST.labels(:test), 0:9);

m(test_X[:,5287]) # Note the 7th index ( corresponding to the digit 6 ) is nearly 1

#decode
onecold(m(test_X[:,5287])) - 1 #minus 1 since we start from 0, but indexing in Julia starts at 1

# Training set accuracy
accuracy(X, Y)

#######
using Flux, Metalhead
using Flux: @epochs
resnet = ResNet().layers

model = Chain(
  resnet[1:end-2],
  x -> reshape(x, 512),
  Dense(512, 13)
)

include("dataloader.jl")

dataset = [gpu.(load_batch(13*5)) for i in 1:10]

model = model |> gpu
opt = ADAM()

loss(x,y) = Flux.Losses.logitcrossentropy(model(x), y)

ps = Flux.params(model[2:end])

@epochs 2 Flux.train!(loss, ps, dataset, opt)

imgs, labels = gpu.(load_batch(13))
display(model(imgs))

labels

using Statistics
using Flux: onehotbatch, onecold, crossentropy, throttle
using Base.Iterators: repeated

accuracy(x, y) = mean(onecold(model(x)) .== onecold(y)) # cute way to find average of correct guesses

X = imgs
Y = labels

X1, Y1 = dataset[1][1], dataset[1][2]
X2, Y2 = dataset[2][1], dataset[2][2]

accuracy(X, Y)
accuracy(X1, Y1)
accuracy(X2, Y2)

max_stars_repo_name: zborffs/Delta

max_issues_repo_name: zborffs/Delta

max_forks_repo_name: zborffs/Delta


text:
@testset "entropy_thresholding" begin
    img = testimage("cameraman")
    edges, histogram = build_histogram(img,256)
    @test_throws ErrorException find_threshold(Entropy(),histogram,edges)
    @test find_threshold(Entropy(),histogram[1:256],edges) ≈ 0.76171875

    img = testimage("mandril_gray")
    edges, histogram = build_histogram(img,256)
    @test find_threshold(Entropy(),histogram[1:256],edges) ≈ 0.408517187461257

end


text:
import Base: show

"""

# Available top level Method

### Method
```julia
*  Sample::Method             : Sampling
*  Optimize::Method           : Optimization
*  Diagnose::Method           : Diagnostics
*  Variational::Method        : Variational Bayes
```
"""
abstract type Method end

const DataDict = Dict{String, Any}

mutable struct Random
  seed::Int64
end
Random(;seed::Number=-1) = Random(seed)

mutable struct Output
  file::String
  diagnostic_file::String
  refresh::Int64
end
Output(;file::String="", diagnostic_file::String="", refresh::Number=100) =
  Output(file, diagnostic_file, refresh)

mutable struct Stanmodel
  name::String
  nchains::Int
  num_warmup::Int
  num_samples::Int
  thin::Int
  id::Int
  model::String
  model_file::String
  monitors::Vector{String}
  data::Vector{DataDict}
  data_file::String
  command::Vector{Base.AbstractCmd}
  method::Method
  random::Random
  init::Vector{DataDict}
  init_file::String
  output::Output
  tmpdir::String
  output_format::Symbol
end

"""
# Method Stanmodel

Create a Stanmodel.

### Constructors
```julia
Stanmodel(
  method=Sample();
  name="noname",
  nchains=4,
  num_warmup=1000,
  num_samples=1000,
  thin=1,
  model="",
  monitors=String[],
  data=DataDict[],
  random=Random(),
  init=DataDict[],
  output=Output(),
  pdir::String=pwd(),
  output_format=:array
)

```
### Required arguments
```julia
* `method::Method`            : See ?Method
```

### Optional arguments
```julia
* `name::String`               : Name for the model
* `nchains::Int`               : Number of chains, if possible execute in parallel
* `num_warmup::Int`            : Number of samples used for num_warmupation
* `num_samples::Int`           : Sample iterations
* `thin::Int`                  : Stan thinning factor
* `model::String`              : Stan program source
* `data::DataDict[]`           : Observed input data as an array of Dicts
* `random::Random`             : Random seed settings
* `init::DataDict[]`           : Initial values for parameters in parameter block
* `output::Output`             : File output options
* `pdir::String`               : Working directory
* `monitors::String[] `        : Filter for variables used in Mamba post-processing
* `output_format::Symbol ` : Specifies the required output format (:array for CmdStan.jl)
```

### CmdStan.jl supports 2 output_format values:
```julia
1. :array                 # Returns an array of draws (default value)
2. :namedarray      # Returns a NamedArray

both return an Array{Float64, 3} with ndraws, nvars, nchains as indices. Other options are available
by `importing` or `using` packages such as StanDataFrames.jl, StanMamba.jl or StanMCMCChain.jl.
See `?CmdStan.convert_a3d`.
```

### Example
```julia
bernoullimodel = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
"

stanmodel = Stanmodel(num_samples=1200, thin=2, name="bernoulli", model=bernoullimodel);
```

### Related help
```julia
?stan                          : Run a Stanmodel
?Sample                        : Sampling settings
?Method                         : List of available methods
?Output                        : Output file settings
?DataDict                      : Input data dictionaries, will be converted to R syntax
?convert_a3d                   : Options for output formats
```
"""
function Stanmodel(
  method=Sample();
  name="noname",
  nchains=4,
  num_warmup=1000,
  num_samples=1000,
  thin=1,
  model="",
  monitors=String[],
  data=DataDict[],
  random=Random(),
  init=DataDict[],
  output=Output(),
  pdir::String=pwd(),
  output_format::Symbol=:array)

  println("=====> $pdir\n")
  cd(pdir)

  tmpdir = joinpath(pdir, "tmp")
  if !isdir(tmpdir)
    mkdir(tmpdir)
  end

  model_file = "$(name).stan"
  if length(model) > 0
    update_model_file(joinpath(tmpdir, "$(name).stan"), strip(model))
  end

  id::Int=0
  data_file::String=""
  init_file::String=""
  cmdarray = fill(``, nchains)

  if num_samples != 1000
    method.num_samples=num_samples
  end

  if num_warmup != 1000
    method.num_warmup=num_warmup
  end

  if thin != 1
    method.thin=thin
  end

  Stanmodel(name, nchains,
    num_warmup, num_samples, thin,
    id, model, model_file, monitors,
    data, data_file, cmdarray, method, random,
    init, init_file, output, tmpdir, output_format);
end

function model_show(io::IO, m::Stanmodel, compact::Bool)
  println("  name =                    \"$(m.name)\"")
  println("  nchains =                 $(m.nchains)")
  println("  num_samples =             $(m.num_samples)")
  println("  num_warmup =                $(m.num_warmup)")
  println("  thin =                    $(m.thin)")
  println("  monitors =                $(m.monitors)")
  println("  model_file =              \"$(m.model_file)\"")
  println("  data_file =               \"$(m.data_file)\"")
  println("  output =                  Output()")
  println("    file =                    \"$(m.output.file)\"")
  println("    diagnostics_file =        \"$(m.output.diagnostic_file)\"")
  println("    refresh =                 $(m.output.refresh)")
  println("  output_format =           :$(m.output_format)")
  if isa(m.method, Sample)
    sample_show(io, m.method, compact)
  elseif isa(m.method, Optimize)
    optimize_show(io, m.method, compact)
  elseif isa(m.method, Variational)
    variational_show(io, m.method, compact)
  else
    diagnose_show(io, m.method, compact)
  end
end

show(io::IO, m::Stanmodel) = model_show(io, m, false)


text:
# Subtypes of Stationary must define the following functions:
# metric(kernel::Stationary) = ::Metric
# cov(k::Stationary, r::Float64) = ::Float64
# grad_kern!

abstract Stationary <: Kernel
abstract StationaryData <: KernelData

distance(k::Stationary, X::Matrix{Float64}) = pairwise(metric(k), X)
distance(k::Stationary, X::Matrix{Float64}, Y::Matrix{Float64}) = pairwise(metric(k), X, Y)
distance(k::Stationary, x::Vector{Float64}, y::Vector{Float64}) = evaluate(metric(k), x, y)

cov(k::Stationary, x::Vector{Float64}, y::Vector{Float64}) = cov(k, distance(k, x, y))

function cov(k::Stationary, x1::Matrix{Float64}, x2::Matrix{Float64})
    nobsv1 = size(x1, 2)
    nobsv2 = size(x2, 2)
    R = distance(k, x1, x2)
    for i in 1:nobsv1, j in 1:nobsv2
        R[i,j] = cov(k, R[i,j])
    end
    return R
end

function cov(k::Stationary, X::Matrix{Float64}, data::StationaryData)
    nobsv = size(X, 2)
    R = distance(k, X, data)
    for i in 1:nobsv, j in 1:i
        @inbounds R[i,j] = cov(k, R[i,j])
        @inbounds R[j,i] = R[i,j]
    end
    return R
end

# Isotropic Kernels

abstract Isotropic <: Stationary

type IsotropicData <: StationaryData
    R::Matrix{Float64}
end

function KernelData(k::Isotropic, X::Matrix{Float64})
     IsotropicData(distance(k, X))
end

distance(k::Isotropic, X::Matrix{Float64}, data::IsotropicData) = copy(data.R)

# StationaryARD Kernels

abstract StationaryARD <: Stationary

type StationaryARDData <: StationaryData
    dist_stack::Array{Float64, 3}
end

# May need to customized in subtypes
function KernelData(k::StationaryARD, X::Matrix{Float64})
    d, nobsv = size(X)
    dist_stack = Array(Float64, nobsv, nobsv, d)
    for i in 1:d
        grad_ls = view(dist_stack, :, :, i)
        pairwise!(grad_ls, SqEuclidean(), X[i,:])
    end
    StationaryARDData(dist_stack)
end

function distance(k::StationaryARD, X::Matrix{Float64}, data::StationaryARDData)
    ### This commented section is slower than recalculating the distance from scratch...
    # nobsv = size(data.dist_stack,1)
    # d = length(k.ℓ2)
    # weighted = broadcast(/, data.dist_stack, reshape(k.ℓ2, (1,1,d)))
    # return reshape(sum(weighted, 3), (nobsv, nobsv))
    return pairwise(metric(k), X)
end


text:
###
# Specialized functions for Newsvendor with Ber(p) demand
# unit holding, backorder cost: b = s/(1-s)
# Because of Bernoulli demand, represent data as mhat = count zeta ==1
# xk = 1 iff p(alpha) > 1-s

# prob xk = 1
nv_prob_k(p, p0, alpha, Nhat_k, s) = ccdf(Binomial(Nhat_k, p), (Nhat_k + alpha) * (1 - s) - alpha * p0)

#vectorized way to simulate counts
#returns vector mhat_k for 1s only.
nv_sim_path(ps, Ns) = @.sum(rand(Bernoulli(ps), Ns))

#E[Z_k(alpha)], i.e. expected true obj at x(alpha)
function exp_nv_obj_k(p, p0, alpha, Nhat_k, s)
    b = s / (1 - s)
    (1 - p - b * p) * nv_prob_k(p, p0, alpha, Nhat_k, s) + b * p
end

#expected value of fullInfo solution
function exp_nv_fullInfo_k(p, s)
	b = s / (1 - s)
	xstar = p > 1 - s ? 1. : 0.
	b * p + (1 - p - b * p) * xstar
end

#Z_k(alpha)..  mhat is # zeta-k == 1
function nv_obj_k(mhat, p, p0, alpha, Nhat_k, s)
	p_alpha = shrink(mhat / Nhat_k, p0, alpha, Nhat_k)
	if p_alpha > 1 - s #xk = 1
		return 1 - p
	else #xk = 0
		return s / (1 - s) * p
	end
	return 0.  #never reached
end

#Computes Nint * ZLoo_k(alpha)
#Nint is the overall intensity for observations
#Assumes lambdas = 1
function nv_loo_k(mhat_k, p0, alpha, Nhat_k, s)
	out = 0.
	#i corresponding to demand = 1
	if mhat_k > 0
		p_alpha = shrink((mhat_k - 1)/(Nhat_k - 1), p0, alpha, Nhat_k - 1)
		#only non-zero if x neq 1
		out += p_alpha > 1 - s ? 0. : mhat_k * s / (1 - s)
	end

	#i corresponding to demand = 0
	if mhat_k < Nhat_k
		p_alpha = shrink(mhat_k / (Nhat_k - 1), p0, alpha, Nhat_k - 1)
		#only non-zero if x  = 1
		out += p_alpha > 1 - s ? (Nhat_k - mhat_k) : 0.
	end
	return out
end

#this is a funny saa-ish term in the Loo-Stability decomposition
#Returns Nint * saa-ish term.
function nv_saa_k(mhat_k, p0, alpha, Nhat_k, s)
	p_alpha = shrink(mhat_k/Nhat_k, p0, alpha, Nhat_k)
	if p_alpha > 1 - s #xk = 1
		return Nhat_k - mhat_k
	else #xk = 0
		return s / (1 - s) * mhat_k  #Note, not scaling by Nhat_k!
	end
	return 0.  #never reached
end

#Optimzes Z(\alpha) over alpha_grid
#returns alphaOR, minimizingAlphaIndex, curveInAlpha
function nv_oracle_alpha(mhats, ps, p0, alpha_grid, Nhats, ss)
	out = map(a-> nv_obj(mhats, ps, p0, a, Nhats, ss), alpha_grid)
	jstar = argmin(out)
	return alpha_grid[jstar], jstar, out
end

#returns p0OR, alphaOR, jOR for alpha
function nv_oracle_both(mhats, ps, p0_grid, alpha_grid, Nhats, ss)
	out = zeros(length(alpha_grid))
	or_perf, p0OR, alphaOR, jOR = Inf, 0., 0., 0
	for p0 in p0_grid
		out[:] = map(a-> nv_obj(mhats, ps, p0, a, Nhats, ss), alpha_grid)
		jstar = argmin(out)
		if out[jstar] < or_perf
			jOR = jstar
			or_perf = out[jstar]
			p0OR = p0
			alphaOR = alpha_grid[jstar]
		end
	end
	return p0OR, alphaOR, jOR
end

#optimizes ZLoo(\alpha)over alpha_grid
#Note curve in alpha = Nintensity * lamavg * ZLoo(alpha)
#returns alpha^LOO, maximizingAlphaIndex, curveInAlpha
function nv_loo_alpha(mhats, p0, alpha_grid, Nhats, ss)
	out = map(a-> nv_loo(mhats, p0, a, Nhats, ss), alpha_grid)
	jstar = argmin(out)
	return alpha_grid[jstar], jstar, out
end

#loo estimate optimizing boht p0 and alpha0
function nv_loo_both(mhats, p0_grid, alpha_grid, Nhats, ss)
	out = zeros(length(alpha_grid))
	loo_perf = Inf
	p0LOO = 0.
	alphaLOO = 0.
	jLOO = -1
	for p0 in p0_grid
		out[:] = map(a-> nv_loo(mhats, p0, a, Nhats, ss), alpha_grid)
		jstar = argmin(out)
		if out[jstar] < loo_perf
			jLOO = jstar
			loo_perf = out[jstar]
			p0LOO = p0
			alphaLOO = alpha_grid[jstar]
		end
	end
	return p0LOO, alphaLOO, jLOO
end

#Optimzes E[Z(\alpha)] over alpha_grid
#returns alphaAP, minimizingAlphaIndex, curveInAlpha
function nv_apriori_alpha(ps, p0, alpha_grid, Nhats, ss)
	out = map(a-> exp_nv_obj(ps, p0, a, Nhats, ss), alpha_grid)
	jstar = argmin(out)
	return alpha_grid[jstar], jstar, out
end

#vector versions assume lam_k = 1 for all k
nv_loo(mhats, p0, alpha, Nhats, ss) =  mean( nv_loo_k.(mhats, p0, alpha, Nhats, ss) )
nv_obj(mhats, ps, p0, alpha, Nhats, ss) = mean( nv_obj_k.(mhats, ps, p0, alpha, Nhats, ss) )
nv_saa(mhats, p0, alpha, Nhats, ss) = mean( nv_saa_k.(mhats, p0, alpha, Nhats, ss) )

exp_nv_obj(ps, p0, alpha, Nhats, ss) = mean( exp_nv_obj_k.(ps, p0, alpha, Nhats, ss) )
nv_var_obj(ps, p0, alpha, Nhats, ss) = mean( nv_var_obj_k.(ps, p0, alpha, Nhats, ss) )
exp_nv_fullInfo(ps, s) = mean(exp_nv_fullInfo_k.(ps, s))

max_issues_repo_name: CausalML/DataPooling

max_forks_repo_name: CausalML/DataPooling


text:
using Documenter

makedocs(
    # modules = [GUI],
    sitename = "GUI",
    pages  = Any[
        "Overview" => "index.md"
        "HTML" => "html.md"
        ]
    )

deploydocs(
    repo = "github.com/JeffreySarnoff/GUI.jl.git",
    target = "build"
)