text:
module nt_signal

using PyPlot

export imageplot, resize_image, load_image, cconvol, upsampling, subsampling, plot_wavelet, snr, perform_wavortho_transf

#######################################
function imageplot(g, str="", a=-1, b=-1, c=-1)
	if size(g,3)==1
		# special for b&w images.
		g = g[:,:,1];
	end
	if c>0
		subplot(a,b,c);
	end
	imshow(g, interpolation="nearest")
	set_cmap("gray")
	axis("off")
	if str!=""
		title(str);
	end
end

#######################################
function load_image(name,N)
	g = imread(name);
	g = resize_image(g, N);
	return g;
end

#######################################
function resize_image(f, N)
	P = size(f,1);
	# add 1 pixel boundary
	g = f;
	g = cat(2, g, g[:,1,:]);
	g = cat(1, g, g[1,:,:]);
	# interpolate
	t = range(1, stop=P, length=N);
	ti = floor(t); tj = ceil(t);
	fi = t-floor(t); fj = 1-fi;
	h = zeros(N,N,size(f,3));
	for s=1:size(f,3)
	    h[:,:,s] = g[ti,ti,s] .* (fj*fj') + g[tj,tj,s] .* (fi*fi') + g[ti,tj,s] .* (fj*fi') + g[tj,ti,s] .* (fi*fj');
	end
	return h;
end

###################################################
function trim_dim(x)
	# remove empty trailing dimensions
	if size(x,3)==1
		x = x[:,:,1];
		if size(x,2)==1
			x = x[:,1];
		end
	end
	return x;
end

###################################################
function subsampling(x,d=1,p=2)

	# downsampling - subsampling along dimension d
	#
	#   y = subsampling(x,d,p);
	#
	#   default is p==2, d==1
	#
	#   Copyright (c) 2009 Gabriel Peyre

	if d==1
	        y = x[1:p:end,:,:];
	elseif d==2
	        y = x[:,1:p:end,:];
	elseif d==3
	        y = x[:,:,1:p:end];
	else
	        error("Not implemented");
	end
	y = trim_dim(y);

	return y;

end # subsampling

###################################################
function upsampling(x,d=1,p=2)

	# upsampling - add p zeros between samples along dimension d
	#
	#   y = upsampling(x,d,p);
	#
	#   default is p==2, d==1
	#
	#   Copyright (c) 2009 Gabriel Peyre

	if d==1
	        y = zeros(p*size(x,1),size(x,2),size(x,3));
	        y[1:p:end,:,:] = x;
	elseif d==2
	        y = zeros(size(x,1),p*size(x,2),size(x,3));
	        y[:,1:p:end,:] = x;
	elseif d==3
	        y = zeros(size(x,1),size(x,2),p*size(x,3));
	        y[:,:,1:p:end] = x;
	else
	        error("Not implemented");
	end
	y = trim_dim(y);

	return y;
end # upsampling

