text:
@model function model_v2(
    θ,
    predict,
	::Type{TV} = Vector{Float64},
	::Type{V}  = Float64;
) where {TV, V}

	@unpack num_observations, num_total_days, num_rt_steps, num_iar_steps, invlink, deathmodel, casemodel = θ
    # If we don't want to predict the future, we only need to compute up-to time-step `num_obs_countries[m]`
    num_time_steps   = predict ? num_total_days : num_observations
	############# 1.) time varying reproduction number

	R0         ~ truncated(Normal(3., 1.), 1.5, 5.)
	R1         ~ truncated(Normal(.8, .1), .5, 1.1)
	σ_rt       ~ truncated(Normal(0.15, .05), 0, .25)
	latent_Rt  ~ RandomWalk(num_rt_steps, σ_rt, invlink(R1))

	Rt = TV(undef, num_time_steps)

	random_walks!(Rt, θ, predict, latent_Rt, R0, σ_rt)

	############# 2.) time varying infection ascertainment rate
	σ_iar      ~ truncated(Normal(0.05, .03), 0, .15)
	iar0       ~ Beta(1,10)
	latent_iar ~ RandomWalk(num_iar_steps, σ_iar, logit(iar0))

	iar = TV(undef, num_time_steps)

	iar!(iar, θ, predict, latent_iar, iar0, σ_iar)

	############ 3.) infection dynamics
	τ  ~ Exponential(1 / 0.03) # `Exponential` has inverse parameterization of the one in Stan
	T  = typeof(τ)
	y  ~ truncated(Exponential(τ),T(1),T(1000))

	newly_infected       = TV(undef, num_time_steps)
	cumulative_infected  = TV(undef, num_time_steps)
	effective_Rt         = TV(undef, num_time_steps)

	# infections!(newly_infected, cumulative_infected, θ, τ, y, Rt)
	infections!(newly_infected, cumulative_infected, effective_Rt, θ, y, Rt)
	########### 4.) derive observables
	μ_i2d ~ truncated(Normal(20., 2.), 10, 25)
	μ_i2c ~ truncated(Normal(4.5, 1.), 2, 10)
	ifr   ~ truncated(Normal(6/1000, 3/1000), 2/1000, 10/1000)
	# ihr   ~ truncated(Normal(1/100,1/100),.1/100,5/100)

	### day of week effect
	weekdayeffect_simplex ~ Dirichlet([7,7,7,7,7,3,1])
	weekdayeffect = TV(undef, 7)
	weekdayeffect!(weekdayeffect, weekdayeffect_simplex)

	### holiday effect
	holidayeffect ~ Beta(1,1)

	########### 3.) observation model
	expected_daily_deaths = TV(undef, num_time_steps)
	expected_daily_cases  = TV(undef, num_time_steps)

	ϕ_c  ~ truncated(Normal(50, 10), 30, Inf)
	ϕ_d  ~ truncated(Normal(40, 10), 10, Inf)

    cases_observation_model = WeekdayHolidayObsModel(
        casemodel, μ_i2c, iar, ϕ_c, expected_daily_cases, weekdayeffect, holidayeffect
    )
    deaths_observation_model = SimpleObsModel(deathmodel, μ_i2d, ifr, ϕ_d, expected_daily_deaths)

	expected!(deaths_observation_model, newly_infected)
	expected!(cases_observation_model, newly_infected)

	# ## 4.3) seroprevalence study
	# serodelay = 16 # 14d self deferral + 16d seroconversion (see Erikstrup et al. https://doi.org/10.1093/cid/ciaa849)
	# expected_seropositive[1:serodelay]    .= zero(V)
	# expected_seropositive[serodelay+1:end] = cumulative_infected[1:end-serodelay]

	# select_seroprev = 2:4
	# seroprev_σ ~ arraydist( InverseGamma2.(seroprev_std[ select_seroprev ], Ref(.1)) )
    # seroprev_μ = expected_seropositive[ seroprev_idx[ select_seroprev ] ]
	# seroprev_mean[select_seroprev] ~ MvNormal(seroprev_μ, seroprev_σ) # select only early observations

	########### 4.) compare model to observations

	ℓ  = zero(V)
    ℓ += logpdf(cases_observation_model, θ.cases)
    ℓ += logpdf(deaths_observation_model, θ.deaths)
	Turing.@addlogprob! ℓ

	return (
		newly_infected        = newly_infected,
		expected_daily_deaths = expected_daily_deaths,
		expected_daily_cases  = expected_daily_cases,
		Rt                    = Rt,
		effective_Rt          = effective_Rt,
		cumulative_infected   = cumulative_infected,
		iar                   = iar
	)
end

# ================================================================================================
#                   Cases + Hospit + Deaths
# =================================================================================================
@model function model_v3(
    θ,
    predict,
	::Type{TV} = Vector{Float64},
	::Type{V}  = Float64;
) where {TV, V}

	@unpack num_observations, num_total_days, num_rt_steps, num_iar_steps, invlink, deathmodel, casemodel, hospitmodel = θ
    # If we don't want to predict the future, we only need to compute up-to time-step `num_obs_countries[m]`
    num_time_steps   = predict ? num_total_days : num_observations
	############# 1.) time varying reproduction number

	R0         ~ truncated(Normal(3., 1.), 1.5, 5.)
	R1         ~ truncated(Normal(.8, .1), .5, 1.1)
	σ_rt       ~ truncated(Normal(0.15, .05), 0, .25)
	latent_Rt  ~ RandomWalk(num_rt_steps, σ_rt, invlink(R1))

	Rt = TV(undef, num_time_steps)

	random_walks!(Rt, θ, predict, latent_Rt, R0, σ_rt)

	############# 2.) time varying infection ascertainment rate
	σ_iar      ~ truncated(Normal(0.05, .03), 0, .15)
	iar0       ~ Beta(1,10)
	latent_iar ~ RandomWalk(num_iar_steps, σ_iar, logit(iar0))

	iar = TV(undef, num_time_steps)

	iar!(iar, θ, predict, latent_iar, iar0, σ_iar)

	############ 3.) infection dynamics
	τ  ~ Exponential(1 / 0.03) # `Exponential` has inverse parameterization of the one in Stan
	T  = typeof(τ)
	y  ~ truncated(Exponential(τ),T(1),T(1000))

	newly_infected       = TV(undef, num_time_steps)
	cumulative_infected  = TV(undef, num_time_steps)
	effective_Rt         = TV(undef, num_time_steps)

	# infections!(newly_infected, cumulative_infected, θ, τ, y, Rt)
	infections!(newly_infected, cumulative_infected, effective_Rt, θ, y, Rt)
	########### 4.) derive observables
	μ_i2c ~ truncated(Normal(4.5, 1.), 2, 10)
	μ_i2h ~ truncated(Normal(12., 3.), 9, 16)
	μ_i2d ~ truncated(Normal(20., 2.), 10, 25)
	ihr   ~ truncated(Normal(1/100,1/100),.1/100,5/100)
	ifr   ~ truncated(Normal(6/1000, 3/1000), 2/1000, 10/1000)

	### day of week effect
	weekdayeffect_simplex ~ Dirichlet([7,7,7,7,7,3,1])
	weekdayeffect = TV(undef, 7)
	weekdayeffect!(weekdayeffect, weekdayeffect_simplex)

	### holiday effect
	holidayeffect ~ Beta(1,1)

	########### 3.) observation model
	expected_daily_cases  = TV(undef, num_time_steps)
	expected_daily_hospit = TV(undef, num_time_steps)
	expected_daily_deaths = TV(undef, num_time_steps)

	ϕ_c  ~ truncated(Normal(50, 10), 30, Inf)
	ϕ_h  ~ truncated(Normal(50, 10), 20, Inf)
	ϕ_d  ~ truncated(Normal(40, 10), 10, Inf)

    cases_observation_model = WeekdayHolidayObsModel(
        casemodel, μ_i2c, iar, ϕ_c, expected_daily_cases, weekdayeffect, holidayeffect
    )
	hospit_observation_model = SimpleObsModel(hospitmodel, μ_i2h, ihr, ϕ_h, expected_daily_hospit)
    deaths_observation_model = SimpleObsModel(deathmodel, μ_i2d, ifr, ϕ_d, expected_daily_deaths)

	expected!(cases_observation_model, newly_infected)
	expected!(hospit_observation_model, newly_infected)
	expected!(deaths_observation_model, newly_infected)

	# ## 4.3) seroprevalence study
	# serodelay = 16 # 14d self deferral + 16d seroconversion (see Erikstrup et al. https://doi.org/10.1093/cid/ciaa849)
	# expected_seropositive[1:serodelay]    .= zero(V)
	# expected_seropositive[serodelay+1:end] = cumulative_infected[1:end-serodelay]

	# select_seroprev = 2:4
	# seroprev_σ ~ arraydist( InverseGamma2.(seroprev_std[ select_seroprev ], Ref(.1)) )
    # seroprev_μ = expected_seropositive[ seroprev_idx[ select_seroprev ] ]
	# seroprev_mean[select_seroprev] ~ MvNormal(seroprev_μ, seroprev_σ) # select only early observations

	########### 4.) compare model to observations

	ℓ  = zero(V)
    ℓ += logpdf(cases_observation_model, θ.cases)
    ℓ += logpdf(hospit_observation_model, θ.hospit)
	ℓ += logpdf(deaths_observation_model, θ.deaths)
	Turing.@addlogprob! ℓ

	return (
		newly_infected        = newly_infected,
		expected_daily_cases  = expected_daily_cases,
		expected_daily_hospit = expected_daily_hospit,
		expected_daily_deaths = expected_daily_deaths,
		Rt                    = Rt,
		effective_Rt          = effective_Rt,
		cumulative_infected   = cumulative_infected,
		iar                   = iar
	)
end
# =============================================================================
#                      model_cases
# ==============================================================================
@model function model_cases(
    θ,
    predict,
	::Type{TV} = Vector{Float64},
	::Type{V}  = Float64;
) where {TV, V}

	@unpack num_observations, num_total_days, num_rt_steps, invlink, deathmodel, casemodel, seromodel, num_covariates = θ
    # If we don't want to predict the future, we only need to compute up-to time-step `num_obs_countries[m]`
    num_time_steps   = predict ? num_total_days : num_observations
	num_time_steps_2 = predict ? num_total_days : deathmodel.stop
	############# 2.) time varying reproduction number

	R0         ~ truncated(Normal(3., 1.), 1.5, 5.)
	R1         ~ truncated(Normal(.8, .1), .5, 1.1)
	σ_rt       ~ truncated(Normal(0.15, .05), 0, .25)
	latent_Rt  ~ RandomWalk(num_rt_steps, σ_rt, invlink(R1))

	Rt = TV(undef, num_time_steps)

	if num_covariates > 0
		effect ~ filldist(Exponential(0.2), num_covariates)
		if θ.semiparametric
			semiparametric!(Rt, θ, predict, latent_Rt, R0, σ_rt, effect)
		else
			mixed!(Rt, θ, predict, latent_Rt, R0, σ_rt, effect)
		end
	else
		random_walks!(Rt, θ, predict, latent_Rt, R0, σ_rt)
	end

	############ 3.) infection dynamics
	τ  ~ Exponential(1 / 0.03) # `Exponential` has inverse parameterization of the one in Stan
	T  = typeof(τ)
	y  ~ truncated(Exponential(τ),T(1),T(1000))

	newly_infected       = TV(undef, num_time_steps)
	cumulative_infected  = TV(undef, num_time_steps)
	effective_Rt         = TV(undef, num_time_steps)

	# infections!(newly_infected, cumulative_infected, θ, τ, y, Rt)
	infections!(newly_infected, cumulative_infected, effective_Rt, θ, y, Rt)
	########### 4.) derive observables
	#1
	# μ_i2d ~ truncated(Normal(20., 4.), 10, 30)
	# μ_i2c ~ truncated(Normal(4.5, 1.), 2, 10)
	#2
	μ_i2d ~ truncated(Normal(20., 2.), 10, 25)
	μ_i2c ~ truncated(Normal(2., 1.), 0, 5)
	iar   ~ Beta(1,10)
	ifr   ~ truncated(Normal(6/1000, 3/1000), 2/1000, 10/1000)
	# ihr   ~ truncated(Normal(1/100,1/100),.1/100,5/100)

	### day of week effect
	weekdayeffect_simplex ~ Dirichlet([7,7,7,7,7,3,1])
	weekdayeffect = TV(undef, 7)
	weekdayeffect!(weekdayeffect, weekdayeffect_simplex)

	### holiday effect
	holidayeffect ~ Beta(1,1)

	########### 3.) observation model
	expected_daily_deaths = TV(undef, num_time_steps_2)
	expected_daily_cases  = TV(undef, num_time_steps)

	ϕ_c  ~ truncated(InverseGamma2(80, .2), 30, Inf)
	ϕ_d  ~ truncated(Normal(40, 10), 10, Inf)

    cases_observation_model = WeekdayHolidayObsModel(
        casemodel, μ_i2c, iar, ϕ_c, expected_daily_cases, weekdayeffect, holidayeffect
    )
    deaths_observation_model = SimpleObsModel(deathmodel, μ_i2d, ifr, ϕ_d, expected_daily_deaths)
	sero_observation_model   = SimpleSeroObsModel2(seromodel, cumulative_infected)

	expected!(cases_observation_model, newly_infected)
	expected!(deaths_observation_model, newly_infected)

	########### 4.) compare model to observations

	ℓ  = zero(V)
    ℓ += logpdf(cases_observation_model, θ.cases)
    ℓ += logpdf(deaths_observation_model, θ.deaths)
	ℓ += logpdf(sero_observation_model, θ.sero)

	Turing.@addlogprob! ℓ

	return (
		newly_infected = newly_infected,
		expected_daily_deaths = expected_daily_deaths,
		expected_daily_cases = expected_daily_cases,
		Rt = Rt,
		effective_Rt = effective_Rt,
		expected_seropos = cumulative_infected
	)
