text:
Base.@kwdef mutable struct AMGXVector <: AMGXObject
    handle::API.AMGX_vector_handle = API.AMGX_vector_handle(C_NULL)
    mode::Union{Mode, Nothing} = nothing
    resources::Union{Resources, Nothing} = nothing
    function AMGXVector(handle::API.AMGX_vector_handle, mode::Union{Mode, Nothing},
                       resources::Union{Resources, Nothing})
        v = new(handle, mode, resources)
        finalizer(warn_not_destroyed_on_finalize, v)
        return v
    end
end
function Base.show(io::IO, mime::MIME"text/plain", v::AMGXVector)
    invoke(show, Tuple{IO, MIME"text/plain", AMGXObject}, io, mime, v)
    v.handle == C_NULL && return
    n, block_dim = vector_get_size(v)
    if n !== 0
        if block_dim == 1
            print(io, " of length $n")
        else
            print(io, " of length $n⋅$block_dim")
        end
    end
end
get_api_destroy_call(::Type{AMGXVector}) = API.AMGX_vector_destroy
function dec_refcount_parents(v::AMGXVector)
    dec_refcount!(v.resources)
    v.resources = nothing
    nothing
end

function create!(v::AMGXVector, res::Resources, mode::Mode)
    vec_handle_ptr = Ref{API.AMGX_vector_handle}()
    amgx_mode = API.AMGX_Mode(Int(mode))
    @checked API.AMGX_vector_create(vec_handle_ptr, res.handle, amgx_mode)
    v.handle = vec_handle_ptr[]
    v.mode = mode
    v.resources = res
    inc_refcount!(v.resources)
    return v
end
AMGXVector(res::Resources, mode::Mode) = create!(AMGXVector(), res, mode)

function upload!(v::AMGXVector, data::VectorOrCuVector; block_dim::Int=1)
    n, rem = divrem(length(data), block_dim)
    if rem != 0
        throw(ArgumentError("vector length ($(length(data))) not an integer multiple of block dimension ($(block_dim))"))
    end
    GC.@preserve data begin
        upload_raw!(v, pointer(data), n, block_dim)
    end
    return v
end

function upload_raw!(v::AMGXVector, data::PtrOrCuPtrUnion{T}, n::Int, block_dim::Int=1) where T <: Union{Float64, Float32}
    # TODO: This should be handled in the ccall
    data = reinterpret(Ptr{T}, data)
    vT = vector_type(v.mode)
    if vT != T
        throw(ArgumentError("inconsistent AMGX vector mode ($vT) with element type of upload ($T)"))
    end
    @checked API.AMGX_vector_upload(v.handle, n, block_dim, data)
    return v
end

function download(v::AMGXVector)
    vT = vector_type(v.mode)
    buffer = Vector{vT}(undef, length(v))
    download!(buffer, v)
    return buffer
end
Base.Vector(v::AMGXVector) = download(v)
Base.Array(v::AMGXVector) = download(v)
CUDA.CuArray(v::AMGXVector) = CUDA.CuVector(v)
function CUDA.CuVector(v::AMGXVector)
    vT = vector_type(v.mode)
    buffer = CUDA.CuVector{vT}(undef, length(v))
    download!(buffer, v)
    return buffer
end
function download!(buffer::VectorOrCuVector, v::AMGXVector)
    m, n = length(buffer), length(v)
    if length(buffer) !== n
        throw(ArgumentError("invalid buffer length, got $m, expected $n"))
    end
    GC.@preserve buffer begin
        download_raw(pointer(buffer), v)
    end
    return buffer
end
Base.copy!(buffer::VectorOrCuVector, v::AMGXVector) = download!(buffer, v)

function download_raw(ptr::PtrOrCuPtrUnion{T}, v::AMGXVector) where T <: Union{Float64, Float32}
    # TODO: This should be handled in the ccall
    ptr = reinterpret(Ptr{T}, ptr)
    vT = vector_type(v.mode)
    if vT != T
        throw(ArgumentError("inconsistent AMGX vector mode ($vT) with element type of download ($T)"))
    end
    @checked API.AMGX_vector_download(v.handle, ptr)
end

function vector_get_size(v::AMGXVector)
    n_ptr, block_dim_ptr = Ref{Cint}(), Ref{Cint}()
    @checked API.AMGX_vector_get_size(v.handle, n_ptr, block_dim_ptr)
    return Int(n_ptr[]), Int(block_dim_ptr[])
end
Base.length(v::AMGXVector) = prod(vector_get_size(v))

function set_zero!(v::AMGXVector, n::Int=vector_get_size(v)[1]; block_dim=vector_get_size(v)[2])
    @checked API.AMGX_vector_set_zero(v.handle, Cint(n), Cint(block_dim))
    return v
end


text:
function init_configuration()
    config_initialized[] && return true
    ret = LibLikwid.init_configuration()
    if ret == 0
        config_initialized[] = true
        _config[] = unsafe_load(LibLikwid.get_configuration())
        _build_jl_config()
        return true
    end
    return false
end

function destroy_configuration()
    config_initialized[] || return false
    ret = LibLikwid.destroy_configuration()
    if ret == 0
        config_initialized[] = false
        return true
    end
    return false
end

function _build_jl_config()
    c = _config[]

    config[] = Likwid_Configuration(
        c.configFileName == C_NULL ? "" : unsafe_string(c.configFileName),
        c.topologyCfgFileName == C_NULL ? "" : unsafe_string(c.topologyCfgFileName),
        c.daemonPath == C_NULL ? "" : unsafe_string(c.daemonPath),
        c.groupPath == C_NULL ? "" : unsafe_string(c.groupPath),
        c.daemonMode,
        c.maxNumThreads,
        c.maxNumNodes,
    )
    return nothing
end

function get_configuration()
    if !config_initialized[]
        init_configuration() || error("Couldn't init configuration.")
    end
    return config[]
end

ext: jl


text:
struct MiniApp
	jg::JobGraph
	fromGUI::Channel{Pair{Symbol,Vector}}
	toGUI::Channel{Pair{Symbol,Any}}
	lastSchedulerTime::Ref{UInt64}
end
MiniApp() = MiniApp(JobGraph(verbose=true), Channel{Pair{Symbol,Vector}}(Inf), Channel{Pair{Symbol,Any}}(Inf), Ref{UInt64}(0))

setvalue(app::MiniApp, value::Array{Any}) = put!(app.fromGUI, :setvalue=>value)

function init(app::MiniApp)
	setvalue(app, Any["samplesimplexmethod", "SA"])
	setvalue(app, Any["loadrowsassamples", "true"])
	setvalue(app, Any["logtransform", "false"])
	setvalue(app, Any["logtransformoffset", "1"])
	setvalue(app, Any["normalizemethod", "Mean=0"])
	setvalue(app, Any["dimreductionmethod", "PMA"])
	setvalue(app, Any["knearestneighbors", "0"])
	setvalue(app, Any["distnearestneighbors", "0"])
	setvalue(app, Any["xaxis", "1"])
	setvalue(app, Any["yaxis", "2"])
	setvalue(app, Any["zaxis", "3"])
	setvalue(app, Any["plotwidth", "1024"])
	setvalue(app, Any["plotheight", "768"])
	setvalue(app, Any["showpoints", "true"])
	setvalue(app, Any["showlines", "true"])
	setvalue(app, Any["showtriangles", "false"])
	setvalue(app, Any["markersize", "4"])
	setvalue(app, Any["linewidth", "1"])
	setvalue(app, Any["triangleopacity", "0.05"])
	setvalue(app, Any["colorby", "Auto"])
	setvalue(app, Any["exportmode", "Variables"])
	setvalue(app, Any["exportsingledim", "1"])
	setvalue(app, Any["exportsinglesort", "Abs"])
	setvalue(app, Any["exportmultipledim", "3"]	)
end

function runall(app::MiniApp)
	while isready(app.fromGUI) || wantstorun(app.jg.scheduler) || isactive(app.jg.scheduler)
		process_step(app.jg, app.fromGUI, app.toGUI, app.lastSchedulerTime; verbosityLevel=2) || return false
	end
	true
end

