text:
using Magnus
using SimpleFields

# Calculation in the dipole approximation, linear
# polarization.
function calc{T<:AbstractFloat}(H₀::AbstractMatrix{T},
                                D::AbstractMatrix,
                                ψ₀::Vector{Complex{T}},
                                field::Field, ndt::Integer;
                                gobble!::Function = Ψ -> (),
                                observe::Function = (Ψ,i,τ,field) -> (),
                                observables::Vector{Symbol} = Symbol[],
                                mode = :cpu,
                                verbose = true,
                                magnus_kwargs...)
    if verbose
        H = H₀ + D
        @printf("Hamiltonian dimensions: %dx%d, sparsity: %03.3f%%\n",
                size(H,1), size(H,2),
                100(1.0-length(nonzeros(H))/length(H)))
    end

    H₀ *= one(Complex{T})
    D *= one(Complex{T})

    N = ceil(Int, ndt*field.tmax)

    f = t -> field(t/field.T)

    obss = Dict{AbstractString,Any}()
    for o in observables
        os = string(o)
        os in keys(observable_types) || error("Unknown observable, $(os)")
        obss[os] = observable_types[os]{typeof(field.tmax)}(N)
    end

    results = integrate(ψ₀, T(field.tmax*field.T), N,
                        H₀, f, D, -im;
                        mode = mode,
                        verbose = verbose,
                        magnus_kwargs...) do Ψ,i,τ
                            gobble!(Ψ)
                            observe(Ψ,i,τ,field)
                            for o in values(obss)
                                o(Ψ,i,τ,field)
                            end
                        end

    Dict("E" => real(diag(H₀)), "psi" => results[:V],
         "milliseconds" => results[:milliseconds],
         "performance" => results[:performance],
         [k => v.v for (k,v) in obss]...)
end

# Calculation in the dipole approximation, linear
# polarization. Integration is done in the eigenbasis of the atomic
# Hamiltonian.
function calc{T<:AbstractFloat,
              M<:AbstractMatrix}(E::Vector{Vector{T}},
                                 V::Vector{M},
                                 D::Operator,
                                 field::Field, ndt::Integer;
                                 observe::Function = (Ψ,i,τ,field) -> (),
                                 observables::Vector{Symbol} = Symbol[],
                                 cutoff::Real = 0,
                                 mask_ratio::Real = 0,
                                 mode = :cpu,
                                 verbose = true,
                                 magnus_kwargs...)
    npartial = length(E)
    H₀,D,gst = hamiltonian(E, V, D, cutoff)
    ψ₀ = zeros(Complex{T}, size(H₀,1))
    ψ₀[gst] = 1.0
    gobble! = Egobbler(H₀, npartial, mask_ratio*cutoff, mode)

    calc(H₀, D, ψ₀, field, ndt;
         gobble! = gobble!,
         observe = observe,
         observables = observables,
         mode = mode, verbose = verbose,
         magnus_kwargs...)
end

export calc


text:
using Test
using RenderingGeometry

function test_precision_conversion(type₁, type₂)
    @eval begin
        dataset = generate_set($type₁)
        converted_dataset = [convert($type₂, pair) for pair ∈ dataset]
        reconverted_dataset = [convert($type₁, pair) for pair ∈ converted_dataset]
        for (pair₁, pair₂, pair₃) ∈ zip(dataset, converted_dataset, reconverted_dataset)
            @test typeof(pair₁) == typeof(pair₃) == $type₁
            @test typeof(pair₂) == $type₂
        end
    end
end

function test_type_casting(type₁, type₂)
    @eval begin
        dataset = generate_set($type₁)
        converted_dataset = [convert($type₂, pair) for pair ∈ dataset]
        for (pair₁, pair₂) ∈ zip(dataset, converted_dataset)
            @test typeof(pair₂) <: $type₂
            for (index₁, index₂) ∈ zip(pair₁, pair₂)
                @test typeof(index₁) == typeof(index₂)
                @test isapprox(index₁, index₂; atol=1)
            end
        end
    end
end

function test_all_pair_conversions()
    @testset "Pair conversion tests" begin
        log_test_batch_begin("Pair conversion")
        @testset "Vector conversion tests: $name₁ → $name₂" for (type₁, name₁) ∈ pair_vector_iterable, (type₂, name₂) ∈ pair_vector_iterable
            test_precision_conversion(type₁, type₂)
            increment_progress()
        end

        @testset "Point conversion tests: $name₁ → $name₂" for (type₁, name₁) ∈ pair_point_iterable, (type₂, name₂) ∈ pair_point_iterable
            test_precision_conversion(type₁, type₂)
            increment_progress()
        end

        @testset "Special cases: $name₁ → $name₂" for ((type₁, name₁), (type₂, name₂)) ∈ zip(pair_vector_iterable, pair_point_iterable)
            test_type_casting(type₁, type₂)
            increment_progress()
        end

        @testset "Special cases: $name₁ → $name₂" for ((type₁, name₁), (type₂, name₂)) ∈ zip(pair_point_iterable, pair_vector_iterable)
            test_type_casting(type₁, type₂)
            increment_progress()
        end
        log_test_batch_end()
    end
end

ext: jl


text:
include("./LatticeQCD.jl")
using Plots
using .LatticeQCD
using Random
using Dates
using JLD
#using GFlops

if length(ARGS) == 0
    error("""
    Use input file:
    Like,
    julia analyze.jl parameters.jl
    """)
end

function runtest()
    plot_plaquette(ARGS[1])
    plot_polyakov(ARGS[1])
    plot_plaq_and_poly(ARGS[1])
end

@time runtest()


text:
function OpenMapRealVector()
    return OpenMapRealVector(())
end

function OpenMapRealVector(arg0::OpenMapRealVector)
    return OpenMapRealVector((OpenMapRealVector,), arg0)
end

function OpenMapRealVector(arg0::RealVector)
    return OpenMapRealVector((RealVector,), arg0)
end

function OpenMapRealVector(arg0::Vector{Double})
    return OpenMapRealVector((Vector{Double},), arg0)
end

function OpenMapRealVector(arg0::Vector{Double}, arg1::jdouble)
    return OpenMapRealVector((Vector{Double}, jdouble), arg0, arg1)
end

function OpenMapRealVector(arg0::Vector{jdouble})
    return OpenMapRealVector((Vector{jdouble},), arg0)
end

function OpenMapRealVector(arg0::Vector{jdouble}, arg1::jdouble)
    return OpenMapRealVector((Vector{jdouble}, jdouble), arg0, arg1)
end

function OpenMapRealVector(arg0::jint)
    return OpenMapRealVector((jint,), arg0)
end

function OpenMapRealVector(arg0::jint, arg1::jdouble)
    return OpenMapRealVector((jint, jdouble), arg0, arg1)
end

function OpenMapRealVector(arg0::jint, arg1::jint)
    return OpenMapRealVector((jint, jint), arg0, arg1)
end

function OpenMapRealVector(arg0::jint, arg1::jint, arg2::jdouble)
    return OpenMapRealVector((jint, jint, jdouble), arg0, arg1, arg2)
end

function add(obj::OpenMapRealVector, arg0::OpenMapRealVector)
    return jcall(obj, "add", OpenMapRealVector, (OpenMapRealVector,), arg0)
end

function add(obj::OpenMapRealVector, arg0::RealVector)
    return jcall(obj, "add", RealVector, (RealVector,), arg0)
end

function add_to_entry(obj::RealVector, arg0::jint, arg1::jdouble)
    return jcall(obj, "addToEntry", void, (jint, jdouble), arg0, arg1)
end

function append(obj::OpenMapRealVector, arg0::OpenMapRealVector)
    return jcall(obj, "append", OpenMapRealVector, (OpenMapRealVector,), arg0)
end

function append(obj::OpenMapRealVector, arg0::RealVector)
    return jcall(obj, "append", OpenMapRealVector, (RealVector,), arg0)
end

function append(obj::OpenMapRealVector, arg0::jdouble)
    return jcall(obj, "append", OpenMapRealVector, (jdouble,), arg0)
end

function combine(obj::RealVector, arg0::jdouble, arg1::jdouble, arg2::RealVector)
    return jcall(obj, "combine", RealVector, (jdouble, jdouble, RealVector), arg0, arg1, arg2)
end

function combine_to_self(obj::RealVector, arg0::jdouble, arg1::jdouble, arg2::RealVector)
    return jcall(obj, "combineToSelf", RealVector, (jdouble, jdouble, RealVector), arg0, arg1, arg2)
end

function copy(obj::OpenMapRealVector)
    return jcall(obj, "copy", OpenMapRealVector, ())
end

function dot_product(obj::RealVector, arg0::RealVector)
    return jcall(obj, "dotProduct", jdouble, (RealVector,), arg0)
end

function ebe_divide(obj::OpenMapRealVector, arg0::RealVector)
    return jcall(obj, "ebeDivide", OpenMapRealVector, (RealVector,), arg0)
end

function ebe_multiply(obj::OpenMapRealVector, arg0::RealVector)
    return jcall(obj, "ebeMultiply", OpenMapRealVector, (RealVector,), arg0)
end