end

# =============================================================================
#                      model_hospit
# ==============================================================================
@model function model_hospit(
    θ,
    predict,
	::Type{TV} = Vector{Float64},
	::Type{V}  = Float64;
) where {TV, V}

	@unpack num_observations, num_total_days, num_rt_steps, invlink, hospitmodel, seromodel, num_covariates = θ
    # If we don't want to predict the future, we only need to compute up-to time-step `num_obs_countries[m]`
    num_time_steps = predict ? num_total_days : num_observations

	############# 2.) time varying reproduction number
	R0         ~ truncated(Normal(3., 1.), 1.5, 5.)
	R1         ~ truncated(Normal(.8, .1), .5, 1.1)
	σ_rt       ~ truncated(Normal(0.15, .05), 0, .25)
	latent_Rt  ~ RandomWalk(num_rt_steps, σ_rt, invlink(R1))
	Rt = TV(undef, num_time_steps)
	if num_covariates > 0
		effect ~ filldist(Exponential(0.2), num_covariates)
		if θ.semiparametric
			semiparametric!(Rt, θ, predict, latent_Rt, R0, σ_rt, effect)
		else
			mixed!(Rt, θ, predict, latent_Rt, R0, σ_rt, effect)
		end
	else
		random_walks!(Rt, θ, predict, latent_Rt, R0, σ_rt)
	end

	############ 3.) infection dynamics
	τ  ~ Exponential(1 / 0.03) # `Exponential` has inverse parameterization of the one in Stan
	T  = typeof(τ)
	y  ~ truncated(Exponential(τ),T(1),T(1000))

	newly_infected       = TV(undef, num_time_steps)
	cumulative_infected  = TV(undef, num_time_steps)
	effective_Rt         = TV(undef, num_time_steps)

	# infections!(newly_infected, cumulative_infected, θ, τ, y, Rt)
	infections!(newly_infected, cumulative_infected, effective_Rt, θ, y, Rt)

	########### 4.) derive observables
	#1
    # μ ~ truncated(Normal(12., 1.), 9, 16)
	#2
	μ ~ truncated(Normal(9., 1.), 7, 15)
	α ~ truncated(Normal(1.8/100,0.5/100),1/100,5/100)
	ϕ ~ truncated(Normal(50, 10), 20, Inf)
	# σ_s = seromodel.std

	expected = TV(undef, num_time_steps)
	hospit_observation_model = SimpleObsModel(hospitmodel, μ, α, ϕ, expected)
	sero_observation_model   = SimpleSeroObsModel2(seromodel, cumulative_infected)

	expected!(hospit_observation_model, newly_infected)

	########### 4.) compare model to observations
	## 4.1) observe hospitalizations
	ℓ  = zero(V)
	ℓ += logpdf(hospit_observation_model, θ.hospit)
	ℓ += logpdf(sero_observation_model, θ.sero)
	Turing.@addlogprob! ℓ

	return (
		newly_infected = newly_infected,
		expected_daily_hospit = expected,
		Rt = Rt,
		effective_Rt = effective_Rt,
		expected_seropos = cumulative_infected
	)
end

# =============================================================================
#                      model_deaths
# ==============================================================================
@model function model_deaths(
    θ,
    predict,
	::Type{TV} = Vector{Float64},
	::Type{V}  = Float64;
) where {TV, V}

	@unpack num_observations, num_total_days, num_rt_steps, invlink, deathmodel, seromodel, num_covariates = θ
    # If we don't want to predict the future, we only need to compute up-to time-step `num_obs_countries[m]`
    num_time_steps = predict ? num_total_days : num_observations

	############# 2.) time varying reproduction number
	R0         ~ truncated(Normal(3., 1.), 1.5, 5.)
	R1         ~ truncated(Normal(.8, .1), .5, 1.1)
	σ_rt       ~ truncated(Normal(0.15, .05), 0, .25)
	latent_Rt  ~ RandomWalk(num_rt_steps, σ_rt, invlink(R1))

	Rt = TV(undef, num_time_steps)

	if num_covariates > 0
		effect ~ filldist(Exponential(0.2), num_covariates)
		if θ.semiparametric
			semiparametric!(Rt, θ, predict, latent_Rt, R0, σ_rt, effect)
		else
			mixed!(Rt, θ, predict, latent_Rt, R0, σ_rt, effect)
		end
	else
		random_walks!(Rt, θ, predict, latent_Rt, R0, σ_rt)
	end

	############ 3.) infection dynamics
	τ  ~ Exponential(1 / 0.03) # `Exponential` has inverse parameterization of the one in Stan
	T  = typeof(τ)
	y  ~ truncated(Exponential(τ),T(1),T(1000))

	newly_infected       = TV(undef, num_time_steps)
	cumulative_infected  = TV(undef, num_time_steps)
	effective_Rt         = TV(undef, num_time_steps)

	# infections!(newly_infected, cumulative_infected, θ, τ, y, Rt)
	infections!(newly_infected, cumulative_infected, effective_Rt, θ, y, Rt)

	########### 4.) derive observables
	μ ~ truncated(Normal(20., 2.), 10, 25)
	α ~ truncated(Normal(6/1000, 3/1000), 1/1000, 10/1000)
	ϕ ~ truncated(Normal(40, 10), 10, Inf)

	expected = TV(undef, num_time_steps)
	deaths_observation_model = SimpleObsModel(deathmodel, μ, α, ϕ, expected)
	sero_observation_model   = SimpleSeroObsModel2(seromodel, cumulative_infected)
	expected!(deaths_observation_model, newly_infected)

	########### 4.) compare model to observations
	## 4.1) observe deaths
	ℓ  = zero(V)
	ℓ += logpdf(deaths_observation_model, θ.deaths)
	ℓ += logpdf(sero_observation_model, θ.sero)
	Turing.@addlogprob! ℓ

	return (
		newly_infected = newly_infected,
		expected_daily_deaths = expected,
		Rt = Rt,
		effective_Rt = effective_Rt,
		expected_seropos = cumulative_infected
	)
end

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
module LBblocks

import MacroTools
export @bblock, @lblock, @sblock

# check if m.s is a declared function or type (necessarily const)
function isdeclaredfunction_or_type(m::Module, s::Symbol)
  m_s = getfield(m, s)
  isstruct_in_main = isa(m_s,Type) && isstructtype(m_s)
  return isconst(m,s) && (isa(m_s,Function) || isstruct_in_main) && nameof(m_s)==s
end

# check if m.s is a module function (necessarily const)
function ismod(m::Module, s::Symbol)
  m_s = getfield(m, s)
  return isconst(m,s) && isa(m_s,Module)
end

# Note: we partition
# names(m::Module) = Ω₁ ∪ Ω₂ ∪ Ω₃ where
# Ω₁ = {non const}
# Ω₂ = {consts which are delcared functions or types or Modules}
# Ω₃ = {consts which are *not* delcared functions or types or Modules}

#  Ω₁ = {non const}
function non_const_vals(m::Module)
    return [s for s in names(m) if !isconst(m,s)]
end

# Ω₁ᶜ
function const_vals(m::Module)
    return [s for s in names(m) if isconst(m,s)]
end

# Ω₂  = {consts which are delcared functions or types or Modules}
function const_delcaredfun_or_mod(m::Module)
    return [s for s in const_vals(m) if (ismod(m,s) || isdeclaredfunction_or_type(m,s))]
end

# Ω₃ = {consts which are *not* delcared functions or types or Modules}
function const_other(m::Module)
    return [s for s in const_vals(m) if !(ismod(m,s) || isdeclaredfunction_or_type(m,s))]
end

lhs(ex::Symbol) = ex
lhs(ex::Expr)   = MacroTools.isexpr(ex, :(=)) ? ex.args[1] : nothing

rhs(ex::Symbol) = ex
rhs(ex::Expr)   = MacroTools.isexpr(ex, :(=)) ? ex.args[end] : nothing

isassigned_vars(ex::Symbol)::Bool = false
isassigned_vars(ex::Expr)::Bool = MacroTools.isexpr(ex, :(=)) ? true : false
# ... which can be used as follows
# let_assigned_expr    = filter(isassigned_vars, [:(a=2), :y, :(b=c), :x])
# let_isolated_symbols = filter(!isassigned_vars, [:(a=2), :y, :(b=c), :x])

"""
```
rtn = @bblock begin
    code body
end
```
Begin block decorator which return anon eval
"""
macro bblock(ex)
    @assert ex.head == :block
    :((() -> $(esc(ex)))())
end

"""
```
rtn = @lblock let args
    code body
end
```
Let block decorator for function eval of code blocks ... need to explicitly pass
all non-const variables as let args
"""
macro lblock(ex)
    @assert ex.head == :let

    # grab first line and body of let
    MacroTools.@capture(
        ex,
        let fline_
            body__
        end
    )

    # split the first line into an array of Symbols
    MacroTools.@capture(fline, vlines__)

    # args from first line of let
    lhs_args    = map(lhs, vlines)
    rhs_args    = map(rhs, vlines)

    # declare local if in (Ω₁ - lhs_args) where Ω₁ = {non const}
    declared_local = Any[s for s in non_const_vals(__module__) if !(s in lhs_args)]

    # @show lhs_args
    # @show rhs_args
    # @show declared_local

    # prepend the let body with :(local declared_local_not_in_arglist...)
    new_body = length(declared_local) == 0 ?
        Expr(:block, map(esc,body)...) :
        Expr(:block, esc(:(local $(declared_local...))) , map(esc,body)...)

    # put the modified body and first line of let back together as a function
    fnm = gensym()
    fnm_ex = MacroTools.combinedef(
        Dict{Symbol,Any}(
            :name        => fnm,
            :args        => map(esc,lhs_args),
            :kwargs      => Any[],
            :body        => new_body,
            :whereparams => ()
        )
    )

    # the return expression defines the function and then calls it on the let args
    return MacroTools.rmlines(quote
        $fnm_ex
        $fnm($(map(esc,rhs_args)...))
    end)

end

"""
```
rtn = @sblock let args
    code body
end
```
The most strict let block decorator ... need to explicitly pass
all variables as let args with the exception of modules and explicitly declared
(non-anon) functions
"""
macro sblock(ex)
    @assert ex.head == :let

    # grab first line and body of let
    MacroTools.@capture(
        ex,
        let fline_
            body__
        end
    )

    # split the first line into an array of Symbols
    MacroTools.@capture(fline, vlines__)

    # args from first line of let
    lhs_args    = map(lhs, vlines)
    rhs_args    = map(rhs, vlines)

    # declare local if in (Ω₁ ∪ Ω₃ - lhs_args) where Ω₁ = {non const}
    # and Ω₃ = {consts which are *not* delcared functions or Modules}
    declared_local1  = Any[s for s in non_const_vals(__module__) if !(s in lhs_args)]
    declared_local2  = Any[s for s in const_other(__module__) if !(s in lhs_args)]
    declared_local   = vcat(declared_local1, declared_local2)

    # @show declared_local

    # prepend the let body with :(local declared_local_not_in_arglist...)
    new_body = length(declared_local) == 0 ?
        Expr(:block, map(esc,body)...) :
        Expr(:block, esc(:(local $(declared_local...))) , map(esc,body)...)

    # put the modified body and first line of let back together as a function
    fnm = gensym()
    fnm_ex = MacroTools.combinedef(
        Dict{Symbol,Any}(
            :name        => fnm,
            :args        => map(esc,lhs_args),
            :kwargs      => Any[],
            :body        => new_body,
            :whereparams => ()
        )
    )

    # the return expression defines the function and then calls it on the let args
    return MacroTools.rmlines(quote
        $fnm_ex
        $fnm($(map(esc,rhs_args)...))
    end)

end

end # module

ext: jl


text:
@testset "datalayout" begin

Context() do ctx
DataLayout("E-p:32:32-f128:128:128") do data
    @test convert(String, data) == "E-p:32:32-f128:128:128"

    @test occursin("E-p:32:32-f128:128:128", sprint(io->show(io,data)))

    @test byteorder(data) == LLVM.API.LLVMBigEndian
    @test pointersize(data) == pointersize(data, 0) == 4

    @test intptr(data) == intptr(data, 0) == LLVM.Int32Type()
    @test intptr(data, ctx) == intptr(data, 0, ctx) == LLVM.Int32Type(ctx)

    @test sizeof(data, LLVM.Int32Type()) == storage_size(data, LLVM.Int32Type()) == abi_size(data, LLVM.Int32Type()) == 4

    @test abi_alignment(data, LLVM.Int32Type()) == frame_alignment(data, LLVM.Int32Type()) == preferred_alignment(data, LLVM.Int32Type()) == 4

    LLVM.Module("SomeModule") do mod
        gv = GlobalVariable(mod, LLVM.Int32Type(), "SomeGlobal")
        @test preferred_alignment(data, gv) == 4

        datalayout!(mod, data)
        @test convert(String, datalayout(mod)) == convert(String, data)
    end

    elem = [LLVM.Int32Type(ctx), LLVM.FloatType(ctx)]
    let st = LLVM.StructType(elem, ctx)
        @test element_at(data, st, 4) == 1
        @test offsetof(data, st, 1) == 4
    end
