text:
# Unit testing for hybrid expressions

using Base.Test
using Devectorize

#data

a = [1., 2., 3., 4.]
b = [3., 4., 5., 6.]
c = [9., 8., 7., 6.]
cv = 12.0
cv_a = [12., 12., 12., 12.]

abc = [a b c]

at = a'
bt = b'
ct = c'
abct = [at; bt; ct]

# test cases

@devec r = sum(a) + 2.
@test isequal(r, sum(a) + 2.)

@devec r = sqrt(sum(sqr(a)))
@test isequal(r, sqrt(sum(a .* a)))

@devec r = sum(abc, 1) + maximum(c) .* cv
@test isequal(r, sum(abc, 1) .+ maximum(c) .* cv)

@devec r = sum(maximum(abc, 2)) + sum(abct, 1) .* (cv .+ maximum(b))
@test isequal(r, sum(maximum(abc, 2)) .+ sum(abct, 1) .* (cv .+ maximum(b)) )


text:
A=[0 2;0 0];
t=(-5,0.5,5);

X1=repmat([t;]',length(t));
X2=repmat([t;],1,length(t));

#Let exp(A*t)=E
E=[1 2.*t;0 1];

using PyPlot

quiver(X1,X2,E,1)
plot(X[:,1],X[:,2],label="Problem_9_1",colour="r")

ext: jl

max_stars_repo_name: akshayhiregoudar/Control-System-Design

max_issues_repo_name: akshayhiregoudar/Control-System-Design

max_forks_repo_name: akshayhiregoudar/Control-System-Design


text:
"""
    Tanh(x)

Hyperbolic tangent activation function.

    Usage:
```julia-repl
Tanh(x)
```
Arguments:
- `x` : Input value.
"""
function Tanh(x::Float64)
    x₁ = tanh(x)
    return x₁
end


text:
module RigidBodyDynamicsDiff

export
    DifferentialCache,
    mass_matrix_differential!,
    dynamics_bias_differential!,
    inverse_dynamics_differential!,
    dynamics_differential!

using RigidBodyDynamics
using ForwardDiff
using LinearAlgebra
using StaticArrays

using RigidBodyDynamics: predsucc, successorid, isdirty, Spatial.colwise, Spatial.hat, Spatial.se3_commutator
using RigidBodyDynamics: update_transforms!, update_bias_accelerations_wrt_world!, update_twists_wrt_world!
using RigidBodyDynamics: update_motion_subspaces!, update_spatial_inertias!, update_crb_inertias!
using RigidBodyDynamics: configuration_index_to_joint_id, velocity_index_to_joint_id, supports
using RigidBodyDynamics: velocity_to_configuration_derivative_jacobian, velocity_to_configuration_derivative_jacobian!
using RigidBodyDynamics.CustomCollections: SegmentedBlockDiagonalMatrix
using ForwardDiff: Dual

# TODO:
# * more frame checks for timederiv methods

function timederiv(H::Transform3D, twist::Twist, ::Type{Tag}) where Tag
    # @framecheck H.from twist.body
    # @framecheck H.to twist.base
    # @framecheck twist.frame twist.base
    R = rotation(H)
    p = translation(H)
    ω = angular(twist)
    v = linear(twist)
    dR = colwise(×, ω, R)
    dp = ω × p + v
    T = promote_type(eltype(dR), eltype(dp))
    dH = hcat(
        vcat(dR, zeros(similar_type(dR, Size(1, 3)))),
        vcat(dp, zeros(similar_type(dR, Size(1, 1))))
    )
    Transform3D(H.from, H.to,
        map(Dual{Tag}, H.mat, dH))
end

function timederiv(J::GeometricJacobian, twist::Twist, ::Type{Tag}) where Tag
    # TODO: more frame checks needed?
    @framecheck J.frame twist.frame
    ω = angular(twist)
    v = linear(twist)
    Jω = angular(J)
    Jv = linear(J)
    dJω = colwise(×, ω, Jω)
    dJv = colwise(×, v, Jω) + colwise(×, ω, Jv)
    GeometricJacobian(J.body, J.base, J.frame,
        map(Dual{Tag}, Jω, dJω),
        map(Dual{Tag}, Jv, dJv))
end

function timederiv(twist1::Twist, twist2::Twist, ::Type{Tag}) where Tag
    @framecheck twist1.frame twist2.frame
    dangular, dlinear = se3_commutator(angular(twist2), linear(twist2), angular(twist1), linear(twist1))
    Twist(twist1.body, twist1.base, twist1.frame,
        map(Dual{Tag}, angular(twist1), dangular),
        map(Dual{Tag}, linear(twist1), dlinear))
end

function timederiv(accel::SpatialAcceleration, twist::Twist, ::Type{Tag}) where Tag
    @framecheck accel.frame twist.frame
    dangular, dlinear = se3_commutator(angular(accel), linear(accel), angular(twist), linear(twist))
    SpatialAcceleration(accel.body, accel.base, accel.frame,
        map(Dual{Tag}, angular(accel), dangular),
        map(Dual{Tag}, linear(accel), dlinear))
end

function timederiv(inertia::SpatialInertia, twist::Twist, ::Type{Tag}) where Tag
    # TODO: more frame checks needed?
    @framecheck inertia.frame twist.frame
    ω = angular(twist)
    v = linear(twist)
    Iω = inertia.moment
    mc = inertia.cross_part
    m = inertia.mass
    mv = m * v
    dmc = mv + ω × mc
    ωhat = hat(ω)
    vhat = hat(v)
    mchat = hat(mc)
    dIω = ωhat * Iω - Iω * ωhat - vhat * mchat - mchat * vhat
    SpatialInertia(inertia.frame,
        map(Dual{Tag}, Iω, dIω),
        map(Dual{Tag}, mc, dmc),
        Dual{Tag}(m, zero(m)))
end

struct DifferentialCache{Tag, T, DS<:MechanismState{Dual{Tag, T, 1}}, DR<:DynamicsResult{Dual{Tag, T, 1}}}
    dualstates::Vector{DS}
    dualresults::Vector{DR}
end

function DifferentialCache{Tag}(mechanism::Mechanism{T}) where {Tag, T}
    if !mapreduce(has_fixed_subspaces, &, tree_joints(mechanism), init=true)
        throw(ArgumentError("Can only handle Mechanisms with fixed motion subspaces."))
    end
    D = ForwardDiff.Dual{Tag, T, 1}
    n = Threads.nthreads()
    dualstates = Vector{typeof(MechanismState{D}(mechanism))}(undef, n)
    dualresults = Vector{typeof(DynamicsResult{D}(mechanism))}(undef, n)
    Threads.@threads for i = 1 : n
        # Create MechanismStates and DynamicsResults in separate threads to avoid the possibility of false sharing.
        id = Threads.threadid()
        dualstate = dualstates[id] = MechanismState{D}(mechanism)
        dualresults[id] = DynamicsResult{D}(mechanism)
        dualstate.q .= NaN
        dualstate.v .= NaN
        dualstate.s .= NaN
    end
    DifferentialCache(dualstates, dualresults)
end

function copy_differential_column!(dest::Matrix, src::AbstractVector{<:Dual}, destcol::Integer)
    deststart = (destcol - 1) * size(dest, 1)
    @inbounds @simd for row in eachindex(src)
        dest[deststart + row] = ForwardDiff.partials(src[row], 1)
    end
end

function copy_differential_column!(dest::Matrix, src::Symmetric{<:Dual}, destcol::Integer)
    deststart = (destcol - 1) * size(dest, 1)
    n = size(src, 1)
    upper = src.uplo === 'u'
    @inbounds for col in Base.OneTo(n)
        @simd for row in col : n
            u_index = (row - 1) * n + col
            l_index = (col - 1) * n + row
            data_index = ifelse(upper, u_index, l_index)
            val = ForwardDiff.partials(src.data[data_index], 1)
            dest[deststart + u_index] = val
            dest[deststart + l_index] = val
        end
    end
end

function dual_state_init!(dualstate::MechanismState{<:Dual{Tag}}, state::MechanismState{T}, v_index::Integer,
        transforms::Bool, motion_subspaces::Bool, inertias::Bool, crb_inertias::Bool,
        twists::Bool, bias_accelerations::Bool) where {Tag, T}
    dualstate.q .= NaN # TODO: actually compute this?
    copyto!(dualstate.v, state.v)
    copyto!(dualstate.s, state.s)
    setdirty!(dualstate)
    k = velocity_index_to_joint_id(state, v_index)
    twist = Twist(state.motion_subspaces.data[v_index], SVector(one(T)))

    # TODO:
    if bias_accelerations
        transforms = true
    end

    if transforms
        dtransforms_to_root = dualstate.transforms_to_root
        @inbounds for i in state.treejointids
            i_successor = successorid(i, state)
            H = transform_to_root(state, i_successor, false)
            if supports(k, i_successor, state)
                dtransforms_to_root[i_successor] = timederiv(H, twist, Tag)
            else
                dtransforms_to_root[i_successor] = H
            end
        end
        dualstate.transforms_to_root.dirty = false
    end

    if motion_subspaces
        dmotion_subspaces = dualstate.motion_subspaces.data
        @inbounds for i in state.treejointids
            i_successor = successorid(i, state)
            for v_index_i in velocity_range(state, i)
                S = state.motion_subspaces.data[v_index_i]
                if supports(k, i_successor, state)
                    dmotion_subspaces[v_index_i] = timederiv(S, twist, Tag)
                else
                    dmotion_subspaces[v_index_i] = S
                end
            end
        end
        dualstate.motion_subspaces.dirty = false
    end

    if inertias
        dinertias = dualstate.inertias
        @inbounds for i in state.treejointids
            i_successor = successorid(i, state)
            I_i = state.inertias[i_successor]
            if supports(k, i_successor, state)
                dinertias[i_successor] = timederiv(I_i, twist, Tag)
            else
                dinertias[i_successor] = I_i
            end
        end
        dualstate.inertias.dirty = false
    end

    if crb_inertias
        dcrb_inertias = dualstate.crb_inertias
        k_successor = successorid(k, state)
        @inbounds for i in state.treejointids
            i_successor = successorid(i, state)
            Ic_i = state.crb_inertias[i_successor]
            if supports(k, i_successor, state) || supports(i, k_successor, state)
                p = max(i, k)
                Ic_p = state.crb_inertias[successorid(p, state)]
                Ic_p_dual = timederiv(Ic_p, twist, Tag)
                Iω = map(Dual{Tag}, Ic_i.moment, map(ForwardDiff.partials, Ic_p_dual.moment))
                mc = map(Dual{Tag}, Ic_i.cross_part, map(ForwardDiff.partials, Ic_p_dual.cross_part))
                m = Dual{Tag}(Ic_i.mass, ForwardDiff.partials(Ic_p_dual.mass))
                dcrb_inertias[i_successor] = SpatialInertia(Ic_i.frame, Iω, mc, m)
            else
                dcrb_inertias[i_successor] = Ic_i
            end
        end
        dualstate.crb_inertias.dirty = false
    end

    if twists
        dtwists = dualstate.twists_wrt_world
        @inbounds for i in state.treejointids
            i_successor = successorid(i, state)
            twist_i = twist_wrt_world(state, i_successor, false)
            if supports(k, i_successor, state)
                relative_twist_dual = timederiv(relative_twist(state, i_successor, first(predsucc(k, state))), twist, Tag)
                ω = map(Dual{Tag}, twist_i.angular, map(ForwardDiff.partials, relative_twist_dual.angular))
                v = map(Dual{Tag}, twist_i.linear, map(ForwardDiff.partials, relative_twist_dual.linear))
                dtwists[i_successor] = Twist(twist_i.body, twist_i.base, twist_i.frame, ω, v)
            else
                dtwists[i_successor] = twist_i
            end
        end
        dualstate.twists_wrt_world.dirty = false
    end

    # if bias_accelerations
    #     dbias_accelerations = dualstate.bias_accelerations_wrt_world
    #     @inbounds for i in state.treejointids
    #         κ_i = state.ancestor_joint_masks[i]
    #         i_successor = successorid(i, state)
    #         accel_i = bias_acceleration(state, i_successor, false)
    #         if κ_i[k]
    #             dbias_accelerations[i_successor] = timederiv(accel_i, twist, Tag)
    #         else
    #             dbias_accelerations[i_successor] = accel_i
    #         end
    #     end
    #     dualstate.bias_accelerations_wrt_world.dirty = false
    # end

    nothing
end

function threaded_differential!(f::F, dest::Matrix, state::MechanismState, cache::DifferentialCache;
        transforms::Bool=false, motion_subspaces::Bool=false, inertias::Bool=false, crb_inertias::Bool=false, twists::Bool=false, bias_accelerations::Bool=false) where F
    let state = state
        transforms && update_transforms!(state)
        motion_subspaces && update_motion_subspaces!(state)
        inertias && update_spatial_inertias!(state)
        crb_inertias && update_crb_inertias!(state)
        twists && update_twists_wrt_world!(state)
        bias_accelerations && update_bias_accelerations_wrt_world!(state)
    end
    let dest = dest, cache = cache, state = state
        nv = num_velocities(state)
        Threads.@threads for v_index in Base.OneTo(nv)
            id = Threads.threadid()
            dualstate = cache.dualstates[id]
            dual_state_init!(dualstate, state, v_index, transforms, motion_subspaces,
                inertias, crb_inertias, twists, bias_accelerations)
            dualresult = cache.dualresults[id]
            copy_differential_column!(dest, f(dualresult, dualstate), v_index)
            nothing
        end
    end
    return dest
end

function mass_matrix_differential!(dest::Matrix, state::MechanismState, cache::DifferentialCache)
    nv = num_velocities(state)
    size(dest) == (nv^2, nv) || throw(DimensionMismatch())
    threaded_differential!(mass_matrix!, dest, state, cache;
        motion_subspaces=true, crb_inertias=true)
end

function dynamics_bias_differential!(dest::Matrix, state::MechanismState, cache::DifferentialCache)
    nv = num_velocities(state)
    size(dest) == (nv, nv) || throw(DimensionMismatch())
    threaded_differential!(dynamics_bias!, dest, state, cache;
        motion_subspaces=true, inertias=true, bias_accelerations=true, twists=true)
end

function inverse_dynamics_differential!(dest::Matrix, state::MechanismState, v̇::SegmentedVector{JointID}, cache::DifferentialCache)
    nv = num_velocities(state)
    size(dest) == (nv, nv) || throw(DimensionMismatch())
    f = let v̇ = v̇
        function (result, state)
            # TODO: abusing the dynamicsbias and totalwrenches fields of DynamicsResult.
            inverse_dynamics!(result.dynamicsbias, result.jointwrenches, result.accelerations, state, v̇, result.totalwrenches)
            result.dynamicsbias
        end
    end
    threaded_differential!(f, dest, state, cache;
        transforms=true, motion_subspaces=true, inertias=true, bias_accelerations=true, twists=true)
end

function dynamics_differential!(dest::Matrix, result::DynamicsResult, state::MechanismState, τ::SegmentedVector{JointID}, cache::DifferentialCache)
    nv = num_velocities(state)
    size(dest) == (nv, nv) || throw(DimensionMismatch())
    dynamics!(result, state, τ)
    inverse_dynamics_differential!(dest, state, result.v̇, cache)
    LinearAlgebra.LAPACK.potrs!(result.massmatrix.uplo, result.L, dest)
    @inbounds map!(-, dest, dest)
    return dest
end

end # module


text:
ndims(a)             # a 的维度
eltype(a)            # a 中元素的类型
length(a)            # a 中元素数目
size(a)              # a 每一维度的大小所构成的元组
sum(a, dims=3)       # 对第三维求和，相似地，该用法可以用于以下函数： mean, std,
                     # prod, minimum, maximum, any, all；
                     # 使用统计函数前需要 using Statistics

max_stars_repo_name: zhaiyusci/Julia-Notes

max_issues_repo_name: zhaiyusci/Julia-Notes

max_forks_repo_name: zhaiyusci/Julia-Notes


text:
using LinearAlgebra
using StaticArrays
using Test

using CLIMA
using CLIMA.Atmos
using CLIMA.GenericCallbacks
using CLIMA.ODESolvers
using CLIMA.Mesh.Filters
using CLIMA.Mesh.Grids
using CLIMA.MoistThermodynamics
using CLIMA.PlanetParameters
using CLIMA.VariableTemplates

const p_ground = MSLP
const T_initial = 255
const domain_height = 30e3

function init_heldsuarez!(bl, state, aux, coords, t)
    global p_ground, T_initial

    FT = eltype(state)

    # TODO: change to altitude function
    r = norm(coords, 2)
    h = r - FT(planet_radius)
    # h = altitude(bl.orientation, aux)

    scale_height = R_d * FT(T_initial) / grav
    p            = FT(p_ground) * exp(-h / scale_height)

    state.ρ  = air_density(FT(T_initial), p)
    state.ρu = SVector{3, FT}(0, 0, 0)
    state.ρe = state.ρ * (internal_energy(FT(T_initial)) + aux.orientation.Φ)

    return nothing
end

function config_heldsuarez(FT, N, resolution)
    global T_initial, domain_height

    ref_state = HydrostaticState(IsothermalProfile(FT(T_initial)), FT(0))
    model = AtmosModel{FT}(AtmosGCMConfiguration;
                           ref_state  = ref_state,
                           turbulence = ConstantViscosityWithDivergence(FT(0)),
                           moisture   = DryModel(),
                           source     = (Gravity(), Coriolis(), held_suarez_forcing!),
                           init_state = init_heldsuarez!)
    config = CLIMA.Atmos_GCM_Configuration("HeldSuarez", N, resolution,
                                           FT(domain_height),
                                           init_heldsuarez!;
                                           model = model)

    return config
end

function held_suarez_forcing!(bl, source, state, diffusive, aux, t::Real)
    global T_initial

    FT = eltype(state)

    ρ     = state.ρ
    ρu    = state.ρu
    ρe    = state.ρe
    coord = aux.coord
    Φ     = aux.orientation.Φ
    e     = ρe / ρ
    u     = ρu / ρ
    e_int = e - u' * u / 2 - Φ
    T     = air_temperature(e_int)

    # Held-Suarez constants
    k_a       = FT(1 / (40 * day))
    k_f       = FT(1 / day)
    k_s       = FT(1 / (4 * day)) # TODO: day is actually seconds per day; should it be named better?
    ΔT_y      = FT(60)
    Δθ_z      = FT(10)
    T_equator = FT(315)
    T_min     = FT(200)

    σ_b          = FT(7 / 10)
    r            = norm(coord, 2)
    @inbounds λ  = atan(coord[2], coord[1])
    @inbounds φ  = asin(coord[3] / r)
    h            = r - FT(planet_radius)
    scale_height = R_d * FT(T_initial) / grav
    σ            = exp(-h / scale_height)

    # TODO: use
    #  p = air_pressure(T, ρ)
    #  σ = p/p0
    exner_p       = σ ^ (R_d / cp_d)
    Δσ            = (σ - σ_b) / (1 - σ_b)
    height_factor = max(0, Δσ)
    T_equil       = (T_equator - ΔT_y * sin(φ) ^ 2 - Δθ_z * log(σ) * cos(φ) ^ 2 ) * exner_p
    T_equil       = max(T_min, T_equil)
    k_T           = k_a + (k_s - k_a) * height_factor * cos(φ) ^ 4
    k_v           = k_f * height_factor

    source.ρu += -k_v * ρu
    source.ρe += -k_T * ρ * cv_d * (T - T_equil)
end

function main()
    CLIMA.init()

    FT = Float64

    # DG polynomial order
    N = 5

    # Domain resolution
    nelem_horz = 6
    nelem_vert = 8
    resolution = (nelem_horz, nelem_vert)

    t0 = FT(0)
    timeend = FT(60) # 400day

    driver_config = config_heldsuarez(FT, N, resolution)
    ode_solver_type = CLIMA.ExplicitSolverType(solver_method=LSRK144NiegemannDiehlBusch)
    solver_config = CLIMA.setup_solver(t0, timeend, driver_config,
                                       ode_solver_type=ode_solver_type)

    # Set up the filter callback
    filterorder = 14
    filter = ExponentialFilter(solver_config.dg.grid, 0, filterorder)
    cbfilter = GenericCallbacks.EveryXSimulationSteps(1) do
        Filters.apply!(solver_config.Q, 1:size(solver_config.Q, 2),
                       solver_config.dg.grid, filter)
        nothing
    end

    result = CLIMA.invoke!(solver_config;
                           user_callbacks=(cbfilter,),
                          check_euclidean_distance=true)
end

main()

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
using WoodburyMatrices
using Test

seed!(123)
n = 5

for elty in (Float32, Float64, ComplexF32, ComplexF64, Int), AMat in (x -> Matrix(Diagonal(x)),)

    elty = Float64

    a = rand(n); B = rand(n,2); D = rand(2,2); v = rand(n)

    if elty == Int
        v = rand(1:100, n)
        a = rand(1:100, n)
        B = rand(1:100, 2, n)
        D = rand(1:100, 2, 2)
    else
        v = convert(Vector{elty}, v)
        a = convert(Vector{elty}, a)
        B = convert(Matrix{elty}, B)
        D = convert(Matrix{elty}, D)
    end

    ε = eps(abs2(float(one(elty))))
    A = AMat(a)

    # Woodbury
    for W in (SymWoodbury(A, B, D), SymWoodbury(A, B[:,1][:], 2.))

        F = Matrix(W)
        @test (2*W)*v ≈ 2*(W*v)
        @test W'*v ≈ W*v
        @test (W'W)*v ≈ Matrix(W)*(Matrix(W)*v)
        @test (W*W)*v ≈ Matrix(W)*(Matrix(W)*v)
        @test (W*W')*v ≈ Matrix(W)*(Matrix(W)*v)
        @test W[1:3,1:3]*v[1:3] ≈ Matrix(W)[1:3,1:3]*v[1:3]
        @test sparse(W) ≈ Matrix(W)
        @test W === W'
        @test W*Matrix(1.0I, n, n) ≈ Matrix(W)
        @test W'*Matrix(1.0I, n, n) ≈ Matrix(W)

        Z = randn(n,n)
        @test Matrix(W*Z) ≈ Matrix(W)*Z

        R = rand(n,n)

        for v = (rand(n, 1), view(rand(n,1), 1:n), view(rand(n,2),1:n,1:2))
            @test (2*W)*v ≈ 2*(W*v)
            @test (W*2)*v ≈ 2*(W*v)
            @test (W'W)*v ≈ Matrix(W)*(Matrix(W)*v)
            @test (W*W)*v ≈ Matrix(W)*(Matrix(W)*v)
            @test (W*W')*v ≈ Matrix(W)*(Matrix(W)*v)
            @test W[1:3,1:3]*v[1:3] ≈ Matrix(W)[1:3,1:3]*v[1:3]
            @test Matrix(WoodburyMatrices.conjm(W, R)) ≈ R*Matrix(W)*R'
            @test Matrix(copy(W)'W)*v ≈ Matrix(W)*(Matrix(W)*v)
            @test Matrix(W + A) ≈ Matrix(W)+Matrix(A)
            @test Matrix(A + W) ≈ Matrix(W)+Matrix(A)
        end

        v = rand(n,1)
        W2 = convert(Woodbury, W)
        @test Matrix(W2) ≈ Matrix(W)

        if elty != Int
            @test inv(W)*v ≈ inv(Matrix(W))*v
            @test W\v ≈ inv(Matrix(W))*v
            @test liftFactor(W)(v) ≈ inv(W)*v
            @test WoodburyMatrices.partialInv(W)[1] ≈ inv(W).B
            @test WoodburyMatrices.partialInv(W)[2] ≈ inv(W).D
            @test det(W) ≈ det(Matrix(W))
        end

    end

end

for elty in (Float32, Float64, ComplexF32, ComplexF64, Int)

    elty = Float64

    a1 = rand(n); B1 = rand(n,2); D1 = rand(2,2); v = rand(n)
    a2 = rand(n); B2 = rand(n,2); D2 = rand(2,2);

    if elty == Int
        v = rand(1:100, n)

        a1 = rand(1:100, n)
        B1 = rand(1:100, 2, n)
        D1 = rand(1:100, 2, 2)

        a2 = rand(1:100, n)
        B2 = rand(1:100, 2, n)
        D2 = rand(1:100, 2, 2)
    else
        v = convert(Vector{elty}, v)

        a1 = convert(Vector{elty}, a1)
        B1 = convert(Matrix{elty}, B1)
        D1 = convert(Matrix{elty}, D1)

        a2 = convert(Vector{elty}, a2)
        B2 = convert(Matrix{elty}, B2)
        D2 = convert(Matrix{elty}, D2)
    end

    ε = eps(abs2(float(one(elty))))

    # Woodbury
    A1 = Matrix(Diagonal((a1)))
    A2 = Matrix(Diagonal((a2)))

    W1 = SymWoodbury(A1, B1, D1)
    W2 = SymWoodbury(A2, B2, D2)

    W1r = SymWoodbury(A1, B1[:,1][:], 2.)
    W2r = SymWoodbury(A2, B2[:,1][:], 3.)

    for (W1, W2) = ((W1,W2), (W1r, W2), (W1, W2r), (W1r,W2r))
        @test (W1 + W2)*v ≈ (Matrix(W1) + Matrix(W2))*v
        @test (Matrix(W1) + W2)*v ≈ (Matrix(W1) + Matrix(W2))*v
        @test (W1 + 2*Matrix(Diagonal((a1))))*v ≈ (Matrix(W1) + Matrix(2*Matrix(Diagonal((a1)))))*v
        @test_throws MethodError W1*W2
    end

end

# Sparse U and D

A = Matrix(Diagonal((rand(n))))
B = sprandn(n,2,1.)
D = sprandn(2,2,1.)
W = SymWoodbury(A, B, D)
v = randn(n)
V = randn(n,1)

@test size(W) == (n,n)
@test size(W,1) == n
@test size(W,2) == n

@test inv(W)*v ≈ inv(Matrix(W))*v
@test (2*W)*v ≈ 2*(W*v)
@test (W'W)*v ≈ Matrix(W)*(Matrix(W)*v)
@test (W*W)*v ≈ Matrix(W)*(Matrix(W)*v)
@test (W*W')*v ≈ Matrix(W)*(Matrix(W)*v)
@test liftFactor(W)(v) ≈ inv(W)*v

@test inv(W)*V ≈ inv(Matrix(W))*V
@test (2*W)*V ≈ 2*(W*V)
@test (W'W)*V ≈ Matrix(W)*(Matrix(W)*V)
@test (W*W)*V ≈ Matrix(W)*(Matrix(W)*V)
@test (W*W')*V ≈ Matrix(W)*(Matrix(W)*V)

# Mismatched sizes
@test_throws DimensionMismatch SymWoodbury(rand(5,5),rand(5,2),rand(2,3))
@test_throws DimensionMismatch SymWoodbury(rand(5,5),rand(5,2),rand(3,3))
@test_throws DimensionMismatch SymWoodbury(rand(5,5),rand(3),1.)


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.

# Julia Array related interface

"""
    similar(x::NDArray; writable, ctx)

Create an `NDArray` with similar shape, data type,
and context with the given one.
Note that the returned `NDArray` is uninitialized.
"""
Base.similar(x::NDArray{T,N}; writable = x.writable, ctx = context(x)) where {T,N} =
  NDArray{T,N}(undef, size(x)...; writable = writable, ctx = ctx)

"""
    zeros([DType], dims, ctx::Context = current_context())
    zeros([DType], dims...)
    zeros(x::NDArray)

Create zero-ed `NDArray` with specific shape and type.
"""
function zeros(::Type{T}, dims::NTuple{N,Int},
               ctx::Context = current_context()) where {N,T<:DType}
  x = NDArray{T}(undef, dims..., ctx = ctx)
  x[:] = zero(T)
  x
end

zeros(::Type{T}, dims::Int...) where {T<:DType} = zeros(T, dims)

zeros(dims::NTuple{N,Int}, ctx::Context = current_context()) where N =
  zeros(MX_float, dims, ctx)
zeros(dims::Int...) = zeros(dims)

zeros(x::NDArray)::typeof(x)      = zeros_like(x)
Base.zeros(x::NDArray)::typeof(x) = zeros_like(x)

"""
    ones([DType], dims, ctx::Context = current_context())
    ones([DType], dims...)
    ones(x::NDArray)

Create an `NDArray` with specific shape & type, and initialize with 1.
"""
function ones(::Type{T}, dims::NTuple{N,Int},
              ctx::Context = current_context()) where {N,T<:DType}
  arr = NDArray{T}(undef, dims..., ctx = ctx)
  arr[:] = one(T)
  arr
end

ones(::Type{T}, dims::Int...) where T<:DType = ones(T, dims)

ones(dims::NTuple{N,Int}, ctx::Context = current_context()) where N =
  ones(MX_float, dims, ctx)
ones(dims::Int...) = ones(dims)

ones(x::NDArray)::typeof(x)      = ones_like(x)
Base.ones(x::NDArray)::typeof(x) = ones_like(x)

import Base: length, ndims

"""
    size(x::NDArray)
    size(x::NDArray, dims)

Get the shape of an `NDArray`. The shape is in Julia's column-major convention.
See also the notes on NDArray shapes [`NDArray`](@ref).
"""
function Base.size(x::NDArray)
  ref_ndim  = Ref{MX_uint}(0)
  ref_shape = Ref{Ptr{MX_uint}}(0)
  @mxcall(:MXNDArrayGetShape, (MX_handle, Ref{MX_uint}, Ref{Ptr{MX_uint}}),
          x, ref_ndim, ref_shape)
  tuple(map(Int, reverse(unsafe_wrap(Array, ref_shape[], ref_ndim[])))...)
end

Base.size(x::NDArray{T,N}, dims::Integer) where {T,N} = (dims > N) ? 1 : size(x)[dims]

"""
    length(x::NDArray)

Get the number of elements in an `NDArray`.
"""
length(x::NDArray) = prod(size(x))

"""
    ndims(x::NDArray)

Get the number of dimensions of an `NDArray`.
Is equivalent to `length(size(arr))`.
"""
ndims(x::NDArray) = ndims(x.handle)

function ndims(x::MX_NDArrayHandle)::Int
  ref_ndim  = Ref{MX_uint}(0)
  ref_shape = Ref{Ptr{MX_uint}}(0)
  @mxcall(:MXNDArrayGetShape, (MX_handle, Ref{MX_uint}, Ref{Ptr{MX_uint}}),
          x, ref_ndim, ref_shape)
  ref_ndim[]
end

"""
    eltype(x::NDArray)

Get the element type of an `NDArray`.
"""
function Base.eltype(x::Union{NDArray,MX_NDArrayHandle})
  dtype_ref = Ref{Cint}(0)
  @mxcall(:MXNDArrayGetDType, (MX_handle, Ptr{Cint}), x, dtype_ref)

  if dtype_ref[] == -1 # x->is_none()
    # TODO: unit test for this branch
    throw(MXError("Eltype of $x is not defined"))
  end

  fromTypeFlag(TypeFlag(dtype_ref[]))
end

@inline _first(x::NDArray) = try_get_shared(x, sync = :read) |> first

Base.first(x::NDArray) = _first(x)

Base.lastindex(x::NDArray) = length(x)

"""
    slice(arr :: NDArray, start:stop)

Create a view into a sub-slice of an `NDArray`. Note only slicing at the slowest
changing dimension is supported. In Julia's column-major perspective, this is the last
dimension. For example, given an `NDArray` of shape (2,3,4), `slice(array, 2:3)` will create
a `NDArray` of shape (2,3,2), sharing the data with the original array. This operation is
used in data parallelization to split mini-batch into sub-batches for different devices.
"""
function slice(arr::NDArray, ::Colon)
  arr
end
function slice(arr::NDArray, slice::UnitRange{Int})
  dim1 = size(arr)[end]
  @assert(1 <= slice.start <= slice.stop <= dim1)
  if slice.start == 1 && slice.stop == dim1
    return arr
  end

  hdr_ref = Ref{MX_handle}(0)
  # note Julia is 1-based, inclusive-inclusive indexing, while C++ is
  # 0-based, inclusive-exclusive indexing. So 1:3 in Julia should
  # translates into 0:3 in C++.
  @mxcall(:MXNDArraySlice, (MX_handle, MX_uint, MX_uint, Ref{MX_handle}),
          arr, slice.start-1, slice.stop, hdr_ref)
  return NDArray(MX_NDArrayHandle(hdr_ref[]), arr.writable)
end

function _at(handle::Union{MX_NDArrayHandle, MX_handle}, idx::Integer)
  h_ref = Ref{MX_handle}(C_NULL)
  @mxcall(:MXNDArrayAt, (MX_handle, MX_uint, Ref{MX_handle}),
          handle, idx, h_ref)
  h_ref[]
end

import Base: setindex!

"""
    setindex!(arr::NDArray, val, idx)

Assign values to an `NDArray`.
The following scenarios are supported

* single value assignment via linear indexing: `arr[42] = 24`

* `arr[:] = val`: whole array assignment, `val` could be a scalar or an array (Julia `Array`
  or `NDArray`) of the same shape.
* `arr[start:stop] = val`: assignment to a *slice*, `val` could be a scalar or an array of
  the same shape to the slice. See also [`slice`](@ref).
"""
function setindex!(arr::NDArray, val::Real, idx::Integer)
  # linear indexing
  @assert arr.writable
  _set_value(out=arr[idx], src=val)
end

function setindex!(arr::NDArray, val::Real, ::Colon)
  @assert arr.writable
  _set_value(out = arr, src = dump_mx_param(val))
end

function setindex!(arr::NDArray, val::Array{T}, ::Colon) where T<:Real
  @assert arr.writable
  copy!(arr, val)
end

function setindex!(arr::NDArray, val::NDArray, ::Colon)
  @assert arr.writable
  copy!(arr, val)
end

function setindex!(arr::NDArray, val::Union{T,Array{T},NDArray},
                   idx::UnitRange{Int}) where T<:Real
  @assert arr.writable
  setindex!(slice(arr, idx), val, Colon())
end

import Base: getindex
"""
    getindex(arr::NDArray, idx)

Shortcut for [`slice`](@ref). A typical use is to write

```julia
  arr[:] += 5
```

which translates into

```julia
  arr[:] = arr[:] + 5
```

which furthur translates into

```julia
  setindex!(getindex(arr, Colon()), 5, Colon())
```

!!! note
    The behavior is quite different from indexing into Julia's `Array`. For example, `arr[2:5]`
    create a **copy** of the sub-array for Julia `Array`, while for `NDArray`, this is
    a *slice* that shares the memory.
"""
getindex(arr::NDArray, ::Colon) = arr

"""
Shortcut for [`slice`](@ref).
**NOTE** the behavior for Julia's built-in index slicing is to create a
copy of the sub-array, while here we simply call `slice`,
which shares the underlying memory.
"""
getindex(arr::NDArray, idx::UnitRange{Int}) = slice(arr, idx)

getindex(arr::NDArray) = _first(arr)

function getindex(arr::NDArray, idx::Integer)
  # linear indexing
  len = length(arr)
  size_ = size(arr)

  if idx <= 0 || idx > len
    throw(BoundsError(
      "attempt to access $(join(size_, 'x')) NDArray at index $(idx)"))
  end

  idx -= 1
  offsets = size_[1:end-1] |> reverse ∘ cumprod ∘ collect
  handle = arr.handle
  for offset ∈ offsets
    handle = _at(handle, idx ÷ offset)
    idx %= offset
  end

  _at(handle, idx) |> MX_NDArrayHandle |> x -> NDArray(x, arr.writable)
end

import Base: copy!, copy, convert, deepcopy

"""
    copy!(dst::Union{NDArray, Array}, src::Union{NDArray, Array})

Copy contents of `src` into `dst`.
"""
function copy!(dst::NDArray, src::NDArray)
  @assert(dst.writable)
  if dst.handle == src.handle
    @warn("Copying an NDArray to itself")
    return
  end

  _copyto(src, out=dst)
  return dst
end

function copy!(dst::Array{T}, src::NDArray{T}) where T<:DType
  @assert size(dst) == size(src)
  @mxcall(:MXNDArraySyncCopyToCPU, (MX_handle, Ptr{Cvoid}, Csize_t),
          src, pointer(dst), length(dst))
  dst
end

copy!(dst::Array{<:Real}, src::NDArray) = copy!(dst, copy(src))
copy!(dst::NDArray, src::AbstractArray) = copy!(dst, collect(src))

function copy!(dst::NDArray{T}, src::Array{<:Real}) where {T}
  @assert dst.writable
  @assert size(dst) == size(src)
  src = convert(Array{T}, src) # this might involve copying
  @mxcall(:MXNDArraySyncCopyFromCPU, (MX_handle, Ptr{Cvoid}, Csize_t),
          dst.handle, pointer(src), length(src))
  dst
end

function copy_ignore_shape!(dst::NDArray{T}, src::Array{<:Real}) where {T}
  @assert dst.writable
  @assert length(dst) == length(src)
  src = convert(Array{T}, src) # this might involve copying
  @mxcall(:MXNDArraySyncCopyFromCPU, (MX_handle, Ptr{Cvoid}, Csize_t),
          dst.handle, pointer(src), length(src))
  dst
end

"""
    copy(arr :: NDArray)
    copy(arr :: NDArray, ctx :: Context)
    copy(arr :: Array, ctx :: Context)

Create a copy of an array. When no `Context` is given, create a Julia `Array`.
Otherwise, create an `NDArray` on the specified context.
"""
copy

# Create copy: NDArray -> Julia Array
copy(x::NDArray{T,D}) where{T,D} = copy!(Array{T,D}(undef, size(x)), x)

# Create copy: NDArray -> NDArray in a given context
copy(x::NDArray{T,D}, ctx::Context) where {T,D} =
  copy!(NDArray{T,D}(_ndarray_alloc(T, size(x), ctx, true)), x)

# Create copy: Julia Array -> NDArray in a given context
copy(x::Array{T}, ctx::Context) where {T<:DType} =
  copy!(NDArray{T}(undef, size(x); ctx = ctx), x)

copy(x::AbstractArray, ctx::Context) =
  copy!(NDArray{eltype(x)}(undef, size(x); ctx = ctx), collect(x))

"""
    convert(::Type{Array{<:Real}}, x::NDArray)

Convert an `NDArray` into a Julia `Array` of specific type.
Data will be copied.
"""
convert(T::Type{Array{<:Real}}, x::NDArray) = convert(T, copy(x))

"""
    deepcopy(arr::NDArray)

Get a deep copy of the data blob in the form of an NDArray of default storage
type. This function blocks. Do not use it in performance critical code.
"""
function deepcopy(arr::NDArray)
  out_ref = Ref{MX_handle}(C_NULL)
  @mxcall(:MXNDArrayGetDataNDArray, (MX_handle, Ref{MX_handle}), arr, out_ref)
  NDArray(MX_NDArrayHandle(out_ref[]))
end

"""
    hcat(x::NDArray...)
"""
Base.hcat(xs::NDArray{T}...) where T = cat(xs..., dims = 2)

"""
    vcat(x::NDArray...)
"""
Base.vcat(xs::NDArray{T}...) where T = cat(xs..., dims = 1)

"""
    cat(xs::NDArray...; dims)

Concate the `NDArray`s which have the same element type along the `dims`.
Building a diagonal matrix is not supported yet.
"""
function Base.cat(xs::NDArray{T}...; dims) where T
  ns = ndims.(xs)
  d = Base.max(dims, maximum(ns))
  xs′ = map(zip(ns, xs)) do i
    n, x = i
    (d > n) ? reshape(x, -2, Base.ones(Int, d - n)...) : x
  end
  concat(xs′..., dim = d - dims)
end

"""
    @inplace

Julia does not support re-definiton of `+=` operator (like `__iadd__` in python),
When one write `a += b`, it gets translated to `a = a+b`. `a+b` will allocate new
memory for the results, and the newly allocated `NDArray` object is then assigned
back to a, while the original contents in a is discarded. This is very inefficient
when we want to do inplace update.

This macro is a simple utility to implement this behavior. Write

```julia
  @mx.inplace a += b
```

will translate into

```julia
  mx.add_to!(a, b)
```

which will do inplace adding of the contents of `b` into `a`.
"""
macro inplace(ex)
  f = if ex.head == :+= || ex.head == :.+=
    :add_to!
  elseif ex.head == :-= || ex.head == :.-=
    :sub_from!
  elseif ex.head == :.*=
    :mul_to!
  elseif ex.head == :./=
    :div_from!
  elseif ex.head == :.%=
    :mod_from!
  else
    error("unsupported inplace translation for $ex")
  end
  Expr(:call, f, esc(ex.args[1]), esc(ex.args[2]))
end

"""
    add_to!(dst::NDArray, args::NDArrayOrReal...)

Add a bunch of arguments into `dst`. Inplace updating.
"""
function add_to!(dst::NDArray, args::NDArrayOrReal...)
  @assert dst.writable
  for arg in args
    if isa(arg, Real)
      _plus_scalar(dst, scalar = arg, out = dst)
    else
      _plus!(dst, arg)
    end
  end
  dst
end

"""
    fill!(arr::NDArray, x)

Create an `NDArray` filled with the value `x`, like `Base.fill!`.
"""
function Base.fill!(arr::NDArray, x)
  arr[:] = x
  arr
end

"""
    fill(x, dims, ctx = current_context())
    fill(x, dims...)

Create an `NDArray` filled with the value `x`, like `Base.fill`.
"""
function fill(x::T, dims::NTuple{N,Integer}, ctx::Context = current_context()) where {T,N}
  arr = NDArray{T}(undef, dims, ctx = ctx)
  arr[:] = x
  arr
end

fill(x, dims::Integer...) = fill(x, dims)

import Base: hypot

broadcasted(::typeof(hypot), x::NDArray{T}, y::NDArray{T}) where {T} =
  _broadcast_hypot(x, y)

"""
Manipulating as Julia Arrays
----------------------------

    @nd_as_jl(captures..., statement)

A convenient macro that allows to operate `NDArray` as Julia Arrays. For example,

```julia
  x = mx.zeros(3,4)
  y = mx.ones(3,4)
  z = mx.zeros((3,4), mx.gpu())

  @mx.nd_as_jl ro=(x,y) rw=z begin
    # now x, y, z are just ordinary Julia Arrays
    z[:,1] = y[:,2]
    z[:,2] = 5
  end
```

Under the hood, the macro convert all the declared captures from `NDArray` into Julia
Arrays, by using `try_get_shared`. And automatically commit the modifications back into
the `NDArray` that is declared as `rw`. This is useful for fast prototyping and when
implement non-critical computations, such as `AbstractEvalMetric`.

!!! note
* Multiple `rw` and / or `ro` capture declaration could be made.
* The macro does **not** check to make sure that `ro` captures are not modified. If the
  original `NDArray` lives in CPU memory, then it is very likely the corresponding
  Julia Array shares data with the `NDArray`, so modifying the Julia Array will also
  modify the underlying `NDArray`.
* More importantly, since the `NDArray` is
  asynchronized, we will wait for *writing* for `rw` variables but wait only for *reading*
  in `ro` variables. If we write into those `ro` variables, **and** if the memory is
  shared, racing condition might happen, and the behavior is undefined.
* When an `NDArray` is declared to be captured as `rw`, its contents is always sync
  back in the end.
* The execution results of the expanded macro is always `nothing`.
* The statements are wrapped in a `let`, thus locally introduced new variables will not be
  available after the statements. So you will need to declare the variables before calling the
  macro if needed.
"""
macro nd_as_jl(m_args...)
  @assert(length(m_args) > 0)
  stmts = m_args[end]
  @assert(isa(stmts, Expr) && stmts.head == :block,
          "The last argument should be a statement block (begin-end); but get $stmts")
  stmts = esc(stmts)

  dclrs  = m_args[1:end-1]
  nd_ro  = []
  nd_rw  = []
  nd_all = []
  for declr in dclrs
    @assert(isa(declr, Expr) && declr.head == :(=) && length(declr.args)==2 && declr.args[1] ∈ (:ro,:rw),
            "Invalid declaration, should be rw=(x,y) or ro=z; but get $declr")

    declr_vars = declr.args[2]
    if isa(declr_vars, Symbol)
      declr_vars = (declr_vars,)
    elseif isa(declr_vars, Expr)
      @assert(declr_vars.head ∈ (:tuple, :vect),
              "Capture declaration should be a variable or a tuple of variables; but got $declr_vars")
      declr_vars = declr_vars.args
    else
      @assert(false, "Capture declaration should be a variable or a tuple of variables; but got $declr_vars")
    end
    for declr_var in declr_vars
      @assert(isa(declr_var, Symbol),
              "Captured ndarrays in ro/rw declaration should be variables, but get $(declr_var)")
    end
    append!(nd_all, [declr_vars...])
    if declr.args[1] == :ro
      append!(nd_ro, [declr_vars...])
    else
      append!(nd_rw, [declr_vars...])
    end
  end

  nd_ro    = map(esc, nd_ro)
  nd_rw    = map(esc, nd_rw)
  nd_all   = map(esc, nd_all)
  rw_origs = [gensym() for _ in nd_rw]

  save_statements  = Expr(:block, [:($v_orig = $v) for (v_orig, v) in zip(rw_origs, nd_rw)]...)
  wait_statements  = Expr(:block, [:(_wait_to_read($v)) for v in nd_ro]...,
                                  [:(_wait_to_write($v)) for v in nd_rw]...)
  clear_statements = Expr(:block, [:($v_orig = nothing) for v_orig in rw_origs]...)
  let_assignments  = Expr(:block, [:($v = try_get_shared($v)) for v in nd_all]...)
  sync_statements  = map(rw_origs, nd_rw) do v_orig, v
    quote
      if !is_shared($v, $v_orig)
        # copy data back if not or no longer sharing data
        copy!($v_orig, $v)
      end
    end
  end
  sync_statements  = Expr(:block, sync_statements...)

  let_statement = Expr(:let, let_assignments, quote
    $stmts
    $sync_statements
  end)
  m_body = quote
    $wait_statements
    $save_statements
    $let_statement
    $clear_statements
    nothing # the final results is always nothing
  end

  m_body
end

# NOTE: internal use only. Accessing pointers on a different device (e.g. accessing GPU
# pointers from CPU) leads to undefined behavior.
import Base.pointer
function pointer(arr :: NDArray)
  pdata = Ref{Ptr{Cvoid}}(0)
  @mxcall(:MXNDArrayGetData, (MX_handle, Ref{Ptr{Cvoid}}), arr, pdata)
  return convert(Ptr{eltype(arr)}, pdata[])
end

_ndsig[:reshape] = :(reshape(x; shape = dim, reverse = !reverse))
@_remap Base.reshape(x::NDArray, dim...; reverse = false) reshape
@_remap Base.reshape(x::NDArray, dim   ; reverse = false) reshape

_nddoc[:expand_dims] =
"""
    expand_dims(x::NDArray, dim)

Insert a new axis into `dim`.

```julia
julia> x
4 mx.NDArray{Float64,1} @ CPU0:
 1.0
 2.0
 3.0
 4.0

julia> mx.expand_dims(x, 1)
1×4 mx.NDArray{Float64,2} @ CPU0:
 1.0  2.0  3.0  4.0

julia> mx.expand_dims(x, 2)
4×1 mx.NDArray{Float64,2} @ CPU0:
 1.0
 2.0
 3.0
 4.0
```
"""
@_remap expand_dims(x::NDArray, dim) expand_dims(x; axis = -dim)

@_remap Base.permutedims(x::NDArray, axes) transpose(x; axes = length(axes) .- tuple(axes...))

_nddoc[:broadcast_to] = """
    broadcast_to(x::NDArray, dims)
    broadcast_to(x::NDArray, dims...)

Broadcasts the input array to a new shape.

In the case of broacasting doesn't work out of box,
you can expand the NDArray first.

```jldoctest
julia> x = mx.ones(2, 3, 4);

julia> y = mx.ones(1, 1, 4);

julia> x .+ mx.broadcast_to(y, 2, 3, 4)
2×3×4 mx.NDArray{Float32,3} @ CPU0:
[:, :, 1] =
 2.0  2.0  2.0
 2.0  2.0  2.0

[:, :, 2] =
 2.0  2.0  2.0
 2.0  2.0  2.0

[:, :, 3] =
 2.0  2.0  2.0
 2.0  2.0  2.0

[:, :, 4] =
 2.0  2.0  2.0
 2.0  2.0  2.0
```
"""
@_remap broadcast_to(x::NDArray, dims)    broadcast_to(x; shape = dims)
@_remap broadcast_to(x::NDArray, dims...) broadcast_to(x; shape = dims)

_nddoc[:broadcast_axis] = _nddoc[:broadcast_axes] = """
    broadcast_axis(x::NDArray, dim, size)
    broadcast_axes(x::NDArray, dim, size)

Broadcasts the input array over particular axis(axes).
Parameter `dim` and `size` could be a scalar, a Tuple or an Array.

`broadcast_axes` is just an alias.

```jldoctest
julia> x
1×2×1 mx.NDArray{Int64,3} @ CPU0:
[:, :, 1] =
 1  2

julia> mx.broadcast_axis(x, 1, 2)
2×2×1 mx.NDArray{Int64,3} @ CPU0:
[:, :, 1] =
 1  2
 1  2

julia> mx.broadcast_axis(x, 3, 2)
1×2×2 mx.NDArray{Int64,3} @ CPU0:
[:, :, 1] =
 1  2

[:, :, 2] =
 1  2
```
"""
@_remap(broadcast_axis(x::NDArray, dim, size),
        broadcast_axis(x; axis = ndims(x) .- dim, size = size))
@_remap(Base.broadcast_axes(x::NDArray, dim, size),
        broadcast_axes(x; axis = ndims(x) .- dim, size = size))

################################################################################
# remapping to solving type unstablility
################################################################################

@_remap _broadcast_hypot(x::NDArray, y::NDArray)  broadcast_hypot(x, y)
@_remap _broadcast_hypot!(x::NDArray, y::NDArray) broadcast_hypot(x, y)


text: sumRand(n::Int = 10) = sum(rand(n))


text:
function testFramebuffer(device, pass)
    lava.@scope begin
        info = lava.attachment2D(device.getPhysicalDevice(), 10, 10, vk.VK_FORMAT_D32_SFLOAT)
        lava.@autodestroy image1 = lava.Image(device, info.handleRef()[], vk.VK_IMAGE_VIEW_TYPE_2D)
        info = lava.attachment2D(device.getPhysicalDevice(), 10, 10, vk.VK_FORMAT_B8G8R8A8_SRGB)
        lava.@autodestroy image2 = lava.Image(device, info.handleRef()[], vk.VK_IMAGE_VIEW_TYPE_2D)
        image1.realizeAttachment()
        image2.realizeAttachment()
        lava.@autodestroy fb = lava.Framebuffer(pass, [image1, image2])
        @test length(fb.mViews) == 2
        @test fb.mViews[1].getImage() === image1 && fb.mViews[2].getImage() === image2
    end
    return true
end


text:
using PkgBenchmark
using HMMBase

# TODO: Add timestamp to result files
# results = benchmarkpkg(parent(pathof(HMMBase)))
results = benchmarkpkg(pwd(), resultfile = "benchmark.json")
export_markdown("benchmark.md", results)

max_stars_repo_name: Shu-Lea-Lai/module-MultiHMMBase


text:
mutable struct GuiData
    # export window
    window::GLFW.Window
    ctx::Ptr{Nothing}
    clear_color::Array{Cfloat,1}
    counter::Cint
    float_slide::Cfloat
    display_w::Int32
    display_h::Int32

    show_warning_dialog::Bool

    buffer::String

    function GuiData()
        o = new()

        o.buffer = ""

        o.window = GLFW.CreateWindow(GUI_WIDTH, GUI_HEIGHT, "Deuron7")
        @assert o.window != C_NULL
        GLFW.MakeContextCurrent(o.window)
        GLFW.SwapInterval(1)  # enable vsync

        o.ctx = CImGui.CreateContext()

        o.display_w, o.display_h = GLFW.GetFramebufferSize(o.window)

        # CImGui.StyleColorsDark()
        CImGui.StyleColorsClassic()
        # CImGui.StyleColorsLight()

        fonts_dir = joinpath(@__DIR__, "../gui/", "fonts")
        fonts = CImGui.GetIO().Fonts
        CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Roboto-Medium.ttf"), 20)

        ImGui_ImplGlfw_InitForOpenGL(o.window, true)
        ImGui_ImplOpenGL3_Init(glsl_version)

        o.clear_color = Cfloat[0.45, 0.55, 0.60, 1.00]
        o.show_warning_dialog = false

        o.counter = 0
        o.float_slide = 0.0

        o
    end
end


text:
using Distributed
#addprocs(3)

@everywhere using SeisRemoveEQ, SeisIO

using Dates, JLD2

#==================================================#

# Use SeisDownload.jl to make jld2 file.
finame = "./dataset/BPnetwork_3days_test.jld2"

MAX_MEM_PER_CPU = 2.0 # [GB] maximum allocated memory for one cpu

#---Kurtosis parameters---#
IsKurtosisRemoval = true

max_edgetaper_duration = 60 * 5 #[s] tapering to avoid edge effect of response removal
kurtosis_timewindow = 60 * 3 #300 # timewindow to compute kurtosis [s]: large number causes longer computational time.
kurtosis_threshold = 3.0 #kurtosis threshold to detect outlier

#---STA/LTA parameters---#
IsSTALTARemoval = true

stalta_longtimewindow = 60*60*24 #[s] long time window for STA/LTA: Note that short time window is same with removal_shorttimewindow
stalta_threshold = 0.8 #STA/LTA threshold to detect earthquake and tremors: For our purpose, this threshold is smaller than ordinal detection

#---EQ removal parameter---#
invert_tukey_α = 0.01 #invert_tukey window (0<α<1: 1.0 for smooth removal but causes incomplete tremors removal)
removal_shorttimewindow = 60 * 3 # timewindow to remove EQ: removing data if this tw contains data more than threshold.
overlap = 60; # overlap of removal window[s]
max_wintaper_duration = 60 * 3 # maximum duration of tukey window taper: τs

#---Plotting variables---#
IsSaveFig = true # save figure to check removal
plot_kurtosis_α = 1.2 # amplitude normalization for plotting kurtosis
plot_boxheight = 1.5
plot_span = 100 # to reduce number of plot point

#---Output file info---#
fodir       = "./dataset"
foname      = "BPnetwork_3days_test_neq" # data is saved at ./dataset/$foname.jld2

#==================================================#

#intended to parallelization
InputDictionary = Dict( "MAX_MEM_PER_CPU"        => float(MAX_MEM_PER_CPU),
                        "finame"                 => finame,
                        "IsKurtosisRemoval"      => IsKurtosisRemoval,
                        "max_edgetaper_duration" => float(max_edgetaper_duration),
                        "kurtosis_timewindow"    => float(kurtosis_timewindow),
                        "kurtosis_threshold"     => float(kurtosis_threshold),
                        "max_wintaper_duration"  => float(max_wintaper_duration),
                        "IsSTALTARemoval"        => IsSTALTARemoval,
                        "stalta_longtimewindow"  => float(stalta_longtimewindow),
                        "stalta_threshold"       => float(stalta_threshold),
                        "invert_tukey_α"         => float(invert_tukey_α),
                        "removal_shorttimewindow"=> float(removal_shorttimewindow),
                        "overlap"                => float(overlap),
                        "plot_kurtosis_α"        => float(plot_kurtosis_α),
                        "plot_boxheight"         => float(plot_boxheight),
                        "plot_span"              => trunc(Int, plot_span),
                        "fodir"                  => fodir,
                        "foname"                 => foname,
                        "IsSaveFig"              => IsSaveFig)

seisremoveEQ(InputDictionary)


text:
export PriorityQueue
export push!, pop!, top, height

mutable struct PriorityQueue{T}
    values::Vector{T}
    select

    size::Int64

    function PriorityQueue{T}(select, capacity::Int64 = 0) where {T}
        local pq = new()

        pq.values = Vector{T}()
        sizehint!(pq.values, capacity)
        pq.select = select
        pq.size = 0

        return pq
    end
end

function Base.:push!(pq::PriorityQueue{T}, value::T) where {T}
    local index = 1 + pq.size

    if length(pq.values) < index
        push!(pq.values, value)
    else
        pq.values[index] = value
    end

    pq.size += 1

    promote!(pq, index)
end

function Base.:pop!(pq::PriorityQueue{T}) where {T}
    @assert pq.size > 0

    local t = top(pq)

    if 1 == pq.size
        pq.size = 0
    else
        local last = pq.values[pq.size]
        pq.values[1] = last
        pq.size -= 1

        demote!(pq, 1)
    end

    return t
end

function top(pq::PriorityQueue{T}) where {T}
    @assert pq.size > 0

    return pq.values[1]
end

function height(pq::PriorityQueue{T}) where {T}
    return 1 + Int64(floor(log2(pq.size)))
end

function promote!(pq::PriorityQueue{T}, index::Int64) where {T}
    @assert index > 0 && index <= pq.size

    if 1 == index
        return
    end

    index -= 1
    local parent = Int64(floor(index / 2))

    while index > 0
        local t = pq.values[1+index]

        if t != pq.select(t, pq.values[1+parent])
            break
        end

        pq.values[1+index] = pq.values[1+parent]
        pq.values[1+parent] = t

        local next = parent
        parent = Int64(floor(index / 2))
        index = next
    end
end

function demote!(pq::PriorityQueue{T}, index::Int64) where {T}
    @assert index > 0 && index <= pq.size

    if pq.size == index
        return
    end

    local value = pq.values[index]

    index -= 1
    while index < pq.size
        local rv
        local right = (1 + index) * 2
        if right < pq.size
            rv = pq.values[1+right]
        end

        local lv
        local left = right - 1
        if left < pq.size
            lv = pq.values[1+left]
        end

        local child = -1
        if right < pq.size && left < pq.size && lv == pq.select(lv, rv)
            child = left
        elseif right < pq.size
            child = right
        elseif left < pq.size
            child = left
        end

        if child < 0 || value == pq.select(value, pq.values[1+child])
            break
        else
            pq.values[1+index] = pq.values[1+child]
            pq.values[1+child] = value
            index = child
        end
    end
end


text:
using Rankers.Storages
using LOLTools: ChampionMasteryV4, SummonerV4, http_action
using LOLToolsDataDragon.DataDragon.en_US

api_key = get(ENV, "RIOT_TOKEN", "")
region = "kr"

summonerName = "DRX Deft"
read_from_file = true # false
summonerdto = SummonerV4.summoner_by_name(api_key, region, summonerName; action=store_with(region, http_action, debug=true))
encryptedSummonerId = summonerdto.id

champion_masteries = ChampionMasteryV4.by_summoner_id(api_key, region, encryptedSummonerId; action=update_with(region, http_action, read_from_file=read_from_file, debug=true))
for mastery in sort(champion_masteries, by = x -> x.championPoints, rev=true)[1:5]
    @info (en_US.Champions[mastery.championId].name, mastery.championPoints, mastery.lastPlayTime)
end

championId = 114
champion_masteries_114 = ChampionMasteryV4.by_summoner_id(api_key, region, encryptedSummonerId, championId; action=update_with(region, http_action, read_from_file=read_from_file, debug=true))
@info champion_masteries_114
using Calendars.Times
@info :lastPlayTime millisecond2datetime(champion_masteries_114.lastPlayTime)


text:
using Distributions
using BlackBoxOptim

require("../solis_wets.jl")
require("../experiment_framework.jl")

gaussianrng = (sigma) -> Normal(0, sigma)
levyrng = (sigma) -> Levy(0, sigma)
cauchyrng = (sigma) -> Cauchy(0, sigma)

ps = BlackBoxOptim.as_fixed_dim_problem_set(BlackBoxOptim.example_problems, 2^5)

sf(p) = begin
  solis_wets(p; max_fail_steps = 3, rnggen = gaussianrng)
end

@time repeated_runs(sf, ps, 25; experiment = "solis_wets_2")


text:
module DiscontinuousFESpacesTests

using Test
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.FESpaces

domain =(0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)

order = 3
grid_topology = get_grid_topology(model)
polytopes = get_polytopes(grid_topology)
reffes = [LagrangianRefFE(Float64,p,order) for p in polytopes]

trian = get_triangulation(model)

V = DiscontinuousFESpace(reffes,trian)
test_single_field_fe_space(V)

U = TrialFESpace(V)

f(x) = sin(pi*x[1])*cos(2*pi*x[2])

fh = interpolate(f,U)

uh = FEFunction(V,rand(num_free_dofs(V)))

#using Gridap.Visualization
#
#writevtk(trian,"trian",nsubcells=40,cellfields=["fh"=>fh, "uh"=>uh])

reffes = [LagrangianRefFE(Float64,p,order,space=:P) for p in polytopes]
V = DiscontinuousFESpace(reffes,trian)
test_single_field_fe_space(V)

end # module


text:
# Command: Get Window Handles
"""
    window_handles(session::Session)::Vector{String}

Switching window will select the current top-level browsing context used as the target for all subsequent commands.
In a tabbed browser, this will typically make the tab containing the browsing context the selected tab.
"""
function window_handles(session::Session)::Vector{String}
    @unpack addr, id = session
    response = HTTP.get(
        "$addr/session/$id/window_handles",
        [("Content-Type" => "application/json; charset=utf-8")],
    )
    @assert response.status == 200
    convert(Vector{String}, JSON3.read(response.body).value)::Vector{String}
end


text:
using PyPlot, Distributed, Random
using Revise, transD_GP
##
img =     transD_GP.Img(
          filename         = "4.2.01.png",
          dx               = 10.0,
          fractrain        = 0.02,
          dec              = 2,
          gausskernelwidth = 7)
##
transD_GP.get_image_data(img,
                   sdmaxfrac = 0.05,
                   ybreak = 1000,
                   takeevery = 4)

transD_GP.plot_image_data(img)
transD_GP.calc_image_RMS(img)


text:
module ReactingToInput
using Rocket, GLFW,ModernGL,Setfield,  ..ReactToScroll,  ..ForDisplayStructs
using  ..TextureManag,DataTypesBasic,  ..ReactOnMouseClickAndDrag,  ..ReactOnKeyboard,  ..DataStructs,  ..StructsManag,  ..DisplayWords
using ..MaskDiffrence, ..KeyboardVisibility, ..OtherKeyboardActions, ..WindowControll, ..ChangePlane
export subscribeGLFWtoActor

"""
adding the data into about openGL and GLFW context to enable proper display of main image and masks
"""
function setUpMainDisplay(mainForDisplayObjects:: forDisplayObjects,actor::SyncActor{Any, ActorWithOpenGlObjects})
    actor.actor.mainForDisplayObjects.stopListening[]=true

    actor.actor.mainForDisplayObjects=mainForDisplayObjects
    actor.actor.mainForDisplayObjects.stopListening[]=false

end#setUpMainDisplay

"""
adding the data needed for text display; also activates appropriate quad for the display
    it also configures texture that is build for text display
"""
function setUpWordsDisplay(textDispObject:: ForWordsDispStruct,actor::SyncActor{Any, ActorWithOpenGlObjects})

    actor.actor.mainForDisplayObjects.stopListening[]=true

     bindAndActivateForText(textDispObject.shader_program_words
    , textDispObject.fragment_shader_words
    ,actor.actor.mainForDisplayObjects.vertex_shader
    ,textDispObject.vbo_words
    ,actor.actor.calcDimsStruct)

    texId =  createTexture(UInt8,actor.actor.calcDimsStruct.textTexturewidthh
                            ,actor.actor.calcDimsStruct.textTextureheightt
                            ,GL_R8UI,GL_UNSIGNED_BYTE)

    textSpec= setproperties(textDispObject.textureSpec,(ID=texId) )

    samplerRef= glGetUniformLocation(textDispObject.shader_program_words, "TextTexture1")

    glUniform1i(samplerRef,length(actor.actor.mainForDisplayObjects.listOfTextSpecifications)+1)
    textDispObjectiNITIALIZED= setproperties(textDispObject,(textureSpec=textSpec) )

    actor.actor.textDispObj=textDispObjectiNITIALIZED
    # now reactivating the main vbo and shader program

    reactivateMainObj(actor.actor.mainForDisplayObjects.shader_program
    , actor.actor.mainForDisplayObjects.vbo
    ,actor.actor.calcDimsStruct    )

    actor.actor.mainForDisplayObjects.stopListening[]=false

end#setUpWordsDisplay

"""
adding the data about 3 dimensional arrays that will be source of data used for scrolling behaviour
onScroll Data - list of tuples where first is the name of the texture that we provided and second is associated data (3 dimensional array of appropriate type)

"""
function setUpForScrollData(onScrollData::FullScrollableDat ,actor::SyncActor{Any, ActorWithOpenGlObjects})
    actor.actor.mainForDisplayObjects.stopListening[]=true

    onScrollData.slicesNumber= getSlicesNumber(onScrollData)
    actor.actor.onScrollData=onScrollData
    #In order to refresh all in case we would change the texture dimensions ...
    ChangePlane.processKeysInfo(Option(onScrollData.dataToScrollDims),actor,KeyboardStruct()  )
      #so  It will precalculate some data and later mouse modification will be swift
      oldd = actor.actor.valueForMasToSet

      actor.actor.valueForMasToSet = valueForMasToSetStruct(value = 0)
      ReactOnMouseClickAndDrag.reactToMouseDrag(MouseStruct(true,false, [CartesianIndex(5,5)]),actor )
      actor.actor.valueForMasToSet = oldd

  actor.actor.mainForDisplayObjects.stopListening[]=false

end#setUpMainDisplay

"""
add data needed for proper calculations of mouse, verticies positions ... etc
"""
function setUpCalcDimsStruct(calcDim::CalcDimsStruct ,actor::SyncActor{Any, ActorWithOpenGlObjects})
    actor.actor.mainForDisplayObjects.stopListening[]=true

    actor.actor.calcDimsStruct =calcDim

    actor.actor.mainForDisplayObjects.stopListening[]=false

end#setUpCalcDimsStruct

"""
sets value we are setting to the  active mask vie mause interaction, in case mask is modifiable
"""
function setUpvalueForMasToSet(valueForMasToSett::valueForMasToSetStruct ,actor::SyncActor{Any, ActorWithOpenGlObjects})
    actor.actor.mainForDisplayObjects.stopListening[]=true

    actor.actor.valueForMasToSet =valueForMasToSett

    updateImagesDisplayed(actor.actor.currentlyDispDat
    , actor.actor.mainForDisplayObjects
    , actor.actor.textDispObj
    , actor.actor.calcDimsStruct ,valueForMasToSett )

    actor.actor.mainForDisplayObjects.stopListening[]=false

end#setUpvalueForMasToSet

"""
enables updating just a single slice that is displayed - do not change what will happen after scrolling
one need to pass data to actor in
struct that holds tuple where first entry is
-vector of tuples whee first entry in tuple is name of texture given in the setup and second is 2 dimensional aray of appropriate type with image data
- Int - second is Int64 - that is marking the screen number to which we wan to set the actor state
"""
function updateSingleImagesDisplayedSetUp(singleSliceDat::SingleSliceDat ,actor::SyncActor{Any, ActorWithOpenGlObjects})
    actor.actor.mainForDisplayObjects.stopListening[]=true
    updateImagesDisplayed(singleSliceDat
                        , actor.actor.mainForDisplayObjects
                        , actor.actor.textDispObj
                        , actor.actor.calcDimsStruct
                        ,actor.actor.valueForMasToSet  )

    actor.actor.currentlyDispDat=singleSliceDat
    actor.actor.currentDisplayedSlice = singleSliceDat.sliceNumber
    actor.actor.isSliceChanged = true # mark for mouse interaction that we changed slice

    actor.actor.mainForDisplayObjects.stopListening[]=false

end #updateSingleImagesDisplayed

"""
configuring actor using multiple dispatch mechanism in order to connect input to proper functions; this is not
encapsulated by a function becouse this is configuration of Rocket and needs to be global
"""

Rocket.on_next!(actor::SyncActor{Any, ActorWithOpenGlObjects}, data::Int64) = reactToScroll(data,actor )
Rocket.on_next!(actor::SyncActor{Any, ActorWithOpenGlObjects}, data:: forDisplayObjects) = setUpMainDisplay(data,actor)
Rocket.on_next!(actor::SyncActor{Any, ActorWithOpenGlObjects}, data:: ForWordsDispStruct) = setUpWordsDisplay(data,actor)
Rocket.on_next!(actor::SyncActor{Any, ActorWithOpenGlObjects}, data::CalcDimsStruct) = setUpCalcDimsStruct(data,actor)
Rocket.on_next!(actor::SyncActor{Any, ActorWithOpenGlObjects}, data::valueForMasToSetStruct) = setUpvalueForMasToSet(data,actor)

Rocket.on_next!(actor::SyncActor{Any, ActorWithOpenGlObjects}, data::FullScrollableDat) = setUpForScrollData(data,actor)
Rocket.on_next!(actor::SyncActor{Any, ActorWithOpenGlObjects}, data::SingleSliceDat) = updateSingleImagesDisplayedSetUp(data,actor)
Rocket.on_next!(actor::SyncActor{Any, ActorWithOpenGlObjects}, data::MouseStruct) = reactToMouseDrag(data,actor)
Rocket.on_next!(actor::SyncActor{Any, ActorWithOpenGlObjects}, data::KeyboardStruct) = reactToKeyboard(data,actor)

Rocket.on_error!(actor::SyncActor{Any, ActorWithOpenGlObjects}, err)      = error(err)
Rocket.on_complete!(actor::SyncActor{Any, ActorWithOpenGlObjects})        = ""

"""
when GLFW context is ready we need to use this  function in order to register GLFW events to Rocket actor - we use subscription for this
    actor - Roctet actor that holds objects needed for display like window etc...
    return list of subscriptions so if we will need it we can unsubscribe
"""
function subscribeGLFWtoActor(actor ::SyncActor{Any, ActorWithOpenGlObjects})

    #controll scrolling
    forDisplayConstants = actor.actor.mainForDisplayObjects

    scrollback=  ReactToScroll.registerMouseScrollFunctions(forDisplayConstants.window,forDisplayConstants.stopListening,actor.actor.isBusy)
    GLFW.SetScrollCallback(forDisplayConstants.window, (a, xoff, yoff) -> scrollback(a, xoff, yoff))

    keyBoardAct = registerKeyboardFunctions(forDisplayConstants.window,forDisplayConstants.stopListening)
    buttonSubs  = registerMouseClickFunctions(forDisplayConstants.window,forDisplayConstants.stopListening,actor.actor.calcDimsStruct,actor.actor.isBusy )

    keyboardSub = subscribe!(keyBoardAct, actor)
    scrollSubscription = subscribe!(scrollback, actor)
    mouseClickSub = subscribe!(buttonSubs, actor)

return [scrollSubscription,mouseClickSub,keyboardSub]

end

end #ReactToGLFWInpuut


text:
using Test
import DynamicPolynomials: @polyvar
# using Revise
using PolynomialBasis

function allequal(v1,v2)
    return all(v1 .≈ v2)
end

points = [-1.0,1.0]
matrix = PolynomialBasis.vandermonde_matrix(points)
testmatrix = [1.0  -1.0
              1.0   1.0]
@test allequal(matrix,testmatrix)

points = [-1.0,0.0,1.0]
matrix = PolynomialBasis.vandermonde_matrix(points)
testmatrix = [1.0  -1.0   1.0
              1.0   0.0   0.0
              1.0   1.0   1.0]
@test allequal(matrix,testmatrix)

points = [0.0]
matrix = PolynomialBasis.vandermonde_matrix(points)
testmatrix = [1.0]
@test allequal(matrix,testmatrix)

coeffs = PolynomialBasis.lagrange_polynomial_coefficients([0.0])
testcoeffs = [1.0]
@test allequal(coeffs,testcoeffs)

coeffs = PolynomialBasis.lagrange_polynomial_coefficients([-1.0,1.0])
testcoeffs = [0.5  0.5
             -0.5  0.5]
@test allequal(coeffs,testcoeffs)

coeffs = PolynomialBasis.lagrange_polynomial_coefficients([-1.0,0.0,1.0])
testcoeffs = [+0.0   1.0   0.0
              -0.5   0.0   0.5
              +0.5  -1.0   0.5]
@test allequal(coeffs,testcoeffs)

@polyvar x
coeffs = [1.0]
poly = PolynomialBasis.polynomial_from_coefficients(x,coeffs)
powers = sort!(vcat(poly.x.Z...))
testpowers = [0]
@test allequal(powers,testpowers)
@test allequal(coeffs,sort(poly.a))

coeffs = [1.0,2.0]
poly = PolynomialBasis.polynomial_from_coefficients(x,coeffs)
powers = sort!(vcat(poly.x.Z...))
testpowers = [0,1]
@test allequal(powers,testpowers)
@test allequal(coeffs,sort(poly.a))

polys = PolynomialBasis.lagrange_polynomials(x,[0.0])
@test length(polys) == 1
@test polys[1] == 1.0

polys = PolynomialBasis.lagrange_polynomials(x,[-1.0,1.0])
@test length(polys) == 2
p1 = 0.5 - 0.5x
p2 = 0.5 + 0.5x
@test allequal([p1,p2],polys)

polys = PolynomialBasis.lagrange_polynomials(x,[-1.0,0.0,1.0])
@test length(polys) == 3
p1 = -0.5x + 0.5x^2
p2 = 1.0 - x^2
p3 = 0.5x + 0.5x^2
@test allequal([p1,p2,p3],polys)


text:
struct FaceToCellGlue{O} <: GridapType
  orientation::Val{O}
  face_to_cell::Vector{Int}
  face_to_lface::Vector{Int8}
  cell_to_ctype::Vector{Int8}
  cell_to_lface_to_pindex::Table{Int8,Vector{Int8},Vector{Int32}}
  ctype_to_lface_to_ftype::Vector{Vector{Int8}}
end

is_oriented(::FaceToCellGlue{O}) where O = O

"""
    struct GenericBoundaryTriangulation{Dc,Dp,Gf,Gc,O} <: BoundaryTriangulation{Dc,Dp}
      face_trian::Gf
      cell_trian::Gc
      # + private fields
    end
"""
struct GenericBoundaryTriangulation{Dc,Dp,Gf,Gc,O,A} <: BoundaryTriangulation{Dc,Dp}
  face_trian::Gf
  cell_trian::Gc
  glue::FaceToCellGlue{O}
  cell_around::A

  function GenericBoundaryTriangulation(
    face_trian::Triangulation,
    cell_trian::Triangulation,
    glue::FaceToCellGlue{O},
    cell_around=1) where O

    Dc = num_cell_dims(face_trian)
    Dp = num_point_dims(face_trian)
    Gf = typeof(face_trian)
    Gc = typeof(cell_trian)
    A = typeof(cell_around)
    new{Dc,Dp,Gf,Gc,O,A}(face_trian, cell_trian, glue, cell_around)
  end

end

function GenericBoundaryTriangulation(
  face_trian::Triangulation,
  cell_trian::Triangulation,
  topo::GridTopology,
  face_to_oldface::Vector{Int},
  icell_around=1)

  D = num_cell_dims(cell_trian)
  oldface_to_cells = get_faces(topo,D-1,D)
  cell_to_oldfaces = get_faces(topo,D,D-1)
  cell_to_lface_to_pindex = Table(get_cell_permutations(topo,D-1))
  orientation = OrientationStyle(topo)

  oldface_to_cell = get_local_item(oldface_to_cells, icell_around)
  oldface_to_lface = find_local_index(oldface_to_cell, cell_to_oldfaces)
  face_to_cell = collect(Int,reindex(oldface_to_cell, face_to_oldface))
  face_to_lface = collect(Int8,reindex(oldface_to_lface, face_to_oldface))

  f = (p)->fill(Int8(UNSET),num_faces(p,D-1))
  ctype_to_lface_to_ftype = map( f, get_reffes(cell_trian) )
  face_to_ftype = get_cell_type(face_trian)
  cell_to_ctype = collect(Int8,get_cell_type(cell_trian))
  _fill_ctype_to_lface_to_ftype!(
    ctype_to_lface_to_ftype,
    face_to_cell,
    face_to_lface,
    face_to_ftype,
    cell_to_ctype)

  glue = FaceToCellGlue(
    orientation,
    face_to_cell,
    face_to_lface,
    cell_to_ctype,
    cell_to_lface_to_pindex,
    ctype_to_lface_to_ftype)

  GenericBoundaryTriangulation(face_trian, cell_trian, glue, icell_around)

end

"""
    GenericBoundaryTriangulation(model::DiscreteModel,face_to_mask::Vector{Bool})
"""
function GenericBoundaryTriangulation(
  model::DiscreteModel,face_to_mask::Vector{Bool},icell_around=1)
  D = num_cell_dims(model)
  topo = get_grid_topology(model)
  oldface_grid = Grid(ReferenceFE{D-1},model)
  cell_grid = Grid(ReferenceFE{D},model)
  face_to_oldface = findall(face_to_mask)
  face_trian = TriangulationPortion(oldface_grid,face_to_oldface)
  GenericBoundaryTriangulation(face_trian,cell_grid,topo,face_to_oldface,icell_around)
end

function _fill_ctype_to_lface_to_ftype!(
  ctype_to_lface_to_ftype,
  face_to_cell,
  face_to_lface,
  face_to_ftype,
  cell_to_ctype)
  for (face, cell) in enumerate(face_to_cell)
    ctype = cell_to_ctype[cell]
    lface = face_to_lface[face]
    ftype = face_to_ftype[face]
    ctype_to_lface_to_ftype[ctype][lface] = ftype
  end
end

# Triangulation Interface

function get_cell_coordinates(trian::GenericBoundaryTriangulation)
  get_cell_coordinates(trian.face_trian)
end

function get_reffes(trian::GenericBoundaryTriangulation)
    get_reffes(trian.face_trian)
end

function get_cell_type(trian::GenericBoundaryTriangulation)
  get_cell_type(trian.face_trian)
end

function get_face_to_cell(trian::GenericBoundaryTriangulation)
  trian.glue.face_to_cell
end

function get_face_to_lface(trian::GenericBoundaryTriangulation)
  trian.glue.face_to_lface
end

function get_volume_triangulation(trian::GenericBoundaryTriangulation)
  trian.cell_trian
end

function get_face_to_face(trian::GenericBoundaryTriangulation)
  trian.face_trian.cell_to_oldcell
end

function get_cell_around(trian::GenericBoundaryTriangulation)
  trian.cell_around
end

# Cell to face map

function get_face_to_cell_map(trian::GenericBoundaryTriangulation)
  face_to_fvertex_to_qcoors = FaceCellCoordinates(trian)
  f = (p)-> get_shapefuns(LagrangianRefFE(Float64,get_polytope(p),1))
  ftype_to_shapefuns = map( f, get_reffes(trian.face_trian) )
  face_to_ftype = collect(Int8, get_cell_type(trian.face_trian))
  face_to_shapefuns = CompressedArray(ftype_to_shapefuns, face_to_ftype)
  lincomb(face_to_shapefuns,face_to_fvertex_to_qcoors)
end

struct FaceCellCoordinates{D,T,O} <: AbstractVector{Vector{Point{D,T}}}
  glue::FaceToCellGlue{O}
  cell_to_ctype::Vector{Int8}
  ctype_to_lvertex_to_qcoords::Vector{Vector{Point{D,T}}}
  ctype_to_lface_to_lvertices::Vector{Vector{Vector{Int}}}
  ctype_to_lface_to_pindex_to_perm::Vector{Vector{Vector{Vector{Int}}}}
  ctype_to_lface_to_pindex_to_qcoords::Vector{Vector{Vector{Vector{Point{D,T}}}}}
  function FaceCellCoordinates(trian::GenericBoundaryTriangulation)

    d = num_cell_dims(trian)
    polytopes = map(get_polytope, get_reffes(trian.cell_trian))
    cell_to_ctype = trian.glue.cell_to_ctype
    ctype_to_lvertex_to_qcoords = map(get_vertex_coordinates, polytopes)
    ctype_to_lface_to_lvertices = map((p)->get_faces(p,d,0), polytopes)
    ctype_to_lface_to_pindex_to_perm = map( (p)->get_face_vertex_permutations(p,d), polytopes)

    P = eltype(eltype(ctype_to_lvertex_to_qcoords))
    D = length(P)
    T = eltype(P)
    O = is_oriented(trian.glue)

    ctype_to_lface_to_pindex_to_qcoords = Vector{Vector{Vector{Point{D,T}}}}[]
    for (ctype, lface_to_pindex_to_perm) in enumerate(ctype_to_lface_to_pindex_to_perm)
      lvertex_to_qcoods = ctype_to_lvertex_to_qcoords[ctype]
      lface_to_pindex_to_qcoords = Vector{Vector{Point{D,T}}}[]
      for (lface, pindex_to_perm) in enumerate(lface_to_pindex_to_perm)
        cfvertex_to_lvertex = ctype_to_lface_to_lvertices[ctype][lface]
        nfvertices = length(cfvertex_to_lvertex)
        pindex_to_qcoords = Vector{Vector{Point{D,T}}}(undef,length(pindex_to_perm))
        for (pindex, cfvertex_to_ffvertex) in enumerate(pindex_to_perm)
          ffvertex_to_qcoords = zeros(Point{D,T},nfvertices)
          for (cfvertex, ffvertex) in enumerate(cfvertex_to_ffvertex)
            lvertex = cfvertex_to_lvertex[cfvertex]
            qcoords = lvertex_to_qcoods[lvertex]
            ffvertex_to_qcoords[ffvertex] = qcoords
          end
          pindex_to_qcoords[pindex] = ffvertex_to_qcoords
        end
        push!(lface_to_pindex_to_qcoords,pindex_to_qcoords)
      end
      push!(ctype_to_lface_to_pindex_to_qcoords,lface_to_pindex_to_qcoords)
    end

    new{D,T,O}(
      trian.glue,
      cell_to_ctype,
      ctype_to_lvertex_to_qcoords,
      ctype_to_lface_to_lvertices,
      ctype_to_lface_to_pindex_to_perm,
      ctype_to_lface_to_pindex_to_qcoords)

  end
end

Base.IndexStyle(::Type{FaceCellCoordinates{D,T,O}}) where {D,T,O} = IndexLinear()

Base.size(a::FaceCellCoordinates) = (length(a.glue.face_to_cell),)

function Base.getindex(a::FaceCellCoordinates,face::Integer)
  cell = a.glue.face_to_cell[face]
  lface = a.glue.face_to_lface[face]
  ctype = a.cell_to_ctype[cell]
  p = a.glue.cell_to_lface_to_pindex.ptrs[cell]-1
  pindex = a.glue.cell_to_lface_to_pindex.data[p+lface]
  a.ctype_to_lface_to_pindex_to_qcoords[ctype][lface][pindex]
end

# Value of face to cell map

function apply_lincomb(ax::CompressedArray,b::FaceCellCoordinates)
  FaceToCellMapValue(ax,b)
end

struct FaceToCellMapValue{D,T,O} <: AbstractVector{Vector{Point{D,T}}}
  b::FaceCellCoordinates{D,T,O}
  ctype_to_lface_to_pindex_to_qpoints::Vector{Vector{Vector{Vector{Point{D,T}}}}}

  function FaceToCellMapValue(ax::CompressedArray,b::FaceCellCoordinates{D,T,O}) where {D,T,O}

    # Precompute the integration point coordinates for all possible cases

    ctype_to_lface_to_ftype = b.glue.ctype_to_lface_to_ftype
    ctype_to_lface_to_lvertices = b.ctype_to_lface_to_lvertices
    ctype_to_lvertex_to_qcoords = b.ctype_to_lvertex_to_qcoords
    ctype_to_lface_to_pindex_to_perm = b.ctype_to_lface_to_pindex_to_perm
    ftype_to_shapefunsvals = ax.values

    # Allocate
    ctype_to_lface_to_pindex_to_qpoints = Vector{Vector{Vector{Point{D,T}}}}[]
    for (ctype, lface_to_pindex_to_perm) in enumerate(ctype_to_lface_to_pindex_to_perm)
      lface_to_pindex_to_qpoints = Vector{Vector{Point{D,T}}}[]
      for (lface, pindex_to_perm) in enumerate(lface_to_pindex_to_perm)
        pindex_to_qpoints = Vector{Vector{Point{D,T}}}(undef,length(pindex_to_perm))
        push!(lface_to_pindex_to_qpoints,pindex_to_qpoints)
      end
      push!(ctype_to_lface_to_pindex_to_qpoints,lface_to_pindex_to_qpoints)
    end

    # Fill
    for (ctype, lface_to_ftype) in enumerate(ctype_to_lface_to_ftype)
      for (lface, ftype) in enumerate(lface_to_ftype)
        if ftype == UNSET
          continue
        end

        cfvertex_to_lvertex = ctype_to_lface_to_lvertices[ctype][lface]
        pindex_to_perm = ctype_to_lface_to_pindex_to_perm[ctype][lface]
        lvertex_to_qcoods = ctype_to_lvertex_to_qcoords[ctype]
        nfvertices = length(cfvertex_to_lvertex)
        ffvertex_to_qcoords = zeros(Point{D,T},nfvertices)
        shapefunsvals = ftype_to_shapefunsvals[ftype]

        for (pindex, cfvertex_to_ffvertex) in enumerate(pindex_to_perm)
          for (cfvertex, ffvertex) in enumerate(cfvertex_to_ffvertex)
            lvertex = cfvertex_to_lvertex[cfvertex]
            qcoords = lvertex_to_qcoods[lvertex]
            ffvertex_to_qcoords[ffvertex] = qcoords
          end
          qpoints = shapefunsvals*ffvertex_to_qcoords
          ctype_to_lface_to_pindex_to_qpoints[ctype][lface][pindex] = qpoints
        end
      end
    end

    new{D,T,O}(b,ctype_to_lface_to_pindex_to_qpoints)

  end
end

Base.size(v::FaceToCellMapValue) = (length(v.b.glue.face_to_cell),)

Base.IndexStyle(::Type{FaceToCellMapValue{C,T,O}}) where {C,T,O} = IndexLinear()

function Base.getindex(a::FaceToCellMapValue,face::Integer)
  cell = a.b.glue.face_to_cell[face]
  lface = a.b.glue.face_to_lface[face]
  ctype = a.b.cell_to_ctype[cell]
  p = a.b.glue.cell_to_lface_to_pindex.ptrs[cell]-1
  pindex = a.b.glue.cell_to_lface_to_pindex.data[p+lface]
  a.ctype_to_lface_to_pindex_to_qpoints[ctype][lface][pindex]
end

function evaluate_field_array(g::CompressedArray{<:Field},fx::FaceToCellMapValue)
  CellShapeFunsAtFaces(g,fx)
end

struct CellShapeFunsAtFaces{D,T,O,V} <: AbstractVector{V}
  b::FaceCellCoordinates{D,T,O}
  ctype_to_lface_to_pindex_to_shfnsvals::Vector{Vector{Vector{V}}}
  function CellShapeFunsAtFaces(g::CompressedArray{<:Field},fx::FaceToCellMapValue{D,T,O}) where {D,T,O}

    # Precompute shapefuns for all cases

    ctype_to_shapefuns = g.values
    face_to_ctype = g.ptrs
    ctype_to_lface_to_pindex_to_perm = fx.b.ctype_to_lface_to_pindex_to_perm
    ctype_to_lface_to_pindex_to_qpoints = fx.ctype_to_lface_to_pindex_to_qpoints
    ctype_to_lface_to_ftype = fx.b.glue.ctype_to_lface_to_ftype

    # Allocate
    sfuns = first(ctype_to_shapefuns)
    V = field_return_type(sfuns,zeros(Point{D,T},1))
    ctype_to_lface_to_pindex_to_shfnsvals = Vector{Vector{V}}[]
    for (ctype, lface_to_pindex_to_perm) in enumerate(ctype_to_lface_to_pindex_to_perm)
      lface_to_pindex_to_shfnsvals = Vector{V}[]
      for (lface, pindex_to_perm) in enumerate(lface_to_pindex_to_perm)
        pindex_to_shfnsvals = Vector{V}(undef,length(pindex_to_perm))
        push!(lface_to_pindex_to_shfnsvals,pindex_to_shfnsvals)
      end
      push!(ctype_to_lface_to_pindex_to_shfnsvals,lface_to_pindex_to_shfnsvals)
    end

    # Fill
    for (ctype, lface_to_ftype) in enumerate(ctype_to_lface_to_ftype)
      for (lface, ftype) in enumerate(lface_to_ftype)
        if ftype == UNSET
          continue
        end
        pindex_to_perm = ctype_to_lface_to_pindex_to_perm[ctype][lface]
        sfuns = ctype_to_shapefuns[ctype]
        for (pindex, cfvertex_to_ffvertex) in enumerate(pindex_to_perm)
          qpoints = ctype_to_lface_to_pindex_to_qpoints[ctype][lface][pindex]
          shfnsvals = evaluate(sfuns,qpoints)
          ctype_to_lface_to_pindex_to_shfnsvals[ctype][lface][pindex] = shfnsvals
        end
      end
    end

    new{D,T,O,V}(fx.b,ctype_to_lface_to_pindex_to_shfnsvals)

  end
end

Base.size(v::CellShapeFunsAtFaces) = (length(v.b.glue.face_to_cell),)

Base.IndexStyle(::Type{CellShapeFunsAtFaces{C,T,O,V}}) where {C,T,O,V} = IndexLinear()

function Base.getindex(a::CellShapeFunsAtFaces,face::Integer)
  cell = a.b.glue.face_to_cell[face]
  lface = a.b.glue.face_to_lface[face]
  ctype = a.b.cell_to_ctype[cell]
  p = a.b.glue.cell_to_lface_to_pindex.ptrs[cell]-1
  pindex = a.b.glue.cell_to_lface_to_pindex.data[p+lface]
  a.ctype_to_lface_to_pindex_to_shfnsvals[ctype][lface][pindex]
end

# Normal vector

function get_normal_vector(trian::GenericBoundaryTriangulation)
  k = NormalVectorValued()
  refn = ReferenceNormal(trian)
  cell_map = restrict(get_cell_map(trian.cell_trian),trian)
  J = gradient(cell_map)
  a = apply(k,J,refn)
  GenericCellField(a,cell_map)
end

struct ReferenceNormal{D,T} <: AbstractVector{Point{D,T}}
  face_to_cell::Vector{Int}
  face_to_lface::Vector{Int8}
  cell_to_ctype::Vector{Int8}
  ctype_to_lface_to_qnormal::Vector{Vector{Point{D,T}}}
end

function ReferenceNormal(trian::GenericBoundaryTriangulation)
  face_to_cell = trian.glue.face_to_cell
  face_to_lface = trian.glue.face_to_lface
  cell_to_ctype = trian.glue.cell_to_ctype
  f = (r) -> get_facet_normals(get_polytope(r))
  ctype_to_lface_to_qnormal = map(f, get_reffes(trian.cell_trian))
  ReferenceNormal(
    face_to_cell,
    face_to_lface,
    cell_to_ctype,
    ctype_to_lface_to_qnormal)
end

Base.size(a::ReferenceNormal) = (length(a.face_to_cell),)

Base.IndexStyle(::Type{ReferenceNormal{D,T}}) where {D,T} = IndexLinear()

function Base.getindex(a::ReferenceNormal,face::Integer)
  cell = a.face_to_cell[face]
  lface = a.face_to_lface[face]
  ctype = a.cell_to_ctype[cell]
  a.ctype_to_lface_to_qnormal[ctype][lface]
end

struct NormalField <: Field end

struct NormalVectorValued <: Kernel end

function apply_kernel!(cache,k::NormalVectorValued,J,nref)
  NormalField()
end

function kernel_evaluate(k::NormalVectorValued,x,J,refn)
  Jx = evaluate_field_array(J,x)
  k = bcast(_map_normal)
  apply(k,Jx,refn)
end

function _map_normal(J::TensorValue{D,D,T},n::VectorValue{D,T}) where {D,T}
  v = inv(J)⋅n
  m = sqrt(inner(v,v))
  if m < eps()
    return zero(n)
  else
    return v/m
  end
end


text:
##########################################################################
# The TriangularHetNetNetwork is a cell network with three sites,
# placed on the vertices of an equilateral triangle. The simulation
# parameters are based on 3GPP Case 1 (TR 25.814 and TR 36.814) cell
# simulation environment. Distance-dependent pathloss and log-normal
# shadow fading is used. The shadow fading is correlated (with correlation
# coefficient 0.5) from each MS to the three BSs, but independent
# between MSs.
#
# The typical transmit power is 46 dBm over 10 MHz (600 subcarriers).
# This leads to 10*log10(10^(46/10)/600) = 18.2 dBm for one 15 kHz subcarrier.

##########################################################################
# Network definition
type TriangularHetNetNetwork{MS_t <: PhysicalMS, BS_t <: PhysicalBS, System_t <: System, PropagationEnvironment_t <: PropagationEnvironment} <: PhysicalNetwork
    MSs::Vector{MS_t}
    BSs::Vector{BS_t}

    system::System_t
    num_picos_per_cell::Int
    num_MSs_per_cell::Int
    propagation_environment::PropagationEnvironment_t
    inter_site_distance::Float64
    pico_centre_distance::Float64
    guard_distance::Float64
    aux_network_params::AuxNetworkParams

    assignment::Assignment
end

# Convenience constructor without network params and assignments
TriangularHetNetNetwork(MSs, BSs, system, num_picos_per_cell, num_MSs_per_cell, propagation_environment, inter_site_distance, pico_centre_distance, guard_distance) =
    TriangularHetNetNetwork(MSs, BSs, system, num_picos_per_cell, num_MSs_per_cell, propagation_environment, inter_site_distance, pico_centre_distance, guard_distance, AuxNetworkParams(), Assignment())

Base.show(io::IO, x::TriangularHetNetNetwork) =
    print(io, "TriangularHetNet(I = $(length(x.BSs)), num_MSs_per_cell = $(x.num_MSs_per_cell), ISD = $(x.inter_site_distance), GD = $(x.guard_distance))")
Base.showcompact(io::IO, x::TriangularHetNetNetwork) =
    print(io, "TriangularHetNet($(length(x.BSs)), $(x.num_MSs_per_cell), $(x.inter_site_distance), $(x.guard_distance))")

# The default parameter values are taken from 3GPP Case 1
# (TR 25.814 and TR 36.814).
function setup_triangularhetnet_network(
    num_picos_per_cell, num_MSs_per_cell, num_MS_antennas, num_BS_antennas;
    system = SinglecarrierSystem(2e9, 15e3),
    propagation_environment = SimpleLargescaleFadingEnvironment(37.6, 15.3, 20, 8),
    inter_site_distance = 500.,
    pico_centre_distance = 100.,
    guard_distance = 35.,
    transmit_power = 10^(18.2/10), transmit_powers = transmit_power*ones(Float64, 3 + 3*num_picos_per_cell),
    BS_antenna_gain_params =
        vcat([SixSector3gppAntennaParams(-deg2rad(90),  deg2rad(35), 23),
              SixSector3gppAntennaParams( deg2rad(30),  deg2rad(35), 23),
              SixSector3gppAntennaParams( deg2rad(150), deg2rad(35), 23)],
             [ OmnidirectionalAntennaParams(0) for idx = 1:3*num_picos_per_cell ]),
    user_priority = 1., user_priorities = user_priority*ones(Float64, 3*num_MSs_per_cell),
    num_streams = 1, num_streamss = num_streams*ones(Int, 3*num_MSs_per_cell),
    MS_antenna_gain_dB = 0., MS_antenna_gains_dB = MS_antenna_gain_dB*ones(Float64, 3*num_MSs_per_cell),
    receiver_noise_figure = 9., receiver_noise_figures = receiver_noise_figure*ones(Float64, 3*num_MSs_per_cell))

    # Consistency check
    isa(num_MS_antennas, Vector) || (num_MS_antennas = num_MS_antennas*ones(Int, 3*num_MSs_per_cell))
    isa(num_BS_antennas, Vector) || (num_BS_antennas = num_BS_antennas*ones(Int, 3 + 3*num_picos_per_cell))

    BSs = PhysicalBS[]

    # Place macrocells
    macro_centre_distance = (inter_site_distance/2)/cos(deg2rad(30))
    for i = 1:3
        push!(BSs, PhysicalBS(num_BS_antennas[i], rotate(Position(0, macro_centre_distance), deg2rad((i-1)*120)), transmit_powers[i], BS_antenna_gain_params[i]))
    end

    # Place picocells
    num_picos = 3*num_picos_per_cell; pico_rot_angle = deg2rad(360/num_picos)
    pico_base_angle = deg2rad(-(num_picos_per_cell-1)*360/num_picos/2)
    for i = 4:(3 + num_picos)
        push!(BSs, PhysicalBS(num_BS_antennas[i], rotate(Position(0, pico_centre_distance), pico_base_angle + (i-4)*pico_rot_angle), transmit_powers[i], BS_antenna_gain_params[i]))
    end

    MSs = [ PhysicalMS(num_MS_antennas[k], Position(0, 0), Velocity(0, 0), user_priorities[k], num_streamss[k], MS_antenna_gains_dB[k], receiver_noise_figures[k], SimpleLargescaleFadingEnvironmentState(zeros(Float64, 3 + 3*num_picos_per_cell), falses(3 + 3*num_picos_per_cell))) for k = 1:3*num_MSs_per_cell ]

    TriangularHetNetNetwork(MSs, BSs, system, num_picos_per_cell, num_MSs_per_cell,
        propagation_environment, inter_site_distance, pico_centre_distance, guard_distance)
end

##########################################################################
# Cell assignment

# The ID cell assignment only assigns MSs to the macro BSs. This network
# is then identical to the TriangularMacroNetwork.
function IDCellAssignment!(channel, network::TriangularHetNetNetwork)
    Kc = network.num_MSs_per_cell
    cell_assignment = Array(Int, 3*Kc)

    for i = 1:3
        cell_assignment[(i-1)*Kc+1:i*Kc] = i
    end

    network.assignment = Assignment(cell_assignment, get_num_BSs(network))

    return AssignmentResults()
end

# This is a greedy scheduler based on the large scale fading realizations.
function LargeScaleFadingCellAssignment!(channel, network::TriangularHetNetNetwork)
    num_MSs_per_cell = network.num_MSs_per_cell; num_MSs = 3*num_MSs_per_cell
    num_picos_per_cell = network.num_picos_per_cell; num_BSs = 3 + 3*num_picos_per_cell

    aux_params = get_aux_assignment_params(network)
    @defaultize_param! aux_params "max_num_MSs_per_BS" 1

    # Scheduling matrix
    cell_assignment_matrix = zeros(Int, num_MSs, num_BSs)

    # User selection metric
    F = (channel.large_scale_fading_factor.^2)*Diagonal(get_transmit_powers(network))
    Fsize = size(F)

    # Do not schedule users in the wrong cell
    for cell = 1:3
        BS_ids = vcat(cell, 3 .+ [(cell-1)*num_picos_per_cell+1:cell*num_picos_per_cell])
        other_cell_MS_ids = setdiff(1:num_MSs, [(cell-1)*num_MSs_per_cell+1:cell*num_MSs_per_cell])
        F[other_cell_MS_ids, BS_ids] = 0.
    end

    # Do greedy scheduling
    while !all(F .== 0.)
        _, idx = findmax(F)
        k, l = ind2sub(Fsize, idx)

        if sum(cell_assignment_matrix[:,l]) < aux_params["max_num_MSs_per_BS"]
            cell_assignment_matrix[k,l] = 1
            F[k,:] = 0.
        else
            F[:,l] = 0.
        end
    end

    network.assignment = Assignment(cell_assignment_matrix)

    return AssignmentResults()
end

##########################################################################
# Simulation functions
function draw_user_drop!{MS_t <: PhysicalMS, BS_t <: PhysicalBS, System_t <: System}(network::TriangularHetNetNetwork{MS_t, BS_t, System_t, SimpleLargescaleFadingEnvironment})
    I = get_num_BSs(network); K = get_num_MSs(network)

    # Shadow fading covariance (correlation coefficient 0.5 between macro cells, no correlation between pico cells)
    correlations = eye(Float64, I, I)
    correlations[1:3, 1:3] = (ones(3,3) + eye(3))/2
    Cov_shadow_sqrtm = network.propagation_environment.shadow_sigma_dB*chol(correlations)

    for k = 1:K
        # Generate user position within standard triangle [0,30] degrees, to
        # the right of the base. Apply guard distance to x coordinate.
        xtri = sqrt(((network.inter_site_distance/2)^2 - network.guard_distance^2)*rand() + network.guard_distance^2)
        ytri = (xtri/sqrt(3))*rand()

        # Flip it over with probability 0.5
        if rand() < 0.5
            pos = Position(xtri, ytri)
        else
            pos = Position(xtri, -ytri)
            pos = rotate(pos, deg2rad(60))
        end

        # Generate rotation angle
        i = @compat Int(div(k - 1, K/3) + 1) # serving BS id
        if i == 1
            θ = deg2rad(240)
        elseif i == 2
            θ = deg2rad(0)
        elseif i == 3
            θ = deg2rad(120)
        end

        # Get final user position
        network.MSs[k].position = network.BSs[i].position + rotate(pos, θ)

        # Shadow fading
        network.MSs[k].propagation_environment_state =
            SimpleLargescaleFadingEnvironmentState(Cov_shadow_sqrtm*randn(I), falses(I)) # LoS not used in this model
    end
end

function draw_channel{MS_t <: PhysicalMS, BS_t <: PhysicalBS}(network::TriangularHetNetNetwork{MS_t, BS_t, SinglecarrierSystem, SimpleLargescaleFadingEnvironment})
    K = get_num_MSs(network); I = get_num_BSs(network)
    Ns = get_num_MS_antennas(network); Ms = get_num_BS_antennas(network)

    coefs = Array(Matrix{Complex128}, K, I)
    large_scale_fading_factor = Array(Float64, K, I)

    distances = get_distances(network)
    angles = get_angles(network)

    for k = 1:K
        for i = 1:I
            # Small scale fading
            coefs[k,i] = (1/sqrt(2))*(randn(Ns[k], Ms[i])
                                 + im*randn(Ns[k], Ms[i]))

            # Pathloss
            pathloss_factor = sqrt(10^(-(network.propagation_environment.pathloss_beta + network.propagation_environment.pathloss_alpha*log10(distances[k,i]))/10))

            # Shadow fading
            shadow_factor = sqrt(10^(network.MSs[k].propagation_environment_state.shadow_realizations_dB[i]/10))

            # Penetration loss
            penetration_loss_factor = sqrt(10^(-network.propagation_environment.penetration_loss_dB/10))

            # BS antenna gain
            bs_antenna_gain = sqrt(get_antenna_gain(network.BSs[i].antenna_params, angles[k,i]))

            # MS antenna gain
            ms_antenna_gain = sqrt(10^(network.MSs[k].antenna_gain_dB/10))

            # Apply scale factors
            large_scale_fading_factor[k,i] = pathloss_factor*shadow_factor*penetration_loss_factor*bs_antenna_gain*ms_antenna_gain
            coefs[k,i] *= large_scale_fading_factor[k,i]
        end
    end

    return SinglecarrierChannel(coefs, Ns, Ms, K, I, large_scale_fading_factor)
end

##########################################################################
# Visualization functions
function plot_network_layout(network::TriangularHetNetNetwork)
    K = get_num_MSs(network)

    fig = PyPlot.figure()
    ax = fig[:add_subplot](1, 1, 1)

    # Equilateral triangle, with cell edges
    cm = mean([ [network.BSs[i].position.x, network.BSs[i].position.y] for i = 1:3 ])
    for i = 1:3
        pos1 = network.BSs[i].position
        pos2 = network.BSs[mod(i, 3) + 1].position
        ax[:plot]([pos1.x, pos2.x], [pos1.y, pos2.y]; color="b", linestyle="-")
        ax[:plot]([mean([pos1.x, pos2.x]), cm[1]], [mean([pos1.y, pos2.y]), cm[2]]; color="b", linestyle="--")
    end

    # Macro BSs
    for i = 1:3
        pos = network.BSs[i].position
        ax[:plot](pos.x, pos.y; marker="x", color="b", markersize=8)
        ax[:text](pos.x, pos.y, "MacroBS $i")

        # Antenna boresight
        spos = Position(network.inter_site_distance/10, 0); epos = rotate(spos, network.BSs[i].antenna_params.bore_sight_angle)
        ax[:plot](pos.x .+ [0, epos.x], pos.y .+ [0, epos.y], color="r", linestyle=":")
    end

    # Pico BSs
    for i = 4:(3 + 3*network.num_picos_per_cell)
        pos = network.BSs[i].position
        ax[:plot](pos.x, pos.y; marker=".", color="b", markersize=6)
        ax[:text](pos.x, pos.y, "PicoBS $i")
    end

    # MSs
    for k = 1:K
        pos = network.MSs[k].position
        ax[:plot](pos.x, pos.y; marker="o", color="r", markersize=6)
        ax[:text](pos.x, pos.y, "MS $k")
    end

    display(fig)
end


text:
# functions to create input files
using ODLCommonTools

"""
  Make an input file from an options dictionary.

  Inputs:
    dict: the options dictionary
    fname: the file name (without extension)

  Outputs:
    fname_ex: file name with extension
"""
function make_input(dict::Dict, fname::AbstractString)

  # add extension if not already present
  if !endswith(fname, ".jl")
    fname_ex = string(fname, ".jl")
    rmfile(fname_ex)
  else
    fname_ex = fname
  end

  f = open(fname_ex, "w")

  println(f, "arg_dict = Dict{Any,Any}(")

  for i in keys(dict)
    show(f, i)
    print(f, " => ")
    show(f, dict[i])
    println(f, ",")
  end

  println(f, ")")
  close(f)

  return fname_ex
end