function equals(obj::OpenMapRealVector, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_dimension(obj::OpenMapRealVector)
    return jcall(obj, "getDimension", jint, ())
end

function get_distance(obj::OpenMapRealVector, arg0::OpenMapRealVector)
    return jcall(obj, "getDistance", jdouble, (OpenMapRealVector,), arg0)
end

function get_distance(obj::OpenMapRealVector, arg0::RealVector)
    return jcall(obj, "getDistance", jdouble, (RealVector,), arg0)
end

function get_entry(obj::OpenMapRealVector, arg0::jint)
    return jcall(obj, "getEntry", jdouble, (jint,), arg0)
end

function get_l1_distance(obj::OpenMapRealVector, arg0::OpenMapRealVector)
    return jcall(obj, "getL1Distance", jdouble, (OpenMapRealVector,), arg0)
end

function get_l1_distance(obj::OpenMapRealVector, arg0::RealVector)
    return jcall(obj, "getL1Distance", jdouble, (RealVector,), arg0)
end

function get_l1_norm(obj::RealVector)
    return jcall(obj, "getL1Norm", jdouble, ())
end

function get_l_inf_distance(obj::OpenMapRealVector, arg0::RealVector)
    return jcall(obj, "getLInfDistance", jdouble, (RealVector,), arg0)
end

function get_l_inf_norm(obj::RealVector)
    return jcall(obj, "getLInfNorm", jdouble, ())
end

function get_norm(obj::RealVector)
    return jcall(obj, "getNorm", jdouble, ())
end

function get_sparsity(obj::OpenMapRealVector)
    return jcall(obj, "getSparsity", jdouble, ())
end

function get_sub_vector(obj::OpenMapRealVector, arg0::jint, arg1::jint)
    return jcall(obj, "getSubVector", OpenMapRealVector, (jint, jint), arg0, arg1)
end

function hash_code(obj::OpenMapRealVector)
    return jcall(obj, "hashCode", jint, ())
end

function is_infinite(obj::OpenMapRealVector)
    return jcall(obj, "isInfinite", jboolean, ())
end

function is_na_n(obj::OpenMapRealVector)
    return jcall(obj, "isNaN", jboolean, ())
end

function map(obj::RealVector, arg0::UnivariateFunction)
    return jcall(obj, "map", RealVector, (UnivariateFunction,), arg0)
end

function map_add(obj::OpenMapRealVector, arg0::jdouble)
    return jcall(obj, "mapAdd", OpenMapRealVector, (jdouble,), arg0)
end

function map_add_to_self(obj::OpenMapRealVector, arg0::jdouble)
    return jcall(obj, "mapAddToSelf", OpenMapRealVector, (jdouble,), arg0)
end

function map_divide_to_self(obj::RealVector, arg0::jdouble)
    return jcall(obj, "mapDivideToSelf", RealVector, (jdouble,), arg0)
end

function map_multiply_to_self(obj::RealVector, arg0::jdouble)
    return jcall(obj, "mapMultiplyToSelf", RealVector, (jdouble,), arg0)
end

function map_subtract_to_self(obj::RealVector, arg0::jdouble)
    return jcall(obj, "mapSubtractToSelf", RealVector, (jdouble,), arg0)
end

function map_to_self(obj::RealVector, arg0::UnivariateFunction)
    return jcall(obj, "mapToSelf", RealVector, (UnivariateFunction,), arg0)
end

function outer_product(obj::RealVector, arg0::RealVector)
    return jcall(obj, "outerProduct", RealMatrix, (RealVector,), arg0)
end

function set(obj::OpenMapRealVector, arg0::jdouble)
    return jcall(obj, "set", void, (jdouble,), arg0)
end

function set_entry(obj::OpenMapRealVector, arg0::jint, arg1::jdouble)
    return jcall(obj, "setEntry", void, (jint, jdouble), arg0, arg1)
end

function set_sub_vector(obj::OpenMapRealVector, arg0::jint, arg1::RealVector)
    return jcall(obj, "setSubVector", void, (jint, RealVector), arg0, arg1)
end

function sparse_iterator(obj::OpenMapRealVector)
    return jcall(obj, "sparseIterator", JIterator, ())
end

function subtract(obj::OpenMapRealVector, arg0::OpenMapRealVector)
    return jcall(obj, "subtract", OpenMapRealVector, (OpenMapRealVector,), arg0)
end

function subtract(obj::OpenMapRealVector, arg0::RealVector)
    return jcall(obj, "subtract", RealVector, (RealVector,), arg0)
end

function to_array(obj::OpenMapRealVector)
    return jcall(obj, "toArray", Vector{jdouble}, ())
end

function unit_vector(obj::OpenMapRealVector)
    return jcall(obj, "unitVector", OpenMapRealVector, ())
end

function unitize(obj::OpenMapRealVector)
    return jcall(obj, "unitize", void, ())
end

function walk_in_default_order(obj::RealVector, arg0::RealVectorChangingVisitor)
    return jcall(obj, "walkInDefaultOrder", jdouble, (RealVectorChangingVisitor,), arg0)
end

function walk_in_default_order(obj::RealVector, arg0::RealVectorChangingVisitor, arg1::jint, arg2::jint)
    return jcall(obj, "walkInDefaultOrder", jdouble, (RealVectorChangingVisitor, jint, jint), arg0, arg1, arg2)
end

function walk_in_default_order(obj::RealVector, arg0::RealVectorPreservingVisitor)
    return jcall(obj, "walkInDefaultOrder", jdouble, (RealVectorPreservingVisitor,), arg0)
end

function walk_in_default_order(obj::RealVector, arg0::RealVectorPreservingVisitor, arg1::jint, arg2::jint)
    return jcall(obj, "walkInDefaultOrder", jdouble, (RealVectorPreservingVisitor, jint, jint), arg0, arg1, arg2)
end

function walk_in_optimized_order(obj::RealVector, arg0::RealVectorChangingVisitor)
    return jcall(obj, "walkInOptimizedOrder", jdouble, (RealVectorChangingVisitor,), arg0)
end

function walk_in_optimized_order(obj::RealVector, arg0::RealVectorChangingVisitor, arg1::jint, arg2::jint)
    return jcall(obj, "walkInOptimizedOrder", jdouble, (RealVectorChangingVisitor, jint, jint), arg0, arg1, arg2)
end

function walk_in_optimized_order(obj::RealVector, arg0::RealVectorPreservingVisitor)
    return jcall(obj, "walkInOptimizedOrder", jdouble, (RealVectorPreservingVisitor,), arg0)
end

function walk_in_optimized_order(obj::RealVector, arg0::RealVectorPreservingVisitor, arg1::jint, arg2::jint)
    return jcall(obj, "walkInOptimizedOrder", jdouble, (RealVectorPreservingVisitor, jint, jint), arg0, arg1, arg2)
end


text:
#! /usr/bin/julia

# Rosetta Code, Spiral matrix

immutable Spiral
    m::Int
    n::Int
    cmax::Int
    dir::Array{Array{Int,1},1}
    bdelta::Array{Array{Int,1},1}
end

function Spiral(m::Int, n::Int)
    cmax = m*n
    dir = Array{Int,1}[[0,1], [1,0], [0,-1], [-1,0]]
    bdelta = Array{Int,1}[[0,0,0,1], [-1,0,0,0],
                          [0,-1,0,0], [0,0,1,0]]
    Spiral(m, n, cmax, dir, bdelta)
end

function spiral(m::Int, n::Int)
    0<m&&0<n || error("The matrix dimensions must be positive.")
    Spiral(m, n)
end
spiral(n::Int) = spiral(n, n)

type SpState
    cnt::Int
    dirdex::Int
    cell::Array{Int,1}
    bounds::Array{Int,1}
end

Base.length(sp::Spiral) = sp.cmax
Base.start(sp::Spiral) = SpState(1, 1, [1,1], [sp.n,sp.m,1,1])
Base.done(sp::Spiral, sps::SpState) = sps.cnt > sp.cmax

function Base.next(sp::Spiral, sps::SpState)
    s = sub2ind((sp.m, sp.n), sps.cell[1], sps.cell[2])
    if sps.cell[rem1(sps.dirdex+1, 2)] == sps.bounds[sps.dirdex]
        sps.bounds += sp.bdelta[sps.dirdex]
        sps.dirdex = rem1(sps.dirdex+1, 4)
    end
    sps.cell += sp.dir[sps.dirdex]
    sps.cnt += 1
    return (s, sps)
end

using Formatting

function width{T<:Integer}(n::T)
    w = ndigits(n)
    n < 0 || return w
    return w + 1
end

function pretty{T<:Integer}(a::Array{T,2}, indent::Int=4)
    lo, hi = extrema(a)
    w = max(width(lo), width(hi))
    id = " "^indent
    fe = FormatExpr(@sprintf(" {:%dd}", w))
    s = id
    nrow = size(a)[1]
    for i in 1:nrow
        for j in a[i,:]
            s *= format(fe, j)
        end
        i != nrow || continue
        s *= "\n"*id
    end
    return s
end

n = 5
println("The n = ", n, " spiral matrix:")
a = zeros(Int, (n, n))
for (i, s) in enumerate(spiral(n))
    a[s] = i-1
end
println(pretty(a))

m = 3
println()
println("Generalize to a non-square matrix (", m, "x", n, "):")
a = zeros(Int, (m, n))
for (i, s) in enumerate(spiral(m, n))
    a[s] = i-1
end
println(pretty(a))

p = primes(10^3)
n = 7
println()
println("An n = ", n, " prime spiral matrix:")
a = zeros(Int, (n, n))
for (i, s) in enumerate(spiral(n))
    a[s] = p[i]
end
println(pretty(a))

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
include("../src/overapprox_nd.jl")
using Test
# tests for overest_nd.jl

@testset "test for find_variables" begin
    @test find_variables(:(x+1)) == [:x]
    @test find_variables(:(log(x + y*z))) == [:x, :y, :z]
end

@testset "test for is_affine" begin
    @test is_affine(:(x+1)) == true
    @test is_affine(:(-x+(2y-z))) == true
    @test is_affine(:(log(x))) == false
    @test is_affine(:(x+ xz)) == true # interprets xz as one variable
    @test is_affine(:(x+ x*z)) == false
    @test is_affine(:(x + y*log(2))) == true
end

@testset "test for is_1d" begin
    @test is_1d(:(x*log(2))) == true
    @test is_1d(:(x^2.5)) == true
end

@testset "test for is_unary" begin
    @test is_unary(:(x+y)) == false
    @test is_unary(:(x^2.5)) == false # this is like pow(x, 2.5), hence not unary. is taken care of tho
    @test is_unary(:(sin(x))) == true
    @test is_unary(:(x*y)) == false
end
# test find_UB (really tests for overapprox_1d.jl)

@testset "test find_affine_range" begin
    @test find_affine_range(:(x + 1), Dict(:x => (0, 1))) == (1,2)
    @test find_affine_range(:(- x + 1), Dict(:x => (0, 1))) == (0,1)
    @test find_affine_range(:((-1 + y) + x), Dict(:x => (0,1), :y => (1,2))) == (0,2)
    @test find_affine_range(:(2*x), Dict(:x => (0,1))) == (0,2)
    @test find_affine_range(:(x*2), Dict(:x => (0,1))) == (0,2)
    @test find_affine_range(:(x*π), Dict(:x => (0,1)))[2] ≈ π
    @test find_affine_range(:(x*(π*log(2)) + π), Dict(:x => (0,1)))[2] ≈ π*log(2) + π
end

@testset "test substitute" begin
    @test substitute!(:(x^2+1), :x, :(y+1)) == :((y+1)^2+1)
end

@testset "test count min max" begin
    @test count_min_max(:(min(x) + max(y - min(x)))) == [2,1]
end

# test upperbound_expr_compositions
# need clarity on precise purpose and signature of the function

@testset "test reduce_args_to_2!" begin
    @test reduce_args_to_2(:(x+y+z)) == :(x+(y+z))
    @test reduce_args_to_2(:(sin(x*y*z))) == :(sin(x*(y*z)))
end

N_VARS = 0

@assert add_var(OverApproximation()) == :v_1

@assert apply_fx(:(x + y), :a) == :(a + y)

@assert is_number(:(log(2)/3))

@assert is_unary(:(sin(x)))
@assert !is_unary(:(x + y))

@assert is_binary(:(x + y))

function affine_tests()
    @assert is_affine(:(x+1))
    @assert is_affine(:(-x+(2y-z)))
    @assert !is_affine(:(log(x)))
    @assert !is_affine(:(x + x*z))
    @assert is_affine(:(x/6))
    @assert is_affine(:(5*x))
    @assert is_affine(:((1 / 6) * x))
    @assert is_affine(:(log(2)*x))
    @assert is_affine(:(-x))
end
affine_tests()

@testset "PWL_tests" begin
    @test is_relu(:(relu(x)))
    @test is_relu(:(max(x, 0)))
    @test is_relu(:(max(0, 5*x + log(2))))
    @test is_relu(:(max(x*y, 0)))
    @test !is_relu(:(max(x, y)))
    @test !is_relu(:(min(0, x)))
    @test find_range(:(max(x, y)), Dict(:x => [4.,5.], :y=>[1., 6.])) == [4., 6.]
    @test find_range(:(max(5*x, y)), Dict(:x => [4.,5.], :y=>[1., 6.])) == [20., 25.]
    @test find_range(:(min(y, -5*x)), Dict(:x => [4.,5.], :y=>[1., 6.])) == [-25., -20]
    @test find_range(:(max(x, 0)), Dict(:x => [-6.7, 8.9])) == [0.0, 8.9]
    @test find_range(:(relu(x)), Dict(:x => [-12.3, 35.])) == [0.0, 35]
end

# function outer_affine_tests()
#     @assert is_outer_affine(:(-sin(x) + sin(y)))
#     @assert is_outer_affine(:( -5*(sin(x) - 3*sin(y)) ) )
# end
# outer_affine_tests()

overapprox(:(sin(x)), Dict(:x=>[0,π/2]))

overapprox(:(sin(sin(x))), Dict(:x=>[0,π/2]))

# improve how this is handled
overapprox(:(sin(sin(z + y))), Dict(:z=>[0,π], :y=>[-π,π]))

overapprox(:(2*x), Dict(:x=>[0,1]))

# This isn't reduced to eval(log(2))*x because whole thing is affine...
overapprox(:(log(2)*x), Dict(:x=>[0.,1.]))

overapprox(:(2*log(x)), Dict(:x=>[1,2]))

expr = :(exp(2*sin(x) + y) - log(6)*z)
init_set = Dict{Symbol, Array{Float64,1}}(:x=>[0,1], :z=>[0,π], :y=>[-π,π])
oa = overapprox(expr, init_set)

#
b1 = OverApproximation()
b1.ranges = Dict(:A=>[1.,2.], :B=>[3.,4.])
b1.nvars = 2
expand_multiplication(:A, :B, b1; ξ=0.1)
# reference: (:(exp(log(C) + log(D)) + (1.0 - 0.1) * D + (3.0 - 0.1) * C + (1.0 - 0.1) * (3.0 - 0.1)), Dict(:A => (1.0, 2.0),:D => (0.1, 1.1),:B => (3.0, 4.0),:C => (0.1, 1.1)), Expr[:(C = (A - 1.0) + 0.1), :(D = (B - 3.0) + 0.1)])

bound_binary_functions(:*, :A, :B, b1)

xxx = 0; yyy = 1; zzz = 2; kkk = -10;
e1 = :(xxx + yyy + zzz + kkk - xxx - yyy - zzz - kkk)
e2 = reduce_args_to_2(e1)
eval(e1) == eval(e2)

overapprox(:(x*y), Dict(:x=>[1,2], :y=>[-10,-9]))

overapprox(:(sin(6) + sin(x)), Dict(:x=>[1,2]))

overapprox(:(sin(6)*sin(x)), Dict(:x=>[1,2]))

overapprox(:(sin(6)*sin(x)*sin(y)), Dict(:x=>[1,2], :y=>[1,2])::Dict{Symbol,Array{Int64,1}})

# TODO: when the whole expression is affine, division by a constant still appears...fix so it's converted to multiplication by 1/c
overapprox(:(x/6), Dict(:x=>[-1,1]))
# should be handled by affine base case

overapprox(:(6/x), Dict(:x=>[1,2]))

overapprox(:(6/(x+y)), Dict(:x=>[2,3], :y=>[-1,2]))

overapprox(:(x/y), Dict(:x=>[2,3], :y=>[1,2]))

overapprox(:(sin(x + y)/6), Dict(:x=>[2,3], :y=>[1,2]))

overapprox(:(sin(x + y)/y), Dict(:x=>[2,3], :y=>[1,2]))

overapprox(:(exp(x^2)), Dict(:x=>[-1,1]))

overapprox(:((x+sin(y))^3), Dict(:x=>[2,3], :y=>[1,2]))

overapprox(:(2^x), Dict(:x=>[2,3]))

overapprox(:(-sin(x+y)), Dict(:x=>[2,3], :y=>[1,2]))

overapprox(:(log(x)), Dict(:x => [1.0, 166.99205596346707]); N=-1)

overapprox(:(relu(x)), Dict(:x => [-5.7, 23.5]); N=1)

overapprox(:(max(x,y)), Dict(:x=>[2,3.], :y=>[1,2.]))

overapprox(:(max(sin(x),5*y - 6)), Dict(:x=>[2,3.], :y=>[1,2.]))

overapprox(:(max(min(x, y),5*y - 6*cos(y))), Dict(:x=>[2,3.], :y=>[1,2.]))


text:
function set_100()
    set = Dict{Part, Int}()

    # Parts 1
    for c in [red, green, yellow, blue]
        push!(set, part(c, 1) => 10)
    end

    # Parts 2
    for c in [blue, yellow]
        push!(set, part(c, 2) => 1)
    end

    return set
end


text:
julia> [1 2 3 ; 4 5 6] * [1 2 ; 3 4 ; 5 6]  # product of a 2x3 by a 3x2
2x2 Array{Int64,2}:
 22  28
 49  64

julia> [1 2 3] * [1,2,3]   # product of a row vector by a column vector
1-element Array{Int64,1}:
 14

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module TuringTutorials

import IOCapture

using IJulia
using InteractiveUtils
using Pkg
using Plots
using Requires
using Weave

const REPO_DIR = string(pkgdir(TuringTutorials))::String

include("cache.jl")
include("build.jl")

export build_folder, tutorial_path, folder2filename
export build, build_and_exit, verify_logs, tutorials, changed_tutorials

# Not building PDF, because it is fragile. Maybe later.
default_build_list = (:script, :html, :github, :notebook)

cssfile = joinpath(@__DIR__, "..", "templates", "skeleton_css.css")
latexfile = joinpath(@__DIR__, "..", "templates", "julia_tex.tpl")

function __init__()
    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("weaveplots.jl")
end

"""
    tutorial_path(folder)

Return the absolute path to a tutorial `folder`.
"""
tutorial_path(folder) = joinpath(REPO_DIR, "tutorials", folder)

function polish_latex(path::String)
    # TODO: Is it maybe better to overload https://github.com/JunoLab/Weave.jl/blob/b5ba227e757520f389a6d6e0f2cacb731eab8b12/src/WeaveMarkdown/markdown.jl#L10-L17
    # and replace the `tex.formula` there? Only negative part is that this of course
    # will affect all of the markdown parsers, which is not necessarily desirable.
    txt = open(f -> read(f, String), path)
    open(path, "w+") do f
        txt = replace(txt, "\\\\\n" => "\\\\\\\\\n")
        write(f, txt)
    end
end

function weave_file(
    folder, file, build_list=default_build_list;
    kwargs...
)
    tmp = joinpath(tutorial_path(folder), file)
    Pkg.activate(dirname(tmp))
    Pkg.instantiate()
    args = Dict{Symbol,String}(:folder => folder, :file => file)
    if :script ∈ build_list
        println("Building Script")
        dir = joinpath(REPO_DIR, "script", folder)
        isdir(dir) || mkpath(dir)
        args[:doctype] = "script"
        tangle(tmp;out_path=dir)
    end
    if :html ∈ build_list
        println("Building HTML")
        dir = joinpath(REPO_DIR, "html", folder)
        isdir(dir) || mkpath(dir)
        args[:doctype] = "html"
        weave(
            tmp, doctype = "md2html", out_path=dir, args=args;
            fig_ext=".svg", css=cssfile, kwargs...
        )
    end
    if :pdf ∈ build_list
        println("Building PDF")
        dir = joinpath(REPO_DIR, "pdf", folder)
        isdir(dir) || mkpath(dir)
        args[:doctype] = "pdf"
        try
            weave(
                tmp, doctype="md2pdf", out_path=dir, args=args;
                template=latexfile, kwargs...
            )
        catch ex
            @warn "PDF generation failed" exception=(ex, catch_backtrace())
        end
    end
    if :github ∈ build_list
        println("Building Github Markdown")
        dir = joinpath(REPO_DIR, "markdown", folder)
        isdir(dir) || mkpath(dir)
        args[:doctype] = "github"
        out_path = weave(tmp,doctype = "github",out_path=dir,args=args; kwargs...)
        polish_latex(out_path)
    end
    if :notebook ∈ build_list
        println("Building Notebook")
        dir = joinpath(REPO_DIR, "notebook", folder)
        isdir(dir) || mkpath(dir)
        args[:doctype] = "notebook"
        Weave.convert_doc(tmp, joinpath(dir, file[1:end-4]*".ipynb"))
    end
end

"""
    tutorials::Vector{String}

Return names of the tutorials.
"""
function tutorials()::Vector{String}
    dirs = readdir(joinpath(REPO_DIR, "tutorials"))
    dirs = filter(!=("test.jmd"), dirs)
    # This DiffEq one has to be done manually, because it takes about 12 hours.
    dirs = filter(!=("10-bayesian-differential-equations"), dirs)
    dirs = filter(!=("99-test"), dirs)
end

function weave_all(build_list=default_build_list; kwargs...)
    for tutorial in tutorials()
        weave_folder(folder, build_list; kwargs...)
    end
end

function weave_md(; kwargs...)
    for tutorial in tutorials()
        weave_folder(folder, (:github,); kwargs...)
    end
end

function weave_folder(
    folder, build_list=default_build_list;
    ext = r"^\.[Jj]md", kwargs...
)
    for file in readdir(tutorial_path(folder))
        try
            # HACK: only weave (j)md files
            if occursin(ext, splitext(file)[2])
                println("Building $(joinpath(folder, file))")
                weave_file(folder, file, build_list; kwargs...)
            else
                println("Skipping $(joinpath(folder, file))")
            end
        catch ex
            rethrow(ex)
        end
    end
end

function tutorial_footer(folder=nothing, file=nothing; remove_homedir=true)
    display("text/markdown", """
        ## Appendix
         This tutorial is part of the TuringTutorials repository, found at: <https://github.com/TuringLang/TuringTutorials>.
        """)
    if folder !== nothing && file !== nothing
        display("text/markdown", """
            To locally run this tutorial, do the following commands:
            ```julia, eval = false
            using TuringTutorials
            TuringTutorials.weave_file("$folder", "$file")
            ```
            """)
    end
    display("text/markdown", "Computer Information:")
    vinfo = sprint(InteractiveUtils.versioninfo)
    display("text/markdown",  """
        ```
        $(vinfo)
        ```
        """)

    ctx = Pkg.API.Context()

    pkg_status = let io = IOBuffer()
        Pkg.status(Pkg.API.Context(); io = io)
        String(take!(io))
    end
    projfile = ctx.env.project_file
    remove_homedir && (projfile = replace(projfile, homedir() => "~"))

    display("text/markdown","""
        Package Information:
        """)

    md = "```\n$(pkg_status)\n```"
    display("text/markdown", md)
end
end # module

max_stars_repo_name: staticfloat/TuringTutorials

max_issues_repo_name: staticfloat/TuringTutorials

max_forks_repo_name: staticfloat/TuringTutorials


text:
using SphericalHarmonics
using SphericalHarmonicModes
using Test
using HCubature
using Aqua
using LegendrePolynomials
using WignerD
using OffsetArrays
using MultiplesOfPi
using ForwardDiff

using SphericalHarmonics: NorthPole,
SouthPole,
allocate_y,
allocate_p,
RealHarmonics,
ComplexHarmonics,
_lmax,
associatedLegendre,
associatedLegendrex,
Unnormalized,
LMnorm,
Orthonormal

@testset "project quality" begin
    Aqua.test_all(SphericalHarmonics, ambiguities = false)
end

@testset "allocate" begin
    lmax = 4
    @test size(allocate_y(lmax)) == size(allocate_y(ComplexF64, lmax))
    @test eltype(allocate_y(Complex{BigFloat}, lmax)) == Complex{BigFloat}
    @test eltype(allocate_y(BigFloat, lmax)) == BigFloat

    @test size(allocate_p(lmax)) == size(allocate_p(Float64, lmax))
    @test eltype(allocate_p(BigFloat, lmax)) == BigFloat
end

@testset "Ylm and Plm explicit" begin

    function explicit_shs_complex_fullrange(θ, φ)
        Y00 = 0.5 * sqrt(1/π)
        Y1m1 = 0.5 * sqrt(3/(2*π))*sin(θ)*exp(-im*φ)
        Y10 = 0.5 * sqrt(3/π)*cos(θ)
        Y11 = -0.5 * sqrt(3/(2*π))*sin(θ)*exp(im*φ)
        Y2m2 = 0.25 * sqrt(15/(2*π))*sin(θ)^2*exp(-2*im*φ)
        Y2m1 = 0.5 * sqrt(15/(2*π))*sin(θ)*cos(θ)*exp(-im*φ)
        Y20 = 0.25 * sqrt(5/π)*(3*cos(θ)^2 - 1)
        Y21 = -0.5 * sqrt(15/(2*π))*sin(θ)*cos(θ)*exp(im*φ)
        Y22 = 0.25 * sqrt(15/(2*π))*sin(θ)^2*exp(2*im*φ)
        Y3m3 = 1/8 * exp(-3 * im * φ) * sqrt(35/π) * sin(θ)^3
        Y3m2 = 1/4 * exp(-2 * im * φ) * sqrt(105/(2*π)) * cos(θ) * sin(θ)^2
        Y3m1 = 1/8 * exp(-im * φ) * sqrt(21/π) * (-1 + 5 * cos(θ)^2) * sin(θ)
        Y30 = 1/4 * sqrt(7/π) * (-3 * cos(θ) + 5 * cos(θ)^3)
        Y31 = -(1/8) * exp(im * φ) * sqrt(21/π) * (-1 + 5 * cos(θ)^2) * sin(θ)
        Y32 = 1/4 * exp(2 * im * φ) * sqrt(105/(2*π)) * cos(θ) * sin(θ)^2
        Y33 = -(1/8) * exp(3 * im * φ) * sqrt(35/π) * sin(θ)^3
        return [Y00, Y1m1, Y10, Y11, Y2m2, Y2m1, Y20, Y21, Y22,
               Y3m3, Y3m2, Y3m1, Y30, Y31, Y32, Y33]
    end

    function explicit_shs_complex_zeroto(θ, φ)
        Y00 = 0.5 * sqrt(1/π)
        Y10 = 0.5 * sqrt(3/π)*cos(θ)
        Y11 = -0.5 * sqrt(3/(2*π))*sin(θ)*exp(im*φ)
        Y20 = 0.25 * sqrt(5/π)*(3*cos(θ)^2 - 1)
        Y21 = -0.5 * sqrt(15/(2*π))*sin(θ)*cos(θ)*exp(im*φ)
        Y22 = 0.25 * sqrt(15/(2*π))*sin(θ)^2*exp(2*im*φ)
        Y30 = 1/4 * sqrt(7/π) * (-3 * cos(θ) + 5 * cos(θ)^3)
        Y31 = -(1/8) * exp(im * φ) * sqrt(21/π) * (-1 + 5 * cos(θ)^2) * sin(θ)
        Y32 = 1/4 * exp(2 * im * φ) * sqrt(105/(2*π)) * cos(θ) * sin(θ)^2
        Y33 = -(1/8) * exp(3 * im * φ) * sqrt(35/π) * sin(θ)^3
        return [Y00, Y10, Y11, Y20, Y21, Y22, Y30, Y31, Y32, Y33]
    end

    function explicit_associated_legendre(θ)
        # test only the unnormalized polynomials
        x = cos(θ)
        P00 = 1
        P10 = x
        P11 = -sin(θ)
        P20 = (3x^2 - 1)/2
        P21 = -3x*sin(θ)
        P22 = 3(1-x^2)

        [P00, P10, P11, P20, P21, P22]
    end

    @testset "associated legendre polynomials" begin
        for θ in LinRange(0, π, 100)
            P = computePlmcostheta(θ, 2, norm = Unnormalized())
            @test P ≈ explicit_associated_legendre(θ)
        end

        @testset "condon-shortley phase" begin
            θ = pi/3
            for normtype in [Unnormalized, LMnorm, Orthonormal]
                P = computePlmcostheta(θ, 2, norm = normtype())
                PnoCS = computePlmcostheta(θ, 2, norm = normtype(csphase = false))
                for (ind, (l,m)) in enumerate(ML(0:2, ZeroTo))
                    @test P[ind] ≈ (-1)^m * PnoCS[ind] atol=1e-14 rtol=1e-14
                end
            end
        end
    end

    @testset "complex harmonics" begin

        for θ in LinRange(0, π, 100), ϕ in LinRange(0, 2π, 200)
            Y = computeYlm(θ, ϕ, 3)
            Yex = explicit_shs_complex_fullrange(θ, ϕ)
            for (ind,mode) in enumerate(ML(0:3))
                @test Y[mode] ≈ Yex[ind]
            end

            Y = computeYlm(θ, ϕ, 3, ZeroTo)
            Yex = explicit_shs_complex_zeroto(θ, ϕ)
            for (ind,mode) in enumerate(ML(0:3, ZeroTo))
                @test Y[mode] ≈ Yex[ind]
            end
        end
    end

    @testset "real harmonics" begin
        function explicit_shs_real_fullrange(θ, φ)
            Y00 = 0.5 * sqrt(1/π)
            Y1m1 = -0.5 * sqrt(3/π)*sin(θ)*sin(φ)
            Y10 = 0.5 * sqrt(3/π)*cos(θ)
            Y11 = -0.5 * sqrt(3/π)*sin(θ)*cos(φ)
            Y2m2 = 0.25 * sqrt(15/π)*sin(θ)^2*sin(2*φ)
            Y2m1 = -0.5 * sqrt(15/π)*sin(θ)*cos(θ)*sin(φ)
            Y20 = 0.25 * sqrt(5/π)*(3*cos(θ)^2 - 1)
            Y21 = -0.5 * sqrt(15/π)*sin(θ)*cos(θ)*cos(φ)
            Y22 = 0.25 * sqrt(15/π)*sin(θ)^2*cos(2*φ)
            Y3m3 = -1/8 * sin(3 * φ) * sqrt(2*35/π) * sin(θ)^3
            Y3m2 = 1/4 * sin(2 * φ) * sqrt(105/π) * cos(θ) * sin(θ)^2
            Y3m1 = -1/8 * sin(φ) * sqrt(2*21/π) * (-1 + 5 * cos(θ)^2) * sin(θ)
            Y30 = 1/4 * sqrt(7/π) * (-3 * cos(θ) + 5 * cos(θ)^3)
            Y31 = -(1/8) * cos(φ) * sqrt(2*21/π) * (-1 + 5 * cos(θ)^2) * sin(θ)
            Y32 = 1/4 * cos(2 * φ) * sqrt(105/π) * cos(θ) * sin(θ)^2
            Y33 = -(1/8) * cos(3 * φ) * sqrt(2*35/π) * sin(θ)^3
            return [Y00, Y1m1, Y10, Y11, Y2m2, Y2m1, Y20, Y21, Y22,
                   Y3m3, Y3m2, Y3m1, Y30, Y31, Y32, Y33]
        end

        function explicit_shs_real_zeroto(θ, φ)
            Y00 = 0.5 * sqrt(1/π)
            Y10 = 0.5 * sqrt(3/π)*cos(θ)
            Y11 = -0.5 * sqrt(3/π)*sin(θ)*cos(φ)
            Y20 = 0.25 * sqrt(5/π)*(3*cos(θ)^2 - 1)
            Y21 = -0.5 * sqrt(15/π)*sin(θ)*cos(θ)*cos(φ)
            Y22 = 0.25 * sqrt(15/π)*sin(θ)^2*cos(2*φ)
            Y30 = 1/4 * sqrt(7/π) * (-3 * cos(θ) + 5 * cos(θ)^3)
            Y31 = -(1/8) * cos(φ) * sqrt(2*21/π) * (-1 + 5 * cos(θ)^2) * sin(θ)
            Y32 = 1/4 * cos(2 * φ) * sqrt(105/π) * cos(θ) * sin(θ)^2
            Y33 = -(1/8) * cos(3 * φ) * sqrt(2*35/π) * sin(θ)^3
            return [Y00, Y10, Y11, Y20, Y21, Y22, Y30, Y31, Y32, Y33]
        end

        for θ in LinRange(0, π, 100), ϕ in LinRange(0, 2π, 100)
            Y = computeYlm(θ, ϕ, 3, FullRange, SphericalHarmonics.RealHarmonics())
            Yex = explicit_shs_real_fullrange(θ, ϕ)
            for (ind,mode) in enumerate(ML(0:3))
                @test Y[mode] ≈ Yex[ind]
            end

            Y = computeYlm(θ, ϕ, 3, ZeroTo, SphericalHarmonics.RealHarmonics())
            Yex = explicit_shs_real_zeroto(θ, ϕ)
            for (ind,mode) in enumerate(ML(0:3, ZeroTo))
                @test Y[mode] ≈ Yex[ind]
            end
        end
    end

    @testset "Pole" begin
        @testset "NorthPole" begin
            Y = computeYlm(NorthPole(), 3)
            Yex = explicit_shs_complex_fullrange(0, 0)
            @test Y ≈ Yex

            Y = computeYlm(NorthPole(), 3, ZeroTo)
            Yex = explicit_shs_complex_zeroto(0, 0)
            @test Y ≈ Yex
        end
        @testset "SouthPole" begin
            Y = computeYlm(SouthPole(), 3)
            Yex = explicit_shs_complex_fullrange(π, 0)
            @test Y ≈ Yex
        end
        @testset "RealHarmonics == ComplexHarmonics" begin
            for x in (NorthPole(), SouthPole())
                YC = computeYlm(x, 4, FullRange, SphericalHarmonics.ComplexHarmonics())
                YR = computeYlm(x, 4, FullRange, SphericalHarmonics.RealHarmonics())
                @test YC == YR

                YC = computeYlm(x, 4, ZeroTo, SphericalHarmonics.ComplexHarmonics())
                YR = computeYlm(x, 4, ZeroTo, SphericalHarmonics.RealHarmonics())
                @test YC == YR
            end
        end
    end
end

@testset "computePlm!" begin
    lmax = 10
    θ = pi/3
    coeff = SphericalHarmonics.compute_coefficients(lmax);
    P = SphericalHarmonics.allocate_p(Float64, lmax);
    computePlmcostheta!(P, θ, lmax, coeff);

    @test P == computePlmcostheta(θ, lmax)

    computePlmcostheta!(P, NorthPole(), lmax, coeff);
    @test P == computePlmcostheta(NorthPole(), lmax)

    @testset "single m" begin
        computePlmcostheta!(P, θ, lmax, coeff);
        P2 = SphericalHarmonics.allocate_p(Float64, lmax);
        P3 = SphericalHarmonics.allocate_p(Float64, lmax);
        computePlmcostheta!(P3, θ, lmax, nothing, coeff);
        @test P3 ≈ P

        for l in 0:lmax, m in 0:l
            computePlmcostheta!(P2, θ, l, m, coeff);
            for l2 in m:l
                @test P2[(l2,m)] == P[(l2,m)]
            end
        end

        @test_throws ArgumentError computePlmcostheta!(P2, θ, 2, 3, coeff)
        @test_throws ArgumentError computePlmcostheta!(P2, θ, 2, -1, coeff)
    end
end

@testset "computePlmx and computePlmcostheta" begin
    θ = pi/3
    @test SphericalHarmonics.computePlmx(cos(θ), 4) ≈ computePlmcostheta(θ, 4)
    @test SphericalHarmonics.computePlmx(cos(θ), lmax = 4) ≈ computePlmcostheta(θ, 4)
    @test SphericalHarmonics.computePlmx(cos(θ), 4) ≈ computePlmcostheta(θ, lmax = 4)
    @test SphericalHarmonics.computePlmx(cos(θ), lmax = 4) ≈ computePlmcostheta(θ, lmax = 4)

    @test SphericalHarmonics.computePlmx(cos(θ), 4, 2) ≈ computePlmcostheta(θ, 4, 2)
    @test SphericalHarmonics.computePlmx(cos(θ), lmax = 4, m = 2) ≈ computePlmcostheta(θ, 4, 2)
    @test SphericalHarmonics.computePlmx(cos(θ), 4, 2) ≈ computePlmcostheta(θ, lmax = 4, m = 2)
    @test SphericalHarmonics.computePlmx(cos(θ), lmax = 4, m = 2) ≈ computePlmcostheta(θ, lmax = 4, m = 2)
end

@testset "computeYlm kwargs" begin
    θ, ϕ = pi/3, pi/3
    lmax = 10

    Y1 = SphericalHarmonics.computeYlm(θ, ϕ, lmax)
    Y2 = SphericalHarmonics.computeYlm(θ, ϕ, lmax = lmax)
    Y3 = SphericalHarmonics.computeYlm(θ, ϕ, lmax = lmax, m_range = FullRange)
    Y4 = SphericalHarmonics.computeYlm(θ, ϕ, lmax = lmax, m_range = FullRange, SHType = SphericalHarmonics.ComplexHarmonics())
    @test Y1 == Y2 == Y3 == Y4

    for m = -lmax:lmax
        Y1 = SphericalHarmonics.computeYlm(θ, ϕ, lmax, m)
        Y2 = SphericalHarmonics.computeYlm(θ, ϕ, lmax = lmax, m = m)
        Y3 = SphericalHarmonics.computeYlm(θ, ϕ, lmax = lmax, m = m, m_range = FullRange)
        Y4 = SphericalHarmonics.computeYlm(θ, ϕ, lmax = lmax, m = m, m_range = FullRange, SHType = SphericalHarmonics.ComplexHarmonics())
        @test Y1 == Y2 == Y3 == Y4
    end

    Y1 = SphericalHarmonics.computeYlm(θ, ϕ, lmax, ZeroTo)
    Y2 = SphericalHarmonics.computeYlm(θ, ϕ, lmax = lmax, m_range = ZeroTo)
    Y3 = SphericalHarmonics.computeYlm(θ, ϕ, lmax = lmax, m_range = ZeroTo, SHType = SphericalHarmonics.ComplexHarmonics())
    @test Y1 == Y2 == Y3

    for m in 0:lmax
        Y1 = SphericalHarmonics.computeYlm(θ, ϕ, lmax, m, ZeroTo)
        Y2 = SphericalHarmonics.computeYlm(θ, ϕ, lmax = lmax, m = m, m_range = ZeroTo)
        Y3 = SphericalHarmonics.computeYlm(θ, ϕ, lmax = lmax, m = m, m_range = ZeroTo, SHType = SphericalHarmonics.ComplexHarmonics())
        @test Y1 == Y2 == Y3
    end

    Y1 = SphericalHarmonics.computeYlm(θ, ϕ, lmax, ZeroTo, SphericalHarmonics.RealHarmonics())
    Y2 = SphericalHarmonics.computeYlm(θ, ϕ, lmax = lmax, m_range = ZeroTo, SHType = SphericalHarmonics.RealHarmonics())
    @test Y1 == Y2
end

@testset "computeYlm!" begin
    θ, ϕ = pi/3, pi/3
    lmax = 10
    P = SphericalHarmonics.computePlmcostheta(θ, lmax)
    Y1 = SphericalHarmonics.allocate_y(ComplexF64, lmax)
    Y2 = SphericalHarmonics.computeYlm(θ, ϕ, lmax)

    SphericalHarmonics.computeYlm!(Y1, θ, ϕ, lmax)
    @test Y1 ≈ Y2

    SphericalHarmonics.computeYlm!(Y1, θ, ϕ; lmax = lmax)
    @test Y1 ≈ Y2

    SphericalHarmonics.computeYlm!(Y1, P, θ, ϕ; lmax = lmax)
    @test Y1 ≈ Y2

    SphericalHarmonics.computeYlm!(Y1, θ, ϕ, lmax, ZeroTo)
    Y2 = SphericalHarmonics.computeYlm(θ, ϕ, lmax, ZeroTo)
    @test Y1[1:length(ML(ZeroTo(lmax),ZeroTo))] ≈ Y2

    @testset "single m" begin
        Y2 = SphericalHarmonics.computeYlm(θ, ϕ, lmax)
        Ym = SphericalHarmonics.allocate_y(ComplexF64, lmax)

        for l in 0:lmax, m in -l:l
            computeYlm!(Ym, P, θ, ϕ, l, m)
            for l2 in abs(m):l
                @test Ym[(l2,m)] == Y2[(l2,m)]
            end
        end

        Y2 = SphericalHarmonics.computeYlm(θ, ϕ; lmax = lmax, SHType = SphericalHarmonics.RealHarmonics())
        Ym = SphericalHarmonics.allocate_y(Float64, lmax)

        for l in 0:lmax, m in -l:l
            computeYlm!(Ym, P, θ, ϕ, l, m, SphericalHarmonicModes.FullRange, SphericalHarmonics.RealHarmonics())
            for l2 in abs(m):l
                @test Ym[(l2,m)] == Y2[(l2,m)]
            end
        end

        @test_throws ArgumentError computeYlm!(Ym, P, θ, ϕ, 2, 3)
        @test_throws ArgumentError computeYlm!(Ym, P, θ, ϕ, 2, -3)
    end
end

@testset "Pole" begin
    @testset "utils" begin
        @test one(NorthPole()) == 1
        @test one(SouthPole()) == 1
        @test zero(NorthPole()) == 0
        @test zero(SouthPole()) == 0

        @test float(NorthPole()) == 0
        @test float(SouthPole()) == float(pi)

        @test Float64(NorthPole()) == 0
        @test Float64(SouthPole()) == float(pi)

        @test promote_rule(SouthPole, Float64) == Float64
        @test promote_rule(SouthPole, BigFloat) == BigFloat
    end

    @testset "North pole" begin
        @testset "Plm" begin
            for norm in [LMnorm(), Unnormalized(), Orthonormal()]
                @test computePlmcostheta(0, 10, norm = norm) ≈ computePlmcostheta(NorthPole(),10, norm = norm)
            end
            for (l,m) in ML(0:4, ZeroTo)
                @test associatedLegendre(0, l, m, norm = SphericalHarmonics.Unnormalized()) ≈ (m == 0) atol=1e-14 rtol=1e-8
            end
            P1 = computePlmcostheta(0, 10, norm = SphericalHarmonics.Unnormalized())
            P2 = computePlmx(1, 10, norm = SphericalHarmonics.Unnormalized())
            for (l,m) in ML(0:10, ZeroTo)
                @test P1[(l,m)] ≈ (m == 0) atol=1e-14 rtol=1e-8
                @test P2[(l,m)] ≈ (m == 0) atol=1e-14 rtol=1e-8
            end
        end
        @testset "computePlmcostheta!" begin
            lmax = 3
            P = SphericalHarmonics.allocate_p(lmax)
            for m = 0:lmax
                computePlmcostheta!(P, NorthPole(), lmax, m)

                P2 = computePlmcostheta(NorthPole(), lmax, m)
                P3 = computePlmcostheta(0, lmax, m)

                for l = m:lmax
                    @test P[(l,m)] == P2[(l,m)]
                    @test isapprox(P[(l,m)], P3[(l,m)], atol=1e-10, rtol=1e-10)
                end
            end
        end
        @testset "computeYlm!" begin
            lmax = 3
            Y = SphericalHarmonics.allocate_y(lmax)
            P = computePlmcostheta(NorthPole(), lmax)
            computeYlm!(Y, P, NorthPole(), 0, lmax)
            Y2 = computeYlm(NorthPole(), 0, lmax)
            @test all(Y ≈ Y2)

            for m in -lmax:lmax
                computeYlm!(Y, P, NorthPole(), 0, lmax, m)
                Y2 = computeYlm(NorthPole(), 0, lmax, m)
                Y3 = computeYlm(0, 0, lmax, m)
                for l in abs(m):lmax
                    @test Y[(l,m)] == Y2[(l,m)]
                    @test isapprox(Y[(l,m)], Y3[(l,m)], atol=1e-10, rtol=1e-10)
                end
            end
        end
        @testset "Ylm" begin
            Y0 = computeYlm(0,0,10,ZeroTo)
            Y0R = computeYlm(0,0,10,ZeroTo,RealHarmonics())
            Y1 = computeYlm(0,0,10)
            @test Y1 ≈ computeYlm(NorthPole(),10)
            @test Y1 ≈ computeYlm(NorthPole(),lmax = 10)
            @test Y1 ≈ computeYlm(NorthPole(),π/2,10)
            @test Y0 ≈ computeYlm(NorthPole(),10,ZeroTo)
            @test Y0 ≈ computeYlm(NorthPole(),π/2,10,ZeroTo)
            @test Y0R ≈ computeYlm(NorthPole(),10,ZeroTo,RealHarmonics())
            @test Y0R ≈ computeYlm(NorthPole(),π/2,10,ZeroTo,RealHarmonics())
        end
        @testset "trignometric functions" begin
            @test cos(NorthPole()) == 1
            @test sec(NorthPole()) == 1
            @test sin(NorthPole()) == 0
        end
    end
    @testset "south pole" begin
        @testset "Plm" begin
            @test computePlmcostheta(π,10) ≈ computePlmcostheta(SouthPole(),10)
        end
        @testset "computePlmcostheta!" begin
            lmax = 3
            P = SphericalHarmonics.allocate_p(lmax)

            for m = 0:lmax
                computePlmcostheta!(P, SouthPole(), lmax, m)

                P2 = computePlmcostheta(SouthPole(), lmax, m)
                P3 = computePlmcostheta(π, lmax, m)
                for l = m:lmax
                    @test P[(l,m)] == P2[(l,m)]
                    @test isapprox(P[(l,m)], P3[(l,m)], atol=1e-10, rtol=1e-10)
                end
            end
        end
        @testset "computeYlm!" begin
            lmax = 3
            Y = SphericalHarmonics.allocate_y(lmax)
            P = computePlmcostheta(SouthPole(), lmax)
            computeYlm!(Y, P, SouthPole(), 0, lmax)
            Y2 = computeYlm(SouthPole(), 0, lmax)
            @test all(Y ≈ Y2)

            for m = -lmax:lmax
                computeYlm!(Y, P, SouthPole(), 0, lmax, m)
                Y2 = computeYlm(SouthPole(), 0, lmax, m)
                Y3 = computeYlm(π, 0, lmax, m)
                for l in abs(m):lmax
                    @test Y[(l,m)] == Y2[(l,m)]
                    @test isapprox(Y[(l,m)], Y2[(l,m)], atol=1e-10, rtol=1e-10)
                end
            end
        end
        @testset "Ylm" begin
            @test computeYlm(π,0,10) ≈ computeYlm(SouthPole(),10)
            @test computeYlm(π,0,10) ≈ computeYlm(SouthPole(),π/2,10)
            @test computeYlm(π,0,10,ZeroTo) ≈ computeYlm(SouthPole(),10,ZeroTo)
            @test computeYlm(π,0,10,ZeroTo) ≈ computeYlm(SouthPole(),π/2,10,ZeroTo)
            @test computeYlm(π,0,10,ZeroTo,RealHarmonics()) ≈ computeYlm(SouthPole(),10,ZeroTo,RealHarmonics())
            @test computeYlm(π,0,10,ZeroTo,RealHarmonics()) ≈ computeYlm(SouthPole(),π/2,10,ZeroTo,RealHarmonics())
        end
        @testset "trignometric functions" begin
            @test cos(SouthPole()) == -1
            @test sec(SouthPole()) == -1
            @test sin(SouthPole()) == 0
        end
    end
end

@testset "precision" begin

    isapproxdefault((x,y); kw...) = isapprox(ComplexF64(x), y; kw...)

    @testset "θ = pi/2" begin
        P = computePlmcostheta(big(pi)/2, lmax = big(10))
        @test P[(1,1)] ≈ -(√(3/big(pi)))/2
        @test P[(2,2)] ≈ (√(15/big(pi)))/4
        @test P[(3,3)] ≈ -(√(35/(2big(pi))))/4
        @test P[(4,4)] ≈ (√(35/big(pi)))*3/16
        @test P[(5,5)] ≈ -(√(77/(2big(pi))))*3/16
        @test P[(6,6)] ≈ (√(3003/(2big(pi))))/32
        @test P[(7,7)] ≈ -(√(715/big(pi)))*3/64
        @test P[(8,8)] ≈ (√(12155/big(pi)))*3/256
        @test P[(9,9)] ≈ -(√(230945/(2big(pi))))/256
        @test P[(10,10)] ≈ (√(969969/(2big(pi))))/512
    end

    θ_range_big = LinRange(0, big(pi), 20)
    ϕ_range_big = LinRange(0, 2big(pi), 5)

    @testset "lmax 1000" begin
        lmax = 1000

        coeff_big = SphericalHarmonics.compute_coefficients(big(lmax))
        coeff_Float64 = SphericalHarmonics.compute_coefficients(lmax)

        Y1 = SphericalHarmonics.allocate_y(Complex{BigFloat}, lmax)
        Y2 = SphericalHarmonics.allocate_y(Complex{Float64}, lmax)

        P1 = SphericalHarmonics.allocate_p(BigFloat, lmax)
        P2 = SphericalHarmonics.allocate_p(Float64, lmax)

        for θ in θ_range_big
            computePlmcostheta!(P1, θ, lmax, coeff_big)
            computePlmcostheta!(P2, Float64(θ), lmax, coeff_Float64)
            for ϕ in ϕ_range_big
                computeYlm!(Y1, P1, θ, ϕ, lmax)
                computeYlm!(Y2, P2, Float64(θ), Float64(ϕ), lmax)
                @test all(x -> isapproxdefault(x, atol=1e-11,  rtol=1e-11), zip(Y1, Y2))
            end
        end
    end
    @testset "lmax 100" begin
        lmax = 100

        coeff_big = SphericalHarmonics.compute_coefficients(big(lmax))
        coeff_Float64 = SphericalHarmonics.compute_coefficients(lmax)

        Y1 = SphericalHarmonics.allocate_y(Complex{BigFloat}, lmax)
        Y2 = SphericalHarmonics.allocate_y(Complex{Float64}, lmax)

        P1 = SphericalHarmonics.allocate_p(BigFloat, lmax)
        P2 = SphericalHarmonics.allocate_p(Float64, lmax)

        for θ in θ_range_big
            computePlmcostheta!(P1, θ, lmax, coeff_big)
            computePlmcostheta!(P2, Float64(θ), lmax, coeff_Float64)
            for ϕ in ϕ_range_big
                computeYlm!(Y1, P1, θ, ϕ, lmax)
                computeYlm!(Y2, P2, Float64(θ), Float64(ϕ), lmax)
                @test all(x -> isapproxdefault(x, atol=1e-13,  rtol=1e-13), zip(Y1, Y2))
            end
        end
    end

    function doublefactorial(l)
        prod(float(i) for i in big(l):-2:1)
    end

    # check against the explicit value for θ = pi/2
    for l in 499:1000
        P = SphericalHarmonics.associatedLegendrex(big(0), l = l, m = l, norm = Unnormalized())
        @test P ≈ (-1)^l * doublefactorial(2l-1)
    end
end

@testset "single mode" begin
    @testset "Associated Legendre Polynomials" begin
        lmax = 20
        coeff = SphericalHarmonics.compute_coefficients(lmax)

        @testset "explicit coeff" begin
            P = SphericalHarmonics.allocate_p(lmax)
            for θ in LinRange(0, pi, 10)
                computePlmcostheta!(P, θ, lmax, coeff)
                for l = 0:lmax, m = 0:l
                    Plm = SphericalHarmonics.associatedLegendre(θ, l, m, coeff)
                    @test isapprox(P[(l,m)], Plm, atol = 1e-14, rtol = 1e-14)
                end
            end

            P = computePlmcostheta(SphericalHarmonics.NorthPole(), lmax = lmax)
            for l = 0:lmax, m = 0:l
                Plm = SphericalHarmonics.associatedLegendre(SphericalHarmonics.NorthPole(), l, m, coeff)
                @test isapprox(P[(l,m)], Plm, atol = 1e-14, rtol = 1e-14)
            end

            P = computePlmcostheta(SphericalHarmonics.SouthPole(), lmax = lmax)
            for l = 0:lmax, m = 0:l
                Plm = SphericalHarmonics.associatedLegendre(SphericalHarmonics.SouthPole(), l, m, coeff)
                @test isapprox(P[(l,m)], Plm, atol = 1e-14, rtol = 1e-14)
            end
        end
        @testset "implicit coeff" begin
            P = SphericalHarmonics.allocate_p(lmax)
            for θ in LinRange(0, pi, 10)
                computePlmcostheta!(P, θ, lmax, coeff)
                for l = 0:lmax, m = 0:l
                    Plm = SphericalHarmonics.associatedLegendre(θ, l, m)
                    Plm2 = SphericalHarmonics.associatedLegendre(θ, l = l, m = m)
                    @test Plm2 == Plm
                    @test isapprox(P[(l,m)], Plm, atol = 1e-14, rtol = 1e-14)
                end
            end

            P = computePlmcostheta(SphericalHarmonics.NorthPole(), lmax = lmax)
            for l = 0:lmax, m = 0:l
                Plm = SphericalHarmonics.associatedLegendre(SphericalHarmonics.NorthPole(), l, m)
                @test isapprox(P[(l,m)], Plm, atol = 1e-14, rtol = 1e-14)
            end

            P = computePlmcostheta(SphericalHarmonics.SouthPole(), lmax = lmax)
            for l = 0:lmax, m = 0:l
                Plm = SphericalHarmonics.associatedLegendre(SphericalHarmonics.SouthPole(), l, m)
                @test isapprox(P[(l,m)], Plm, atol = 1e-14, rtol = 1e-14)
            end
        end
    end
    @testset "SphericalHarmonics" begin
        lmax = 20
        coeff = SphericalHarmonics.compute_coefficients(lmax);

        P = SphericalHarmonics.allocate_p(lmax)
        Y = SphericalHarmonics.allocate_y(lmax)
        Yreal = SphericalHarmonics.allocate_y(Float64, lmax)
        @testset "ComplexHarmonics" begin
            @testset "explicit coeff" begin
                for θ in LinRange(0, pi, 10)
                    computePlmcostheta!(P, θ, lmax, coeff)
                    for ϕ in LinRange(0, 2pi, 10)[1:end-1]
                        computeYlm!(Y, P, θ, ϕ, lmax = lmax)
                        for l = 0:lmax, m = -l:l
                            Ylm = SphericalHarmonics.sphericalharmonic(θ, ϕ, l, m,
                                SphericalHarmonics.ComplexHarmonics(), coeff)

                            Ylm2 = SphericalHarmonics.sphericalharmonic(θ, ϕ; l = l, m = m,
                                SHType = SphericalHarmonics.ComplexHarmonics(), coeff = coeff)

                            @test Ylm == Ylm2

                            @test isapprox(Y[(l,m)], Ylm, atol = 1e-14, rtol = 1e-14)
                        end
                    end
                end
            end
            @testset "implicit coeff" begin
                for θ in LinRange(0, pi, 10)
                    computePlmcostheta!(P, θ, lmax, coeff)
                    for ϕ in LinRange(0, 2pi, 10)[1:end-1]
                        computeYlm!(Y, P, θ, ϕ, lmax = lmax)
                        for l = 0:lmax, m = -l:l
                            Ylm = SphericalHarmonics.sphericalharmonic(θ, ϕ, l, m,
                                SphericalHarmonics.ComplexHarmonics())

                            Ylm2 = SphericalHarmonics.sphericalharmonic(θ, ϕ; l = l, m = m,
                                SHType = SphericalHarmonics.ComplexHarmonics())

                            @test Ylm == Ylm2

                            @test isapprox(Y[(l,m)], Ylm, atol = 1e-14, rtol = 1e-14)
                        end
                    end
                end
            end
        end
        @testset "RealHarmonics" begin
            @testset "explicit coeff" begin
                for θ in LinRange(0, pi, 10)
                    computePlmcostheta!(P, θ, lmax, coeff)
                    for ϕ in LinRange(0, 2pi, 10)[1:end-1]
                        computeYlm!(Yreal, P, θ, ϕ, lmax = lmax, SHType = SphericalHarmonics.RealHarmonics())
                        for l = 0:lmax, m = -l:l
                            Ylm = SphericalHarmonics.sphericalharmonic(θ, ϕ, l, m,
                                SphericalHarmonics.RealHarmonics(), coeff)

                            Ylm2 = SphericalHarmonics.sphericalharmonic(θ, ϕ; l = l, m = m,
                                SHType = SphericalHarmonics.RealHarmonics(), coeff = coeff)

                            @test Ylm == Ylm2

                            @test isapprox(Yreal[(l,m)], Ylm, atol = 1e-14, rtol = 1e-8)
                        end
                    end
                end
            end
            @testset "implicit coeff" begin
                for θ in LinRange(0, pi, 10)
                    computePlmcostheta!(P, θ, lmax, coeff)
                    for ϕ in LinRange(0, 2pi, 10)[1:end-1]
                        computeYlm!(Yreal, P, θ, ϕ, lmax = lmax, SHType = SphericalHarmonics.RealHarmonics())
                        for l = 0:lmax, m = -l:l
                            Ylm = SphericalHarmonics.sphericalharmonic(θ, ϕ, l, m,
                                SphericalHarmonics.RealHarmonics())

                            Ylm2 = SphericalHarmonics.sphericalharmonic(θ, ϕ; l = l, m = m,
                                SHType = SphericalHarmonics.RealHarmonics())

                            @test Ylm == Ylm2

                            @test isapprox(Yreal[(l,m)], Ylm, atol = 1e-14, rtol = 1e-8)
                        end
                    end
                end
            end
        end
    end
end

@testset "Orthonormality" begin

    @testset "Associated Legendre Polynomials" begin
        shnorm(::SphericalHarmonics.LMnorm) = 1/pi
        shnorm(::SphericalHarmonics.Orthonormal) = 1
        function testnorm(f, norm = SphericalHarmonics.LMnorm())
            I, E = hcubature(f, [0], [π]);
            @test I ≈ shnorm(norm)
        end
        function testorthogonal(f)
            I, E = hcubature(f, [0], [π], atol=1e-10);
            @test isapprox(abs(I), 0, atol = max(abs(E), 1e-10))
        end
        lmax = 10
        coeff = SphericalHarmonics.compute_coefficients(lmax);
        P = SphericalHarmonics.allocate_p(lmax);

        @testset "Normalization" begin
            function f!(x, P, l, m, norm = SphericalHarmonics.LMnorm())
                θ = x[1]
                computePlmcostheta!(P, θ, l, m, coeff, norm = norm)
                sin(θ) * P[(l,m)]^2
            end

            for norm in [SphericalHarmonics.LMnorm(), SphericalHarmonics.Orthonormal()]
                for m = 0:lmax, l = m:lmax
                    testnorm(x -> f!(x, P, l, m, norm), norm)
                end
            end
        end

        @testset "Orthogonality" begin
            @testset "same m, different l" begin
                function f!(x, P, k, l, m, norm = SphericalHarmonics.LMnorm())
                    θ = x[1]
                    computePlmcostheta!(P, θ, max(k, l), m, coeff, norm = norm)
                    sin(θ) * P[(k,m)] * P[(l,m)]
                end

                for norm in [SphericalHarmonics.LMnorm(), SphericalHarmonics.Orthonormal()]
                    for m = 0:lmax, l = m:lmax, k = m:lmax
                        k == l && continue
                        testorthogonal(x -> f!(x, P, k, l, m, norm))
                    end
                end
            end

            @testset "same l, different m" begin
                function f!(x, P, l, m1, m2)
                    θ = x[1]
                    computePlmcostheta!(P, θ, l, coeff)
                    P[(l,m1)] * P[(l,m2)] / sin(θ)
                end

                for l = 0:lmax, m1 in 0:l, m2 in m1+1:l
                    testorthogonal(x -> f!(x, P, l, m1, m2))
                end
            end
        end
    end

    @testset "SphericalHarmonics" begin
        function testnorm(f)
            I, E = hcubature(f, [0, 0], [π, 2π]);
            @test all(x->isapprox(x, 1), I)
        end
        function testortho(f)
            I, E = hcubature(f, [0, 0], [π, 2π], atol=1e-10);
            @test isapprox(abs(I), 0, atol = max(abs(E), 1e-10))
        end

        @testset "Normalization" begin
            lmax = 10
            @testset "Complex harmonics" begin
                f(θϕ) = sin(θϕ[1]) .* abs2.(computeYlm(θϕ[1], θϕ[2], lmax = lmax))
                testnorm(f)
            end
            @testset "Real harmonics" begin
                f(θϕ) = sin(θϕ[1]) .* abs2.(computeYlm(θϕ[1], θϕ[2], lmax = lmax, SHType = SphericalHarmonics.RealHarmonics()))
                testnorm(f)
            end
        end
        @testset "Orthogonality" begin
            lmax = 3
            coeff = SphericalHarmonics.compute_coefficients(lmax)

            @testset "Complex harmonics" begin
                Y = SphericalHarmonics.allocate_y(lmax)
                P = SphericalHarmonics.allocate_p(lmax)

                function f!(θϕ, Y, P, l, m1, m2)
                    θ, ϕ = θϕ
                    computePlmcostheta!(P, θ, l, coeff)
                    computeYlm!(Y, P, θ, ϕ, l)
                    sin(θ) * Y[(l,m1)] * conj(Y[(l,m2)])
                end

                for l in 0:lmax, m1 in -l:0, m2 in m1+1:l
                    testortho(x->f!(x, Y, P, l, m1, m2))
                end
            end

            @testset "Real harmonics" begin
                Y = SphericalHarmonics.allocate_y(Float64, lmax)
                P = SphericalHarmonics.allocate_p(lmax)

                function f!(θϕ, Y, P, l, m1, m2)
                    θ, ϕ = θϕ
                    computePlmcostheta!(P, θ, l, coeff)
                    computeYlm!(Y, P, θ, ϕ, l, SphericalHarmonicModes.FullRange, SphericalHarmonics.RealHarmonics())
                    sin(θ) * Y[(l,m1)] * Y[(l,m2)]
                end

                for l in 0:lmax, m1 in -l:0, m2 in m1+1:l
                    testortho(x->f!(x, Y, P, l, m1, m2))
                end
            end
        end
    end
end

@testset "cache" begin
    @testset "constructors" begin
        S1 = @inferred SphericalHarmonics.cache(3);
        S2 = @inferred SphericalHarmonics.cache(3, FullRange);
        S3 = @inferred SphericalHarmonics.cache(3, FullRange, ComplexHarmonics());
        S4 = @inferred SphericalHarmonics.cache(Float64, 3, FullRange, ComplexHarmonics());
        S5 = @inferred SphericalHarmonics.cache(3, ComplexHarmonics());
        S6 = @inferred SphericalHarmonics.cache(Float64, 3, ComplexHarmonics());
        @test typeof(S1) == typeof(S2) == typeof(S3) == typeof(S4) == typeof(S5) == typeof(S6)
    end
    @testset "FullRange, ComplexHarmonics" begin
        S = SphericalHarmonics.cache(3);
        @test _lmax(S) == 3
        θ, ϕ = pi/3, pi/4
        SphericalHarmonics.computePlmx!(S, cos(θ), 1)
        @test S.P.lmax == 1
        P2 = SphericalHarmonics.computePlmx(cos(θ), lmax = 1)
        @test S.P[1:length(P2)] == P2
        SphericalHarmonics.computePlmx!(S, cos(2θ))
        @test S.P.lmax == _lmax(S)
        @test S.P == SphericalHarmonics.computePlmx(cos(2θ), lmax = 3)
        computePlmcostheta!(S, θ, 3)
        @test _lmax(S) == 3
        @test S.P == computePlmcostheta(θ, lmax = 3)
        # This should be a no-op
        computePlmcostheta!(S, θ, 3)
        @test S.P == computePlmcostheta(θ, lmax = 3)
        computePlmcostheta!(S, 2θ)
        @test S.P ≈ computePlmcostheta(2θ, lmax = 3)

        computePlmcostheta!(S, θ)
        computeYlm!(S, θ, ϕ, 3)
        @test S.Y == computeYlm(θ, ϕ, 3)
        computePlmcostheta!(S, θ, 4)
        @test _lmax(S) == 4
        @test S.P == computePlmcostheta(θ, lmax = 4)
        computeYlm!(S, θ, ϕ, 4)
        @test S.Y == computeYlm(θ, ϕ, 4)
        computeYlm!(S, θ, ϕ)
        @test S.Y == computeYlm(θ, ϕ, 4)
        computeYlm!(S, θ, ϕ, 2)
        Y2 = computeYlm(θ, ϕ, 2)
        @test S.Y[1:length(Y2)] == Y2

        computePlmcostheta!(S, NorthPole(), 4)
        @test S.P.cosθ == cos(NorthPole())
        @test S.P == computePlmcostheta(NorthPole(), 4)
    end
    @testset "ZeroTo, ComplexHarmonics" begin
        S = SphericalHarmonics.cache(3, m_range = ZeroTo);
        @test _lmax(S) == 3
        θ, ϕ = pi/3, pi/4
        SphericalHarmonics.computePlmx!(S, cos(θ), 3)
        @test S.P == SphericalHarmonics.computePlmx(cos(θ), lmax = 3)
        computePlmcostheta!(S, θ, 3)
        @test S.P == computePlmcostheta(θ, lmax = 3)
        # This should be a no-op
        computePlmcostheta!(S, θ, 3)
        @test S.P == computePlmcostheta(θ, lmax = 3)
        computeYlm!(S, θ, ϕ, 3)
        @test S.Y == computeYlm(θ, ϕ, lmax = 3, m_range = ZeroTo)
        computePlmcostheta!(S, θ, 4)
        @test _lmax(S) == 4
        @test S.P == computePlmcostheta(θ, lmax = 4)
        computeYlm!(S, θ, ϕ, 4)
        @test S.Y == computeYlm(θ, ϕ, lmax = 4, m_range = ZeroTo)
        computeYlm!(S, θ, ϕ)
        @test S.Y == computeYlm(θ, ϕ, lmax = 4, m_range = ZeroTo)
        computeYlm!(S, θ, ϕ, 2)
        Y2 = computeYlm(θ, ϕ, lmax = 2, m_range = ZeroTo)
        @test S.Y[1:length(Y2)] == Y2

        computePlmcostheta!(S, NorthPole(), 4)
        @test S.P.cosθ == cos(NorthPole())
        @test S.P == computePlmcostheta(NorthPole(), 4)
    end
    @testset "FullRange, RealHarmonics" begin
        S = SphericalHarmonics.cache(3, SHType = RealHarmonics());
        @test _lmax(S) == 3
        θ, ϕ = pi/3, pi/4
        SphericalHarmonics.computePlmx!(S, cos(θ), 3)
        @test S.P == SphericalHarmonics.computePlmx(cos(θ), lmax = 3)
        # This should be a no-op
        SphericalHarmonics.computePlmx!(S, cos(θ), 3)
        @test S.P == SphericalHarmonics.computePlmx(cos(θ), lmax = 3)
        computePlmx!(S, cos(θ), 4)
        @test _lmax(S) == 4
        @test S.P == computePlmx(cos(θ), lmax = 4)

        S = SphericalHarmonics.cache(3, SHType = RealHarmonics());
        computePlmcostheta!(S, θ, 3)
        @test S.P == computePlmcostheta(θ, lmax = 3)
        # This should be a no-op
        computePlmcostheta!(S, θ, 3)
        @test S.P == computePlmcostheta(θ, lmax = 3)
        computeYlm!(S, θ, ϕ, 3)
        @test S.Y == computeYlm(θ, ϕ, lmax = 3, SHType = RealHarmonics())

        computePlmcostheta!(S, θ, 4)
        @test _lmax(S) == 4
        @test S.P == computePlmcostheta(θ, lmax = 4)
        computeYlm!(S, θ, ϕ, 4)
        @test S.Y == computeYlm(θ, ϕ, lmax = 4, SHType = RealHarmonics())
        computeYlm!(S, θ, ϕ)
        @test S.Y == computeYlm(θ, ϕ, lmax = 4, SHType = RealHarmonics())
        computeYlm!(S, θ, ϕ, 2)
        Y2 = computeYlm(θ, ϕ, lmax = 2, SHType = RealHarmonics())
        @test S.Y[1:length(Y2)] == Y2

        computePlmcostheta!(S, NorthPole(), 4)
        @test S.P.cosθ == cos(NorthPole())
        @test S.P == computePlmcostheta(NorthPole(), 4)
    end
    @testset "ZeroTo, RealHarmonics" begin
        S = SphericalHarmonics.cache(3, m_range = ZeroTo, SHType = RealHarmonics());
        @test _lmax(S) == 3
        θ, ϕ = pi/3, pi/4
        SphericalHarmonics.computePlmx!(S, cos(θ), 3)
        @test S.P == SphericalHarmonics.computePlmx(cos(θ), lmax = 3)
        computePlmcostheta!(S, θ, 3)
        @test S.P == computePlmcostheta(θ, lmax = 3)
        # This should be a no-op
        computePlmcostheta!(S, θ, 3)
        @test S.P == computePlmcostheta(θ, lmax = 3)
        computeYlm!(S, θ, ϕ, 3)
        @test S.Y == computeYlm(θ, ϕ, lmax = 3, m_range = ZeroTo, SHType = RealHarmonics())
        computePlmcostheta!(S, θ, 4)
        @test _lmax(S) == 4
        @test S.P == computePlmcostheta(θ, lmax = 4)
        computeYlm!(S, θ, ϕ, 4)
        @test S.Y == computeYlm(θ, ϕ, lmax = 4, m_range = ZeroTo, SHType = RealHarmonics())
        computeYlm!(S, θ, ϕ)
        @test S.Y == computeYlm(θ, ϕ, lmax = 4, m_range = ZeroTo, SHType = RealHarmonics())
        computeYlm!(S, θ, ϕ, 2)
        Y2 = computeYlm(θ, ϕ, lmax = 2, m_range = ZeroTo, SHType = RealHarmonics())
        @test S.Y[1:length(Y2)] == Y2

        computePlmcostheta!(S, NorthPole(), 4)
        @test S.P.cosθ == cos(NorthPole())
        @test S.P == computePlmcostheta(NorthPole(), 4)
    end
    @testset "show" begin
        S = SphericalHarmonics.cache(3);
        io = IOBuffer()
        show(io, S)
        s = String(take!(io))
        s_exp = "$(SphericalHarmonics.SphericalHarmonicsCache)(Float64, 3, $(SphericalHarmonicModes.FullRange), $(SphericalHarmonics.ComplexHarmonics)())"
        @test s == s_exp
        summary(io, S.P)
        s = String(take!(io))
        s_exp = "10-element normalized AssociatedLegendrePolynomials{Float64} for lmax = 3 (uninitialized)"
        @test s == s_exp
        θ = 0
        computePlmcostheta!(S, θ, 3)
        summary(io, S.P)
        s = String(take!(io))
        s_exp = "10-element normalized AssociatedLegendrePolynomials{Float64} for lmax = 3 and cosθ = 1"
        @test s == s_exp

        S = SphericalHarmonics.cache(BigFloat, 3);
        summary(io, S.P)
        s = String(take!(io))
        s_exp = "10-element normalized AssociatedLegendrePolynomials{BigFloat} for lmax = 3 (uninitialized)"
    end
    @testset "accessor methods" begin
        S = SphericalHarmonics.cache(Float64, 3);
        @test @inferred SphericalHarmonics.eltypeP(S) == Float64
        @test @inferred SphericalHarmonics.eltypeY(S) == ComplexF64
        @test SphericalHarmonics.getP(S) === S.P
        @inferred SphericalHarmonics.getP(S)
        @test SphericalHarmonics.getY(S) === S.Y
        @inferred SphericalHarmonics.getY(S)
    end
    @testset "deepcopy" begin
        S = SphericalHarmonics.cache(3);
        θ, ϕ = pi/3, pi/4
        SphericalHarmonics.computePlmcostheta!(S, θ)
        SphericalHarmonics.computeYlm!(S, θ, ϕ)
        Sd = deepcopy(S)
        @test Sd.C == S.C
        @test Sd.P == S.P
        @test Sd.Y == S.Y
        @test Sd.SHType == S.SHType
    end
end

@testset "Parity" begin
    lmax = 10
    S1 = SphericalHarmonics.cache(lmax)
    S2 = SphericalHarmonics.cache(lmax)
    for θ in LinRange(0, pi, 10)
        P1 = computePlmcostheta!(S1, θ)
        P2 = computePlmcostheta!(S2, pi - θ)
        for l in 0:lmax, m in 0:l
            @test isapprox(P2[(l,m)], (-1)^(l + m) * P1[(l,m)], atol = 1e-13, rtol = 1e-13)
        end
        for ϕ in LinRange(0, 2pi, 10)
            Y1 = computeYlm!(S1, θ, ϕ)
            Y2 = computeYlm!(S2, pi - θ, pi + ϕ)
            for l in 0:lmax, m in -l:l
                @test isapprox(Y2[(l,m)], (-1)^l * Y1[(l,m)], atol = 1e-13, rtol = 1e-13)
            end
        end
    end
end

@testset "m = 0 Legendre" begin
    # for m = 0, the spherical harmonics are normalized Legendre Polynomials and are purely real
    lmax = 100
    @testset "ComplexHarmonics" begin
        S = SphericalHarmonics.cache(lmax)
        LP = zeros(0:lmax)
        for θ in LinRange(0, pi, 10)
            LegendrePolynomials.collectPl!(LP, cos(θ), lmax = lmax)
            computePlmcostheta!(S, θ)
            for ϕ in LinRange(0, 2pi, 10)
                computeYlm!(S, θ, ϕ)
                for l in 0:lmax
                    @test isapprox(imag(S.Y[(l,0)]), 0, atol = 1e-13)
                    @test isapprox(S.Y[(l,0)], √((2l+1)/4pi) * LP[l], atol = 1e-13, rtol = 1e-13)
                end
            end
        end
    end
    @testset "RealHarmonics" begin
        S = SphericalHarmonics.cache(lmax, SHType = RealHarmonics())
        LP = zeros(0:lmax)
        for θ in LinRange(0, pi, 10)
            LegendrePolynomials.collectPl!(LP, cos(θ), lmax = lmax)
            computePlmcostheta!(S, θ)
            for ϕ in LinRange(0, 2pi, 10)
                computeYlm!(S, θ, ϕ)
                for l in 0:lmax
                    @test isapprox(imag(S.Y[(l,0)]), 0, atol = 1e-13)
                    @test isapprox(S.Y[(l,0)], √((2l+1)/4pi) * LP[l], atol = 1e-13, rtol = 1e-13)
                end
            end
        end
    end
end

@testset "rotation of coordinates" begin
    # Ylm(θ, ϕ) = Yl0(0, 0) D^l_{0,-m}(0, θ, ϕ)
    lmax = 10
    S = SphericalHarmonics.cache(lmax);
    Dvec = OffsetArray([zeros(ComplexF64, 2l+1, 2l+1) for l in 0:lmax], 0:lmax);
    Jyvec = OffsetArray([zeros(ComplexF64, 2l+1, 2l+1) for l in 0:lmax], 0:lmax);
    @testset "NorthPole" begin
        Yl0NP = computeYlm(NorthPole(), 0, lmax = lmax)
        for θ in LinRange(0, pi, 10)
            computePlmcostheta!(S, θ)
            for ϕ in LinRange(0, 2pi, 10)
                Ylmθϕ = computeYlm!(S, θ, ϕ)
                for l in 0:lmax
                    Dp = wignerD!(Dvec[l], l, 0, θ, ϕ, Jyvec[l])
                    D = OffsetArray(Dp, -l:l, -l:l)
                    for m in -l:l
                        @test isapprox(Ylmθϕ[(l,m)], Yl0NP[(l,0)] * D[0,-m], atol = 1e-13, rtol = 1e-8)
                    end
                end
            end
        end
    end
    @testset "arbitrary point" begin
        # We represent an arbitrary point in two coordinate systems
        # We assume that it has coordinates (θ1,ϕ1) in the first frame and (θ2,ϕ2) in the second frame
        # The frames may be shown to be related through S2 = Rz(ϕ1)Ry(θ1-θ2)Rz(-ϕ2) S1,
        # where Ri(ω) represents a rotation about the axis i by an angle ω.
        # We note that the rotation operator may be represented as D = exp(-iαJz)exp(-iβJy)exp(-iγJz)
        # Comparing, we obtain α = ϕ1, β = θ1-θ2, γ = -ϕ2
        # The Spherical harmonics transform through Ylm(n2) = ∑_m′ conj(D^l_{m,m′}(-γ,-β,-α)) Ylm(n1)
        # Substituting, we obtain the matrix conj(D^l_{m,m′}(ϕ2,θ2-θ1,-ϕ1))
        # We note that choosing (θ1,ϕ1) = (0,0) in this gives us the matrix conj(D^l_{m,m′}(ϕ2,θ2,0)),
        # which, using symmetries, reduces to the previous case with the point at the north pole

        θ1, ϕ1 = pi/6, pi/3

        Y1lm = computeYlm(θ1, ϕ1, lmax = lmax)
        for θ2 in LinRange(0, pi, 10)
            computePlmcostheta!(S, θ2)
            for ϕ2 in LinRange(0, 2pi, 10)
                Ylmθϕ = computeYlm!(S, θ2, ϕ2)
                for l in 0:lmax
                    Dp = wignerD!(Dvec[l], l, ϕ2, θ2-θ1, -ϕ1, Jyvec[l])
                    D = OffsetArray(Dp, -l:l, -l:l)
                    for m in -l:l
                        rotY1lm = sum(conj(D[m, m′]) * Y1lm[(l,m′)] for m′ = -l:l)
                        @test isapprox(Ylmθϕ[(l,m)], rotY1lm, atol = 1e-13, rtol = 1e-8)
                    end
                end
            end
        end
    end
end

@testset "negative theta" begin
    lmax = 10
    S1 = SphericalHarmonics.cache(lmax);
    S2 = SphericalHarmonics.cache(lmax);
    for θ in LinRange(0, pi, 10)
        computePlmcostheta!(S1, θ)
        θ2 = -θ
        computePlmcostheta!(S2, θ2)
        for ϕ in LinRange(0, 2pi, 10)
            Ylmθϕ = computeYlm!(S1, θ, ϕ)
            Ylmθ2ϕ = computeYlm!(S2, θ2, ϕ)
            for l in 0:lmax, m in -l:l
                @test isapprox(Ylmθ2ϕ[(l,m)], (-1)^m * Ylmθϕ[(l,m)], atol = 1e-13, rtol = 1e-8)
            end
        end
    end
end

@testset "theta ± nπ" begin
    lmax = 10
    S1 = SphericalHarmonics.cache(lmax);
    S2 = SphericalHarmonics.cache(lmax);
    S3 = SphericalHarmonics.cache(lmax);
    for θ in LinRange(0, pi, 10)
        computePlmcostheta!(S1, θ)
        θ2 = pi + θ
        computePlmcostheta!(S2, θ2)
        θ3 = 2pi + θ
        computePlmcostheta!(S3, θ3)
        for ϕ in LinRange(0, 2pi, 10)
            Ylmθϕ = computeYlm!(S1, θ, ϕ)
            Ylmθ2ϕ = computeYlm!(S2, θ2, ϕ)
            Ylmθ3ϕ = computeYlm!(S3, θ3, ϕ)
            for l in 0:lmax, m in -l:l
                @test isapprox(Ylmθ2ϕ[(l,m)], (-1)^l * Ylmθϕ[(l,m)], atol = 1e-13, rtol = 1e-8)
                @test isapprox(Ylmθ3ϕ[(l,m)], Ylmθϕ[(l,m)], atol = 1e-13, rtol = 1e-8)
            end
        end
    end
end

@testset "derivatives of Plm" begin
    for l in 1:20, m in l:-2:0
        # the polynomials for even l+m are even, so their derivatives are odd
        # we check that these go to zero exactly at θ = pi/2
        dP = ForwardDiff.derivative(θ -> associatedLegendre(θ, l = l, m = m), Pi/2)
        @test dP == 0
        dP = ForwardDiff.derivative(x -> associatedLegendrex(x, l = l, m = m), 0)
        @test dP == 0
    end
end


text:
const TEST_RAW_S = [
    1.16e3, 1.31e3, 1.57e3, 1.76e3, 1.94e3,
    2.31e3, 2.57e3, 3.11e3, 3.38e3, 3.77e3,
    4.58e3, 5.06e3, 6.22e3, 7.51e3, 8.78e3,
    9.68e3, 1.10e4, 1.20e4, 1.27e4, 1.40e4,
    1.57e4, 1.62e4, 1.69e4, 1.77e4, 1.86e4,
    1.98e4, 2.12e4, 2.26e4, 2.36e4, 2.55e4,
    2.65e4, 2.82e4, 2.97e4, 3.12e4, 3.20e4,
    3.37e4, 3.41e4, 3.50e4, 3.64e4, 3.71e4,
    3.78e4, 3.83e4, 3.87e4
]
const TEST_S = [
    1.49e2, 2.99e2, 5.65e2, 7.57e2, 9.29e2,
    1.30e3, 1.56e3, 2.10e3, 2.37e3, 2.76e3,
    3.57e3, 4.05e3, 5.21e3, 6.51e3, 7.77e3,
    8.67e3, 9.95e3, 1.10e4, 1.17e4, 1.30e4,
    1.47e4, 1.52e4, 1.59e4, 1.67e4, 1.76e4,
    1.88e4, 2.02e4, 2.16e4, 2.26e4, 2.45e4,
    2.55e4, 2.72e4, 2.87e4, 3.02e4, 3.10e4,
    3.26e4, 3.31e4, 3.40e4, 3.53e4, 3.61e4,
    3.68e4, 3.73e4, 3.77e4
]
const TEST_RMS = [
    1.22e1, 1.63e1, 2.19e1, 2.52e1, 2.79e1,
    3.32e1, 3.65e1, 4.27e1, 4.57e1, 4.99e1,
    5.71e1, 6.14e1, 6.99e1, 7.84e1, 8.67e1,
    9.20e1, 9.88e1, 1.03e2, 1.08e2, 1.15e2,
    1.22e2, 1.26e2, 1.27e2, 1.31e2, 1.36e2,
    1.39e2, 1.44e2, 1.51e2, 1.55e2, 1.54e2,
    1.58e2, 1.63e2, 1.64e2, 1.64e2, 1.66e2,
    1.60e2, 1.59e2, 1.61e2, 1.51e2, 1.49e2,
    1.38e2, 1.30e2, 1.11e2
]

"""
    testdata()

Loads a tuple of the signal and total noise from the Appendix A table A.1 from Janesick (2007).[^1] This data is generated by a CMOS imager with V/e- nonlinearity. Data units are in DN.

# Examples

```jldoctest
julia> data = testdata();

julia> size(data.S)
(43,)
```

[^1]: James R. Janesick, 2007, ["Photon Transfer"](https://www.spiedigitallibrary.org/ebooks/PM/Photon-Transfer/eISBN-9780819478382/10.1117/3.725073?SSO=1), SPIE
"""
testdata() = (S=TEST_S, total_noise=TEST_RMS, raw_S=TEST_RAW_S)


text:
include("analysis_functions.jl")
include("plot_fcns.jl")

for dos_file in string.("DOSes/",readdir("./DOSes"))
    material = dos_file[7:end-7]
    println(material)
    eq_dos, average_eq_dos, min_E_eq, max_E_eq = get_dos(dos_file)

    for Ef in [0.0, 1.0, -1.0, -2.0]
        dos_f, average_dos, min_E, max_E = get_dos(dos_file; Ef=Ef)
        plot_comparison(dos_f, 5, min_E, max_E, average_dos, plot_title=string("E_f = ", Ef, " eV"))
        savefig(string("figs/rate_compare/", material,"_Ef",Ef,".png"))
    end

    # and plot the DOS too
    E_range = range(min_E_eq, max_E_eq-0.01, length=200)
    plot(E_range, eq_dos.(E_range), label="DOS", xlabel="energy", ylabel="DOS")
    savefig(string("figs/DOS/", material, ".png"))
end


text:
# DiscontinuousGalerkinScheme
# Based on:
# Hesthaven, Warburton, Nodal Discontinuous Galerkin Methods Algorithms, Analysis and applications

mutable struct DiscontinuousGalerkinScheme{T,vType} <: AbstractFEAlgorithm
  basis::PolynomialBasis
  riemann_solver::Function
  max_w_speed :: Function
  D::AbstractArray{T,2}
  Ma::AbstractArray{T,2}    #Mass Matrix on reference element
  S::AbstractArray{T,2}     #Stiffness matrix on reference element
  ν::vType                  #Max wave speed
end

"DiscontinuousGalerkinScheme constructor"
function DiscontinuousGalerkinScheme(basis::PolynomialBasis, riemann_solver;max_w_speed = nothing)
    if max_w_speed == nothing
        max_w_speed = maxfluxρ
    end
    D = basis.dφ*basis.invφ
    Ma = inv(basis.φ*basis.φ')
  DiscontinuousGalerkinScheme(basis, riemann_solver, max_w_speed, D, Ma, Ma*D, 0.0)
end

"Reconstruc solution from basis space"
function reconstruct_u(u::AbstractArray{T,2}, φ::AbstractArray{T2,2}, NC::Int) where {T, T2}
  uh = myblock(φ,NC)*u
  NN = size(φ,1); Nx = size(u,2)
  uₕ = zeros(T, NN*Nx,NC)
  for j in 1:NC
    uₕ[:,j] = uh[(j-1)*NN+1:j*NN,:][:]
  end
  return uₕ
end

"Flat x nodes on u matrix"
function flat_u(u::AbstractArray{T,2}, order::Int, NC::Int) where {T}
  NN = order + 1
  uh = zeros(T, size(u,2)*NN, NC)
  for j in 1:NC
    uh[:,j] = u[(j-1)*NN+1:j*NN,:][:]
  end
  return uh
end

"Get cell face values for each variable"
function get_dg_face_values(u::AbstractArray{T,2}, order::Int, NC::Int) where {T}
    NN = order + 1
    #TODO: number of faces depend on dimension
    us = zeros(T, 2*NC, size(u,2))
    for j in 1:size(u,2)
      for k in 1:NC
          us[(2*k-1):k*2,j] = u[[1+NN*(k-1),NN+NN*(k-1)],j]
      end
    end
    return us
end

"Update dt based on CFL condition"
function update_dt(alg::DiscontinuousGalerkinScheme,u::AbstractArray{T2,2},Flux,
    CFL,mesh::Uniform1DFVMesh) where {T2}
    alg.ν = alg.max_w_speed(u, Flux)
    dx = maximum(cell_volumes(mesh))
    dx * CFL*min(abs(alg.basis.nodes[1]-alg.basis.nodes[2])) / alg.ν
end

"Compute right hand side for time integration"
function residual!(H::AbstractArray{T,2}, u::AbstractArray{T,2}, basis::PolynomialBasis, mesh::Uniform1DFVMesh, alg::DiscontinuousGalerkinScheme, f, riemann_solver, NC, ::Type{Val{false}}) where {T}
    NN = basis.order + 1

    us = get_dg_face_values(u, basis.order, NC)
    #Apply boundary conditions TODO: Other boundary types
    us = apply_boundary(us, mesh)
    q = zeros(u)
    F = zeros(u)
    ur=us[1:2:end,:]
    ul=us[2:2:end,:]
    for i = 1:numcells(mesh)
        # Evaluate edge fluxes
        q[NN:NN:end,i] = riemann_solver(ul[:,i+1],ur[:,i+2], f, alg.ν)
        q[1:NN:end,i] = -riemann_solver(ul[:,i],ur[:,i+1], f, alg.ν)
        # Integrate interior fluxes ∫f(uₕ)φ'(ξ)dξ
        for k = 1:NN
            F[k:NN:end,i] = f(u[k:NN:end,i])
        end
    end

    # Compute right hand size
    ru = myblock(alg.S',NC)*F - q
    h = maximum(cell_volumes(mesh))
    H[:,:] = (h/2*myblock(alg.Ma,NC))\ru;
    if isleftdirichlet(mesh); H[:,1] = 0.0; end
   if isrightdirichlet(mesh); H[:,end] = 0.0; end
end

"Apply boundary conditions on scalar problems"
function apply_boundary(u::AbstractArray{T,2}, mesh::AbstractFVMesh1D) where {T}
  #Add ghost cells
  uh = hcat(u[:,1],u,u[:,end])
  if isleftperiodic(mesh)
      uh[:,1] = uh[:,end-1]
  elseif  isleftzeroflux(mesh)
      uh[:,1] = uh[:,2]
  end
  if isrightperiodic(mesh)
      uh[:,end] = uh[:,2]
  elseif  isrightzeroflux(mesh)
      uh[:,end] = uh[:,end-1]
  end
  return uh
end

@inline function fluxρ(uj::Vector,f)
  maximum(abs,eigvals(f(Val{:jac}, uj)))
end

"build a block diagonal matrix by repeating a matrix N times"
function myblock(A::AbstractArray{T,2},N::Int) where {T}
  M = size(A,1)
  Q = size(A,2)
  B = zeros(T,M*N,Q*N)
  for i = 1:N
    B[(i-1)*M+1:i*M,(i-1)*Q+1:i*Q] = A
  end
  B
end

"Solve scalar 1D conservation laws problems with DG Scheme"
function solve(
  prob::AbstractConservationLawProblem,
  alg::AbstractFEAlgorithm;
  TimeIntegrator::OrdinaryDiffEqAlgorithm = SSPRK22(),use_threads = false, kwargs...)

  # Unpack some useful variables
  @unpack basis, riemann_solver = alg
  #Unroll some important constants
  @unpack tspan,f,f0, mesh = prob

  N = numcells(mesh)
  NC = prob.numvars
  NN = basis.order+1
  #Assign Initial values
  xg = zeros(NN,N)
  for i in 1:N
    b = cell_faces(mesh)[i+1]; a=cell_faces(mesh)[i]
    xg[:,i] = reference_to_interval(basis.nodes,(a,b))
  end

  u0 = zeros(eltype(f0(cell_faces(mesh)[1])),NN*NC, N)
  for i = 1:N
      for k = 1:NN
        u0[k:NN:NN*NC,i] = f0(xg[k,i])
      end
  end
  #Time loop
  #First dt
  u0ₕ = flat_u(u0, basis.order, NC)
  dt = update_dt(alg, u0ₕ, f, prob.CFL, mesh)
  # Setup time integrator
  semidiscretef(du,u,p,t) = residual!(du, u, basis, mesh, alg, f, riemann_solver, NC, Val{use_threads})
  ode_prob = ODEProblem(semidiscretef, u0, prob.tspan)
  timeIntegrator = init(ode_prob, TimeIntegrator;dt=dt, kwargs...)
  @inbounds for i in timeIntegrator
    uₕ = flat_u(timeIntegrator.u, basis.order,NC)
    dt = update_dt(alg, uₕ, f, prob.CFL, mesh)
    set_proposed_dt!(timeIntegrator, dt)
  end
  if timeIntegrator.sol.t[end] != prob.tspan[end]
    savevalues!(timeIntegrator)
  end
  return build_solution(timeIntegrator.sol,xg,basis,prob, NC)
end

max_stars_repo_name: Paulms/DiscontinuousGalerkin

max_issues_repo_name: Paulms/DiscontinuousGalerkin

max_forks_repo_name: Paulms/DiscontinuousGalerkin


text:
@testset "HW1" begin
    using DMUStudent.HW1
    @test 1 <= fx([2]) <= 4
    @test 1 <= fy([2]) <= 4
    @test 1 <= fy([1, 2, 4]) <= 7

    f(x, y) = sqrt(x^2 + y^2)
    results = HW1.evaluate(f, "test1234@colorado.edu")
    @test results.score == 1

    results = HW1.evaluate(f)
    @test results.score == 1

    g(x, y) = x+y
    results = HW1.evaluate(g, "test1234@colorado.edu")
    @test results.score == 0
end


text:
# This file was generated, do not modify it. # hide
params(krb) |> pprint


text:
import AtomicLevels: AbstractOrbital, symmetry

"""
    Conjugate(orbital)

Type representing the conjugation of an `orbital`.

# Examples

```jldoctest
julia> Conjugate(:a)
:a†
```
"""
struct Conjugate{O}
    orbital::O
end

function Base.show(io::IO, co::Conjugate{O}) where O
    show(io, co.orbital)
    write(io, "†")
end

"""
    conj(o::AbstractOrbital)

Convenience function to conjugate an `AbstractOrbital`.

# Examples

```jldoctest
julia> conj(o"1s")
1s†
```
"""
Base.conj(o::O) where {O<:AbstractOrbital} = Conjugate(o)

"""
    conj(o::Conjugate)

Convenience function to unconjugate a conjugated orbital.

# Examples

```jldoctest
julia> conj(Conjugate(:a))
:a
```
"""
Base.conj(co::Conjugate{O}) where O = co.orbital

AtomicLevels.symmetry(co::Conjugate{O}) where O = symmetry(co.orbital)

export Conjugate


text:
using Wavelets
using PyPlot

n = 2^11
x0 = testfunction(n,"Doppler")
x = x0 + 0.05*randn(n)
y = denoise(x,TI=true)
f, ax = plt.subplots(3,1, sharex=true)

cc=0.6
ax[1][:plot](x0,"b")
ax[1][:set_title]("original",fontsize=11)
ax[2][:plot](x,"g")
ax[2][:set_title]("noisy",fontsize=11)
ax[3][:plot](y,"r")
ax[3][:set_title]("denoised (TI VisuShrink)",fontsize=11)
for i=1:length(ax)
    ax[i][:set_xlim]([0,n])
    ax[i][:set_ylim]([-cc,cc])
end
f[:set_size_inches](6,5)

savefig("denoise_doppler.png")


text:
#
# SurfaceNets in Julia
#
# Ported from the Javascript implementation by Mikola Lysenko (C) 2012
# https://github.com/mikolalysenko/mikolalysenko.github.com/blob/master/Isosurface/js/surfacenets.js
# MIT License
#
# Based on: S.F. Gibson, "Constrained Elastic Surface Nets". (1998) MERL Tech Report.
#

#Look up Table
include("lut/sn.jl")

"""
Generate a mesh using naive surface nets.
This takes the center of mass of the voxel as the vertex for each cube.
"""
function surface_nets(data::Vector{T}, dims,eps,scale,origin) where {T}

    vertices = Point{3,T}[]
    faces = Face{4,Int}[]

    sizehint!(vertices,ceil(Int,maximum(dims)^2))
    sizehint!(faces,ceil(Int,maximum(dims)^2))

    n = 0
    R = Array{Int}([1, (dims[1]+1), (dims[1]+1)*(dims[2]+1)])
    buf_no = 1

    buffer = fill(zero(Int),R[3]*2)

    #March over the voxel grid
    zi = 0
    @inbounds while zi<dims[3]-1

        # m is the pointer into the buffer we are going to use.
        # This is slightly obtuse because javascript does not have good support for packed data structures, so we must use typed arrays :(
        # The contents of the buffer will be the indices of the vertices on the previous x/y slice of the volume
        m = 1 + (dims[1]+1) * (1 + buf_no * (dims[2]+1))

        yi=0
        @inbounds while yi<dims[2]-1

            xi=0
            @inbounds while xi < dims[1]-1

                # Read in 8 field values around this vertex and store them in an array
                # Also calculate 8-bit mask, like in marching cubes, so we can speed up sign checks later
                @inbounds grid = (data[n+1],
                                  data[n+2],
                                  data[n+dims[1]+1],
                                  data[n+dims[1]+2],
                                  data[n+dims[1]*2+1 + dims[1]*(dims[2]-2)],
                                  data[n+dims[1]*2+2 + dims[1]*(dims[2]-2)],
                                  data[n+dims[1]*3+1 + dims[1]*(dims[2]-2)],
                                  data[n+dims[1]*3+2 + dims[1]*(dims[2]-2)])

                mask = _get_cubeindex(grid, 0)

                # Check for early termination if cell does not intersect boundary
                if mask == 0x00 || mask == 0xff
                    xi += 1
                    n += 1
                    m += 1
                    continue
                end

                #Sum up edge intersections
                edge_mask = sn_edge_table[mask]

                _sn_add_verts!(xi, yi, zi, vertices, grid, edge_mask, buffer, m, scale, origin, eps, T, Val(true))

                #Now we need to add faces together, to do this we just loop over 3 basis components
                x = (xi,yi,zi)
                for i=0:2
                    #The first three entries of the edge_mask count the crossings along the edge
                    if (edge_mask & (1<<i)) == 0
                        continue
                    end

                    # i = axes we are point along.  iu, iv = orthogonal axes
                    iu = (i+1)%3
                    iv = (i+2)%3

                    #If we are on a boundary, skip it
                    if (x[iu+1] == 0 || x[iv+1] == 0)
                        continue
                    end

                    #Otherwise, look up adjacent edges in buffer
                    du = R[iu+1]
                    dv = R[iv+1]

                    #Remember to flip orientation depending on the sign of the corner.
                    if (mask & 0x01) != 0x00
                        push!(faces,Face{4,Int}(buffer[m+1]+1, buffer[m-du+1]+1, buffer[m-du-dv+1]+1, buffer[m-dv+1]+1));
                    else
                        push!(faces,Face{4,Int}(buffer[m+1]+1, buffer[m-dv+1]+1, buffer[m-du-dv+1]+1, buffer[m-du+1]+1));
                    end
                end
                xi += 1
                n += 1
                m += 1
            end
            yi += 1
            n += 1
            m += 2
        end
        zi += 1
        n+=dims[1]
        buf_no = xor(buf_no,1)
        R[3]=-R[3]
    end
    #All done!  Return the result

    vertices, faces # faces are quads, indexed to vertices
end

"""
Generate a mesh using naive surface nets.
This takes the center of mass of the voxel as the vertex for each cube.
"""
function surface_nets(f::Function, dims::NTuple{3,Int},eps,scale,origin)

    # TODO
    T = Float64

    vertices = Point{3,T}[]
    faces = Face{4,Int}[]

    sizehint!(vertices,ceil(Int,maximum(dims)^2))
    sizehint!(faces,ceil(Int,maximum(dims)^2))

    n = 0
    R = Array{Int}([1, (dims[1]+1), (dims[1]+1)*(dims[2]+1)])
    buf_no = 1

    buffer = fill(zero(Int),R[3]*2)

    grid = Vector{Float64}(undef,8)

    #March over the voxel grid
    zi = 0
    @inbounds while zi<dims[3]-1

        # m is the pointer into the buffer we are going to use.
        # This is slightly obtuse because javascript does not have good support for packed data structures, so we must use typed arrays :(
        # The contents of the buffer will be the indices of the vertices on the previous x/y slice of the volume
        m = 1 + (dims[1]+1) * (1 + buf_no * (dims[2]+1))

        yi=0
        @inbounds while yi<dims[2]-1

            xi=0
            @inbounds while xi < dims[1]-1

                # Read in 8 field values around this vertex and store them in an array
                points = (Point{3,Float64}(xi,yi,zi).* scale + origin,
                          Point{3,Float64}(xi+1,yi,zi).* scale + origin,
                          Point{3,Float64}(xi,yi+1,zi).* scale + origin,
                          Point{3,Float64}(xi+1,yi+1,zi).* scale + origin,
                          Point{3,Float64}(xi,yi,zi+1).* scale + origin,
                          Point{3,Float64}(xi+1,yi,zi+1).* scale + origin,
                          Point{3,Float64}(xi,yi+1,zi+1).* scale + origin,
                          Point{3,Float64}(xi+1,yi+1,zi+1).* scale + origin)

                if xi == 0
                    for i = 1:8
                        grid[i] = f(points[i])
                    end
                else
                    grid[1] = grid[2]
                    grid[2] = f(points[2])
                    grid[3] = grid[4]
                    grid[4] = f(points[4])
                    grid[5] = grid[6]
                    grid[6] = f(points[6])
                    grid[7] = grid[8]
                    grid[8] = f(points[8])
                end

                # Also calculate 8-bit mask, like in marching cubes, so we can speed up sign checks later
                mask = _get_cubeindex(grid, 0)

                # Check for early termination if cell does not intersect boundary
                if mask == 0x00 || mask == 0xff
                    xi += 1
                    n += 1
                    m += 1
                    continue
                end

                #Sum up edge intersections
                edge_mask = sn_edge_table[mask]

                # add vertices
                _sn_add_verts!(xi, yi, zi, vertices, grid, edge_mask, buffer, m, scale, origin, eps, T, Val(false))

                #Now we need to add faces together, to do this we just loop over 3 basis components
                x = (xi,yi,zi)
                for i=0:2
                    #The first three entries of the edge_mask count the crossings along the edge
                    if (edge_mask & (1<<i)) == 0
                        continue
                    end

                    # i = axes we are point along.  iu, iv = orthogonal axes
                    iu = (i+1)%3
                    iv = (i+2)%3

                    #If we are on a boundary, skip it
                    if (x[iu+1] == 0 || x[iv+1] == 0)
                        continue
                    end

                    #Otherwise, look up adjacent edges in buffer
                    du = R[iu+1]
                    dv = R[iv+1]

                    #Remember to flip orientation depending on the sign of the corner.
                    if (mask & 0x01) != 0x00
                        push!(faces,Face{4,Int}(buffer[m+1]+1, buffer[m-du+1]+1, buffer[m-du-dv+1]+1, buffer[m-dv+1]+1));
                    else
                        push!(faces,Face{4,Int}(buffer[m+1]+1, buffer[m-dv+1]+1, buffer[m-du-dv+1]+1, buffer[m-du+1]+1));
                    end
                end
                xi += 1
                n += 1
                m += 1
            end
            yi += 1
            n += 1
            m += 2
        end
        zi += 1
        n+=dims[1]
        buf_no = xor(buf_no,1)
        R[3]=-R[3]
    end
    #All done!  Return the result

    vertices, faces # faces are quads, indexed to vertices
end

@inline function _sn_add_verts!(xi, yi, zi, vertices, grid, edge_mask, buffer, m, scale, origin, eps, T, translate_pt)
    v = Point{3,T}(zero(T),zero(T),zero(T))
    e_count = 0

    #For every edge of the cube...
    @inbounds for i=0:11

        #Use edge mask to check if it is crossed
        if (edge_mask & (1<<i)) == 0
            continue
        end

        #If it did, increment number of edge crossings
        e_count += 1

        #Now find the point of intersection
        e0 = cube_edges[(i<<1)+1]       #Unpack vertices
        e1 = cube_edges[(i<<1)+2]
        g0 = grid[e0]                 #Unpack grid values
        g1 = grid[e1]
        t  = g0 - g1                 #Compute point of intersection
        if abs(t) > eps
            t = g0 / t
        else
            continue
        end

        #Interpolate vertices and add up intersections (this can be done without multiplying)
        # TODO lut table change may have made this incorrect
        # in which case e0=e0-1 and e1=e1-1
        xj, yj, zj = zero(T), zero(T), zero(T)
        a = e0 & 1
        b = e1 & 1
        (a != 0) && (xj += one(T))
        (a != b) && (xj += (a != 0 ? - t : t))
        a = e0 & 2
        b = e1 & 2
        (a != 0) && (yj += one(T))
        (a != b) && (yj += (a != 0 ? - t : t))
        a = e0 & 4
        b = e1 & 4
        (a != 0) && (zj += 1.0)
        (a != b) && (zj += (a != 0 ? - t : t))
        v += Point{3,T}(xj,yj,zj)

    end # edge check

    #Now we just average the edge intersections and add them to coordinate
    s = 1.0 / e_count
    if typeof(translate_pt) == Val{true}
        @inbounds v = (Point{3,T}(xi,yi,zi)  + s .* v) .* scale + origin
    else
        @inbounds v = (Point{3,T}(xi,yi,zi) + s .* v)# * scale[i] + origin[i]
    end

    #Add vertex to buffer, store pointer to vertex index in buffer
    buffer[m+1] = length(vertices)
    push!(vertices, v)
end

#
# Constructors
#

struct NaiveSurfaceNets{T} <: AbstractMeshingAlgorithm
    iso::T
    eps::T
    reduceverts::Bool
end

NaiveSurfaceNets(;iso::T1=0.0, eps::T2=1e-3, reduceverts::Bool=true) where {T1, T2} = NaiveSurfaceNets{promote_type(T1, T2)}(iso, eps, reduceverts)
NaiveSurfaceNets(iso) = NaiveSurfaceNets(iso=iso)
NaiveSurfaceNets(iso,eps) = NaiveSurfaceNets(iso=iso,eps=eps)

function (::Type{MT})(sdf::SignedDistanceField, method::NaiveSurfaceNets) where {MT <: AbstractMesh}
    bounds = sdf.bounds
    orig = origin(bounds)
    w = widths(bounds)
    scale = w ./ Point(size(sdf) .- 1)  # subtract 1 because an SDF with N points per side has N-1 cells

    d = vec(sdf.data)

    # Run iso surface additions here as to not
    # penalize surface net inner loops, and we are using a copy anyway
    if method.iso != 0.0
        for i = eachindex(d)
            d[i] -= method.iso
        end
    end

    vts, fcs = surface_nets(d,
                            size(sdf.data),
                            method.eps,
                            scale,
                            orig)
    MT(vts, fcs)::MT
end

function (::Type{MT})(f::Function, bounds::HyperRectangle, size::NTuple{3,Int}, method::NaiveSurfaceNets) where {MT <: AbstractMesh}
    orig = origin(bounds)
    w = widths(bounds)
    scale = w ./ Point(size .- 1)  # subtract 1 because an SDF with N points per side has N-1 cells

    # TODO ISO val

    vts, fcs = surface_nets(f,
                            size,
                            method.eps,
                            scale,
                            orig)
    MT(vts, fcs)::MT
end

function (::Type{MT})(f::Function, bounds::HyperRectangle, method::NaiveSurfaceNets;size::NTuple{3,Int}=(128,128,128)) where {MT <: AbstractMesh}
    MT(f,bounds,size,method)
end


text:
using ArgParse

function parse_commandline()
    s = ArgParseSettings(description = """This program is used to generate profile lines""")
    @add_arg_table! s begin
        """--nproc_old"""
        help = "number of slices of the old mesh"
        arg_type = Int64
        required = true
        """--old_mesh_dir"""
        help = "directory holds proc*_reg1_solver_data.bin"
        arg_type = String
        required = true
        """--old_model_dir"""
        help = "directory holds proc*_reg1_<model_tag>.bin"
        arg_type = String
        required = true
        """--model_tags"""
        help = "comma delimited string, e.g. vsv,vsh,rho"
        arg_type = String
        required = true
        """--output_file"""
        help = "output directory for interpolated model files"
        arg_type = String
        required = true
        """--region"""
        help = "lon1/lat1/lon2/lat2/dep1/dep2, if dep1==dep2, cut horizontally, otherwise cut vertically"
        arg_type = String
        required = true
        """--npts"""
        help = "hnpts/vnpts if cut vertically, lon_npts/lat_npts if cut horizontally, should be interger"
        arg_type = String
        required = true
    end
    return parse_args(s)
end

max_stars_repo_name: ziyixi/seisflow

max_issues_repo_name: ziyixi/seisflow

max_forks_repo_name: ziyixi/seisflow