end
end

end


text:
"""
    add(network::QNetwork, channel::QChannel)

Add a channel to the network
"""
function add(network::QNetwork, channel::QChannel)
    push!(network.channels, channel)
end

function update(channel::QChannel, old_time::Float64, new_time::Float64)
end

"""
    distance(src::QNode, dest::QNode)

Cartesian distance between two nodes
"""
function distance(src::QNode, dest::QNode)
    v = src.location
    w = dest.location
    return sqrt((v.x - w.x)^2 + (v.y - w.y)^2 + (v.z - w.z)^2)
end

"""
    mutable struct BasicChannel <: QChannel

The default Channel type. Costs are assumed to be exponential with length
"""
mutable struct BasicChannel <: QChannel
    name::String
    costs::Dict
    src::QNode
    dest::QNode
    length::Float64
    active::Bool
    directed::Bool

    """
    Initalise a Basic Channel with unit costs
    """
    function BasicChannel(src::QNode, dest::QNode, ;exp_cost::Bool=false)
        tmpchannel = new("", unit_costvector(), src, dest, distance(src, dest), true, false)
        if exp_cost == true
            tmpchannel.costs = cost(tmpchannel)
        end
        return tmpchannel
    end

    function BasicChannel(name::String, src::QNode, dest::QNode, exp_cost::Bool=false)
        tmpchannel = new(name, unit_costvector(), src, dest, distance(src, dest), true, false)
        if exp_cost == true
            tmpchannel.costs = cost(tmpchannel)
        end
        return tmpchannel
    end

    function BasicChannel(name::String, src::QNode, dest::QNode, edge_cost::Dict)
        tmpchannel = new(name, edge_cost, src, dest, distance(src, dest), true, false)
        return tmpchannel
    end

    #BasicChannel(src::QNode, dest::QNode) = new("", unit_costvector(), src, dest, distance(src, dest), true, false)
    #BasicChannel(name::String, src::QNode, dest::QNode) = new(string(name), unit_costvector(), src, dest, distance(src, dest), true, false)
end

function cost(channel::BasicChannel)
    d = channel.length
    β = 0.001
    loss = P_to_dB(exp(-β * d))
    Z = Z_to_dB((1 + exp(-β * d))/2)
    costVector = Dict([("loss",loss), ("Z",Z)])
    return costVector
end

mutable struct AirChannel <: QChannel
    name::String
    costs::Dict
    src::QNode
    dest::QNode
    length::Float64
    active::Bool
    directed::Bool

    function AirChannel(src::QNode, dest::QNode)
        tmpchannel = new("", unit_costvector(), src, dest, distance(src, dest), true, false)
        tmpchannel.costs = cost(tmpchannel)
        return tmpchannel
    end

    function AirChannel(name::String, src::QNode, dest::QNode)
        tmpchannel = new(name, unit_costvector(), src, dest, distance(src, dest), true, false)
        tmpchannel.costs = cost(tmpchannel)
        return tmpchannel
    end
end

function cost(channel::AirChannel)
    """
    Line integral for effective density

    / L'
    |     rho(x * sin(theta)) dx
    / 0
    """

    v = channel.src.location
    w = channel.dest.location
    L = channel.length
    # sin(theta)
    st = abs(v.z - w.z)/L
    # atmosphere function
    ρ = expatmosphere
    f(x) = ρ(x*st)
    # Effective atmospheric depth
    d = quadgk(f, 0, L)[1]

    β = 10e-5
    d₀ = 10e7

    # Calculate the decibelic forms of loss and Z
    P = exp(-β*d)*(d₀)^2/(d + d₀)^2
    Z = (1+exp(-β*d))/2

    # Put them in a cost vector and return
    costVector = Dict([("loss",P_to_dB(P)), ("Z",Z_to_dB(Z))])
    return costVector
end

function update(channel::AirChannel, old_time::Float64, new_time::Float64)
    channel.length = distance(channel.src, channel.dest)
    channel.costs = cost(channel)
end

#     BasicMemory(node::QNode) = new(node, unit_costvector())
#     BasicMemory(node::QNode, costs::Dict) = new(node, costs)
# end

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
# calculate required particle num for kld sampling

module ParticleNumKld
  using Plots, Distributions
  pyplot()

  function num(epsilon, delta, bin_num)
    return ceil(quantile(Chisq(bin_num-1), 1-delta)/(2*epsilon))
  end

  function main(is_test=false)
    bs = 2:1:9
    n = [num(0.1, 0.01, b) for b in bs]
    p1 = plot(bs, n, title="bin: 2-10")

    bs = 2:1:9999
    n = [num(0.1, 0.01, b) for b in bs]
    p2 = plot(bs, n, title="bin: 2-100000")

    plot(p1, p2, layout=(1, 2), size=(750, 300), legend=false)

    if is_test == false
      save_path = joinpath(split(@__FILE__, "src")[1], "src/localization/particle_filter/kld_sampling/particle_num_kld.png")
      savefig(save_path)
    end
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
# # Closest Correlation Matrix
#
# We consider the problem of finding the closest correlation matrix $X$ to a given random matrix $C$.
# With closest correlation matrix we mean a positive semidefinite matrix with ones on the diagonal.
# The problem is given by:
# $$
# \begin{array}{ll} \text{minimize} &  \frac{1}{2}||X - C||_F^2\\
# \text{subject to} &  X_{ii} = 1, \quad i=1,\dots,n \\
#                   &  X \succeq 0.
# \end{array}
# $$
# Notice that we use `JuMP` to model the problem. `COSMO` is chosen as the backend solver. And COSMO-specific settings are passed using
# the `optimizer_with_attributes()` function.

using COSMO, JuMP, LinearAlgebra, SparseArrays, Test, Random

#-
rng = Random.MersenneTwister(12345);
## create a random test matrix C
n = 8;
C = -1 .+ rand(rng, n, n) .* 2;
c = vec(C);

# Define problem in `JuMP`:
q = -vec(C);
r = 0.5 * vec(C)' * vec(C);
m = JuMP.Model(optimizer_with_attributes(COSMO.Optimizer, "verbose" => true));
@variable(m, X[1:n, 1:n], PSD);
x = vec(X);
@objective(m, Min, 0.5 * x' * x  + q' * x + r);
for i = 1:n
  @constraint(m, X[i, i] == 1.);
end

# Solve the `JuMP` model with `COSMO` and query the solution `X_sol`:
status = JuMP.optimize!(m);
obj_val = JuMP.objective_value(m);
X_sol = JuMP.value.(X);

# Double check result against known solution:
known_opt_val = 12.5406
known_solution =  [
  1.0         0.732562   -0.319491   -0.359985   -0.287543   -0.15578     0.0264044  -0.271438;
  0.732562    1.0         0.0913246  -0.0386357   0.299199   -0.122733    0.126612   -0.187489;
 -0.319491    0.0913246   1.0        -0.0863377   0.432948    0.461783   -0.248641   -0.395299;
 -0.359985   -0.0386357  -0.0863377   1.0         0.503379    0.250601    0.141151    0.286088;
 -0.287543    0.299199    0.432948    0.503379    1.0        -0.0875199   0.137518    0.0262425;
 -0.15578    -0.122733    0.461783    0.250601   -0.0875199   1.0        -0.731556    0.0841783;
  0.0264044   0.126612   -0.248641    0.141151    0.137518   -0.731556    1.0        -0.436274;
 -0.271438   -0.187489   -0.395299    0.286088    0.0262425   0.0841783  -0.436274    1.0  ];
@test isapprox(obj_val, known_opt_val , atol=1e-3)

#-
@test norm(X_sol - known_solution, Inf) < 1e-3


text:
module TestUtils

using Test
using TMLE
using MLJ
using StableRNGs
using Distributions

LinearBinaryClassifier = @load LinearBinaryClassifier pkg=GLM verbosity=0
LinearRegressor = @load LinearRegressor pkg=MLJLinearModels verbosity=0
LogisticClassifier = @load LogisticClassifier pkg=MLJLinearModels verbosity=0

@testset "Test expected_value & maybelogit" begin
    n = 100
    X = MLJ.table(rand(n, 3))

    # Probabilistic Classifier
    y = categorical(rand([0, 1], n))
    mach = machine(ConstantClassifier(), X, y)
    fit!(mach; verbosity=0)
    proba = mach.fitresult[2][2]
    ŷ = MLJ.predict(mach)
    expectation = TMLE.expected_value(ŷ, typeof(mach.model), target_scitype(mach.model))
    @test expectation == repeat([proba], n)
    @test TMLE.maybelogit(expectation, typeof(mach.model), target_scitype(mach.model)) == TMLE.logit(expectation)

    # Probabilistic Regressor
    y = rand(n)
    mach = machine(ConstantRegressor(), X, y)
    fit!(mach; verbosity=0)
    ŷ = MLJ.predict(mach)
    expectation = TMLE.expected_value(ŷ, typeof(mach.model), target_scitype(mach.model))
    @test expectation ≈ repeat([mean(y)], n) atol=1e-10
    @test TMLE.maybelogit(expectation, typeof(mach.model), target_scitype(mach.model)) == expectation

    # Deterministic Regressor
    mach = machine(LinearRegressor(), X, y)
    fit!(mach; verbosity=0)
    ŷ = MLJ.predict(mach)
    expectation = TMLE.expected_value(ŷ, typeof(mach.model), target_scitype(mach.model))
    @test expectation == ŷ
    @test TMLE.maybelogit(expectation, typeof(mach.model), target_scitype(mach.model)) == expectation

end

@testset "Test interaction_combinations" begin
    # With 1 treatment variable
    query = (t₁=["a", "b"],)
    combinations = TMLE.interaction_combinations(query)
    expected_comb = [(t₁ = "a",),
                     (t₁ = "b",)]
    for (i, comb) in enumerate(combinations)
        @test comb == expected_comb[i]
    end

    # With 3 treatment variables
    query = (t₁=["a", "b"], t₂ = ["c", "d"], t₃ = [true, false])
    combinations = TMLE.interaction_combinations(query)
    expected_comb = [(t₁ = "a", t₂ = "c", t₃ = true),
                     (t₁ = "b", t₂ = "c", t₃ = true),
                     (t₁ = "a", t₂ = "d", t₃ = true),
                     (t₁ = "b", t₂ = "d", t₃ = true),
                     (t₁ = "a", t₂ = "c", t₃ = false),
                     (t₁ = "b", t₂ = "c", t₃ = false),
                     (t₁ = "a", t₂ = "d", t₃ = false),
                     (t₁ = "b", t₂ = "d", t₃ = false)]
    for (i, comb) in enumerate(combinations)
        @test comb == expected_comb[i]
    end
end

@testset "Test indicator_fns" begin
    # For 1 treatment variable
    query = (t₁=["a", "b"],)
    indicators = TMLE.indicator_fns(query)
    @test indicators == Dict(
        (t₁ = "b",) => -1,
        (t₁ = "a",) => 1
    )

    # For 2 treatment variables:
    # The "case" treatment is: (true, false)
    query = (t₁=[true, false], t₂ = [false, true])
    indicators = TMLE.indicator_fns(query)
    @test indicators == Dict(
        (t₁ = 1, t₂ = 0) => 1,
        (t₁ = 0, t₂ = 0) => -1,
        (t₁ = 1, t₂ = 1) => -1,
        (t₁ = 0, t₂ = 1) => 1
    )

    # For 3 treatment variables:
    # The "case" treatment is: (a, c, true)
    query = (t₁=["a", "b"], t₂ = ["c", "d"], t₃ = [true, false])
    indicators = TMLE.indicator_fns(query)
    @test indicators == Dict(
        (t₁ = "b", t₂ = "c", t₃ = 1) => -1, # (-1)^{2}=1
        (t₁ = "a", t₂ = "c", t₃ = 1) => 1,  # (-1)^{0}=1
        (t₁ = "b", t₂ = "d", t₃ = 0) => -1, # (-1)^{3}=-1
        (t₁ = "b", t₂ = "c", t₃ = 0) => 1,  # (-1)^{2}=1
        (t₁ = "a", t₂ = "d", t₃ = 1) => -1, # (-1)^{1}=-1
        (t₁ = "a", t₂ = "c", t₃ = 0) => -1, # (-1)^{1}=-1
        (t₁ = "a", t₂ = "d", t₃ = 0) => 1,  # (-1)^{2}=1
        (t₁ = "b", t₂ = "d", t₃ = 1) => 1   # (-1)^{2}=1
    )
end

@testset "Test adapt" begin
    T = (a=1,)
    @test TMLE.adapt(T) == 1

    T = (a=1, b=2)
    @test TMLE.adapt(T) == T