###################################################
function cconvol(x,h,d=1)

	# cconvol - spatial domain circular convolution
	#
	#   y = cconv(x,h,d);
	#
	#   Circular convolution along dimension |d|.
	#
	#   If p=length(h) is odd then h((p+1)/2) corresponds to the zero position.
	#   If p is even, h(p/2) corresponds to the zero position.
	#
	#   Copyright (c) 2009 Gabriel Peyre

	if d==2
	    return ( cconvol( x',h, 1 )' );
	elseif d==3
		# permute do not exists in python
	    return permute( cconvol( permute(x,[3, 2, 1]),h), [3 2 1]);
	end

	p = length(h);
	if mod(p,2)==0
	    pc = p/2;
	else
	    pc = (p+1)/2;
	end

	y = zeros(size(x));
	for i=1:length(h)
	    y = y + h[i]*circshift(x,i-pc);
	end
	return y;

end # cconvol

#####################################################
function plot_wavelet(fW, Jmin=0)
    #    plot_wavelet - plot wavelets coefficients.
    #    U = plot_wavelet(fW, Jmin):
    #    Copyright (c) 2014 Gabriel Peyre

    function rescaleWav(A)
        v = maximum( abs(A[:]) );
        B = copy(A)
        if v > 0
            B = .5 + .5 * A / v;
        end
        return B;
    end
    function rescale(x,a=0,b=1)
		m = minimum(x[:]);
		M = maximum(x[:]);
		if M-m<1e-10
		    y = x;
		else
		    y = (b-a) * (x-m)/(M-m) + a;
		end
		return y;
	end ## end of rescale
    ##
    n = size(fW,1);
    Jmax = log2(n) - 1;
    U = copy(fW);
    for j = Jmax:-1:Jmin
    	L = 1:2^j; H = (2^j)+1:2^(j+1);
        U[L,H] = rescaleWav( U[L,H] );
        U[H,L] = rescaleWav( U[H,L] );
        U[H,H] = rescaleWav( U[H,H] );
    end
    # coarse scale
    L = 1:2^Jmin;
    U[L,L] = rescale(U[L,L]);
    # plot underlying image
    imageplot(U);
    # display crosses
    for j=Jmax:-1:Jmin
        plot([0, 2^(j+1)], [2^j, 2^j], "r");
        plot([2^j, 2^j], [0, 2^(j + 1)], "r");
    end
    # display box
    plot([0, n], [0, 0], "r");
    plot([0, n], [n, n], "r");
    plot([0, 0], [0, n], "r");
    plot([n, n], [0, n], "r");
    return U;

end # plot_wavelet

#########################################
function subselectdim(f,sel,d)
	g = [];
	if d==1
	    g = f[sel,:,:];
	elseif d==2
	    g = f[:,sel,:];
	elseif d==3
	    g = f[:,:,sel];
	end
	return trim_dim(g);
end

#########################################
function subassign(f,sel,g)
	d = ndims(f);
	if d==1
	    f[sel] = g;
	elseif d==2
	    f[sel,sel] = g;
	elseif d==3
	    f[sel,sel,sel] = g;
	end
	return f;
end # subassign

#########################################
function subselect(f,sel)
	d = ndims(f)
    if d==1
        return f[sel];
    elseif d==2
        return f[sel,sel];
    elseif d==3
        return f[sel,sel,sel];
	end
	return [];
end

##############################################################
function perform_wavortho_transf(f0,Jmin,dir,h)

	# perform_wavortho_transf - compute orthogonal wavelet transform
	#
	#   fw = perform_wavortho_transf[f,Jmin,dir,options);
	#
	#   You can give the filter in options.h.
	#
	#   Works in arbitrary dimension.
	#
	#   Copyright (c) 2009 Gabriel Peyre

	n = size(f0,1);
	Jmax = log2(n)-1;

	g = [0, h[length(h):-1:2]] .* (-1).^(1:length(h));

	f = copy(f0)
	if dir==1
	    ### FORWARD ###
	    for j=Jmax:-1:Jmin
	        sel = 1:2^(j+1);
	        a = subselect(f,sel);
	        for d=1:ndims(f)
	            a = cat(d, subsampling(cconvol(a,h,d),d), subsampling(cconvol(a,g,d),d) );
	        end
	        f = subassign(f,sel,a);
	    end
	else
	    ### FORWARD ###
	    for j=Jmin:Jmax
	        sel = 1:2^(j+1);
	        a = subselect(f,sel);
	        for d=1:ndims(f)
	            w = subselectdim(a,2^j+1:2^(j+1),d);
	            a = subselectdim(a,1:2^j,d);
	            a = cconvol(upsampling(a,d),reverse(h),d) + cconvol(upsampling(w,d),reverse(g),d);
	        end
	        f = subassign(f,sel,a);
	    end
	end
	return f;
end # perform_wavortho_transf

#####################################################
function snr(x,y)

	# snr - signal to noise ratio
	#
	#   v = snr(x,y);
	#
	# v = 20*log10( norm(x[:]) / norm(x[:]-y[:]) )
	#
	#   x is the original clean signal (reference).
	#   y is the denoised signal.
	#
	#   Copyright (c) 2008 Gabriel Peyre

	return 20*log10(norm(x[:])/norm(x[:]-y[:]));

end # snr

end ## module nt_signal


text:
using ActiveContours
using Test
using TestSetExtensions
using LazySets
using LazySets: Ellipsoid
using LinearAlgebra
using LinearAlgebra: Diagonal

ext: jl


text:
# this file used something like:

# julia-0.4 -e 'using Go; d = undump("evol-1.dump")[1]; replay_direct(d, "514/10000  13/25   cb91577f4ddb175e:1   bt e7d7a5f40ef149a5:2    3 sc  356 mv  24 sp  337 st")'

# where the quoted string is cut+paste from the log

# here, a plays first (black)
function play(a::Vector{UInt8}, b::Vector{UInt8}, board_size, max_moves, display_move, display_pass)
    passed, t, h = 0, black, [Position{board_size}()]
    # white must pass last
    while (passed < 2 || t == white) && h[end].stats.moves < max_moves
        played = false
        for m in moves(a, h, t)
            p = Position(h[end])
            move!(p, t, m...)
            if length(filter(x -> x.board == p.board, h)) > 0
                println("repeated position")
            else
                push!(h, p)
                display_move(p, m...)
                played = true
                break
            end
        end
        if played
            passed = 0
        else
            passed += 1
            push!(h, pass!(t, Position(h[end])))
            display_pass(h[end], t)
        end
        a, b, t = b, a, other(t)
    end
    h[end]
end

null_display(p, x, y) = nothing
null_display(p, t) = nothing

function board_display(p, x, y)
    t = point(p, x, y)
    @printf("\n\n%d (%3.1f): %s at (%d,%d)\n", p.stats.moves, p.stats.score, t == black ? "black" : "white", x, y)
    println(join(fmtpoint(p.board), "\n"))
end

function pass_display(p, t)
    @printf("\n\n%d (%3.1f): %s pass\n", p.stats.moves, p.stats.score, t == black ? "black" : "white")
end

# for replay, pick a game from the log.  perhaps
# 1000/1000   6/20   ceaa93b88f34d14f:1   bt fb9b4e8b26cd7702:3    9 sc   81 mv   4 sp  70 st
# and calculate the seed.  in this case (1000-1)*20+6
# the read in the population and replay:
# > d = undump("evol-1.dump");
# > replay_direct(d["ceaa93b88f34d14f"], d["fb9b4e8b26cd7702"], 9, 81, 20*999+6)

function replay_direct(a::Vector{UInt8}, b::Vector{UInt8}, board_size, max_moves)
    p = play(a, b, board_size, max_moves, board_display, pass_display)
    println("\n\nfinal position:")
    println(p)
    p
end

function replay_direct(d::Dict{AbstractString, Vector{UInt8}}, line::AbstractString; board_size=19, max_moves=1000)
# 50/10000  10/25 ! 8414ac2df622accd:38  bt c340be963d5fc19f:5   14 sc  100 mv   5 sp  67 st
# 50/10000  11/25   c340be963d5fc19f:6   bt 6b34dd72635a19e6:19   2 sc   45 mv   6 sp  45 st
    p = r"^\s*(?P<i>\d+)/(?P<n>\d+)\s+(?P<j>\d+)/(?P<m>\d+) (?P<surprise>(?:!| )) (?P<a>[a-f0-9]+):\d+\s+(?:>|=|<)\s+(?P<b>[a-f0-9]+):\d+.*$"
    m = match(p, line)
    println(m.captures)
    i, n, j, m, surprise, a, b = m.captures
    seed = (parse(Int, i)-1) * parse(Int, m) + parse(Int, j)
    println("$a v $b (seed $(seed))")
    replay_direct(d[a], d[b], board_size, max_moves)
end

function replay_pair(path, a, b; board_size=19, max_moves=1000)
    d = undump(path)[1]
    replay_direct(d[a], d[b], board_size, max_moves)
end

function replay_latest(; log_path="evol-2.log", dump_path="evol-2.dump",
                       a=1, b=2, n=100, fraction=0.9,
                       board_size=19, max_moves=1000)
    events = parse_log(log_path, dump_path, n, fraction)
    popn = Vector{Vector{Individual}}()
    for e in events
        expand!(popn, e)
    end
    named = undump(dump_path)[1]
    ea, eb = named[popn[end][a].tag], named[popn[end][b].tag]
    if ea > eb
        ea, eb = eb, ea
        a, b = b, a
    end
    println("black: $(name(ea)):$(a); white: $(name(eb)):$(b)\n")
    replay_direct(ea, eb, board_size, max_moves)
end


text:
# Squared Exponential Function with ARD

@doc """
# Description
Constructor for the ARD Squared Exponential kernel (covariance)

k(x,x') = σ²exp(-(x-x')ᵀL⁻²(x-x')/2), where L = diag(ℓ₁,ℓ₂,...)
# Arguments:
* `ll::Vector{Float64}`: Log of the length scale ℓ
* `lσ::Float64`: Log of the signal standard deviation σ
""" ->
type SEArd <: StationaryARD
    ℓ2::Vector{Float64}      # Log of Length scale
    σ2::Float64              # Log of Signal std
    SEArd(ll::Vector{Float64}, lσ::Float64) = new(exp(2.0*ll),exp(2.0*lσ))
end

function set_params!(se::SEArd, hyp::Vector{Float64})
    length(hyp) == num_params(se) || throw(ArgumentError("SEArd only has $(num_params(se)) parameters"))
    d = length(se.ℓ2)
    se.ℓ2 = exp(2.0*hyp[1:d])
    se.σ2 = exp(2.0*hyp[d+1])
end

get_params(se::SEArd) = [log(se.ℓ2)/2.0; log(se.σ2)/2.0]
get_param_names(k::SEArd) = [get_param_names(k.ℓ2, :ll); :lσ]
num_params(se::SEArd) = length(se.ℓ2) + 1

metric(se::SEArd) = WeightedSqEuclidean(1.0./(se.ℓ2))
cov(se::SEArd, r::Float64) = se.σ2*exp(-0.5*r)

function grad_kern(se::SEArd, x::Vector{Float64}, y::Vector{Float64})
    r = distance(se, x, y)
    exp_r = exp(-0.5*r)
    wdiff = ((x-y).^2)./se.ℓ2

    g1   = se.σ2.*wdiff*exp_r   #dK_d(log ℓ)
    g2 = 2.0*se.σ2*exp_r        #dK_d(log σ)

    return [g1; g2]
end

function grad_stack!(stack::AbstractArray, se::SEArd, X::Matrix{Float64}, data::StationaryARDData)
    d = size(X,1)
    ck = cov(se, X, data)
    broadcast!(*, view(stack, :, :, 1:d), data.dist_stack, reshape(1.0./se.ℓ2, (1,1,d)), ck)
    stack[:,:, d+1] = 2.0 * ck
    return stack
end


text:
using JSON

function read_ils_table(file::String, jsonFile::String)
    @info "reading ILS file $file"
    if isfile(file) && isfile(jsonFile)
        jsonDict     = JSON.parsefile(jsonFile);
        ils_json     = jsonDict["ILS"]
        ilsFile      = Dataset(file);
        ils_grid     = getNC_var(ilsFile, ils_json["ils_grid"])
        ils_response = getNC_var(ilsFile, ils_json["ils_response"])
        dispersion   = getNC_var(ilsFile, ils_json["dispersion"])
        @info "ILS table size: " size(ils_response)
        close(ilsFile);
        return ils_grid, ils_response, dispersion
    else
        @error "ILS files don't exist" file isfile(file) jsonFile isfile(jsonFile)
        return nothing, nothing
    end
end

function getNC_var(fin, path)
    try
        loc = split(path, r"/")
        # println(loc)
        if length(loc) == 1
            return fin[path].var
        elseif length(loc) > 1
            gr = []
            for i in 1:length(loc) - 1
                if i == 1
                    gr = fin.group[loc[i]]
                else
                    gr = gr.group[loc[i]]
                end
            end
            return gr[loc[end]].var
        end
    catch e
        @error e
        println("Error in getNC_var ", path)
        return nothing
    end
end


text:
using GenomicAnnotations
using BioSequences
using Test

@testset "GenomicAnnotations" begin
    @testset "readgbk" begin
        s = "     gene            1..1"
        @test GenomicAnnotations.parseposition(s) == ("gene", Locus(1:1, '+'))
        s = "     gene            complement(order(3300..4037,4047..4052))"
        @test GenomicAnnotations.parseposition(s) == ("gene", Locus(3300:4052, '-', true, true, UnitRange{Int}[3300:4037, 4047:4052], false))

        chrs = readgbk("example.gbk")
        @test length(chrs) == 2
        @test chrs[2].name == "plasmid1"
    end

    chr = readgbk("example.gbk")[1]

    @testset "Extended methods" begin
        @test length(chr.genes[1]) == length(sequence(chr.genes[1]))
    end

    @testset "Gene properties" begin
        @test length(propertynames(chr.genes[1])) == 14
        @test chr.genes[2].locus_tag == "tag01"
        @test (chr.genes[2].locus_tag = "tag01") == "tag01"
        @test begin
            chr.genes[2].test = "New column"
            chr.genes[2].test == "New column" && all(ismissing, chr.genes[[1,3,4,5,6]].test)
        end
        @test get(chr.genes[1], :locus_tag, "") == ""
        @test get(chr.genes[2], :locus_tag, "") == "tag01"
        @test get(chr.genes[2], :locustag, "") == ""
        @test begin
            GenomicAnnotations.pushproperty!(chr.genes[2], :db_xref, "GI:123")
            chr.genes[2].db_xref == ["GI:1293614", "GI:123"]
        end
        @test GenomicAnnotations.vectorise(Union{Missing,Int}[1,1,1]) == [[1],[1],[1]]
    end

    @testset "Iteration" begin
        @test length([g.locus_tag for g in chr.genes]) == 7
    end

    @testset "Adding/removing genes" begin
        addgene!(chr, "CDS", Locus(300:390, '+'), locus_tag = "tag04")
        @test chr.genes[end].locus_tag == "tag04"
        delete!(chr.genes[end])
        @test chr.genes[end].locus_tag == "reg01"
    end

    @testset "@genes" begin
        @test @genes(chr, :feature == "CDS") == chr.genes[[2,4,6]]
        @test @genes(chr, iscomplement(gene)) == chr.genes[[5,6,7]]
        @test @genes(chr, :feature == "CDS", !iscomplement(gene)) == chr.genes[[2,4]]
        @test @genes(chr, length(gene) < 300)[1] == chr.genes[2]
        @test length(@genes(chr, get(gene, :locus_tag, "") == "")) == 3
    end

    @testset "Broadcast" begin
        # Broadcasted assignment on existing property
        chr.genes[3:4].gene .= "AXL2P"
        @test all(chr.genes[3:4].gene .== "AXL2P")
        # Broadcasted assignment on previously missing property
        chr.genes[3:4].newproperty .= true
        @test all(chr.genes[3:4].newproperty .== true)
        # Broadcasted assignment with @genes
        @genes(chr, :feature == "gene").newproperty .= false
        @test all(chr.genes[[3,5]].newproperty .== false)
    end

    @testset "Locus" begin
        locus = Locus(1:1, '.', true, true, UnitRange{Int}[], false)
        @test Locus() == locus
        @test chr.genes[2].locus < chr.genes[4].locus
        @test chr.genes[2].locus == chr.genes[2].locus
        @test iscomplement(chr.genes[2]) == false
        @test iscomplement(chr.genes[5]) == true
    end

    seq = dna"atgtccatatacaacggtatctccacctcaggtttagatctcaacaacggaaccattgccgacatgagacagttaggtatcgtcgagagttacaagctaaaacgagcagtagtcagctctgcatctgaagccgctgaagttctactaagggtggataacatcatccgtgcaagaccaagaaccgccaatagacaacatatgtaa"
    @test sequence(chr.genes[2]) == seq
    @test length(chr.genes[2]) == length(seq)

    @testset "Chromosome" begin
        chr = Chromosome()
        @test chr.name == ""
        @test chr.sequence == dna""
        @test chr.header == ""
        @test chr.genes == Gene[]
        @test names(chr.genedata) == [:feature, :locus]
    end
end


text:
abstract type AbstractPolynomialVectorBasis{PT<:MP.AbstractPolynomialLike, PV<:AbstractVector{PT}} <: AbstractPolynomialBasis end

Base.length(basis::AbstractPolynomialVectorBasis) = length(basis.polynomials)
Base.copy(basis::AbstractPolynomialVectorBasis) = typeof(basis)(copy(basis.polynomials))

MP.nvariables(basis::AbstractPolynomialVectorBasis) = MP.nvariables(basis.polynomials)
MP.variables(basis::AbstractPolynomialVectorBasis) = MP.variables(basis.polynomials)
MP.monomialtype(::Type{<:AbstractPolynomialVectorBasis{PT}}) where PT = MP.monomialtype(PT)

empty_basis(B::Type{<:AbstractPolynomialVectorBasis{PT, Vector{PT}}}) where PT = B(PT[])
function MP.polynomialtype(basis::AbstractPolynomialVectorBasis{PT}, T::Type) where PT
    C = MP.coefficienttype(PT)
    U = MA.promote_operation(*, C, T)
    V = MA.promote_operation(+, U, U)
    return MP.polynomialtype(PT, V)
end
function MP.polynomial(f::Function, basis::AbstractPolynomialVectorBasis)
    return MP.polynomial(mapreduce(
        ip -> f(ip[1]) * ip[2], MA.add!, enumerate(basis.polynomials)))
end

function MP.polynomial(Q::AbstractMatrix, basis::AbstractPolynomialVectorBasis,
                       T::Type)
    n = length(basis)
    @assert size(Q) == (n, n)
    return MP.polynomial(mapreduce(row -> adjoint(basis.polynomials[row]) *
        mapreduce(col -> Q[row, col] * basis.polynomials[col], MA.add!, 1:n),
        MA.add!, 1:n), T)
end

"""
    struct FixedPolynomialBasis{PT<:MP.AbstractPolynomialLike, PV<:AbstractVector{PT}} <: AbstractPolynomialBasis
        polynomials::PV
    end

Polynomial basis with the polynomials of the vector `polynomials`.
For instance, `FixedPolynomialBasis([1, x, 2x^2-1, 4x^3-3x])` is the Chebyshev
polynomial basis for cubic polynomials in the variable `x`.
"""
struct FixedPolynomialBasis{PT<:MP.AbstractPolynomialLike, PV<:AbstractVector{PT}} <: AbstractPolynomialVectorBasis{PT, PV}
    polynomials::PV
end


text:
export SolverParameters, SolverState, Solver
export setup_coffee_lounge, add_coffee_break, solve
export load_snapshot
export make_solver_parameters
export init_solve, onestep_solve, finalize_solve
import Base.Meta: quot

@compat abstract type SolverMethod end # An enum type to identify the solver in validation functions

@compat abstract type InternalSolverState end # All the state a solver needs to update an iteration
const SolverParameters = Dict{Symbol,Any}

immutable Solver{T<:SolverMethod}
  method        :: T
  params        :: SolverParameters
  coffee_lounge :: Any # forward declaration
end

Solver{T}(method::T, params::SolverParameters) = begin
    validate_parameters(method, params)
    Solver(method, params, CoffeeLounge())
end

type SolverState{T<:InternalSolverState}
  iter               :: Int
  obj_val            :: Float64
  losses             :: Dict
  internal           :: T
end
function copy_solver_state!{T<:InternalSolverState}(dst::SolverState{T}, src::SolverState{T})
  dst.iter = src.iter
  dst.obj_val = src.obj_val
  dst.losses = src.losses
  dst.internal = src.internal
end

SolverState{T<:InternalSolverState}(internal::T) = SolverState{T}(0, Inf, Dict(), internal)

@compat abstract type SolverStateSnapshot end # Just the serializable part of the solver state, for snapshot files

function make_solver_parameters(;kwargs...)
    localDict = @compat Dict(:max_iter => Inf,
           :regu_coef => 0.0005,
           :load_from => "")
	merge(localDict,
          SolverParameters(kwargs))
end

function validate_parameters(params::SolverParameters, args...)
    if params[:regu_coef] < 0
        error("regu_coef must be non-negative")
    end
    for a in args
        if !haskey(params, a)
            error("Must provide parameter $a")
        end
    end
end

############################################################
#  API functions to be implemented by each solver instance
############################################################

function get_statistic(state::SolverState, name::AbstractString)
  if haskey(state.losses, name)
    return state.losses[name]
  else
    return NaN
  end
end

function get_statistic(state::SolverState, name::Symbol)
    if name in [:obj_val, :iter]
      @eval $(quot(state)).$name
    else
        @eval $(quot(state.internal)).$name
    end
end

function format_statistic(state::SolverState, name)
    format_statistic(state, name, get_statistic(state, name))
end

function format_statistic(state::SolverState, name, value)
    if name == :iter
        @sprintf("%s=%06d", name, value)
    else
        @sprintf("%s=%.8f", name, value)
    end
end

validate_parameters(method::SolverMethod, params::SolverParameters) = begin # should turn validate method-specific parameters
    validate_parameters(params)
    error("solver_state is not implemented for the method type $(typeof(method)), and params type $(typeof(params))")
end

function snapshot(state::SolverState) # should turn a SolverState{T} into a SolverStateSnapshot instance
    error("snapshot is not implemented for the state type $(typeof(state))")
end

function solver_state(net::Net, snapshot::SolverStateSnapshot) # should use a SolverStateSnapshot instance to create a SolverState{T}
    error("solver_state is not implemented for the net type $(typeof(net)), and snapshot type $(typeof(snapshot))")
end

function solver_state(solver::SolverMethod, net::Net, params::SolverParameters) # should create SolverState{T} from SolverParameters dictionary
    error("solver_state is not implemented for the solver type $(typeof(solver)), net type $(typeof(net)), and params type $(typeof(params))")
end

function update{T}(solver::Solver{T}, net::Net, state::SolverState) # should do one iteration of update
    error("update is not implemented for the solver type $(typeof(solver)), net type $(typeof(net)), and state type $(typeof(state))")
end
function shutdown(state::SolverState) # should shutdown the solver
    error("shutdown is not implemented for the $(typeof(state)) type")
end

############################################################
# General utilities that could be used by all solvers
############################################################
function load_snapshot(net::Net, path::AbstractString="", state=nothing)
  if isempty(path)
    return state
  end

  if endswith(path, ".hdf5") || endswith(path, ".h5")
    # load from HDF5 file, possibly exported from caffe, but training
    # from the beginning (iteration 0) as the solver state is not saved
    # in a HDF5 file
    if isfile(path)
      @info("Loading existing model from $path")
      h5open(path) do file
        load_network(file, net)
      end
    end
    return state
  else
    if endswith(path, ".jld")
      # load from some specific JLD sanpshot, the solver state is also
      # recovered
      filename = path
    else
      # automatically load from the latest snapshot in a directory
      filename = ""
      if isdir(path)
        # load the latest snapshot from the directory
        snapshots = glob(path, r"^snapshot-[0-9]+\.jld", sort_by=:mtime)
        if length(snapshots) > 0
          filename = joinpath(path, snapshots[end])
        end
      end
    end

    if !isempty(filename) && isfile(filename)
      @info("Loading existing model from $filename")
      jldopen(filename) do file
        load_network(file, net)
        return solver_state(net, read(file, SOLVER_STATE_KEY))
      end
    else
      return state
    end
  end
end

function stop_condition_satisfied(solver::Solver, state::SolverState, net::Net)
  if state.iter >= solver.params[:max_iter]
    return true
  end
  return false
end

############################################################
# General Solver Loop
############################################################

function init_solve(solver::Solver, net::Net)
    @debug("#DEBUG Checking network topology for back-propagation")
    check_bp_topology(net)

    state = solver_state(solver.method, net, solver.params)
    state = load_snapshot(net, solver.params[:load_from], state)

    # we init network AFTER loading. If the parameters are loaded from file, the
    # initializers will be automatically set to NullInitializer
    init(net)
    state.obj_val = forward(net, solver.params[:regu_coef])

    @debug("#DEBUG Initializing coffee breaks")
    setup(solver.coffee_lounge, state, net)

    # coffee break for iteration 0, before everything starts
    check_coffee_break(solver.coffee_lounge, state, net)

    return state
end

function finalize_solve(solver::Solver, net::Net, state::SolverState)
    shutdown(solver.coffee_lounge, net)
    shutdown(state)
end

function onestep_solve(solver::Solver, net::Net, state::SolverState)
    state.iter += 1
    layer_states = updatable_layer_states(net)

    backward(net, solver.params[:regu_coef])
    update(solver, net, state)

    # apply weight constraints
    for layer_state in layer_states
      for param in layer_state.parameters
        cons_every = param.constraint.every_n_iter
        if cons_every > 0 && state.iter % cons_every == 0
          constrain!(net.backend, param.constraint, param.blob)
        end
      end
    end

    state.obj_val = forward(net, solver.params[:regu_coef])

    state.losses = Dict()
    for i = 1:length(net.layers)
      if has_loss(net.layers[i])
        state.losses[net.layers[i].name] = net.states[i].loss
      end
    end
    check_coffee_break(solver.coffee_lounge, state, net)
end

function solve(solver::Solver, net::Net)
    state = init_solve(solver, net)
    do_solve_loop(solver, net, state)
    finalize_solve(solver, net, state)
end

function do_solve_loop(solver::Solver, net::Net, state::SolverState)
  @debug("#DEBUG Entering solver loop")
  while !stop_condition_satisfied(solver, state, net)
    onestep_solve(solver,net,state)

    if stop_condition_satisfied(solver, state, net)
      break
    end
  end
  return state
end

trainable_layers(net::Net) = filter(i -> has_param(net.layers[i]) && !is_frozen(net.states[i]),
                                    1:length(net.layers))

updatable_layer_states(net::Net) = [net.states[i] for i in trainable_layers(net)]


text:
using BipolarSphericalHarmonics
using Test
using BipolarSphericalHarmonics: kronindex, cache, monopolarharmonics, GSHPerm, VSHPerm

using LegendrePolynomials
using SphericalHarmonics
using VectorSphericalHarmonics
using VectorSphericalHarmonics: basisconversionmatrix
using SphericalHarmonicModes

using SphericalHarmonicArrays
using SphericalHarmonics: NorthPole, SouthPole
import WignerSymbols
using OffsetArrays
using LinearAlgebra
using StaticArrays
using Rotations
using WignerD

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

isapproxdefault(x, y) = isapprox(x, y, atol = 1e-14, rtol = 1e-8)
isapproxdefault(x::Nothing, y::Nothing) = true

cosχ(θ1, ϕ1, θ2, ϕ2) = cos(θ1)cos(θ2) + sin(θ1)sin(θ2)cos(ϕ1-ϕ2)
∂ϕ₂cosχ(θ1, ϕ1, θ2, ϕ2) = sin(θ1)sin(θ2)sin(ϕ1-ϕ2)

@testset "Clebsch-Gordan" begin

    # Computes the Clebsch-Gordan coefficient C_{l₁m₁l₂m₂}^{lm} for all valid l
    function clebschgordan(l1::Integer, m1::Integer, l2::Integer, m::Integer)
        m₂ = m - m1
        w3j = zeros(Cdouble, l1 + l2 + 1)
        jmin, jmax = extrema(BipolarSphericalHarmonics.degreerange(l1, l2, m))
        C = zeros(jmin:jmax)
        BipolarSphericalHarmonics.clebschgordan!(C, w3j, l1, m1,l2, m)
        return C
    end

    @testset "allocating" begin
        CG = clebschgordan(1,1,1,0)

        @test CG[0] ≈ WignerSymbols.clebschgordan(1,1,1,-1,0,0) ≈ 1/√3
        @test CG[1] ≈ WignerSymbols.clebschgordan(1,1,1,-1,1,0) ≈ 1/√2
        @test CG[2] ≈ WignerSymbols.clebschgordan(1,1,1,-1,2,0) ≈ 1/√6

        CG = clebschgordan(1,-1,1,0)

        @test CG[0] ≈ WignerSymbols.clebschgordan(1,-1,1,1,0,0) ≈ 1/√3
        @test CG[1] ≈ WignerSymbols.clebschgordan(1,-1,1,1,1,0) ≈ -1/√2
        @test CG[2] ≈ WignerSymbols.clebschgordan(1,-1,1,1,2,0) ≈ 1/√6

        CG = clebschgordan(1,0,1,0)

        @test CG[0] ≈ WignerSymbols.clebschgordan(1,0,1,0,0,0) ≈ -1/√3
        @test CG[1] ≈ WignerSymbols.clebschgordan(1,0,1,0,1,0) ≈ 0
        @test CG[2] ≈ WignerSymbols.clebschgordan(1,0,1,0,2,0) ≈ √(2/3)

        for l1 in 0:10, l2 in 0:10
            j = l1 + l2
            m = j
            CG = clebschgordan(l1,l1,l2,m)
            @test CG[j] ≈ 1
        end
    end

    @testset "prime factorization" begin
        C1 = zeros(200)
        C2 = zero(C1)
        C3 = zero(C1)
        W1 = zero(C1)
        W2 = zero(C2)
        W3 = zero(C3)
        BipolarSphericalHarmonics.clebschgordan!(C1, W1, 50, 0, 50, 0, primefactorization_cutoff = 0)
        BipolarSphericalHarmonics.clebschgordan!(C2, W2, 50, 0, 50, 0, primefactorization_cutoff = 200)
        BipolarSphericalHarmonics.clebschgordan!(C3, W3, 50, 0, 50, 0, primefactorization_cutoff = 50)
        @test C1 ≈ C2 ≈ C3
        C1 .= 0
        C2 .= 0
        C3 .= 0
        BipolarSphericalHarmonics.clebschgordan!(C1, W1, 50, 50, 50, 50, primefactorization_cutoff = 0)
        BipolarSphericalHarmonics.clebschgordan!(C2, W2, 50, 50, 50, 50, primefactorization_cutoff = 200)
        BipolarSphericalHarmonics.clebschgordan!(C3, W3, 50, 50, 50, 50, primefactorization_cutoff = 50)
        @test C1 ≈ C2 ≈ C3
    end
end

function _biposh(m, l1, l2, θ1, ϕ1, θ2, ϕ2, B = cache(SH(), max(l1, l2)))
    jrange = BipolarSphericalHarmonics.degreerange(l1, l2, m)
    Y = zeros(ComplexF64, LM(jrange, m))
    computePlmcostheta!(B.S1, θ1, l1)
    Y1 = computeYlm!(B.S1, θ1, ϕ1, l1)
    computePlmcostheta!(B.S2, θ2, l2)
    Y2 = computeYlm!(B.S2, θ2, ϕ2, l2)
    for m1 in -l1:l1
        m2 = m - m1
        abs(m2) <= l2 || continue
        BipolarSphericalHarmonics.clebschgordan!(B.C, B.W3j, l1, m1, l2, m)
        C = BipolarSphericalHarmonics.C_view_j1j2m_j(B.C, l1, l2, m, jrange)
        Y1_l1m1 = Y1[(l1, m1)]
        Y2_l2m2 = Y2[(l2, m2)]
        for ind in eachindex(Y, C)
            Y[ind] += C[ind] * Y1_l1m1 * Y2_l2m2
        end
    end
    return Y
end

@testset "kronindex" begin
    @testset "GSH" begin
        r = -1:1
        for (ind, (i, j)) in enumerate(Base.product(r,r))
            @test kronindex(GSH(), j, i) == ind
        end

        @testset "values" begin
            θ1, ϕ1, θ2, ϕ2 = pi/2, 0, pi/2, pi/3
            j, m, j1, j2 = 0, 0, 2, 2
            B1 = biposh(SH(), θ1, ϕ1, θ2, ϕ2, j, m, j1, j2)
            BG = biposh(GSH(), θ1, ϕ1, θ2, ϕ2, j, m, j1, j2)
            @test BG[kronindex(GSH(), 0, 0)] ≈ B1 atol=1e-14 rtol=1e-8
        end
    end

    @testset "VSH" begin
        struct TupMerge{N}
           t :: NTuple{N,Int}
        end
        Base.Tuple(t::TupMerge) = t.t
        Base.:(*)(a::TupMerge, b::TupMerge) = TupMerge((a.t..., b.t...))
        S = [TupMerge((i,j)) for i in -1:1, j in -1:1];
        SS = kron(S, S);
        for (ind,t) in zip(CartesianIndices(SS), SS)
            ttup = Tuple(t)
            @test kronindex(VSH(PB(), HelicityCovariant()), ttup...) == ind
            @test kronindex(VSH(PB(), SphericalCovariant()), ttup...) == ind
            indt11, indt12, indt21, indt22 = ttup
            indt11s = indt11 + 2
            indt21s = indt21 + 2
            ttups = (indt11s, indt12, indt21s, indt22)
            @test kronindex(VSH(PB(), Polar()), ttups...) == ind
            @test kronindex(VSH(PB(), Cartesian()), ttups...) == ind
        end

        # For Hansen, H⁻¹ⱼₘ(̂n) = Yⱼₘ(̂n)̂n, therefore its radial component is Yⱼₘ(̂n)
        # The radial basis vector ̂n is the same as the HelicityCovariant basis vector χ₀
        # The Hansen vector H⁻¹ⱼₘ(̂n) is equivalent to the PB vector P⁰ⱼₘ(̂n)
        @testset "values" begin
            θ1, ϕ1, θ2, ϕ2 = pi/2, 0, pi/2, pi/3
            j, m, j1, j2 = 0, 0, 2, 2
            B1 = biposh(SH(), θ1, ϕ1, θ2, ϕ2, j, m, j1, j2)

            BV = biposh(VSH(Hansen(), Polar()), θ1, ϕ1, θ2, ϕ2, j, m, j1, j2)
            @test BV[kronindex(VSH(Hansen(), Polar()), 1, -1, 1, -1)] ≈ B1 atol=1e-14 rtol=1e-8
            BV = biposh(VSH(Hansen(), HelicityCovariant()), θ1, ϕ1, θ2, ϕ2, j, m, j1, j2)
            @test BV[kronindex(VSH(Hansen(), HelicityCovariant()), 0, -1, 0, -1)] ≈ B1 atol=1e-14 rtol=1e-8

            BV = biposh(VSH(PB(), Polar()), θ1, ϕ1, θ2, ϕ2, j, m, j1, j2)
            @test BV[kronindex(VSH(PB(), Polar()), 1, 0, 1, 0)] ≈ B1 atol=1e-14 rtol=1e-8
            BV = biposh(VSH(PB(), HelicityCovariant()), θ1, ϕ1, θ2, ϕ2, j, m, j1, j2)
            @test BV[kronindex(VSH(PB(), HelicityCovariant()), 0, 0, 0, 0)] ≈ B1 atol=1e-14 rtol=1e-8
        end
    end

    @testset "broadcasting" begin
        @testset "GSH" begin
            for r1_st in -1:1, r1_end in -1:1
                r1 = r1_st:r1_end
                for i2 in -1:1
                    ki = kronindex.(GSH(), r1, i2)
                    @test ki == kronindex.(GSH(), collect(r1), collect(i2))
                    @test length(ki) == length(r1)
                    for (ind, i1) in enumerate(r1)
                        @test ki[ind] == kronindex(GSH(), i1, i2)
                    end
                end
            end
            for r2_st in -1:1, r2_end in -1:1
                r2 = r2_st:r2_end
                for i1 in -1:1
                    ki = kronindex.(GSH(), i1, r2)
                    @test ki == kronindex.(GSH(), collect(i1), collect(r2))
                    @test length(ki) == length(r2)
                    for (ind, i2) in enumerate(r2)
                        @test ki[ind] == kronindex(GSH(), i1, i2)
                    end
                end
            end
            for r2_st in -1:1, r2_end in -1:1
                r2 = r2_st:r2_end
                for r1_st in -1:1, r1_end in -1:1
                    r1 = r1_st:r1_end
                    if length(r1) != 1 && length(r2) != 1 && length(r1) != length(r2)
                        continue
                    end
                    @inferred Broadcast.broadcast(kronindex, Ref(GSH()), r1, r2)
                    ki = kronindex.(GSH(), r1, r2)
                    @test ki == kronindex.(GSH(), collect(r1), collect(r2))
                    els = ((x,y) -> (x,y)).(r1,r2)
                    for (ind, (i1, i2)) in enumerate(els)
                        @test ki[ind] == kronindex(GSH(), i1, i2)
                    end
                    if length(r2) == 1
                        @test length(ki) == length(r1)
                        for (ind, i1) in enumerate(r1)
                            @test ki[ind] == kronindex.(GSH(), i1, r2)[1]
                        end
                    end
                    if length(r1) == 1
                        @test length(ki) == length(r2)
                        for (ind, i2) in enumerate(r2)
                            @test ki[ind] == kronindex.(GSH(), r1, i2)[1]
                        end
                    end
                    if length(r1) == 1 && length(r2) == 1
                        @test length(ki) == 1
                        @test ki[1] == kronindex(GSH(), first(r1), first(r2))
                    end
                    if length(r1) == 0 || length(r2) == 0
                        @test length(ki) == 0
                    end
                end
            end
        end
    end
end

@testset "monopolarharmonics" begin
    lmax = 3
    @testset "SH" begin
        B = cache(SH(), lmax);
        YSH1 = computeYlm(pi/2, 0, lmax = lmax)
        YSH2 = computeYlm(pi/3, pi/4, lmax = lmax)
        BipolarSphericalHarmonics.monopolarharmonics!(B, pi/2, 0, pi/3, pi/4);
        Y1, Y2 = monopolarharmonics(B);
        for ((j,m), Y) in zip(first(SphericalHarmonicArrays.modes(Y1)), Y1)
            @test isapproxdefault(Y, YSH1[(j,m)])
        end
        for ((j,m), Y) in zip(first(SphericalHarmonicArrays.modes(Y2)), Y2)
            @test isapproxdefault(Y, YSH2[(j,m)])
        end
        BipolarSphericalHarmonics.monopolarharmonics!(B, pi/2, 0, pi/3, pi/4, ((1,2),(2,1)));
        for ((j,m), Y) in zip(first(SphericalHarmonicArrays.modes(Y1)), Y1)
            @test isapproxdefault(Y, YSH1[(j,m)])
        end
        for ((j,m), Y) in zip(first(SphericalHarmonicArrays.modes(Y2)), Y2)
            @test isapproxdefault(Y, YSH2[(j,m)])
        end
    end
    @testset "GSH" begin
        B = cache(GSH(), lmax);
        BipolarSphericalHarmonics.monopolarharmonics!(B, pi/2, 0, pi/3, pi/4);
        Y1, Y2 = monopolarharmonics(B);
        for ((j,m), Y) in zip(first(SphericalHarmonicArrays.modes(Y1)), Y1)
            @test isapproxdefault(Y, genspharm(j, m, pi/2, 0))
        end
        for ((j,m), Y) in zip(first(SphericalHarmonicArrays.modes(Y2)), Y2)
            @test isapproxdefault(Y, genspharm(j, m, pi/3, pi/4))
        end
        BipolarSphericalHarmonics.monopolarharmonics!(B, pi/2, 0, pi/3, pi/4, ((1,2),(2,1)));
        for ((j,m), Y) in zip(first(SphericalHarmonicArrays.modes(Y1)), Y1)
            @test isapproxdefault(Y, genspharm(j, m, pi/2, 0))
        end
        for ((j,m), Y) in zip(first(SphericalHarmonicArrays.modes(Y2)), Y2)
            @test isapproxdefault(Y, genspharm(j, m, pi/3, pi/4))
        end
    end
    @testset "VSH" begin
        B = cache(VSH(PB(), Polar()), lmax);
        BipolarSphericalHarmonics.monopolarharmonics!(B, pi/2, 0, pi/3, pi/4);
        Y1, Y2 = monopolarharmonics(B);
        for ((j,m), Y) in zip(first(SphericalHarmonicArrays.modes(Y1)), Y1)
            @test isapproxdefault(Y, vshbasis(PB(), Polar(), j, m, pi/2, 0))
        end
        for ((j,m), Y) in zip(first(SphericalHarmonicArrays.modes(Y2)), Y2)
            @test isapproxdefault(Y, vshbasis(PB(), Polar(), j, m, pi/3, pi/4))
        end
        BipolarSphericalHarmonics.monopolarharmonics!(B, pi/2, 0, pi/3, pi/4, ((1,2),(2,1)));
        for ((j,m), Y) in zip(first(SphericalHarmonicArrays.modes(Y1)), Y1)
            @test isapproxdefault(Y, vshbasis(PB(), Polar(), j, m, pi/2, 0))
        end
        for ((j,m), Y) in zip(first(SphericalHarmonicArrays.modes(Y2)), Y2)
            @test isapproxdefault(Y, vshbasis(PB(), Polar(), j, m, pi/3, pi/4))
        end
    end
end

@testset "NorthPole" begin
    θ1 = NorthPole(); ϕ1 = 0;
    l1 = 1; l2 = 1;
    B = BipolarSphericalHarmonics.cache(BipolarSphericalHarmonics.SH(), max(l1,l2))
    Y12 = zeros(ComplexF64, l1 + l2 + 1);

    m = -2;
    jrange = BipolarSphericalHarmonics.degreerange(l1, l2, m)
    modes = LM(jrange, m)
    Y = SHArray((@view Y12[1:length(modes)]), modes)
    for θ2 in LinRange(0, pi, 10), ϕ2 in LinRange(0, 2pi, 10)
        BipolarSphericalHarmonics.monopolarharmonics!(B, θ1, ϕ1, θ2, ϕ2)
        BipolarSphericalHarmonics.biposh!(Y12, B, θ1, ϕ1, θ2, ϕ2, :, m, l1, l2)
        @test all(iszero, Y)
    end

    m = -1;
    jrange = BipolarSphericalHarmonics.degreerange(l1, l2, m)
    modes = LM(jrange, m)
    Y = SHArray((@view Y12[1:length(modes)]), modes)
    for θ2 in LinRange(0, pi, 10), ϕ2 in LinRange(0, 2pi, 10)
        BipolarSphericalHarmonics.monopolarharmonics!(B, θ1, ϕ1, θ2, ϕ2)
        BipolarSphericalHarmonics.biposh!(Y12, B, θ1, ϕ1, θ2, ϕ2, :, m, l1, l2)
        Y2 = _biposh(m, l1, l2, θ1, ϕ1, θ2, ϕ2, B)
        @test isapproxdefault(Y[(1,-1)], 3/8pi*sin(θ2)*cis(-ϕ2))
        @test isapproxdefault(Y[(2,-1)], 3/8pi*sin(θ2)*cis(-ϕ2))
        @test isapproxdefault(Y, Y2)

        BipolarSphericalHarmonics.monopolarharmonics!(B, θ2, ϕ2, θ1, ϕ1)
        BipolarSphericalHarmonics.biposh!(Y12, B, θ2, ϕ2, θ1, ϕ1, :, m, l1, l2)
        Y2 = _biposh(m, l1, l2, θ2, ϕ2, θ1, ϕ1, B)
        @test isapproxdefault(Y, Y2)
    end

    m = 0;
    jrange = BipolarSphericalHarmonics.degreerange(l1, l2, m)
    modes = LM(jrange, m)
    Y = SHArray((@view Y12[1:length(modes)]), modes)
    for θ2 in LinRange(0, pi, 10), ϕ2 in LinRange(0, 2pi, 10)
        BipolarSphericalHarmonics.monopolarharmonics!(B, θ1, ϕ1, θ2, ϕ2)
        BipolarSphericalHarmonics.biposh!(Y12, B, θ1, ϕ1, θ2, ϕ2, :, m, l1, l2)
        Y2 = _biposh(m, l1, l2, θ1, ϕ1, θ2, ϕ2, B)
        @test isapproxdefault(Y[(0,0)], -√3/4pi * cos(θ2))
        @test isapproxdefault(Y[(1,0)], 0)
        @test isapproxdefault(Y[(2,0)], √(3/2)/2pi * cos(θ2))
        @test isapproxdefault(Y, Y2)

        BipolarSphericalHarmonics.monopolarharmonics!(B, θ2, ϕ2, θ1, ϕ1)
        BipolarSphericalHarmonics.biposh!(Y12, B, θ2, ϕ2, θ1, ϕ1, :, m, l1, l2)
        Y2 = _biposh(m, l1, l2, θ2, ϕ2, θ1, ϕ1, B)
        @test isapproxdefault(Y, Y2)
    end

    m = 1;
    jrange = BipolarSphericalHarmonics.degreerange(l1, l2, m)
    modes = LM(jrange, m)
    Y = SHArray((@view Y12[1:length(modes)]), modes)
    for θ2 in LinRange(0, pi, 10), ϕ2 in LinRange(0, 2pi, 10)
        BipolarSphericalHarmonics.monopolarharmonics!(B, θ1, ϕ1, θ2, ϕ2)
        BipolarSphericalHarmonics.biposh!(Y12, B, θ1, ϕ1, θ2, ϕ2, :, m, l1, l2)
        Y2 = _biposh(m, l1, l2, θ1, ϕ1, θ2, ϕ2, B)
        @test isapproxdefault(Y[(1,1)], 3/8pi*sin(θ2)*cis(ϕ2))
        @test isapproxdefault(Y[(2,1)], -3/8pi*sin(θ2)*cis(ϕ2))
        @test isapproxdefault(Y, Y2)

        BipolarSphericalHarmonics.monopolarharmonics!(B, θ2, ϕ2, θ1, ϕ1)
        BipolarSphericalHarmonics.biposh!(Y12, B, θ2, ϕ2, θ1, ϕ1, :, m, l1, l2)
        Y2 = _biposh(m, l1, l2, θ2, ϕ2, θ1, ϕ1, B)
        @test isapproxdefault(Y, Y2)
    end

    m = 2;
    jrange = BipolarSphericalHarmonics.degreerange(l1, l2, m)
    modes = LM(jrange, m)
    Y = SHArray((@view Y12[1:length(modes)]), modes)
    for θ2 in LinRange(0, pi, 10), ϕ2 in LinRange(0, 2pi, 10)
        BipolarSphericalHarmonics.monopolarharmonics!(B, θ1, ϕ1, θ2, ϕ2)
        BipolarSphericalHarmonics.biposh!(Y12, B, θ1, ϕ1, θ2, ϕ2, :, m, l1, l2)
        @test all(iszero, Y)
    end

    # check NorthPole and 0
    lmax = 5
    B = BipolarSphericalHarmonics.cache(BipolarSphericalHarmonics.SH(), lmax)
    for l1 in 0:lmax, l2 in 0:lmax
        BipolarSphericalHarmonics.monopolarharmonics!(B, 0, 0, 0, 0)
        B1 = biposh(B, 0, 0, 0, 0, :, :, l1, l2)
        BipolarSphericalHarmonics.monopolarharmonics!(B, NorthPole(), 0, 0, 0)
        B2 = biposh(B, NorthPole(), 0, 0, 0, :, :, l1, l2)
        BipolarSphericalHarmonics.monopolarharmonics!(B, 0, 0, NorthPole(), 0)
        B3 = biposh(B, 0, 0, NorthPole(), 0, :, :, l1, l2)
        BipolarSphericalHarmonics.monopolarharmonics!(B, NorthPole(), 0, NorthPole(), 0)
        B4 = biposh(B, NorthPole(), 0, NorthPole(), 0, :, :, l1, l2)
        @test isapproxdefault(B1, B2)
        @test isapproxdefault(B1, B3)
        @test isapproxdefault(B1, B4)
    end
end

@testset "l1l2modes" begin
    l1range = 1:2; l2range = 1:2
    θ1, ϕ1 = pi/3, pi/4
    θ2, ϕ2 = pi/2, pi/3
    for l1l2modes = Any[Iterators.product(l1range, l2range),
            L2L1Triangle(0:1, 1, 0:1), L1L2Triangle(0:1, 1, 0:1)]
        Y = biposh(SH(), θ1, ϕ1, θ2, ϕ2, :, :, l1l2modes)
        for (ind, (l1, l2)) in enumerate(l1l2modes)
            Y2 = biposh(SH(), θ1, ϕ1, θ2, ϕ2, :, :, l1, l2)
            @test isapproxdefault(Y[ind], Y2)
        end
        Y = biposh(SH(), θ1, ϕ1, θ2, ϕ2, :, 0, l1l2modes)
        for (ind, (l1, l2)) in enumerate(l1l2modes)
            Y2 = biposh(SH(), θ1, ϕ1, θ2, ϕ2, :, 0, l1, l2)
            @test isapproxdefault(Y[ind], Y2)
        end
    end

    l1range = 2:2; l2range = 2:2
    jrange = 0:4
    θ1, ϕ1 = pi/3, pi/4
    θ2, ϕ2 = pi/2, pi/3
    l1l2modes = Iterators.product(l1range, l2range)
    Y = biposh(SH(), θ1, ϕ1, θ2, ϕ2, jrange, :, l1l2modes)
    for (ind, (l1, l2)) in enumerate(l1l2modes)
        Y2 = biposh(SH(), θ1, ϕ1, θ2, ϕ2, jrange, :, l1, l2)
        @test isapproxdefault(Y[ind], Y2)
    end
    Y = biposh(SH(), θ1, ϕ1, θ2, ϕ2, jrange, 0, l1l2modes)
    for (ind, (l1, l2)) in enumerate(l1l2modes)
        Y2 = biposh(SH(), θ1, ϕ1, θ2, ϕ2, jrange, 0, l1, l2)
        @test isapproxdefault(Y[ind], Y2)
    end

    jₒjₛ_allmodes = L2L1Triangle(1:50, 2, 1:52)
    Y12 = biposh(SH(), θ1, ϕ1, θ2, ϕ2, 2, 0, jₒjₛ_allmodes)
    for (ind,(j2,j1)) in enumerate(jₒjₛ_allmodes)
        @test Y12[ind] ≈ biposh(SH(), θ1, ϕ1, θ2, ϕ2, 2, 0, j2, j1)
    end

    jₒjₛ_allmodes = L1L2Triangle(1:50, 2, 1:52)
    Y12 = biposh(SH(), θ1, ϕ1, θ2, ϕ2, 2, 0, jₒjₛ_allmodes)
    for (ind,(j1,j2)) in enumerate(jₒjₛ_allmodes)
        @test Y12[ind] ≈ biposh(SH(), θ1, ϕ1, θ2, ϕ2, 2, 0, j1, j2)
    end
end

@testset "Y^{11}_{00}" begin
    θ1, ϕ1 = pi*rand(), 2pi*rand()
    θ2, ϕ2 = pi*rand(), 2pi*rand()

    @testset "SH" begin
        B = biposh(SH(), θ1, ϕ1, θ2, ϕ2, 0, 0, 1, 1)
        @test B ≈ -√3/4π * cosχ(θ1, ϕ1, θ2, ϕ2)
    end

    @testset "GSH" begin
        B = biposh(GSH(), θ1, ϕ1, θ2, ϕ2, 0, 0, 1, 1)
        @test B[kronindex(GSH(), 0, 0)] ≈ -√3/4π * cosχ(θ1, ϕ1, θ2, ϕ2)
        for α1 = -1:1, α2 = -1:1
            @test isapprox(B[kronindex(GSH(), -α1, -α2)], (-1)^(α1 + α2) * conj(B[kronindex(GSH(), α1, α2)]),
                atol = 1e-14, rtol = 1e-8)
        end

        # We check the values at the Equator
        B = biposh(GSH(), pi/2, pi/4, pi/2, pi/2, 0, 0, 1, 1)
        B_exp = vec(transpose(
                [√3/8π*(-1 + cosh(im*(-pi/4)))    √(3/2)/4π*sinh(im*(-pi/4))     √3/8π*(1 + cosh(im*(-pi/4)))
                -√(3/2)/4π*sinh(im*(-pi/4))       -√3/4π*cosh(im*(-pi/4))       -√(3/2)/4π*sinh(im*(-pi/4))
                √3/8π*(1 + cosh(im*(-pi/4)))     √(3/2)/4π*sinh(im*(-pi/4))      √3/8π*(-1 + cosh(im*(-pi/4)))]))

        @test isapproxdefault(B, B_exp)
    end
end

@testset "Y^{ll}_{00}" begin
    θ1, ϕ1 = pi*rand(), 2pi*rand()
    θ2, ϕ2 = pi*rand(), 2pi*rand()

    lmax = 10
    P = collectPl(cosχ(θ1, ϕ1, θ2, ϕ2), lmax = lmax)
    Yll_00 = [P[l]*(-1)^l * √(2l+1)/4π for l in 1:lmax]
    @testset "SH" begin
        YB_00 = [biposh(SH(), θ1, ϕ1, θ2, ϕ2, 0, 0, l, l) for l in 1:lmax]
        @test Yll_00 ≈ YB_00
    end

    @testset "GSH" begin
        for l in 1:lmax
            B = biposh(GSH(), θ1, ϕ1, θ2, ϕ2, 0, 0, l, l)
            for α1 = -1:1, α2 = -1:1
                @test isapprox(B[kronindex(GSH(), -α1, -α2)], (-1)^(α1 + α2) * conj(B[kronindex(GSH(), α1, α2)]),
                    atol = 1e-14, rtol = 1e-8)
            end
            @test B[kronindex(GSH(), 0, 0)] ≈ Yll_00[l]
        end
    end
end

@testset "Y^{j1j1}_{10}" begin
    θ1, ϕ1 = pi*rand(), 2pi*rand()
    θ2, ϕ2 = pi*rand(), 2pi*rand()

    lmax = 10
    dP = collectdnPl(cosχ(θ1, ϕ1, θ2, ϕ2), lmax = lmax, n = 1)
    Yll_10 = [dP[l]*im*(-1)^l * √(3*(2l+1)/(l*(l+1)))/4π * ∂ϕ₂cosχ(θ1, ϕ1, θ2, ϕ2) for l in 1:lmax]

    @testset "SH" begin
        YB_10_n1n2 = zeros(ComplexF64, 1:lmax)
        YB_10_n2n1 = similar(YB_10_n1n2)
        for l in 1:lmax
            YB_10_n1n2[l] = biposh(SH(), θ1, ϕ1, θ2, ϕ2, 1, 0, l, l)
            YB_10_n2n1[l] = biposh(SH(), θ2, ϕ2, θ1, ϕ1, 1, 0, l, l)
        end

        @test Yll_10 ≈ YB_10_n1n2
        @test YB_10_n1n2 ≈ -YB_10_n2n1

        @test biposh(SH(), θ1, ϕ1, θ2, ϕ2, LM(1, 0), 1, 1)[1] == biposh(SH(), θ1, ϕ1, θ2, ϕ2, 1, 0, 1, 1)
    end

    @testset "GSH" begin
        for l in 1:lmax
            B = biposh(GSH(), θ1, ϕ1, θ2, ϕ2, 1, 0, l, l)
            for α1 = -1:1, α2 = -1:1
                @test isapprox(B[kronindex(GSH(), -α1, -α2)], (-1)^(α1 + α2 + 1) * conj(B[kronindex(GSH(), α1, α2)]),
                    atol = 1e-14, rtol = 1e-8)
            end
            @test B[kronindex(GSH(), 0, 0)] ≈ Yll_10[l]
        end

        # We check the values for specific j1
        B = biposh(GSH(), pi/2, pi/4, pi/2, pi/2, 1, 0, 1, 1)
        B_exp = vec(transpose(
                [3/(8*√2π)*sinh(im*(-pi/4))    3/8π*cosh(im*(-pi/4))     3/(8*√2π)*sinh(im*(-pi/4))
                -3/8π*cosh(im*(-pi/4))       -3/(4*√2π)*sinh(im*(-pi/4))      -3/8π*cosh(im*(-pi/4))
                3/(8*√2π)*sinh(im*(-pi/4))     3/8π*cosh(im*(-pi/4))      3/(8*√2π)*sinh(im*(-pi/4))]))

        @test isapproxdefault(B, B_exp)
    end
end

# test symmetries
@testset "Y^{j1j1}_{j0}" begin
    θ1, ϕ1 = pi*rand(), 2pi*rand()
    θ2, ϕ2 = pi*rand(), 2pi*rand()

    lmax = 5; jmax = 5
    @testset "SH" begin
        for j in 2:jmax, l in div(j,2)+1:lmax
            B = biposh(SH(), θ1, ϕ1, θ2, ϕ2, j, 0, l, l)
            if iseven(j)
                @test isapprox(imag(B), 0, atol = 1e-14)
            else
                @test isapprox(real(B), 0, atol = 1e-14)
            end
        end
    end
    @testset "GSH" begin
        for j in 2:jmax, l in div(j,2)+1:lmax
            B = biposh(GSH(), θ1, ϕ1, θ2, ϕ2, j, 0, l, l)
            for α1 = -1:1, α2 = -1:1
                @test isapprox(B[kronindex(GSH(), -α1, -α2)], (-1)^(α1 + α2 + j) * conj(B[kronindex(GSH(), α1, α2)]),
                    atol = 1e-14, rtol = 1e-8)
            end
            if iseven(j)
                @test isapprox(imag(B[kronindex(GSH(), 0, 0)]), 0, atol = 1e-14)
            else
                @test isapprox(real(B[kronindex(GSH(), 0, 0)]), 0, atol = 1e-14)
            end
        end

        # test values
        B = biposh(GSH(), pi/2, pi/4, pi/2, pi/2, 2, 0, 1, 1)
        B_exp = vec(transpose(
                [√(3/2)/8π*(2 + cosh(im*(-pi/4)))    √3/8π*sinh(im*(-pi/4))       √(3/2)/8π*(-2 + cosh(im*(-pi/4)))
                 -√3/8π*sinh(im*(-pi/4))            -√(3/2)/4π*cosh(im*(-pi/4))   -√3/8π*sinh(im*(-pi/4))
                 √(3/2)/8π*(-2 + cosh(im*(-pi/4)))   √3/8π*sinh(im*(-pi/4))       √(3/2)/8π*(2 + cosh(im*(-pi/4))) ]))

        @test isapproxdefault(B, B_exp)
    end
end

@testset "Y^{j1j2}_{jm}" begin
    θ1, ϕ1 = pi*rand(), 2pi*rand()
    θ2, ϕ2 = pi*rand(), 2pi*rand()

    lmax = 2;

    function testmultiplemodes(C, θ1, ϕ1, θ2, ϕ2, lmax)
        for l1 in 0:lmax, l2 in 0:lmax
            B = biposh(C, θ1, ϕ1, θ2, ϕ2, :, :, l1, l2)
            if l1 + l2 > 0
                Bm = biposh(C, θ1, ϕ1, θ2, ϕ2, :, 1:(l1+l2), l1, l2)
            end
            for j in abs(l1 - l2):l1 + l2
                Bj = biposh(C, θ1, ϕ1, θ2, ϕ2, j, :, l1, l2)
                Bjneg = biposh(C, θ1, ϕ1, θ2, ϕ2, j, -j:-min(1,j), l1, l2)
                Bjpos = biposh(C, θ1, ϕ1, θ2, ϕ2, j, min(1,j):j, l1, l2)
                for m in -j:j
                    Bjm = biposh(C, θ1, ϕ1, θ2, ϕ2, j, m, l1, l2)
                    @test B[(j,m)] ≈ Bjm
                    @test Bj[(j,m)] ≈ Bjm
                    B2jm = biposh(C, θ1, ϕ1, θ2, ϕ2, j, m, ((l1, l2),))
                    @test B2jm[1] ≈ Bjm

                    B3jm = biposh(C, θ1, ϕ1, θ2, ϕ2, LM(j, m), l1, l2)
                    @test B3jm[(j,m)] ≈ Bjm
                    B4jm = biposh(C, θ1, ϕ1, θ2, ϕ2, LM(j, m), ((l1, l2),))
                    @test B4jm[1][(j,m)] ≈ Bjm

                    if m < 0
                        @test isapproxdefault(Bjneg[(j,m)], Bjm)
                    elseif m > 0
                        @test isapproxdefault(Bjpos[(j,m)], Bjm)
                        @test isapproxdefault(Bm[(j,m)], Bjm)
                    end
                end
            end
        end
    end

    @testset "SH" begin
        C = monopolarharmonics(SH(), θ1, ϕ1, θ2, ϕ2, lmax);

        @test biposh(C, θ1, ϕ1, θ2, ϕ2, 3, 3, 1, 1) === nothing

        @testset "multiple vs one" begin
            testmultiplemodes(C, θ1, ϕ1, θ2, ϕ2, lmax)
        end
    end

    @testset "GSH" begin
        C = monopolarharmonics(GSH(), θ1, ϕ1, θ2, ϕ2, lmax);
        for l1 in 0:lmax, l2 in 0:lmax, j in abs(l1 - l2):l1 + l2, m in 0:j
            Bm = biposh(C, θ1, ϕ1, θ2, ϕ2, j, m, l1, l2)
            Bminm = biposh(C, θ1, ϕ1, θ2, ϕ2, j, -m, l1, l2)
            for α1 = -1:1, α2 = -1:1
                @test isapprox(Bm[kronindex(GSH(), -α1, -α2)], (-1)^(α1 + α2 + j + m + l1 + l2) * conj(Bminm[kronindex(GSH(), α1, α2)]),
                    atol = 1e-14, rtol = 1e-8)
            end
        end

        @test biposh(C, θ1, ϕ1, θ2, ϕ2, 3, 3, 1, 1) === nothing

        @testset "multiple vs one" begin
            testmultiplemodes(C, θ1, ϕ1, θ2, ϕ2, lmax)
        end
    end

    @testset "VSH" begin
        @testset "PB" begin
            CGSH = monopolarharmonics(GSH(), θ1, ϕ1, θ2, ϕ2, lmax);
            @testset "PB HelicityCovariant" begin
                CVSH = monopolarharmonics(VSH(PB(), HelicityCovariant()), θ1, ϕ1, θ2, ϕ2, lmax);
                for l1 in 0:lmax, l2 in 0:lmax
                    G = biposh(CGSH, θ1, ϕ1, θ2, ϕ2, :, :, l1, l2)
                    V = biposh(CVSH, θ1, ϕ1, θ2, ϕ2, :, :, l1, l2)
                    for j in abs(l1 - l2):l1 + l2, m in -j:j
                        @test diag(V[(j,m)]) ≈ G[(j,m)]
                    end
                end
            end
        end
        for YT in [PB(), Hansen(), Irreducible()], B in [SphericalCovariant(), HelicityCovariant(), Polar(), Cartesian()]
            @testset "$(typeof(YT)) $(typeof(B))" begin
                CVSH = monopolarharmonics(VSH(YT, B), θ1, ϕ1, θ2, ϕ2, lmax);
                testmultiplemodes(CVSH, θ1, ϕ1, θ2, ϕ2, lmax)
                @test biposh(CVSH, θ1, ϕ1, θ2, ϕ2, 3, 3, 1, 1) === nothing
            end
        end
    end
end

@testset "swap points" begin
    θ1, ϕ1 = pi*rand(), 2pi*rand()
    θ2, ϕ2 = pi*rand(), 2pi*rand()

    jmax = 2;

    @testset "biposh" begin
        function testswap(SHT, θ1, ϕ1, θ2, ϕ2, jmax, perm)
            C12 = monopolarharmonics(SHT, θ1, ϕ1, θ2, ϕ2, jmax);
            C21 = monopolarharmonics(SHT, θ2, ϕ2, θ1, ϕ1, jmax);
            testswap(C12, C21, θ1, ϕ1, θ2, ϕ2, jmax, perm)
        end

        function testswap(C12, C21, θ1, ϕ1, θ2, ϕ2, jmax, perm)
            for j1 in 0:jmax, j2 in 0:jmax
                for j in abs(j1 - j2):j1+j2, m in -j:j
                    B12 = biposh(C12, θ1, ϕ1, θ2, ϕ2, j, m, j1, j2);
                    B21 = biposh(C21, θ2, ϕ2, θ1, ϕ1, j, m, j2, j1);
                    @test begin
                        res = isapproxdefault(B21, (-1)^(j1 + j2 - j) * B12[perm])
                        if !res
                            @show j1, j2, j, m
                        end
                        res
                    end
                end
            end
        end

        @testset "SH" begin
            tstart = time()
            testswap(SH(), θ1, ϕ1, θ2, ϕ2, jmax, 1)
            tend = time()
            @info "Finished testing SH() in $(round(tend - tstart, sigdigits = 1)) seconds"
        end

        @testset "GSH" begin
            tstart = time()
            testswap(GSH(), θ1, ϕ1, θ2, ϕ2, jmax, GSHPerm)
            tend = time()
            @info "Finished testing GSH() in $(round(tend - tstart, sigdigits = 1)) seconds"
        end

        @testset "VSH" begin
            for YT in [PB(), Hansen(), Irreducible()], B in [Polar(), Cartesian(), HelicityCovariant(), SphericalCovariant()]
                @testset "$(typeof(YT)) $(typeof(B))" begin
                    tstart = time()
                    testswap(VSH(YT, B), θ1, ϕ1, θ2, ϕ2, jmax, VSHPerm)
                    tend = time()
                    @info "Finished testing $YT $B in $(round(tend - tstart, sigdigits = 1)) seconds"
                end
            end
        end
    end

    @testset "biposh_flippoints" begin
        j1j2modes = L2L1Triangle(1:jmax, jmax, 1:jmax)

        function testswap(SHT, θ1, ϕ1, θ2, ϕ2, j1j2modes)
            B12 = monopolarharmonics(SHT, θ1, ϕ1, θ2, ϕ2, j1j2modes);
            B21 = monopolarharmonics(SHT, θ2, ϕ2, θ1, ϕ1, j1j2modes);
            Y12_alljm, Y21_alljm = biposh_flippoints(B12, B21, θ1, ϕ1, θ2, ϕ2, :, :, j1j2modes);
            for (indj1j2, (j1, j2)) in enumerate(j1j2modes)
                Y12_j1j2_alljm = biposh(B12, θ1, ϕ1, θ2, ϕ2, :, :, j1, j2)
                @test isapproxdefault(Y12_alljm[indj1j2], Y12_j1j2_alljm)
                Y21_j1j2_alljm = biposh(B21, θ2, ϕ2, θ1, ϕ1, :, :, j1, j2)
                @test isapproxdefault(Y21_alljm[indj1j2], Y21_j1j2_alljm)
            end
            j, m = 1, 1
            Y12_jm, Y21_jm = biposh_flippoints(B12, B21, θ1, ϕ1, θ2, ϕ2, j, m, j1j2modes);
            for (indj1j2, (j1, j2)) in enumerate(j1j2modes)
                Y12_j1j2_jm = biposh(B12, θ1, ϕ1, θ2, ϕ2, j, m, j1, j2)
                if Y12_j1j2_jm === nothing
                    @test Y12_jm[indj1j2] === nothing
                else
                    @test isapproxdefault(Y12_jm[indj1j2], Y12_j1j2_jm)
                end
                Y21_j1j2_jm = biposh(B21, θ2, ϕ2, θ1, ϕ1, j, m, j1, j2)
                if Y21_j1j2_jm === nothing
                    @test Y21_jm[indj1j2] === nothing
                else
                    @test isapproxdefault(Y21_jm[indj1j2], Y21_j1j2_jm)
                end
            end

            j1, j2 = 2, 2
            Y12, Y21 = biposh_flippoints(B12, B21, θ1, ϕ1, θ2, ϕ2, j, m, j1, j2);
            Y12_j1j2 = biposh(B12, θ1, ϕ1, θ2, ϕ2, j, m, j1, j2)
            Y21_j1j2 = biposh(B21, θ2, ϕ2, θ1, ϕ1, j, m, j1, j2)
            @test isapproxdefault(Y12, Y12_j1j2)
            @test isapproxdefault(Y21, Y21_j1j2)

            @test biposh_flippoints(B12, B21, θ1, ϕ1, θ2, ϕ2, 3, 3, 1, 1) === nothing
            @test biposh_flippoints(B12, B21, θ1, ϕ1, θ2, ϕ2, LM(3, 3), 1, 1) === nothing
            @test biposh_flippoints(B12, B21, θ1, ϕ1, θ2, ϕ2, 3, 3, L2L1Triangle(1:1, 0)) == ([nothing], [nothing])
            @test biposh_flippoints(B12, B21, θ1, ϕ1, θ2, ϕ2, LM(3, 3), L2L1Triangle(1:1, 0)) == ([nothing], [nothing])
            return nothing
        end

        for SHT in [SH(), GSH()]
            @testset "$(typeof(SHT))" begin
                tstart = time()
                testswap(SHT, θ1, ϕ1, θ2, ϕ2, j1j2modes)
                Y12_jmj1j2 = biposh_flippoints(SHT, θ1, ϕ1, θ2, ϕ2, 1, 1, 1, 1)[1]
                Y12_jmj1j2modes = biposh_flippoints(SHT, θ1, ϕ1, θ2, ϕ2, LM(1, 1), 1, 1)[1][1]
                @test isapproxdefault(Y12_jmj1j2, Y12_jmj1j2)
                tend = time()
                @info "Finished testing $SHT in $(round(tend - tstart, sigdigits = 1)) seconds"
            end
        end
        @testset "VSH" begin
            for YT in [PB(), Hansen(), Irreducible()], B in [Polar(), Cartesian(), HelicityCovariant(), SphericalCovariant()]
                @testset "$(typeof(YT)) $(typeof(B))" begin
                    tstart = time()
                    testswap(VSH(YT, B), θ1, ϕ1, θ2, ϕ2, j1j2modes)
                    tend = time()
                    @info "Finished testing $YT $B in $(round(tend - tstart, sigdigits = 1)) seconds"
                end
            end
        end
    end
end

@testset "Rotation of coordinates" begin
    # define some helper functions to rotate the points
    cartvec(θ, ϕ) = SVector{3}(sin(θ)cos(ϕ), sin(θ)sin(ϕ), cos(θ))
    function polcoords(v)
        # This is not unique at the poles, as (θ=0, ϕ) map to the north pole for all ϕ
        vunit = normalize(v)
        if vunit[3] == 1
            # north pole
            return SVector{2}(promote(oftype(float(vunit[3]), 0), 0))
        elseif vunit[3] == -1
            # south pole
            return SVector{2}(promote(oftype(float(vunit[3]), pi), 0))
        end
        θ = acos(vunit[3])
        ϕ = mod2pi(atan(vunit[2], vunit[1]))
        SVector{2}(promote(θ, ϕ))
    end
    for θ1 in LinRange(0, pi, 10), ϕ1 in LinRange(0, 2pi, 10)
        @test polcoords(cartvec(θ1, ϕ1)) ≈ SVector{2}(θ1, abs(cos(θ1)) == 1 ? 0 : ϕ1)
    end
    # A rotation that takes (θ1, ϕ1) to (θ1′, ϕ1′). The rotation is not unique
    rotn1n1′(θ1, ϕ1, θ1′, ϕ1′) = RotZYZ(ϕ1′, θ1′ - θ1, -ϕ1)
    for θ1 in LinRange(0, pi, 4), ϕ1 in LinRange(0, 2pi, 4), θ1′ in LinRange(0, pi, 4), ϕ1′ in LinRange(0, 2pi, 4)
        R = rotn1n1′(θ1, ϕ1, θ1′, ϕ1′)
        @test det(R) ≈ 1
        @test R'R ≈ I
        @test R * cartvec(θ1, ϕ1) ≈ cartvec(θ1′, ϕ1′)
    end
    # A rotation about the axis n1′ that takes n2′′ to n2′
    function rotn1′n2n2′(n1′, n2′′, n2′)
        n1′x, n1′y, n1′z = n1′
        a = n2′′ - n1′
        b = n2′ - n1′
        if (isapproxdefault(cross(n1′, a), zero(SVector{3})) || isapproxdefault(cross(n1′, b), zero(SVector{3})))
            return one(SMatrix{3,3,eltype(n1′)})
        end
        sinω = b ⋅ cross(n1′, a) / norm(cross(n1′, a))^2
        cosω = b ⋅ (a - (n1′⋅a)n1′) / norm(a - (n1′⋅a)n1′)^2
        n1′C = SMatrix{3,3}(0, n1′z, -n1′y, -n1′z, 0, n1′x, n1′y, -n1′x, 0)
        n1′n1′ = n1′ * n1′'
        R = cosω*I + sinω * n1′C + (1 - cosω)*n1′n1′
    end
    for Δθ in LinRange(0, pi, 4), ϕ1 in LinRange(0, 2pi, 4), ϕ2 in LinRange(0, 2pi, 4)
        n1′ = cartvec(0, 0)
        n2′′ = cartvec(Δθ, ϕ1)
        n2′ = cartvec(Δθ, ϕ2)
        R = rotn1′n2n2′(n1′, n2′′, n2′)
        @test det(R) ≈ 1
        @test R'R ≈ I
        @test R * n1′ ≈ n1′
        @test R * n2′′ ≈ n2′
    end
    # a rotation that takes (θ1, ϕ1) and (θ2, ϕ2) to (θ1′, ϕ1′) and (θ2′, ϕ2′) for pairs of points that are equispaced in angle
    function rotn1n2n1′n2′(θ1, ϕ1, θ2, ϕ2, θ1′, ϕ1′, θ2′, ϕ2′)
        local n1, n2
        n1 = cartvec(θ1, ϕ1)
        n2 = cartvec(θ2, ϕ2)
        n1′ = cartvec(θ1′, ϕ1′)
        n2′ = cartvec(θ2′, ϕ2′)
        n1 ⋅ n2 ≈ n1′ ⋅ n2′ || throw(ArgumentError("Points are not related by a rotation"))
        R11′ = rotn1n1′(θ1, ϕ1, θ1′, ϕ1′)
        n2′′ = R11′ * n2
        R2′′2′ = rotn1′n2n2′(n1′, n2′′, n2′)
        R = RotZYZ(R2′′2′ * R11′)
    end
    for θ1 in LinRange(0, pi, 4), θ2 in LinRange(0, pi, 4), ϕ1′ in LinRange(0, 2pi, 4)
        Δθ21 = θ2 - θ1
        ϕ1 = ϕ2 = 0.0
        n1 = cartvec(θ1, ϕ1)
        n2 = cartvec(θ2, ϕ2)

        θ1′ = θ2′ = pi/2
        Δϕ2′1′ = Δθ21
        ϕ2′ = ϕ1′ + Δϕ2′1′
        n1′ = cartvec(θ1′, ϕ1′)
        n2′ = cartvec(θ2′, ϕ2′)
        @test dot(n1, n2) ≈ dot(n1′, n2′)

        R = rotn1n2n1′n2′(θ1, ϕ1, θ2, ϕ2, θ1′, ϕ1′, θ2′, ϕ2′)
        @test det(R) ≈ 1
        @test R'R ≈ I
        @test R * n1 ≈ n1′
        @test R * n2 ≈ n2′

        # Test by rotating the points n2 and n2′
        # We rotate n2 about n1 and n2′ about n1′
        # This preserves the angular distance between the pairs
        for ω1 in LinRange(0, 2pi, 6), ω2 in LinRange(0, 2pi, 6)
            n2_rot = AngleAxis(ω1, n1...) * n2
            θ2_rot, ϕ2_rot = polcoords(n2_rot)
            @test dot(n1, n2_rot) ≈ dot(n1, n2)

            n2′_rot = AngleAxis(ω2, n1′...) * n2′
            @test dot(n1′, n2′_rot) ≈ dot(n1′, n2′)
            θ2′_rot, ϕ2′_rot = polcoords(n2′_rot)

            R = rotn1n2n1′n2′(θ1, ϕ1, θ2_rot, ϕ2_rot, θ1′, ϕ1′, θ2′_rot, ϕ2′_rot)
            @test det(R) ≈ 1
            @test R'R ≈ I
            @test R * n1 ≈ n1′
            @test R * n2_rot ≈ n2′_rot
        end
    end

    lmax = 5
    modes = ML(0:lmax)
    Dvec = OffsetArray([zeros(ComplexF64, 2l+1, 2l+1) for l in 0:2lmax], 0:2lmax);

    θ1, ϕ1 = pi/4, 0; n1 = cartvec(θ1, ϕ1);
    θ2, ϕ2 = pi/3, pi/3; n2 = cartvec(θ2, ϕ2);

    @testset "SH" begin
        B12 = cache(SH(), lmax);
        B1′2′ = cache(SH(), lmax);
        BipolarSphericalHarmonics.monopolarharmonics!(B12, θ1, ϕ1, θ2, ϕ2);
        for θ1′ in LinRange(0, pi, 6), ϕ1′ in LinRange(0, 2pi, 6)
            # start by choosing the rotation that maps n1 to n1′
            R11′ = rotn1n1′(θ1, ϕ1, θ1′, ϕ1′)
            n2′ = R11′ * cartvec(θ2, ϕ2)
            θ2′, ϕ2′ = polcoords(n2′)
            n1′ = cartvec(θ1′, ϕ1′)
            BipolarSphericalHarmonics.monopolarharmonics!(B1′2′, θ1′, ϕ1′, θ2′, ϕ2′);
            @test R11′ * n1 ≈ n1′
            @test R11′ * n2 ≈ n2′
            for j in 0:2lmax
                Dp = wignerD!(Dvec[j], j, ϕ1, θ1 - θ1′, -ϕ1′);
                D = OffsetArray(Dp, -j:j, -j:j)
                for l1 in 0:lmax, l2 in abs(j - l1):min(lmax, l1 + j)
                    Y  = biposh(B12, θ1, ϕ1, θ2, ϕ2, j, :, l1, l2)
                    Y′ = biposh(B1′2′, θ1′, ϕ1′, θ2′, ϕ2′, j, :, l1, l2)
                    Y_rot = parent([sum(D[m′, m] * Y[(j, m′)] for m′ in -j:j) for m in axes(D, 2)])
                    @test isapproxdefault(Y′, Y_rot)
                end
            end

            # next rotate (θ2′,ϕ2′) to a new position
            for ω2 in LinRange(0, 2pi, 4)[2:end]
                n2_rot = AngleAxis(ω2, n1...) * n2
                θ2_rot, ϕ2_rot = polcoords(n2_rot)
                @test dot(n1, n2_rot) ≈ dot(n1, n2)
                BipolarSphericalHarmonics.monopolarharmonics!(B12, θ1, ϕ1, θ2_rot, ϕ2_rot);

                for ω2′ in LinRange(0, 2pi, 4)[2:end]
                    n2′_rot = AngleAxis(ω2′, n1′...) * n2′
                    θ2′_rot, ϕ2′_rot = polcoords(n2′_rot)
                    @test dot(n1′, n2′_rot) ≈ dot(n1′, n2′) ≈ dot(n1, n2)

                    BipolarSphericalHarmonics.monopolarharmonics!(B1′2′, θ1′, ϕ1′, θ2′_rot, ϕ2′_rot);
                    R = rotn1n2n1′n2′(θ1, ϕ1, θ2_rot, ϕ2_rot, θ1′, ϕ1′, θ2′_rot, ϕ2′_rot)

                    @test det(R) ≈ 1
                    @test R'R ≈ I
                    @test R * n1 ≈ n1′
                    @test R * n2_rot ≈ n2′_rot
                    # The passive rotation for the frame is the inverse of R11′
                    RSS′ = inv(R)
                    α, β, γ = RSS′.theta1, RSS′.theta2, RSS′.theta3
                    for j in 0:lmax
                        # The rotation that transforms n1 to n1′ is the inverse of the one that transforms the harmoncics
                        Dp = wignerD!(Dvec[j], j, α, β, γ);
                        D = OffsetArray(Dp, -j:j, -j:j)
                        for l1 in 0:lmax, l2 in abs(j - l1):min(lmax, l1 + j)
                            Y  = biposh(B12, θ1, ϕ1, θ2_rot, ϕ2_rot, j, :, l1, l2)
                            Y′ = biposh(B1′2′, θ1′, ϕ1′, θ2′_rot, ϕ2′_rot, j, :, l1, l2)
                            Y_rot = parent([sum(D[m′, m] * Y[(j, m′)] for m′ in -j:j) for m in axes(D, 2)])
                            @test isapproxdefault(Y′, Y_rot)
                        end
                    end
                end
            end
        end
    end

    @testset "GSH" begin
        B12 = cache(GSH(), lmax);
        B1′2′ = cache(GSH(), lmax);
        B12_rot = cache(GSH(), lmax);
        B1′2′_rot = cache(GSH(), lmax);
        BipolarSphericalHarmonics.monopolarharmonics!(B12, θ1, ϕ1, θ2, ϕ2);
        Un1 = basisconversionmatrix(Cartesian(), HelicityCovariant(), θ1, ϕ1);
        Un2 = basisconversionmatrix(Cartesian(), HelicityCovariant(), θ2, ϕ2);
        x̂_n1 = Un1 * SVector{3}(1,0,0)
        x̂_n2 = Un2 * SVector{3}(1,0,0)
        x̂x̂_n1n2 = kron(x̂_n1, x̂_n2)

        for θ1′ in LinRange(0, pi, 6), ϕ1′ in LinRange(0, 2pi, 6)
            # start by choosing the rotation that maps n1 to n1′
            n1′ = cartvec(θ1′, ϕ1′)
            R11′ = rotn1n1′(θ1, ϕ1, θ1′, ϕ1′)
            n2′ = R11′ * n2
            θ2′, ϕ2′ = polcoords(n2′)
            @test dot(n1, n2) ≈ dot(n1′, n2′)
            @test R11′ * n1 ≈ n1′
            @test R11′ * n2 ≈ n2′
            BipolarSphericalHarmonics.monopolarharmonics!(B1′2′, θ1′, ϕ1′, θ2′, ϕ2′);
            # The passive rotation for the frame is the inverse of R11′
            RSS′ = inv(R11′)
            α, β, γ = RSS′.theta1, RSS′.theta2, RSS′.theta3
            U′n1 = basisconversionmatrix(Cartesian(), HelicityCovariant(), θ1′, ϕ1′);
            U′n2 = basisconversionmatrix(Cartesian(), HelicityCovariant(), θ2′, ϕ2′);
            M′1 = U′n1 * R11′ * Un1'
            M′2 = U′n2 * R11′ * Un2'
            M′ = kron(M′1, M′2)
            # For HelicityCovariant both M′1 and M′2, and consequently their kronecker product, are diagonal
            @test M′ ≈ Diagonal(M′)

            x̂x̂_n1′n2′ = kron(M′1 * x̂_n1, M′2 * x̂_n2)
            R⁻¹x̂R⁻¹x̂_n1n2 = kron(inv(M′1) * x̂_n1, inv(M′2) * x̂_n2)

            for j in 0:2lmax
                Dp = wignerD!(Dvec[j], j, α, β, γ);
                D = OffsetArray(Dp, -j:j, -j:j);
                for l1 in 0:lmax, l2 in abs(j - l1):min(lmax, l1 + j)
                    for m in -j:j
                        DY = sum(D[m′, m] * biposh(B12, θ1, ϕ1, θ2, ϕ2, j, m′, l1, l2) for m′ in -j:j)
                        Y_rot = M′ * DY
                        Y′ = biposh(B1′2′, θ1′, ϕ1′, θ2′, ϕ2′, j, m, l1, l2)
                        @test isapproxdefault(Y′, Y_rot)
                        # x′x′ : Bjm′(n1′,n2′) = M1xM2x : M1M2 ∑Dlmm′ Bjm′(n1,n2) = xx : ∑Dlmm′ Bjm′(n1,n2)
                        xxY12_rot = x̂x̂_n1n2' * DY
                        xxY1′2′ = x̂x̂_n1′n2′' * Y′
                        @test isapproxdefault(xxY1′2′, xxY12_rot)
                        # xx : Bjm′(n1′,n2′) = xx : M1M2 ∑Dlmm′ Bjm′(n1,n2) = M1⁻¹xM2⁻¹x : ∑Dlmm′ Bjm′(n1,n2)
                        xxY12_rot = R⁻¹x̂R⁻¹x̂_n1n2' * DY
                        xxY1′2′ = x̂x̂_n1n2' * Y′
                        @test isapproxdefault(xxY1′2′, xxY12_rot)
                    end
                end
            end

            # next rotate (θ2,ϕ2) and (θ2′,ϕ2′) to new positions
            for ω2 in LinRange(0, 2pi, 4)[2:end]
                R = AngleAxis(ω2, n1...)
                n2_rot = R * n2
                θ2_rot, ϕ2_rot = polcoords(n2_rot)
                @test dot(n1, n2_rot) ≈ dot(n1, n2)
                BipolarSphericalHarmonics.monopolarharmonics!(B12_rot, θ1, ϕ1, θ2_rot, ϕ2_rot);

                RSS′ = inv(R) # passive rotation that maps n2 to n2_rot
                Un2_rot = basisconversionmatrix(Cartesian(), HelicityCovariant(), θ2_rot, ϕ2_rot);
                M22_rot = Un2_rot * RSS′' * Un2'

                x̂_n2_rot = M22_rot * SVector{3}(1,0,0)
                x̂x̂_n1n2_rot = kron(x̂_n1, x̂_n2_rot)

                for ω2′ in LinRange(0, 2pi, 4)[2:end]
                    n2′_rot = AngleAxis(ω2′, n1′...) * n2′
                    θ2′_rot, ϕ2′_rot = polcoords(n2′_rot)
                    @test dot(n1′, n2′_rot) ≈ dot(n1′, n2′) ≈ dot(n1, n2)
                    R = rotn1n2n1′n2′(θ1, ϕ1, θ2_rot, ϕ2_rot, θ1′, ϕ1′, θ2′_rot, ϕ2′_rot)

                    @test det(R) ≈ 1
                    @test R'R ≈ I
                    @test R * n1 ≈ n1′
                    @test R * n2_rot ≈ n2′_rot
                    RSS′ = inv(R)
                    α, β, γ = RSS′.theta1, RSS′.theta2, RSS′.theta3
                    U′n2_rot = basisconversionmatrix(Cartesian(), HelicityCovariant(), θ2′_rot, ϕ2′_rot);
                    M′1 = U′n1 * R * Un1'
                    M′2 = U′n2_rot * R * Un2_rot'
                    M′ = kron(M′1, M′2)
                    @test M′ ≈ Diagonal(M′)

                    x̂_n1′ = M′1 * x̂_n1
                    x̂_n2′_rot = M′2 * x̂_n2_rot
                    x̂x̂_n1′n2′_rot = kron(x̂_n1′, x̂_n2′_rot)

                    BipolarSphericalHarmonics.monopolarharmonics!(B1′2′_rot, θ1′, ϕ1′, θ2′_rot, ϕ2′_rot);

                    for j in 0:lmax
                        Dp = wignerD!(Dvec[j], j, α, β, γ);
                        D = OffsetArray(Dp, -j:j, -j:j)
                        for l1 in 0:lmax, l2 in abs(j - l1):min(lmax, l1 + j), m in -j:j
                            DY = sum(D[m′, m] * biposh(B12_rot, θ1, ϕ1, θ2_rot, ϕ2_rot, j, m′, l1, l2) for m′ in -j:j)
                            RY12_rot = M′ * DY
                            Y1′2′_rot = biposh(B1′2′_rot, θ1′, ϕ1′, θ2′_rot, ϕ2′_rot, j, m, l1, l2)
                            @test isapproxdefault(Y1′2′_rot, RY12_rot)
                            # x′x′ : Bjm′(n1′,n2′) = M1xM2x : M1M2 ∑Dlmm′ Bjm′(n1,n2) = xx : ∑Dlmm′ Bjm′(n1,n2)
                            RxxY12_rot = x̂x̂_n1n2_rot' * DY
                            xxY1′2′_rot = x̂x̂_n1′n2′_rot' * Y1′2′_rot
                            @test isapproxdefault(xxY1′2′_rot, RxxY12_rot)
                        end
                    end
                end
            end
        end

        # Compare Equator with the prime meridian (as both are great circles)
        ϕ1_Mer = ϕ2_Mer = 0
        θ1_Equator = θ2_Equator = pi/2
        for θ1_Mer in LinRange(0, pi, 4), θ2_Mer in LinRange(0, pi, 4),
            Δθ = θ2_Mer - θ1_Mer
            BipolarSphericalHarmonics.monopolarharmonics!(B12, θ1_Mer, ϕ1_Mer, θ2_Mer, ϕ2_Mer);

            Un1 = basisconversionmatrix(Cartesian(), HelicityCovariant(), θ1_Mer, ϕ1_Mer);
            Un2 = basisconversionmatrix(Cartesian(), HelicityCovariant(), θ2_Mer, ϕ2_Mer);

            for ϕ1_Equator in LinRange(0, pi, 4)
                Δϕ = Δθ
                ϕ2_Equator = ϕ1_Equator + Δϕ
                U′n1 = basisconversionmatrix(Cartesian(), HelicityCovariant(), θ1_Equator, ϕ1_Equator);
                U′n2 = basisconversionmatrix(Cartesian(), HelicityCovariant(), θ2_Equator, ϕ2_Equator);
                BipolarSphericalHarmonics.monopolarharmonics!(B1′2′, θ1_Equator, ϕ1_Equator, θ2_Equator, ϕ2_Equator);

                R = rotn1n2n1′n2′(θ1_Mer, ϕ1_Mer, θ2_Mer, ϕ2_Mer, θ1_Equator, ϕ1_Equator, θ2_Equator, ϕ2_Equator)
                RSS′ = inv(R)
                α, β, γ = RSS′.theta1, RSS′.theta2, RSS′.theta3
                M′1 = U′n1 * RSS′' * Un1'
                M′2 = U′n2 * RSS′' * Un2'
                M′ = kron(M′1, M′2)
                @test M′ ≈ Diagonal(M′)

                for j in 0:2lmax
                    Dp = wignerD!(Dvec[j], j, α, β, γ);
                    D = OffsetArray(Dp, -j:j, -j:j);
                    for l1 in 0:lmax, l2 in abs(j - l1):min(lmax, l1 + j), m in -j:j
                        Y_rot = M′ * sum(D[m′, m] * biposh(B12, θ1_Mer, ϕ1_Mer, θ2_Mer, ϕ2_Mer, j, m′, l1, l2) for m′ in -j:j)
                        Y′ = biposh(B1′2′, θ1_Equator, ϕ1_Equator, θ2_Equator, ϕ2_Equator, j, m, l1, l2)
                        @test isapproxdefault(Y′, Y_rot)
                    end
                end
            end
        end
    end

    @testset "VSH" begin
        for YT in [Irreducible(), PB(), Hansen()], B in [Cartesian(), Polar(), SphericalCovariant(), HelicityCovariant()]
            tstart = time()
            B12 = cache(VSH(YT, B), lmax);
            B12_rot = cache(VSH(YT, B), lmax);
            B1′2′ = cache(VSH(YT, B), lmax);
            B1′2′_rot = cache(VSH(YT, B), lmax);
            Un1 = basisconversionmatrix(Cartesian(), B, θ1, ϕ1);
            Un2 = basisconversionmatrix(Cartesian(), B, θ2, ϕ2);
            BipolarSphericalHarmonics.monopolarharmonics!(B12, θ1, ϕ1, θ2, ϕ2);
            x̂_n1 = Un1 * SVector{3}(1,0,0)
            x̂_n2 = Un2 * SVector{3}(1,0,0)
            x̂x̂_n1n2 = kron(x̂_n1, x̂_n2)
            for θ1′ in LinRange(0, pi, 4), ϕ1′ in LinRange(0, 2pi, 5)
                # start by choosing the rotation that maps n1 to n1′
                n1′ = cartvec(θ1′, ϕ1′)
                R11′ = rotn1n1′(θ1, ϕ1, θ1′, ϕ1′)
                n2′ = R11′ * n2
                θ2′, ϕ2′ = polcoords(n2′)
                @test dot(n1, n2) ≈ dot(n1′, n2′)
                @test R11′ * n1 ≈ n1′
                @test R11′ * n2 ≈ n2′

                BipolarSphericalHarmonics.monopolarharmonics!(B1′2′, θ1′, ϕ1′, θ2′, ϕ2′);
                # The passive rotation for the frame is the inverse of R11′
                RSS′ = inv(R11′)
                α, β, γ = RSS′.theta1, RSS′.theta2, RSS′.theta3
                U′n1 = basisconversionmatrix(Cartesian(), B, θ1′, ϕ1′);
                U′n2 = basisconversionmatrix(Cartesian(), B, θ2′, ϕ2′);
                M′1 = U′n1 * R11′ * Un1'
                M′2 = U′n2 * R11′ * Un2'
                M′ = kron(M′1, M′2)

                x̂_n1′ = M′1 * x̂_n1
                x̂_n2′ = M′2 * x̂_n2
                x̂x̂_n1′n2′ = kron(x̂_n1′, x̂_n2′)
                R⁻¹x̂R⁻¹x̂_n1n2 = kron(inv(M′1) * x̂_n1, inv(M′2) * x̂_n2)

                for j in 0:lmax
                    Dp = wignerD!(Dvec[j], j, α, β, γ);
                    D = OffsetArray(Dp, -j:j, -j:j);
                    for l1 in 0:lmax, l2 in abs(j - l1):min(lmax, l1 + j)
                        for m in -j:j
                            DY = sum(D[m′, m] * biposh(B12, θ1, ϕ1, θ2, ϕ2, j, m′, l1, l2) for m′ in -j:j)
                            RY = M′ * DY
                            Y′ = biposh(B1′2′, θ1′, ϕ1′, θ2′, ϕ2′, j, m, l1, l2)
                            @test isapproxdefault(Y′, RY)

                            # x′x′ : Bjm′(n1′,n2′) = M1xM2x : M1M2 ∑Dlmm′ Bjm′(n1,n2) = xx : ∑Dlmm′ Bjm′(n1,n2)
                            RxxY12 = x̂x̂_n1n2' * DY
                            xxY1′2′ = x̂x̂_n1′n2′' * Y′
                            @test isapproxdefault(xxY1′2′, RxxY12)

                            # xx : Bjm′(n1′,n2′) = xx : M1M2 ∑Dlmm′ Bjm′(n1,n2) = M1⁻¹xM2⁻¹x : ∑Dlmm′ Bjm′(n1,n2)
                            xxY12_rot = R⁻¹x̂R⁻¹x̂_n1n2' * DY
                            xxY1′2′ = x̂x̂_n1n2' * Y′
                            @test isapproxdefault(xxY1′2′, xxY12_rot)
                        end
                    end
                end

                # next rotate (θ2,ϕ2) and (θ2′,ϕ2′) to new positions
                for ω2 in LinRange(0, pi, 4)[2:end]
                    R = AngleAxis(ω2, n1...)
                    n2_rot = R * n2
                    θ2_rot, ϕ2_rot = polcoords(n2_rot)
                    @test dot(n1, n2_rot) ≈ dot(n1, n2)
                    BipolarSphericalHarmonics.monopolarharmonics!(B12_rot, θ1, ϕ1, θ2_rot, ϕ2_rot);
                    Un2_rot = basisconversionmatrix(Cartesian(), B, θ2_rot, ϕ2_rot);

                    RSS′ = inv(R) # passive rotation that maps n2 to n2_rot
                    M22_rot = Un2_rot * R * Un2'

                    x̂_n2_rot = M22_rot * SVector{3}(1,0,0)
                    x̂x̂_n1n2_rot = kron(x̂_n1, x̂_n2_rot)

                    for ω2′ in LinRange(0, pi, 4)[2:end]
                        n2′_rot = AngleAxis(ω2′, n1′...) * n2′
                        θ2′_rot, ϕ2′_rot = polcoords(n2′_rot)
                        @test dot(n1′, n2′_rot) ≈ dot(n1′, n2′) ≈ dot(n1, n2)
                        R = rotn1n2n1′n2′(θ1, ϕ1, θ2_rot, ϕ2_rot, θ1′, ϕ1′, θ2′_rot, ϕ2′_rot)
                        @test det(R) ≈ 1
                        @test R'R ≈ I
                        @test R * n1 ≈ n1′
                        @test R * n2_rot ≈ n2′_rot

                        RSS′ = inv(R)
                        α, β, γ = RSS′.theta1, RSS′.theta2, RSS′.theta3
                        U′n2_rot = basisconversionmatrix(Cartesian(), B, θ2′_rot, ϕ2′_rot);
                        M′1 = U′n1 * R * Un1'
                        M′2 = U′n2_rot * R * Un2_rot'
                        M′ = kron(M′1, M′2)

                        x̂x̂_n1′n2′_rot = kron(M′1 * x̂_n1, M′2 * x̂_n2_rot)
                        R⁻¹x̂R⁻¹x̂_n1n2_rot = kron(inv(M′1) * x̂_n1, inv(M′2) * x̂_n2_rot)

                        BipolarSphericalHarmonics.monopolarharmonics!(B1′2′_rot, θ1′, ϕ1′, θ2′_rot, ϕ2′_rot);

                        for j in 0:lmax
                            Dp = wignerD!(Dvec[j], j, α, β, γ);
                            D = OffsetArray(Dp, -j:j, -j:j)
                            for l1 in 0:2:lmax, l2 in abs(j - l1):2:min(lmax, l1 + j), m in -j:2:j
                                DY = sum(D[m′, m] * biposh(B12_rot, θ1, ϕ1, θ2_rot, ϕ2_rot, j, m′, l1, l2) for m′ in -j:j)
                                RY12_rot = M′ * DY
                                Y1′2′_rot = biposh(B1′2′_rot, θ1′, ϕ1′, θ2′_rot, ϕ2′_rot, j, m, l1, l2)
                                @test isapproxdefault(Y1′2′_rot, RY12_rot)

                                # x′x′ : Bjm′(n1′,n2′) = M1xM2x : M1M2 ∑Dlmm′ Bjm′(n1,n2) = xx : ∑Dlmm′ Bjm′(n1,n2)
                                xxY1′2′_rot = x̂x̂_n1′n2′_rot' * Y1′2′_rot
                                RxxY12_rot = x̂x̂_n1n2_rot' * DY
                                @test isapproxdefault(xxY1′2′_rot, RxxY12_rot)

                                # xx : Bjm′(n1′,n2′) = xx : M1M2 ∑Dlmm′ Bjm′(n1,n2) = M1⁻¹xM2⁻¹x : ∑Dlmm′ Bjm′(n1,n2)
                                xxY1′2′_rot = x̂x̂_n1n2_rot' * Y1′2′_rot
                                R⁻¹xR⁻¹xY12_rot = R⁻¹x̂R⁻¹x̂_n1n2_rot' * DY
                                @test isapproxdefault(xxY1′2′_rot, R⁻¹xR⁻¹xY12_rot)
                            end
                        end
                    end
                end
            end

            # Compare Equator with the prime meridian (as both are great circles)
            ϕ1_Mer = ϕ2_Mer = 0
            θ1_Equator = θ2_Equator = pi/2
            for θ1_Mer in LinRange(0, pi, 4), θ2_Mer in LinRange(0, pi, 4),
                Δθ = θ2_Mer - θ1_Mer
                BipolarSphericalHarmonics.monopolarharmonics!(B12, θ1_Mer, ϕ1_Mer, θ2_Mer, ϕ2_Mer);

                Un1 = basisconversionmatrix(Cartesian(), B, θ1_Mer, ϕ1_Mer);
                Un2 = basisconversionmatrix(Cartesian(), B, θ2_Mer, ϕ2_Mer);

                for ϕ1_Equator in LinRange(0, pi, 4)
                    Δϕ = Δθ
                    ϕ2_Equator = ϕ1_Equator + Δϕ
                    U′n1 = basisconversionmatrix(Cartesian(), B, θ1_Equator, ϕ1_Equator);
                    U′n2 = basisconversionmatrix(Cartesian(), B, θ2_Equator, ϕ2_Equator);
                    BipolarSphericalHarmonics.monopolarharmonics!(B1′2′, θ1_Equator, ϕ1_Equator, θ2_Equator, ϕ2_Equator);

                    R = rotn1n2n1′n2′(θ1_Mer, ϕ1_Mer, θ2_Mer, ϕ2_Mer, θ1_Equator, ϕ1_Equator, θ2_Equator, ϕ2_Equator)
                    RSS′ = inv(R)
                    α, β, γ = RSS′.theta1, RSS′.theta2, RSS′.theta3
                    M′1 = U′n1 * R * Un1'
                    M′2 = U′n2 * R * Un2'
                    M′ = kron(M′1, M′2)

                    for j in 0:lmax
                        Dp = wignerD!(Dvec[j], j, α, β, γ);
                        D = OffsetArray(Dp, -j:j, -j:j);
                        for l1 in 0:2:lmax, l2 in abs(j - l1):2:min(lmax, l1 + j), m in -j:2:j
                            RY_mer = M′ * sum(D[m′, m] * biposh(B12, θ1_Mer, ϕ1_Mer, θ2_Mer, ϕ2_Mer, j, m′, l1, l2) for m′ in -j:j)
                            Y_Equator = biposh(B1′2′, θ1_Equator, ϕ1_Equator, θ2_Equator, ϕ2_Equator, j, m, l1, l2)
                            @test isapproxdefault(Y_Equator, RY_mer)
                        end
                    end
                end
            end

            tend = time()
            @info "Finished testing $YT $B in $(round(tend - tstart, sigdigits = 1)) seconds"
        end
    end
end

@testset "allocation" begin
    θ1, ϕ1, θ2, ϕ2 = pi*rand(), 2pi*rand(), pi*rand(), 2pi*rand()

    function test_alloc(SHT, θ1, ϕ1, θ2, ϕ2)
        B12 = monopolarharmonics(SHT, θ1, ϕ1, θ2, ϕ2, 1, 1)
        Y12 = zeros(BipolarSphericalHarmonics.eltypeY(B12), 1)
        biposh!(Y12, B12, θ1, ϕ1, θ2, ϕ2, 1, 1, 1, 1)
        @test (@allocated biposh!(Y12, B12, θ1, ϕ1, θ2, ϕ2, 1, 1, 1, 1)) == 0
    end

    test_alloc(SH(), θ1, ϕ1, θ2, ϕ2)
    test_alloc(GSH(), θ1, ϕ1, θ2, ϕ2)
    for YT in [PB(), Hansen(), Irreducible()], B in [Polar(), Cartesian(), HelicityCovariant(), SphericalCovariant()]
        test_alloc(VSH(YT, B), θ1, ϕ1, θ2, ϕ2)
    end
end


text:
#=
P05 (*) Reverse a list.
=#

function my_reverse(coll)
    # With built-in `reverse` function
    reverse(coll)
end

function reverse_with_reversed(coll, reversed)
    # With tail recursion
    if isempty(coll)
        reversed
    else
        insert!(reversed, 1, coll[1])
        reverse_with_reversed(coll[2:end], reversed)
    end
end

function my_reverse_2(coll)
    # Reversed collection initially empty
    reverse_with_reversed(coll, [])
end

@assert my_reverse([]) == []
@assert my_reverse([:a, :b, :c, :d]) == [:d, :c, :b, :a]
@assert my_reverse(Vector{Symbol}[[:a, :b], [:c, :d]]) == Vector{Symbol}[[:c, :d], [:a, :b]]

@assert my_reverse_2([]) == []
@assert my_reverse_2([:a, :b, :c, :d]) == [:d, :c, :b, :a]
my_reverse_2(Vector{Symbol}[[:a, :b], [:c, :d]])
@assert my_reverse_2(Vector{Symbol}[[:a, :b], [:c, :d]]) == Vector{Symbol}[[:c, :d], [:a, :b]]

println("Tests passed: JL-05.jl")


text:
# Lazy broadcasted dictionary

struct BroadcastedDictionary{I, T, F, Data <: Tuple} <: AbstractDictionary{I, T}
    f::F
    data::Data
    sharetokens::Bool
end

_f(d::BroadcastedDictionary) = getfield(d, :f)
_data(d::BroadcastedDictionary) = getfield(d, :data)

@propagate_inbounds function BroadcastedDictionary(f, data)
    dicts = _dicts(data...)
    sharetokens = _sharetokens(dicts...)
    I = keytype(dicts[1])
    Ts = Base.Broadcast.eltypes(data)
    T = Core.Compiler.return_type(f, Ts)

    return BroadcastedDictionary{I, T, typeof(f), typeof(data)}(f, data, sharetokens)
end

@inline Base.keys(d::BroadcastedDictionary) = _keys(_data(d)...)

@propagate_inbounds function Base.getindex(d::BroadcastedDictionary{I}, i::I) where {I}
    if istokenizable(d)
        t = gettoken(d, i)
        return _f(d)(_gettokenvalue(t, _data(d)...)...)
    else
        return _f(d)(_getindex(i, _data(d)...)...)
    end
end

function Base.isassigned(d::BroadcastedDictionary{I}, i::I) where {I}
    return _isassigned(i, _data(d)...)
end

istokenizable(d::BroadcastedDictionary) = getfield(d, :sharetokens)
function tokens(d::BroadcastedDictionary)
    _tokens(_data(d)...)
end

@propagate_inbounds function gettoken(d::BroadcastedDictionary{I}, i::I) where {I}
    return gettoken(_tokens(_data(d)...), i)
end

function istokenassigned(d::BroadcastedDictionary, t)
    return _istokenassigned(t, _data(d)...)
end

@propagate_inbounds function gettokenvalue(d::BroadcastedDictionary, t)
    return _f(d)(_gettokenvalue(t, _data(d)...)...)
end

@inline function Base.similar(d::BroadcastedDictionary, ::Type{T}) where {T}
    return similar(_dicts(d.data...)[1], T)
end

@inline _dicts(d::AbstractDictionary, ds...) = (d, _dicts(ds...)...)
@inline _dicts(d, ds...) = (_dicts(ds...)...,)
_dicts() = ()

@inline _tokens(d::AbstractDictionary, ds...) = tokens(d)
@inline _tokens(d, ds...) = _tokens(ds...)

@inline _keys(d::AbstractDictionary, ds...) = keys(d)
@inline _keys(d, ds...) = _keys(ds...)

_sharetokens(d) = true
@propagate_inbounds function _sharetokens(d, d2, ds...)
    if sharetokens(d, d2)
        return _sharetokens(d, ds...)
    else
        @boundscheck if !isequal(keys(d), keys(d2))
            throw(IndexError("Indices do not match"))
        end
        return false
    end
end

_getindex(i) = ()
@propagate_inbounds _getindex(i, d::AbstractDictionary, ds...) = (d[i], _getindex(i, ds...)...)
@propagate_inbounds _getindex(i, d::Ref, ds...) = (d[], _getindex(i, ds...)...)
@propagate_inbounds _getindex(i, d, ds...) = (d[CartesianIndex()], _getindex(i, ds...)...)

_gettokenvalue(t) = ()
@propagate_inbounds _gettokenvalue(t, d::AbstractDictionary, ds...) = (gettokenvalue(d, t), _gettokenvalue(t, ds...)...)
@propagate_inbounds _gettokenvalue(t, d::Ref, ds...) = (d[], _gettokenvalue(t, ds...)...)
@propagate_inbounds _gettokenvalue(t, d, ds...) = (d[CartesianIndex()], _gettokenvalue(t, ds...)...)

_isassigned(i) = true
_isassigned(i, d, ds...) = _isassigned(i, ds...)
function _isassigned(i, d::AbstractDictionary, ds...)
    if isassigned(d, i)
        return _isassigned(i, ds...)
    else
        return false
    end
end

_istokenassigned(t) = true
_istokenassigned(t, d, ds...) = _istokenassigned(t, ds...)
function _istokenassigned(t, d::AbstractDictionary, ds...)
    if istokenassigned(d, t)
        return _istokenassigned(t, ds...)
    else
        return false
    end
end

## Hook into `Base.Broadcast` machinery

Base.Broadcast.broadcastable(d::AbstractDictionary) = d

struct DictionaryStyle <: Base.Broadcast.BroadcastStyle
end

Base.Broadcast.BroadcastStyle(::Type{<:AbstractDictionary}) = DictionaryStyle()
Base.Broadcast.BroadcastStyle(::DictionaryStyle, ::Base.Broadcast.AbstractArrayStyle{0}) = DictionaryStyle()

# Now we overload `broadcasted` directly
function Base.Broadcast.broadcasted(::DictionaryStyle, f, args...)
    return BroadcastedDictionary(f, args)
end

Base.Broadcast.materialize(d::BroadcastedDictionary) = copy(d)
Base.Broadcast.materialize!(out::AbstractDictionary, d::BroadcastedDictionary) = copyto!(out, d)


text:
export KSP, KSPSetOperators, KSPSetFromOptions, KSPSolve, KSPSolveTranspose, KSPSetUp, KSP_NULL

export KSPGetConvergedReason, PetscView, KSPSetType, KSPGetType, KSPSetTolerances, KSPGetTolerances, KSPSetInitialGuessNonzero, KSPGetInitialGuessNonzero, KSPGetResidualNorm

type KSP
  pobj::Ptr{Void}

  function KSP(comm::MPI_Comm)  # constructor for KSP owned by the user
      ptr = Array(Ptr{Void}, 1)
      ierr = ccall((:KSPCreate,petsc),PetscErrorCode,(comm_type, Ptr{Void}),comm, ptr)
      @assert(ierr == 0)
      obj = new(ptr[1])
#      finalizer(obj, PetscDestroy)
      return obj
  end

  function KSP(ptr::Ptr{Void})  # constructor for KSP *not* owned by the user
    return new(ptr)
  end
end

global const KSP_NULL = KSP(C_NULL)

function PetscDestroy(ksp::KSP)
    err = ccall((:KSPDestroy,petsc),PetscErrorCode,(Ptr{Ptr{Void}},),&ksp.pobj)
    println("KSPDestroy called")
#    sleep(5)
end

function PetscView(ksp::KSP, viewer)
    ccall((:KSPView,petsc),PetscErrorCode,(Ptr{Void}, Int64), ksp.pobj, viewer)
end

function KSPSetOperators(ksp::KSP,Amat::PetscMat,Pmat::PetscMat)
   err = ccall((:KSPSetOperators,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void},Ptr{Void}), ksp.pobj, Amat.pobj, Pmat.pobj)
