text:
function construct_device!(ps_m::CanonicalModel,
                           device::Type{R},
                           device_formulation::Type{D},
                           system_formulation::Type{S},
                           sys::PSY.System,
                           time_range::UnitRange{Int64};
                           kwargs...) where {R <: PSY.RenewableGen,
                                             D <: AbstractRenewableDispatchForm,
                                             S <: PM.AbstractPowerFormulation}

    if !isa(sys.generators.renewable, Nothing)

        parameters = get(kwargs, :parameters, true)

        fixed_resources = [fs for fs in sys.generators.renewable if isa(fs,PSY.RenewableFix)]

        controllable_resources = [fs for fs in sys.generators.renewable if !isa(fs,PSY.RenewableFix)]

        if !isempty(controllable_resources)

            #Variables
            activepower_variables(ps_m, controllable_resources, time_range);

            reactivepower_variables(ps_m, controllable_resources, time_range);

            #Constraints
            activepower_constraints(ps_m, controllable_resources, device_formulation, system_formulation, time_range, parameters)

            reactivepower_constraints(ps_m, controllable_resources, device_formulation, system_formulation, time_range)

            #Cost Function
            cost_function(ps_m, controllable_resources, device_formulation, system_formulation)

        else
            @warn("The Data Doesn't Contain Controllable Renewable Resources, Consider Changing the Device Formulation to RenewableFixed")
        end

        #add to expression

        if !isempty(fixed_resources)
            nodal_expression(ps_m, fixed_resources, system_formulation, time_range, parameters)
        end

    end

    return

end

function construct_device!(ps_m::CanonicalModel,
                           device::Type{R},
                           device_formulation::Type{D},
                           system_formulation::Type{S},
                           sys::PSY.System,
                           time_range::UnitRange{Int64};
                           kwargs...) where {R <: PSY.RenewableGen,
                                             D <: AbstractRenewableDispatchForm,
                                             S <: PM.AbstractActivePowerFormulation}

    if !isa(sys.generators.renewable, Nothing)

        parameters = get(kwargs, :parameters, true)

        fixed_resources = [fs for fs in sys.generators.renewable if isa(fs,PSY.RenewableFix)]

        controllable_resources = [fs for fs in sys.generators.renewable if !isa(fs,PSY.RenewableFix)]

        if !isempty(controllable_resources)

            #Variables
            activepower_variables(ps_m, controllable_resources, time_range)

            #Constraints
            activepower_constraints(ps_m, controllable_resources, device_formulation, system_formulation, time_range, parameters)

            #Cost Function
            cost_function(ps_m, controllable_resources, device_formulation, system_formulation)

        else
            @warn("The Data Doesn't Contain Controllable Renewable Resources, Consider Changing the Device Formulation to RenewableFixed")

        end

        #add to expression

        if !isempty(fixed_resources)
            nodal_expression(ps_m, fixed_resources, system_formulation, time_range, parameters)
        end

    end

    return

end

function construct_device!(ps_m::CanonicalModel,
                           device::Type{R},
                           device_formulation::Type{PSI.RenewableFixed},
                           system_formulation::Type{S},
                           sys::PSY.System,
                           time_range::UnitRange{Int64};
                           kwargs...) where {R <: PSY.RenewableGen,
                                             S <: PM.AbstractPowerFormulation}

    if !isa(sys.generators.renewable, Nothing)

        parameters = get(kwargs, :parameters, true)

        nodal_expression(ps_m, sys.generators.renewable, system_formulation, time_range, parameters)
    end

    return

end


text:
import Oceananigans.Grids: required_halo_size

struct ScalarDiffusivity{TD, F, N, K} <: AbstractScalarDiffusivity{TD, F}
    ν :: N
    κ :: K

    function ScalarDiffusivity{TD, F}(ν::N, κ::K) where {TD, F, N, K}
        return new{TD, F, N, K}(ν, κ)
    end
end

"""
ScalarDiffusivity(time_disc=ExplicitTimeDiscretization,
                 formulation=ThreeDimensionalFormulation, FT=Float64;
                 ν=0, κ=0,
                 discrete_form = false)

Returns parameters for an isotropic diffusivity model with viscosity `ν`
and thermal diffusivities `κ` for each tracer field in `tracers`
`ν` and the fields of `κ` may be constants, arrays, fields or
- functions of `(x, y, z, t)` if `discrete_diffusivity = false`
- functions of `(i, j, k, grid, LX, LY, LZ)` with `LX`, `LY` and `LZ` are either `Face()` or `Center()` if
  `discrete_form = true`.

`κ` may be a `NamedTuple` with fields corresponding to each tracer, or a
single number to be a applied to all tracers.
"""

  VerticalScalarDiffusivity(time_disc=ExplicitTimeDiscretization(), FT::DataType=Float64; kwargs...) = ScalarDiffusivity(time_disc, VerticalFormulation(), FT; kwargs...)
HorizontalScalarDiffusivity(time_disc=ExplicitTimeDiscretization(), FT::DataType=Float64; kwargs...) = ScalarDiffusivity(time_disc, HorizontalFormulation(), FT; kwargs...)

# Aliases that allow specify the floating type, assuming that the discretization is Explicit in time
          ScalarDiffusivity(FT::DataType; kwargs...) = ScalarDiffusivity(ExplicitTimeDiscretization(), ThreeDimensionalFormulation(), FT; kwargs...)
  VerticalScalarDiffusivity(FT::DataType; kwargs...) = ScalarDiffusivity(ExplicitTimeDiscretization(), VerticalFormulation(), FT; kwargs...)
HorizontalScalarDiffusivity(FT::DataType; kwargs...) = ScalarDiffusivity(ExplicitTimeDiscretization(), HorizontalFormulation(), FT; kwargs...)

function ScalarDiffusivity(time_disc=ExplicitTimeDiscretization(),
                           formulation=ThreeDimensionalFormulation(), FT=Float64;
                           ν=0, κ=0,
                           discrete_form = false)

    if formulation == HorizontalFormulation() && time_discretization == VerticallyImplicitTimeDiscretization()
        throw(ArgumentError("VerticallyImplicitTimeDiscretization is only supported for `HorizontalFormulation` or `ThreeDimensionalFormulation`"))
    end
    κ = convert_diffusivity(FT, κ, Val(discrete_form))
    ν = convert_diffusivity(FT, ν, Val(discrete_form))
    return ScalarDiffusivity{typeof(time_disc), typeof(formulation)}(ν, κ)
end

required_halo_size(closure::ScalarDiffusivity) = 1

function with_tracers(tracers, closure::ScalarDiffusivity{TD, F}) where {TD, F}
    κ = tracer_diffusivities(tracers, closure.κ)
    return ScalarDiffusivity{TD, F}(closure.ν, κ)
end

@inline viscosity(closure::ScalarDiffusivity, args...) = closure.ν
@inline diffusivity(closure::ScalarDiffusivity, ::Val{tracer_index}, args...) where tracer_index = closure.κ[tracer_index]

calculate_diffusivities!(diffusivities, ::ScalarDiffusivity, args...) = nothing

Base.show(io::IO, closure::ScalarDiffusivity{TD, F})  where {TD, F}=
    print(io, "ScalarDiffusivity:\n",
              "ν=$(closure.ν), κ=$(closure.κ)\n",
              "time discretization: $TD\n",
              "formulation: $F")

ext: jl


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    ClassificationTask(x, y)

A classification task consists of finding a function `f` such that `y ~ f(x)`
for all training examples `(x,y)` with `y` a categorical variable.
"""
struct ClassificationTask{N} <: AbstractLearningTask
  features::NTuple{N,Symbol}
  label::Symbol
end

ClassificationTask(x::Symbol, y::Symbol) = ClassificationTask{1}((x,), y)

ClassificationTask(x::AbstractVector{Symbol}, y::Symbol) =
  ClassificationTask{length(x)}(Tuple(x), y)

issupervised(task::ClassificationTask) = true

"""
    features(task)

Return the features of the classification `task`.
"""
features(task::ClassificationTask) = task.features

"""
    label(task)

Return the label of the classification `task`.
"""
label(task::ClassificationTask) = task.label

# ------------
# IO methods
# ------------
function Base.show(io::IO, task::ClassificationTask)
  x = task.features
  y = task.label
  lhs = length(x) > 1 ? "("*join(x, ", ")*")" : "$(x[1])"
  print(io, "Classification $lhs → $y")
end


text:
"""
    Predicted state in here
"""
function predict_output(filter, x, k)
    H = get_H(filter, k)
    y_pred = H*x
end
"""
    Aka propagation. Time  update. Take the filtered state and propagate.
"""
function time_update(model, x::MvNormal, k)
    F = jacobian_x(model, x, k)
    v = get_process_noise(model, x, k) # Zero mean RV

    x_pred = F_k*x + v

    return x_pred
end
function measurement_update(model, y, x_pred, k)
    y_pred = predict_output(model, x_pred, k)
    e = y - y_pred # Innovation

    P_pred = cov(x_pred)
    S = cov(e)

    H = jacobian_x(model, k)
    K = P_pred*H'/S

    x = mean(x_pred) + K*mean(e)
    P = (I - K*H)*P_pred
    # Trick  to make the  covariance matrix symmetric
    return MvNormal(x, 0.5*(P + P')), e
end


text:
export
    VectorizedNGSIMEnv,
    reset,
    step,
    observation_space_spec,
    action_space_spec,
    obs_names,
    render,
    vectorized,
    num_envs

#=
Description:
    NGSIM env that plays NGSIM trajectories, allowing the agent to take the place
    of one of the vehicles in the trajectory
=#
type VectorizedNGSIMEnv <: Env
    n_envs::Int
    envs::Vector{NGSIMEnv}
    x::Array{Float64}
    r::Vector{Float64}
    dones::Vector{Bool}

    function VectorizedNGSIMEnv(
            params::Dict;
            n_envs::Int = 10,
            reclength::Int = 10,
            Δt::Float64 = .1,
            primesteps::Int = 50,
            H::Int = 50,
            terminate_on_collision::Bool = true,
            terminate_on_off_road::Bool = true,
            render_params::Dict = Dict("zoom"=>5., "viz_dir"=>"/tmp"))
        param_keys = keys(params)
        @assert in("trajectory_filepaths", param_keys)

        # optionally overwrite defaults
        n_envs = get(params, "n_envs", n_envs)
        reclength = get(params, "reclength", reclength)
        primesteps = get(params, "primesteps", primesteps)
        H = get(params, "H", H)
        for (k,v) in get(params, "render_params", render_params)
            render_params[k] = v
        end
        terminate_on_collision = get(params, "terminate_on_collision", terminate_on_collision)
        terminate_on_off_road = get(params, "terminate_on_off_road", terminate_on_off_road)

        # load trajdatas
        trajdatas, trajinfos, roadways = load_ngsim_trajdatas(
            params["trajectory_filepaths"],
            minlength=primesteps + H
        )

        # build envs
        envs = Vector{NGSIMEnv}(n_envs)
        for i in 1:n_envs
            envs[i] = NGSIMEnv(
                params,
                trajdatas=trajdatas,
                trajinfos=trajinfos,
                roadways=roadways,
                reclength=reclength,
                Δt=Δt,
                primesteps=primesteps,
                H=H,
                terminate_on_collision=terminate_on_collision,
                terminate_on_off_road=terminate_on_off_road,
                render_params=render_params
            )
        end
        n_features = length(obs_names(envs[1]))
        x = zeros(n_envs, n_features)
        r = zeros(n_envs)
        dones = Vector{Bool}(n_envs)

        return new(n_envs, envs, x, r, dones)
    end
end
function reset(env::VectorizedNGSIMEnv, dones::Vector{Bool} = fill!(Vector{Bool}(env.n_envs), true))
    for i in 1:env.n_envs
        if dones[i]
            env.x[i, :] = reset(env.envs[i])
        end
    end
    return deepcopy(env.x)
end
function Base.step(env::VectorizedNGSIMEnv, actions::Array{Float64})
    infos = Vector{Dict}(env.n_envs)
    for i in 1:env.n_envs
            env.x[i, :], env.r[i], env.dones[i], infos[i] = step(env.envs[i], actions[i, :])
    end
    # vectorized sampler does not call reset on the environment
    # but expects the environment to handle resetting, so do that here
    # note: this mutates env.x in order to return the correct obs when resetting
    reset(env, env.dones)
    return deepcopy(env.x), env.r, env.dones, stack_tensor_dict_list(infos)
end

observation_space_spec(env::VectorizedNGSIMEnv) = observation_space_spec(env.envs[1])
action_space_spec(env::VectorizedNGSIMEnv) = action_space_spec(env.envs[1])
obs_names(env::VectorizedNGSIMEnv) = obs_names(env.envs[1])
render(env::VectorizedNGSIMEnv) = render(env.envs[1])
vectorized(env::VectorizedNGSIMEnv) = true
num_envs(env::VectorizedNGSIMEnv) = env.n_envs


text:
export
    BP,
    infer

"""
    abstract type BP <: InstantAlgorithm