end

@testset "Test compute_covariate" begin
    # First case: 1 categorical variable
    # Using a trivial classifier
    # that outputs the proportions of of the classes
    T = (t₁ = categorical(["a", "b", "c", "a", "a", "b", "a"]),)
    W = MLJ.table(rand(7, 3))

    Gmach = machine(ConstantClassifier(),
                    W,
                    TMLE.adapt(T))
    fit!(Gmach, verbosity=0)

    indicators = TMLE.indicator_fns((t₁=["a", "b"],))

    cov = TMLE.compute_covariate(Gmach, W, T, indicators)
    @test cov == [1.75,
                 -3.5,
                 0.0,
                 1.75,
                 1.75,
                 -3.5,
                 1.75]

    # Second case: 2 binary variables
    # Using a trivial classifier
    # that outputs the proportions of of the classes
    T = (t₁ = categorical([1, 0, 0, 1, 1, 1, 0]),
         t₂ = categorical([1, 1, 1, 1, 1, 0, 0]))
    W = MLJ.table(rand(7, 3))

    Gmach = machine(FullCategoricalJoint(ConstantClassifier()),
                    W,
                    T)
    fit!(Gmach, verbosity=0)

    indicators = TMLE.indicator_fns((t₁=[1, 0], t₂ = [1, 0]))

    cov = TMLE.compute_covariate(Gmach, W, T, indicators)
    @test cov == [2.3333333333333335,
                 -3.5,
                 -3.5,
                 2.3333333333333335,
                 2.3333333333333335,
                 -7.0,
                 7.0]

    # Third case: 3 mixed categorical variables
    # Using a trivial classifier
    # that outputs the proportions of of the classes
    T = (t₁ = categorical(["a", "a", "b", "b", "c", "b", "b"]),
         t₂ = categorical([3, 2, 1, 1, 2, 2, 2]),
         t₃ = categorical([true, false, true, false, false, false, false]))
    W = MLJ.table(rand(7, 3))

    Gmach = machine(FullCategoricalJoint(ConstantClassifier()),
                    W,
                    T)
    fit!(Gmach, verbosity=0)

    indicators = TMLE.indicator_fns((t₁=["a", "b"], t₂ = [1, 2], t₃ = [true, false]))

    cov = TMLE.compute_covariate(Gmach, W, T, indicators)
    @test cov == [0,
                  7.0,
                 -7,
                  7,
                  0,
                 -3.5,
                 -3.5]
end

@testset "Test compute_offset" begin
    n = 10
    X = rand(n, 3)

    # When Y is binary
    y = categorical([1, 1, 1, 1, 0, 0, 0, 0, 0, 0])
    mach = machine(ConstantClassifier(), MLJ.table(X), y)
    fit!(mach, verbosity=0)
    # Should be equal to logit(Ê[Y|X])= logit(4/10) = -0.4054651081081643
    @test TMLE.compute_offset(mach, X) == repeat([-0.4054651081081643], n)

    # When Y is continuous
    y = [1., 2., 3, 4, 5, 6, 7, 8, 9, 10]
    mach = machine(MLJ.DeterministicConstantRegressor(), MLJ.table(X), y)
    fit!(mach, verbosity=0)
    # Should be equal to Ê[Y|X] = 5.5
    @test TMLE.compute_offset(mach, X) == repeat([5.5], n)

end

@testset "Test logit" begin
    @test TMLE.logit([0.4, 0.8, 0.2]) ≈ [
        -0.40546510810,
        1.38629436112,
        -1.38629436112
    ]
    @test TMLE.logit([1, 0]) == [Inf, -Inf]
end

@testset "Test compute_fluctuation" begin
    n = 10
    W = (w₁=ones(n),)
    T = (t₁=categorical([1, 1, 1, 0, 0, 0, 0, 1, 1, 0]),
         t₂=categorical([0, 0, 1, 0, 0, 1, 1, 1, 0, 0]))
    y = categorical([1, 1, 0, 0, 1, 0 , 1, 0, 0, 0])
    query = (t₁=[1, 0], t₂ = [1, 0])
    indicators = TMLE.indicator_fns(query)

    # Fit encoder
    Hmach = machine(OneHotEncoder(features=[:t₁, :t₂], drop_last=true), T)
    fit!(Hmach, verbosity=0)
    Thot = transform(Hmach)
    # Fit Q̅
    X = merge(Thot, W)
    Q̅mach = machine(ConstantClassifier(), X, y)
    fit!(Q̅mach, verbosity=0)
    # Fit G
    Gmach = machine(FullCategoricalJoint(ConstantClassifier()), W, T)
    fit!(Gmach, verbosity=0)
    # Fit Fluctuation
    offset = TMLE.compute_offset(Q̅mach, X)
    covariate = TMLE.compute_covariate(Gmach, W, T, indicators)
    Xfluct = (covariate=covariate, offset=offset)
    Fmach = machine(LinearBinaryClassifier(fit_intercept=false, offsetcol=:offset), Xfluct, y)
    fit!(Fmach, verbosity=0)

    # We are using constant classifiers
    # The offset is equal to: -0.40546510810 all the time
    # The covariate is ≈ [-3.3, -3.3, 5., 3.3, 3.3, -5., -5., 5., -3.3, 3.3]
    # The fluctuation value is equal to Ê[Y|W, T] where Ê is computed via Fmach
    # The coefficient for the fluctuation seems to be -0.222835 here
    expected_mean(cov) = TMLE.expit(-0.40546510 - 0.22283549318*cov)
    # Let's look at the different counterfactual treatments
    # T₁₁: cov=5.
    T₁₁ = (t₁=categorical(ones(n), levels=levels(T[1])), t₂=categorical(ones(n), levels=levels(T[2])))
    fluct = TMLE.compute_fluctuation(Fmach, Q̅mach, Gmach, indicators, W, T₁₁, X)
    @test fluct ≈ repeat([expected_mean(5.)], n) atol=1e-5
    # T₁₀: cov=-3.333333
    T₁₀ = (t₁=categorical(ones(n), levels=[0, 1]), t₂=categorical(zeros(n), levels=[0, 1]))
    fluct = TMLE.compute_fluctuation(Fmach, Q̅mach, Gmach, indicators, W, T₁₀, X)
    @test fluct ≈ repeat([expected_mean(-3.333333)], n) atol=1e-5
    # T₀₁: cov=-5.
    T₀₁ = (t₁=categorical(zeros(n), levels=[0, 1]), t₂=categorical(ones(n), levels=[0, 1]))
    fluct = TMLE.compute_fluctuation(Fmach, Q̅mach, Gmach, indicators, W, T₀₁, X)
    @test fluct ≈ repeat([expected_mean(-5.)], n) atol=1e-5
    # T₀₀: cov=3.333333
    T₀₀ = (t₁=categorical(zeros(n), levels=[0, 1]), t₂=categorical(zeros(n), levels=[0, 1]))
    fluct = TMLE.compute_fluctuation(Fmach, Q̅mach, Gmach, indicators, W, T₀₀, X)
    @test fluct ≈ repeat([expected_mean(3.333333)], n) atol=1e-5

end

@testset "Test log_over_threshold" begin
    covariate = source([4, 2, 3])
    @test TMLE.log_over_threshold(covariate, 0.4)() == [1, 3]

    # End to end in the fit process
    n = 10000
    rng = StableRNG(123)
    T = (t=categorical(rand(rng, Bernoulli(0.001), n)),)
    W = MLJ.table(rand(rng, n, 3))
    y = rand(rng, n)

    query = (t=[true, false],)

    Q̅ = LinearRegressor()
    G = LogisticClassifier()
    tmle = TMLEstimator(Q̅, G, query)

    mach = machine(tmle, T, W, y)
    fit!(mach, verbosity=0)
    @test length(report(mach).extreme_propensity_idx) == 12
end

end;

true

max_stars_repo_name: olivierlabayle/TMLE

max_issues_repo_name: olivierlabayle/TMLE

max_forks_repo_name: olivierlabayle/TMLE


text:
const ComponentsByType = Dict{DataType, Dict{String, <:InfrastructureSystemsComponent}}

struct Components
    data::ComponentsByType
    time_series_storage::TimeSeriesStorage
    validation_descriptors::Vector
end

function Components(
    time_series_storage::TimeSeriesStorage,
    validation_descriptors = nothing,
)
    if isnothing(validation_descriptors)
        validation_descriptors = Vector()
    end

    return Components(ComponentsByType(), time_series_storage, validation_descriptors)
end

function serialize(components::Components)
    # time_series_storage and validation_descriptors are serialized elsewhere.
    return [serialize(x) for y in values(components.data) for x in values(y)]
end

"""
Add a component.

Throws ArgumentError if the component's name is already stored for its concrete type.

Throws InvalidRange if any of the component's field values are outside of defined valid
range.
"""
function add_component!(
    components::Components,
    component::T;
    skip_validation = false,
    deserialization_in_progress = false,
) where {T <: InfrastructureSystemsComponent}
    component_name = get_name(component)
    if !isconcretetype(T)
        throw(ArgumentError("add_component! only accepts concrete types"))
    end

    if !haskey(components.data, T)
        components.data[T] = Dict{String, T}()
    elseif haskey(components.data[T], component_name)
        throw(ArgumentError("$(component_name) is already stored for type $T"))
    end

    !skip_validation && check_component(components, component)

    if !deserialization_in_progress && has_time_series(component)
        throw(ArgumentError("cannot add a component with time_series: $component"))
    end

    set_time_series_storage!(component, components.time_series_storage)
    components.data[T][component_name] = component
    return
end

function check_components(components::Components)
    for component in iterate_components(components)
        check_component(component)
    end
end

function check_component(components::Components, comp::InfrastructureSystemsComponent)
    if !isempty(components.validation_descriptors) && !validate_fields(components, comp)
        throw(InvalidRange("Invalid value"))
    end

    if !validate_struct(comp)
        throw(InvalidValue("Invalid value for $(comp)"))
    end
end

"""
Removes all components from the system.
"""
function clear_components!(components::Components)
    for type_ in collect(keys(components.data))
        remove_components!(type_, components)
    end
end

"""
Remove all components of type T.

Throws ArgumentError if the type is not stored.
"""
function remove_components!(
    ::Type{T},
    components::Components,
) where {T <: InfrastructureSystemsComponent}
    if !haskey(components.data, T)
        throw(ArgumentError("component $T is not stored"))
    end

    _components = pop!(components.data, T)
    for component in values(_components)
        prepare_for_removal!(component)
    end

    @debug "Removed all components of type" T
    return values(_components)
end

"""
Remove a component by its value.

Throws ArgumentError if the component is not stored.
"""
function remove_component!(
    components::Components,
    component::T,
) where {T <: InfrastructureSystemsComponent}
    return _remove_component!(T, components, get_name(component))
end

"""
Remove a component by its name.

Throws ArgumentError if the component is not stored.
"""
function remove_component!(
    ::Type{T},
    components::Components,
    name::AbstractString,
) where {T <: InfrastructureSystemsComponent}
    return _remove_component!(T, components, name)
end

function _remove_component!(
    ::Type{T},
    components::Components,
    name::AbstractString,
) where {T <: InfrastructureSystemsComponent}
    if !haskey(components.data, T)
        throw(ArgumentError("component $T is not stored"))
    end

    if !haskey(components.data[T], name)
        throw(ArgumentError("component $T name=$name is not stored"))
    end

    component = pop!(components.data[T], name)
    if isempty(components.data[T])
        pop!(components.data, T)
    end
    prepare_for_removal!(component)
    @debug "Removed component" T name
    return component
end

"""
Check to see if a component exists.
"""
function has_component(
    ::Type{T},
    components::Components,
    name::AbstractString,
) where {T <: InfrastructureSystemsComponent}
    !isconcretetype(T) && return !isempty(get_components_by_name(T, components, name))
    !haskey(components.data, T) && return false
    return haskey(components.data[T], name)
end

"""
Get the component of type T with name. Returns nothing if no component matches. If T is an abstract
type then the names of components across all subtypes of T must be unique.

See [`get_components_by_name`](@ref) for abstract types with non-unique names across subtypes.

Throws ArgumentError if T is not a concrete type and there is more than one component with
    requested name
"""
function get_component(
    ::Type{T},
    components::Components,
    name::AbstractString,
)::Union{T, Nothing} where {T <: InfrastructureSystemsComponent}
    if !isconcretetype(T)
        components = get_components_by_name(T, components, name)
        if length(components) > 1
            throw(
                ArgumentError(
                    "More than one abstract component of type $T with name $name in the system. Operation can't continue",
                ),
            )
        end
        return isempty(components) ? nothing : first(components)
    end

    if !haskey(components.data, T)
        @debug "components of type $T are not stored"
        return nothing
    end

    return get(components.data[T], name, nothing)
end

"""
Get the components of abstract type T with name. Note that
InfrastructureSystems enforces unique names on each concrete type but not
across concrete types.

See [`get_component`](@ref) if the concrete type is known.

Throws ArgumentError if T is not an abstract type.
"""
function get_components_by_name(
    ::Type{T},
    components::Components,
    name::AbstractString,
)::Vector{T} where {T <: InfrastructureSystemsComponent}
    if isconcretetype(T)
        throw(ArgumentError("get_components_by_name does not support concrete types: $T"))
    end

    _components = Vector{T}()
    for key in keys(components.data)
        if key <: T
            component = get_component(key, components, name)
            if !isnothing(component)
                push!(_components, component)
            end
        end
    end

    return _components