function factorizationcmp(F1,F2)
	@test size(F1.U)==size(F2.U)
	@test size(F1.S)==size(F2.S)
	@test size(F1.V)==size(F2.V)
	@test size(F1.Vt)==size(F2.Vt)
	@test F1.S ≈ F2.S
	sgn = sign.(diag(F1.U'F2.U))
	@test F1.U.*sgn' ≈ F2.U
	@test F1.V.*sgn' ≈ F2.V
	@test F1.Vt.*sgn ≈ F2.Vt
	nothing
end


text:
using BenchmarkTools
using LinearAlgebra
using CUDA
using Test

free_parameters = 1000

#C = rand(ComplexF32, free_parameters, free_parameters)
C = rand(free_parameters, free_parameters)
C += C'
C_GPU = CuArray(convert(Matrix{Float32}, C))

@btime val, vec = eigen(C)
@btime CUDA.@sync val_GPU, vec_GPU = CUDA.CUSOLVER.syevd!('V', 'U', C_GPU)

#val, vec = eigen(C)
#val_GPU, vec_GPU = CUDA.CUSOLVER.syevd!('V', 'U', C_GPU)

#@test val ≈ Array(val_GPU)

println("Finished")


text:
using DoubleExponentialFormulas
using LinearAlgebra: norm
using Test

quades32 = QuadES(Float32)
quades64 = QuadES(Float64)
quadesBF = QuadES(BigFloat)

# Test Base.show works without error
let
    io = IOBuffer()
    @test (Base.show(io, MIME("text/plain"), quades32); true)
    @test (Base.show(io, MIME("text/plain"), quades64); true)
    @test (Base.show(io, MIME("text/plain"), quadesBF); true)
end

# Test [0, ∞] interval (with Exp-Sinh quadrature)
let
    f(x::AbstractFloat) = exp(-x)
    expect = 1

    I, E = quades32(f)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quades64(f)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quadesBF(f, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test [0, ∞] interval (with Exp-Sinh quadrature)
let
    f(x::AbstractFloat) = -exp(-x)*log(x)
    expect = BigFloat("5.77215664901532860606512090082e-1")  # γ constant

    I, E = quades32(f)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quades64(f)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quadesBF(f, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test [0, ∞] interval (with Exp-Sinh quadrature)
let
    f(x::AbstractFloat) = 2/(1 + x^2)
    expect = π

    I, E = quades32(f)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quades64(f)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quadesBF(f, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test [0, ∞] interval (with Exp-Sinh quadrature)
# Dirichlet integral
let
    f(x::AbstractFloat) = 2*sin(x)/x
    expect = π

    I, E = quades32(f)
    @test I isa Float32
    @test_broken I ≈ expect
    @test_skip E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quades64(f)
    @test I isa Float64
    @test_broken I ≈ expect
    @test_skip E ≤ sqrt(eps(typeof(I)))*norm(I)

    # FIXME: Something wrong. Infinite loop?
    # rtol = 1e-17
    # I, E = quadesBF(f, rtol=rtol)
    # @test I isa BigFloat
    # @test_broken isapprox(I, expect, rtol=10rtol)
    # @test_skip E ≤ rtol*norm(I)
end

# Test `f(xk)*wk = 0*Inf` case
let
    f(x::AbstractFloat) = x*exp(-x^2)
    expect = 1/BigFloat(2)

    I, E = quades32(f)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quades64(f)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quadesBF(f, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test non-scalar output with QuadES
let
    f(x::AbstractFloat) = [exp(-x), 2exp(-x)]
    expect = [1, 2]

    I, E = quades32(f)
    @test eltype(I) == Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    I, E = quades64(f)
    @test eltype(I) == Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    rtol = 1e-30
    I, E = quadesBF(f, rtol=rtol)
    @test eltype(I) == BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end


text:
using Plots: length
using Plots,DifferentialEquations,Measures

############################################### CSTR FUNCTION #################################################################
function cstr(du,u,p,t)

    Ca,T = u # States
    Ea,k₀,V,ρ,Cₚ,ΔH,UA,q,Caf,Tf,Tc = p # Parameters
    k(x) = k₀*exp(-Ea/8.314/x) # kinetic constant

    du[1]  = dCa = (q/V)*(Caf-Ca) - k(T)*Ca # derivatie of Ca
    du[2]  = dT = (q/V)*(Tf-T) + (-ΔH/ρ/Cₚ)*k(T)*Ca + (UA/V/ρ/Cₚ)*(Tc-T) # Derivative of T

end
##############################################################################################################################

############################################### PARAMETERS ###################################################################
Ea = 72750          # activation energy J/gmol
k₀ = 7.2e10         # Arrhenius rate constant 1/min
V = 100.0           # Volume [L]
ρ = 1000.0          # Density [g/L]
Cₚ = 0.239          # Heat capacity [J/g/K]
ΔH = -5.0e4         # Enthalpy of reaction [J/mol]
UA = 5.0e4          # Heat transfer [J/min/K]
q = 100.0           # Flowrate [L/min]
Caf = 1.0           # Inlet feed concentration [mol/L]
Tf = 350.0          # Inlet feed temperature [K]
Tc = 280.0          # Coolant temperature [K]
##############################################################################################################################

############################################### INITIAL COND #################################################################
Ca0 = 0.5
T0 = 350.0
##############################################################################################################################

############################################### SIMULATION ###################################################################
u0 = [Ca0 , T0] # Initial Condition Array
p = [Ea,k₀,V,ρ,Cₚ,ΔH,UA,q,Caf,Tf,Tc] # Parameters Array
tspan = (0.0 , 10.0) # Creating time span to solve the ODE from 0 to 10 min
prob = ODEProblem(cstr,u0,tspan,p) # Creating the ODEProblem
sol = solve(prob,Tsit5(),saveat=0.1) # Solving the problem and saving the results at each 0.1 min
##############################################################################################################################

############################################### PLOTTING SOLUTION ###############################################################
gr(size=(800,450)) # Backend choice
t = 0:0.1:10 # Time array
T_react = sol[2,:] # Getting Temperature solution vector
Ca_react = sol[1,:] # Getting Concentration solution vector
p_Ca = plot(t,Ca_react,title="Reactor's Concentration",ylabel="Ca (mol/L)",
            label=""); # Creating a plot to Concentration

p_T = plot(t,T_react,title="Reactor's Temperature",ylabel="Temperature (K)",
            xlabel="Time (min)",label=""); # Creating a plot to Temperature

transient_plot = plot(p_Ca,p_T,layout=(2,1)) # Adding both plots to a grid layout of 2 rows and 1 column
savefig(transient_plot,"Images/transient_plot.png") # Saving the transient_plot

phase_plot = plot(sol,vars=(1,2),ylabel="Temperature (K)",xlabel="Concentration (mol/L)",
                    title="Phase Space",label="") # Creating a plot from the States (Ca,T) a.k.a Phase Plot
savefig(phase_plot,"Images/phase_plot.png") # Saving the phase_plot

sub = plot(p_Ca,p_T,phase_plot,layout=@layout([[a;b] c]),left_margin=5mm)
savefig(sub,"Images/all_in_one.png") # Saving all plots in one
##############################################################################################################################

############################################### ANIMATION PLOT ###############################################################
# This function simulates the CSTR for a giving Coolant Temerature.
function sim_cstr(Tc)
    u0 = [Ca0 , T0] # Initial Condition Array
    p = [Ea,k₀,V,ρ,Cₚ,ΔH,UA,q,Caf,Tf,Tc] # Parameters Array
    tspan = (0.0 , 10.0) # Creating time span to solve the ODE from 0 to 10 min
    prob = ODEProblem(cstr,u0,tspan,p) # Creating the ODEProblem
    sol = solve(prob,Tsit5(),saveat=0.1) # Solving the problem and saving the results at each 0.1 min
    return sol
end

Tcs = 280:0.5:350 # Array of Coolant Temperatures
sols = sim_cstr.(Tcs); # Solving for each Coolant Temperature in Tcs

# Creating animation using the animate macro
anim = @animate for i ∈ 1:length(Tcs)
    # Plots each solution in the sols array (these makes a new plot for each solution, thus older solutions don't appear)
    plot(
        plot(sols[i].t,sols[i][1,:],title="Reactor's Concentration (Tc = "*string(Tcs[i])*"K)",ylabel="Ca (mol/L)",label=:none),
        plot(sols[i].t,sols[i][2,:],title="Reactor's Temperature (Tc = "*string(Tcs[i])*"K)",ylabel="Temperature (K)",xlabel="Time (min)",label=:none),
        plot(sols[i],vars=(1,2),ylabel="Temperature (K)",xlabel="Concentration (mol/L)",title="Phase Space (Tc = "*string(Tcs[i])*"K)",label=:none),
        layout=@layout([[a;b] c]),left_margin=10mm
    )
    scatter!((Ca0,T0),sp=3,label="Initial Condition",markercolor=:black) # Adds a point in the phase plot regarding the initial condition
end
gif(anim, "Images/cstr_coolant_effect.gif", fps = 10) # create a gif using the anim variable where we've stored our frames
##############################################################################################################################

function derivatives_cstr(Ca,T,Tc)

    k(x) = k₀*exp(-Ea/8.314/x) # kinetic constant

    dCa = (q/V)*(Caf-Ca) - k(T)*Ca # derivatie of Ca
    dT = (q/V)*(Tf-T) + (-ΔH/ρ/Cₚ)*k(T)*Ca + (UA/V/ρ/Cₚ)*(Tc-T) # Derivative of T

    return [dCa,dT]

end

meshgrid(x,y) = (repeat(x, outer=length(y)),repeat(y, inner=length(x)))
Ca_test,T_test = meshgrid(0:0.1:1,280:10:350)
dv = derivatives_cstr.(Ca_test,T_test,305)
dca,dt = hcat(dv...)'[:,1],hcat(dv...)'[:,2]
scale = 1
dca,dt = dca./sqrt.(dca.^2+dt.^2),dt./sqrt.(dca.^2+dt.^2)
quiver(Ca_test,T_test,quiver=scale.*(dca,dt),color=:inferno)
scatter!([Ca0],[T0])


text:
# Functions for contracting two arrays, both for general StridedArrays as well as
# for SymArrays. we only have a few specialized functions that we have needed up
# to now, but carefully optimize each of them

using TensorOperations
using LinearAlgebra
using Requires

# Array[i]*SymArray[(i,j,k)]
# indices 1, 2, and 3 are exchangeable here
function contract(A::StridedVector{T},S::SymArray{(3,),U},n::Union{Val{1},Val{2},Val{3}}) where {T,U}
    TU = promote_type(T,U)
    @assert size(S,1) == length(A)
    res = SymArray{(2,),TU}(size(S,1),size(S,2))
    contract!(res,A,S,n)
end

# We know that $j\leq k$ (because $R$ is itself exchange symmetric)
# \begin{align}
# R_{jk} &= \sum_{i=1}^N g_i S_{ijk}
# \end{align}

# Matrix elements represented by $S_{ijk}$:
# \begin{equation}
# \begin{cases}
# S_{ijk}, S_{ikj}, S_{jik}, S_{jki}, S_{kij}, S_{kji} & i<j<k\\
# S_{ijk}, S_{jik}, S_{jki} & i<j=k\\
# S_{ijk}, S_{ikj}, S_{kij} & i=j<k\\
# S_{ijk} & i=j=k
# \end{cases}
# \end{equation}
# We only need to take the contributions that show up for each $R_{jk}$

# Array[i]*SymArray[(i,j,k)]
# indices 1, 2, and 3 are exchangeable here
function contract!(res::SymArray{(2,),TU}, A::StridedVector{T}, S::SymArray{(3,),U}, n::Union{Val{1},Val{2},Val{3}}) where {T,U,TU}
    # only loop over S once, and put all the values where they should go
    # R[j,k] = sum_i A[i] B[i,j,k]
    # S[i,j,k] with i<=j<=k represents the 6 (not always distinct) terms: Bijk, Bikj, Bjik, Bjki, Bkij, Bkji
    # since R[j,k] is also exchange symmetric, we only need to calculate j<=k
    # commented out terms below are the ones where j>k
    @assert size(S,1) == length(A)
    @assert size(S,1) == size(res,1)
    res.data .= 0
    @inbounds for (v,inds) in zip(S.data,CartesianIndices(S))
        i,j,k = Tuple(inds)
        res[j,k] += v*A[i]
        if i==j==k
            continue
        elseif i==j # i=j < k
            # exch. terms: Sikj and Skij
            # res[k,j] += v*A[i] i2>i3
            res[i,j] += v*A[k]
        elseif j==k # i < j=k
            # exch. terms Sjik, Sjki
            res[i,k] += v*A[j]
            # res[k,i] += v*A[j] i3>i2
        else # i<j<k
            # exch. terms Sikj, Sjik, Sjki, Skij, Skji
            # res[k,j] += v*A[i] i3>i2
            res[i,k] += v*A[j]
            # res[k,i] += v*A[j] i3>i2
            res[i,j] += v*A[k]
            # res[j,i] += v*A[k] i3>i2
        end
    end
    res
end

# Array[k]*SymArray[(i,j),k]
function contract(A::StridedVector{T},S::SymArray{(2,1),U},n::Val{3}) where {T,U}
    TU = promote_type(T,U)
    sumsize = length(A)
    @assert sumsize == size(S,3)
    res = SymArray{(2,),TU}(size(S,1),size(S,2))
    contract!(res,A,S,n)
end

# Array[k]*SymArray[(i,j),k]
function contract!(res::SymArray{(2,),TU},A::StridedVector{T},S::SymArray{(2,1),U},::Val{3}) where {T,U,TU}
    # use that S[(i,j),k] == S[I,k] (i.e., the two symmetric indices act like a "big" index)
    mul!(res.data,reshape(S.data,:,length(A)),A)
    res
end

# Array[i]*SymArray[(i,j),k)]
# since indices 1 and 2 are exchangeable here, use this
function contract(A::StridedVector{T},S::SymArray{(2,1),U},n::Union{Val{1},Val{2}}) where {T,U}
    TU = promote_type(T,U)
    @assert size(S,1) == length(A)
    # the result is a normal 2D array
    res = Array{TU,2}(undef,size(S,2),size(S,3))
    contract!(res,A,S,n)
end

# Array[i]*SymArray[(i,j),k]
# since indices 1 and 2 are exchangeable here, use this
function contract!(res::StridedArray{TU,2},A::StridedVector{T},S::SymArray{(2,1),U},::Union{Val{1},Val{2}}) where {T,U,TU}
    # only loop over S once, and put all the values where they should go
    @assert size(A,1) == size(S,1)
    @assert size(res,1) == size(S,1)
    @assert size(res,2) == size(S,3)
    res .= zero(TU)
    @inbounds for (v,inds) in zip(S.data,CartesianIndices(S))
        i1,i2,i3 = Tuple(inds)
        res[i2,i3] += v*A[i1]
        # if i1 != i2, we have to add the equal contribution from S[i2,i1,i3]
        if i1 != i2
            res[i1,i3] += v*A[i2]
        end
    end
    res
end

# Array[i]*SymArray[(i,j)]
# this is symmetric in i1 and i2
function contract(A::StridedVector{T},S::SymArray{(2,),U},n::Union{Val{1},Val{2}}) where {T,U}
    TU = promote_type(T,U)
    @assert size(S,1) == length(A)
    # the result is a normal 1D vector
    res = Vector{TU}(undef,size(S,2))
    contract!(res,A,S,n)
end

# Array[i]*SymArray[(i,j)]
# this is symmetric in i1 and i2
function contract!(res::StridedVector{TU},A::StridedVector{T},S::SymArray{(2,),U},::Union{Val{1},Val{2}}) where {T,U,TU}
@assert size(A,1) == size(S,1)
    @assert size(res,1) == size(S,1)
    res .= zero(TU)
    # only loop over S once, and put all the values where they should go
    for (v,inds) in zip(S.data,CartesianIndices(S))
        i1,i2 = Tuple(inds)
        res[i2] += v*A[i1]
        # if i1 != i2, we have to add the contribution from S[i2,i1]
        if i1 != i2
            res[i1] += v*A[i2]
        end
    end
    res
end

# Array[i_n]*Array[i1,i2,i3,...,iN]
function contract(A::StridedVector{T},B::StridedArray{U,N},::Val{n}) where {T,U,N,n}
    TU = promote_type(T,U)
    @assert 1 <= n <= N

    resdims = size(B)[1:N .!= n]
    res = similar(B,TU,resdims)

    A = convert(AbstractArray{TU},A)
    B = convert(AbstractArray{TU},B)
    contract!(res,A,B,Val{n}())
end

mygemv!(args...) = BLAS.gemv!(args...)

function _contract_middle!(res,A,B)
    @inbounds for k=1:size(B,3)
        mul!(@view(res[:,k]), @view(B[:,:,k]), A)
    end
end

function __init__()
    @require CuArrays = "3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
        mygemv!(tA,alpha,A::CuArrays.CuArray,args...) = CuArrays.CUBLAS.gemv!(tA,alpha,A,args...)
        _contract_middle!(res::CuArrays.CuArray,A,B) = error("please install CuTensorOperations.jl to allow this operation with CuArrays")
    end
    @require CuTensorOperations = "01a401d6-3e09-11e9-13be-6f1bb8674acf" begin
        _contract_middle!(res::CuArrays.CuArray,A,B) = (@tensor res[i,k] = B[i,j,k] * A[j])
    end
end

# Array[i_n]*Array[i1,i2,i3,...,iN]
function contract!(res::StridedArray{TU},A::StridedVector{TU},B::StridedArray{TU,N},::Val{n}) where {TU,N,n}
    nsum = length(A)
    @assert size(B,n) == nsum
    @assert ndims(res)+1 == ndims(B)
    ii = 0
    for jj = 1:ndims(B)
        jj==n && continue
        ii += 1
        @assert size(B,jj) == size(res,ii)
    end

    if n==1      # A[i]*B[i,...]
        mygemv!('T',one(TU),reshape(B,nsum,:),A,zero(TU),vec(res))
    elseif n==N  # B[...,i]*A[i]
        mygemv!('N',one(TU),reshape(B,:,nsum),A,zero(TU),vec(res))
    else
        rightsize = prod(size(B,i) for i=n+1:N)
        Br = reshape(B,:,nsum,rightsize)
        resr = reshape(res,:,rightsize)
        _contract_middle!(resr,A,Br)
    end
    res
end


text:
using Random

gsize = 100
reads = 400
genome = [0 for n in 1:gsize]
for i = 1:reads
	genome[rand(1:gsize)] += 1
end
println(genome)

max_stars_repo_name: KorfLab/julie

max_issues_repo_name: KorfLab/julie

max_forks_repo_name: KorfLab/julie


text:
using JLD2
using VUMPS:SquareVUMPSRuntime, ALCtoAC

function observable(model, Ni, Nj, atype, folder, D, χ, tol=1e-10, maxiter=10)
    observable_log = folder*"/$(model)_$(Ni)x$(Nj)/D$(D)_χ$(χ)_observable.log"
    # if isfile(observable_log)
    #     println("load observable from $(observable_log)")
    #     f = open(observable_log, "r" )
    #     occ,doubleocc = parse.(Float64,split(readline(f), "   "))
    #     close(f)
    # else
        ipeps, key = init_ipeps(model; Ni = Ni, Nj = Nj, atype = atype, folder = folder, D=D, χ=χ, tol=tol, maxiter= maxiter)
        folder, model, Ni, Nj, atype, D, χ, tol, maxiter = key
        T = reshape([parity_conserving(ipeps[:,:,:,:,:,ABBA(i)]) for i = 1:Ni*Nj], (Ni, Nj))
        b = reshape([zeros(D^2,D^2,D^2,D^2) for i = 1:Ni*Nj], (Ni, Nj))
        # b = reshape([permutedims(bulk(T[i]),(2,3,4,1)) for i = 1:Ni*Nj], (Ni, Nj))
        op = reshape([permutedims(bulkop(T[i]),(2,3,4,1,5,6)) for i = 1:Ni*Nj], (Ni, Nj))

        chkp_file_obs = folder*"/$(model)_$(Ni)x$(Nj)/obs_D$(D^2)_χ$(χ).jld2"
        FLo, FRo = load(chkp_file_obs)["env"]
        chkp_file_up = folder*"/$(model)_$(Ni)x$(Nj)/up_D$(D^2)_χ$(χ).jld2"
        rtup = SquareVUMPSRuntime(b, chkp_file_up, χ)
        FLu, FRu, ALu, ARu, Cu = rtup.FL, rtup.FR, rtup.AL, rtup.AR, rtup.C
        chkp_file_down = folder*"/$(model)_$(Ni)x$(Nj)/down_D$(D^2)_χ$(χ).jld2"
        rtdown = SquareVUMPSRuntime(b, chkp_file_down, χ)
        ALd,ARd,Cd = rtdown.AL,rtdown.AR,rtdown.C

        ACu = ALCtoAC(ALu,Cu)
        ACd = ALCtoAC(ALd,Cd)

        hocc = atype([0.0 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 2])
        hdoubleocc = atype([0.0 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 1])
        Nzup = atype([0.0 0 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 1])
        Nzdn = atype([0.0 0 0 0; 0 0 0 0; 0 0 1 0; 0 0 0 1])
        Nxup = atype([0.0 0 0 0; 0 1  1 0; 0  1 1 0; 0 0 0 2]./2)
        Nxdn = atype([0.0 0 0 0; 0 1 -1 0; 0 -1 1 0; 0 1 1 2]./2)
        Nyup = atype([0.0 0 0 0; 0 1 -1im 0; 0 1im 1 0; 0 0 0 2]./2)
        Nydn = atype([0.0 0 0 0; 0 1 1im 0; 0 -1im 1 0; 0 1 1 2]./2)
        U = atype([1 0 0 0;0 0 1 0;0 -1 0 0;0 0 0 1])
        occ = 0
        doubleocc = 0
        for j = 1:Nj, i = 1:Ni
            println("==========$i $j==========")
            ir = Ni + 1 - i
            ρ = ein"(((adf,abc),dgebpq),fgh),ceh -> pq"(FLo[i,j],ACu[i,j],op[i,j],ACd[ir,j],FRo[i,j])
            # if (i,j) in [(2,1),(1,2)]
            #     ρ = U' * ρ * U
            # end
            Occ = ein"pq,pq -> "(ρ,hocc)
            DoubleOcc = ein"pq,pq -> "(ρ,hdoubleocc)
            NNzup = ein"pq,pq -> "(ρ,Nzup)
            NNzdn = ein"pq,pq -> "(ρ,Nzdn)
            NNxup = ein"pq,pq -> "(ρ,Nxup)
            NNxdn = ein"pq,pq -> "(ρ,Nxdn)
            NNyup = ein"pq,pq -> "(ρ,Nyup)
            NNydn = ein"pq,pq -> "(ρ,Nydn)
            n = ein"pp -> "(ρ)
            occ += Array(Occ)[]/Array(n)[]
            doubleocc += Array(DoubleOcc)[]/Array(n)[]
            println("N = $(Array(Occ)[]/Array(n)[])")
            println("DN = $(Array(DoubleOcc)[]/Array(n)[])")
            println("Nz↑ = $(Array(NNzup)[]/Array(n)[])")
            println("Nz↓ = $(Array(NNzdn)[]/Array(n)[])")
            println("Nx↑ = $(Array(NNxup)[]/Array(n)[])")
            println("Nx↓ = $(Array(NNxdn)[]/Array(n)[])")
            println("Ny↑ = $(Array(NNyup)[]/Array(n)[])")
            println("Ny↓ = $(Array(NNydn)[]/Array(n)[])")
            println("{$(real(Array(NNzup)[]/Array(n)[])),$(real(Array(NNzdn)[]/Array(n)[])),$(real(Array(NNxup)[]/Array(n)[])),$(real(Array(NNxdn)[]/Array(n)[])),$(real(Array(NNyup)[]/Array(n)[])),$(real(Array(NNydn)[]/Array(n)[]))}")
        end

        occ = real(occ)/Ni/Nj
        doubleocc = real(doubleocc)/Ni/Nj
        message = "$occ   $doubleocc\n"
        logfile = open(observable_log, "a")
        write(logfile, message)
        close(logfile)
    # end
    return occ, doubleocc
end


text:
#   Unit examples of the FourierAnalysis Package for julia language
#   v 0.0.1 - last update 24th of September 2019
#
#   MIT License
#   Copyright (c) 2019, Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This example show how to compute filter banks.

#   ~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~  #
#                                                                             #
#   ~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~  #

####################
# test Filters

using Plots, DSP, FourierAnalysis
x=sinusoidal(2, 8, 128, 512, 0)
f, Y=filterbank(x, 128)
plot(x)
plot!((Y[:, 7]))

# Show amplitude distorsions at the edges and the linear phase response
plot!((x-Y[:, 7]))

f, Y=filterbank(x, 128; filtkind=Butterworth(4))

x2=randn(512)+sinusoidal(2, 8, 128, 512, 0)
f, Y=filterbank(x2, 128)
plot(x)
plot!(x2)
plot!((Y[:, 7]))
plot!((x-Y[:, 7]))


text:
@everywhere begin
    import AARMBEMC.AARMBEMCheat.heattransfer
end

heattransfer(ARGS)


text:
module SphericalHarmonics

using IrrationalConstants
using SphericalHarmonicModes
using SphericalHarmonicArrays
using Printf
using StaticArrays
using Setfield
using SpecialFunctions: loggamma

export computeYlm, computeYlm!
export computePlmcostheta, computePlmcostheta!
export computePlmx, computePlmx!

using Base: @propagate_inbounds

include("irrationals.jl")

abstract type HarmonicType end
struct RealHarmonics <: HarmonicType end
struct ComplexHarmonics <: HarmonicType end

"""
    SphericalHarmonicsCache

Preallocate arrays of associated Legendre polynomials and spherical harmonics.
Such an object may be allocated using [`cache`](@ref).
"""
mutable struct SphericalHarmonicsCache{T, C, PLM<:AbstractVector{T}, YLM<:AbstractVector, SHT}
    C :: C
    P :: PLM
    Y :: YLM
    SHType :: SHT
end

struct AssociatedLegendrePolynomials{T, PLM<:AbstractVector{T}} <: AbstractVector{T}
    cosθ :: T
    lmax :: Int
    P :: PLM
    initialized :: Bool
end
AssociatedLegendrePolynomials(cosθ, lmax, P::AbstractVector, initialized::Bool) =
    AssociatedLegendrePolynomials(convert(eltype(P), cosθ), Int(lmax), P, initialized)
Base.parent(P::AssociatedLegendrePolynomials) = P.P
Base.size(P::AssociatedLegendrePolynomials) = size(parent(P))
Base.axes(P::AssociatedLegendrePolynomials) = axes(parent(P))
Base.IndexStyle(::Type{<:AssociatedLegendrePolynomials{<:Any, PLM}}) where {PLM} = IndexStyle(PLM)
@propagate_inbounds Base.getindex(P::AssociatedLegendrePolynomials, I...) = getindex(parent(P), I...)
function Base.summary(io::IO, P::AssociatedLegendrePolynomials)
    print(io, "$(length(P))-element normalized AssociatedLegendrePolynomials{$(eltype(P))} for lmax = $(Int(P.lmax))")
    if P.initialized
        print(io, " and cosθ = ")
        @printf io "%.4g" P.cosθ
    else
        print(io, " (uninitialized)")
    end
end

function SphericalHarmonicsCache(T::Type, lmax::Int, ::Type{m_range}, SHType) where {m_range}
    C = compute_coefficients(T, lmax)
    P = AssociatedLegendrePolynomials(T(0), lmax, allocate_p(T, lmax), false)
    Y = allocate_y(eltypeY(T, SHType), lmax, m_range)
    SphericalHarmonicsCache(C, P, Y, SHType)
end
SphericalHarmonicsCache(lmax::Int, T::Type, SHType) = SphericalHarmonicsCache(Float64, lmax, T, SHType)

function Base.show(io::IO, S::SphericalHarmonicsCache)
    print(io, "$SphericalHarmonicsCache($(eltypeP(S)), $(_lmax(S)), $(_mrange_basetype(S)), $(S.SHType))")
end

"""
    cache([T::Type = Float64], lmax, [m_range = SphericalHarmonicModes.FullRange], [SHType = SphericalHarmonics.ComplexHarmonics()])

Allocate arrays to evaluate associated Legendre polynomials and spherical harmonics.
The returned object may be passed to [`computePlmcostheta!`](@ref) and [`computeYlm!`](@ref).
The coefficients are cached and need not be recomputed.

# Examples
```jldoctest cache
julia> S = SphericalHarmonics.cache(1);

julia> computePlmcostheta!(S, pi/3, 1)
3-element normalized AssociatedLegendrePolynomials{Float64} for lmax = 1 and cosθ = 0.5:
  0.3989422804014327
  0.34549414947133555
 -0.4231421876608172

julia> computeYlm!(S, pi/3, pi/4, 1)
4-element SHArray(::Vector{ComplexF64}, (ML(0:1, -1:1),)):
   0.2820947917738782 + 0.0im
  0.21157109383040865 - 0.2115710938304086im
  0.24430125595146002 + 0.0im
 -0.21157109383040865 - 0.2115710938304086im
```

Choosing a new `lmax` in `computePlmcostheta!` expands the cache if necessary.

```jldoctest cache
julia> computePlmcostheta!(S, pi/3, 2)
6-element normalized AssociatedLegendrePolynomials{Float64} for lmax = 2 and cosθ = 0.5:
  0.3989422804014327
  0.34549414947133555
 -0.4231421876608172
 -0.11150775725954817
 -0.4730873478787801
  0.40970566147202964
```
"""
cache(args...; m_range = FullRange, SHType = ComplexHarmonics()) = cache(args..., m_range, SHType)
cache(::Type{T}, lmax::Int, ::Type{m_range}, SHType::HarmonicType = ComplexHarmonics()) where {T,m_range} = SphericalHarmonicsCache(T, lmax, m_range, SHType)
cache(::Type{T}, lmax::Int, SHType::HarmonicType = ComplexHarmonics()) where {T} = SphericalHarmonicsCache(T, lmax, FullRange, SHType)
cache(lmax::Int, ::Type{m_range}, SHType::HarmonicType = ComplexHarmonics()) where {m_range} = SphericalHarmonicsCache(lmax, m_range, SHType)
cache(lmax::Int, SHType::HarmonicType = ComplexHarmonics()) = SphericalHarmonicsCache(lmax, FullRange, SHType)

const SHMRange = Union{Type{FullRange}, Type{ZeroTo}}

sizeP(maxDegree::Int) = div((maxDegree + 1) * (maxDegree + 2), 2)
sizeY(maxDegree::Int, ::Type{FullRange} = FullRange) = (maxDegree + 1) * (maxDegree + 1)
sizeY(maxDegree::Int, ::Type{ZeroTo}) = sizeP(maxDegree)

"""
    SphericalHarmonics.allocate_p([T::Type = Float64], lmax::Integer)

Allocate an array large enough to store an entire set of Associated Legendre
Polynomials ``\\bar{P}_ℓ^m(x)`` of maximum degree ``ℓ``.
"""
allocate_p(T::Type, lmax::Integer) = SHArray{T,1}(undef, (ML(ZeroTo(lmax), ZeroTo),))
allocate_p(lmax::Integer) = allocate_p(Float64, lmax)

"""
    SphericalHarmonics.allocate_y([T::Type = ComplexF64], lmax::Integer)

Allocate an array large enough to store an entire set of spherical harmonics
``Y_{ℓ,m}(θ,ϕ)`` of maximum degree ``ℓ``.
"""
allocate_y(T::Type, lmax::Integer, m_range = FullRange) = SHArray{T,1}(undef, (ML(ZeroTo(lmax), m_range),))
allocate_y(lmax::Integer, m_range = FullRange) = allocate_y(ComplexF64, lmax, m_range)

# This method uses the somewhat complicated promote_type mechanism instead of complex(R), as
# some types such as Dual numbers currently do not have complex(R) defined for them.
# We effectively evaluate typeof(zero(R) * im) instead.
eltypeY(::Type{R}, ::ComplexHarmonics) where {R} = promote_type(R, complex(Bool))
eltypeY(::Type{R}, ::RealHarmonics) where {R} = R

# define accessor methods that may be used by wrappers
eltypeP(S::SphericalHarmonicsCache{T}) where {T} = T
eltypeY(S::SphericalHarmonicsCache{<:Any,<:Any,<:Any,Y}) where {Y} = eltype(Y)
getP(S::SphericalHarmonicsCache) = S.P
getY(S::SphericalHarmonicsCache) = S.Y

function _mrange_basetype(S::SphericalHarmonicsCache)
    Y = getY(S)
    modes = first(SphericalHarmonicArrays.shmodes(Y))
    _mrange_basetype(m_range(modes))
end
_mrange_basetype(::FullRange) = FullRange
_mrange_basetype(::ZeroTo) = ZeroTo

function _lmax(S::SphericalHarmonicsCache)
    modes = first(SphericalHarmonicArrays.shmodes(S.C))
    maximum(l_range(modes))
end

@doc raw"""
    SphericalHarmonics.compute_coefficients(lmax)

Precompute coefficients ``a_ℓ^m`` and ``b_ℓ^m`` for all ``2 ≤ ℓ ≤ ℓ_\mathrm{max}`` and ``0 ≤ m ≤ ℓ-2``.

    SphericalHarmonics.compute_coefficients(lmax, m)

Precompute coefficients ``a_ℓ^m`` and ``b_ℓ^m`` for all ``|m| + 2 ≤ ℓ ≤ ℓ_\mathrm{max}`` and the specified ``m``.
"""
compute_coefficients(lmax::Integer) = compute_coefficients(typeof(sqrt(lmax)), lmax)
compute_coefficients(lmax::Integer, ::Nothing) = compute_coefficients(lmax)
compute_coefficients(lmax::Integer, m::Integer) = compute_coefficients(typeof(sqrt(lmax)), lmax, m)
function compute_coefficients(T::Type, lmax::Integer)
    @assert lmax >= 0 "degree must be non-negative"

    shmodes = ML(ZeroTo(Int(lmax)), ZeroTo)
    A = zeros(SVector{2,T}, length(shmodes))
    coeff = SHArray(A, (shmodes,))

    @inbounds _compute_coefficients!(coeff, 2, Int(lmax))
    return coeff
end
compute_coefficients(T::Type, lmax::Integer, ::Nothing) = compute_coefficients(T, lmax)
function compute_coefficients(T::Type, lmax::Integer, m::Integer)
    @assert lmax >= 0 "degree must be non-negative"
    @assert m >= 0 "m must be non-negative"

    shmodes = ML(ZeroTo(Int(lmax)), SingleValuedRange(m))
    coeff = zeros(SVector{2,T}, shmodes)

    @inbounds _compute_coefficients!(coeff, abs(Int(m)) + 2, Int(lmax), abs(Int(m)), abs(Int(m)))
    return coeff
end
function compute_coefficients!(S::SphericalHarmonicsCache, lmax::Integer)
    @assert lmax >= 0 "degree must be >= 0"
    if lmax > _lmax(S)
        shmodesP = ML(ZeroTo(Int(lmax)), ZeroTo)
        shmodesY = ML(ZeroTo(Int(lmax)), _mrange_basetype(S))
        A = parent(S.C)
        resize!(A, length(shmodesP))
        coeff = SHArray(A, (shmodesP,))

        @inbounds _compute_coefficients!(coeff, Int(_lmax(S))+1, Int(lmax))

        S.C = coeff
        P_len = length(S.P)
        P_new = resize!(parent(parent(S.P)), length(shmodesP))
        P_new[P_len+1:end] .= zero(eltype(P_new))
        P_new_SA = SHArray(P_new, (shmodesP,))
        P = S.P
        P = @set P.P = P_new_SA
        S.P = P

        Y_len = length(S.Y)
        Y_new = resize!(parent(S.Y), length(shmodesY))
        Y_new[Y_len+1:end] .= zero(eltype(Y_new))
        S.Y = SHArray(Y_new, (shmodesY,))
    end
    return S.C
end

@propagate_inbounds function _compute_coefficients!(coeff, lmin, lmax, m_min = 0, m_max = lmax - 2)
    T = eltype(eltype(coeff))
    for l in lmin:lmax
        # pre-compute certain terms to improve performance
        invlm1 = 1/T(l-1)
        invl = 1/T(l)
        Anum = 4 - invl^2
        Bden = 1/√(4 - invlm1^2)
        for m in max(0, m_min):min(l-2, m_max)
            coeff[(l, m)] = SVector{2}(√(Anum / (1 - (m * invl)^2)), -√(1 - (m * invlm1)^2) * Bden)
        end
    end
end

function checksize(sz, minsize)
    @assert sz >= minsize "array needs to have a minimum size of $minsize, received size $sz"
end
function checksizesP(P, lmax)
    lmax >= 0 || throw(ArgumentError("lmax = $lmax does not correspond to a valid mode"))
    checksize(length(P), sizeP(Int(lmax)))
end

@propagate_inbounds readcoeffs(coeff::SHArray, T::Type, l, m) = coeff[(l,m)]

function readcoeffs(coeff::Nothing, T::Type, l, m)
    @assert abs(m) <= l - 2 "m must be <= l - 2"
    invl = 1/T(l)
    invlm1 = 1/T(l-1)
    alm = √(T(4 - invl^2)/T(1 - (m * invl)^2))
    blm = -√(T(1 - (m * invlm1)^2 )/T(4 - invlm1^2))
    return alm, blm
end

######################################################################################
# Associated Legendre polynomials
######################################################################################

_promotetype(costheta, sintheta, coeff::Nothing) = float(promote_type(typeof(costheta), typeof(sintheta)))
_promotetype(costheta, sintheta, coeff::AbstractArray) = float(reduce(promote_type, (typeof(costheta), typeof(sintheta), eltype(eltype(coeff)))))

_viewmodes(P, modes) = @view P[(firstindex(P) - 1) .+ (1:length(modes))]

function _wrapSHArray(P, modes)
    SHArray(_viewmodes(P, modes), (modes,))
end
function _wrapSHArray(P, lmax, m_range)
    modes = ML(ZeroTo(Int(lmax)), m_range)
    _wrapSHArray(P, modes)
end

abstract type PLMnorm end
struct Orthonormal <: PLMnorm cs::Bool end
struct Unnormalized <: PLMnorm cs::Bool end
struct LMnorm <: PLMnorm cs::Bool end # Limpanuparb-Milthorpe norm as used in the paper

(::Type{T})(; csphase = true) where {T<:PLMnorm} = T(csphase)

# Condon-Shortey phase is included by default. This may be disabled by passing
# the flag to the normalization specifier
csphase(norm, m) = !norm.cs ? (-1)^m : 1
function invnorm(l, m, norm::Unnormalized)
    f = exp(-(log(2l+1) - log2π + loggamma(l-m+1) - loggamma(l+m+1))/2)
    f * csphase(norm, m)
end
invnorm(l, m, norm::Orthonormal) = (f = sqrtπ; f * csphase(norm, m))
invnorm(l, m, norm::LMnorm) = (f = true; f * csphase(norm, m))

m_l(modes, l, ::Nothing) =  m_range(modes, l)
m_l(modes, l, m) =  intersect(m_range(modes, l), m)
function normalize!(P, norm, l, m = nothing)
    norm isa LMnorm && norm.cs && return P # short-circuit the standard case
    modes = LM(l, m === nothing ? ZeroTo : m)
    Plm = _wrapSHArray(P, maximum(l), ZeroTo)
    T = eltype(Plm)
    for l in l_range(modes)
        for m in m_l(modes, l, m)
            Plm[(l,m)] *= invnorm(T(l), T(m), norm)
        end
    end
    return P
end

@propagate_inbounds function _computePlmcostheta!(P, costheta, sintheta, l_range,
    coeff::Union{AbstractArray,Nothing} = nothing; norm::PLMnorm = LMnorm())

    lmin, lmax = map(Int, extrema(l_range))
    @assert lmin >= 0 "degree must be non-negative"

    T = _promotetype(costheta, sintheta, coeff)
    Plm = _wrapSHArray(P, lmax, ZeroTo)

    if lmin == 0
        Plm[(0, 0)] = invsqrt2π
    end

    if (lmax > 0)
        P11 = -(_sqrt3by4pi * sintheta)
        temp = T(P11)
        if lmin <= 1
            Plm[(1, 0)] = _sqrt3by2pi * costheta
            Plm[(1, 1)] = P11
        end

        # compute temp without losing precision
        for l in 2:lmin - 1
            temp *= -√(T(2l + 1)/2l) * sintheta
        end

        for l in max(2, lmin):lmax
            for m in 0:(l-2)
                alm, blm = readcoeffs(coeff, T, l, m)
                Plm[(l,m)] = alm * (costheta * Plm[(l - 1, m)] + blm * Plm[(l - 2, m)])
            end
            Plm[(l, l - 1)] = costheta * √(T(2l + 1)) * temp
            temp *= -√(T(2l + 1)/2l) * sintheta
            Plm[(l, l)] = temp
        end
    end

    normalize!(P, norm, l_range)
    return P
end

@propagate_inbounds function _computePlmcostheta!(P, costheta, sintheta, l_range, m::Integer,
    coeff::Union{AbstractArray,Nothing} = nothing; norm::PLMnorm = LMnorm())

    lmin, lmax = map(Int, extrema(l_range))
    @assert lmin == 0 "minimum degree must be zero"
    # @assert lmin >= 0 "degree must be non-negative"
    0 <= m <= lmax || throw(ArgumentError("m = $m does not satisfy 0 ≤ m ≤ lmax = $lmax"))

    T = _promotetype(costheta, sintheta, coeff)

    Plm = _wrapSHArray(P, lmax, ZeroTo)

    if m == 0 && lmin == 0
        Plm[(0, 0)] = invsqrt2π
    end
    temp = T(invsqrt2π)

    if (lmax > 0)
        P11 = -(_sqrt3by4pi * sintheta)
        if lmin <= 1 && m > 0
            temp = T(P11)
            Plm[(1, 1)] = temp
        end

        # This is necessary for lmin > 2
        # for l in 2:lmin - 1
        #     temp *= -√(T(2l + 1)/2l) * sintheta
        # end

        # Compute Pmm using recursion over l and m
        for mi in max(lmin, 2):Int(m)
            temp *= -√(T(2mi + 1)/2mi) * sintheta
            Plm[(mi, mi)] = temp
        end

        if m == lmax
            normalize!(Plm, norm, l_range, m)
            return P
        end

        Plm[(m + 1, m)] = √(T(2m + 3)) * costheta * temp

        # Compute Plm using recursion over l at a fixed m
        for l in max(lmin, Int(m) + 2):lmax
            alm, blm = readcoeffs(coeff, T, l, m)
            Plm[(l,m)] = alm * (costheta * Plm[(l - 1, m)] + blm * Plm[(l - 2, m)])
        end
    end

    normalize!(Plm, norm, l_range, m)
    return P
end

@propagate_inbounds function _computePlmcostheta!(P, costheta, sintheta, l_range,
        m::Nothing, coeff::Union{AbstractArray,Nothing}; norm::PLMnorm = LMnorm())
    _computePlmcostheta!(P, costheta, sintheta, l_range, coeff; norm = norm)
end

@propagate_inbounds function Plmrecursion(l, m, costheta, Pmm, Pmp1m, coeff, T)
    Plm2_m, Plm1_m = Pmm, Pmp1m

    # Compute Plm using recursion over l at a fixed m, starting from l = m + 2
    for li in Int(m) + 2:Int(l)
        alm, blm = readcoeffs(coeff, T, li, m)

        Pl_m = alm * (costheta * Plm1_m + blm * Plm2_m)

        Plm2_m, Plm1_m = Plm1_m, Pl_m
    end

    return Plm1_m
end

@propagate_inbounds function _computePlmcostheta(costheta, sintheta, l, m, coeff; norm::PLMnorm = LMnorm())
    0 <= m <= l || throw(ArgumentError("m = $m does not satisfy 0 ≤ m ≤ l = $l"))

    T = _promotetype(costheta, sintheta, coeff)

    Pmm = T(invsqrt2π)

    if m > 0
        Pmm = -T(_sqrt3by4pi * sintheta)
    end

    # Compute Pmm using recursion over l and m
    for mi in 2:Int(m)
        Pmm = -√(1 + 1 / 2T(mi)) * sintheta * Pmm
    end

    if m == l
        Pmm *= invnorm(T(l), T(m), norm)
        return Pmm
    end

    Pmp1m = √(T(2m + 3)) * costheta * Pmm

    # Recursion at a constant m to compute Pl,m from Pm,m and Pm+1,m
    Plm = Plmrecursion(l, m, costheta, Pmm, Pmp1m, coeff, T)
    Plm *= invnorm(T(l), T(m), norm)
    return Plm
end

function zeroP!(P, lmax)
    checksizesP(P, lmax)
    N = sizeP(Int(lmax))
    Pv = @view P[firstindex(P) .+ (0:N-1)]
    fill!(Pv, zero(eltype(Pv)))
    return P
end

"""
    computePlmx!(P::AbstractVector, x, lmax::Integer, coeff::AbstractMatrix; [norm = SphericalHarmonics.LMnorm()])

Compute an entire set of normalized Associated Legendre Polynomials ``\\bar{P}_ℓ^m(x)``
using the given coefficients, and store in the array `P`.
The matrix `coeff` may be computed using [`compute_coefficients`](@ref).

The argument `x` needs to lie in ``-1 ≤ x ≤ 1``. The function implicitly assumes that
``x = \\cos\\theta`` where ``0 ≤ \\theta ≤ π``.

The keyword argument `norm` may be used to specify the how the polynomials are normalized.
See [`computePlmcostheta`](@ref) for the possible normalization options.

    computePlmx!(P::AbstractVector, x, lmax::Integer, m::Integer, coeff::AbstractMatrix)

Compute the set of normalized Associated Legendre Polynomials ``\\bar{P}_ℓ^m(x)`` for for the specified ``m``
and all ``ℓ`` lying in ``|m| ≤ ℓ ≤ ℓ_\\mathrm{max}`` .
"""
function computePlmx!(P::AbstractVector, x, lmax::Integer, args...; norm::PLMnorm = LMnorm())
    @assert lmax >= 0 "degree must be non-negative"
    -1 <= x <= 1 || throw(DomainError("x", "The argument to associated Legendre polynomials must satisfy -1 <= x <= 1"))
    zeroP!(P, lmax)
    _computePlmx_range!(P, x, 0:lmax, args...; norm = norm)
    return P
end

"""
    computePlmx!(S::SphericalHarmonicsCache, x[, lmax::Integer]; [norm = SphericalHarmonics.LMnorm()])

Compute an entire set of normalized Associated Legendre Polynomials ``\\bar{P}_ℓ^m(x)``
using the pre-computed coefficients in `S`, and store the result in `S`. If `lmax` is not provided,
the value of `lmax` for which coefficients have been computed in `S` is used.

The keyword argument `norm` may be used to specify the how the polynomials are normalized.
See [`computePlmcostheta`](@ref) for the possible normalization options.
"""
function computePlmx!(S::SphericalHarmonicsCache, x, lmax::Integer = _lmax(S); norm::PLMnorm = LMnorm())
    compute_coefficients!(S, lmax)
    P = computePlmx!(S.P, x, lmax, S.C, norm = norm)
    S.P = P
    return S.P
end
function computePlmx!(P::AssociatedLegendrePolynomials, x, lmax::Integer, coeff = nothing; kw...)
    @assert lmax >= 0 "degree must be non-negative"
    -1 <= x <= 1 || throw(DomainError("x", "The argument to associated Legendre polynomials must satisfy -1 <= x <= 1"))
    if P.initialized && P.cosθ == x && P.lmax >= lmax
        return P
    elseif P.initialized && P.cosθ == x && P.lmax < lmax
        _computePlmx_range!(parent(P), x, P.lmax+1:lmax, coeff; kw...)
    else
        computePlmx!(parent(P), x, lmax, coeff; kw...)
        P = @set P.cosθ = x
    end
    P = @set P.lmax = lmax
    P = @set P.initialized = true
    return P
end
function _computePlmx_range!(P::AbstractVector, x, l_range::AbstractUnitRange{<:Integer}, args...; norm::PLMnorm = LMnorm())
    lmin, lmax = map(Int, extrema(l_range))
    cosθ, sinθ  = promote(x, √(1-x^2))
    @inbounds _computePlmcostheta!(P, cosθ, sinθ, lmin:lmax, args...; norm = norm)
    return P
end

"""
    computePlmx(x; lmax::Integer, [m::Integer], [norm = SphericalHarmonics.LMnorm()])
    computePlmx(x, lmax::Integer[, m::Integer]; [norm = SphericalHarmonics.LMnorm()])

Compute an entire set of normalized Associated Legendre Polynomials ``\\bar{P}_ℓ^m(x)`` where
``0 ≤ ℓ ≤ ℓ_\\mathrm{max}`` and ``0 ≤ m ≤ ℓ``. If `m` is provided then only the polynomials for that
azimuthal order are computed.

The argument `x` needs to lie in ``-1 ≤ x ≤ 1``. The function implicitly assumes that
``x = \\cos\\theta`` where ``0 ≤ \\theta ≤ π``.

The keyword argument `norm` may be used to specify the how the polynomials are normalized.
See [`computePlmcostheta`](@ref) for the possible normalization options.
"""
computePlmx(x; lmax::Integer, m::Union{Integer, Nothing} = nothing, norm::PLMnorm = LMnorm()) = computePlmx(x, lmax, m; norm = norm)
function computePlmx(x, lmax::Integer, m::Union{Integer, Nothing} = nothing; norm::PLMnorm = LMnorm())
    P = allocate_p(float(typeof(x)), lmax)
    coeff = compute_coefficients(lmax, m)
    _applymaybeallm(computePlmx!, P, x, lmax, m, coeff; norm = norm)
    return AssociatedLegendrePolynomials(x, lmax, P, true)
end

"""
    computePlmcostheta!(P::AbstractVector, θ, lmax::Integer, coeff; [norm = SphericalHarmonics.LMnorm()])
    computePlmcostheta!(P::AbstractVector, θ::SphericalHarmonics.Pole, lmax::Integer; [norm = SphericalHarmonics.LMnorm()])

Compute an entire set of normalized Associated Legendre Polynomials ``\\bar{P}_ℓ^m(\\cos θ)``
using the given coefficients, and store in the array `P`. The matrix `coeff` may be computed
using [`compute_coefficients`](@ref).

The keyword argument `norm` may be used to specify the how the polynomials are normalized.
See [`computePlmcostheta`](@ref) for the possible normalization options.

    computePlmcostheta!(P::AbstractVector, θ, lmax::Integer, m::Integer, coeff)

Compute the Associated Legendre Polynomials ``\\bar{P}_ℓ^m(\\cos θ)`` for for the specified ``m``
and all ``ℓ`` lying in ``|m| ≤ ℓ ≤ ℓ_\\mathrm{max}``. The array `P` needs to be large enough to hold all the polynomials
for ``0 ≤ ℓ ≤ ℓ_\\mathrm{max}`` and ``0 ≤ m ≤ ℓ``, as the recursive evaluation requires the computation of extra elements.
Pre-existing values in `P` may be overwritten, even for azimuthal orders not equal to ``m``.
"""
function computePlmcostheta!(P::AbstractVector, θ, lmax, args...; norm::PLMnorm = LMnorm())
    zeroP!(P, lmax)
    _computePlmcostheta_range!(P, θ, 0:lmax, args...; norm = norm)
    return P
end
function _computePlmcostheta_range!(P::AbstractVector, θ, l_range::AbstractUnitRange{<:Integer}, args...; norm::PLMnorm = LMnorm())
    lmin, lmax = map(Int, extrema(l_range))
    cosθ, sinθ  = promote(cos(θ), sin(θ))
    @inbounds _computePlmcostheta!(P, cosθ, sinθ, lmin:lmax, args...; norm = norm)
    return P
end

function computePlmcostheta!(P::AbstractVector, θ::Pole, lmax, args...; norm::PLMnorm = LMnorm())
    zeroP!(P, lmax)
    _computePlmcostheta_range!(P, θ, 0:lmax, args...; norm = norm)
    return P
end
function _computePlmcostheta_range!(P::AbstractVector, θ::Pole, l_range::AbstractUnitRange{<:Integer}, m::Integer, args...; kw...)
    lmin, lmax = map(Int, extrema(l_range))
    if !iszero(m)
        return P
    end
    _computePlmcostheta_m0_range!(P, θ, lmin:lmax; kw...)
    return P
end
function _computePlmcostheta_range!(P::AbstractVector, θ::Pole, l_range::AbstractUnitRange{<:Integer}, args...; kw...)
    _computePlmcostheta_m0_range!(P, θ, l_range; kw...)
    return P
end

function _computePlmcostheta_m0_range!(P::AbstractVector{T}, ::NorthPole, l_range::AbstractUnitRange{<:Integer}; norm::PLMnorm = LMnorm()) where {T}
    lmin, lmax = map(Int, extrema(l_range))
    Plm = _wrapSHArray(P, lmax, ZeroTo)
    @inbounds for l in lmin:lmax
        Plm[(l, 0)] = invsqrt2π * √(T(2l + 1))
    end
    normalize!(Plm, norm, l_range, 0)
    return P
end

function _computePlmcostheta_m0_range!(P::AbstractVector{T}, ::SouthPole, l_range::AbstractUnitRange{<:Integer}; norm::PLMnorm = LMnorm()) where {T}
    lmin, lmax = map(Int, extrema(l_range))
    Plm = _wrapSHArray(P, lmax, ZeroTo)
    phase = 1
    @inbounds for l in lmin:lmax
        Plm[(l, 0)] = phase * invsqrt2π * √(T(2l + 1))
        phase *= -1
    end
    normalize!(Plm, norm, l_range, 0)
    return P
end

function computePlmcostheta!(P::AssociatedLegendrePolynomials, θ::Pole, lmax, args...; kw...)
    P = _computePlmcostheta_alp!(P, θ, lmax, args...; kw...)
    return P
end
function computePlmcostheta!(P::AssociatedLegendrePolynomials, θ, lmax, args...; kw...)
    P = _computePlmcostheta_alp!(P, θ, lmax, args...; kw...)
    return P
end
function _computePlmcostheta_alp!(P::AssociatedLegendrePolynomials, θ, lmax, args...; kw...)
    if P.initialized && P.cosθ == cos(θ) && P.lmax >= lmax
        return P
    elseif P.initialized && P.cosθ == cos(θ) && P.lmax < lmax
        _computePlmcostheta_range!(parent(P), θ, P.lmax+1:lmax, args...; kw...)
    else
        computePlmcostheta!(parent(P), θ, lmax, args...; kw...)
        P = @set P.cosθ = cos(θ)
    end
    P = @set P.lmax = lmax
    P = @set P.initialized = true
    return P
end

"""
    computePlmcostheta!(S::SphericalHarmonicsCache, θ, [lmax::Integer]; [norm = SphericalHarmonics.LMnorm()])

Compute an entire set of normalized Associated Legendre Polynomials ``\\bar{P}_ℓ^m(\\cos θ)``
using the pre-computed coefficients in `S`, and store the result in `S`. If `lmax` is not provided,
the value of `lmax` for which coefficients have been computed in `S` is used.

The keyword argument `norm` may be used to specify the how the polynomials are normalized.
See [`computePlmcostheta`](@ref) for the possible normalization options.
"""
function computePlmcostheta!(S::SphericalHarmonicsCache, θ, lmax::Integer = _lmax(S); kw...)
    compute_coefficients!(S, lmax)
    P = computePlmcostheta!(S.P, θ, lmax, S.C; kw...)
    S.P = P
    return P
end

"""
    computePlmcostheta(θ; lmax::Integer, [m::Integer], [norm = SphericalHarmonics.LMnorm()])
    computePlmcostheta(θ, lmax::Integer[, m::Integer]; [norm = SphericalHarmonics.LMnorm()])

Compute an entire set of normalized Associated Legendre Polynomials ``\\bar{P}_ℓ^m(\\cos θ)`` where
``0 ≤ ℓ ≤ ℓ_\\mathrm{max}`` and ``0 ≤ m ≤ ℓ`` for colatitude ``\\theta``. If `m` is provided then only the
polynomials corresponding to the specified `m` are computed.

The polynomials are normalized as

```math
\\bar{P}_{\\ell}^m = \\sqrt{\\frac{(2\\ell + 1)(\\ell-m)!}{2\\pi (\\ell+m)!}} P_{\\ell m},
```

where ``P_{\\ell m}`` are the standard
[associated Legendre polynomials](https://en.wikipedia.org/wiki/Associated_Legendre_polynomials#Alternative_notations),
and are defined in terms of Legendre polynomials ``P_\\ell(x)`` as

```math
P_{\\ell m}\\left(x\\right)=(-1)^m \\left(1-x^{2}\\right)^{m/2}\\frac{d^{m}}{dx^{m}}P_{\\ell}\\left(x\\right).
```

The normalized polynomials ``\\bar{P}_{\\ell}^m`` satisfy

```math
\\int_{-1}^{1} dx\\,\\left| \\bar{P}_{\\ell}^m(x) \\right|^2 = \\frac{1}{\\pi}
```

A different normalization may be chosen by specifying the keyword argument `norm`, with optionally the
Condon-Shortley phase disabled by passing the keyword argument `csphase` to the constructor of the
normalization specifier. The possible normalization options are:

* `SphericalHarmonics.LMnorm([; csphase = true])`: the default normalization described above
* `SphericalHarmonics.Orthonormal([; csphase = true])`: Orthonormal polynomials that are defined as
```math
\\tilde{P}_{\\ell}^m = \\sqrt{\\frac{(2\\ell + 1)(\\ell-m)!}{2(\\ell+m)!}} P_{\\ell m} =
\\sqrt{\\pi} \\bar{P}_{\\ell m},
```
and satisfy
```math
\\int_{-1}^{1} \\tilde{P}_{\\ell m}(x) \\tilde{P}_{k m}(x) dx = \\delta_{ℓk}
```
* `SphericalHarmonics.Unnormalized([; csphase = true])`: The polynomials ``P_{ℓm}`` that satisfy ``P_{ℓm}(1)=\\delta_{m0}``
within numerical precision bounds, as well as

```math
\\int_{-1}^{1} P_{\\ell m}(x) P_{k m}(x) dx = \\frac{2(\\ell+m)!}{(2\\ell+1)(\\ell-m)!}\\delta_{ℓk}
```

In each case setting `csphase = false` will lead to an overall factor of ``(-1)^m`` being multiplied
to the polynomials.

Returns an `AbstractVector` that may be indexed using `(ℓ,m)` pairs aside from the
canonical indexing with `Int`s.

The precision of the result may be increased by using arbitrary-precision arguments.

# Examples
```jldoctest
julia> P = computePlmcostheta(pi/2, lmax = 1)
3-element normalized AssociatedLegendrePolynomials{Float64} for lmax = 1 and cosθ = 6.123e-17:
  0.3989422804014327
  4.231083042742082e-17
 -0.4886025119029199

julia> P[(0,0)]
0.3989422804014327

julia> P = computePlmcostheta(big(pi)/2, lmax = 1)
3-element normalized AssociatedLegendrePolynomials{BigFloat} for lmax = 1 and cosθ = 5.485e-78:
  0.3989422804014326779399460599343818684758586311649346576659258296706579258993008
  3.789785583114350800838137317730900078444216599640987847808409161681770236721676e-78
 -0.4886025119029199215863846228383470045758856081942277021382431574458410003616367
```
"""
computePlmcostheta(θ; lmax::Integer, m::Union{Integer, Nothing} = nothing, norm::PLMnorm = LMnorm()) = computePlmcostheta(θ, lmax, m; norm = norm)

_applymaybeallm(f, P, θ, lmax, m::Nothing, coeff...; kw...) = f(P, θ, lmax, coeff...; kw...)
_applymaybeallm(f, P, θ, lmax, m, coeff...; kw...) = f(P, θ, lmax, m, coeff...; kw...)

_maybecomputecoeff(lmax, θ::Pole, m) = nothing
_maybecomputecoeff(lmax, θ, m) = compute_coefficients(promote_type(Float64, float(typeof(θ))), Int(lmax), m)

function computePlmcostheta(θ, lmax::Integer, m::Union{Integer, Nothing} = nothing; norm::PLMnorm = LMnorm())
    coeff = _maybecomputecoeff(lmax, θ, m)
    P = allocate_p(float(typeof(θ)), lmax)
    _applymaybeallm(computePlmcostheta!, P, θ, lmax, m, coeff; norm = norm)
    return AssociatedLegendrePolynomials(cos(θ), lmax, P, true)
end

@doc raw"""
    SphericalHarmonics.associatedLegendrex(x; l::Integer, m::Integer, [coeff = nothing], [norm = SphericalHarmonics.LMnorm()])

Evaluate the normalized associated Legendre polynomial ``\bar{P}_ℓ^m(x)``.
Optionally a matrix of coefficients returned by [`compute_coefficients`](@ref) may be provided.

The keyword argument `norm` may be used to specify the how the polynomials are normalized.
See [`computePlmcostheta`](@ref) for the possible normalization options.
"""
associatedLegendrex(x; l::Integer, m::Integer, coeff = nothing, norm::PLMnorm = LMnorm()) =
    associatedLegendrex(x, l, m, coeff; norm = norm)
function associatedLegendrex(x, l::Integer, m::Integer, coeff = nothing; norm::PLMnorm = LMnorm())
    _computePlmcostheta(x, √(1-x^2), l, m, coeff; norm = norm)
end

@doc raw"""
    SphericalHarmonics.associatedLegendre(θ; l::Integer, m::Integer, [coeff = nothing], [norm = SphericalHarmonics.LMnorm()])

Evaluate the normalized associated Legendre polynomial ``\bar{P}_ℓ^m(\cos \theta)``.
Optionally a matrix of coefficients returned by [`compute_coefficients`](@ref) may be provided.

The keyword argument `norm` may be used to specify the how the polynomials are normalized.
See [`computePlmcostheta`](@ref) for the possible normalization options.
"""
associatedLegendre(θ; l::Integer, m::Integer, coeff = nothing, norm::PLMnorm = LMnorm()) = associatedLegendre(θ, l, m, coeff; norm = norm)
function associatedLegendre(θ, l::Integer, m::Integer, coeff = nothing; norm::PLMnorm = LMnorm())
    _computePlmcostheta(cos(θ), sin(θ), l, m, coeff; norm = norm)
end

_Plm_pole(T, l, m, norm) = invsqrt2π * √(T(2l + 1)) * invnorm(T(l), T(m), norm)

function associatedLegendre(::NorthPole, l::Integer, m::Integer, coeff = nothing; norm::PLMnorm = LMnorm())
    T = promote_type(Float64, float(typeof(l)))
    if m != 0
        return zero(T)
    end
    _Plm_pole(T, l, m, norm)
end

function associatedLegendre(::SouthPole, l::Integer, m::Integer, coeff = nothing; norm::PLMnorm = LMnorm())
    T = promote_type(Float64, float(typeof(l)))
    if m != 0
        return zero(T)
    end
    (-1)^Int(l) * _Plm_pole(T, l, m, norm)
end

######################################################################################
# Spherical Harmonics
######################################################################################

function phase(::RealHarmonics, ::Type{FullRange}, m, ϕ, norm, CSphase)
    S, C = sincos(abs(Int(m))*ϕ)
    C, S
end

function phase(::ComplexHarmonics, ::Type{FullRange}, m, ϕ, norm, CSphase)
    ep = cis(Int(m)*ϕ) * norm
    em = CSphase * conj(ep)
    ep, em
end

phase(::RealHarmonics, ::Type{ZeroTo}, m, ϕ, norm, CSphase) = cos(m*ϕ)
phase(::ComplexHarmonics, ::Type{ZeroTo}, m, ϕ, norm, CSphase) = cis(m*ϕ) * norm

@propagate_inbounds function fill_m_maybenegm!(Y, P, lmax, m, ϕ, CSphase, ::Type{FullRange}, SHType)
    m >= 0 || throw(ArgumentError("m must be ≥ 0"))
    @assert lmax >= 0 "degree must be non-negative"
    phasempos, phasemneg = phase(SHType, FullRange, m, ϕ, invsqrt2, CSphase)

    PS = _wrapSHArray(P, lmax, ZeroTo)
    Ylm = _wrapSHArray(Y, lmax, FullRange)

    for l in m:Int(lmax)
        Plm = PS[(l, m)]
        Ylm[(l, -m)] = phasemneg * Plm
        Ylm[(l,  m)] = phasempos * Plm
    end
    return Y
end

@propagate_inbounds function fill_m_maybenegm!(Y, P, lmax, m, ϕ, CSphase, ::Type{ZeroTo}, SHType)
    m >= 0 || throw(ArgumentError("m must be ≥ 0"))
    @assert lmax >= 0 "degree must be non-negative"
    phasem = phase(SHType, ZeroTo, m, ϕ, invsqrt2, CSphase)

    PS = _wrapSHArray(P, lmax, ZeroTo)
    Ylm = _wrapSHArray(Y, lmax, ZeroTo)

    for l in abs(Int(m)):Int(lmax)
        Ylm[(l,m)] = phasem * PS[(l,m)]
    end
    return Y
end

@propagate_inbounds function fill_m!(Y, P, lmax, m, ϕ, CSphase, ::Type{FullRange}, SHType::ComplexHarmonics)
    @assert lmax >= 0 "degree must be non-negative"
    phasempos, phasemneg = phase(SHType, FullRange, m, ϕ, invsqrt2, CSphase)

    PS = _wrapSHArray(P, lmax, ZeroTo)
    Ylm = _wrapSHArray(Y, lmax, FullRange)

    for l in abs(Int(m)):Int(lmax)
        Plm = PS[(l, abs(m))]
        if m >= 0
            Ylm[(l, m)] = phasempos * Plm
        else
            Ylm[(l, m)] = (-1)^m * phasempos * Plm
        end
    end
    return Y
end

@propagate_inbounds function fill_m!(Y, P, lmax, m, ϕ, CSphase, ::Type{FullRange}, SHType::RealHarmonics)
    @assert lmax >= 0 "degree must be non-negative"
    C, S = phase(SHType, FullRange, m, ϕ, invsqrt2, CSphase)

    PS = _wrapSHArray(P, lmax, ZeroTo)
    Ylm = _wrapSHArray(Y, lmax, FullRange)

    for l in abs(Int(m)):Int(lmax)
        Plm = PS[(l, abs(m))]
        if m >= 0
            Ylm[(l, m)] = C * Plm
        else
            Ylm[(l, m)] = S * Plm
        end
    end
    return Y
end

@propagate_inbounds function fill_m!(Y, P, lmax, m, ϕ, CSphase, ::Type{ZeroTo}, SHType)
    m >= 0 || throw(ArgumentError("m must be ≥ 0"))
    @assert lmax >= 0 "degree must be non-negative"

    phasem = phase(SHType, ZeroTo, m, ϕ, invsqrt2, CSphase)
    PS = _wrapSHArray(P, lmax, ZeroTo)
    Ylm = _wrapSHArray(Y, lmax, ZeroTo)

    for l in Int(m):Int(lmax)
        Ylm[(l, m)] = phasem * P[(l, m)]
    end
    return Y
end

@propagate_inbounds function computeYlm_maybeallm!(Y, P, ϕ, lmax, ::Nothing, m_range, SHType)
    @assert lmax >= 0 "degree must be non-negative"
    PS = _wrapSHArray(P, lmax, ZeroTo)
    Ylm = _wrapSHArray(Y, lmax, m_range)
    for l in 0:Int(lmax)
        Ylm[(l, 0)] = PS[(l, 0)] * invsqrt2
    end

    CSphase = 1
    for m in 1:Int(lmax)
        CSphase *= -1
        fill_m_maybenegm!(Y, P, lmax, m, ϕ, CSphase, m_range, SHType)
    end
    return Y
end
@propagate_inbounds function computeYlm_maybeallm!(Y, P, ϕ, lmax, m::Integer, m_range, SHType)
    @assert lmax >= 0 "degree must be non-negative"
    -lmax <= m <= lmax || throw(ArgumentError("m = $m does not satisfy 0 ≤ |m| ≤ lmax = $lmax"))

    PS = _wrapSHArray(P, lmax, ZeroTo)
    Ylm = _wrapSHArray(Y, lmax, m_range)
    if iszero(m)
        for l in 0:Int(lmax)
            Ylm[(l, 0)] = PS[(l, 0)] * invsqrt2
        end
    else
        fill_m!(Y, P, lmax, m, ϕ, (-1)^m, m_range, SHType)
    end
    return Y
end

_maybeabs(::Nothing) = nothing
_maybeabs(m::Integer) = abs(m)

function checksize_zeroY!(Y, P, lmax, m_range)
    checksize(length(P), sizeP(Int(lmax)))
    Ny = sizeY(Int(lmax), m_range)
    checksize(length(Y), Ny)
    Yv = @view Y[firstindex(Y) .+ (0:Ny-1)]
    fill!(Yv, zero(eltype(Yv)))
    return nothing
end

"""
    computeYlm!(Y::AbstractVector, P::AbstractVector, θ, ϕ; lmax::Integer, [m::Integer], [m_range = SphericalHarmonics.FullRange], [SHType = SphericalHarmonics.ComplexHarmonics()])
    computeYlm!(Y::AbstractVector, P::AbstractVector, θ, ϕ, lmax::Integer, [m::Integer, [m_range = SphericalHarmonics.FullRange, [SHType = SphericalHarmonics.ComplexHarmonics()]]])

Compute an entire set of spherical harmonics ``Y_{ℓ,m}(θ,ϕ)``
using the precomputed associated Legendre Polynomials ``\\bar{P}_ℓ^m(\\cos θ)``,
and store in the array `Y`. The array `P` may be computed using [`computePlmcostheta`](@ref).

The optional argument `m_range` decides if the spherical harmonics for negative `m` values are computed.
By default the functions for all values of `m` are evaluated. Setting `m_range` to `SphericalHarmonics.ZeroTo` would result
in only functions for `m ≥ 0` being evaluated. Providing `m` would override this, in which case only the polynomials
corresponding to the azimuthal order `m` would be computed.

The optional argument `SHType` may be used to choose between real and complex harmonics.
To compute real spherical harmonics, set this to `RealHarmonics()`.

!!! note
    This function assumes that the associated Legendre Polynomials have been pre-computed, and does not perform any
    check on the values of `P`.
"""
function computeYlm!(Y::AbstractVector, P::AbstractVector, θ,
    ϕ, lmax::Integer, m::Union{Integer,Nothing} = nothing,
    m_range::SHMRange = FullRange, SHType::HarmonicType = ComplexHarmonics())

    checksize_zeroY!(Y, P, lmax, m_range)
    @inbounds computeYlm_maybeallm!(Y, P, ϕ, lmax, m, m_range, SHType)
    return Y
end
function computeYlm!(Y::AbstractVector, P::AbstractVector, θ,
    ϕ, lmax::Integer, m_range::SHMRange, SHType::HarmonicType = ComplexHarmonics())

    computeYlm!(Y, P, θ, ϕ, lmax, nothing, m_range, SHType)
    return Y
end

function computeYlm!(Y::AbstractVector, P::AbstractVector, θ,
    ϕ; lmax::Integer, m::Union{Integer,Nothing} = nothing,
    m_range::SHMRange = FullRange, SHType::HarmonicType = ComplexHarmonics())

    computeYlm!(Y, P, θ, ϕ, lmax, m, m_range, SHType)
    return Y
end

function _computeYlm_poles!(Y, P, ϕ, lmax, m_range)
    PS = _wrapSHArray(P, lmax, ZeroTo)
    Ylm = _wrapSHArray(Y, lmax, m_range)

    @inbounds for l in 0:lmax
        Ylm[(l, 0)] = PS[(l, 0)] * invsqrt2
    end
    return Y
end

function computeYlm!(Y::AbstractVector, P::AbstractVector, θ::Pole,
    ϕ, lmax::Integer, m_range::SHMRange = FullRange,
    SHType::HarmonicType = ComplexHarmonics())

    checksize_zeroY!(Y, P, lmax, m_range)
    _computeYlm_poles!(Y, P, ϕ, lmax, m_range)
    return Y
end
function computeYlm!(Y::AbstractVector, P::AbstractVector, θ::Pole,
    ϕ, lmax::Integer, m::Integer, m_range::SHMRange = FullRange,
    SHType::HarmonicType = ComplexHarmonics())

    checksize_zeroY!(Y, P, lmax, m_range)
    if iszero(m)
        _computeYlm_poles!(Y, P, ϕ, lmax, m_range)
    end
    return Y
end

"""
    computeYlm!(S::SphericalHarmonicsCache, θ, ϕ, [lmax::Integer])

Compute an entire set of spherical harmonics ``Y_{ℓ,m}(θ,ϕ)`` for ``0 ≤ ℓ ≤ ℓ_\\mathrm{max}`` using the
pre-computed associated Legendre polynomials saved in `S`, and store the result in `S`. If `lmax` is not provided,
the value of `lmax` for which associated Legendre polynomials have been computed in `S` is used.

!!! note
    This function assumes that the associated Legendre polynomials have been pre-computed, and does not perform
    any check on their values. In general `computeYlm!(S::SphericalHarmonicsCache, θ, ϕ, lmax)` should only be
    called after a preceeding call to `computePlmcostheta!(S, θ, lmax)` in order to obtain meaningful results.
"""
function computeYlm!(S::SphericalHarmonicsCache, θ, ϕ, lmax::Integer = _lmax(S))
    @assert lmax <= _lmax(S) "Plm for lmax = $lmax is not available, please run computePlmcostheta!(S, θ, lmax) first"
    !S.P.initialized && throw(ArgumentError("please run computePlmcostheta!(S, θ, lmax) first"))
    computeYlm!(S.Y, S.P, θ, ϕ, lmax, nothing, _mrange_basetype(S), S.SHType)
    return S.Y
end

"""
    computeYlm!(Y::AbstractVector, θ, ϕ; lmax::Integer, [m::Integer] [m_range = SphericalHarmonics.FullRange], [SHType = SphericalHarmonics.ComplexHarmonics()])

Compute an entire set of spherical harmonics ``Y_{ℓ,m}(θ,ϕ)`` for ``0 ≤ ℓ ≤ ℓ_\\mathrm{max}``,
and store them in the array `Y`.

The optional argument `m_range` decides if the spherical harmonics for negative `m` values are computed.
By default the functions for all values of `m` are evaluated. Setting `m_range` to `SphericalHarmonics.ZeroTo` would result
in only functions for `m ≥ 0` being evaluated. Providing `m` would override this, in which case only the polynomials
corresponding to the azimuthal order `m` would be computed.

The optional argument `SHType` may be used to choose between real and complex harmonics.
To compute real spherical harmonics, set this to `SphericalHarmonics.RealHarmonics()`.
"""
function computeYlm!(Y::AbstractVector, θ::Number, ϕ, lmax::Integer, m::Union{Integer,Nothing} = nothing,
    m_range::SHMRange = FullRange, SHType::HarmonicType = ComplexHarmonics())
    P = computePlmcostheta(θ, lmax, _maybeabs(m))
    computeYlm!(Y, P, θ, ϕ, lmax, m, m_range, SHType)
    return Y
end
function computeYlm!(Y::AbstractVector, θ::Number, ϕ, lmax::Integer, m_range::SHMRange,
    SHType::HarmonicType = ComplexHarmonics())
    computeYlm!(Y, θ, ϕ, lmax, nothing, m_range, SHType)
    return Y
end

function computeYlm!(Y::AbstractVector, θ::Number, ϕ; lmax::Integer, m::Union{Integer,Nothing} = nothing,
    m_range::SHMRange = FullRange, SHType::HarmonicType = ComplexHarmonics())
    computeYlm!(Y, θ, ϕ, lmax, m, m_range, SHType)
    return Y
end

function _computeYlm(P, θ, ϕ, lmax, m, m_range, SHType)
    T = eltypeY(promote_type(eltype(P), typeof(ϕ)), SHType)
    Y = allocate_y(T, lmax, m_range)
    computeYlm!(Y, P, θ, ϕ, lmax, m, m_range, SHType)
    return Y
end

@doc raw"""
    computeYlm(θ, ϕ; lmax::Integer, [m::Integer], [m_range = SphericalHarmonics.FullRange], [SHType = SphericalHarmonics.ComplexHarmonics()])
    computeYlm(θ::SphericalHarmonics.Pole; lmax::Integer, [m_range = SphericalHarmonics.FullRange], [SHType = SphericalHarmonics.ComplexHarmonics()])

Compute an entire set of spherical harmonics ``Y_{ℓ,m}(θ,ϕ)`` for
``0 ≤ ℓ ≤ ℓ_\mathrm{max}`` and ``-ℓ ≤ m ≤ ℓ``, for colatitude ``\theta`` and
azimuth ``\phi``. If ``m`` is provided then only the polynomials for the specified ``m`` are computed.

Returns an `AbstractVector` that may be indexed using `(l,m)` pairs aside from the canonical
indexing with `Int`s.

The optional argument `m_range` decides if the spherical harmonics for negative `m` values are computed.
By default the functions for all values of `m` are evaluated. Setting `m_range` to `SphericalHarmonics.ZeroTo` would result
in only functions for `m ≥ 0` being evaluated.

The optional argument `SHType` may be used to choose between real and complex harmonics.
To compute real spherical harmonics, set this to `SphericalHarmonics.RealHarmonics()`. The complex harmonics
are defined as

```math
Y_{\ell,m}\left(\theta,\phi\right)=\begin{cases}
\frac{1}{\sqrt{2}}\bar{P}_{\ell}^{m}\left(\cos\theta\right)\exp\left(im\phi\right), & m\geq0,\\
\left(-1\right)^{m}Y_{\ell,-m}^{*}\left(\theta,\phi\right), & m<0.
\end{cases}
```

This definition corresponds to [Laplace sphercial harmonics](https://en.wikipedia.org/wiki/Spherical_harmonics#Conventions),
whereas the quantum mechanical definition prepends a Condon-Shortley phase on the harmonics.

The real spherical harmonics are defined as

```math
Y_{\ell,m}\left(\theta,\phi\right)=\begin{cases}
\bar{P}_{\ell}^{\left|m\right|}\left(\cos\theta\right)\sin\left|m\right|\phi, & m<0,\\
\bar{P}_{\ell}^{0}\left(\cos\theta\right)/\sqrt{2}, & m=0,\\
\bar{P}_{\ell}^{m}\left(\cos\theta\right)\cos m\phi, & m>0.
\end{cases}
```

The precision of the result may be increased by using arbitrary-precision arguments.
"""*
"""
# Examples
```jldoctest
julia> Y = computeYlm(pi/2, 0, lmax = 1)
4-element SHArray(::$(Array{Complex{Float64},1}), (ML(0:1, -1:1),)):
     0.2820947917738782 + 0.0im
     0.3454941494713355 - 0.0im
 2.9918275112863375e-17 + 0.0im
    -0.3454941494713355 - 0.0im

julia> Y[(1,-1)] # index using (l,m)
0.3454941494713355 - 0.0im

julia> Y = computeYlm(big(pi)/2, big(0), lmax = big(1)) # Arbitrary precision
4-element SHArray(::$(Array{Complex{BigFloat},1}), (ML(0:1, -1:1),)):
    0.2820947917738781434740397257803862929220253146644994284220428608553212342207478 + 0.0im
    0.3454941494713354792652446460318896831393773703262433134867073548945156550201567 - 0.0im
 2.679783085063171668225419916118067917387251852939708540164955895366691604430101e-78 + 0.0im
   -0.3454941494713354792652446460318896831393773703262433134867073548945156550201567 - 0.0im

julia> computeYlm(SphericalHarmonics.NorthPole(), 0, lmax = 1)
4-element SHArray(::$(Array{Complex{Float64},1}), (ML(0:1, -1:1),)):
 0.2820947917738782 + 0.0im
               -0.0 + 0.0im
   0.48860251190292 + 0.0im
                0.0 + 0.0im

julia> Y = computeYlm(pi/2, pi/3, lmax = 1, m_range = SphericalHarmonics.ZeroTo, SHType = SphericalHarmonics.RealHarmonics())
3-element SHArray(::$(Array{Float64,1}), (ML(0:1, 0:1),)):
  0.2820947917738782
  2.9918275112863375e-17
 -0.24430125595146002
```
"""
function computeYlm(θ, ϕ, lmax::Integer, m::Union{Integer,Nothing} = nothing,
    m_range::SHMRange = FullRange, SHType::HarmonicType = ComplexHarmonics())

    T = float(promote_type(typeof(θ), typeof(ϕ)))
    P = computePlmcostheta(T(θ), lmax, _maybeabs(m))
    Y = _computeYlm(P, T(θ), T(ϕ), lmax, m, m_range, SHType)
    return Y
end
function computeYlm(θ, ϕ, lmax::Integer,
    m_range::SHMRange, SHType::HarmonicType = ComplexHarmonics())

    computeYlm(θ, ϕ, lmax, nothing, m_range, SHType)
end

function computeYlm(θ::Pole, ϕ, lmax::Integer, m::Union{Integer,Nothing} = nothing,
    m_range::SHMRange = FullRange,
    SHType::HarmonicType = ComplexHarmonics())

    P = computePlmcostheta(θ, lmax, _maybeabs(m))
    Y = _computeYlm(P, θ, ϕ, lmax, m, m_range, SHType)
    return Y
end
function computeYlm(θ::Pole, ϕ, lmax::Integer, m_range::SHMRange,
    SHType::HarmonicType = ComplexHarmonics())

    computeYlm(θ, ϕ, lmax, nothing, m_range, SHType)
end

function computeYlm(θ::Pole, lmax::Integer, m_range::SHMRange = FullRange,
    SHType::HarmonicType = ComplexHarmonics())
    computeYlm(θ, 0, lmax, nothing, m_range, SHType)
end

function computeYlm(θ, ϕ; lmax::Integer, m::Union{Integer,Nothing} = nothing,
    m_range::SHMRange = FullRange, SHType::HarmonicType = ComplexHarmonics())

    computeYlm(θ, ϕ, lmax, m, m_range, SHType)
end
function computeYlm(θ::Pole; lmax::Integer, m::Union{Integer,Nothing} = nothing,
    m_range::SHMRange = FullRange, SHType::HarmonicType = ComplexHarmonics())
    computeYlm(θ, 0, lmax, m, m_range, SHType)
end

"""
    SphericalHarmonics.sphericalharmonic(θ, ϕ; l, m, [SHType = ComplexHarmonics()], [coeff = nothing])

Evaluate the spherical harmonic ``Y_{ℓ,m}(θ, ϕ)``. The flag `SHType` sets the type of the harmonic computed,
and setting this to `RealHarmonics()` would evaluate real spherical harmonics. Optionally a precomputed
matrix of coefficients returned by [`compute_coefficients`](@ref) may be provided.

# Example
```jldoctest
julia> SphericalHarmonics.sphericalharmonic(π/3, π/3, l = 500, m = 250)
-0.18910100312194328 - 0.32753254516944075im

julia> SphericalHarmonics.sphericalharmonic(π/3, π/3, l = 500, m = 250, SHType = SphericalHarmonics.RealHarmonics())
-0.26742920327340913
```
"""
function sphericalharmonic(θ, ϕ; l::Integer, m::Integer,
    SHType::HarmonicType = ComplexHarmonics(), coeff = nothing)

    sphericalharmonic(θ, ϕ, l, m, SHType, coeff)
end
function sphericalharmonic(θ, ϕ, l::Integer, m::Integer, SHType::HarmonicType = ComplexHarmonics(),
    coeff = nothing)

    P = associatedLegendre(θ, l, abs(m), coeff)
    if m == 0
        return invsqrt2 * P
    end
    phasepos, phaseneg = phase(SHType, FullRange, abs(m), ϕ, invsqrt2, (-1)^Int(m))
    norm = m >= 0 ? phasepos : phaseneg
    norm * P
end

include("precompile.jl")

end


text:
using Test
using LinearAlgebra
using Distributed
include(joinpath(Sys.BINDIR, "..", "share", "julia", "test", "testenv.jl"))
addprocs_with_testenv(3)
@test nworkers() == 3

@everywhere using Logging
for w in workers()
    # Do not log on worker nodes
    remotecall(()->global_logger(NullLogger()),w)
end

@everywhere using StochasticPrograms
using LShapedSolvers
using JuMP
using GLPKMathProgInterface
using Ipopt

τ = 1e-5
reference_solver = GLPKSolverLP()
qp_solver = IpoptSolver(print_level=0)
dlsolvers = [(LShapedSolver(reference_solver, distributed=true, log=false),"L-Shaped"),
             (LShapedSolver(qp_solver, subsolver=reference_solver, crash=Crash.EVP(), regularization = :rd, distributed = true, autotune=true, log=false),"RD L-Shaped"),
             (LShapedSolver(reference_solver, crash=Crash.EVP(), regularization = :rd, distributed = true, autotune=true, log=false, linearize=true),"Linearized RD L-Shaped"),
             (LShapedSolver(reference_solver, crash=Crash.EVP(), regularization = :tr, distributed = true, autotune=true, log=false),"TR L-Shaped"),
             (LShapedSolver(reference_solver, projectionsolver=qp_solver, regularization = :lv, distributed = true, log=false),"Leveled L-Shaped"),
             (LShapedSolver(reference_solver, regularization = :lv, distributed = true, log=false, linearize=true),"Linearized Leveled L-Shaped")]
lsolvers = [(LShapedSolver(reference_solver,log=false),"L-Shaped"),
            (LShapedSolver(qp_solver, subsolver=reference_solver, crash=Crash.EVP(), regularization = :rd, autotune=true, log=false),"RD L-Shaped"),
            (LShapedSolver(reference_solver, crash=Crash.EVP(), regularization = :rd, autotune=true, log=false, linearize=true),"Linearized RD L-Shaped"),
            (LShapedSolver(reference_solver, crash=Crash.EVP(), regularization = :tr, autotune=true, log=false),"TR L-Shaped"),
            (LShapedSolver(reference_solver, projectionsolver=qp_solver, regularization = :lv, log=false),"Leveled L-Shaped"),
            (LShapedSolver(reference_solver, regularization = :lv, log=false, linearize=true),"Linearized Leveled L-Shaped")]
problems = Vector{Tuple{<:StochasticProgram,String}}()
@info "Loading test problems..."
@info "Loading simple..."
include("simple.jl")
@info "Loading farmer..."
include("farmer.jl")
@info "Loading infeasible..."
include("infeasible.jl")

@testset "Distributed solvers" begin
    @testset "Simple problems" begin
        @info "Test problems loaded. Starting test sequence."
        @testset "Distributed $lsname Solver with Distributed Data: $name" for (lsolver,lsname) in dlsolvers, (sp,name) in problems
            optimize!(sp, solver=reference_solver)
            x̄ = optimal_decision(sp)
            Q̄ = optimal_value(sp)
            optimize!(sp, solver=lsolver)
            @test abs(optimal_value(sp) - Q̄)/(1e-10+abs(Q̄)) <= τ
            @test norm(optimal_decision(sp) - x̄)/(1e-10+norm(x̄)) <= sqrt(τ)
        end
        @testset "Distributed Bundled $lsname Solver: $name" for (lsolver,lsname) in dlsolvers, (sp,name) in problems
            optimize!(sp, solver=reference_solver)
            x̄ = optimal_decision(sp)
            Q̄ = optimal_value(sp)
            add_params!(lsolver, bundle=2)
            optimize!(sp, solver=lsolver)
            @test abs(optimal_value(sp) - Q̄)/(1e-10+abs(Q̄)) <= τ
            @test norm(optimal_decision(sp) - x̄)/(1e-10+norm(x̄)) <= sqrt(τ)
        end
        @testset "Distributed Bundled Single Node $lsname Solver: $name" for (lsolver,lsname) in dlsolvers, (sp,name) in problems
            sp_onenode = copy(sp)
            add_scenarios!(sp_onenode, scenarios(sp), workers()[1])
            optimize!(sp_onenode, solver=reference_solver)
            x̄ = optimal_decision(sp_onenode)
            Q̄ = optimal_value(sp_onenode)
            add_params!(lsolver, bundle=2)
            optimize!(sp, solver=lsolver)
            @test abs(optimal_value(sp_onenode) - Q̄)/(1e-10+abs(Q̄)) <= τ
            @test norm(optimal_decision(sp_onenode) - x̄)/(1e-10+norm(x̄)) <= sqrt(τ)
        end
        @testset "Distributed $lsname Solver: $name" for (lsolver,lsname) in dlsolvers, (sp,name) in problems
            sp_nondist = copy(sp, procs = [1])
            add_scenarios!(sp_nondist, scenarios(sp))
            optimize!(sp_nondist, solver=reference_solver)
            x̄ = optimal_decision(sp_nondist)
            Q̄ = optimal_value(sp_nondist)
            optimize!(sp_nondist, solver=lsolver)
            @test abs(optimal_value(sp_nondist) - Q̄)/(1e-10+abs(Q̄)) <= τ
            @test norm(optimal_decision(sp_nondist) - x̄)/(1e-10+norm(x̄)) <= sqrt(τ)
        end
        @testset "$lsname Solver with Distributed Data: $name" for (lsolver,lsname) in lsolvers, (sp,name) in problems
            optimize!(sp, solver=reference_solver)
            x̄ = optimal_decision(sp)
            Q̄ = optimal_value(sp)
            with_logger(NullLogger()) do
                optimize!(sp, solver=lsolver)
            end
            @test abs(optimal_value(sp) - Q̄)/(1e-10+abs(Q̄)) <= τ
            @test norm(optimal_decision(sp) - x̄)/(1e-10+norm(x̄)) <= sqrt(τ)
        end
    end
    @testset "Infeasible problem" begin
        @testset "$lsname Solver: Feasibility cuts" for (lsolver,lsname) in dlsolvers
            optimize!(infeasible, solver=reference_solver)
            x̄ = optimal_decision(infeasible)
            Q̄ = optimal_value(infeasible)
            with_logger(NullLogger()) do
                @test optimize!(infeasible, solver=lsolver) == :Infeasible
            end
            add_params!(lsolver, checkfeas=true)
            optimize!(infeasible, solver=lsolver)
            @test abs(optimal_value(infeasible) - Q̄)/(1e-10+abs(Q̄)) <= τ
            @test norm(optimal_decision(infeasible) - x̄)/(1e-10+norm(x̄)) <= sqrt(τ)
        end
        @testset "Bundled $lsname Solver: Feasibility cuts" for (lsolver,lsname) in dlsolvers
            optimize!(infeasible, solver=reference_solver)
            x̄ = optimal_decision(infeasible)
            Q̄ = optimal_value(infeasible)
            add_params!(lsolver, checkfeas=false,bundle=2)
            with_logger(NullLogger()) do
                @test optimize!(infeasible, solver=lsolver) == :Infeasible
            end
            add_params!(lsolver, checkfeas=true, bundle=2)
            optimize!(infeasible, solver=lsolver)
            @test abs(optimal_value(infeasible) - Q̄)/(1e-10+abs(Q̄)) <= τ
            @test norm(optimal_decision(infeasible) - x̄)/(1e-10+norm(x̄)) <= sqrt(τ)
        end
    end
end


text:
module Ping

using Discord

function on_message_create(client::Client, event::MessageCreate)
    if event.message.content == "ping"
       reply(client, event.message, "Pong!")
    elseif event.message.content == "pong"
       reply(client, event.message, "Ping!")
    end
end

function main()
    client = Client(ENV["DISCORD_TOKEN"])
    open(client)
    add_handler!(client, MessageCreate, on_message_create; tag=:ping)
    return client
end

end

if abspath(PROGRAM_FILE) == @__FILE__
    client = Ping.main()
    wait(client)
end


text:
# TODO: fix width limiting
# TODO: introducing vertical scrolling so nothing breaks when the current Tree is
#       higher than the terminal has lines
using REPL

include("utils.jl")

export @ishow

mutable struct Tree{F}
    head
    children_gen::F
    children::Vector

    materialized::Bool
    expanded::Bool
    options::Vector{String}

    pagesize::Int
    pageoffset::Int

    selected
    lastHeight::Int

    cursor
end

function Tree(head, children::Vector{<:Any})
    Tree(head, () -> (), children, true, false, String[], length(children), 0, nothing, 0, 0)
end

function Tree(head, f::Function)
    Tree(head, f, Any[], false, false, String[], length(f()), 0, nothing, 0, 0)
end

function children(t::Tree)
    if t.materialized
        t.children
    else
        tmp = t.children_gen()
        t.children = tmp
        t.materialized = true
        tmp
    end
end

toggle(t::Tree) = (t.expanded = !t.expanded)

showmethod(T) = which(show, (IO, T))

getfield′(x, f) = isdefined(x, f) ? getfield(x, f) : Text("#undef")

function defaultrepr(x; smethod=false)
    if smethod && showmethod(typeof(x)) ≠ showmethod(Any)
        b = IOBuffer()
        print(b, Text(io -> show(IOContext(io, :limit => true), MIME"text/plain"(), x)))
        Text(String(take!(b)))
    else
        fields = fieldnames(typeof(x))

        if isempty(fields)
            Tree(string(typeof(x), "()"), [])
        else
            Tree(string(typeof(x)),
                 [Tree(string(f), [defaultrepr(getfield′(x, f), smethod=true)]) for f in fields])
        end
    end
end

function generateTreeView(x)
    header = Text(sprint(treelabel, x, MIME"text/plain"()))

    numberofnodes(x) == 0 && return Tree(header, [])

    genchildren = function ()
      children = Any[]
      for i in 1:numberofnodes(x)
        node = treenode(x, i)

        cheader = sprint(nodelabel, x, i, MIME"text/plain"())
        if isempty(cheader)
          push!(children, hastreeview(node) ? generateTreeView(node) : node)
        elseif node === nothing
          push!(children, Text(cheader))
        else
          # would be nicer if this were a SubTree, but works fine as is
          push!(children, Tree(Text(cheader),
                               () -> [(hastreeview(node) ?
                                       generateTreeView(node) :
                                       Text(io -> show(IOContext(io, :limit => true), MIME"text/plain"(), node)))]))
      end
      end
      children
    end

    return Tree(header, genchildren)
end

macro ishow(x)
    :(ishow($(esc(x))))
end

function ishow(x)
    request(hastreeview(x) ? generateTreeView(x) : defaultrepr(x))
end

const INDENTSIZE = 3

indent(level) = " "^INDENTSIZE

function printTreeChild(buf::IOBuffer, child::Tree, cursor, term_width::Int; level::Int = 0)
    cur = cursor == -1
    symbol = length(children(child)) > 0 ? child.expanded ? "▼" : "▶" : " "

    cur ? print(buf, "[$symbol] ") : print(buf, " $symbol  ")
    if child.expanded
        # print Tree with additional nesting, but without an active cursor
        # init=true assures that the Tree printing doesn't mess with anything
        cursor = printMenu(buf, child, cursor; init=true, level = level)
    else
        # only print header
        tb = IOBuffer()
        print(tb, child.head)
        s = String(take!(tb))
        print(buf, join(limitLineLength([s], term_width-(2*level + 10)), '\n'))
    end

    cursor
end

function limitLineLength(strs, term_width)
    outstrs = String[]
    for str in strs
        if length(str) >= term_width
            while length(str) >= term_width
                push!(outstrs, str[1:term_width])
                str = str[max(term_width, 1):end]
            end
        else
            push!(outstrs, str)
        end
    end
    outstrs
end

function writeChild(buf::IOBuffer, t::Tree, idx::Int, cursor, term_width::Int; level::Int = 0)
    tmpbuf = IOBuffer()

    child = children(t)[idx]

    cursor -= 1
    cur = cursor == -1
    if child isa Tree
        cursor = printTreeChild(tmpbuf, child, cursor, term_width, level = level)
    else
        # if there's a specially designed show method we fall back to that
        if showmethod(typeof(child)) ≠ showmethod(Any)
            cur ? print(buf, "[ ] ") : print(buf, "    ")
            b = IOBuffer()
            print(b, Text(io -> show(IOContext(io, :limit => true), MIME"text/plain"(), child)))
            s = join(limitLineLength(split(String(take!(b)), '\n'), term_width-(2*level + 10)), "\n"*indent(level))
            print(tmpbuf, s)
        else
            d = defaultrepr(child)
            if d isa Tree
                cursor = printTreeChild(tmpbuf, d, cursor, term_width, level = level)
            else
                b = IOBuffer()
                print(b, d)
                s = join(limitLineLength(split(String(take!(b)), '\n'), term_width-(2*level + 10)), "\n"*indent(level))
                print(tmpbuf, s)
            end
        end
    end

    str = String(take!(tmpbuf))

    str = join(split(str, '\n'), "\n"*indent(level))

    print(buf, str)

    cursor
end

function printMenu(out, m::Tree, cursor; init::Bool=false, level=0)
    buf = IOBuffer()

    if init
        m.pageoffset = 0
    else
        # move cursor to beginning of current menu
        print(buf, "\x1b[999D\x1b[$(m.lastHeight)A")
        # clear display until end of screen
        print(buf, "\x1b[0J")
    end

    term_width = REPL.Terminals.width(terminal)

    # print header
    tb = IOBuffer()
    print(tb, m.head)
    println(buf, join(limitLineLength([String(take!(tb))], term_width-(2*level + 10)), '\n'))
    cs = children(m)

    for i in 1:length(cs)
        print(buf, "\x1b[2K")

        cursor = writeChild(buf, m, i, cursor, term_width, level=level+1)

        # dont print an \r\n on the last line
        i != (m.pagesize+m.pageoffset) && print(buf, "\r\n")
    end

    str = String(take!(buf))

    m.lastHeight = count(c -> c == '\n', str)

    print(out, str)
    cursor
end

cancel(t::Tree) = nothing

request(m::Tree) = request(terminal, m)

function request(term::REPL.Terminals.TTYTerminal, m::Tree)
    global mem
    mem = Dict{Any, Any}()

    cursor = 0

    printMenu(term.out_stream, m, cursor, init=true)

    raw_mode_enabled = enableRawMode(term)
    raw_mode_enabled && print(term.out_stream, "\x1b[?25l") # hide the cursor
    try
        while true
            c = readKey(term.in_stream)

            currentItem, _ = findItem(m, cursor)

            if c == Int(ARROW_UP)
                if cursor > 0
                    cursor -= 1
                end
            elseif c == Int(ARROW_DOWN)
                cursor += 1
            elseif c == 13 # <enter>
                # will break if pick returns true
                currentItem isa Tree && toggle(currentItem)
            elseif c == UInt32('q')
                cancel(m)
                break
            elseif c == 3 # ctrl-c
                cancel(m)
                break
            else
                # will break if keypress returns true
                keypress(m, c) && break
            end

            printMenu(term.out_stream, m, cursor)
        end
    finally
        # always disable raw mode even even if there is an
        #  exception in the above loop
        if raw_mode_enabled
            print(term.out_stream, "\x1b[?25h") # unhide cursor
            disableRawMode(term)
        end
    end
    println(term.out_stream)

    return m.selected
end

function findItem(t::Tree, cursor)
    i = nothing
    for c in children(t)
        if cursor == 0
            return c, cursor
        end

        cursor -= 1

        if c isa Tree && c.expanded
            i, cursor = findItem(c, cursor)
        end

        if i ≠ nothing
            return i, cursor
        end
    end
    return i, cursor
end

function toggleall(t::Tree, expand)
    for child in children(t)
        if child isa Tree
            child.expanded = expand
            toggleall(child, expand)
        end
    end
end

function keypress(t::Tree, key::UInt32)
    if key == UInt32('e') || key == UInt32('E')
        toggleall(t, true)
    elseif key == UInt32('c') || key == UInt32('C')
        toggleall(t, false)
    end
    false # don't break
end