text:
export SplitExplicitLSRK3nSolver

using KernelAbstractions
using KernelAbstractions.Extras: @unroll
using StaticArrays
using ...SystemSolvers
using ...MPIStateArrays: array_device, realview
using ...GenericCallbacks
#using Printf

using ...ODESolvers:
    AbstractODESolver, LowStorageRungeKutta3N, update!, updatedt!, getdt
import ...ODESolvers: dostep!

using ...BalanceLaws:
    tendency_from_slow_to_fast!,
    cummulate_fast_solution!,
    reconcile_from_fast_to_slow!

LSRK3N = LowStorageRungeKutta3N

@doc """
    SplitExplicitLSRK3nSolver(slow_solver, fast_solver; dt, t0 = 0, coupled = true)

This is a time stepping object for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q_fast} = f_fast(Q_fast, Q_slow, t)
  \\dot{Q_slow} = f_slow(Q_slow, Q_fast, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

This method performs an operator splitting to timestep the Sea-Surface elevation
and vertically averaged horizontal velocity of the model at a faster rate than
the full model, using LowStorageRungeKutta3N time-stepping.

""" SplitExplicitLSRK3nSolver
mutable struct SplitExplicitLSRK3nSolver{SS, FS, RT, MSA, MSB} <:
               AbstractODESolver
    "slow solver"
    slow_solver::SS
    "fast solver"
    fast_solver::FS
    "time step"
    dt::RT
    "time"
    t::RT
    "elapsed time steps"
    steps::Int
    "storage for transfer tendency"
    dQ2fast::MSA
    "saving original fast state"
    S_fast::MSB

    function SplitExplicitLSRK3nSolver(
        slow_solver::LSRK3N,
        fast_solver,
        Q = nothing;
        dt = getdt(slow_solver),
        t0 = slow_solver.t,
    ) where {AT <: AbstractArray}
        SS = typeof(slow_solver)
        FS = typeof(fast_solver)
        RT = real(eltype(slow_solver.dQ))

        dQ2fast = similar(slow_solver.dQ)
        dQ2fast .= -0.0
        #- Warning: Number of fast-solution to save (here only 1, in S_fast) should be as
        #   large as number of non zero Butcher Coeff (including Weight "b") below 2 diagonal,
        #   i.e., with ns= Number of Stages and a(ns+1,:) = b(:), all non zero a(i,j)|_{i > j+1}.
        #  Saving only 1 fast-solution workd for LS3NRK33Heuns.
        S_fast = similar(fast_solver.dQ)
        S_fast .= -0.0
        MSA = typeof(dQ2fast)
        MSB = typeof(S_fast)
        return new{SS, FS, RT, MSA, MSB}(
            slow_solver,
            fast_solver,
            RT(dt),
            RT(t0),
            0,
            dQ2fast,
            S_fast,
        )
    end
end

function dostep!(
    Qvec,
    split::SplitExplicitLSRK3nSolver{SS},
    param,
    time::Real,
) where {SS <: LSRK3N}
    slow = split.slow_solver
    fast = split.fast_solver

    Qslow = Qvec.slow
    Qfast = Qvec.fast

    dQslow = slow.dQ
    dRslow = slow.dR
    dQ2fast = split.dQ2fast
    S_fast = split.S_fast

    slow_bl = slow.rhs!.balance_law
    fast_bl = fast.rhs!.balance_law

    rv_Q = realview(Qslow)
    rv_dQ = realview(dQslow)
    rv_dR = realview(dRslow)
    groupsize = 256

    slow_dt = getdt(slow)
    fast_dt_in = getdt(fast)
    rkA = slow.RKA
    rkB = slow.RKB
    rkC = slow.RKC
    rkW = slow.RKW
    nStages = length(rkC)

    rv_dR .= -0
    for s in 1:nStages
        # Current slow state time
        slow_stage_time = time + rkC[s] * slow_dt
        # @printf("-- main dostep! stage s=%3i , t= %10.2f\n",s,slow_stage_time)

        # Initialize fast model and set time-step and number of substeps we need
        fast_steps = [0 0 0]
        FT = typeof(slow_dt)
        fast_time_rec = [fast_dt_in FT(0) FT(0)]
        set_fast_for_stepping!(
            slow_bl,
            fast_bl,
            fast.rhs!,
            Qfast,
            S_fast,
            slow_dt,
            rkC,
            rkW,
            s,
            nStages,
            fast_time_rec,
            fast_steps,
        )
        # Initialize tentency adjustment before evaluation of slow mode
        initialize_adjustment!(
            slow_bl,
            fast_bl,
            slow.rhs!,
            fast.rhs!,
            Qslow,
            Qfast,
        )

        # Evaluate the slow mode
        # --> save tendency for the fast
        slow.rhs!(dQ2fast, Qslow, param, slow_stage_time, increment = false)

        # vertically integrate slow tendency to advance fast equation
        # and use vertical mean for slow model (negative source)
        # ---> work with dQ2fast as input
        tendency_from_slow_to_fast!(
            slow_bl,
            fast_bl,
            slow.rhs!,
            fast.rhs!,
            Qslow,
            Qfast,
            dQ2fast,
        )

        # Compute (and RK update) slow tendency
        slow.rhs!(dQslow, Qslow, param, slow_stage_time, increment = true)

        # Update (RK-stage) slow state
        event = Event(array_device(Qslow))
        event = update!(array_device(Qslow), groupsize)(
            rv_dQ,
            rv_dR,
            rv_Q,
            rkA[s % nStages + 1, 1],
            rkA[s % nStages + 1, 2],
            rkB[s, 1],
            rkB[s, 2],
            slow_dt,
            nothing,
            nothing,
            nothing;
            ndrange = length(rv_Q),
            dependencies = (event,),
        )
        wait(array_device(Qslow), event)

        # Determine number of substeps we need
        fast_dt = fast_time_rec[1]
        nsubsteps = fast_steps[3]

        updatedt!(fast, fast_dt)

        for substep in 1:nsubsteps
            fast_time = time + fast_time_rec[3] + (substep - 1) * fast_dt
            # @printf("-- main dostep! substep=%3i , t= %10.2f\n",substep,fast_time)
            dostep!(Qfast, fast, param, fast_time)

            # cumulate fast solution
            cummulate_fast_solution!(
                fast_bl,
                fast.rhs!,
                Qfast,
                fast_time,
                fast_dt,
                substep,
                fast_steps,
                fast_time_rec,
            )
        end

        # reconcile slow equation using fast equation
        reconcile_from_fast_to_slow!(
            slow_bl,
            fast_bl,
            slow.rhs!,
            fast.rhs!,
            Qslow,
            Qfast,
            fast_time_rec;
            lastStage = (s == nStages),
        )

    end

    # reset fast time-step to original value
    updatedt!(fast, fast_dt_in)

    # now do implicit mixing step
    nImplSteps = slow_bl.numImplSteps
    if nImplSteps > 0
        # 1. get implicit mising model, model state variable array and solver handles
        ivdc_dg = slow.rhs!.modeldata.ivdc_dg
        ivdc_bl = ivdc_dg.balance_law
        ivdc_Q = slow.rhs!.modeldata.ivdc_Q
        ivdc_solver = slow.rhs!.modeldata.ivdc_bgm_solver
        # ivdc_solver_dt = getdt(ivdc_solver) # would work if solver time-step was set
        # FT = typeof(slow_dt)
        # ivdc_solver_dt = slow_dt / FT(nImplSteps) # just recompute time-step
        ivdc_solver_dt = ivdc_bl.parent_om.ivdc_dt
        # println("ivdc_solver_dt = ",ivdc_solver_dt )
        # 2. setup start RHS, initial guess and values for computing mixing coeff
        ivdc_Q.θ .= Qslow.θ
        ivdc_RHS = slow.rhs!.modeldata.ivdc_RHS
        ivdc_RHS.θ .= Qslow.θ
        ivdc_RHS.θ .= ivdc_RHS.θ ./ ivdc_solver_dt
        ivdc_dg.state_auxiliary.θ_init .= ivdc_Q.θ
        # 3. Invoke iterative solver

        println("BEFORE maximum(ivdc_Q.θ[:]): ", maximum(ivdc_Q.realdata[:]))
        println("BEFORE minimum(ivdc_Q.θ[:]): ", minimum(ivdc_Q.realdata[:]))

        lm!(y, x) = ivdc_dg(y, x, nothing, 0; increment = false)
        solve_tot = 0
        iter_tot = 0
        for i in 1:nImplSteps
            solve_time = @elapsed iters =
                linearsolve!(lm!, nothing, ivdc_solver, ivdc_Q, ivdc_RHS)
            solve_tot = solve_tot + solve_time
            iter_tot = iter_tot + iters
            # Set new RHS and initial values
            ivdc_RHS.θ .= ivdc_Q.θ ./ ivdc_solver_dt
            ivdc_dg.state_auxiliary.θ_init .= ivdc_Q.θ
        end
        println("solver iters, time: ", iter_tot, ", ", solve_tot)

        println("AFTER  maximum(ivdc_Q.θ[:]): ", maximum(ivdc_Q.realdata[:]))
        println("AFTER  minimum(ivdc_Q.θ[:]): ", minimum(ivdc_Q.realdata[:]))

        # exit()
        # Now update
        Qslow.θ .= ivdc_Q.θ

    end

    return nothing
end


text:
@system ThermalTime begin
    T: temperature ~ track(u"°C", override)
    Δt(context.clock.step): timestep ~ preserve(u"hr")
    ΔT: magnitude ~ track
    r(ΔT, Δt): rate => ΔT / Δt ~ track(u"hr^-1")
end

@system GrowingDegree(ThermalTime) begin
    Tb: base_temperature ~ preserve(u"°C", extern, parameter)
    To: optimal_temperature ~ preserve(u"°C", optional, extern, parameter)
    Tx: maximum_temperature ~ preserve(u"°C", optional, extern, parameter)

    ΔT(T, Tb, To, Tx): magnitude => begin
        T = !isnothing(To) ? min(T, To) : T
        T = !isnothing(Tx) && T >= Tx ? Tb : T
        T - Tb
    end ~ track(u"K", min=0)
    r(ΔT, Δt): rate => ΔT / Δt ~ track(u"K/hr")
end

@system BetaFunction(ThermalTime) begin
    Tn: minimum_temperature => 0 ~ preserve(u"°C", extern, parameter)
    To: optimal_temperature ~ preserve(u"°C", extern, parameter)
    Tx: maximum_temperature ~ preserve(u"°C", extern, parameter)
    β: beta => 1 ~ preserve(parameter)

    ΔT(T, Tn, To, Tx, β): magnitude => begin
        # beta function, See Yin et al. (1995), Ag For Meteorol., Yan and Hunt (1999) AnnBot, SK
        if (Tn < T < Tx) && (Tn < To < Tx)
            Ton = To - Tn
            Txo = Tx - To
            f = (T - Tn) / Ton
            g = (Tx - T) / Txo
            α = β * (Ton / Txo)
            f^α * g^β
        else
            0.
        end
    end ~ track
end

@system Q10Function(ThermalTime) begin
    To: optimal_temperature ~ preserve(u"°C", extern, parameter)
    Q10 => 2 ~ preserve(extern, parameter)

    #FIXME: Q10 isn't actually a thermal function like others (not a rate, check unit)
    ΔT(T, To, Q10): magnitude => begin
        Q10^((T - To) / 10u"K")
    end ~ track
end

export GrowingDegree, BetaFunction, Q10Function

ext: jl


text:
#Capacity Market Clearing Module

# Capacity demand curve creator
function create_capacity_demand_curve(
    input_file::String,
    system_peak_load::Float64,
)

# Gather parameter data
capacity_demand_params = read_data(input_file)[1, :]
fpr = capacity_demand_params["Forecast Pool Req"] #Forecast Pool Req
rel_req = system_peak_load * fpr  # Reliability Requirement
irm = capacity_demand_params["IRM"] # Installed Reserve Margin
irm_perc_points = parse.(Float64, split(capacity_demand_params["IRM perc points"], ";")) #Installed Reserve Margin percentage points

net_CONE = capacity_demand_params["Net CONE per day"] * 365  # Net CONE
net_CONE_perc_points = parse.(Float64, split(capacity_demand_params["Net CONE perc points"], ";")) #Net CONE percentage points

@assert length(net_CONE_perc_points) == length(irm_perc_points)

gross_CONE_value = capacity_demand_params["Gross CONE per day"] * 365   # Gross CONE
gross_CONE_points = zeros(length(net_CONE_perc_points))
gross_CONE_points[1] = gross_CONE_value

max_clear = capacity_demand_params["Max Clear"]

# Construct demand curve
demand_curve_break_points = zeros(length(net_CONE_perc_points) + 2)
demand_curve_price_points = zeros(length(net_CONE_perc_points) + 2)

demand_curve_break_points[1] = 0.0
demand_curve_price_points[1] = net_CONE*net_CONE_perc_points[1]

demand_curve_break_points[end] = rel_req * max_clear
demand_curve_price_points[end] = 0.0

# Demand curve points based on PJM capacity market
for i in 1:(length(net_CONE_perc_points))
demand_curve_break_points[i + 1] = rel_req * ( 1 + irm + irm_perc_points[i]) / ( 1 + irm)
demand_curve_price_points[i + 1] = max(net_CONE * net_CONE_perc_points[i], gross_CONE_points[i])
end

capacity_demand_curve = CapacityMarket(demand_curve_break_points, demand_curve_price_points)
return capacity_demand_curve
end