end

function KSPSetFromOptions(ksp::KSP)
    ccall((:KSPSetFromOptions,petsc),PetscErrorCode,(Ptr{Void},),ksp.pobj)
end

function KSPSolve(ksp::KSP, b::PetscVec, x::PetscVec)
    err = ccall((:KSPSolve,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void},Ptr{Void}), ksp.pobj, b.pobj, x.pobj)
end

function KSPSolveTranspose(ksp::KSP, b::PetscVec, x::PetscVec)
    err = ccall((:KSPSolveTranspose,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void},Ptr{Void}), ksp.pobj, b.pobj, x.pobj)
end

function KSPSetUp(ksp::KSP)
    err = ccall((:KSPSetUp,petsc),PetscErrorCode,(Ptr{Void},), ksp.pobj)
end

function KSPGetConvergedReason(ksp::KSP)
    reason = Array(KSPConvergedReason, 1)
    ccall((:KSPGetConvergedReason,petsc),PetscErrorCode,(Ptr{Void},Ptr{KSPConvergedReason}), ksp.pobj, reason)

    return reason[1]
end

function KSPSetType(ksp::KSP, ksptype::KSPType)
    ccall((:KSPSetType,petsc),PetscErrorCode,(Ptr{Void},Cstring), ksp.pobj, ksptype)