Belief Propagation algorithm
"""
abstract type BP <: InstantAlgorithm end

function answer(::Marginal, ::BP, runtime::Runtime, instance::VariableInstance, bounds = false)
    if bounds
        error("BP cannot provide bounded answers")
    end
    return get_belief(runtime, instance)
end

function infer(algorithm::BP, runtime::InstantRuntime,
    evidence::Dict{Symbol, Score} = Dict{Symbol, Score}(),
    interventions::Dict{Symbol, Dist} = Dict{Symbol, Dist}(),
    placeholder_beliefs = get_placeholder_beliefs(runtime))
    net = get_network(runtime)
    ensure_all!(runtime)
    order = topsort(get_initial_graph(net))
    set_ranges!(runtime, evidence, algorithm.default_range_size, 1, order, placeholder_beliefs)
    for ph in get_placeholders(net)
        if !(ph.name in keys(placeholder_beliefs))
            error("Placeholder ", ph.name, " does not have a belief")
        end
        pi = placeholder_beliefs[ph.name]
        for ch in get_children(net, ph)
            set_message!(runtime, ph, ch, :pi_message, pi)
        end
    end
    for (n,e) in evidence
        v = get_node(net, n)
        inst = current_instance(runtime, v)
        post_evidence!(runtime, inst, e)
    end
    for (n,i) in interventions
        v = get_node(net, n)
        inst = current_instance(runtime, v)
        post_intervention!(runtime, inst, i)
    end
    run_bp(algorithm, runtime)
end

# Run a backward step of BP for a single variable.
function _backstep(runtime, var, ranges)
    network = get_network(runtime)
    inst = current_instance(runtime, var)
    evidence = get_evidence(runtime, inst)
    sf = get_sfunc(inst)
    chs = get_children(get_network(runtime), var)
    msgs = collect_messages(runtime, chs, var, :lambda_message)
    O = output_type(sf)
    range::Vector{<:O} = [convert(output_type(sf), r) for r in ranges[var.name]]

    if isempty(msgs)
        ilams = Score{output_type(sf)}[]
        if !isnothing(evidence)
            push!(ilams, evidence)
        end
        lam = compute_lambda(sf, range, ilams)
    else
        incoming_lams = Score[]
        for m in msgs
            m1::Score = m
            push!(incoming_lams, m1)
        end
        if !isnothing(evidence)
            push!(incoming_lams, evidence)
        end
        lam = compute_lambda(sf, range, incoming_lams)
    end

    set_value!(runtime, inst, :lambda, lam)
    pars = get_parents(network, var)
    incoming_pis::Vector{Dist} = collect_messages(runtime, pars, var, :pi_message)
    parranges = [ranges[p.name] for p in pars]
    outgoing_lams = outgoing_lambdas(sf, lam, range, tuple(parranges...), tuple(incoming_pis...))
    distribute_messages!(runtime, var, pars, :lambda_message, outgoing_lams)
end

# Run a forward step of BP for a single variable
# The procedure is slightly different for the first pass, so use the firstpass flag
function _forwardstep(runtime, var, ranges, firstpass)
    network = get_network(runtime)
    inst = current_instance(runtime, var)
    sf = get_sfunc(inst)
    pars = get_parents(network, var)
    O = output_type(sf)
    range = convert(Vector{O}, ranges[var.name])
    parranges = [ranges[p.name] for p in pars]
    incoming_pis = tuple(collect_messages(runtime, pars, var, :pi_message)...)
    if has_intervention(runtime, inst)
        intervention = get_intervention(runtime, inst)
        pi = intervention
    elseif firstpass && has_evidence(runtime, inst)
        # In the first pass, evidence is treated like lambda messages from elsewhere
        evidence = get_evidence(runtime, inst)
        pi = Cat(range, [get_score(evidence, x) for x in range])
    else
        pi = compute_pi(sf, range, tuple(parranges...), incoming_pis)
    end
    set_value!(runtime, inst, :pi, pi)

    chs = get_children(network, var)
    if firstpass
        outpis = fill(pi, length(chs))
    else
        lam = get_value(runtime, inst, :lambda)
        bel = compute_bel(sf, range, pi, lam)
        post_belief!(runtime, inst, bel)
        incoming_lams::Vector{Score} = collect_messages(runtime, chs, var, :lambda_message)
        outpis = outgoing_pis(sf, range, bel, incoming_lams)
    end
    distribute_messages!(runtime, var, chs, :pi_message, outpis)
end


text:
using HopfNormalForm

Flutter_opt.optimize_NLstiff(3.862,13.956)
Flutter_opt.optimize_NLstiff2(4.538,19.9976)


text:
export map_components
using Base.Iterators: product

@with_kw struct MultiPriorTracking{C} <: Tracking
  cohere::C
  time_constants::Array{typeof(1.0s)}
  time_constant_bias::Vector{Float64}
  source_priors::AxisArray
  freq_prior
  max_sources::Int = 4
  min_norm::Float64 = Inf
  normalize::Bool = false
end
function Tracking(C,::Val{:multi_prior};
                  time_constants=[4s],
                  time_constant_bias=zeros(length(time_constants)),
                  source_prior_sds,
                  source_prior_strengths=[1.0s],
                  source_prior_sd_bias=zeros(length(source_prior_sds)),
                  source_prior_strength_bias=
                    zeros(length(source_prior_strengths)),
                  freq_ridge=0.0, scale_ridge=0.0,
                  ridge_threshold=0.05,
                  freq_prior_N = 2, freq_prior_bias = 0,
                  normalize=false,min_norm=Inf,
                  params...)
  checknorm(normalize,min_norm)
  source_prior_Ns = asseconds.(source_prior_strengths)./Δt(C)
  pr = product(zip(source_prior_sds, source_prior_sd_bias),
               zip(source_prior_Ns, source_prior_strength_bias))

  source_priors = if iszero(freq_ridge) && iszero(scale_ridge)
    vals = [(isonorm(sd, N, (size(C,2),size(C,3))), sd_b + N_b)
                            for ((sd,sd_b),(N,N_b)) in pr]
    prior_params = map(x -> map(a -> a[1],x),pr)
    AxisArray(vec(vals), Axis{:params}(vec(prior_params)))
  else
    vals = [(ridgenorm(sd, N, (size(C,2),size(C,3)),
                       freq=freq_ridge,scale=scale_ridge,
                       threshold=ridge_threshold), sd_b + N_b)
            for ((sd,sd_b),(N,N_b)) in pr]
    prior_params = map(x -> map(a -> a[1],x),pr)
    AxisArray(vec(vals), Axis{:params}(vec(prior_params)))
  end

  freq_prior = freqprior(freq_prior_bias,freq_prior_N)
  MultiPriorTracking(;cohere=getmeta(C), source_priors=source_priors,
                     freq_prior=freq_prior,
                     time_constants=asseconds.(time_constants),
                     time_constant_bias=time_constant_bias,
                     normalize=normalize,min_norm=min_norm,
                     params...)
end

function expand_params(params::MultiPriorTracking)
  pr = params.source_priors

  AxisArray([(PriorTracking(params.cohere,tc,prior,params.freq_prior,
                            params.max_sources,params.min_norm), tb + pb)
             for (tc,tb) in zip(params.time_constants,
                                params.time_constant_bias)
             for (prior,pb) in pr],
            Axis{:params}([(tc,sd,N) for tc in params.time_constants
                           for (sd,N) in axisvalues(params.source_priors)[1]]))
end

function nitr(C::Coherence,params::MultiPriorTracking)
  ntimes(C) * length(params.time_constants) * length(params.source_priors)
end

function track(C::Coherence,params::MultiPriorTracking,progressbar=true,
               progress=track_progress(progressbar,nitr(C,params),"multi-prior"))

  C_ = prepare_coherence(C,params.min_norm)
  all_params = expand_params(params)
  S = Array{SourceTracking}(undef,size(all_params,1))
  lp = Array{Array{Float64}}(undef,size(all_params,1))

  # in my tests, using @Threads does not seem to help hereo This might be
  # related to a now fixed bug in Julia
  # (https://github.com/JuliaLang/julia/issues/17395) or it could because the
  # operation is memory bound; in either case, I haven't gone back to test this
  # and see if threading works any better now
  #=@Threads.threads=# for (i,(p,bias)) in collect(enumerate(all_params))
    S[i], lp[i] = track(C_,p,true,progress)
    lp[i] .+= bias
  end

  (AxisArray(S, AxisArrays.axes(all_params,1)),
   AxisArray(hcat(lp...), AxisArrays.axes(C,1), AxisArrays.axes(all_params,1)))
end

function map_components(fn,tracks::AxisArray{<:SourceTracking},
                        tracks_lp::AxisArray{<:Float64};
                        window=500ms,step=250ms)
  windows = windowing(tracks[1],length=window,step=step)

  result = map(enumerate(windows)) do (i,ixs)
    best_track = argmax(dropdims(mean(Array(tracks_lp[ixs,:]),dims=1),dims=1))
    fn(tracks[best_track][Axis{:time}(ixs)])
  end

  AxisArray(result,AxisArrays.axes(windows,Axis{:time}))
end

function mask(sp::CorticalSpecrtralTemporalResponses.AuditorySpectrogram,
              tracks::AxisArray{<:SourceTracking},
              tracks_lp::AxisArray{<:Float64},
              settings;progressbar=false,kwds...)
  scales = settings["scales"]["values"] .* cycoct
  freql,freqh = settings["rates"]["freq_limits"] .* Hz

  cr = cortical(sp,scales=scales,progressbar=progressbar)
  cr = cr[:,:,freql .. freqh]

  mask(cr,tracks,tracks_lp;progressbar=progressbar,kwds...)
end

function mask(cr::MetaUnion{AxisArray},
              tracks::AxisArray{<:SourceTracking},
              tracks_lp::AxisArray{<:Float64},
              order=1;window=500ms,mask_wait=0,
              delta=250ms,progressbar=false)

  @assert axisdim(cr,Axis{:time}) == 1
  @assert axisdim(cr,Axis{:scale}) == 2
  @assert axisdim(cr,Axis{:freq}) == 3
  @assert size(cr)[2:3] == size(tracks[1])[1:2] "Dimension mismatch"

  windows = windowing(tracks[1],length=window,step=delta)

  progress = progressbar ? Progress(length(windows),"Masking: ") : nothing
  mask_helper(cr,tracks,tracks_lp,order,windows,mask_wait,progress)
end

function mask_helper(cr,tracks,tracks_lp,order,windows,mask_wait,progress)
  y = zero(cr)
  norm = similar(cr,real(eltype(cr)))
  norm .= zero(eltype(norm))

  cohere_windows =
    collect(windowing(cr,tracks[1].cohere))

  for (i,ixs) = enumerate(windows)
    best_track = argmax(dropdims(mean(Array(tracks_lp[ixs,:]),dims=1),dims=1))
    components = view(tracks[best_track],:,:,:,ixs)
    sorting = sortperm(component_means(components),rev=true)
    c = sorting[order]

    for (ti,t) in enumerate(ixs)
      if i <= mask_wait
        y[cohere_windows[t],:,:] .+= 1
        norm[cohere_windows[t],:,:] .+= 1
      else
        resh = reshape(view(components,:,:,c,ti),1,size(y,2),size(y,3))
        y[cohere_windows[t],:,:] .+= resh
        norm[cohere_windows[t],:,:] .+= sqrt.(mean(x -> x^2,resh))
      end
    end

    next!(progress)
  end
  y ./= max.(1,norm)
  y ./= maximum(abs,y)
  y .*= cr

  y
end

max_stars_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_issues_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_forks_repo_name: JuliaTagBot/AuditoryBistabilityLE


text:
using SimpleNeRF
using Test

@testset "SimpleNeRF.jl" begin
    # Write your tests here.
    @testset "data" begin
        include("test_data.jl")
    end
    @testset "model" begin
        include("test_model.jl")
    end
    @testset "render" begin
        include("test_render.jl")
    end
end


text:
# Bayesian Logistic Regression using the Polya-Gamma augmentation.
#

using LinearAlgebra
using StatsFuns: log1pexp

export BinaryLogisticRegression
export LogisticRegression
export predict
export reorder

# TO REMOVE
export regressors
export encode
export pgaugmentation

#######################################################################
# Helper functions
#

function regressors(hasbias::Bool, X::Matrix{T}) where T <: AbstractFloat
    if hasbias
        # We add an extra dimension always set to "1" to incorporate a bias
        # parameter to the model.
        return vcat(X, ones(T, 1, size(X, 2)))
    end
    return X
end

# Mean of PolyaGamma distributions
function pgmean(c::Union{Vector{T}, Matrix{T}}; b::Real = 1) where T <: AbstractFloat
    retval = (b ./ (2 * c)) .* tanh.(c ./ 2)

    # When c = 0 the mean is not defined but can be extended by
    # continuity observing that lim_{x => 0} (e^(x) - 1) / x = 0             !
    # which lead to the mean = b / 4
    idxs = isnan.(retval)
    retval[idxs] .= b / 4
    return retval
end

# Log-normalizer of PolyaGamma distributions
function pglognorm(c::Union{Vector{T}, Matrix{T}}; b::Real = 1) where T <: AbstractFloat
    return -b .* (log1pexp.(c) .- log(2) .- c ./ 2)
end

# Compute the necessary results of the PolyaGamma augmentation:
#   1. mean of the ω variables
#   2. KL divergence between the variational factors q(ω) and p(ω)
function pgaugmentation(Ψ²::Union{Vector{T}, Matrix{T}}; b::Real = 1) where T <: AbstractFloat
    c = sqrt.(Ψ²)
    E_ω = pgmean(c)
    kl = -.5 .* (c.^2) .* E_ω .- pglognorm(c)
    E_ω, kl
end

# For the purpose of operating with matrices the function encodes
# the vector of classes `z` as:
#   * onehot: one-hot encoding for the classes, where the class
#     'n_classes' is encoded as all 0's
#   * N_up_to: for each sample of class c, row i is 1 if i>=c and
#     0 otherwise. Note that if c=n_classes, all rows are 0's
#
# Note: K is the number of classes
function encode(z::Vector{T}, K::Int) where T <: Int
    onehot = zeros(T, K - 1, length(z))
    for i = 1:K-1 onehot[i, z .== i] .= 1 end

    N_up_to = zeros(T, K - 1, length(z))
    for i = 1:K-1 N_up_to[i, z .>= i] .= 1 end

    onehot, N_up_to
end

#######################################################################
# Binary logistic regression

struct BinaryLogisticRegression <: Model
    β::ConjugateParameter{<:Normal}
    hasbias::Bool

    function BinaryLogisticRegression(β::ConjugateParameter{<:Normal}, hasbias::Bool)
        model = new(β, hasbias)
        β.stats = data -> stats_β(model, data...)
        return model
    end
end

function BinaryLogisticRegression(μ₀::Vector{T}, Σ₀::Matrix{T},
                                  μ::Vector{T}, Σ::Matrix{T},
                                  hasbias::Bool) where T <: AbstractFloat
    BinaryLogisticRegression(ConjugateParameter(Normal(μ₀, Σ₀), Normal(μ, Σ)),
                             hasbias)
end

function BinaryLogisticRegression(T = Float64; inputdim::Integer, initσ::Real = 0.0,
                                  pseudocounts::Real = 1.0, hasbias::Bool = true)
    dim = hasbias ? inputdim + 1 : inputdim

    # Prior parameters
    μ₀ = zeros(T, dim)
    Σ₀ = T(1/pseudocounts) * Matrix{T}(I, dim, dim)

    # Initial paremeterization of the posterior
    μ = μ₀ .+ T(initσ) .* randn(T)
    Σ = deepcopy(Σ₀)

    BinaryLogisticRegression(μ₀, Σ₀, μ, Σ, hasbias)
end

function (model::BinaryLogisticRegression)(
    X::Matrix{T1},
    z::Vector{T2},
) where {T1 <: AbstractFloat, T2 <: Real}

    X̂ = regressors(model.hasbias, X)

    # Expectation of the parameters β
    μ, Σ = stdparam(model.β.posterior)
    E_β = μ
    E_ββᵀ = Σ + μ * μ'

    # Pre-compute the terms of the lower-bound
    ψ = E_βᵀX̂ = X̂' * E_β
    ψ² = dropdims(sum(X̂ .* (E_ββᵀ * X̂), dims = 1), dims = 1)

    # PolyaGamma augmentation
    E_ω, kl_ω = pgaugmentation(ψ²)

    (z .- 0.5) .* ψ .-.5 .* E_ω .* ψ² .- log(2) .- kl_ω
end

function stats_β(model::BinaryLogisticRegression,
    X::Matrix{T1},
    z::Vector{T2},
) where {T1 <: AbstractFloat, T2 <: Real}

    X̂ = regressors(model.hasbias, X)

    # Expectation of the parameters β
    μ, Σ = stdparam(model.β.posterior)
    E_β = μ
    E_ββᵀ = Σ + μ * μ'

    # Polya-Gamma augmentation: compute the prior/posterior over ω
    ψ² = dropdims(sum(X̂ .* (E_ββᵀ * X̂), dims = 1), dims = 1)

    # PolyaGamma augmentation
    E_ω, kl_ω = pgaugmentation(ψ²)

    # 1st order statistics
    zstats = (z .- 0.5)
    s1 = zeros(T1, size(X̂, 1))
    for i = 1:length(z)
        s1 .+= X̂[:, i] * zstats[i]
    end

    # 2nd order statistics.
    s2 = zeros(T1, size(X̂, 1), size(X̂, 1))
    for i = 1:length(z)
        x̂ᵢ = X̂[:, i]
        s2 .+= (E_ω[i] * x̂ᵢ) * x̂ᵢ'
    end
    s2 .*= -.5

    vcat(s1, vec(s2))
end

# Predict the classes using the Maximum A Posteriori (MAP) parameters
function predict_map(model::BinaryLogisticRegression, X::Matrix{T}) where T <: AbstractFloat

    X̂ = regressors(model.hasbias, X)
    μ, _ = stdparam(model.β.posterior)
    y = X̂' * μ
    1 ./ (1 .+ exp.(-y))
end

# Predict the classes the Posteriori (MAP) parameters
function predict_marginal(model::BinaryLogisticRegression,
                 X::Matrix{T}; a::T = 0.368) where T <: AbstractFloat
    # Following https://arxiv.org/pdf/1703.00091.pdf, we use the
    # approximation:
    # ⟨σ(ψ)⟩ ≈ σ(μ / √(1 + a * σ²))
    # where a = 0.368
    # In Bishop 4.5.2 a very similar approximation uses a = π / 8 ≈ 0.392...

    X̂ = regressors(model.hasbias, X)
    μ, Σ = stdparam(model.β.posterior)
    E_ββᵀ = Σ + μ * μ'

    ψ_μ = X̂' * μ
    ψ² = dropdims(sum(X̂ .* (E_ββᵀ * X̂), dims = 1), dims = 1)
    ψ_σ² = ψ² .- (ψ_μ .^ 2)

    y = ψ_μ ./ sqrt.(1 .+ a * ψ_σ²)
    1 ./ (1 .+ exp.(-y))
end

function predict(model::BinaryLogisticRegression, X::Matrix{T};
                 marginalize::Bool = true) where T <: AbstractFloat
    if marginalize
        return predict_marginal(model, X)
    end
    return predict_map(model, X)
end

#######################################################################
# K-class Logistic Regression (LR) based on the stick-breaking process
#
# References
# ----------
# [1] Dependent Multinomial Modes Made Easy: Stick-Breaing with the
#     Polya-Gamma Augmentation (https://arxiv.org/pdf/1506.05843.pdf)

struct LogisticRegression{K} <: Model
    stickbreaking::Vector{BinaryLogisticRegression}
    ordering::Dict{Int, Int}

    function LogisticRegression(stickbreaking::Vector{BinaryLogisticRegression},
                                hasbias::Bool)
        K = length(stickbreaking) + 1
        model = new{K}(stickbreaking, Dict(i => i for i in 1:K))

        # We override the stats
        for (k, blr) in enumerate(stickbreaking)
            blr.β.stats = data -> stats_β(model, k, data...)
        end

        return model
    end
end

function LogisticRegression(μ₀::Vector{T}, Σ₀::Matrix{T}, μ::Vector{T},
                            Σ::Matrix{T}, hasbias::Bool; nclasses::Integer
                           ) where T <: AbstractFloat
    stickbreaking = Vector{BinaryLogisticRegression}()
    for i = 1:nclasses - 1
        push!(stickbreaking,
              BinaryLogisticRegression(μ₀, Σ₀, deepcopy(μ), deepcopy(Σ), hasbias))
    end
    LogisticRegression(stickbreaking, hasbias)
end

function LogisticRegression(T = Float64; inputdim::Integer, nclasses::Integer,
                            initσ::Real = 0.0, pseudocounts::Real = 1.0,
                            hasbias::Bool = true)
    dim = hasbias ? inputdim + 1 : inputdim

    # Prior parameters
    μ₀ = zeros(T, dim)
    Σ₀ = T(1/pseudocounts) * Matrix{T}(I, dim, dim)

    # Initial paremeterization of the posterior
    μ = μ₀ .+ T(initσ) .* randn(T)
    Σ = deepcopy(Σ₀)

    LogisticRegression(μ₀, Σ₀, μ, Σ, hasbias, nclasses = nclasses)
end

getconjugateparams(model::LogisticRegression) = [blr.β for blr in model.stickbreaking]

function reorder(
    model::LogisticRegression{K},
    X::Matrix{T1},
    c::Vector{T2}
) where {K, T1 <: AbstractFloat, T2 <: Integer}

    # Compute the order-dependent terms of the ELBO
    ln1_ν = ones(T1, K-1, length(c))
    for (i, blr) in enumerate(model.stickbreaking)
        X̂ = regressors(blr.hasbias, X)
        μ, Σ = stdparam(blr.β.posterior)
        E_β = μ
        E_ββᵀ = Σ + μ * μ'

        ψ = X̂' * E_β
        ψ² = dropdims(sum(X̂ .* (E_ββᵀ * X̂), dims = 1), dims = 1)

        E_ω, kl_ω = pgaugmentation(ψ²)

        ln1_ν[i, :] = .5 * (-ψ .- E_ω .* ψ² .- kl_ω)
    end

    # Sum over the frame dimension. The result is a KxK-1 matrix.
    M = zeros(K, K-1)
    for i in 1:3 M[i, :] = sum(ln1_ν .* (c .== i)', dims = 2) end

    idxs = sortperm(sum(M, dims = 2)[:, 1], rev = false)
    for (i, j) in enumerate(idxs)
        model.ordering[i] = j
    end
end

function (model::LogisticRegression{K})(
    X::Matrix{T1},
    c::Vector{T2}
) where {K, T1 <: AbstractFloat, T2 <: Integer}

    # Re-map the labels to the "ideal" ordering
    c = [model.ordering[cₙ] for cₙ in c]

    llh = zeros(T1, length(c))
    for i = 1:K-1
        idxs = c .≥ i
        zᵢ= Int.(c[idxs] .== i)
        llh[idxs] .+= model.stickbreaking[i](X[:, idxs], zᵢ)
    end
    return llh
end

function stats_β(
    model::LogisticRegression,
    i::Integer,
    X::Matrix{T1},
    c::Vector{T2}
) where {T1 <: AbstractFloat, T2 <: Real}

    # Re-map the labels to the "ideal" ordering
    c = [model.ordering[cₙ] for cₙ in c]

    idxs = c .≥ i
    zᵢ= Int.(c[idxs] .== i)
    stats_β(model.stickbreaking[i], X[:, idxs], zᵢ)
end

function predict(
    model::LogisticRegression{K},
    X::Matrix{T};
    marginalize::Bool = true
) where {K, T <: AbstractFloat}

    retval = zeros(T, K, size(X, 2))
    for k in 1:K
        residual = dropdims(sum(retval[1:k, :], dims = 1), dims = 1)
        if k < K
            retval[k, :] = predict(model.stickbreaking[k], X,
                                   marginalize = marginalize) .* (1 .- residual)
        else
            retval[k, :] = (1 .- residual)
        end
    end

    # Reverse "ideal" ordering of the target labels
    rordering = Dict(value => key for (key, value) in model.ordering)
    idxs = [rordering[i] for i in 1:K]

    retval[idxs, :]
end

max_stars_repo_name: BUTSpeechFIT/BayesianModels

max_issues_repo_name: BUTSpeechFIT/BayesianModels

max_forks_repo_name: BUTSpeechFIT/BayesianModels


text:
# common utilities

## convenient type alias
#
#  These types signficantly reduces the need of using
#  type parameters in functions (which are often just
#  for the purpose of restricting the arrays to real)
#
# These could be removed when the Base supports
# covariant type notation, i.e. AbstractVector{<:Real}
#

@compat RealArray{T<:Real,N} = AbstractArray{T,N}
@compat RealVector{T<:Real} = AbstractArray{T,1}
@compat RealMatrix{T<:Real} = AbstractArray{T,2}

@compat IntegerArray{T<:Integer,N} = AbstractArray{T,N}
@compat IntegerVector{T<:Integer} = AbstractArray{T,1}
@compat IntegerMatrix{T<:Integer} = AbstractArray{T,2}

@compat const RealFP = Union{Float32, Float64}

## conversion from real to fp types

@compat fptype{T<:Union{Float32,Bool,Int8,UInt8,Int16,UInt16}}(::Type{T}) = Float32
@compat fptype{T<:Union{Float64,Int32,UInt32,Int64,UInt64,Int128,UInt128}}(::Type{T}) = Float64
fptype(::Type{Complex64}) = Complex64
fptype(::Type{Complex128}) = Complex128


text:
function noop_loader(args...; kwargs...)
    @warn "Noop data loader, exiting gracefully..."
    exit()
end


text:
# Take Ni "intermediate" time steps at a time before printing a progress
# statement and updating the time step.
Ni = 1000

function print_progress(simulation)
    model = simulation.model
    i, t = model.clock.iteration, model.clock.time

    progress = 100 * (model.clock.time / end_time)

    umax = maximum(abs, model.velocities.u.data.parent)
    vmax = maximum(abs, model.velocities.v.data.parent)
    wmax = maximum(abs, model.velocities.w.data.parent)
    meanb = mean(view(model.tracers.b.data, 2:Nx-1, 2:Ny-1, 2:Nz-1))

    @printf("[%05.2f%%] i: %d, t: %.2e days, umax: (%6.3e, %6.3e, %6.3e) m/s, CFL: %6.4e, next Δt: %.1e s\n",
            progress, i, t / day, umax, vmax, wmax, cfl(model), Δt_wizard.Δt)
    println(" ")
    @printf("The mean value of b is %.2e", meanb)
    println(" ")
end

max_stars_repo_name: sandreza/Mesoscale

max_issues_repo_name: sandreza/Mesoscale

max_forks_repo_name: sandreza/Mesoscale


text:
using GeoAcceleratedArrays
using Test

@testset "GeoAcceleratedArrays.jl" begin
    # Accelerate 1000 points
    points = rand(1, 1000)
    acc_points = accelerate(points, KDTreeIndex)

    center = [0.0]
    radius = 0.1
    sphere = HyperSphere(center, radius)

    indices = findall(in(sphere), acc_points)
    in_sphere = acc_points[:, indices]
    @test all(in_sphere .<= radius)
end


text:
include("../examples/smiley_examples.jl")

using Test

using NumberIntervals, Turmeric
using .SmileyExample22, .SmileyExample52, .SmileyExample54, .SmileyExample55

const tol = 1e-6

#NOTE: Bisection method performs badly in all examples
for method in (Krawczyk(), Newton())

@testset "$method: $(SmileyExample22.title)" begin
    region = NumberInterval.(SmileyExample22.region)
    roots_found, rest_regions = roots(SmileyExample22.f, region, method, tol)
    @test length(roots_found) == 8
    @test length(rest_regions) == 0
    # no reference data for roots given
end

for example in ifelse(method == Newton(),
    (SmileyExample52, SmileyExample54, ), # SmileyExample55 takes too long with Newton
    (SmileyExample52, SmileyExample54, SmileyExample55),
)
    @testset "$method: $(example.title)" begin
        region = NumberInterval.(example.region)
        roots_found, rest_regions = roots(example.f, region, method, tol)
        @test length(roots_found) == length(example.known_roots)
        @test length(rest_regions) == 0
        for rf in roots_found
            # check there is exactly one known root for each found root
            @test sum(!any(isempty.(rk .∩ rf)) for rk in example.known_roots) == 1
        end
    end
end

end # method loop


text:
# Default values for unspecified objects
baseval{T <: String}(s::Type{T}) = ""
baseval(x::Any) = zero(x)

# iterate over vals array to find the most generic non-NA type
function typeloop(vals)
    toptype = None
    for i = 1:length(vals)
        if vals[i] != :NA
            toptype = promote_type(toptype, typeof(vals[i]))
        end
    end
    if !method_exists(baseval, (toptype, ))
        error("No baseval exists for type: $(toptype)")
    end
    return toptype
end

# Enumerate a BitArray. This is faster than using enumerate()
macro bitenumerate(ba, i, x, code)
    esc(quote
        $i = 1
        for $x in $ba
            $code
            $i += 1
        end
    end)
end


text:
import Base.getindex

export Chromosome, getindex, print_chromosome, random_chromosome

type Chromosome
    params::Parameters
    inputs::Vector{InputNode}
    interiors::Matrix{InteriorNode}
    outputs::Vector{OutputNode}
    active_set::Bool
    number_active_nodes::Integer
end

function Chromosome(p::Parameters)
    inputs = Array(InputNode, p.numinputs)
    interiors = Array(InteriorNode, p.numlevels, p.numperlevel)
    outputs = Array(OutputNode, p.numoutputs)
    fitness = 0.0
    active_set = false
    number_active_nodes = 0

    return Chromosome(p, inputs, interiors, outputs, active_set, number_active_nodes)
end

function getindex(c::Chromosome, level::Integer, index::Integer)
    if level == 0
        return c.inputs[index]
    end

    if level > c.params.numlevels
        return c.outputs[index]
    end

    return c.interiors[level, index]
end

# Prints the chromosome in a compact text format (on one line).  Active notes
# are indicated with a "+" and inactive notes with a "*".  If the chromosome
# has not been evaluated, the active field is set to "?" since the status of
# the nodes (with the exception of output nodes) is unknown.  If active_only is
# true, then only the active notes are shown.
function print_chromosome(c::Chromosome, active_only::Bool)
    active_set = c.active_set

    # Input nodes
    for i = 1:length(c.inputs)
        if c.inputs[i].active || !active_only
            active = c.inputs[i].active ? "+" : (active_set ? "*" : "?")
            print("[in", i, active, "] ")
        end
    end

    # Interior nodes
    for i = 1:c.params.numlevels
        for j = 1:c.params.numperlevel
            if c.interiors[i,j].active || !active_only
                active = c.interiors[i,j].active ? "+" : (active_set ? "*" : "?")
                print("[")
                for k = 1:length(c.interiors[i,j].inputs)
                    if c.interiors[i,j].inputs[k][1] == 0
                        print("(in", c.interiors[i,j].inputs[k][2], ")")
                    else
                        print(c.interiors[i,j].inputs[k])
                    end
                end
                print("\"", c.interiors[i, j].func.name, "\"", (i, j), active, "] ")
            end
        end
    end

    # Output nodes
    for i = 1:length(c.outputs)
        active = c.outputs[i].active ? "+" : (active_set ? "*" : "?")
        print("[", c.outputs[i].input, "out", i, active, "] ")
    end
    println()
    return
end

function print_chromosome(c::Chromosome)
    print_chromosome(c, false)
end

function first_in_level(p::Parameters, level::Integer)
    if level == 0
        first = 1
    else
        first = p.numinputs + 1 + (level - 1) * p.numperlevel
    end
    return first
end

function last_in_level(p::Parameters, level::Integer)
    if level == 0
        last = p.numinputs
    else
        last = p.numinputs + level * p.numperlevel
    end
    return last
end

function random_node_position(p::Parameters, minlevel::Integer, maxlevel::Integer)
    first = first_in_level(p, minlevel)
    last = last_in_level(p, maxlevel)

    i = rand(first:last)

    if i <= p.numinputs
        level = 0
        index = i
    else
        level = div(i - p.numinputs - 1, p.numperlevel) + 1
        index = rem(i - p.numinputs - 1, p.numperlevel) + 1
    end

    return (level, index)
end

function random_chromosome(p::Parameters )
    return random_chromosome(p, p.funcs)
end

function random_chromosome(p::Parameters, funcs::Vector{Func})
    c = Chromosome(p)

    for index = 1:length(c.inputs)
        c.inputs[index] = InputNode(index)
    end

    for level = 1:p.numlevels
        minlevel = max(level - p.numlevelsback, 0)
        maxlevel = level - 1
        for index = 1:p.numperlevel
            func = funcs[rand(1:end)]
            inputs = Array(NodePosition, func.arity)
            for i = 1:func.arity
                (input_level, input_index) = random_node_position(p, minlevel, maxlevel)
                inputs[i] = (input_level, input_index)
            end
            c.interiors[level, index] = InteriorNode(func, inputs)
        end
    end

    minlevel = p.numlevels - p.numlevelsback + 1
    maxlevel = p.numlevels
    for i = 1:length(c.outputs)
        (level, index) = random_node_position(p, minlevel, maxlevel)
        c.outputs[i] = OutputNode((level, index))
    end
    return c
end


text:
module DOT

using GraphIO.ParserCombinator.Parsers
using Graphs
using Graphs: AbstractGraphFormat

import Graphs: loadgraph, loadgraphs, savegraph

export DOTFormat

struct DOTFormat <: AbstractGraphFormat end

function savedot(io::IO, g::Graphs.AbstractGraph, gname::String = "")
    isdir = Graphs.is_directed(g)
    println(io,(isdir ? "digraph " : "graph ") * gname * " {")
    for i in Graphs.vertices(g)
         println(io,"\t" * string(i))
    end
    if isdir
        for u in Graphs.vertices(g)
            out_nbrs = Graphs.outneighbors(g, u)
            length(out_nbrs) == 0 && continue
            println(io, "\t" * string(u) * " -> {" * join(out_nbrs,',') * "}")
        end
    else
        for e in Graphs.edges(g)
            source = string(Graphs.src(e))
            dest = string(Graphs.dst(e))
            println(io, "\t" * source * " -- " * dest)
        end
    end
    println(io,"}")
    return 1
end

function savedot_mult(io::IO, graphs::Dict)
    ng = 0
    for (gname, g) in graphs
        ng += savedot(io, g, gname)
    end
    return ng
end

function _dot_read_one_graph(pg::Parsers.DOT.Graph)
    isdir = pg.directed
    nvg = length(Parsers.DOT.nodes(pg))
    nodedict = Dict(zip(collect(Parsers.DOT.nodes(pg)), 1:nvg))
    if isdir
        g = Graphs.DiGraph(nvg)
    else
        g = Graphs.Graph(nvg)
    end
    for es in Parsers.DOT.edges(pg)
        s = nodedict[es[1]]
        d = nodedict[es[2]]
        add_edge!(g, s, d)
    end
    return g
end

_name(pg::Parsers.DOT.Graph) =
    pg.id !== nothing ? pg.id.id :
    Parsers.DOT.StringID(pg.directed ? "digraph" : "graph")

function loaddot(io::IO, gname::String)
    p = Parsers.DOT.parse_dot(read(io, String))
    for pg in p
        _name(pg) == gname && return _dot_read_one_graph(pg)
    end
    error("Graph $gname not found")
end

function loaddot_mult(io::IO)
    p = Parsers.DOT.parse_dot(read(io, String))
    graphs = Dict{String,AbstractGraph}()

    for pg in p
        graphs[_name(pg)] = _dot_read_one_graph(pg)
    end
    return graphs
end

loadgraph(io::IO, gname::String, ::DOTFormat) = loaddot(io, gname)
loadgraphs(io::IO, ::DOTFormat) = loaddot_mult(io)
savegraph(io::IO, g::AbstractGraph, gname::String, ::DOTFormat) = savedot(io, g, gname)
savegraph(io::IO, d::Dict, ::DOTFormat) = savedot_mult(io, d)

end #module


text:
module HonlyHelperTurret

using ..Ahorn, Maple

@mapdef Entity "HonlyHelper/Turret" Turret(x::Integer, y::Integer, turretID::String = "TurretID", aimTime::Number = 2.0, cooldownTime::Number = 2.0, randomCooldownMultiplier::Number = 0.0, desiredBulletSpeed::Number = 2000.0, accelerationMultiplier::Number = 1.0)

const placements = Ahorn.PlacementDict(
    "Turret (HonlyHelper, WIP)" => Ahorn.EntityPlacement(
        Turret
    )
)

sprite = "objects/HonlyHelper/Turret/idle00.png"

function Ahorn.selection(entity::Turret)
    x, y = Ahorn.position(entity)
    return Ahorn.Rectangle(x - 16, y - 16, 24, 24)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Turret, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, -4, -4)

end

max_stars_repo_name: HollyMagala/HonlyHelper

max_issues_repo_name: HollyMagala/HonlyHelper

max_forks_repo_name: HollyMagala/HonlyHelper


text:
using Documenter, Faker

makedocs(format = Documenter.HTML(),
         modules = [Faker],
         sitename = "Faker.jl"
         )
deploydocs(
	repo = "github.com/codeneomatrix/Faker.jl.git"
	#devbranch = "master",
    #devurl = "dev",
    #versions = ["stable" => "v^", "v#.#", "devurl" => devurl]
    )


text:
#=
rect.jl
=#

#using ImagePhantoms #: Object, Object2d

export Rect
export Square
export phantom, radon, spectrum

"""
    Rect <: AbstractShape2