"""
#Construct VRR Curve
VRR_curve_points = [[(rel_req*(1+irm+irm_perc_points[i])/(1+irm)+ee_add_back),
max(net_CONE*net_CONE_perc_points[i],gross_CONE[i])] for i =1:length(irm_perc_points)];
pushfirst!(VRR_curve_points,[0, max(net_CONE*net_CONE_perc_points[1],gross_CONE[1])])
push!(VRR_curve_points,[rel_req+25000, 0])
"""
function capacity_market_clearing(
    demand_curve::CapacityMarket,
    supply_curve::Vector{Vector{Union{String, Float64}}},
    optimizer::JuMP.MOI.OptimizerWithAttributes,
)

"""
Each Element of the Supply Curve:
[1] - Project Name
[2] - Project Size
[3] - Project Capacity Bid
[4] - Project Derating Factor
"""

demand_segmentsize, demand_segmentgrad, demand_pricepoints = makecapacitydemand(demand_curve)

#Number of segments of capacity supply and demand curves
n_demand_seg = length(demand_segmentsize);
n_supply_seg = length(supply_curve);

#----------Capacity Market Clearing Problem---------------------------------------------------------------------------------
cap_mkt = JuMP.Model(optimizer);

#Define the variables
JuMP.@variables(cap_mkt, begin
Q_supply[s=1:n_supply_seg] >= 0 # Quantity of cleared capacity supply offers
Q_demand[d=1:n_demand_seg] >= 0 # Quantity of cleared capacity demand from the demand curve
end)

#Functions----------------------------------------------------------------------------------------

#Cost of procuring capacity supply for each segment
supply_cost(Q_seg, s) = supply_curve[s][3] * Q_seg;

#Welfare from meeting the capacity resource requirement for each segment
demand_welfare(Q_seg, d) = demand_pricepoints[d] * Q_seg + 0.5 * demand_segmentgrad[d] * (Q_seg^2);

#Expressions---------------------------------------------------------------------------------------

#Totoal Cost of procuring capacity supply
JuMP.@expression(cap_mkt, total_supply_cost, sum(supply_cost(Q_supply[s], s) for s = 1:n_supply_seg))

#Total Welfare from meeting the capacity resource requirement
JuMP.@expression(cap_mkt, total_dem_welfare, sum(demand_welfare(Q_demand[d], d) for d = 1:n_demand_seg))

#Constraints--------------------------------------------------------------------------------------

#Cleared capacity supply limit for each segment
JuMP.@constraint(cap_mkt, [s=1:n_supply_seg], Q_supply[s] <= supply_curve[s][2])

#Cleared demand limit for each segment
JuMP.@constraint(cap_mkt, [d=1:n_demand_seg], Q_demand[d] <= demand_segmentsize[d])

#Total cleared capacity supply should meet the total cleared demand
JuMP.@constraint(cap_mkt, mkt_clear, sum(Q_supply[s] * supply_curve[s][4] for s in 1:n_supply_seg) == sum(Q_demand))

#Define Objective Function - Social Welfare Maxmization
JuMP.@objective(cap_mkt, Max, total_dem_welfare - total_supply_cost)

println("Actual Capacity Market Clearing:")
JuMP.optimize!(cap_mkt)
println(JuMP.termination_status(cap_mkt))
println(JuMP.objective_value(cap_mkt))

#Capacity Market Clearing Price is the shadow variable of the capacity balance constraint
cap_price = JuMP.dual(mkt_clear)
cap_accepted_bid = Dict(supply_curve[s][1] => JuMP.value.(Q_supply[s]) / supply_curve[s][2]  for s in 1:n_supply_seg)
#------------------------------------------------------------------------------------------------
return cap_price, cap_accepted_bid

end


text:
#
# test the linear program solver on a simple problem
#

include("../src/homogeneous_algorithm/homogeneous_algorithm.jl")
include_print("testing_tools.jl")

#simple_tests();

file_name = "Problems/QAP8.mat";

A, b, c = get_netlib_problem(file_name);
println("Solving ", file_name, " with the homogeneous algorithm")
println(size(A,2), " variables and ", size(A,1), " constraints")
println("Non-zeros: ", length(nonzeros(A)))
settings = class_settings();

println("=================== Linear system solver is julia LU (normal system) =============")
settings.linear_system_solver = linear_solver_JULIA();
settings.linear_system_solver.options.sym = 0;
test_homogeneous_algorithm(A, b, c, settings);

println("=================== Linear system solver is julia LU factor =============")
settings.linear_system_solver = linear_solver_julia();
test_homogeneous_algorithm(A, b, c, settings);

println("=================== Linear system solver is MUMPS LU ====================")
settings.linear_system_solver = linear_solver_MUMPS();
settings.linear_system_solver.options.sym = 0;
test_homogeneous_algorithm(A, b, c, settings);

println("=================== Linear system solver is MUMPS LDL ===================")
#settings.linear_system_solver = linear_solver_MUMPS();
#settings.linear_system_solver.options.sym = 2;
#test_homogeneous_algorithm(A, b, c, settings);

println("=================== Linear system solver is MUMPS chol (normal system) ==================")
settings.linear_system_solver = linear_solver_MUMPS();
settings.linear_system_solver.options.normal = true;
settings.linear_system_solver.options.sym = 1;
test_homogeneous_algorithm(A, b, c, settings);

println("=================== Linear system solver is matlab ldl ==================")
settings.linear_system_solver = linear_solver_matlab();
settings.linear_system_solver.options.normal = false;
settings.linear_system_solver.options.sym = 2;
test_homogeneous_algorithm(A, b, c, settings);

using Mosek
println("=========================================================================")
println("Calls Mosek")
solve_with_JuMP(A, b, c, MosekSolver());

using Ipopt
println("=========================================================================")
println("Calls IPOPT")
solve_with_JuMP(A, b, c, IpoptSolver());

max_forks_repo_name: ohinder/Everest


text:
#
# Aliases for Unitful dimension singletons
#

const MassParameter = MassParameter
const Velocity = Unitful.Velocity
const Time = Unitful.Time
const Mass = Unitful.Mass


text:
# TODO: gravity_j2, qEI, date2Y, all the Julian date stuff
# TODO: igrf12

# Point-mass gravity model
function point_mass_gravity(r, mu)
    return (-mu / norm(r)^3) * r
end


text:
module Halite

@enum DockedStatus UNDOCKED=0 DOCKING=1 DOCKED=2 UNDOCKING=3

module Constants
    #: Max number of units of distance a ship can travel in a turn
    const MAX_SPEED = 7
    #: Radius of a ship
    const SHIP_RADIUS = 0.5
    #: Starting health of ship, also its max
    const MAX_SHIP_HEALTH = 255
    #: Starting health of ship, also its max
    const BASE_SHIP_HEALTH = 255
    #: Weapon cooldown period
    const WEAPON_COOLDOWN = 1
    #: Weapon damage radius
    const WEAPON_RADIUS = 5.0
    #: Weapon damage
    const WEAPON_DAMAGE = 64
    #: Radius in which explosions affect other entities
    const EXPLOSION_RADIUS = 10.0
    #: Distance from the edge of the planet at which ships can try to dock
    const DOCK_RADIUS = 4.0
    #: Number of turns it takes to dock a ship
    const DOCK_TURNS = 5
    #: Number of production units per turn contributed by each docked ship
    const BASE_PRODUCTIVITY = 6
    #: Distance from the planets edge at which new ships are created
    const SPAWN_RADIUS = 2.0
end

include("entity.jl")
include("player.jl")
include("game_map.jl")
include("game.jl")

export Game, Constants, DockedStatus
export start_game, thrust, dock, undock
export update_map, all_players, all_planets, all_ships, all_docked_ships, get_player, get_me, get_planet, get_docked_ship
export isdocked, isowned, isfull, can_dock
export navigate, send_command_queue, closest_point_to, nearest_unoccupied_planet

end

max_stars_repo_name: xinnosuke/Halite-II

max_issues_repo_name: xinnosuke/Halite-II

max_forks_repo_name: xinnosuke/Halite-II


text:
struct BandedQR{T} <: Factorization{T}
    H::Matrix{T}  # Represents the orthogonal matrix Q
    R::BandedMatrix{T}
end

size(QR::BandedQR,k) = size(QR.R,k)

struct BandedQ{T} <: AbstractMatrix{T}
    H::Matrix{T}
    m::Int   # number of rows/cols
end

size(A::BandedQ) = (A.m, A.m)
size(A::BandedQ,i::Integer) = i <= 0 ? error("dimension out of range") :
                                i == 1 ? A.m :
                                i == 2 ? A.m : 1

transpose(A::BandedQ) = Transpose(A)
adjoint(A::BandedQ) = Adjoint(A)