end

function KSPGetType(ksp::KSP)
    ksptype = Array(Ptr{UInt8}, 1)
    ccall((:KSPGetType,petsc),PetscErrorCode,(Ptr{Void}, Ptr{Ptr{UInt8}}), ksp.pobj, ksptype)

    return bytestring(ksptype[1])
end

function KSPSetTolerances(ksp::KSP,rtol::PetscReal, abstol::PetscReal, dtol::PetscReal, maxits::PetscInt)

    ccall((:KSPSetTolerances,petsc),PetscErrorCode,(Ptr{Void}, PetscReal, PetscReal, PetscReal, PetscInt), ksp.pobj, rtol, abstol, dtol, maxits)
end

function KSPGetTolerances(ksp::KSP)
    rtol = Array(PetscReal,1)
    abstol = Array(PetscReal,1)
    dtol = Array(PetscReal,1)
    maxits = Array(PetscInt, 1)

    ccall((:KSPGetTolerances,petsc),PetscErrorCode,(KSP,Ptr{PetscReal},Ptr{PetscReal},Ptr{PetscReal},Ptr{PetscInt}), ksp, rtol, abstol, dtol, maxits)

    return rtol[1], abstol[1], dtol[1], maxits[1]
end

function KSPSetInitialGuessNonzero(ksp::KSP, flg::PetscBool)
    ccall((:KSPSetInitialGuessNonzero,petsc),PetscErrorCode,(Ptr{Void},PetscBool), ksp.pobj, flg)