"""
struct Rect <: AbstractShape2 end

# constructors

"""
    Rect(cx, cy, wx=1, wy=wx, ϕ=0, value::Number=1)
    Rect(center::NTuple{2,RealU}, width::NTuple{2,RealU}=(1,1), ϕ::RealU=0, v=1)
    Rect([6-vector])
Construct `Rect` object from parameters;
here `width` is the full-width.
"""
function Rect(
    cx::RealU,
    cy::RealU,
    wx::RealU = oneunit(cx),
    wy::RealU = wx,
    ϕ::RealU = 0,
    value::Number = 1,
)
    (cx, cy, wx, wy) = promote(cx, cy, wx, wy)
    Object(Rect(), (cx,cy), (wx,wy), ϕ, value)
end

function Rect(
    center::NTuple{2,RealU},
    width::NTuple{2,RealU} = (1,1) .* oneunit(center[1]),
    ϕ::RealU = 0,
    value::Number = 1,
)
    Rect(center..., width..., ϕ, value)
end

function Rect(v::AbstractVector{<:Number})
    length(v) == 6 || throw(ArgumentError("$v wrong length"))
    Rect(v...)
end

# squares as a special case

"""
    Square(x,y,w,v=1) (square of width `w` centered at `(x,y)`)
    Square((x,y), w=1, v=1) ditto
    Square([4-vector]) ditto
    Square(w, v=1) centered at origin