*(At::Transpose{T,BandedQ{T}}, B::AbstractVecOrMat{T}) where {T<:Real} = transpose(At)'*B
mul!(Y, At::Transpose{T,BandedQ{T}}, B::AbstractVecOrMat{T}) where {T<:Real} =
    mul!(Y, transpose(At)', B)

*(Ac::Adjoint{T,BandedQ{T}}, B::AbstractVector{V}) where {T<:Real,V<:Complex} =
    Ac * real(B) + im* Ac * imag(B)
*(Ac::Adjoint{T,BandedQ{T}}, B::AbstractVector{V}) where {T<:Real,V<:Real} =
    Ac * convert(Vector{T}, B)
*(Ac::Adjoint{T,BandedQ{T}}, B::AbstractVector) where {T<:Complex} =
    Ac * convert(Vector{T}, B)
*(Ac::Adjoint{T,BandedQ{T}}, B::AbstractVector{T}) where {T<:Real} =
    mul!(similar(B), Ac, B)
*(Ac::Adjoint{T,BandedQ{T}}, B::AbstractVector{T}) where {T<:Complex} =
    mul!(similar(B), Ac, B)

function mul!(Y::Vector{T}, Ac::Adjoint{T,BandedQ{T}}, B::Vector{T}) where {T<:BlasFloat}
    A = parent(Ac)
    if length(Y) != size(A,1) || length(B) != size(A,2)
        throw(DimensionMismatch("Matrices have wrong dimensions"))
    end

    H=A.H
    m=A.m

    M=size(H,1)
    b=pointer(B)
    y=pointer(Y)
    h=pointer(H)
    st=stride(H,2)

    sz=sizeof(T)

    BLAS.blascopy!(m,B,1,Y,1)

    for k=1:min(size(H,2),m-M+1)
        wp=h+sz*st*(k-1)
        yp=y+sz*(k-1)

        dt=dot(M,wp,1,yp,1)
        axpy!(M,-2*dt,wp,1,yp,1)
    end

    for k=m-M+2:size(H,2)
        p=k-m+M-1

        wp=h+sz*st*(k-1)
        yp=y+sz*(k-1)

        dt=dot(M-p,wp,1,yp,1)
        axpy!(M-p,-2*dt,wp,1,yp,1)
    end
    Y
end

# Each householder is symmetyric, this just reverses the order of application
function mul!(Y::Vector{T}, A::BandedQ{T}, B::Vector{T}) where {T<:BlasFloat}
    if length(Y) != size(A,1) || length(B) != size(A,2)
        throw(DimensionMismatch("Matrices have wrong dimensions"))
    end

    H=A.H
    m=A.m

    M=size(H,1)
    b=pointer(B)
    y=pointer(Y)
    h=pointer(H)
    st=stride(H,2)

    sz=sizeof(T)

    BLAS.blascopy!(m,B,1,Y,1)

    for k=size(H,2):-1:m-M+2
        p=k-m+M-1

        wp=h+sz*st*(k-1)
        yp=y+sz*(k-1)

        dt=dot(M-p,wp,1,yp,1)
        axpy!(M-p,-2*dt,wp,1,yp,1)
    end

    for k=min(size(H,2),m-M+1):-1:1
        wp=h+sz*st*(k-1)
        yp=y+sz*(k-1)

        dt=dot(M,wp,1,yp,1)
        axpy!(M,-2*dt,wp,1,yp,1)
    end

    Y
end

function mul!(Y::Matrix, A::BandedQ, B::Matrix)
    for j=1:size(A,2)
        Y[:,j]=A*B[:,j]
    end
    Y
end

function mul!(Y::Matrix,A::Adjoint{T,BandedQ{T}},B::Matrix) where T
    for j=1:size(A,2)
        Y[:,j]=A'*B[:,j]
    end
    Y
end

function (*)(A::BandedQ{T}, x::AbstractVector{S}) where {T,S}
    TS = promote_type(T, S)
    mul!(similar(x,TS,size(A,1)),A,x)
end

Base.IndexStyle(::Type{BandedQ{T}}) where {T} = IndexCartesian()
Base.getindex(A::BandedQ,k::Int,j::Int) = (A*eltype(A)[zeros(j-1);1.0;zeros(size(A,2)-j)])[k]

function Base.getindex(QR::BandedQR,d::Symbol)
    d == :Q && return BandedQ(QR.H,size(QR,1))
    d == :R && return QR.R

    throw(KeyError(d))
end

function getproperty(QR::BandedQR, d::Symbol)
    if d == :R
        return getfield(QR,:R)
    elseif d == :Q
        return BandedQ(QR.H,size(QR,1))
    else
        getfield(QR, d)
    end
end
# iteration for destructuring into components
iterate(S::BandedQR) = (S.Q, Val(:R))
iterate(S::BandedQR, ::Val{:R}) = (S.R, Val(:done))
iterate(S::BandedQR, ::Val{:done}) = nothing

function qr(A::BandedMatrix)
    R=BandedMatrix(Zeros{eltype(A)}(size(A)), (A.l,A.l+A.u))
    R.data[A.l+1:end,:]=A.data
    banded_qr!(R)
end

flipsign(x,y) = Base.flipsign(x,y)
flipsign(x::BigFloat,y::BigFloat) = sign(y)==1 ? x : (-x)
flipsign(x,y::Complex) = y==0 ? x : x*sign(y)

function banded_qr!(R::BandedMatrix{T}) where T
    M=R.l+1   # number of diag+subdiagonal bands
    m,n=size(R)
    W=Matrix{T}(undef, M, (n<m ? n : m-1))
    w=pointer(W)
    r=pointer(R.data)
    sz=sizeof(T)
    st=stride(R.data,2)
    stw=stride(W,2)

    for k=1:min(size(R,1)-R.l,n)
        v=r+sz*(R.u + (k-1)*st)    # diagonal entry
        wp=w+stw*sz*(k-1)          # k-th column of W
        BLAS.blascopy!(M,v,1,wp,1)
        W[1,k]+= flipsign(BLAS.nrm2(M,wp,1),W[1,k])
        normalize!(M,wp)

        for j=k:min(k+R.u,n)
            v=r+sz*(R.u + (k-1)*st + (j-k)*(st-1))
            dt=dot(M,wp,1,v,1)
            axpy!(M,-2*dt,wp,1,v,1)
        end
    end

    for k=m-R.l+1:(n<m ? n : m-1)
        p=k-m+R.l
        v=r+sz*(R.u + (k-1)*st)    # diagonal entry
        wp=w+stw*sz*(k-1)          # k-th column of W
        BLAS.blascopy!(M-p,v,1,wp,1)
        W[1,k]+= flipsign(BLAS.nrm2(M-p,wp,1),W[1,k])
        normalize!(M-p,wp)

        for j=k:min(k+R.u,n)
            v=r+sz*(R.u + (k-1)*st + (j-k)*(st-1))
            dt=dot(M-p,wp,1,v,1)
            axpy!(M-p,-2*dt,wp,1,v,1)
        end
    end

    BandedQR(W, _BandedMatrix(R.data[1:R.u+1,:],m,0,R.u))
end


text:
module Sass

const deps_file = normpath(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
if isfile(deps_file)
    include(deps_file)
else
    error("Sass not properly installed. Please run `] build Sass` in the REPL")
end
include("version.jl")
include("context.jl")
include("options.jl")
include("compiler.jl")
include("julian_api.jl")

examplefolder = normpath(joinpath(@__DIR__, "..", "test", "css"))

end # module


text:
function solve(
    i::SpanningTreeInstance{T, Maximise},
    ::DynamicProgramming,
) where {T}
    # A dynamic-programming stance on Prim's algorithm.

    # Get an indexed list of *all* edges, sorted by increasing value.
    sorted_edges = Pair{Edge, Float64}[]
    for e in edges(i.graph)
        push!(sorted_edges, e => reward(i, e))
    end
    sort!(sorted_edges, by=(p) -> -p[2])

    # Fill the DP table, considering edges one by one (from the best to the
    # poorest), Prim-like.
    # Index: considering edges up to that index (included).
    V = Vector{Float64}(undef, ne(i.graph))
    S = Vector{Vector{Edge{T}}}(undef, ne(i.graph))

    first_edge, first_edge_value = sorted_edges[1]
    V[1] = first_edge_value
    S[1] = Edge{T}[first_edge]

    ld = LoopDetector(nv(i.graph))
    visit_edge(ld, src(first_edge), dst(first_edge))

    # Dynamic part.
    for i in 2:ne(i.graph)
        # Can edges_sorted[i] be taken?
        edge, edge_value = sorted_edges[i]

        if edge_would_create_loop(ld, src(edge), dst(edge)) # Don't take i.
            V[i] = V[i - 1]
            S[i] = S[i - 1]
        else # Take i.
            V[i] = V[i - 1] + edge_value
            S[i] = copy(S[i - 1])
            push!(S[i], edge)

            visit_edge(ld, src(edge), dst(edge))
        end
    end

    return SpanningTreeSolution(i, S[end])
end


text:
using Base.Broadcast: broadcastable, Broadcasted
using Swizzles.ExtrudedArrays
using Swizzles.ArrayifiedArrays

foo(x) = x

@testset "ExtrudedArrays" begin
    for arg in ([1, 2, 3.0], [1, 2, 3], [], [1 2; 3 4], transpose([1, 2]))

        x = foo.(arg)
        y = foo.(ExtrudedArray(arg))
        @test x == y
        @test typeof(x) == typeof(y)

        x = arg .+ arg
        y = ExtrudedArray(arg) .+ arg
        @test x == y
        @test typeof(x) == typeof(y)

        x = arg .+ arg
        y = ExtrudedArray(arg) .+ ExtrudedArray(arg)
        @test x == y
        @test typeof(x) == typeof(y)
    end

    bc  = Broadcasted(+, ((1, 2, 3), identity.([1, 2, 3]')))
    bc′ = Broadcasted(+, ((1, 2, 3), ExtrudedArray(identity.([1, 2, 3]'), (extrude, keep))))
    @test typeof(lift_keeps(bc)) == typeof(bc′)
    @test keeps(bc) == (keep, true)

    bc  = Broadcasted(+, ((1, 2, 3), identity.([1]')))
    bc′ = Broadcasted(+, ((1, 2, 3), ExtrudedArray(identity.([1]'), (extrude, extrude))))
    @test typeof(lift_keeps(bc)) == typeof(bc′)
    @test keeps(bc) == (keep, false)

    bc  = Broadcasted(+, ((1, 2, 3), [1, 2, 3]))
    bc′ = Broadcasted(+, ((1, 2, 3), ExtrudedArray([1, 2, 3], (keep,))))
    @test typeof(lift_keeps(bc)) == typeof(bc′)
    @test keeps(bc) == (keep,)

    bc  = Broadcasted(+, ((1,), [1, 2, 3]))
    bc′ = Broadcasted(+, ((1,), ExtrudedArray([1, 2, 3], (keep,))))
    @test typeof(lift_keeps(bc)) == typeof(bc′)
    @test keeps(bc) == (true,)

    bc  = Broadcasted(+, (1, [1, 2, 3]))
    bc′ = Broadcasted(+, (1, ExtrudedArray([1, 2, 3], (keep,))))
    @test typeof(lift_keeps(bc)) == typeof(bc′)
    @test keeps(bc) == (true,)

    bc  = Broadcasted(+, ((1, 2, 3), [1]))
    bc′ = Broadcasted(+, ((1, 2, 3), ExtrudedArray([1], (extrude,))))
    @test typeof(lift_keeps(bc)) == typeof(bc′)
    @test keeps(bc) == (keep,)

    bc  = Broadcasted(+, ((1, 2, 3), 1))
    bc′ = Broadcasted(+, ((1, 2, 3), 1))
    @test typeof(lift_keeps(bc)) == typeof(bc′)
    @test keeps(bc) == (keep,)

    bc  = Broadcasted(+, ((1,), 1))
    bc′ = Broadcasted(+, ((1,), 1))
    @test typeof(lift_keeps(bc)) == typeof(bc′)
    @test keeps(bc) == (extrude,)

    bc  = Delay().(Swizzle(+, (2, 1)).(Broadcasted(+, ((1, 2, 3), identity.([1, 2, 3]')))))
    bc′ = Delay().(Swizzle(+, (2, 1)).(Broadcasted(+, ((1, 2, 3), ExtrudedArray(identity.([1, 2, 3]'), (extrude, keep))))))
    @test typeof(lift_keeps(bc)) == typeof(bc′)
    @test keeps(bc) == (true, keep)

    bc  = Delay().(Swizzle(+, (2, 1)).(Broadcasted(+, ((1,), identity.([1, 2, 3]')))))
    bc′ = Delay().(Swizzle(+, (2, 1)).(Broadcasted(+, ((1,), ExtrudedArray(identity.([1, 2, 3]'), (extrude, keep))))))
    @test typeof(lift_keeps(bc)) == typeof(bc′)
    @test keeps(bc) == (true, false)

    a  = [1, 2, 3]'
    a′ = ExtrudedArray([1, 2, 3], (keep,))'
    @test typeof(lift_keeps(a)) == typeof(a′)
    @test keeps(a) == (extrude, true)
end


text:
import Base: getindex, lastindex

@doc doc"""
    load_chains(filename; burnin=0, burnin_chunks=0, thin=1, join=false, unbatch=true)

Load a single chain or multiple parallel chains which were written to
a file by [`sample_joint`](@ref).

Keyword arguments:

* `burnin` — Remove this many samples from the start of each chain, or
  if negative, keep only this many samples at the end of each chain.
* `burnin_chunks` — Same as burnin, but in terms of chain "chunks"
  stored in the chain file, rather than in terms of samples.
* `thin` — If `thin` is an integer, thin the chain by this factor. If
  `thin == :hasmaps`, return only samples which have maps saved. If
  thin is a `Function`, filter the chain by this function (e.g.
  `thin=haskey(:g)` on Julia 1.5+)
* `unbatch` — If true, [unbatch](@ref) the chains if they are batched.
* `join` — If true, concatenate all the chains together.
* `skip_missing_chunks` — Skip missing chunks in the chain instead of
  terminating the chain there.

The object returned by this function is a `Chain` or `Chains` object,
which simply wraps an `Array` of `Dicts` or an `Array` of `Array` of
`Dicts`, respectively (each sample is a `Dict`). The wrapper object
has some extra indexing properties for convenience:

* It can be indexed as if it were a single multidimensional object,
  e.g. `chains[1,:,:accept]` would return the `:accept` key of all
  samples in the first chain.
* Leading colons can be dropped, i.e. `chains[:,:,:accept]` is the
  same as `chains[:accept]`.
* If some samples are missing a particular key, `missing` is returned
  for those samples insted of an error.
* The recursion goes arbitrarily deep into the objects it finds. E.g.,
  since sampled parameters are stored in a `NamedTuple` like
  `(Aϕ=1.3,)` in the `θ` key of each sample `Dict`, you can do
  `chain[:θ,:Aϕ]` to get all `Aϕ` samples as a vector.

"""
function load_chains(filename; burnin=0, thin=1, join=false, unbatch=true, dropmaps=false, burnin_chunks=0)
    chains = jldopen(filename) do io
        ks = keys(io)
        chunk_ks = sort([k for k in ks if startswith(k,"chunks_")], by=k->parse(Int,k[8:end]))
        chunk_ks = chunk_ks[burnin_chunks>=0 ? (burnin_chunks+1:end) : (end+burnin_chunks+1:end)]
        @showprogress for (isfirst,k) in flagfirst(chunk_ks)
            if isfirst
                chains = read(io,k)
            else
                append!.(chains, read(io,k))
            end
            if dropmaps
                map!(chain -> filter!.(((_,v),)->!(v isa Field), chain), chains, chains)
            end
        end
        chains
    end
    if thin isa Int
        chains = [chain[burnin>=0 ? ((1+burnin):thin:end) : (end+(1+burnin):thin:end)] for chain in chains]
    elseif thin == :hasmaps
        chains = [[samp for samp in chain[(1+burnin):end] if :ϕ in keys(samp)] for chain in chains]
    elseif thin isa Function
        chains = [filter(thin,chain) for chain in chains]
    else
        error("`thin` should be an Int, :hasmaps, or a filter function")
    end
    chains = wrap_chains(chains)
    if unbatch
        chains = CMBLensing.unbatch(chains)
    end
    return join ? Chain(reduce(vcat, chains)) : chains
end

# multiple parallel chains
struct Chains{T} <: AbstractVector{T}
    chains :: Vector{T}
end
Chains(chains::Chains) = chains
getindex(c::Chains, k::Symbol, ks...) = _getindex(c.chains, :, :, k, ks...)
getindex(c::Chains, k,         ks...) = _getindex(c.chains, k, ks...)
size(c::Chains) = size(c.chains)
function Base.print_array(io::IO, cs::Chains; indent="  ")
    for c in cs
        print(io, "  ")
        Base.summary(io, c)
        println(io)
        Base.print_array(io, c, indent="    ")
    end
end
function lastindex(c::Chains, d)
    d==1 ? lastindex(c.chains) :
    d==2 ? lastindex(c.chains[1]) :
    error("`end` only valid in first or second dim of Chains")
end

# single chain
struct Chain{T} <: AbstractVector{T}
    chain :: Vector{T}
end
getindex(c::Chain, k::Symbol,              ks...) = _getindex(c.chain, :, k, ks...)
getindex(c::Chain, k,                      ks...) = _getindex(c.chain, k, ks...)
getindex(c::Chain, k::Union{Colon,AbstractRange}) = Chain(getindex(c.chain, k))
lastindex(c::Chain) = lastindex(c.chain)
lastindex(c::Chain, d) = d==1 ? lastindex(c.chain) : error("`end` only valid in first dim of Chain")
size(c::Chain) = size(c.chain)
function Base.print_array(io::IO, c::Chain; indent="  ")
    _,cols = displaysize(io)
    for k in keys(c[end])
        str = string("$(indent)$(k) => ", repr(c[k]; context=(:limit => true)))
        println(io, Base._truncate_at_width_or_chars(str, cols))
    end
end

# recurse
_getindex(x, k,                             k2, ks...) = _getindex(_getindex(x, k), k2, ks...)
_getindex(x, k::Union{Colon,AbstractRange}, k2, ks...) = broadcast(y -> _getindex(y, k2, ks...),  x[k])

# base cases
_getindex(x::Union{Dict,NamedTuple}, k::Symbol) = haskey(x,k) ? getindex(x, k) : missing
_getindex(x,                         k) = getindex(x, k)

wrap_chains(chains::Vector{<:Vector}) = Chains(Chain.(chains))
wrap_chains(chain::Vector) = Chain(chain)

# batching
@doc doc"""
    unbatch(chain::Chain)

Convert a chain of batch-length-`D` fields to `D` chains of unbatched fields.
"""
function unbatch(chain::Chain)
    D = batch_length(chain[end][:lnP])
    (D==1) && return chain
    Chains(map(1:D) do I
        Chain(map(chain) do samp
            Dict(map(collect(samp)) do (k,v)
                if v isa Union{BatchedReal,FlatField}
                    k => batch_index(v, I)
                elseif v isa NamedTuple{<:Any, <:NTuple{<:Any,<:Union{BatchedReal,BatchedVal}}}
                    k => map(x -> (x isa Union{BatchedReal,BatchedVal} ? batch_index(x,I) : x), v)
                elseif v isa AbstractArray{<:Union{BatchedReal,BatchedVal}}
                    k => batch_index.(v, I)
                else
                    k => v
                end
            end)
        end)
    end)
end

@doc doc"""
    unbatch(chains::Chains)

Expand each chain in this `Chains` object by unbatching it.
"""
unbatch(chains::Chains) = Chains(mapreduce(unbatch, vcat, chains))

@doc doc"""
    mean_std_and_errors(samples; N_bootstrap=10000)

Get the mean and standard deviation of a set of correlated `samples` from a
chain where the error on the mean and standard deviation is estimated with
bootstrap resampling using the calculated "effective sample size" of the chain.
"""
function mean_std_and_errors(samples; N_bootstrap=10000, N_in_paren=2, tol=50)

    Neff = round(Int, length(samples) / @ondemand(PyCall.pyimport)(:emcee).autocorr.integrated_time(samples; tol)[1])

    μ = mean(samples)
    σ = std(samples)

    SEμ = std([mean(samples[rand(1:end, Neff)]) for i=1:N_bootstrap])
    SEσ = std([ std(samples[rand(1:end, Neff)]) for i=1:N_bootstrap])

    "$(paren_errors(μ, SEμ; N_in_paren=N_in_paren)) ± $(paren_errors(σ, SEσ; N_in_paren=N_in_paren))"

end

@doc doc"""
    paren_errors(μ, σ; N_in_paren=2)

Get a string represntation of `μ ± σ` in "parenthesis" format, e.g. `1.234 ±
0.012` becomes `1.234(12)`.
"""
function paren_errors(μ, σ; N_in_paren=2)
    N = round(Int, floor(log10(1/σ))) + N_in_paren
    fmt = "%.$(N)f"
    @ondemand(Formatting.sprintf1)(fmt, μ)*"($(round(Int,σ*10^N)))"
end

# an N-dimensional getdist Kernel Density Estimate
struct GetDistKDE{N}
    kde
end
(k::GetDistKDE)(x...) = k.kde(x...)

"""
    kde(samples::AbstractVector; [boundary=(min,max), normalize="integral" or "max"])
    kde(samples::AbstractMatrix; [boundary=[(min1,max1),(min2,max2)], normalize="integral" or "max", smooth_scale_2D])

Return a Kernel Density Estimate for a set of 1D or 2D samples. The
return object is a function which can be evaluated anywhere to compute
the PDF. If provided, `boundary` specifies a hard upper/lower bound
for the 1 or 2 or parameters, `normalize` specifies whether to
normalize the PDF to unit integral or unit maximum, and
`smooth_scale_2D` specifies how much smoothing to do for the 2D case.

Based on Python [GetDist](https://getdist.readthedocs.io/en/latest/intro.html),
which must be installed.
"""
function kde(samples::AbstractVector; boundary=(nothing,nothing), normalize="integral")
    getdist = @ondemand(PyCall.pyimport)("getdist")
    getdist.chains.print_load_details = false
    kde = getdist.MCSamples(;
        samples, weights=nothing, names=["x"], ranges=Dict("x"=>boundary)
    )
    GetDistKDE{1}(kde.get1DDensity(0).normalize(normalize))
end

function kde(samples::AbstractMatrix; boundary=((nothing,nothing),(nothing,nothing)), normalize="integral", smooth_scale_2D=nothing)
    if size(samples,1) == 2
        samples = samples'
    elseif size(samples,2) != 2
        error("KDE only supports 1 or 2 dimensional samples.")
    end
    getdist = @ondemand(PyCall.pyimport)("getdist")
    getdist.chains.print_load_details = false
    kde = getdist.MCSamples(;
        samples, weights=nothing, names=["x","y"], ranges=Dict("x"=>boundary[1], "y"=>boundary[2])
    )
    density_kwargs = isnothing(smooth_scale_2D) ? () : (;smooth_scale_2D)
    GetDistKDE{2}(kde.get2DDensity(0, 1; density_kwargs...).normalize(normalize))
end


text:
const BlockSparseMatrix{ElT,StoreT,IndsT} = BlockSparseTensor{ElT,2,StoreT,IndsT}
const DiagBlockSparseMatrix{ElT,StoreT,IndsT} = DiagBlockSparseTensor{ElT,2,StoreT,IndsT}
const DiagMatrix{ElT,StoreT,IndsT} = DiagTensor{ElT,2,StoreT,IndsT}

function _truncated_blockdim(S::DiagMatrix,
                             docut::Float64;
                             singular_values=false,
                             truncate=true)
  !truncate && return diaglength(S)
  newdim = 0
	val = singular_values ? getdiagindex(S,newdim+1)^2 : getdiagindex(S,newdim+1)
  while newdim+1 ≤ diaglength(S) && val > docut
    newdim += 1
    if newdim+1 ≤ diaglength(S)
      val = singular_values ? getdiagindex(S,newdim+1)^2 :
             getdiagindex(S,newdim+1)
    end
  end
  return newdim
end

"""
svd(T::BlockSparseTensor{<:Number,2}; kwargs...)

svd of an order-2 BlockSparseTensor.

This function assumes that there is one block
per row/column, otherwise it fails.
This assumption makes it so the result can be
computed from the dense svds of seperate blocks.
"""
function LinearAlgebra.svd(T::BlockSparseMatrix{ElT};
                           kwargs...) where {ElT}
  alg::String = get(kwargs, :alg, "recursive")

  truncate = haskey(kwargs, :maxdim) || haskey(kwargs, :cutoff)

  Us = Vector{BlockSparseMatrix{ElT}}(undef, nnzblocks(T))
  Ss = Vector{DiagBlockSparseMatrix{real(ElT)}}(undef, nnzblocks(T))
  Vs = Vector{BlockSparseMatrix{ElT}}(undef, nnzblocks(T))

  # Sorted eigenvalues
  d = Vector{real(ElT)}()

  for n in 1:nnzblocks(T)
    b = nzblock(T, n)
    blockT = blockview(T, n)
    Ub, Sb, Vb = svd(blockT; alg = alg)
    Us[n] = Ub
    Ss[n] = Sb
    Vs[n] = Vb
    append!(d, vector(diag(Sb)))
  end

  # Square the singular values to get
  # the eigenvalues
  d .= d .^ 2
  sort!(d; rev=true)

  # Get the list of blocks of T
  # that are not dropped
  nzblocksT = nzblocks(T)

  dropblocks = Int[]
  if truncate
    truncerr,docut = truncate!(d; kwargs...)
    for n in 1:nnzblocks(T)
      blockdim = _truncated_blockdim(Ss[n],
                                     docut;
                                     singular_values=true,
                                     truncate=truncate)
      if blockdim == 0
        push!(dropblocks,n)
      else
        Strunc = tensor(Diag(store(Ss[n])[1:blockdim]),
                        (blockdim,blockdim))
        Ss[n] = Strunc
        Us[n] = copy(Us[n][1:dim(Us[n],1),1:blockdim])
        Vs[n] = copy(Vs[n][1:dim(Vs[n],1),1:blockdim])
      end
    end
    deleteat!(Ss,dropblocks)
    deleteat!(Us,dropblocks)
    deleteat!(Vs,dropblocks)
    deleteat!(nzblocksT,dropblocks)
    else
    truncerr,docut = 0.0,0.0
  end

  # The number of blocks of T remaining
  nnzblocksT = nnzblocks(T)-length(dropblocks)

  #
  # Put the blocks into U,S,V
  #

  nb1_lt_nb2 = (nblocks(T)[1] < nblocks(T)[2] ||
                (nblocks(T)[1] == nblocks(T)[2] &&
                 dim(T,1) < dim(T,2)))

  if nb1_lt_nb2
    uind = sim(ind(T,1))
  else
    uind = sim(ind(T,2))
  end

  deleteat!(uind,dropblocks)

  # uind may have too many blocks
  if nblocks(uind) > nnzblocksT
    resize!(uind,nnzblocksT)
  end

  for n in 1:nnzblocksT
    setblockdim!(uind, minimum(dims(Ss[n])), n)
  end

  if dir(uind) != dir(inds(T)[1])
    uind = dag(uind)
  end
  indsU = setindex(inds(T), dag(uind), 2)

  vind = sim(uind)
  if dir(vind) != dir(inds(T)[2])
    vind = dag(vind)
  end
  indsV = setindex(inds(T),dag(vind),1)
  indsV = permute(indsV,(2,1))

  indsS = setindex(inds(T),uind,1)
  indsS = setindex(indsS,vind,2)

  nzblocksU = Vector{Block{2}}(undef,nnzblocksT)
  nzblocksS = Vector{Block{2}}(undef,nnzblocksT)
  nzblocksV = Vector{Block{2}}(undef,nnzblocksT)

  for n in 1:nnzblocksT
    blockT = nzblocksT[n]

    blockU = (blockT[1],n)
    nzblocksU[n] = blockU

    blockS = (n,n)
    nzblocksS[n] = blockS

    blockV = (blockT[2],n)
    nzblocksV[n] = blockV
  end

  U = BlockSparseTensor(undef,nzblocksU,indsU)
  V = BlockSparseTensor(undef,nzblocksV,indsV)
  S = DiagBlockSparseTensor(undef,nzblocksS,indsS)

  for n in 1:nnzblocksT
    Ub,Sb,Vb = Us[n],Ss[n],Vs[n]

    blockU = nzblocksU[n]
    blockV = nzblocksV[n]
    blockS = nzblocksS[n]

    blockview(U,blockU) .= Ub
    blockview(V,blockV) .= Vb

    blockviewS = blockview(S,blockS)
    for i in 1:diaglength(Sb)
      setdiagindex!(blockviewS,getdiagindex(Sb,i),i)
    end
  end

  return U,S,V,Spectrum(d,truncerr)
end

function LinearAlgebra.eigen(T::Hermitian{ElT,<:BlockSparseMatrix{ElT}};
                             kwargs...) where {ElT<:Union{Real,Complex}}
  truncate = haskey(kwargs,:maxdim) || haskey(kwargs,:cutoff)

  Us = Vector{BlockSparseMatrix{ElT}}(undef,nnzblocks(T))
  Ds = Vector{DiagBlockSparseMatrix{real(ElT)}}(undef,nnzblocks(T))

  # Sorted eigenvalues
  d = Vector{real(ElT)}()

  for n in 1:nnzblocks(T)
    b = nzblock(T, n)
    blockT = blockview(T, n)
    Ub,Db = eigen(blockT)
    Us[n] = Ub
    Ds[n] = Db
    append!(d,vector(diag(Db)))
  end

  dropblocks = Int[]
  sort!(d; rev=true, by=abs)
  if truncate
    truncerr,docut = truncate!(d; kwargs...)
    for n in 1:nnzblocks(T)
      blockdim = _truncated_blockdim(Ds[n],docut)
      if blockdim == 0
        push!(dropblocks,n)
      else
        Dtrunc = tensor(Diag(store(Ds[n])[1:blockdim]),
                        (blockdim,blockdim))
        Ds[n] = Dtrunc
        Us[n] = copy(Us[n][1:dim(Us[n],1),1:blockdim])
      end
    end
    deleteat!(Ds,dropblocks)
    deleteat!(Us,dropblocks)
  else
    truncerr = 0.0
  end

  # Get the list of blocks of T
  # that are not dropped
  nzblocksT = nzblocks(T)
  deleteat!(nzblocksT,dropblocks)

  # The number of blocks of T remaining
  nnzblocksT = nnzblocks(T)-length(dropblocks)

  #
  # Put the blocks into U,D
  #

  nb1_lt_nb2 = (nblocks(T)[1] < nblocks(T)[2] || (nblocks(T)[1] == nblocks(T)[2] && dim(T,1) < dim(T,2)))

  if nb1_lt_nb2
    uind = sim(ind(T,1))
  else
    uind = sim(ind(T,2))
  end

  deleteat!(uind,dropblocks)

  # uind may have too many blocks
  if nblocks(uind) > nnzblocksT
    resize!(uind,nnzblocksT)
  end

  for n in 1:nnzblocksT
    setblockdim!(uind,minimum(dims(Ds[n])),n)
  end

  if dir(uind) != dir(inds(T)[1])
    uind = dag(uind)
  end
  indsU = setindex(inds(T),dag(uind),2)

  vind = sim(uind)
  if dir(vind) != dir(inds(T)[2])
    vind = dag(vind)
  end
  indsV = setindex(inds(T),dag(vind),1)
  indsV = permute(indsV,(2,1))

  indsD = setindex(inds(T),uind,1)
  indsD = setindex(indsD,vind,2)

  nzblocksU = Vector{Block{2}}(undef,nnzblocksT)
  nzblocksD = Vector{Block{2}}(undef,nnzblocksT)

  for n in 1:nnzblocksT
    blockT = nzblocksT[n]

    blockU = (blockT[1],n)
    nzblocksU[n] = blockU

    blockD = (n,n)
    nzblocksD[n] = blockD
  end

  U = BlockSparseTensor(ElT,undef,nzblocksU,indsU)
  D = DiagBlockSparseTensor(ElT,undef,nzblocksD,indsD)

  for n in 1:nnzblocksT
    Ub,Db = Us[n],Ds[n]

    blockU = nzblocksU[n]
    blockD = nzblocksD[n]

    blockview(U,blockU) .= Ub

    blockviewD = blockview(D,blockD)
    for i in 1:diaglength(Db)
      setdiagindex!(blockviewD,getdiagindex(Db,i),i)
    end
  end

  return U,D,Spectrum(d,truncerr)
end

function LinearAlgebra.eigen(T::BlockSparseMatrix{ElT};
                             kwargs...) where {ElT}
  truncate = haskey(kwargs,:maxdim) || haskey(kwargs,:cutoff)

  if truncate
    error("Truncation is not currently supported by non-Hermitian eigendecomposition")
  end

  Us = Vector{BlockSparseMatrix{complex(ElT)}}(undef,nnzblocks(T))
  Ds = Vector{DiagBlockSparseMatrix{complex(ElT)}}(undef,nnzblocks(T))

  # Sorted eigenvalues
  d = Vector{real(ElT)}()

  for n in 1:nnzblocks(T)
    b = nzblock(T, n)
    blockT = blockview(T, n)
    Ub,Db = eigen(blockT)
    Us[n] = complex(Ub)
    Ds[n] = complex(Db)
    append!(d,abs.(vector(diag(Db))))
  end
  sort!(d; rev=true)
  truncerr = 0.0

  # Get the list of blocks of T
  # that are not dropped
  nzblocksT = nzblocks(T)

  # The number of blocks of T remaining
  nnzblocksT = nnzblocks(T)

  #
  # Put the blocks into U,D
  #

  nb1_lt_nb2 = (nblocks(T)[1] < nblocks(T)[2] || (nblocks(T)[1] == nblocks(T)[2] && dim(T,1) < dim(T,2)))

  if nb1_lt_nb2
    uind = sim(ind(T,1))
  else
    uind = sim(ind(T,2))
  end

  # uind may have too many blocks
  if nblocks(uind) > nnzblocksT
    resize!(uind,nnzblocksT)
  end

  for n in 1:nnzblocksT
    setblockdim!(uind,minimum(dims(Ds[n])),n)
  end

  if dir(uind) != dir(inds(T)[1])
    uind = dag(uind)
  end
  indsU = setindex(inds(T),dag(uind),2)

  vind = sim(uind)
  if dir(vind) != dir(inds(T)[2])
    vind = dag(vind)
  end
  indsV = setindex(inds(T),dag(vind),1)
  indsV = permute(indsV,(2,1))

  indsD = setindex(inds(T),uind,1)
  indsD = setindex(indsD,vind,2)

  nzblocksU = Vector{Block{2}}(undef,nnzblocksT)
  nzblocksD = Vector{Block{2}}(undef,nnzblocksT)

  for n in 1:nnzblocksT
    blockT = nzblocksT[n]

    blockU = (blockT[1],n)
    nzblocksU[n] = blockU

    blockD = (n,n)
    nzblocksD[n] = blockD
  end

  U = BlockSparseTensor(complex(ElT),undef,nzblocksU,indsU)
  D = DiagBlockSparseTensor(complex(ElT),undef,nzblocksD,indsD)

  for n in 1:nnzblocksT
    Ub,Db = Us[n],Ds[n]

    blockU = nzblocksU[n]
    blockD = nzblocksD[n]

    blockview(U,blockU) .= Ub

    blockviewD = blockview(D,blockD)
    for i in 1:diaglength(Db)
      setdiagindex!(blockviewD,getdiagindex(Db,i),i)
    end
  end

  return U,D,Spectrum(d,truncerr)
end


text:
module MathOptInterfaceMosek

import MathOptInterface
using Mosek

immutable MosekSolver <: MathOptInterface.AbstractSolver
  options
end

MosekSolver(;kwargs...) = MosekSolver(kwargs)
export MosekSolver

include("LinkedInts.jl")

const DEBUG = false

mosek_block_type_unallocated = 0
mosek_block_type_zero   = 1
mosek_block_type_nonneg = 2
mosek_block_type_nonpos = 3
mosek_block_type_range  = 4
mosek_block_type_qcone  = 5
mosek_block_type_rqcone = 6
mosek_block_type_psd    = 7
mosek_block_type_integer = 8

problemtype_linear    = 0
problemtype_conic     = 1
problemtype_quadratic = 2

import MathOptInterface

boundflag_lower = 0x1
boundflag_upper = 0x2
boundflag_cone  = 0x4
boundflag_int   = 0x8
boundflag_all   = 0x0f

# Mapping of all constraint types to its index
struct ConstraintMap
    x_lessthan           :: Dict{UInt64,Int} # (SingleVariable,LessThan) -> constraint number
    x_greaterthan        :: Dict{UInt64,Int} # (SingleVariable,GreaterThan) -> constraint number
    x_equalto            :: Dict{UInt64,Int} # (SingleVariable,EqualTo) -> constraint number
    x_interval           :: Dict{UInt64,Int} # (SingleVariable,Interval) -> constraint number
    x_nonpositives       :: Dict{UInt64,Int} # (SingleVariable,Nonpositives) -> constraint number
    x_nonnegatives       :: Dict{UInt64,Int} # (SingleVariable,Nonnegatives) -> constraint number
    x_binary             :: Dict{UInt64,Int} # (SingleVariable,ZeroOne) -> constraint number
    x_integer            :: Dict{UInt64,Int} # (SingleVariable,Integer) -> constraint number

    xs_nonpositives      :: Dict{UInt64,Int} # (VectorOfVariables,Nonpositives) -> constraint number
    xs_nonnegatives      :: Dict{UInt64,Int} # (VectorOfVariables,Nonnegatives) -> constraint number
    xs_zeros             :: Dict{UInt64,Int} # (VectorOfVariables,Zeros) -> constraint number
    xs_reals             :: Dict{UInt64,Int} # (VectorOfVariables,Reals) -> constraint number
    xs_qcone             :: Dict{UInt64,Int} # (VectorOfVariables,SecondOrderCone) -> constraint number
    xs_rqcone            :: Dict{UInt64,Int} # (VectorOfVariables,RotatedSecondOrderCone) -> constraint number
    xs_psdconetriangle   :: Dict{UInt64,Int} # (VectorOfVariables,PositiveSemidefiniteConeTriangle) -> constraint number

    axb_lessthan         :: Dict{UInt64,Int} # (ScalarAffineFunction,LessThan) -> constraint number
    axb_greaterthan      :: Dict{UInt64,Int} # (ScalarAffineFunction,GreaterThan) -> constraint number
    axb_equalto          :: Dict{UInt64,Int} # (ScalarAffineFunction,EqualTo) -> constraint number
    axb_interval         :: Dict{UInt64,Int} # (ScalarAffineFunction,Interval) -> constraint number
    axb_nonpositives     :: Dict{UInt64,Int} # (ScalarAffineFunction,Nonpositives) -> constraint number
    axb_nonnegatives     :: Dict{UInt64,Int} # (ScalarAffineFunction,Nonnegatives) -> constraint number
    axb_binary           :: Dict{UInt64,Int} # (ScalarAffineFunction,ZeroOne) -> constraint number
    axb_integer          :: Dict{UInt64,Int} # (ScalarAffineFunction,Integer) -> constraint number

    axbs_nonpositives    :: Dict{UInt64,Int} # (VectorAffineFunction,Nonpositives) -> constraint number
    axbs_nonnegatives    :: Dict{UInt64,Int} # (VectorAffineFunction,Nonnegatives) -> constraint number
    axbs_zeros           :: Dict{UInt64,Int} # (VectorAffineFunction,Zeros) -> constraint number
    axbs_reals           :: Dict{UInt64,Int} # (VectorAffineFunction,Reals) -> constraint number
    axbs_qcone           :: Dict{UInt64,Int} # (VectorAffineFunction,SecondOrderCone) -> constraint number
    axbs_rqcone          :: Dict{UInt64,Int} # (VectorAffineFunction,RotatedSecondOrderCone) -> constraint number
    axbs_psdconetriangle :: Dict{UInt64,Int} # (VectorAffineFunction,PositiveSemidefiniteConeTriangle) -> constraint number
end

ConstraintMap() = ConstraintMap([Dict{UInt64,Int}() for i in 1:30]...)
select(cm::ConstraintMap,::Type{MathOptInterface.SingleVariable},               ::Type{MathOptInterface.LessThan{Float64}}) =                cm.x_lessthan
select(cm::ConstraintMap,::Type{MathOptInterface.SingleVariable},               ::Type{MathOptInterface.GreaterThan{Float64}}) =             cm.x_greaterthan
select(cm::ConstraintMap,::Type{MathOptInterface.SingleVariable},               ::Type{MathOptInterface.EqualTo{Float64}}) =                 cm.x_equalto
select(cm::ConstraintMap,::Type{MathOptInterface.SingleVariable},               ::Type{MathOptInterface.Interval{Float64}}) =                cm.x_interval
select(cm::ConstraintMap,::Type{MathOptInterface.SingleVariable},               ::Type{MathOptInterface.Nonpositives}) =                     cm.x_nonpositives
select(cm::ConstraintMap,::Type{MathOptInterface.SingleVariable},               ::Type{MathOptInterface.Nonnegatives}) =                     cm.x_nonnegatives
select(cm::ConstraintMap,::Type{MathOptInterface.SingleVariable},               ::Type{MathOptInterface.ZeroOne}) =                          cm.x_binary
select(cm::ConstraintMap,::Type{MathOptInterface.SingleVariable},               ::Type{MathOptInterface.Integer}) =                          cm.x_integer
select(cm::ConstraintMap,::Type{MathOptInterface.VectorOfVariables},            ::Type{MathOptInterface.Nonpositives}) =                     cm.xs_nonpositives
select(cm::ConstraintMap,::Type{MathOptInterface.VectorOfVariables},            ::Type{MathOptInterface.Nonnegatives}) =                     cm.xs_nonnegatives
select(cm::ConstraintMap,::Type{MathOptInterface.VectorOfVariables},            ::Type{MathOptInterface.Zeros}) =                            cm.xs_zeros
select(cm::ConstraintMap,::Type{MathOptInterface.VectorOfVariables},            ::Type{MathOptInterface.Reals}) =                            cm.xs_reals
select(cm::ConstraintMap,::Type{MathOptInterface.VectorOfVariables},            ::Type{MathOptInterface.SecondOrderCone}) =                  cm.xs_qcone
select(cm::ConstraintMap,::Type{MathOptInterface.VectorOfVariables},            ::Type{MathOptInterface.RotatedSecondOrderCone}) =           cm.xs_rqcone
select(cm::ConstraintMap,::Type{MathOptInterface.VectorOfVariables},            ::Type{MathOptInterface.PositiveSemidefiniteConeTriangle}) = cm.xs_psdconetriangle
select(cm::ConstraintMap,::Type{MathOptInterface.ScalarAffineFunction{Float64}},::Type{MathOptInterface.LessThan{Float64}}) =                cm.axb_lessthan
select(cm::ConstraintMap,::Type{MathOptInterface.ScalarAffineFunction{Float64}},::Type{MathOptInterface.GreaterThan{Float64}}) =             cm.axb_greaterthan
select(cm::ConstraintMap,::Type{MathOptInterface.ScalarAffineFunction{Float64}},::Type{MathOptInterface.EqualTo{Float64}}) =                 cm.axb_equalto
select(cm::ConstraintMap,::Type{MathOptInterface.ScalarAffineFunction{Float64}},::Type{MathOptInterface.Interval{Float64}}) =                cm.axb_interval
select(cm::ConstraintMap,::Type{MathOptInterface.ScalarAffineFunction{Float64}},::Type{MathOptInterface.ZeroOne}) =                          cm.axb_binary
select(cm::ConstraintMap,::Type{MathOptInterface.ScalarAffineFunction{Float64}},::Type{MathOptInterface.Integer}) =                          cm.axb_integer
select(cm::ConstraintMap,::Type{MathOptInterface.ScalarAffineFunction{Float64}},::Type{MathOptInterface.Nonpositives}) =                     cm.axb_nonpositives
select(cm::ConstraintMap,::Type{MathOptInterface.ScalarAffineFunction{Float64}},::Type{MathOptInterface.Nonnegatives}) =                     cm.axb_nonnegatives
select(cm::ConstraintMap,::Type{MathOptInterface.VectorAffineFunction{Float64}},::Type{MathOptInterface.Nonpositives}) =                     cm.axbs_nonpositives
select(cm::ConstraintMap,::Type{MathOptInterface.VectorAffineFunction{Float64}},::Type{MathOptInterface.Nonnegatives}) =                     cm.axbs_nonnegatives
select(cm::ConstraintMap,::Type{MathOptInterface.VectorAffineFunction{Float64}},::Type{MathOptInterface.Zeros}) =                            cm.axbs_zeros
select(cm::ConstraintMap,::Type{MathOptInterface.VectorAffineFunction{Float64}},::Type{MathOptInterface.Reals}) =                            cm.axbs_reals
select(cm::ConstraintMap,::Type{MathOptInterface.VectorAffineFunction{Float64}},::Type{MathOptInterface.SecondOrderCone}) =                  cm.axbs_qcone
select(cm::ConstraintMap,::Type{MathOptInterface.VectorAffineFunction{Float64}},::Type{MathOptInterface.RotatedSecondOrderCone}) =           cm.axbs_rqcone
select(cm::ConstraintMap,::Type{MathOptInterface.VectorAffineFunction{Float64}},::Type{MathOptInterface.PositiveSemidefiniteConeTriangle}) = cm.axbs_psdconetriangle

Base.getindex{F,D}(cm::ConstraintMap,r :: MathOptInterface.ConstraintReference{F,D}) = select(cm,F,D)[r.value]

struct MosekSolution
    whichsol :: Soltype
    solsta   :: Solsta
    prosta   :: Prosta

    xxstatus :: Vector{Stakey}
    xx       :: Vector{Float64}
    slx      :: Vector{Float64}
    sux      :: Vector{Float64}
    snx      :: Vector{Float64}

    cstatus  :: Vector{Stakey}
    xc       :: Vector{Float64}
    slc      :: Vector{Float64}
    suc      :: Vector{Float64}
    y        :: Vector{Float64}
end

"""
    MosekModel <: MathOptInterface.AbstractModel

Linear variables and constraint can be deleted. For some reason MOSEK
does not support deleting PSD variables.

Note also that adding variables and constraints will permanently add
some (currently between 1 and 3) Int64s that a `delete!` will not
remove. This ensures that References (Variable and constraint) that
are deleted are thereafter invalid.
"""
mutable struct MosekModel  <: MathOptInterface.AbstractSolverInstance
    task :: Mosek.MSKtask

    """
    Number of variables explicitly created by user
    """
    publicnumvar :: Int

    """
    """
    constrmap :: ConstraintMap

    """
    The total length of `x_block` matches the number of variables in
    the underlying task, and the number of blocks corresponds to the
    number variables allocated in the Model.
    """
    x_block      :: LinkedInts

    """
    One entry per scalar variable in the task indicating which bound
    types are defined
    """
    x_boundflags :: Vector{Int}

    """
    One entry per scalar variable in the task defining the number of
    variable constraints imposed on that variable. It is not allowed
    to delete a variable without deleting the constraints on it first
    (to get around the problem of deleting roots in conic constraints).
    """
    x_numxc :: Vector{Int}

    """
    One entry per variable-constraint
    """
    xc_block     :: LinkedInts

    """
    One entry per variable-constraint block indicating which bound
    types it defines. The values are binary ORed `boundflag_...` values.
    """
    xc_bounds    :: Vector{UInt8}
    """
    One entry per variable-constraint block indicating which cone it belongs to, 0 if none.
    """
    xc_coneid    :: Vector{Int}

    """
    One entry per scalar variable-constraint, defining which native
    variables the bound block covers.
    """
    xc_idxs      :: Vector{Int}

    ###########################
    """
    One scalar entry per constraint in the underlying task. One block
    per constraint allocated in the Model.
    """
    c_block :: LinkedInts

    """
    One entry per allocated scalar constraint. Defines the fixed term
    on the left-hand for each scalar constraint.
    """
    c_constant :: Vector{Float64}

    """
    One entry per allocated scalar constraint.
    Each element is either
    - 0: meaning: no slack, when domain is defined directly as a bound,
    - positive: a `x_block` reference, e.g. for qcones, or
    - negative: Negated index of a PSD variable in the underlying task
    """
    c_block_slack   :: Vector{Int}
    """
    One entry per variable-constraint block indicating which cone it belongs to, 0 if none.
    """
    c_coneid   :: Vector{Int}

    ###########################
    conecounter :: Int

    ###########################
    trm :: Rescode
    solutions :: Vector{MosekSolution}
end

function MathOptInterface.SolverInstance(solver::MosekSolver)
    t = maketask()
    try
        be_quiet = false
        for (option,val) in solver.options
            parname = string(option)
            if parname == "QUIET"
                be_quiet = be_quiet || convert(Bool,val)
            elseif startswith(parname, "MSK_IPAR_")
                Mosek.putnaintparam(t, parname, convert(Integer, val))
            elseif startswith(parname, "MSK_DPAR_")
                Mosek.putnadouparam(t, parname, convert(AbstractFloat, val))
            elseif startswith(parname, "MSK_SPAR_")
                Mosek.putnastrparam(t, parname, convert(AbstractString, val))
            elseif isa(val, Integer)
                parname = "MSK_IPAR_$parname"
                Mosek.putnaintparam(t, parname, val)
            elseif isa(val, AbstractFloat)
                parname = "MSK_DPAR_$parname"
                Mosek.putnadouparam(t, parname, val)
            elseif isa(val, AbstractString)
                parname = "MSK_SPAR_$parname"
                Mosek.putnastrparam(t, parname, val)
            else
                error("Value $val for parameter $option has unrecognized type")
            end
        end
        if ! be_quiet
            Mosek.putstreamfunc(t,Mosek.MSK_STREAM_LOG,printstream)
        end
        MosekModel(t,# task
                   0, # public numvar
                   ConstraintMap(), # public constraints
                   LinkedInts(),# c_block
                   Int[], # x_boundflags
                   Int[], # x_boundflags
                   LinkedInts(), # xc_block
                   UInt8[], # xc_bounds
                   Int[], # xc_coneid
                   Int[], # xc_idxs
                   LinkedInts(), # c_block
                   Float64[], # c_constant
                   Int[], # c_block_slack
                   Int[], # c_coneid
                   0, # cone counter
                   Mosek.MSK_RES_OK,
                   MosekSolution[]) # trm
    catch
        deletetask(m.task)
        rethrow()
    end
end

function MathOptInterface.free!(m::MosekModel)
    deletetask(m.task)
end

function MathOptInterface.optimize!(m::MosekModel)
    m.trm = optimize(m.task)
    m.solutions = MosekSolution[]
    if solutiondef(m.task,MSK_SOL_ITG)

        push!(m.solutions,
              MosekSolution(MSK_SOL_ITG,
                            getsolsta(m.task,MSK_SOL_ITG),
                            getprosta(m.task,MSK_SOL_ITG),
                            getskx(m.task,MSK_SOL_ITG),
                            getxx(m.task,MSK_SOL_ITG),
                            Float64[],
                            Float64[],
                            Float64[],
                            getskc(m.task,MSK_SOL_ITG),
                            getxc(m.task,MSK_SOL_ITG),
                            Float64[],
                            Float64[],
                            Float64[]))
    end
    if solutiondef(m.task,MSK_SOL_BAS)
        push!(m.solutions,
              MosekSolution(MSK_SOL_BAS,
                            getsolsta(m.task,MSK_SOL_BAS),
                            getprosta(m.task,MSK_SOL_BAS),
                            getskx(m.task,MSK_SOL_BAS),
                            getxx(m.task,MSK_SOL_BAS),
                            getslx(m.task,MSK_SOL_BAS),
                            getsux(m.task,MSK_SOL_BAS),
                            Float64[],
                            getskc(m.task,MSK_SOL_BAS),
                            getxc(m.task,MSK_SOL_BAS),
                            getslc(m.task,MSK_SOL_BAS),
                            getsuc(m.task,MSK_SOL_BAS),
                            gety(m.task,MSK_SOL_BAS)))
    end
    if solutiondef(m.task,MSK_SOL_ITR)
        push!(m.solutions,
              MosekSolution(MSK_SOL_ITR,
                            getsolsta(m.task,MSK_SOL_ITR),
                            getprosta(m.task,MSK_SOL_ITR),
                            getskx(m.task,MSK_SOL_ITR),
                            getxx(m.task,MSK_SOL_ITR),
                            getslx(m.task,MSK_SOL_ITR),
                            getsux(m.task,MSK_SOL_ITR),
                            getsnx(m.task,MSK_SOL_ITR),
                            getskc(m.task,MSK_SOL_ITR),
                            getxc(m.task,MSK_SOL_ITR),
                            getslc(m.task,MSK_SOL_ITR),
                            getsuc(m.task,MSK_SOL_ITR),
                            gety(m.task,MSK_SOL_ITR)))
    end
end

function MathOptInterface.writeproblem(m::MosekModel, filename :: String)
    putintparam(m.task,MSK_IPAR_OPF_WRITE_SOLUTIONS, MSK_ON)
    writedata(m.task,filename)
end

# For linear objectives we accept:
# EITER affine left-hand side and ranged, unbounded, half-open, fixed (equality), PSD or SOC domains
# OR affine and quadratic left-hand side, and ranged, unbounded, half-open, fixed (equality) domains (quadratic constraints must be unbounded or half-open)
#
# For non-quadratic problems we allow binary and integer variables (but not constraints)
function supportsconstraints(m::MosekSolver, constraint_types) :: Bool
    for (fun,dom) in constraint_types
        if  fun in [MathOptInterface.ScalarAffineFunction{Float64},
                    MathOptInterface.SingleVariable,
                    MathOptInterface.VectorAffineFunction{Float64}] &&
            dom in [MathOptInterface.Zeros,
                    MathOptInterface.Reals,
                    MathOptInterface.Nonnegatives,
                    MathOptInterface.Nonpositives,
                    MathOptInterface.GreaterThan{Float64},
                    MathOptInterface.LessThan{Float64},
                    MathOptInterface.EqualTo{Float64},
                    MathOptInterface.Interval{Float64},
                    MathOptInterface.SecondOrderCone,
                    MathOptInterface.RotatedSecondOrderCone,
                    MathOptInterface.PositiveSemidefiniteConeTriangle,
                    MathOptInterface.PositiveSemidefiniteConeScaled ]
            # ok
        elseif dom in [MathOptInterface.ZeroOne,
                       MathOptInterface.Integer] &&
                           fun in [MathOptInterface.SingleVariable,
                                   MathOptInterface.VectorOfVariables]
            # ok
        else
            return false
        end
    end
    true
end

MathOptInterface.supportsproblem(m::MosekSolver, ::Type{MathOptInterface.SingleVariable},                constraint_types) :: Bool = supportsconstraints(m,constraint_types)
MathOptInterface.supportsproblem(m::MosekSolver, ::Type{MathOptInterface.ScalarAffineFunction{Float64}}, constraint_types) :: Bool = supportsconstraints(m,constraint_types)
MathOptInterface.supportsproblem{F}(m::MosekSolver, ::Type{F}, constraint_types) :: Bool = false

ref2id(ref :: MathOptInterface.VariableReference) :: Int =
    if ref.value & 1 == 0
        Int(ref.value >> 1)
    else
        - Int(ref.value >> 1)
    end

ref2id(ref :: MathOptInterface.ConstraintReference) :: Int =
    if ref.value & 1 == 0
        Int(ref.value >> 1)
    else
        - Int(ref.value >> 1)
    end

id2vref(id :: Int) :: MathOptInterface.VariableReference =
    if id < 0
        MathOptInterface.VariableReference((UInt64(-id) << 1) | 1)
    else
        MathOptInterface.VariableReference(UInt64(id) << 1)
    end

include("objective.jl")
include("variable.jl")
include("constraint.jl")
include("attributes.jl")

export MosekModel

end # module


text:
module TestFirstLast

using Mimi
using Test

import Mimi: time_labels, set_first_last!

#
# Define some Components
#

@defcomp grosseconomy begin
	YGROSS	= Variable(index=[time])	# Gross output
	K	    = Variable(index=[time])	# Capital
	l	    = Parameter(index=[time])	# Labor
	tfp	    = Parameter(index=[time])	# Total factor productivity
	s	    = Parameter(index=[time])	# Savings rate
	depk	= Parameter()			    # Depreciation rate on capital - Note that it has no time index
	k0	    = Parameter()			    # Initial level of capital
	share	= Parameter()			    # Capital share

	function run_timestep(p, v, d, t)
		if is_first(t)
			v.K[t] 	= p.k0
		else
			v.K[t] 	= (1 - p.depk)^5 * v.K[t-1] + v.YGROSS[t-1] * p.s[t-1] * 5
		end
		v.YGROSS[t] = p.tfp[t] * v.K[t]^p.share * p.l[t]^(1-p.share)
	end
end

@defcomp emissions begin
	E 	    = Variable(index=[time])	# Total greenhouse gas emissions
	sigma	= Parameter(index=[time])	# Emissions output ratio
	YGROSS	= Parameter(index=[time])	# Gross output - Note that YGROSS is now a parameter

	function run_timestep(p, v, d, t)
		v.E[t] = p.YGROSS[t] * p.sigma[t]	# Note the p. in front of YGROSS
	end
end

#
# Test using first and last for one (the second) component
#

m = Model()
set_dimension!(m, :time, collect(2015:5:2110)) # 20 timesteps
add_comp!(m, grosseconomy)
add_comp!(m, emissions, first = 2020, last = 2105)

# check that the attributes of the ModelDef and ComponentDef(s) have been set_dimension
# as expected
@test collect(2015:5:2110) == time_labels(m.md) == [keys(m.md.namespace[:emissions].dim_dict[:time])...] == [keys(m.md.namespace[:grosseconomy].dim_dict[:time])...]
@test m.md.first == m.md.namespace[:grosseconomy].first
@test m.md.last == m.md.namespace[:grosseconomy].last
@test m.md.namespace[:emissions].first == 2020
@test m.md.namespace[:emissions].last == 2105

# update parameters for the grosseconomy component
update_param!(m, :grosseconomy, :l, [(1. + 0.015)^t *6404 for t in 1:20])
update_param!(m, :grosseconomy, :tfp, [(1 + 0.065)^t * 3.57 for t in 1:20])
update_param!(m, :grosseconomy, :s, ones(20).* 0.22)
update_param!(m, :grosseconomy, :depk, 0.1)
update_param!(m, :grosseconomy, :k0, 130.)
update_param!(m, :grosseconomy, :share, 0.3)

# update and connect parameters for the emissions component
@test_throws ErrorException update_param!(m, :emissions, :sigma, [(1. - 0.05)^t *0.58 for t in 1:19]) # the parameter needs to be length of model
update_param!(m, :emissions, :sigma, [(1. - 0.05)^t *0.58 for t in 1:20])
connect_param!(m, :emissions, :YGROSS, :grosseconomy, :YGROSS)

run(m)

# test that there are missing values in :emissions variables outside of the component's
# run period, and no missing values in the :grosseconomy variables
@test ismissing(m[:emissions, :E][1])
@test ismissing(m[:emissions, :E][20])
@test sum(ismissing.(m[:emissions, :E][2:19])) == 0
@test sum(ismissing.(m[:grosseconomy, :l])) == 0

# change the model dimension (widen it)
set_dimension!(m, :time, collect(2015:5:2115))

# check that the first, last, and time have been updated properly for both the
# ModelDef and ComponentDef(s)
@test collect(2015:5:2115) == time_labels(m.md) == [keys(m.md.namespace[:emissions].dim_dict[:time])...] == [keys(m.md.namespace[:grosseconomy].dim_dict[:time])...]
@test m.md.namespace[:grosseconomy].first == 2015 # same as original model dim
@test m.md.namespace[:grosseconomy].last == 2110 # same as original model dim
@test m.md.namespace[:emissions].first == 2020 # explicitly set
@test m.md.namespace[:emissions].last == 2105 # explicitly set

# reset any parameters that have a time dimension
update_param!(m, :grosseconomy, :l, [(1. + 0.015)^t *6404 for t in 1:21])
update_param!(m, :grosseconomy, :tfp, [(1 + 0.065)^t * 3.57 for t in 1:21])
update_param!(m, :grosseconomy, :s, ones(21).* 0.22)
update_param!(m, :emissions, :sigma, [(1. - 0.05)^t *0.58 for t in 1:21])

run(m)

# test that there are missing values in :emissions variables outside of the component's
# run period, and no missing values in the :grosseconomy variables
@test ismissing(m[:emissions, :E][1])
@test all(ismissing, m[:emissions, :E][20:21])
@test sum(ismissing.(m[:emissions, :E][2:19])) == 0
@test sum(ismissing.(m[:grosseconomy, :l])) == 0

#
# Test bounds - both start late - with is_first()
#

m = Model()
set_dimension!(m, :time, collect(2015:5:2110)) # 20 timesteps
add_comp!(m, grosseconomy, first = 2020)
add_comp!(m, emissions, first = 2020)

# update parameters for the grosseconomy component
update_param!(m, :grosseconomy, :l, [(1. + 0.015)^t *6404 for t in 1:20])
update_param!(m, :grosseconomy, :tfp, [(1 + 0.065)^t * 3.57 for t in 1:20])
update_param!(m, :grosseconomy, :s, ones(20).* 0.22)
update_param!(m, :grosseconomy, :depk, 0.1)
update_param!(m, :grosseconomy, :k0, 130.)
update_param!(m, :grosseconomy, :share, 0.3)

# update and connect the parameters for the emissions component
update_param!(m, :emissions, :sigma, [(1. - 0.05)^t *0.58 for t in 1:20])
connect_param!(m, :emissions, :YGROSS, :grosseconomy, :YGROSS)

run(m)

# neither component should have a value for the first timestep
@test ismissing(m[:emissions, :E][1])
@test ismissing(m[:grosseconomy, :YGROSS][1])

# Test bounds - both end early
#

m = Model()
set_dimension!(m, :time, collect(2015:5:2110)) # 20 timesteps
add_comp!(m, grosseconomy, last = 2105)
add_comp!(m, emissions, last = 2105)

# update parameters for the grosseconomy component
update_param!(m, :grosseconomy, :l, [(1. + 0.015)^t *6404 for t in 1:20])
update_param!(m, :grosseconomy, :tfp, [(1 + 0.065)^t * 3.57 for t in 1:20])
update_param!(m, :grosseconomy, :s, ones(20).* 0.22)
update_param!(m, :grosseconomy, :depk, 0.1)
update_param!(m, :grosseconomy, :k0, 130.)
update_param!(m, :grosseconomy, :share, 0.3)

# update and connect parameters for the emissions component
update_param!(m, :emissions, :sigma, [(1. - 0.05)^t *0.58 for t in 1:20])
connect_param!(m, :emissions, :YGROSS, :grosseconomy, :YGROSS)

run(m)

@test ismissing(m[:emissions, :E][20])
@test ismissing(m[:grosseconomy, :YGROSS][20])

#
# Test bounds - components starting or ending before/after the model
#

m = Model()
set_dimension!(m, :time, collect(2015:5:2110)) # 20 timesteps

# components cannot start before or end after the model's time dimension
@test_throws ErrorException add_comp!(m, grosseconomy, first = 2000)
add_comp!(m, grosseconomy)
@test_throws ErrorException add_comp!(m, emissions, last = 2120)

#
# Test is_first and is_last
#

# is_first and is_last should be equivalent to comparing the timestep to the first
# and last years provided, and t.t should be 1 when the year is first i.e. the
# first year the component actually runs for.

@defcomp MyComp begin
	a 	= Variable(index=[time])
	function run_timestep(p, v, d, t)
		if t == TimestepValue(2) v.a[t] = -999
		elseif t == TimestepValue(9) v.a[t] = 999
		else v.a[t] = t.t
		end
	end
end

@defcomp MyComp2 begin
	a 	= Variable(index=[time])
	function run_timestep(p, v, d, t)
		if is_first(t) v.a[t] = -999
		elseif is_last(t) v.a[t] = 999
		else v.a[t] = t.t
		end
	end
end

m = Model()
set_dimension!(m, :time, collect(1:10))
add_comp!(m, MyComp, first = 2, last = 9)
run(m)

m2 = Model()
set_dimension!(m2, :time, collect(1:10))
add_comp!(m2, MyComp, first = 2, last = 9)
run(m2)

@test ismissing(m[:MyComp, :a][1])
@test ismissing(m[:MyComp, :a][end])
@test (m[:MyComp, :a])[2:9] == (m2[:MyComp, :a])[2:9] == [-999., 2., 3., 4., 5., 6., 7., 999.]

#
# TimestepIndex and TimestepValue
#

# TimestepIndex: Test Equality - should match up with t.t
@defcomp MyComp begin
	a 	= Variable(index=[time])
	b 	= Variable(index=[time])
	function run_timestep(p, v, d, t)
		v.a[t] = t == TimestepIndex(t.t)
		v.b[t] = 1.
	end
end

m = Model()
set_dimension!(m, :time, collect(1:15))
add_comp!(m, MyComp, first = 5, last = 10)
run(m)

for i in collect(1:15)
	@test m[:MyComp, :a][i] === m[:MyComp, :b][i]
end

# TimestepValue: Test Equality - should match up with the time index
@defcomp MyComp begin
	a 	= Variable(index=[time])
	b 	= Variable(index=[time])
	function run_timestep(p, v, d, t)
		v.a[t] = t == TimestepValue(t.t + 4)
		v.b[t] = 1.
	end
end

m = Model()
set_dimension!(m, :time, collect(1:15))
add_comp!(m, MyComp, first = 5, last = 10)
run(m)

for i in collect(1:15)
	@test m[:MyComp, :a][i] === m[:MyComp, :b][i]
end

# TimestepIndex: Test that Get and Set Index are Relative to Component, not Model
@defcomp MyComp begin
	a 	= Variable(index=[time])
	b 	= Variable(index=[time])
	function run_timestep(p, v, d, t)
		v.a[t] = t.t
		v.b[TimestepIndex(t.t)] = t.t
	end
end

m = Model()
set_dimension!(m, :time, collect(1:15))
add_comp!(m, MyComp, first = 5, last = 10)
run(m)
for (i, element) in enumerate(m[:MyComp, :a])
	@test element === m[:MyComp, :b][i]
end

@defcomp MyComp begin
	a 	= Variable(index=[time])
	function run_timestep(p, v, d, t)
		if t == TimestepIndex(1)
			v.a[TimestepIndex(1)] = 1
		else
			v.a[t] = 0
		end
	end
end

for year in collect(1995:1999)
	m = Model()
	set_dimension!(m, :time, collect(1995:2000))
	add_comp!(m, MyComp, first = year)
	run(m)
	idx = year - 1995 + 1
	@test m[:MyComp, :a][idx] == 1.0
end

# TimestepValue: Test that Get and Set Index are Relative to Component, not Model
@defcomp MyComp begin
	a 	= Variable(index=[time])
	function run_timestep(p, v, d, t)
		if t == TimestepValue(1999)
			v.a[TimestepValue(1999)] = 1
		else
			v.a[t] = 0
		end
	end
end

for year in collect(1995:1999)
	m = Model()
	set_dimension!(m, :time, collect(1995:2000))
	add_comp!(m, MyComp, first = year)
	run(m)
	@test m[:MyComp, :a][5] == 1.0
end

# test composite

@defcomp Comp1 begin
    par_1_1 = Parameter(index=[time])      # external input
    var_1_1 = Variable(index=[time])       # computed
    foo = Parameter()

    function run_timestep(p, v, d, t)
        v.var_1_1[t] = p.par_1_1[t]
    end
end

@defcomp Comp2 begin
    par_2_1 = Parameter(index=[time])      # connected to Comp1.var_1_1
    par_2_2 = Parameter(index=[time])      # external input
    var_2_1 = Variable(index=[time])       # computed
    foo = Parameter()

    function run_timestep(p, v, d, t)
        v.var_2_1[t] = p.par_2_1[t] + p.foo * p.par_2_2[t]
    end
end

@defcomp Comp3 begin
    par_3_1 = Parameter(index=[time])      # connected to Comp2.var_2_1
    var_3_1 = Variable(index=[time])       # external output
    foo = Parameter(default=30)

    function run_timestep(p, v, d, t)
        # @info "Comp3 run_timestep"
        v.var_3_1[t] = p.par_3_1[t] * 2
    end
end

@defcomp Comp4 begin
    par_4_1 = Parameter(index=[time])      # connected to Comp2.var_2_1
    var_4_1 = Variable(index=[time])        # external output
    foo = Parameter(default=300)

    function run_timestep(p, v, d, t)
        # @info "Comp4 run_timestep"
        v.var_4_1[t] = p.par_4_1[t] * 2
    end
end

@defcomposite A begin
    Component(Comp1)
    Component(Comp2)

    foo1 = Parameter(Comp1.foo)
    foo2 = Parameter(Comp2.foo)

    var_2_1 = Variable(Comp2.var_2_1)

    connect(Comp2.par_2_1, Comp1.var_1_1)
    connect(Comp2.par_2_2, Comp1.var_1_1)
end

@defcomposite B begin
    Component(Comp3)
    Component(Comp4)

    foo3 = Parameter(Comp3.foo)
    foo4 = Parameter(Comp4.foo)

    var_3_1 = Variable(Comp3.var_3_1)
end

@defcomposite top begin
    Component(A)

    fooA1 = Parameter(A.foo1)
    fooA2 = Parameter(A.foo2)

    Component(B, first = 2010, last = 2015)
    foo3 = Parameter(B.foo3)
    foo4 = Parameter(B.foo4)

    var_3_1 = Variable(B.var_3_1)

    connect(B.par_3_1, A.var_2_1)
    connect(B.par_4_1, B.var_3_1)
end

# We have created the following composite structure:
#
#          top
#        /    \
#       A       B
#     /  \     /  \
#    1    2   3    4

m = Model()
set_dimension!(m, :time, 2005:2020)
add_comp!(m, top, nameof(top))

update_param!(m, :top, :fooA1, 1)
update_param!(m, :top, :fooA2, 2)
update_param!(m, :top, :foo3, 10)
update_param!(m, :top, :foo4, 20)
update_param!(m, :top, :par_1_1, collect(1:length(time_labels(m))))

run(m)

# check that first and last moved through properly
@test m.md.namespace[:top].first == m.md.namespace[:top].namespace[:A].first == m.md.namespace[:top].namespace[:A].namespace[:Comp1].first== 2005
@test m.md.namespace[:top].last == m.md.namespace[:top].namespace[:A].last ==m.md.namespace[:top].namespace[:A].namespace[:Comp1].last == 2020

@test m.md.namespace[:top].namespace[:B].first == m.md.namespace[:top].namespace[:B].namespace[:Comp3].first== 2010
@test m.md.namespace[:top].namespace[:B].last ==m.md.namespace[:top].namespace[:B].namespace[:Comp3].last == 2015

 #
 # Test set_first_last! function
 #

 m = Model()
 set_dimension!(m, :time, collect(2015:5:2110)) # 20 timesteps
 add_comp!(m, grosseconomy)
 add_comp!(m, emissions)

 # check that the attributes of the ModelDef and ComponentDef(s) have been set
 # as expected
 @test collect(2015:5:2110) == time_labels(m.md) == [keys(m.md.namespace[:emissions].dim_dict[:time])...] == [keys(m.md.namespace[:grosseconomy].dim_dict[:time])...]
 @test m.md.first == m.md.namespace[:grosseconomy].first == m.md.namespace[:emissions].first
 @test m.md.last == m.md.namespace[:grosseconomy].last == m.md.namespace[:emissions].last

 # now set the emissions first and last and check
 set_first_last!(m, :emissions, first = 2020, last = 2105)

 @test collect(2015:5:2110) == time_labels(m.md) == [keys(m.md.namespace[:emissions].dim_dict[:time])...] == [keys(m.md.namespace[:grosseconomy].dim_dict[:time])...]
 @test m.md.first == m.md.namespace[:grosseconomy].first
 @test m.md.last == m.md.namespace[:grosseconomy].last
 @test m.md.namespace[:emissions].first == 2020
 @test m.md.namespace[:emissions].last == 2105

 # check warnings
 @test_throws ErrorException set_first_last!(m, :grosseconomy, first = 2000) # too early
 @test_throws ErrorException set_first_last!(m, :grosseconomy, last = 3000)

end #module


text:
__precompile__(true)

module SphericalHarmonics

export *

"""
	sizeP(maxDegree)

Return the size of the set of Associated Legendre Polynomials ``P_l^m(x)`` of
degree less than or equal to the given maximum degree
"""
sizeP(maxDegree) = div((maxDegree + 1) * (maxDegree + 2), 2)

"""
	sizeY(maxDegree)

Return the size of the set of real spherical harmonics ``Y_{l,m}(θ,φ)`` of
degree less than or equal to the given maximum degree
"""
sizeY(maxDegree) = (maxDegree +1 ) * (maxDegree + 1)

"""
	index_p(l,m)

Return the index into a flat array of Associated Legendre Polynomials ``P_l^m``
for the given indices ``(l,m)``.
``P_l^m`` are stored in l-major order i.e. [P(0,0), [P(1,0), P(1,1), P(2,0), ...]
"""
index_p(l,m) = m + div(l*(l+1), 2) + 1

"""
	index_y(l,m)

Return the index into a flat array of real spherical harmonics ``Y_{l,m}``
for the given indices ``(l,m)``.
``Y_{l,m}`` are stored in l-major order i.e.
[Y(0,0), [Y(1,-1), Y(1,0), Y(1,1), Y(2,-2), ...]
"""
index_y(l,m) = m + l + (l*l) + 1

type ALPCoefficients
	A::Array{Float64}
	B::Array{Float64}
	function ALPCoefficients(maxDegree::Int)
		A = Array{Float64}(sizeP(maxDegree))
		B = Array{Float64}(sizeP(maxDegree))
		new(A, B)
	end
end

"""
	compute_coefficients(L)

Precompute coefficients ``a_l^m`` and ``b_l^m`` for all l <= L, m <= l
"""
function compute_coefficients(L::Int)
	coeff = ALPCoefficients(L)
	for l in 2:L
		ls = l*l
		lm1s = (l-1) * (l-1)
		for m in 0:(l-2)
			ms = m * m
			coeff.A[index_p(l, m)] = sqrt((4 * ls - 1.) / (ls - ms))
			coeff.B[index_p(l, m)] = -sqrt((lm1s - ms) / (4 * lm1s - 1.))
		end
	end
	return coeff
end

"""
	compute_coefficients(L)

Create an array large enough to store an entire set of Associated Legendre
Polynomials ``P_l^m(x)`` of maximum degree L.
"""
function allocate_p(L::Int)
	P = Array{Float64}(sizeP(L))
end

"""
	compute_p(L, x, coeff, P)

Compute an entire set of Associated Legendre Polynomials ``P_l^m(x)``
using the given coefficients, and store in the array P.
"""
function compute_p(L::Int, x::Float64, coeff::ALPCoefficients,
	P::Array{Float64,1})
	@assert length(coeff.A) >= sizeP(L)
	@assert length(coeff.B) >= sizeP(L)
	@assert length(P) >= sizeP(L)

	sintheta = sqrt(1. - x * x)
	temp = 0.39894228040143267794 # = sqrt(0.5/M_PI)
	P[index_p(0, 0)] = temp

	if (L > 0)
		const SQRT3 = 1.7320508075688772935
		P[index_p(1, 0)] = x * SQRT3 * temp
		const SQRT3DIV2 = -1.2247448713915890491
		temp = SQRT3DIV2 * sintheta * temp
		P[index_p(1, 1)] = temp

		for l in 2:L
			for m in 0:(l-2)
				P[index_p(l, m)] = coeff.A[index_p(l, m)] *(X * P[index_p(l - 1, m)]
						     + coeff.B[index_p(l, m)] * P[index_p(l - 2, m)])
			end
			P[index_p(l, l - 1)] = x * sqrt(2 * (l - 1) + 3) * temp
			temp = -sqrt(1.0 + 0.5 / l) * sintheta * temp
			P[index_p(l, l)] = temp
		end
	end
	return P
end

"""
	compute_p(L, x)

Compute an entire set of Associated Legendre Polynomials ``P_l^m(x)`` where
``0 ≤ l ≤ L`` and ``0 ≤ m ≤ l``. Assumes ``|x| ≤ 1``.
"""
function compute_p(L::Int, x::Float64)
	P = Array{Float64}(sizeP(L))
	coeff = compute_coefficients(L)
	compute_p(L, x, coeff, P)
	return P
end

"""
	compute_y(L, P, φ, Y)

Compute an entire set of real spherical harmonics ``Y_{l,m}(θ,φ)``
using the given Associated Legendre Polynomials ``P_l^m(cos θ)``
and store in array Y
"""
function compute_y(L::Int, P::Array{Float64,1}, phi::Float64,
	Y::Array{Float64,1})
	@assert length(P) >= sizeP(L)
	@assert length(Y) >= sizeY(L)

	const SQRT2 = 1.41421356237309504880
	for l in 0:L
		Y[index_y(l, 0)] = P[index_p(l, 0)] * 0.5 * SQRT2
	end

	# NR2 5.5.4-5.5.5
	c1 = 1.0; c2 = cos(phi)
	s1 = 0.0; s2 = -sin(phi)
	tc = 2.0 * c2
	for m in 1:L
		s = tc * s1 - s2
		c = tc * c1 - c2
		s2 = s1
		s1 = s
		c2 = c1
		c1 = c
		for l in m:L
			Y[index_y(l, -m)] = P[index_p(l, m)] * s
			Y[index_y(l, m)] = P[index_p(l, m)] * c
		end
	end
end

"""
	compute_y(L, x, φ)

Compute an entire set of real spherical harmonics ``Y_{l,m}(θ,φ)`` for
``x = cos θ`` where ``0 ≤ l ≤ L`` and ``-l ≤ m ≤ l``.
"""
function compute_y(L::Int, x::Float64, phi::Float64)
	P = Array{Float64}(sizeP(L))
	coeff = compute_coefficients(L)
	compute_p(L, x, coeff, P)
	Y = Array{Float64}(sizeY(L))
	compute_y(L, x, phi, P, Y)
	return Y
end

end


text:
# constants.jl

#################### Helper function(s) ####################

function bytemask(num_quarters::Integer, qoffset::Integer)
    mask = 0b11_11_11_11

    for j in qoffset:(qoffset + num_quarters - 1)
        mask -= 0b11 << 2*j
    end

    mask
end

function flipsnp(x::UInt8)
    x = (~x) & 0b00_00_00_11
    (x << 1 | x >> 1) & 0b00_00_00_11
end

#################### Constants ####################

### .bed file constants
# Note: this is big endian [0x6c, 0x1b]
const plinkmagic = 0b0110_1100_0001_1011
const modes = Dict(0b0000_0001 => :SNPmajor,
                   0b0000_0000 => :SNPminor)

### BED format to Raw format mapping
const NA_byte = 0b11
const quarterstohuman = (0b10, NA_byte, 0b01, 0b00)

const flipBEDbytes = tuple([(flipsnp(snp4) << 6) | (flipsnp(snp3) << 4) | (flipsnp(snp2) << 2) | flipsnp(snp1) for
                            snp4 in 0b00:0b11 for
                            snp3 in 0b00:0b11 for
                            snp2 in 0b00:0b11 for
                            snp1 in 0b00:0b11]...)

"""
    const bytetoquarters::Tuple{Vararg{Tuple{UInt8,UInt8,UInt8,UInt8}, 256}}

A constant array storing all 256 bytes split into 4 quarters. Note
that the mapping is given by `bytetoquarters[BEDbyte + 1]` because of
julia's 1-indexing.

# Example
```julia
julia> BEDbyte = 0b11101001
0xe9
julia> BEDMatrices.bytetoquarters[BEDbyte + 1]
4-element Array{UInt8,1}:
 0x03
 0x01
 0x01
 0x00
```

This can be understood as breaking the byte into quarters, in reverse
order: `0b01`, `0b10`, `0b10`, `0b11`, and then changing from BED
format to `rawformat` with `0x03` representing missing value.

"""
const bytetoquarters = tuple([(quarterstohuman[(byte & 0b00_00_00_11) + 1],
                               quarterstohuman[(byte & 0b00_00_11_00) >> 2 + 1],
                               quarterstohuman[(byte & 0b00_11_00_00) >> 4 + 1],
                               quarterstohuman[(byte & 0b11_00_00_00) >> 6 + 1]) for
                              byte in 0x00:0xff]...)

"""
    const defaultbytemap::Matrix{Tuple{Vararg{UInt8, 4}}}

This is a matrix of `size(defaultbytemap) == (256, 2)`. The first
columne is the same as `bytetoquarters`: `defaultbytemap[:, false + 1]
== collect(bytetoquarters)`, see docstring for `bytetoquarters`. The
second column corresponds to doing a major--minor allele flip.

"""
const defaultbytemap = hcat([(quarterstohuman[(byte & 0b00_00_00_11) + 1],
                              quarterstohuman[(byte & 0b00_00_11_00) >> 2 + 1],
                              quarterstohuman[(byte & 0b00_11_00_00) >> 4 + 1],
                              quarterstohuman[(byte & 0b11_00_00_00) >> 6 + 1]) for
                             byte in 0x00:0xff],
                            [(quarterstohuman[(byte & 0b00_00_00_11) + 1],
                              quarterstohuman[(byte & 0b00_00_11_00) >> 2 + 1],
                              quarterstohuman[(byte & 0b00_11_00_00) >> 4 + 1],
                              quarterstohuman[(byte & 0b11_00_00_00) >> 6 + 1]) for
                             byte in flipBEDbytes])

### BED byte to RAW math
const onebyte = 0b10_10_10_10  # vector of ones in RAW format

const partialbytemasks = [bytemask(numq, qoff) for numq in 1:4, qoff in 0:3]  # tuple?
# const partialbytemasks = [BEDMatrices.bytemask(numq, qoff) for numq in 1:4, qoff in 0:3]

const inversebytemap = Dict(bytetoquarters[b + 1] => b for b in 0x0:0xff)
const natozeromap = [inversebytemap[map(q -> q == NA_byte ? 0b00 : q, bytetoquarters[b + 1])] for b in 0x0:0xff]  # tuple?
const na_sup_map = [inversebytemap[map(q -> q == NA_byte ? 0b01 : 0b00, bytetoquarters[b + 1])] for b in 0x0:0xff]  # tuple?

const nacountmap = [count(q -> q === NA_byte, bytetoquarters[b + 1]) for b in 0x0:0xff]  # tuple?
const distributionmap = tuple(map(b -> (count(q -> q === 0b00, bytetoquarters[b + 1]),
                                        count(q -> q === 0b01, bytetoquarters[b + 1]),
                                        count(q -> q === 0b10, bytetoquarters[b + 1]),
                                        count(q -> q === 0b11, bytetoquarters[b + 1])), 0x0:0xff)...)

const hasNAmap = tuple([c > 0 for c in nacountmap]...)

const bytebytemulttable = [dot(collect(bytetoquarters[natozeromap[b1 + 1] + 1]),
                               collect(bytetoquarters[natozeromap[b2 + 1] + 1])) for b1 in 0x0:0xff, b2 in 0x0:0xff]


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license
module AbstractLogic

##############################################################################
##
## Dependencies
##
##############################################################################
using Markdown, StatsBase, Crayons, ReplMaker

##############################################################################
##
## Exported methods and types (in addition to everything reexported above)
##
##############################################################################

export abstractlogic,
       activehistory,
       checkfeasible,
       discover,
       dependenton,
       LogicalCombo,
       logicalparse,
       search,
       showfeasible,
       dashboard!,
       dashboard,
       expand,
       help,
       nfeasible,
       percievedfeasible,
       replerror,
       activehistory,
       replset,
       sessionhistory,
       setcompare

##############################################################################
##
## Load main files
##
##############################################################################

include("main/LogicalCombo.jl")

include("main/feasible.jl")

include("main/counter.jl")

include("main/logicalparse.jl")

include("main/expand.jl")

include("main/usefulfunctions.jl")

include("main/definelogicalset.jl")

include("main/dependenton.jl")

include("main/operatoreval.jl")

include("main/superoperator.jl")

include("main/metaoperator.jl")

include("main/operatorspawn.jl")

include("main/search.jl")

include("main/checkfeasible.jl")

include("main/showfeasible.jl")

include("main/setcompare.jl")

include("main/help.jl")

include("main/discover.jl")

##############################################################################
##
## Load REPL
##
##############################################################################

include("repl/History.jl")

include("repl/ALback_next.jl")

include("repl/ALcheck.jl")

include("repl/ALclear.jl")

include("repl/ALcompare.jl")

include("repl/ALhistory.jl")

include("repl/ALimport_export.jl")

include("repl/ALparse.jl")

include("repl/ALpreserve_restore.jl")

include("repl/ALrange.jl")

include("repl/ALsearch.jl")

include("repl/ALshow.jl")

include("repl/replcalls.jl")

include("repl/testcall.jl")

include("repl/repl.jl")

isdefined(Base, :active_repl) && include("repl/replmaker.jl")

end # module AbstractLogic


text:
base_dir = dirname(dirname(pathof(PowerSystems)))

sys = PowerSystems.parse_standard_files(joinpath(MATPOWER_DIR, "case5_re.m"))

@testset "Check bus index" begin
    @test sort([b.number for b in collect(get_components(Bus, sys))]) == [1, 2, 3, 4, 10]
    @test sort(collect(Set([b.arc.from.number for
                b in collect(get_components(Branch,sys))]))) == [1, 2, 3, 4]
    @test sort(collect(Set([b.arc.to.number for
                b in collect(get_components(Branch,sys))]))) == [2, 3, 4, 10]

    # TODO: add test for loadzones testing MAPPING_BUSNUMBER2INDEX

end