end

function KSPGetInitialGuessNonzero(ksp::KSP)
    flg_arr = Array(PetscBool, 1)
    ccall((:KSPGetInitialGuessNonzero,petsc),PetscErrorCode,(Ptr{Void}, Ptr{PetscBool}), ksp.pobj, flg_arr)
    return flg_arr[1]
end

function KSPGetResidualNorm(ksp::KSP)
   rnorm = Array(PetscReal, 1)
   ccall((:KSPGetResidualNorm,petsc),PetscErrorCode,(Ptr{Void}, Ptr{PetscReal}),ksp.pobj, rnorm)

    return rnorm[1]
end

### new function
# not tested
function KSPSetReusePreconditioner(arg1::KSP,arg2::PetscBool)
    ccall((:KSPSetReusePreconditioner,petsc),PetscErrorCode,(KSP,PetscBool),arg1,arg2)
end


text:
# Code derived from metabolicEP (https://github.com/anna-pa-m/Metabolic-EP)

struct EPMatT0{T<:AbstractFloat} <: AbstractEPMat
    Σd::Matrix{T}
    Σi::Matrix{T}
    G::Matrix{T}
    lb::Vector{T}
    ub::Vector{T}
    vd::Vector{T}
    vi::Vector{T}
    Y::Vector{T}
    usperm::Vector{Int} # unsorted permutation (see echelonize)