Construct `Square` objects as special cases of `Rect` objects.
"""
Square(cx::RealU, cy::RealU, w::RealU, v::Number = 1) =
    Rect(cx, cy, w, w, 0, v)
Square(center::NTuple{2,RealU}, w::RealU = oneunit(center[1]), v::Number = 1) =
    Square(center..., w, v)
Square(w::RealU, v::Number = 1) = Square((zero(w),zero(w)), w, v)

function Square(v::AbstractVector{<:Number})
    length(v) == 4 || throw(ArgumentError("$v wrong length"))
    Square(v...)
end

# helper

"""
    trapezoid(t::RealU, t1, t2, t3, t4)
Unit-height trapezoid with breakpoints `t1`, `t2`, `t3`, `t4`.
"""
function trapezoid(t::RealU, t1::RealU, t2::RealU, t3::RealU, t4::RealU)
    (t, t1, t2, t3, t4) = promote(t, t1, t2, t3, t4)
    T = eltype(t)
    if t1 < t < t2
        return (t - t1)/(t2 - t1)
    elseif t2 <= t <= t3
        return one(T)
    elseif t3 < t < t4
        return (t4 - t)/(t4 - t3)
    end
    return 0*one(T) # zero(T) has units for Unitful :(
end

# methods

"""
    phantom(ob::Object2d{Rect})