end

"""
Returns an iterator of components. T can be concrete or abstract.
Call collect on the result if an array is desired.

# Arguments
- `T`: component type
- `components::Components`: Components of the sytem
- `filter_func::Union{Nothing, Function} = nothing`: Optional function that accepts a component
   of type T and returns a Bool. Apply this function to each component and only return components
   where the result is true.

See also: [`iterate_components`](@ref)
"""
function get_components(
    ::Type{T},
    components::Components,
    filter_func::Union{Nothing, Function} = nothing,
) where {T <: InfrastructureSystemsComponent}
    if isconcretetype(T)
        _components = get(components.data, T, nothing)
        if !isnothing(filter_func) && !isnothing(_components)
            _filter_func = x -> filter_func(x.second)
            _components = values(filter(_filter_func, _components))
        end
        if isnothing(_components)
            iter = FlattenIteratorWrapper(T, Vector{Base.ValueIterator}([]))
        else
            iter =
                FlattenIteratorWrapper(T, Vector{Base.ValueIterator}([values(_components)]))
        end
    else
        types = [x for x in keys(components.data) if x <: T]
        if isnothing(filter_func)
            _components = [values(components.data[x]) for x in types]
        else
            _filter_func = x -> filter_func(x.second)
            _components = [values(filter(_filter_func, components.data[x])) for x in types]
        end
        iter = FlattenIteratorWrapper(T, _components)
    end

    @assert_op eltype(iter) == T
    return iter
end

"""
Iterates over all components.

# Examples
```Julia
for component in iterate_components(obj)
    @show component
end
```

See also: [`get_components`](@ref)
"""
function iterate_components(components::Components)
    Channel() do channel
        for comp_dict in values(components.data)
            for component in values(comp_dict)
                put!(channel, component)
            end
        end
    end
end

function iterate_components_with_time_series(components::Components)
    Channel() do channel
        for comp_dict in values(components.data)
            for component in values(comp_dict)
                if has_time_series(component)
                    put!(channel, component)
                end
            end
        end
    end
end

function get_num_components(components::Components)
    count = 0
    for components in values(components.data)
        count += length(components)
    end
    return count
end

function clear_time_series!(components::Components)
    for component in iterate_components_with_time_series(components)
        clear_time_series!(component)
    end
end


text:
# Calculate intersections with surface(s)

@inline point(r::Ray{T}, t::T) where T <: AbstractFloat = r.origin .+ t .* r.dir # equivalent to C++'s ray.at()

"""Equivalent to `hit_record.set_face_normal()`"""
@inline @fastmath function ray_to_HitRecord(t::T, p, outward_n⃗, r_dir::Vec3{T}, mat::Material{T})::Union{HitRecord,Nothing} where T
	front_face = r_dir ⋅ outward_n⃗ < 0
	n⃗ = front_face ? outward_n⃗ : -outward_n⃗
	HitRecord(t,p,n⃗,front_face,mat)
end

@inline @fastmath function hit(s::Sphere{T}, r::Ray{T}, tmin::T, tmax::T)::Union{HitRecord,Nothing} where T
    oc = r.origin - s.center
    #a = r.dir ⋅ r.dir # unnecessary since `r.dir` is normalized
	a = 1
	half_b = oc ⋅ r.dir
    c = oc⋅oc - s.radius^2
    discriminant = half_b^2 - a*c
	if discriminant < 0 return nothing end # no hit!
	sqrtd = √discriminant

	# Find the nearest root that lies in the acceptable range
	root = (-half_b - sqrtd) / a
	if root < tmin || tmax < root
		root = (-half_b + sqrtd) / a
		if root < tmin || tmax < root
			return nothing # no hit!
		end
	end

	t = root
	p = point(r, t)
	n⃗ = (p - s.center) / s.radius
	return ray_to_HitRecord(t, p, n⃗, r.dir, s.mat)
end

"""Find closest hit between `Ray r` and a list of Hittable objects `h`, within distance `tmin` < `tmax`"""
@inline function hit(hittables::HittableList, r::Ray{T}, tmin::T, tmax::T)::Union{HitRecord,Nothing} where T
    closest = tmax # closest t so far
    best_rec::Union{HitRecord,Nothing} = nothing # by default, no hit
    @inbounds for i in eachindex(hittables) # @paulmelis reported gave him a 4X speedup?!
		h = hittables[i]
        rec = hit(h, r, tmin, closest)
        if rec !== nothing
            best_rec = rec
            closest = best_rec.t # i.e. ignore any further hit > this one's.
        end
    end
    best_rec
end


text:
using CayleyDickson
using Base.Test: @test, @test_throws

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigInt()
    l = Exo1Quaternion(Exo1Real(a, b))
    r = Exo1Quaternion(a, 0, 0, 0, b, 0, 0, 0)
    l == r
end

@test begin
    a = CayleyDickson.randomBigFloat()
    isreal(Exo1Quaternion(a))
end

@test begin
    a = 1
    !isreal(Exo1Quaternion(a, a, a, a))
end

@test begin
    z = random(Exo1Quaternion{Int})
    z == +(z)
end

@test begin
    a = CayleyDickson.randomBigFloat()
    real(Exo1Quaternion(a)) == a
end

@test begin
    io = IOBuffer()
    show(io, Exo1Quaternion(1.0, 2, 3, 4, 5, 6, 7, 8))
    l = String(take!(io))
    r = "[1: 1.0, i: 2.0, j: 3.0, ij: 4.0, a: 5.0, ia: 6.0, ja: 7.0, (ij)a: 8.0]"
    l == r
end

@test begin
    io = IOBuffer()
    show(io, Exo1Quaternion(1.0, 2, 3, 4, 5, 6, 7))
    l = String(take!(io))
    r = "[1: 1.0, i: 2.0, j: 3.0, ij: 4.0, a: 5.0, ia: 6.0, ja: 7.0, (ij)a: 0.0]"
    l == r
end

@test begin
    io = IOBuffer()
    show(io, Exo1Quaternion(1.0, 2, 3, 4, 5, 6))
    l = String(take!(io))
    r = "[1: 1.0, i: 2.0, j: 3.0, ij: 4.0, a: 5.0, ia: 6.0, ja: 0.0, (ij)a: 0.0]"
    l == r
end

@test begin
    io = IOBuffer()
    show(io, Exo1Quaternion(1.0, 2, 3, 4, 5))
    l = String(take!(io))
    r = "[1: 1.0, i: 2.0, j: 3.0, ij: 4.0, a: 5.0, ia: 0.0, ja: 0.0, (ij)a: 0.0]"
    l == r
end

@test begin
    io = IOBuffer()
    show(io, Exo1Quaternion(1.0, 2, 3, 4))
    l = String(take!(io))
    r = "[1: 1.0, i: 2.0, j: 3.0, ij: 4.0, a: 0.0, ia: 0.0, ja: 0.0, (ij)a: 0.0]"
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigFloat()
    c = CayleyDickson.randomBigInt()
    d = CayleyDickson.randomBigFloat()
    l = Exo1Quaternion(Quaternion(a, b, c, d))
    r = Exo1Quaternion(a, b, c, d)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigFloat()
    l = Exo1Quaternion(Binion(a, b))
    r = Exo1Quaternion(a, b)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigFloat()
    c = CayleyDickson.randomBigInt()
    l = Exo1Quaternion(a, b, c)
    r = Exo1Quaternion(a, b, c, 0)
    l == r
end