end

function EPMatT0(K::AbstractArray{T,2}, Y::Vector{T},
        lb::Vector{T}, ub::Vector{T}) where T <: Real
    M, N = size(K)
    M > N && @warn("numeber of rows M=$M larger than number of cols N=$N")
    _, usperm, EK, EY = echelonize(Matrix(K),Y)
    Mech, Nech = size(EK)
    length(EY) == Mech || error("vector size incompatible with matrix")
    # sortperm ci is the inverse permutation that sends me back to the original model permutation
    return EPMatT0(
        #=Σd=#  zeros(T,Mech,Mech),
        #=Σi=#  zeros(T,Nech-Mech,Nech-Mech),
        #=G=#   copy(EK[1:Mech,Mech+1:Nech]),
        #=lb=#  lb[usperm],
        #=ub=#  ub[usperm],
        #=vd=#  zeros(T,Mech),
        #=vi=#  zeros(T,Nech-Mech),
        #=Y=#   EY,
        #=idx=# usperm
    )
end


text:
#
# Powell.jl --
#
# Mike Powell's derivative free optimization methods for Julia.
#
# -----------------------------------------------------------------------------
#
# This file is part of OptimPack.jl which is licensed under the MIT
# "Expat" License:
#
# Copyright (C) 2015-2019, Éric Thiébaut <https://github.com/emmt/OptimPack.jl>.
#