Returns function of `(x,y)` for making image.
"""
phantom(ob::Object2d{Rect}) = (x,y) ->
    ob.value * (maximum(abs, coords(ob, x, y)) ≤ 0.5)

"""
    radon_rect(r, ϕ, cx, cy, wx, wy, θ)
Radon transform at `(r,ϕ)` of rectangle (involves a trapezoid).
"""
function radon_rect(r, ϕ, cx, cy, wx, wy, θ)
    (sinϕ, cosϕ) = sincos(ϕ)
    r -= cx * cosϕ + cy * sinϕ # Radon translation property
    (sinϕ, cosϕ) = sincos(ϕ - θ) # Radon rotation property
    xmax = wx * abs(cosϕ)
    ymax = wy * abs(sinϕ)
    lmax = wx * wy / max(xmax, ymax)
    dmax = (xmax + ymax) / 2
    dbreak = abs(xmax - ymax) / 2
    return lmax * trapezoid(r, -dmax, -dbreak, dbreak, dmax)
end

"""
    radon(ob::Object2d{Rect})
Returns function of `(r,ϕ)` for making a sinogram.
"""
radon(ob::Object2d{Rect}) = (r,ϕ) -> ob.value *
    radon_rect(r, ϕ, ob.center..., ob.width..., ob.angle[1])

function spectrum_rect(fx, fy, cx, cy, wx, wy, θ)
    (kx, ky) = rotate2d(fx, fy, θ) # rect is rotated first, then translated
    return wx * sinc(kx * wx) * exp(-2im*π*fx*cx) *
           wy * sinc(ky * wy) * exp(-2im*π*fy*cy)
end

"""
    spectrum(ob::Object2d{Rect})