@test begin
    x = random(Exo1Quaternion{BigFloat})
    y = random(Exo1Quaternion{BigFloat})
    x != y
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(Exo1Quaternion{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = random(Quaternion{BigInt})
    z = random(Exo1Quaternion{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(Exo1Quaternion{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    a = random(Quaternion{BigInt})
    z = random(Exo1Quaternion{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    x = random(Exo1Quaternion{BigInt})
    y = random(Exo1Quaternion{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(Exo1Quaternion{BigInt})
    y = random(Exo1Quaternion{BigInt})
    l = x * y
    r = y * x
    l != r
end

@test begin
    x = random(Exo1Quaternion{BigInt})
    y = random(Exo1Quaternion{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(Exo1Quaternion{BigInt})
    y = random(Exo1Quaternion{BigInt})
    z = random(Exo1Quaternion{BigInt})
    l = (x + y) + z
    r = x + (y + z)
    l == r
end

@test begin
    x = random(Exo1Quaternion{BigInt})
    y = random(Exo1Quaternion{BigInt})
    z = random(Exo1Quaternion{BigInt})
    l = (x * y) * z
    r = x * (y * z)
    l != r
end

@test begin
    x = random(Exo1Quaternion{BigInt})
    y = random(Exo1Quaternion{BigInt})
    l = alternatorL(x, y)
    r = zero(Exo1Quaternion{BigInt})
    l == r
end

@test begin
    x = random(Exo1Quaternion{BigInt})
    y = random(Exo1Quaternion{BigInt})
    l = alternatorR(x, y)
    r = zero(Exo1Quaternion{BigInt})
    l == r
end

@test begin
    x = random(Exo1Quaternion{BigInt})
    conj(conj(x)) == x
end

@test begin
    x = random(Exo1Quaternion{Rational{BigInt}})
    inv(inv(x)) == x
end

@test begin
    x = random(Exo1Quaternion{BigInt})
    y = random(Exo1Quaternion{BigInt})
    l = conj(x * y)
    r = conj(y) * conj(x)
    l == r
end

@test begin
    x = random(Exo1Quaternion{BigInt})
    y = random(Exo1Quaternion{BigInt})
    l = conj(x + y)
    r = conj(x) + conj(y)
    l == r
end

@test begin
    x = random(Exo1Quaternion{BigInt})
    y = random(Exo1Quaternion{BigInt})
    l = conj(x - y)
    r = conj(x) - conj(y)
    l == r
end

@test begin
    x = random(Exo1Quaternion{Rational{BigInt}})
    y = random(Exo1Quaternion{Rational{BigInt}})
    l = inv(x * y)
    r = inv(y) * inv(x)
    l == r
end

@test begin
    x = random(Exo1Quaternion{Rational{BigInt}})
    y = random(Exo1Quaternion{Rational{BigInt}})
    l = inv(x / y)
    r = y / x
    l == r
end

@test begin
    x = random(Exo1Quaternion{Rational{BigInt}})
    y = random(Exo1Quaternion{Rational{BigInt}})
    l = inv(x \ y)
    r = y \ x
    l == r
end

@test begin
    x = random(Exo1Quaternion{BigInt})
    y = random(Exo1Quaternion{BigInt})
    z = random(Exo1Quaternion{BigInt})
    l = x * (y + z)
    r = (x * y) + (x * z)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(Exo1Quaternion{BigInt})
    y = random(Exo1Quaternion{BigInt})
    l = a * (x + y)
    r = (a * x) + (a * y)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(Exo1Quaternion{BigInt})
    l = a * x
    r = x * a
    l == r
end

@test begin
    a = random(Quaternion{BigInt})
    x = random(Exo1Quaternion{BigInt})
    l = commutator(a, x)
    r = zero(Exo1Quaternion{BigInt})
    l != r
end

@test begin
    x = random(Exo1Quaternion{BigInt})
    y = random(Exo1Quaternion{BigInt})
    z = random(Exo1Quaternion{BigInt})
    l = (x + y) * z
    r = (x * z) + (y * z)
    l == r
end

@test begin
    x = random(Exo1Quaternion{BigInt})
    y = random(Exo1Quaternion{BigInt})
    z = random(Exo1Quaternion{BigInt})
    l = x * (y - z)
    r = (x * y) - (x * z)
    l == r
end

@test begin
    x = random(Exo1Quaternion{BigInt})
    y = random(Exo1Quaternion{BigInt})
    z = random(Exo1Quaternion{BigInt})
    l = (x - y) * z
    r = (x * z) - (y * z)
    l == r
end

@test begin
    x = random(Exo1Quaternion{BigInt})
    quadrance(x) > big(0)
end

@test begin
    x = random(Exo1Quaternion{BigInt})
    y = random(Exo1Quaternion{BigInt})
    l = quadrance(x * y)
    r = quadrance(x) * quadrance(y)
    l == r
end

@test_throws ErrorException begin
    inv(Exo1Quaternion(0, 0.0))
end

@test_throws ErrorException begin
    random(Exo1Quaternion{Int}) / 0
end

@test_throws ErrorException begin
    0 \ random(Exo1Quaternion{Int})
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(Exo1Quaternion{Rational{BigInt}})
    l = z / a
    r = a \ z
    l == r
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(Exo1Quaternion{Rational{BigInt}})
    l = a / z
    r = z \ a
    l == r
end

@test begin
    z = random(Exo1Quaternion{BigInt})
    z == z + zero(Exo1Quaternion{BigInt})
end

@test begin
    z = random(Exo1Quaternion{BigInt})
    z == z * one(Exo1Quaternion{BigInt})
end

@test begin
    z = random(Exo1Quaternion{BigInt})
    z == z + zero(z)
end

@test begin
    z = random(Exo1Quaternion{BigInt})
    z == z * one(z)
end


text:
################################################################################

"Feed on hidden Markov model object."
function feed!(self::HMM, frame::Int64, d::Array{Float64, 2}; params::HMMParams)
  for ι ∈ eachindex(self.model)
    plus = euclDist(self.data[ι], d[frame, :])
    for ε ∈ eachindex(self.model)
      lpen = copy(params.penalty)
      if (ε == ι) lpen = 0 end
      if self.model[ε][frame + 1] < 0. || self.model[ι][frame] + plus + lpen < self.model[ε][frame + 1]
        self.model[ε][frame + 1] = self.model[ι][frame] + plus + lpen
      end
    end
  end
end

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


text:
using Compat.Test
using DataStructures
const IntSet = DataStructures.IntSet
using Primes, Compat
using Compat.Unicode: lowercase

@test isempty(detect_ambiguities(Base, Core, DataStructures))

tests = ["int_set",
         "deque",
         "circ_deque",
         "sorted_containers",
         "stack_and_queue",
         "accumulator",
         "classified_collections",
         "disjoint_set",
         "binheap",
         "mutable_binheap",
         "default_dict",
         "ordered_dict",
         "ordered_set",
         "trie",
         "list",
         "multi_dict",
         "circular_buffer",
         "sorting",
         "priority_queue"
        ]

if length(ARGS) > 0
    tests = ARGS
end

@testset "DataStructures" begin

for t in tests
    fp = joinpath(dirname(@__FILE__), "test_$t.jl")
    println("$fp ...")
    include(fp)
end

end # @testset


text:
include("module.jl")

import .Example

Example.fn()

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox


text:
module ProjAssistant

    import DrWatson
    import Logging
    import Logging: SimpleLogger, global_logger, with_logger
    import Serialization: serialize, deserialize
    import FileIO
    import Requires: @require
    import Pkg
    using FilesTreeTools
    using ExtractMacro
    using Base.Threads
    using DataFileNames

    include("utils.jl")
    include("_io_print.jl")
    include("io_global_conf.jl")
    include("_save_load.jl")
    include("scache_lcache.jl")
    include("sdat_ldat.jl")
    include("group_files.jl")
    include("gen_sub_proj.jl")
    include("gen_top_proj.jl")
    include("create_proj_dirs.jl")
    include("sfig_sgif.jl")
    include("sglob_lglob_macros_ex.jl")
    include("sglob_lglob.jl")
    include("proj_functions.jl")
    include("fileid.jl")
    include("quickactivate.jl")

    export parentproj, topproj, istop_proj
    export projdir, projname, set_projdir
    export devdir, datdir, srcdir, plotsdir, scriptsdir, papersdir
    export cachedir, procdir, rawdir
    export lsdevdir, lsdatdir, lssrcdir, lsplotsdir, lsscriptsdir,
           lspapersdir, lsprocdir, lsrawdir, lscachedir
    export sdat, sprocdat, srawdat, sfig, sgif
    export ldat, lprocdat, lrawdat, lfig
    export cfname, is_cfname, scache, lcache, delcache
    export globfile, sglob, lglob, delglob, @lglob, @sglob
    export gen_top_proj, @gen_top_proj
    export gen_sub_proj, @gen_sub_proj
    export create_proj_dirs, @create_proj_dirs
    export @fileid
    export quickactivate, @quickactivate

    # Re-exports
    export dfname, parse_dfname
    export @extract
    export walkdown, filtertree

    function __init__()
        _init_globals()

        @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
            import ImgTools
        end
    end

end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

module TestIO

using MXNet
using Base.Test

using ..Main: rand_dims

function test_mnist()
  info("IO::MNIST")
  filenames = mx.get_mnist_ubyte()

  batch_size = 10
  mnist_provider = mx.MNISTProvider(image=filenames[:train_data],
                                    label=filenames[:train_label],
                                    batch_size=batch_size, silent=true, shuffle=false)
  data_spec = mx.provide_data(mnist_provider)
  label_spec = mx.provide_label(mnist_provider)
  @test data_spec == [(:data, (28,28,1,batch_size))]
  @test label_spec == [(:softmax_label, (batch_size,))]

  n_batch = 0
  for batch in mnist_provider
    if n_batch == 0
      data_array  = mx.empty(28,28,1,batch_size)
      label_array = mx.empty(batch_size)
      # have to use "for i=1:1" to get over the legacy "feature" of using
      # [ ] to do concatenation in Julia
      data_targets = [[(1:batch_size, data_array)] for i = 1:1]
      label_targets = [[(1:batch_size, label_array)] for i = 1:1]

      mx.load_data!(mnist_provider, batch, data_targets)
      mx.load_label!(mnist_provider, batch, label_targets)

      true_labels = [5,0,4,1,9,2,1,3,1,4] # the first 10 labels in MNIST train
      got_labels  = Int[copy(label_array)...]
      @test true_labels == got_labels
    end

    n_batch += 1
  end

  @test n_batch == 60000 / batch_size
end

function test_arrays_impl(data::Vector, label::Vector, provider::mx.ArrayDataProvider)
  data = convert(Vector{Array{Float64}}, data)
  label = convert(Vector{Array{Float64}}, label)

  sample_count = size(data[1])[end]
  batch_size   = mx.get_batch_size(provider)
  idx_all = 1:batch_size:sample_count

  for (d1, (_, d2)) in zip(data, mx.provide_data(provider))
    @test size(d1)[1:end-1] == d2[1:end-1]
    @test batch_size == d2[end]
  end
  for (d1, (_, d2)) in zip(label, mx.provide_label(provider))
    @test size(d1)[1:end-1] == d2[1:end-1]
    @test batch_size == d2[end]
  end

  info("IO::Array::#data=$(length(data)),#label=$(length(label)),batch_size=$batch_size")
  for (idx, batch) in zip(idx_all, provider)
    data_batch = [x[[Colon() for i=1:ndims(x)-1]..., idx:min(idx+batch_size-1,sample_count)] for x in data]
    data_get   = mx.get_data(provider, batch)

    for (d_real, d_get) in zip(data_batch, data_get)
      @test d_real ≈ copy(d_get)[[1:n for n in size(d_real)]...]
      @test mx.count_samples(provider, batch) == size(d_real)[end]
    end
  end
end

function test_arrays()
  sample_count = 15
  batch_size   = 4
  dims_data    = [rand_dims()..., sample_count]
  data         = rand(dims_data...)
  provider     = mx.ArrayDataProvider(data, batch_size=batch_size)
  test_arrays_impl(Array[data], [], provider)

  dims_label   = [rand_dims()..., sample_count]
  label        = rand(dims_label...)
  provider     = mx.ArrayDataProvider(data, label, batch_size=batch_size)
  test_arrays_impl(Array[data], Array[label], provider)

  provider     = mx.ArrayDataProvider(:data=>data, :my_label=>label, batch_size=batch_size)
  test_arrays_impl(Array[data], Array[label], provider)

  dims_data2   = [rand_dims()..., sample_count]
  data2        = rand(dims_data2...)
  provider     = mx.ArrayDataProvider((:data=>data, :data2=>data2), label, batch_size=batch_size)
  test_arrays_impl(Array[data,data2], Array[label], provider)
end

function test_arrays_shuffle()
  info("IO::Array::shuffle")

  sample_count = 15
  batch_size   = 4
  data         = rand(mx.MX_float, 1, sample_count)
  label        = collect(1:sample_count)
  provider     = mx.ArrayDataProvider(data, :index => label, batch_size=batch_size, shuffle=true)

  idx_all      = 1:batch_size:sample_count
  data_got     = similar(data)
  label_got    = similar(label)
  for (idx, batch) in zip(idx_all, provider)
    data_batch  = mx.get(provider, batch, :data)
    label_batch = mx.get(provider, batch, :index)
    ns_batch    = mx.count_samples(provider, batch)
    data_got[idx:idx+ns_batch-1]  = copy(data_batch)[1:ns_batch]
    label_got[idx:idx+ns_batch-1] = copy(label_batch)[1:ns_batch]
  end

  @test label_got != label
  @test sort(label_got) == label
  @test size(data_got) == size(data[:, Int[label_got...]])
  @test data_got ≈ data[:, Int[label_got...]]
end

@testset "IO Test" begin
  test_arrays_shuffle()
  test_arrays()
  test_mnist()
end

end

max_stars_repo_name: Vikas-kum/incubator-mxnet

max_issues_repo_name: Vikas-kum/incubator-mxnet

max_forks_repo_name: Vikas-kum/incubator-mxnet


text:
using Documenter, DEMC

makedocs(;
    modules=[DEMC],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/chrished/DEMC.jl/blob/{commit}{path}#L{line}",
    sitename="DEMC.jl",
    authors="christoph hedtrich"
)

deploydocs(;
    repo="github.com/chrished/DEMC.jl",
)


text:
using CmdStan, StanMCMCChain, Test, Statistics

ProjDir = joinpath(dirname(@__FILE__), "..", "examples", "Bernoulli")
cd(ProjDir) do

  isdir("tmp") && rm("tmp", recursive=true);

  include("../examples/Bernoulli/bernoulli.jl")

  @test 0.1 <  mean(chains.value[:, 8, :] ) < 0.6

  isdir("tmp") && rm("tmp", recursive=true);

end # cd


text:
using Test
using MPI
using Ipopt
using ExaPF
using JuMP
using ProxAL
using DelimitedFiles, Printf
using LazyArtifacts

const DATA_DIR = joinpath(artifact"ExaData", "ExaData")
case = "case9"
T = 3
ramp_scale = 0.5
load_scale = 1.0
maxρ = 0.1
quad_penalty = 0.1
rtol = 1e-4

# Load data
case_file = joinpath(DATA_DIR, "$(case).m")
load_file = joinpath(DATA_DIR, "mp_demand", "$(case)_oneweek_168")

@testset "Block Model Backends" begin
    # ctgs_arr = deepcopy(rawdata.ctgs_arr)

    # Model/formulation settings
    modelinfo = ModelInfo()
    modelinfo.num_time_periods = T
    modelinfo.load_scale = load_scale
    modelinfo.ramp_scale = ramp_scale
    modelinfo.allow_obj_gencost = true
    modelinfo.allow_constr_infeas = false

    # Algorithm settings
    algparams = AlgParams()
    algparams.verbose = 0
    algparams.mode = :coldstart
    algparams.ρ_t = maxρ
    algparams.ρ_c = maxρ
    algparams.τ = 3.0*maxρ

    modelinfo.case_name = case
    algparams.optimizer = optimizer_with_attributes(
        Ipopt.Optimizer,
        "print_level" => Int64(algparams.verbose > 0)*5,
    )

    K = 0

    # Instantiate primal and dual buffers
    nlp = ProxALEvaluator(case_file, load_file, modelinfo, algparams, JuMPBackend(), nothing)
    primal = ProxAL.OPFPrimalSolution(nlp)
    dual = ProxAL.OPFDualSolution(nlp)

    modelinfo_local = deepcopy(nlp.modelinfo)
    modelinfo_local.num_time_periods = 1
    blkid = 1

    @testset "Timestep $t" for t in 1:T
        opfdata_c = ProxAL.opf_loaddata(nlp.rawdata;
                            time_horizon_start = t,
                            time_horizon_end = t,
                            load_scale = load_scale,
                            ramp_scale = ramp_scale,
                            corr_scale = modelinfo.corr_scale)

        local solution, n
        @testset "JuMP Block backend" begin
            blockmodel = ProxAL.JuMPBlockBackend(blkid, opfdata_c, nlp.rawdata, algparams, modelinfo_local, t, 1, T)
            ProxAL.init!(blockmodel, nlp.algparams)

            ProxAL.set_objective!(blockmodel, nlp.algparams, primal, dual)
            n = JuMP.num_variables(blockmodel.model)
            x0 = zeros(n)
            solution = ProxAL.optimize!(blockmodel, x0, nlp.algparams)
            @test solution.status ∈ ProxAL.MOI_OPTIMAL_STATUSES
        end
        obj_jump = solution.minimum
        pg_jump = solution.pg
        slack_jump = solution.st

        # @testset "ExaPF Block backend" begin
        #     # TODO: currently, we need to build directly ExaPF object
        #     # with rawdata, as ExaPF is dealing only with struct of arrays objects.
        #     blockmodel = ProxAL.ExaBlockBackend(blkid, opfdata_c, nlp.rawdata, algparams, modelinfo_local, t, 1, T)
        #     ProxAL.init!(blockmodel, nlp.algparams)
        #     ProxAL.set_objective!(blockmodel, nlp.algparams, primal, dual)

        #     # Better to set nothing than 0 to avoid convergence issue
        #     # in the powerflow solver.
        #     x0 = nothing

        #     # Set up optimizer
        #     algparams.gpu_optimizer = optimizer_with_attributes(
        #         Ipopt.Optimizer,
        #         "print_level" => 0,
        #         "limited_memory_max_history" => 50,
        #         "hessian_approximation" => "limited-memory",
        #         "derivative_test" => "first-order",
        #         "tol" => 1e-6,
        #     )

        #     solution = ProxAL.optimize!(blockmodel, x0, algparams)
        #     @test solution.status ∈ ProxAL.MOI_OPTIMAL_STATUSES
        #     obj_exa = solution.minimum
        #     pg_exa = solution.pg
        #     slack_exa = solution.st
        #     @test obj_jump ≈ obj_exa
        #     @test pg_jump ≈ pg_exa rtol=1e-6
        #     if t > 1  # slack could be of any value for t == 1
        #         @test slack_jump ≈ slack_exa rtol=1e-5
        #     end
        # end

        @testset "ExaTron BlockModel" begin
            blockmodel = ProxAL.TronBlockBackend(
                blkid, opfdata_c, nlp.rawdata, algparams, modelinfo_local, t, 1, T;
            )
            ProxAL.init!(blockmodel, nlp.algparams)
            ProxAL.set_objective!(blockmodel, nlp.algparams, primal, dual)

            # Test optimization
            x0 = nothing
            solution = ProxAL.optimize!(blockmodel, x0, nlp.algparams)
            @test solution.status ∈ ProxAL.MOI_OPTIMAL_STATUSES
            obj_tron   = solution.minimum
            pg_tron    = solution.pg
            slack_tron = solution.st
            # TODO: implement ProxAL objective in ExaTron
            @test_broken obj_jump ≈ obj_tron
            @test pg_jump ≈ pg_tron rtol=1e-2
            if t > 1  # slack could be of any value for t == 1
                @test slack_jump ≈ slack_tron rtol=1e-2
            end
        end
    end

    @testset "OPFBlocks" begin
        blocks = ProxAL.OPFBlocks(
            nlp.opfdata, nlp.rawdata;
            modelinfo=nlp.modelinfo,
            algparams=nlp.algparams,
            comm=nothing
        )

        @test length(blocks.blkModel) == nlp.modelinfo.num_time_periods
    end
end

max_stars_repo_name: exanauts/proxALM

max_issues_repo_name: exanauts/proxALM

max_forks_repo_name: exanauts/proxALM


text:
#=
        ScalarAffineFunction -in- LessThan
        ScalarAffineFunction -in- GreaterThan
        ScalarAffineFunction -in- EqualTo
        ScalarAffineFunction -in- Interval
=#
constrdict(m::LinQuadOptimizer, ::LCI{LE})  = cmap(m).less_than
constrdict(m::LinQuadOptimizer, ::LCI{GE})  = cmap(m).greater_than
constrdict(m::LinQuadOptimizer, ::LCI{EQ})  = cmap(m).equal_to
constrdict(m::LinQuadOptimizer, ::LCI{IV})  = cmap(m).interval

function MOI.addconstraint!(m::LinQuadOptimizer, func::Linear, set::T) where T <: LinSets
    cfunc = MOIU.canonical(func)
    addlinearconstraint!(m, cfunc, set)
    m.last_constraint_reference += 1
    ref = LCI{T}(m.last_constraint_reference)
    dict = constrdict(m, ref)
    dict[ref] = get_number_linear_constraints(m)
    push!(m.constraint_primal_solution, NaN)
    push!(m.constraint_dual_solution, NaN)
    push!(m.constraint_constant, func.constant)
    return ref
end

function addlinearconstraint!(m::LinQuadOptimizer, func::Linear, set::S) where S <: Union{LE, GE, EQ}
    addlinearconstraint!(m, func, backend_type(m,set), _getrhs(set))
end

function addlinearconstraint!(m::LinQuadOptimizer, func::Linear, set::IV)
    add_ranged_constraints!(m, [1], getcol.(m, func.variables), func.coefficients, [set.lower], [set.upper])
end

function addlinearconstraint!(m::LinQuadOptimizer, func::Linear, sense::Cchar, rhs)
    if abs(func.constant) > eps(Float64)
        warn("Constant in scalar function moved into set.")
    end
    add_linear_constraints!(m, [1], getcol.(m, func.variables), func.coefficients, [sense], [rhs - func.constant])
end

#=
    Add linear constraints (plural)
=#

function MOI.addconstraints!(m::LinQuadOptimizer, func::Vector{Linear}, set::Vector{S}) where S <: LinSets
    # canonicalize
    cfunc = MOIU.canonical.(func)

    @assert length(cfunc) == length(set)
    numrows = get_number_linear_constraints(m)
    addlinearconstraints!(m, cfunc, set)
    crefs = Vector{LCI{S}}(length(cfunc))
    for i in 1:length(cfunc)
        m.last_constraint_reference += 1
        ref = LCI{S}(m.last_constraint_reference)
        dict = constrdict(m, ref)
        dict[ref] = numrows + i
        push!(m.constraint_primal_solution, NaN)
        push!(m.constraint_dual_solution, NaN)
        push!(m.constraint_constant, cfunc[i].constant)
        crefs[i] = ref
    end
    return crefs
end

function addlinearconstraints!(m::LinQuadOptimizer, func::Vector{Linear}, set::Vector{S}) where S <: LinSets
    addlinearconstraints!(m, func, backend_type.(m,set), [_getrhs(s) for s in set])
end

function addlinearconstraints!(m::LinQuadOptimizer, func::Vector{Linear}, set::Vector{IV})
    # loop through once to get number of non-zeros and to move rhs across
    lowerbounds = [s.lower for s in set]
    upperbounds = [s.upper for s in set]
    nnz = 0
    for (i, f) in enumerate(func)
        if abs(f.constant) > eps(Float64)
            warn("Constant in scalar function moved into set.")
            lowerbounds[i] -= f.constant
            upperbounds[i] -= f.constant
        end
        nnz += length(f.coefficients)
    end
    row_starts     = Vector{Int}(length(func))  # index of start of each row
    column_indices = Vector{Int}(nnz)           # flattened columns for each function
    coefficients   = Vector{Float64}(nnz)       # corresponding non-zeros
    i = 1
    for (fi, f) in enumerate(func)
        row_starts[fi] = i
        for (var, coef) in zip(f.variables, f.coefficients)
            column_indices[i] = getcol(m, var)
            coefficients[i]   = coef
            i += 1
        end
    end
    add_ranged_constraints!(m, row_starts, column_indices, coefficients, lowerbounds, upperbounds)
end

function addlinearconstraints!(m::LinQuadOptimizer, func::Vector{Linear}, sense::Vector{Cchar}, rhs::Vector{Float64})
    # loop through once to get number of non-zeros and to move rhs across
    nnz = 0
    for (i, f) in enumerate(func)
        if abs(f.constant) > eps(Float64)
            warn("Constant in scalar function moved into set.")
            rhs[i] -= f.constant
        end
        nnz += length(f.coefficients)
    end
    rowbegins = Vector{Int}(length(func))   # index of start of each row
    column_indices = Vector{Int}(nnz)       # flattened columns for each function
    nnz_vals = Vector{Float64}(nnz)         # corresponding non-zeros
    cnt = 1
    for (fi, f) in enumerate(func)
        rowbegins[fi] = cnt
        for (var, coef) in zip(f.variables, f.coefficients)
            column_indices[cnt] = getcol(m, var)
            nnz_vals[cnt] = coef
            cnt += 1
        end
    end
    add_linear_constraints!(m, rowbegins, column_indices, nnz_vals, sense, rhs)
end

#=
    Constraint set of Linear function
=#

MOI.canget(m::LinQuadOptimizer, ::MOI.ConstraintSet, ::Type{LCI{S}}) where S <: Union{LE, GE, EQ} = true
function MOI.get(m::LinQuadOptimizer, ::MOI.ConstraintSet, c::LCI{S}) where S <: Union{LE, GE, EQ}
    rhs = get_rhs(m, m[c])
    S(rhs+m.constraint_constant[m[c]])
end

MOI.canget(m::LinQuadOptimizer, ::MOI.ConstraintSet, ::Type{LCI{IV}}) = false
# TODO(odow): get constraint sets for ranged constraints.
# function MOI.get(m::LinQuadOptimizer, ::MOI.ConstraintSet, c::LCI{IV})
#     ???
#     IV(lowerbound+m.constraint_constant[m[c]], upperbound + m.constraint_constant[m[c]])
# end

#=
    Constraint function of Linear function
=#

MOI.canget(m::LinQuadOptimizer, ::MOI.ConstraintFunction, ::Type{<:LCI{<: LinSets}}) = true
function MOI.get(m::LinQuadOptimizer, ::MOI.ConstraintFunction, c::LCI{<: LinSets})
    # TODO more efficiently
    colidx, coefs = get_linear_constraint(m, m[c])
    Linear(m.variable_references[colidx+1], coefs, -m.constraint_constant[m[c]])
end

#=
    Scalar Coefficient Change of Linear Constraint
=#

MOI.canmodifyconstraint(m::LinQuadOptimizer, c::LCI{<: LinSets}, ::Type{MOI.ScalarCoefficientChange{Float64}}) = true
function MOI.modifyconstraint!(m::LinQuadOptimizer, c::LCI{<: LinSets}, chg::MOI.ScalarCoefficientChange{Float64})
    col = m.variable_mapping[chg.variable]
    change_coefficient!(m, m[c], col, chg.new_coefficient)
end

#=
    Change RHS of linear constraint without modifying sense
=#

MOI.canmodifyconstraint(m::LinQuadOptimizer, c::LCI{S}, ::Type{S}) where S <: Union{LE, GE, EQ} = true
function MOI.modifyconstraint!(m::LinQuadOptimizer, c::LCI{S}, newset::S) where S <: Union{LE, GE, EQ}
    # the column 0 (or -1 in 0-index) is the rhs.
    change_coefficient!(m, m[c], 0, _getrhs(newset))
end

MOI.canmodifyconstraint(m::LinQuadOptimizer, c::LCI{IV}, ::Type{IV}) = true
function MOI.modifyconstraint!(m::LinQuadOptimizer, c::LCI{IV}, set::IV)
    modify_ranged_constraints!(m, [m[c]], [set.lower], [set.upper])
end

#=
    Delete a linear constraint
=#

MOI.candelete(m::LinQuadOptimizer, c::LCI{<: LinSets}) = MOI.isvalid(m, c)
function MOI.delete!(m::LinQuadOptimizer, c::LCI{<: LinSets})
    deleteconstraintname!(m, c)
    dict = constrdict(m, c)
    row = dict[c]
    delete_linear_constraints!(m, row, row)
    deleteat!(m.constraint_primal_solution, row)
    deleteat!(m.constraint_dual_solution, row)
    deleteat!(m.constraint_constant, row)
    # shift all the other references
    shift_references_after_delete_affine!(m, row)
    delete!(dict, c)
end

#=
    Transform scalar constraint
=#

function MOI.cantransformconstraint(m::LinQuadOptimizer, ref::LCI{S}, newset::S) where S
    false
end
function MOI.transformconstraint!(m::LinQuadOptimizer, ref::LCI{S}, newset::S) where S
    error("Cannot transform constraint of same set. use `modifyconstraint!` instead.")
end

function MOI.cantransformconstraint(m::LinQuadOptimizer, ref::LCI{S1}, newset::S2) where S1 where S2 <: Union{LE, GE, EQ}
    true
end
function MOI.transformconstraint!(m::LinQuadOptimizer, ref::LCI{S1}, newset::S2) where S1 where S2 <: Union{LE, GE, EQ}
    dict = constrdict(m, ref)
    row = dict[ref]
    change_linear_constraint_sense!(m, [row], [backend_type(m,newset)])
    m.last_constraint_reference += 1
    ref2 = LCI{S2}(m.last_constraint_reference)
    dict2 = constrdict(m, ref2)
    dict2[ref2] = row
    delete!(dict, ref)
    deleteconstraintname!(m, ref)
    return ref2
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: sso_oidc
using AWS.Compat
using AWS.UUIDs

"""
    create_token(client_id, client_secret, device_code, grant_type)
    create_token(client_id, client_secret, device_code, grant_type, params::Dict{String,<:Any})

Creates and returns an access token for the authorized client. The access token issued will
be used to fetch short-term credentials for the assigned roles in the AWS account.

# Arguments
- `client_id`: The unique identifier string for each client. This value should come from
  the persisted result of the RegisterClient API.
- `client_secret`: A secret string generated for the client. This value should come from
  the persisted result of the RegisterClient API.
- `device_code`: Used only when calling this API for the device code grant type. This
  short-term code is used to identify this authentication attempt. This should come from an
  in-memory reference to the result of the StartDeviceAuthorization API.
- `grant_type`: Supports grant types for authorization code, refresh token, and device code
  request.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"code"`: The authorization code received from the authorization service. This parameter
  is required to perform an authorization grant request to get access to a token.
- `"redirectUri"`: The location of the application that will receive the authorization
  code. Users authorize the service to send the request to this location.
- `"refreshToken"`: The token used to obtain an access token in the event that the access
  token is invalid or expired. This token is not issued by the service.
- `"scope"`: The list of scopes that is defined by the client. Upon authorization, this
  list is used to restrict permissions when granting an access token.
"""
function create_token(
    clientId,
    clientSecret,
    deviceCode,
    grantType;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return sso_oidc(
        "POST",
        "/token",
        Dict{String,Any}(
            "clientId" => clientId,
            "clientSecret" => clientSecret,
            "deviceCode" => deviceCode,
            "grantType" => grantType,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_token(
    clientId,
    clientSecret,
    deviceCode,
    grantType,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return sso_oidc(
        "POST",
        "/token",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "clientId" => clientId,
                    "clientSecret" => clientSecret,
                    "deviceCode" => deviceCode,
                    "grantType" => grantType,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    register_client(client_name, client_type)
    register_client(client_name, client_type, params::Dict{String,<:Any})

Registers a client with AWS SSO. This allows clients to initiate device authorization. The
output should be persisted for reuse through many authentication requests.

# Arguments
- `client_name`: The friendly name of the client.
- `client_type`: The type of client. The service supports only public as a client type.
  Anything other than public will be rejected by the service.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"scopes"`: The list of scopes that are defined by the client. Upon authorization, this
  list is used to restrict permissions when granting an access token.
"""
function register_client(
    clientName, clientType; aws_config::AbstractAWSConfig=global_aws_config()
)
    return sso_oidc(
        "POST",
        "/client/register",
        Dict{String,Any}("clientName" => clientName, "clientType" => clientType);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function register_client(
    clientName,
    clientType,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return sso_oidc(
        "POST",
        "/client/register",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("clientName" => clientName, "clientType" => clientType),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    start_device_authorization(client_id, client_secret, start_url)
    start_device_authorization(client_id, client_secret, start_url, params::Dict{String,<:Any})

Initiates device authorization by requesting a pair of verification codes from the
authorization service.

# Arguments
- `client_id`: The unique identifier string for the client that is registered with AWS SSO.
  This value should come from the persisted result of the RegisterClient API operation.
- `client_secret`: A secret string that is generated for the client. This value should come
  from the persisted result of the RegisterClient API operation.
- `start_url`: The URL for the AWS SSO user portal. For more information, see Using the
  User Portal in the AWS Single Sign-On User Guide.

"""
function start_device_authorization(
    clientId, clientSecret, startUrl; aws_config::AbstractAWSConfig=global_aws_config()
)
    return sso_oidc(
        "POST",
        "/device_authorization",
        Dict{String,Any}(
            "clientId" => clientId, "clientSecret" => clientSecret, "startUrl" => startUrl
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function start_device_authorization(
    clientId,
    clientSecret,
    startUrl,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return sso_oidc(
        "POST",
        "/device_authorization",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "clientId" => clientId,
                    "clientSecret" => clientSecret,
                    "startUrl" => startUrl,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
export ft_setflowcontrol

export FT_FLOW_NONE, FT_FLOW_RTS_CTS, FT_FLOW_DTR_DSR, FT_FLOW_XON_XOFF

#Flow Control (see ft_setflowcontrol)
const FT_FLOW_NONE = 0x0000
const FT_FLOW_RTS_CTS = 0x0100
const FT_FLOW_DTR_DSR = 0x0200
const FT_FLOW_XON_XOFF = 0x0400

function ft_setflowcontrol(ft_handle::Culong, flowcontrol::Integer,
                           xon::UInt8 = 0x11, xoff::UInt8 = 0x13)
  @assert flowcontrol == FT_FLOW_NONE ||
          flowcontrol == FT_FLOW_RTS_CTS ||
          flowcontrol == FT_FLOW_DTR_DSR ||
          flowcontrol == FT_FLOW_XON_XOFF "invalid flow control"
  ft_status = ccall((:FT_SetFlowControl, d2xx),
                     Cuint,
                     (Culong, Cushort, Cuchar, Cuchar),
                     ft_handle, flowcontrol, xon, xoff)
  checkstatus(ft_status)
  return nothing
end


text:
"""
    RecursiveApply

This module contains operators to recurse over nested `Tuple`s or `NamedTuple`s.

To extend to another type `T`, define `RecursiveApply.rmap(fn, args::T...)`
"""
module RecursiveApply

export ⊞, ⊠, ⊟

"""
    rmap(fn, X...)

Recursively apply `fn` to each element of `X`
"""
rmap(fn::F, X) where {F} = fn(X)
rmap(fn::F, X, Y) where {F} = fn(X, Y)
rmap(fn::F, X::Tuple) where {F} = map(x -> rmap(fn, x), X)
rmap(fn::F, X::Tuple, Y::Tuple) where {F} = map((x, y) -> rmap(fn, x, y), X, Y)
rmap(fn::F, X::NamedTuple{names}) where {F, names} =
    NamedTuple{names}(rmap(fn, Tuple(X)))
rmap(fn::F, X::NamedTuple{names}, Y::NamedTuple{names}) where {F, names} =
    NamedTuple{names}(rmap(fn, Tuple(X), Tuple(Y)))

"""
    rmaptype(fn, T)

The return type of `rmap(fn, X::T)`.
"""
rmaptype(fn::F, ::Type{T}) where {F, T} = fn(T)
rmaptype(fn::F, ::Type{T}) where {F, T <: Tuple} =
    Tuple{map(fn, tuple(T.parameters...))...}
rmaptype(
    fn::F,
    ::Type{T},
) where {F, T <: NamedTuple{names, tup}} where {names, tup} =
    NamedTuple{names, rmaptype(fn, tup)}

"""
    rmul(w, X)
    w ⊠ X

Recursively scale each element of `X` by `w`.
"""
rmul(w::Number, X) = rmap(x -> w * x, X)
rmul(X, w::Number) = rmap(x -> x * w, X)
rmul(w1::Number, w2::Number) = w1 * w2
const ⊠ = rmul

"""
    radd(X, Y)
    X ⊞ Y

Recursively add elements of `X` and `Y`.
"""
radd(X, Y) = rmap(+, X, Y)
const ⊞ = radd

"""
    rsub(X, Y)
    X ⊟ Y

Recursively subtract elements of `Y` from `X`.
"""
rsub(X) = rmap(-, X)
rsub(X, Y) = rmap(-, X, Y)
const ⊟ = rsub

rdiv(X, w::Number) = rmap(x -> x / w, X)

"""
    rmuladd(w, X, Y)

Recursively add elements of `w * X + Y`.
"""
rmuladd(w::Number, X, Y) = rmap((x, y) -> muladd(w, x, y), X, Y)
rmuladd(X, w::Number, Y) = rmap((x, y) -> muladd(x, w, y), X, Y)
rmuladd(w::Number, x::Number, y::Number) = muladd(w, x, y)

"""
    rmatmul1(W, S, i, j)

Recursive matrix product along the 1st dimension of `S`. Equivalent to:

    mapreduce(⊠, ⊞, W[i,:], S[:,j])

"""
function rmatmul1(W, S, i, j)
    Nq = size(W, 2)
    r = W[i, 1] ⊠ S[1, j]
    for ii in 2:Nq
        r = rmuladd(W[i, ii], S[ii, j], r)
    end
    return r
end

"""
    rmatmul2(W, S, i, j)

Recursive matrix product along the 2nd dimension `S`. Equivalent to:

    mapreduce(⊠, ⊞, W[j,:], S[i, :])

"""
function rmatmul2(W, S, i, j)
    Nq = size(W, 2)
    r = W[j, 1] ⊠ S[i, 1]
    for jj in 2:Nq
        r = rmuladd(W[j, jj], S[i, jj], r)
    end
    return r
end

end # module


text:
# Multinest Dumper Routine
function multinest_default_dumper(physlive::Matrix{T},posterior::Matrix{T},
                          paramconstr::Matrix{T}, maxloglike::T,
                          logz::T, inslogz::T, logzerr::T) where T <: AbstractFloat
end

function multinest_dumper(nsamples_::Ptr{Cint},
                          nlive_::Ptr{Cint},
                          npar_::Ptr{Cint},
                          physlive_::Ptr{Ptr{Cdouble}},
                          posterior_::Ptr{Ptr{Cdouble}},
                          paramconstr_::Ptr{Ptr{Cdouble}},
                          maxloglike_::Ptr{Cdouble},
                          logZ_::Ptr{Cdouble},
                          insLogZ_::Ptr{Cdouble},
                          logZerr_::Ptr{Cdouble},
                          nested_::Ptr{Nothing})

    nsamples = unsafe_load(nsamples_)
    nlive = unsafe_load(nlive_)
    npar = unsafe_load(npar_)
    physlive = unsafe_wrap(Array,unsafe_load(physlive_),(nlive, Int32(npar+1)))
    posterior = unsafe_wrap(Array,unsafe_load(posterior_),(nsamples, Int32(npar+2)))
    paramconstr = unsafe_wrap(Array,unsafe_load(paramconstr_),(npar, Int32(4)))
    maxloglike = unsafe_load(maxloglike_)
    logZ = unsafe_load(logZ_)
    insLogZ = unsafe_load(insLogZ_)
    logZerr = unsafe_load(logZerr_)
    nested = unsafe_pointer_to_objref(nested_)::MultiNest
    nested.dumper(physlive, posterior, paramconstr, maxloglike, logZ, insLogZ, logZerr)
    return
end

# Multinest Specific Settings
struct MultiNest
    loglike::Function
    root::String
    ins::Cint
    mmodal::Cint
    ceff::Cint
    efr::Cdouble
    mode_ztol::Cdouble
    ndims::Cint
    npar::Cint
    nclspar::Cint
    updint::Cint
    maxmodes::Cint
    pwrap::Vector{Cint}
    seed::Cint
    outfile::Cint
    initmpi::Cint
    logzero::Cdouble
    maxiter::Cint
    dumper::Function
end

function MultiNest(loglike::Function, nparams::Int, name::String;
                   prior::Function = identity,
                   importance_sampling::Bool = true,
                   multimodal::Bool = true,
                   max_modes::Int = 100,
                   mode_ztol::Float64 = -1e90,
                   nderived::Int = 0,
                   ncluster_dims::Int = 0,
                   mode_npar::Int = (nparams - nderived),
                   constant_efficiency::Bool = false,
                   efficiency_rate::AbstractFloat = 0.8,
                   seed::Int = -1,
                   update_every::Int = 1000,
                   write_output_files::Bool = true,
                   work_dir::String = ".",
                   init_mpi::Bool = true,
                   periodic::Vector{Bool} = fill(false, nparams - nderived),
                   maxiter::Int = 0,
                   dumper::Function = multinest_default_dumper)

    ndims = nparams - nderived

    # Verify Inputs
    nparams > 0 || error("nparams must be greater than zero")
    nparams > nderived || error("nderived must be less than nparams")
    ncluster_dims < ndims || error("ncluster_dims must be less than (nparams - nderived)")
    0 < efficiency_rate <= 1 || error("efficiency_rate must be in the range (0,1]")
    max_modes > 0 || error("max_modes must be greater than 0")
    length(periodic) == ndims || error("size(periodic) != nparams - nderived")
    maxiter >= 0 || error("max_iter must be non-negative")
    update_every > 0 || error("update_every must be non-negative")
    seed >= -1 || error("seed must be positive integer")

    MultiNest(x->loglike(prior(x)), rpad(joinpath(work_dir,name*"_"), 100, ' '), Cint(importance_sampling),
              Cint(multimodal), Cint(constant_efficiency), Cdouble(efficiency_rate),
              Cdouble(mode_ztol), Cint(ndims), Cint(nparams), Cint(ncluster_dims),
              Cint(update_every), Cint(max_modes), Cint.(periodic), Cint(seed),
              Cint(write_output_files), Cint(init_mpi), Cdouble(nextfloat(-Inf)),
              Cint(maxiter), dumper)

end

function multinest_loglike(cube_::Ptr{Cdouble}, ndim_::Ptr{Cint},
                           npar_::Ptr{Cint}, lnew_::Ptr{Cdouble},
                           nested_::Ptr{Nothing})
    ndim = unsafe_load(ndim_)
    npar = unsafe_load(npar_)
    cube = unsafe_wrap(Array, cube_, npar)
    nested = unsafe_pointer_to_objref(nested_)::MultiNest
    lnew = nested.loglike(cube)
    unsafe_store!(lnew_, lnew)
    return
end

function run(S::MultiNest; nlive::Int=25*S.ndims, ztol=0.1, resume=false, verbosity=0)

    loglike = @cfunction(multinest_loglike, Nothing, (Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Nothing}))

    dumper = @cfunction(multinest_dumper, Nothing, (Ptr{Cint}, Ptr{Cint}, Ptr{Cint},
                                         Ptr{Ptr{Cdouble}}, Ptr{Ptr{Cdouble}},
                                         Ptr{Ptr{Cdouble}}, Ptr{Cdouble},
                                         Ptr{Cdouble}, Ptr{Cdouble},Ptr{Cdouble},
                                         Ptr{Nothing}))

    root_dir = dirname(S.root)
    println(root_dir)
    isdir(root_dir) || mkdir(root_dir)

    ccall(multinest_ptr, Nothing, ( Ref{Cint}, Ref{Cint}, Ref{Cint},
          Ref{Cint}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cint}, Ref{Cint},
          Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cdouble}, Ptr{UInt8},
          Ref{Cint}, Ptr{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint},
          Ref{Cdouble}, Ref{Cint}, Ptr{Nothing}, Ptr{Nothing}, Any ),
          S.ins, S.mmodal, S.ceff, nlive, ztol, S.efr,
          S.ndims, S.npar, S.nclspar, S.maxmodes, S.updint, S.mode_ztol,
          S.root, S.seed, S.pwrap, verbosity, resume, S.outfile, S.initmpi,
          S.logzero, S.maxiter, loglike, dumper, S)
end