module Powell

export
    iterate,
    restart,
    getstatus,
    getreason,
    getradius,
    getncalls,
    getlastf,
    Cobyla,
    cobyla,
    cobyla!,
    Newuoa,
    newuoa,
    newuoa!,
    Bobyqa,
    bobyqa,
    bobyqa!

using Compat

import Base: ==, iterate

# The dynamic libraries implementing the methods.
import ..libcobyla
import ..libbobyqa
import ..libnewuoa

abstract type AbstractContext end
abstract type AbstractStatus end

==(a::T, b::T) where {T<:AbstractStatus} = (a._code == b._code)
==(a::AbstractStatus, b::AbstractStatus) = false

"""
The `iterate(ctx, ...)` method performs the next iteration of the reverse
communication associated with the context `ctx`.  Other arguments depend on the
type of algorithm.

For **COBYLA** algorithm, the next iteration is performed by:

    iterate(ctx, f, x, c) -> status

or

    iterate(ctx, f, x) -> status

on entry, the workspace status must be `COBYLA_ITERATE`, `f` and `c` are the
function value and the constraints at `x`, the latter can be omitted if there
are no constraints.  On exit, the returned value (the new workspace status) is:
`COBYLA_ITERATE` if a new trial point has been stored in `x` and if user is
requested to compute the function value and the constraints on the new point;
`COBYLA_SUCCESS` if algorithm has converged and `x` has been set with the
variables at the solution (the corresponding function value can be retrieved
with `getlastf`); anything else indicates an error (see `getreason`
for an explanatory message).

For **NEWUOA** algorithm, the next iteration is performed by:

    iterate(ctx, f, x) -> status

on entry, the wokspace status must be `NEWUOA_ITERATE`, `f` is the function
value at `x`.  On exit, the returned value (the new wokspace status) is:
`NEWUOA_ITERATE` if a new trial point has been stored in `x` and if user is
requested to compute the function value for the new point; `NEWUOA_SUCCESS` if
algorithm has converged; anything else indicates an error (see `getreason` for
an explanatory message).

"""
function iterate end