Returns function of ``(f_x,f_y)`` for the spectrum (2D Fourier transform).
"""
spectrum(ob::Object2d{Rect}) = (fx,fy) -> ob.value *
    spectrum_rect(fx, fy, ob.center..., ob.width..., ob.angle[1])


text:
export V3DExport, V3DEvents, V3DExportSource, V3DEventsSource

export writeeventsfile

struct V3DExport; end
struct V3DEvents; end

const V3DExportSource = Source{V3DExport}
const V3DEventsSource = Source{V3DEvents}

function DatasetManager.readsource(s::V3DExportSource; kwargs...)
    return matread(sourcepath(s))
end

function DatasetManager.readsegment(seg::Segment{V3DExportSource};
    events::Union{Nothing,Vector}=nothing,
    series::Union{Nothing,Vector}=nothing
)
    file = matopen(sourcepath(seg.source))
    fs = only(read(file, "FRAME_RATE"))
    start, finish = seg.start, seg.finish

    # TODO: Check that seg.finish is <= the (time) length of the source

    outevents = Dict{String,Vector{Float64}}()
    if events !== nothing
        for event in events
            name = String(event)
            if exists(file, name)
                tmp = only(read(file, name))
                if tmp isa AbstractArray
                    outevents[name] = vec(tmp)
                else
                    outevents[name] = [tmp]
                end
                startidx = findfirst(>=(start), outevents[name])
                if startidx == 0
                    # @warn "no $event events during given start and end times"
                    delete!(outevents, name)
                    break
                end
                # Might throw, FIXME, see above
                finidx = isnothing(finish) ? lastindex(outevents[name]) : findlast(<=(finish), outevents[name])

                # Shift events to be index accurate for data subsection
                outevents[name] = outevents[name][startidx:finidx] .- start .+ (1/fs)
                sort!(outevents[name])
            else
                # @warn "Requested event $event does not exist in source data"
            end
        end
    end

    outseries = Dict{String,Matrix{Float64}}()
    if series !== nothing
        for t in series
            name = String(t)
            if exists(file, name)
                tmp = only(read(file, string(t)))
                len = size(tmp, 1)
                startidx = round(Int, start*fs)
                startidx += (startidx === 0) ? 1 : 0
                finidx = isnothing(finish) ? len : min(len, round(Int, finish*fs))
                outseries[name] = tmp[startidx:finidx, :]
            else
                @warn "Requested time series $t does not exist in source data"
            end
        end
    end

    close(file)

    return (fs, outevents, outseries)
end

function any_or_bust(vec)
    return (vec isa Vector{Union{}}) ? [] : vec
end

function DatasetManager.readsource(s::V3DEventsSource; kwargs...)
    events = CSV.File(sourcepath(s); header=2, skipto=6, drop=[1], kwargs...)

    # @show identity(collect(identity(collect(Tables.getcolumn(events, :RFO)))))
    return Dict(string(name) => sort(any_or_bust(collect(skipmissing(Tables.getcolumn(events, name)))))
        for name in Tables.columnnames(events))
end

function deletefirst!(a::Vector, n::Integer)
    zero(n) < n < length(a) && Base._deletebeg!(a, n)
    return a
end

function DatasetManager.readsegment(seg::Segment{V3DEventsSource}; kwargs...)
    src = readsource(seg.source; kwargs...)

    if !isnothing(seg.finish)
        foreach(src) do (_, v)
            l = searchsortedlast(v, seg.finish)
            resize!(v, l)
        end
    end

    if !isnothing(seg.start)
        foreach(src) do (_, v)
            f = searchsortedfirst(v, seg.start)
            if f ≤ length(v)
                deletefirst!(v, f-1)
            else
                empty!(v)
            end
            v .-= seg.start
        end
    end

    return src
end

function DatasetManager.generatesource(trial, src::V3DEventsSource, deps; genfunc, kwargs...)
    mkpath(dirname(sourcepath(src)))

    events = genfunc(trial; kwargs...)
    writeeventsfile(sourcepath(src), events)

    return src
end

function writeeventsfile(
    fn::String, events)
    writeeventsfile(fn, pairs(events)...)
end

function writeeventsfile(
    fn::String, events::Vararg{Union{Tuple{K,<:AbstractVector},Pair{K,<:AbstractVector}}}
) where K <: Union{Symbol, String}
    header = fill("", (5,1))
    header[5,1] = "ITEM"
    numevents = maximum(broadcast(x -> size(x[2],1), events))
    eventdata = string.(collect(1:numevents))

    file = splitext(basename(fn))[1]*".c3d"

    for event in events
        header = [ header [ file,
                            string(event[1]),
                            "EVENT_LABEL",
                            "ORIGINAL",
                            "X" ] ]
        tmp = fill("", (numevents))
        tmp[axes(event[2], 1)] .= string.(event[2])
        eventdata = [ eventdata tmp ]
    end

    mkpath(dirname(fn))
    CSV.write(fn, Tables.table([header; eventdata]); delim='\t', writeheader=false)

    return nothing
end


text:
"""
    `function readsong(fname::String, starttime::Real, endtime::Real, dir::String=pwd())`
Given a file name, a working directory, a start time and an end time, read a segment of a song from a file

(note: This is primarily used as a helper function to the `readsongs_*` family of functions)
"""
function readsong(fname::String, starttime::Real, endtime::Real, dir::String=pwd())
    return wavread(dir*"/"*fname, Int(round(samplerate*starttime)):Int(round(samplerate*endtime)))[1][:, 1]
end

"""
    `function readsongs(calls::DataFrame; dir::String=pwd(), filenamecol::Symbol=:file, starttimecol::Symbol=:start_time, endtimecol::Symbol=:end_time)`

Given a DataFrame `calls`, read out each call in the DataFrame into an array of Float64 and return an array of all the arrays.

Parameters
----------
 - `calls`: A DataFrame, containing the call definitions. It must have the three following columns
 - `dir`: The directory where the files are located. Defaults to the current working directory
 - `filenamecol`: The name of the column in the DataFrame containing file names
 - `starttimecol`: The name of the column containing start times of calls
 - `endtimecol`: The name of the column containing end times of calls
Example:
--------
Given the following DataFrame:

        | Row |   Name   | start | end |
    x = |_____|__________|_______|_____|
        |  1  | rec1.wav |  0.02 | 1.6 |
        |  2  | rec1.wav |  1.9  | 3.7 |
        |  3  | rec2.wav |  3.5  | 5.1 |

The call `calls = readsongs(x, filenamecol=:Name, starttimecol=:start, endtimecol=:end)` will then
return an array containing the sound data in the file rec1.wav from .02 to 1.6 seconds, the sound data
in rec1.wav from 1.9 to 3.7 seconds, and the sound data in rec2.wav from 3.5 to 5.1 seconds.
"""

function readsongs(calls::DataFrame; dir::String=pwd(), filenamecol::Symbol=:file, starttimecol::Symbol=:start_time, endtimecol::Symbol=:end_time)
    #Check for sanity of calls
    if !haskey(calls, filenamecol)
	ArgumentError("Argument `calls` must have a column named '$filenamecol'")
    end
    if !haskey(calls, starttimecol)
	ArgumentError("Argument `calls` must have a column named '$starttimecol'")
    end
    if !haskey(calls, endtimecol)
	ArgumentError("Argument `calls` must have a column named '$endtimecol'")
    end

    nSongs = size(calls, 1)
    songs = Array{Array{Float64, 1}, 1}(nSongs)

@showprogress    for i in 1:nSongs
	songs[i] = readsong(calls[filenamecol][i], calls[starttimecol][i], calls[endtimecol][i], dir)
    end

    return songs
end

"""
    `function readsongs_dict(calls::DataFrame; dir::String=pwd(), filenamecol::Symbol=:file, starttimecol::Symbol=:start_time, endtimecol::Symbol=:end_time)`
Like `readsongs` but instead of an array, returns a dict with keys of the form 'filename_starttime-endtime'
"""
function readsongs_dict(calls::DataFrame; dir::String=pwd(), filenamecol::Symbol=:file, starttimecol::Symbol=:start_time, endtimecol::Symbol=:end_time)
    #Check for sanity of calls
    if !haskey(calls, filenamecol)
	ArgumentError("Argument `calls` must have a column named '$filenamecol'")
    end
    if !haskey(calls, starttimecol)
	ArgumentError("Argument `calls` must have a column named '$starttimecol'")
    end
    if !haskey(calls, endtimecol)
	ArgumentError("Argument `calls` must have a column named '$endtimecol'")
    end

    nSongs = size(calls, 1)
    songs = Dict{String, Array{Float64, 1}}()

@showprogress    for i in 1:nSongs
	name = calls[filenamecol][i]*"_"*"$(calls[starttimecol][i])"*"-"*"$(calls[endtimecol][i])"
	songs[name] = readsong(calls[filenamecol][i], calls[starttimecol][i], calls[endtimecol][i], dir)
    end

    return songs
end

"""
    `function readsongs_string(calls::DataFrame; dir::String=pwd(), filenamecol::Symbol=:file, starttimecol::Symbol=:start_time, endtimecol::Symbol=:end_time)`
Like readsongs_dict, but only returns an array containing what would be the keys of the dictionary
"""
function readsongs_string(calls::DataFrame; dir::String=pwd(), filenamecol::Symbol=:file, starttimecol::Symbol=:start_time, endtimecol::Symbol=:end_time)
    #Check for sanity of calls
    if !haskey(calls, filenamecol)
	ArgumentError("Argument `calls` must have a column named '$filenamecol'")
    end
    if !haskey(calls, starttimecol)
	ArgumentError("Argument `calls` must have a column named '$starttimecol'")
    end
    if !haskey(calls, endtimecol)
	ArgumentError("Argument `calls` must have a column named '$endtimecol'")
    end

    nSongs = size(calls, 1)
    songs = Array{String, 1}(nSongs)

    for i in 1:nSongs
	name = calls[filenamecol][i]*"_"*"$(calls[starttimecol][i])"*"-"*"$(calls[endtimecol][i])"
	songs[i] = name
    end

    return songs
end

export readsong
export readsongs
export readsongs_dict
export readsongs_string


text:
module ReadFeatures

using LightGraphs: AbstractGraph, nv, ne,
                   Δ, δ, density, degree,
                   diameter, radius,
                   global_clustering_coefficient,
                   k_core
using Statistics
using StatsBase
using DataStructures

include("../utils/misc.jl")
include("../utils/graph.jl")
include("readfeatures.jl")

end # module

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
@testset "control comparison" begin
    original_network = "data/nordic_fault/base/fault_4032_4044.m"
    # We choose to put the new network data in the same folder because
    # we also want to use some data that is already there. Each upgrade
    # gets its own folder.
    new_data_directory = "results/inputs/nordic_fault/"
    standard_frequency = 50.0
    lfac_branch_upgrades = [27,33]

    VariableFrequencyOPF.upgrade_branches(
        original_network,
        new_data_directory,
        standard_frequency,
        indices=lfac_branch_upgrades
    )

    # Minimize generation in areas 2 and 3 with plots including
    # the power flow between areas 1 and 2
    objective = "areagen"
    gen_areas = [2,3]
    area_transfer = [1,2]
    # We want to plot the results of all branch upgrades,
    # so we set `enum_branch` to true
    enum_branches = true

    VariableFrequencyOPF.control_comparison(
        new_data_directory,
        "results/nordic_fault/",
        objective,
        gen_areas=gen_areas,
        area_transfer=area_transfer,
        enum_branches=true,
        results_folders=[],
        start_vals=Dict("sn"=>Dict(2=>50))
    )

    results = JSON.parsefile("results/nordic_fault/results_dict_allplots.json")
    results_dict = results[1]

    subnet1 = findall(results_dict["subnet"]["br27"].==1)[1]
    subnet2 = 3-subnet1

    @test results_dict["status"]["br27"] == "LOCALLY_SOLVED"
    @test results_dict["status"]["br33"] == "LOCALLY_SOLVED"
    @test isapprox(results_dict["cost"]["br27"][subnet1], 4.4015e6; rtol = 1e-2)
    @test isapprox(results_dict["cost"]["br27"][subnet2], 0.0; rtol = 1e-2)
    @test isapprox(results_dict["frequency (Hz)"]["br27"][subnet1], 50.0; atol = 1e0)
    # @test isapprox(results_dict["frequency (Hz)"]["br27"][subnet2], 20.92; atol = 1e0)
    @test isapprox(results_dict["cost"]["br33"][subnet1], 4.3645e6; rtol = 1e-2)
    @test isapprox(results_dict["cost"]["br33"][subnet2], 0.0; rtol = 1e-2)
    @test isapprox(results_dict["frequency (Hz)"]["br33"][subnet1], 50.0; atol = 1e0)
    # @test isapprox(results_dict["frequency (Hz)"]["br33"][subnet2], 22.43; atol = 1e0)

    @test isfile("results/nordic_fault/plots/total/total generation.pdf")
    @test isfile("results/nordic_fault/plots/subnet1/subnet 1 line loss.pdf")
    @test isfile("results/nordic_fault/plots/subnet2/subnet 2 frequency (Hz).pdf")

    try
        rm("results",force=true,recursive=true)
    catch IOError
        println("Unable to delete the results files generated by this test.")
    end
end

@testset "hvdc comparison" begin
    original_network = "data/nordic_fault/base/fault_4032_4044.m"
    # We choose to put the new network data in the same folder because
    # we also want to use some data that is already there. Each upgrade
    # gets its own folder.
    new_data_directory = "results/inputs/nordic_fault/"
    standard_frequency = 50.0
    lfac_branch_upgrades = [27,33]

    VariableFrequencyOPF.upgrade_branches(
        original_network,
        new_data_directory,
        standard_frequency,
        indices=lfac_branch_upgrades,
        noloss=true
    )

    # Minimize generation in areas 2 and 3 with plots including
    # the power flow between areas 1 and 2
    objective = "areagen"
    gen_areas = [2,3]
    area_transfer = [1,2]
    # We want to plot the results of all branch upgrades,
    # so we set `enum_branch` to true
    enum_branches = true

    VariableFrequencyOPF.hvdc_comparison(
        new_data_directory,
        "results/nordic_fault/",
        objective,
        gen_areas=gen_areas,
        area_transfer=area_transfer,
        enum_branches=true,
        series_labels=[
           "LFAC",
           "k_cond=1,k_ins=sqrt(2)",
           "k_cond=2_3,k_ins=sqrt(2)",
           "k_cond=1,k_ins=1",
           "k_cond=2_3,k_ins=1"
           ],
        k_cond=[1,2/3,1,2/3],
        k_ins=[sqrt(2),sqrt(2),1,1],
        dc_subnet=-1,
    )

    results = JSON.parsefile("results/nordic_fault/results_dict_allplots.json")
    results_dict = results[2]

    subnet1 = findall(results_dict["subnet"]["br27"].==1)[1]
    subnet2 = 3-subnet1

    @test results_dict["status"]["br27"] == "LOCALLY_SOLVED"
    @test results_dict["status"]["br33"] == "LOCALLY_SOLVED"
    @test isapprox(results_dict["cost"]["br27"][subnet1], 4.3564e6; rtol = 1e-2)
    @test isapprox(results_dict["cost"]["br27"][subnet2], 0.0; rtol = 1e-2)
    @test isapprox(results_dict["frequency (Hz)"]["br27"][subnet1], 50.0; atol = 1e0)
    @test isapprox(results_dict["frequency (Hz)"]["br27"][subnet2], 0.0; atol = 1e0)
    @test isapprox(results_dict["cost"]["br33"][subnet1], 4.3105e6; rtol = 1e-2)
    @test isapprox(results_dict["cost"]["br33"][subnet2], 0.0; rtol = 1e-2)
    @test isapprox(results_dict["frequency (Hz)"]["br33"][subnet1], 50.0; atol = 1e0)
    @test isapprox(results_dict["frequency (Hz)"]["br33"][subnet2], 0.0; atol = 1e0)

    @test isfile("results/nordic_fault/plots/total/total generation.pdf")
    @test isfile("results/nordic_fault/plots/subnet1/subnet 1 line loss.pdf")
    @test isfile("results/nordic_fault/plots/subnet2/subnet 2 frequency (Hz).pdf")

    try
        rm("results",force=true,recursive=true)
    catch IOError
        println("Unable to delete the results files generated by this test.")
    end
end


text:
export setupsymmetry
export WithoutSymmetry, PlanarSymmetry, PlanarAzimuthalSymmetry, AzimuthalSymmetry, RadialSymmetry

"""
    AbstractSetupSymmetry

An abstract types which dictates the symmetry of the setup. That is, the symmetry shared between the incident wave and the shape of the material.
"""
abstract type AbstractSetupSymmetry{Dim} end

"""
    WithoutSymmetry

A type used to describe materials and incident waves which share no common symmetry. This will lead to the most general regular wave expansion for the eignvectors.
"""
struct WithoutSymmetry{Dim} <: AbstractSetupSymmetry{Dim} end

"""
An incident plane-wave and halfspace material will result in all fields being plane-waves.
"""
abstract type AbstractPlanarSymmetry{Dim} <: AbstractSetupSymmetry{Dim} end

"""
    PlanarSymmetry

A type used to describe materials and incident waves that both share a planar symmetry.
"""
struct PlanarSymmetry{Dim} <: AbstractPlanarSymmetry{Dim} end

"""
For spatial dimension > 2, we can consider problems that have azimuthal symmetry. For example, a plane-wave incident on a sphere.
"""
abstract type AbstractAzimuthalSymmetry{Dim} <: AbstractSetupSymmetry{Dim} end

"""
    AzimuthalSymmetry

A type used to describe materials and incident waves in 3 spatial dimensions that share a symmetry about the azimuth.
"""
struct AzimuthalSymmetry{Dim} <: AbstractAzimuthalSymmetry{Dim} end
AzimuthalSymmetry() = AzimuthalSymmetry{3}()

"""
    RadialSymmetry

A type used to describe materials and incident waves in that are both radially symmetric. That is, the material is a sphere, and the incident wave is radially symmetric around the center of this spherical material.
"""
struct RadialSymmetry{Dim} <: AbstractAzimuthalSymmetry{Dim} end

"""
    PlanarAzimuthalSymmetry

A type used to describe a materail and incident wave which have both [`PlanarSymmetry`](@ref) and [`AzimuthalSymmetry`](@ref).
"""
struct PlanarAzimuthalSymmetry{Dim} <: AbstractPlanarSymmetry{Dim} end
PlanarAzimuthalSymmetry() = PlanarAzimuthalSymmetry{3}()

# """Extract the dimension of the space that this physical property lives in"""
# dim(p::AbstractSetupSymmetry{Dim}) where {Dim} = Dim

"""
    Specie

Represents a set of particles which are all the same. The type of particle is given by `Specie.particle` and the volume fraction this specie occupies is given by `Specie.volume_fraction`.

We can use `Specie.numberofparticles` to specify the number of particles, otherwise for an infinite `Specie.numberofparticles = Inf`.

The minimum distance between any two particles will equal `outer_radius(Specie) * Specie.exclusion_distance`.
"""
struct Specie{T<:AbstractFloat,Dim,P<:AbstractParticle{T,Dim}}
    particle::P
    volume_fraction::T
    exclusion_distance::T
end

# Convenience constructor which does not require explicit types/parameters
function Specie(p::AbstractParticle{T,Dim}; number_density::T = 0.0, volume_fraction::T = number_density * volume(p), exclusion_distance::T = 1.005) where {Dim,T<:AbstractFloat}
    Specie{T,Dim,typeof(p)}(p,volume_fraction,exclusion_distance)
end

function Specie(medium::P,s::S; kws...) where {Dim,T,P<:PhysicalMedium{T,Dim},S<:Shape{T,Dim}}
    Specie(Particle(medium, s); kws...)
end

function Specie(medium::P, radius::T; kws...) where {T,P<:PhysicalMedium{T}}
    Specie(Particle(medium, radius); kws...)
end

# Shorthand for all Vectors of species
Species{T<:AbstractFloat,Dim,P} = Vector{S} where S<:Specie{T,Dim,P}

"Returns the volume fraction of the specie."
volume_fraction(s::Specie) = s.volume_fraction
volume_fraction(ss::Species) = sum(volume_fraction.(ss))

import MultipleScattering.volume
volume(s::Specie) = volume(s.particle)

import MultipleScattering.outer_radius

"""
    number_density(s::Specie)

Gives the number of particles per unit volume. Note this is given exactly by `N / V` where `V` is the volume of the region containing the origins of all particles. For consistency, [`volume_fraction`](@ref) is given by `N * volume(s) / V`.
"""
number_density(s::Specie) = s.volume_fraction / volume(s)
# number_density(s::Specie{T,2}) where {T} = s.volume_fraction / (outer_radius(s.particle)^2 * pi)
# number_density(s::Specie{T,3}) where {T} = s.volume_fraction / (T(4/3) * outer_radius(s.particle)^3 * pi)
number_density(ss::Species) = sum(number_density.(ss))

outer_radius(s::Specie) = outer_radius(s.particle)

import MultipleScattering.get_t_matrices
import MultipleScattering.t_matrix

get_t_matrices(medium::PhysicalMedium, species::Vector{S}, ω::AbstractFloat, Nh::Integer) where S<:Specie = get_t_matrices(medium, [s.particle for s in species], ω, Nh)

t_matrix(s::Specie, medium::PhysicalMedium, ω::AbstractFloat, order::Integer) = t_matrix(s.particle, medium, ω, order)

"""
    Material(region::Shape, species::Species [, numberofparticles = Inf])