"""

    restart(ctx) -> status

restarts the reverse communication algorithm associated with the context `ctx`
using the same parameters.  The return value is the new status of the
algorithm, see `getstatus` for details.

"""
function restart end

"""

    getstatus(ctx) -> status

get the current status of the reverse communication algorithm associated with
the context `ctx`.  Possible values are:

* for **COBYLA**: `COBYLA_ITERATE`, if user is requested to compute `f(x)` and
  `c(x)`; `COBYLA_SUCCESS`, if algorithm has converged;

* for **NEWUOA**: `NEWUOA_ITERATE`, if user is requested to compute `f(x)`;
  `NEWUOA_SUCCESS`, if algorithm has converged;

Anything else indicates an error (see `getreason` for an explanatory message).

"""
function getstatus end

"""

    getreason(ctx) -> msg

or

    getreason(status) -> msg

get an explanatory message about the current status of the reverse
communication algorithm associated with the context `ctx` or with the status
returned by an optimization method of by `getstatus(ctx)`.

"""
function getreason end

getreason(ctx::AbstractContext) = getreason(getstatus(ctx))

"""

    getlastf(ctx) -> fx

get the last function value in the reverse communication algorithm associated
with the context `ctx`.  Upon convergence of `iterate`, this value corresponds
to the function at the solution; otherwise, this value corresponds to the
previous set of variables.

"""
function getlastf end

"""

    getncalls(ctx) -> nevals

get the current number of function evaluations in the reverse communication
algorithm associated with the context `ctx`.  Result is -1 if something is
wrong, nonnegative otherwise.

"""
function getncalls end

"""

    getradius(ctx) -> rho

get the current size of the trust region of the reverse communication algorithm
associated with the context `ctx`.  Result is 0 if algorithm not yet started
(before first iteration), -1 if something is wrong, strictly positive
otherwise.

"""
function getradius end

"""

```julia
grow!(x, n) -> x
```

grows vector `x` so that it has at least `n` elements, does nothing if `x` is
large enough.  Argument `x` is returned.

See also [`resize!`](@ref).

"""
grow!(x::Vector, n::Integer) = begin
    length(x) < n && resize!(x, n)
    return x
end

include("newuoa.jl")
import .Newuoa: newuoa, newuoa!

include("cobyla.jl")
import .Cobyla: cobyla, cobyla!

include("bobyqa.jl")
import .Bobyqa: bobyqa, bobyqa!

end # module Powell


text:
using GeometricIntegrators
using ChargedParticleDynamics.GuidingCenter4d.TokamakFastLoop
using GeometricExamples

const Δt     = 2.
const ntime  = 25000
const nloop  = 200
const nsave  = 25
const nplot  = 10
const run_id = "poincare_1st_fast_firk_dt2"

include("../guiding_center_4d_settings_firk.jl")

tableau_list = get_tableau_list_firk()

pinv = guiding_center_4d_ode_poincare_invariant_1st(Δt, nloop, ntime, nsave)

include("guiding_center_4d_fast_poincare_invariant_1st.jl")