Creates a material filled with [`Specie`](@ref)'s inside `region`.
"""
struct Material{Dim,S<:Shape,Sps<:Species}
    shape::S
    species::Sps
    numberofparticles::Number
    # Enforce that the Dims and Types are all the same
    function Material{Dim,S,Sps}(shape::S,species::Sps,numberofparticles::Number = Inf) where {T,Dim,S<:Shape{T,Dim},Sps<:Species{T,Dim}}
        new{Dim,S,Sps}(shape,species,numberofparticles)
    end
end

# Convenience constructor which does not require explicit types/parameters
function Material(shape::S,species::Sps) where {T,Dim,S<:Shape{T,Dim},Sps<:Species{T,Dim}}

    V = volume(shape)
    numberofparticles = round(sum(
        volume_fraction(s) * V / volume(s)
    for s in species))

    Material{Dim,S,Sps}(shape,species,numberofparticles)
end

function Material(shape::S,specie::Sp) where {T,Dim,S<:Shape{T,Dim},Sp<:Specie{T,Dim}}
    Material(shape,[specie])
end

import MultipleScattering.PhysicalMedium

PhysicalMedium(s::Specie) = typeof(s.particle.medium)
PhysicalMedium(m::Material) = PhysicalMedium(m.species[1])

"""
    setupsymmetry(source::AbstractSource, material::Material)

Returns the shared symmetries between the `source` and `materail`.
"""
setupsymmetry(source::AbstractSource, material::Material{Dim}) where Dim = WithoutSymmetry{Dim}()

function setupsymmetry(source::PlaneSource{T,3,1}, material::Material{3,Sphere{T,3}};
        basis_order::Int = 4, ω::T = 0.9) where T

    ls, ms = spherical_harmonics_indices(basis_order)

    gs = regular_spherical_coefficients(source)(basis_order,origin(material.shape),ω)
    azimuthal_symmetry = norm((ms[n] != 0) ? gs[n] : zero(T) for n in eachindex(gs)) ≈ zero(T)

    return if azimuthal_symmetry
        AzimuthalSymmetry{3}()
    else
        WithoutSymmetry{3}()
    end
end

function setupsymmetry(psource::PlaneSource{T,Dim}, material::Material{Dim,S}) where {T<:AbstractFloat, Dim, S<:Union{Halfspace{T,Dim},Plate{T,Dim}}}

    hv = material.shape.normal
    kv = psource.direction

    if abs(dot(hv, kv)^2) ≈ abs(dot(hv, hv) * dot(kv, kv))
        # for direct incidence
        return PlanarAzimuthalSymmetry{Dim}()
    else
        return PlanarSymmetry{Dim}()
    end
end


text:
####################################
# Lighting

@field Shininess = 16f0

@composed type Shading
    Ambient::Vec3f0
    Specular::Vec3f0
    Shininess::Float32
end

@composed type SceneUniforms
    Projection::Mat4f0
    View::Mat4f0
    ProjectionView::Mat4f0
    Resolution::Vec2f0
    EyePosition::Vec3f0
end

function blinnphong{NV, T}(V::Vec{NV, T}, N, L, color, shading, light)
    lambertian = max(dot(L, N), 0f0)
    half_direction = normalize(L .+ V)
    specular_angle = max(dot(half_direction, N), 0.0)
    specular = specular_angle ^ shading.shininess
    surf_color = (lambertian * color) .+ (specular * shading.specular)
    return light.ambient .+ surf_color
end

flat_shading(V, N, L, color, shading, light) = color

immutable Vert2Frag
    color::Vec4f0
    position::Vec3f0
    normal::Vec3f0
    uv::Vec2f0
    lightdir::Vec3f0
end

immutable Vertex2Geom
    uvrect::Vec4f0
    color::Vec4f0
    rect::Vec4f0
end

abstract type AbstractVertex end

immutable VertexN <: AbstractVertex
    position::Vec3f0
    normal::Vec3f0
end
immutable VertexNUV <: AbstractVertex
    position::Vec3f0
    normal::Vec3f0
    uv::Vec2f0
end

abstract type AbstractLineVertex <: AbstractVertex end

@composed immutable LineVertex <: AbstractLineVertex
    Position
    Thickness::Float32
    Color::Vec4f0
end

abstract type Particle <: AbstractVertex end

immutable VertexCS{N, T} <: Particle
    position::Vec{N, T}
    color::Vec4f0
    scale::Vec2f0
end
immutable Sprite{N, T} <: Particle
    position::Point{N, T}
    offset::Vec{2, T}
    scale::Vec{2, T}
    uv::Vec{4, T}
    color::Vec{4, T}
end

# since we don't have traits in that sense, we just define the unions ourselves
# when codegen and the transpiler improves, we could just use FixedDicts with haskey
# which compile away.
const VertexWithColor = Union{Sprite, VertexCS, LineVertex}
const VertexWithThickness = Union{LineVertex}

const Sprite2 = Sprite{2, Float32}
const Sprite3 = Sprite{3, Float32}

to_vec4(v::Vec3f0) = Vec4f0(v[1], v[2], v[3], 1f0)
to_vec4(v::Vec2f0) = Vec4f0(v[1], v[2], 0f0, 1f0)

getuvrect(x::AbstractVertex) = Vec4f0(-0.5f0, -0.5f0, 0.5f0, 0.5f0)
getuvrect(x::Sprite) = x.uv

getuv(x::VertexNUV) = x.uv
getuv(x::AbstractVertex) = Vec2f0(0)

getcolor(v::AbstractVertex, uniforms) = uniforms.color
getcolor(v::VertexWithColor, uniforms) = v.color
getstrokecolor(x) = x.strokecolor
getglowcolor(x) = x.glowcolor
getscale(x) = x.scale

getposition(x::Vec4f0) = x
getposition(x::AbstractVertex) = x.position

getnormal(x) = x.normal

getthickness(x::VertexWithThickness, uniforms) = x.thickness
getthickness(x, uniforms) = uniforms.thickness

function aastep{T}(threshold1::T, value)
    afwidth = norm(Vec2f0(dFdx(value), dFdy(value))) * T(1.05);
    smoothstep(threshold1 - afwidth, threshold1 + afwidth, value)
end

function aastep{T}(threshold1::T, threshold2::T, value::T)
    afwidth = norm(Vec2f0(dFdx(value), dFdy(value))) * T(1.05);
    return (
        smoothstep(threshold1 - afwidth, threshold1 + afwidth, value) -
        smoothstep(threshold2 - afwidth, threshold2 + afwidth, value)
    )
end


text:
using Test, Suppressor, GNSSDecoder

include("test_data.jl")
include("bin_2_dec.jl")
include("find_preamble.jl")
include("decode_words.jl")
include("parity_check.jl")


text:
"""
    struct SubVector{T,A<:AbstractVector{T}} <: AbstractVector{T}
      vector::A
      pini::Int
      pend::Int
    end
"""
struct SubVector{T,A<:AbstractVector{T}} <: AbstractVector{T}
  vector::A
  pini::Int
  pend::Int
end

Base.size(a::SubVector) = (1+a.pend-a.pini,)

@propagate_inbounds function Base.getindex(a::SubVector,i::Integer)
  a.vector[i+a.pini-1]
end

@propagate_inbounds function Base.setindex!(a::SubVector,v,i::Integer)
  a.vector[i+a.pini-1] = v
end

Base.IndexStyle(::Type{SubVector{T,A}}) where {T,A} = IndexLinear()

Base.unaliascopy(A::SubVector) = typeof(A)(unaliascopy(A.vector), A.pini, A.pend)


text:
__precompile__()
module LightML

using Gadfly
using DataFrames
using ForwardDiff
using Distributions
using PyCall
using PyPlot
using DataStructures
using Distances
using Clustering

@pyimport sklearn.datasets as dat

export

    test_LinearRegression,
    test_LogisticRegression,

    test_ClassificationTree,
    test_RegressionTree,

    test_GDA,
    test_HMM,

    test_kneast_regression,
    test_kneast_classification,

    test_label_propagation,

    test_LDA,
    test_LDA_reduction,

    test_naive,

    test_NeuralNetwork,

    test_svm,

    test_GaussianMixture,

    test_kmeans_random,
    test_kmeans_speed,

    test_PCA,

    test_spec_cluster,

    test_Adaboost,
    test_BoostingTree,
    test_randomForest,

    make_cla,
    make_reg,
    make_digits,
    make_blo,
    make_iris,

    train!,
    predict,

    RegressionTree,
    ClassificationTree,
    randomForest,
    Adaboost,
    LinearRegression,
    LogisticRegression,
    GDA,
    BoostingTree,
    KnnRegression,
    KnnClassifier,
    label_propagation,
    show_example,
    LDA,
    plot_in_2d,
    NaiveBayes,
    NeuralNetwork,
    svm,
    GaussianMixture,
    Kmeans,
    PCA,
    transform,
    spec_clustering

Features = Union{String, Real}

#Supervised_learning

include("supervised_learning/baseRegression.jl")
include("supervised_learning/decisionTree.jl")
include("supervised_learning/gaussianDiscriminantAnalysis.jl")
include("supervised_learning/hiddenMarkovModel.jl")
include("supervised_learning/kNearestNeighbors.jl")
include("supervised_learning/labelPropagation.jl")
include("supervised_learning/linearDiscriminantAnalysis.jl")
include("supervised_learning/naivdBayes.jl")
include("supervised_learning/neuralNetwork_bp.jl")
include("supervised_learning/support_vector_machine.jl")
include("supervised_learning/adaboost.jl")
include("supervised_learning/randomForests.jl")
include("supervised_learning/xgboost.jl")
include("supervised_learning/GradientBoostingTree.jl")

#Unsupervised_learning

include("unsupervised_learning/gaussianMixtureModel.jl")
include("unsupervised_learning/kMeans.jl")
include("unsupervised_learning/principalComponentAnalysis.jl")
include("unsupervised_learning/spectralCluster.jl")
include("unsupervised_learning/largeScaleSpectralClustering.jl")

#Utils

include("utils/utils.jl")
include("utils/demo.jl")

end


text:
using ChaosTools
using DynamicalSystemsBase
using Test

@testset "Tipping points" begin
    xg = yg = range(-4, 4, length = 100)
    ds = Systems.magnetic_pendulum(d=0.2, α=0.2, ω=0.8, N=3)
    basins_before, = basins_of_attraction((xg, yg), ds; idxs=1:2,
    diffeq=(;:reltol => 1e-9), show_progress = false)
    ds = Systems.magnetic_pendulum(d=0.2, α=0.2, ω=0.8, N=3, γs = [1.0, 1.0, 0.05])
    basins_after, = basins_of_attraction((xg, yg), ds; idxs = 1:2,
    diffeq=(;:reltol => 1e-9), show_progress = false)
    @testset "Basin fractions" begin
        fs = basin_fractions(basins_before)
        @test sum(values(fs)) ≈ 1
        @test all(0.32 .≤ values(fs) .≤ 0.34)
        fs_after = basin_fractions(basins_after)
        @test length(fs_after) == 2
    end
    @testset "tipping probabilities" begin
        P = tipping_probabilities(basins_before, basins_after)
        @test size(P) == (3,2)
        @test all(0.495 .≤ P[1, :] .≤ 0.505) # must be ≈ 50%
        for i in 1:3; @test sum(P[i, :]) ≈ 1; end
    end
end