text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing

    Base.precompile(Tuple{typeof(LoopVectorization._avx_loopset),Core.SimpleVector,Core.SimpleVector,Core.SimpleVector,Core.SimpleVector,Core.SimpleVector,Any})
    Base.precompile(Tuple{typeof(LoopVectorization.add_broadcast!),LoopVectorization.LoopSet,Symbol,Symbol,Vector{Symbol},Type{Adjoint{Float64, Matrix{Float64}}},Int64})
    Base.precompile(Tuple{typeof(LoopVectorization.add_broadcast!),LoopVectorization.LoopSet,Symbol,Symbol,Vector{Symbol},Type{LowDimArray{(false,), Float32, 3, Array{Float32, 3}}},Int64})
    Base.precompile(Tuple{typeof(LoopVectorization.add_broadcast!),LoopVectorization.LoopSet,Symbol,Symbol,Vector{Symbol},Type{SubArray{Float32, 2, Array{Float32, 3}, Tuple{Int64, Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}}, true}},Int64})
    Base.precompile(Tuple{typeof(LoopVectorization.add_broadcast!),LoopVectorization.LoopSet,Symbol,Symbol,Vector{Symbol},Type{SubArray{Float64, 2, Array{Float64, 3}, Tuple{Int64, Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}}, true}},Int64})
    Base.precompile(Tuple{typeof(LoopVectorization.add_broadcast!),LoopVectorization.LoopSet,Symbol,Symbol,Vector{Symbol},Type{SubArray{Int32, 2, Array{Int32, 3}, Tuple{Int64, Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}}, true}},Int64})
    Base.precompile(Tuple{typeof(LoopVectorization.add_broadcast!),LoopVectorization.LoopSet,Symbol,Symbol,Vector{Symbol},Type{SubArray{Int64, 2, Array{Int64, 3}, Tuple{Int64, Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}}, true}},Int64})
    Base.precompile(Tuple{typeof(LoopVectorization.add_compute!),LoopVectorization.LoopSet,Symbol,Expr,Int64,Int64,Nothing})
    Base.precompile(Tuple{typeof(LoopVectorization.add_constant!),LoopVectorization.LoopSet,Float64,Vector{Symbol},Symbol,Int64})
    Base.precompile(Tuple{typeof(LoopVectorization.add_if!),LoopVectorization.LoopSet,Symbol,Expr,Int64,Int64,LoopVectorization.ArrayReferenceMetaPosition})
    Base.precompile(Tuple{typeof(avx_body), LoopSet, Tuple{Int8,Int8,Int8,Int}})
    Base.precompile(Tuple{typeof(avx_loopset), Vector{Instruction}, Vector{OperationStruct}, Vector{ArrayRefStruct}, Core.SimpleVector, Core.SimpleVector, Core.SimpleVector, Any})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedBitPointer{2, 1, 0, (1, 2), Tuple{ArrayInterface.StaticInt{1}, Int64}, Tuple{ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float32, 2, 1, 0, (1, 2), Tuple{ArrayInterface.StaticInt{4}, ArrayInterface.StaticInt{292}}, Tuple{ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float32, 2, 1, 0, (1, 2), Tuple{ArrayInterface.StaticInt{4}, ArrayInterface.StaticInt{308}}, Tuple{ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float32, 3, 1, 0, (1, 2, 3), Tuple{ArrayInterface.StaticInt{4}, Int64, Int64}, Tuple{ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float32, 5, -1, 0, (2, 3, 4, 5, 6), NTuple{5, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float32, 5, 1, 0, (1, 2, 3, 4, 5), Tuple{ArrayInterface.StaticInt{4}, Int64, Int64, Int64, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float32, 5, 1, 0, (1, 2, 3, 4, 6), Tuple{ArrayInterface.StaticInt{4}, Int64, Int64, Int64, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float32, 5, 1, 0, (1, 2, 3, 5, 6), Tuple{ArrayInterface.StaticInt{4}, Int64, Int64, Int64, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float32, 5, 1, 0, (1, 2, 4, 5, 6), Tuple{ArrayInterface.StaticInt{4}, Int64, Int64, Int64, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float32, 5, 1, 0, (1, 3, 4, 5, 6), Tuple{ArrayInterface.StaticInt{4}, Int64, Int64, Int64, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float32, 6, 1, 0, (1, 2, 3, 4, 5, 6), Tuple{ArrayInterface.StaticInt{4}, Int64, Int64, Int64, Int64, Int64}, NTuple{6, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 1, -1, -1, (3,), Tuple{Int64}, Tuple{ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{ArrayInterface.StaticInt{8}}, Tuple{ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 2, -1, -1, (2, 3), Tuple{Int64, Int64}, Tuple{ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 2, 1, 0, (1, 2), Tuple{ArrayInterface.StaticInt{8}, ArrayInterface.StaticInt{584}}, Tuple{ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 2, 2, 0, (2, 1), Tuple{ArrayInterface.StaticInt{584}, ArrayInterface.StaticInt{8}}, Tuple{ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 2, 2, 0, (2, 1), Tuple{Int64, ArrayInterface.StaticInt{8}}, Tuple{ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 2, 2, 0, (3, 1), Tuple{Int64, ArrayInterface.StaticInt{8}}, Tuple{ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 3, 1, 0, (1, 2, 3), Tuple{ArrayInterface.StaticInt{8}, Int64, Int64}, Tuple{ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 3, 2, 0, (2, 1, 4), Tuple{Int64, ArrayInterface.StaticInt{8}, Int64}, Tuple{ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 4, -1, 0, (2, 3, 4, 5), NTuple{4, Int64}, NTuple{4, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 4, 1, 0, (1, 2, 3, 4), Tuple{ArrayInterface.StaticInt{8}, Int64, Int64, Int64}, NTuple{4, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 4, 1, 0, (1, 2, 3, 5), Tuple{ArrayInterface.StaticInt{8}, Int64, Int64, Int64}, NTuple{4, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 4, 1, 0, (1, 2, 4, 5), Tuple{ArrayInterface.StaticInt{8}, Int64, Int64, Int64}, NTuple{4, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 4, 1, 0, (1, 3, 4, 5), Tuple{ArrayInterface.StaticInt{8}, Int64, Int64, Int64}, NTuple{4, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 4, 2, 0, (2, 1, 4, 3), Tuple{Int64, ArrayInterface.StaticInt{8}, Int64, Int64}, NTuple{4, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 4, 2, 0, (3, 1, 4, 2), Tuple{Int64, ArrayInterface.StaticInt{8}, Int64, Int64}, NTuple{4, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 5, -1, 0, (2, 3, 4, 5, 6), NTuple{5, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 5, 1, 0, (1, 2, 3, 4, 5), Tuple{ArrayInterface.StaticInt{8}, Int64, Int64, Int64, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 5, 1, 0, (1, 2, 3, 4, 6), Tuple{ArrayInterface.StaticInt{8}, Int64, Int64, Int64, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 5, 1, 0, (1, 2, 3, 5, 6), Tuple{ArrayInterface.StaticInt{8}, Int64, Int64, Int64, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 5, 1, 0, (1, 2, 4, 5, 6), Tuple{ArrayInterface.StaticInt{8}, Int64, Int64, Int64, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 5, 1, 0, (1, 3, 4, 5, 6), Tuple{ArrayInterface.StaticInt{8}, Int64, Int64, Int64, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Float64, 6, 1, 0, (1, 2, 3, 4, 5, 6), Tuple{ArrayInterface.StaticInt{8}, Int64, Int64, Int64, Int64, Int64}, NTuple{6, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Int32, 2, 2, 0, (2, 1), Tuple{Int64, ArrayInterface.StaticInt{4}}, Tuple{ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Int32, 3, 1, 0, (1, 2, 3), Tuple{ArrayInterface.StaticInt{4}, Int64, Int64}, Tuple{ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Int32, 5, -1, 0, (2, 3, 4, 5, 6), NTuple{5, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Int32, 5, 1, 0, (1, 2, 3, 4, 5), Tuple{ArrayInterface.StaticInt{4}, Int64, Int64, Int64, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Int32, 5, 1, 0, (1, 2, 3, 4, 6), Tuple{ArrayInterface.StaticInt{4}, Int64, Int64, Int64, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Int32, 5, 1, 0, (1, 2, 3, 5, 6), Tuple{ArrayInterface.StaticInt{4}, Int64, Int64, Int64, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Int32, 5, 1, 0, (1, 2, 4, 5, 6), Tuple{ArrayInterface.StaticInt{4}, Int64, Int64, Int64, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Int32, 5, 1, 0, (1, 3, 4, 5, 6), Tuple{ArrayInterface.StaticInt{4}, Int64, Int64, Int64, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Int32, 6, 1, 0, (1, 2, 3, 4, 5, 6), Tuple{ArrayInterface.StaticInt{4}, Int64, Int64, Int64, Int64, Int64}, NTuple{6, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Int64, 2, 2, 0, (2, 1), Tuple{ArrayInterface.StaticInt{552}, ArrayInterface.StaticInt{8}}, Tuple{ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Int64, 3, 1, 0, (1, 2, 3), Tuple{ArrayInterface.StaticInt{8}, Int64, Int64}, Tuple{ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Int64, 5, -1, 0, (2, 3, 4, 5, 6), NTuple{5, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Int64, 5, 1, 0, (1, 2, 3, 4, 5), Tuple{ArrayInterface.StaticInt{8}, Int64, Int64, Int64, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Int64, 5, 1, 0, (1, 2, 3, 4, 6), Tuple{ArrayInterface.StaticInt{8}, Int64, Int64, Int64, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Int64, 5, 1, 0, (1, 2, 3, 5, 6), Tuple{ArrayInterface.StaticInt{8}, Int64, Int64, Int64, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Int64, 5, 1, 0, (1, 2, 4, 5, 6), Tuple{ArrayInterface.StaticInt{8}, Int64, Int64, Int64, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Int64, 5, 1, 0, (1, 3, 4, 5, 6), Tuple{ArrayInterface.StaticInt{8}, Int64, Int64, Int64, Int64}, NTuple{5, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.StridedPointer{Int64, 6, 1, 0, (1, 2, 3, 4, 5, 6), Tuple{ArrayInterface.StaticInt{8}, Int64, Int64, Int64, Int64, Int64}, NTuple{6, ArrayInterface.StaticInt{1}}}},Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.cost_vec_buf),LoopVectorization.LoopSet})
    Base.precompile(Tuple{typeof(LoopVectorization.evaluate_cost_tile),LoopVectorization.LoopSet,Vector{Symbol},LoopVectorization.UnrollSymbols})
    Base.precompile(Tuple{typeof(LoopVectorization.pushop!),LoopVectorization.LoopSet,LoopVectorization.Operation,Symbol})
    Base.precompile(Tuple{typeof(LoopVectorization.recursive_muladd_search!),Expr,Vector{Any},Bool,Bool})
    Base.precompile(Tuple{typeof(LoopVectorization.solve_unroll),SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true},SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true},Float64,Float64,Int64,Int64})
    Base.precompile(Tuple{typeof(LoopVectorization.stride_penalty),LoopVectorization.LoopSet,Vector{Symbol}})
    Base.precompile(Tuple{typeof(LoopVectorization.substitute_broadcast),Expr,Symbol})
    Base.precompile(Tuple{typeof(println),IOBuffer,Vector{LoopVectorization.Operation}})
    Base.precompile(Tuple{typeof(show),IOBuffer,Tuple{Int64, LoopVectorization.NumberType}})
end


text:
const _rainbowColors = RGBA{Float64}[colorant"purple", colorant"blue", colorant"green", colorant"orange", colorant"red"]
const _testColors = RGBA{Float64}[colorant"darkblue", colorant"blueviolet",  colorant"darkcyan",colorant"green",
                     darken(colorant"yellow",0.3), colorant"orange", darken(colorant"red",0.2)]

const _misc_color_lib = ColorLibrary(IdDict(:default => :sequential, :sequential => :heat, :diverging => :bluesreds), IdDict(
    :reds         => RGBA{Float64}[colorant"lightpink", colorant"darkred"],
    :greens       => RGBA{Float64}[colorant"lightgreen", colorant"darkgreen"],
    :redsblues    => RGBA{Float64}[colorant"darkred", RGB(0.8,0.85,0.8), colorant"darkblue"],
    :bluesreds    => RGBA{Float64}[colorant"darkblue", RGB(0.8,0.85,0.8), colorant"darkred"],
    :heat         => RGBA{Float64}[colorant"lightyellow", colorant"orange", colorant"darkred"],
    :grays        => RGBA{Float64}[RGB(.05,.05,.05),RGB(.95,.95,.95)],
    :rainbow      => _rainbowColors,
    :lightrainbow => map(lighten, _rainbowColors),
    :darkrainbow  => map(darken, _rainbowColors),
    :darktest     => _testColors,
    :lighttest    => map(c -> lighten(c, 0.3), _testColors),
))

register_color_library(:misc, _misc_color_lib)

register_gradient_colors(:cividis, sample_evenly(
[
    RGB(0.0000,0.1262,0.3015),
    RGB(0.0000,0.1292,0.3077),
    RGB(0.0000,0.1321,0.3142),
    RGB(0.0000,0.1350,0.3205),
    RGB(0.0000,0.1379,0.3269),
    RGB(0.0000,0.1408,0.3334),
    RGB(0.0000,0.1437,0.3400),
    RGB(0.0000,0.1465,0.3467),
    RGB(0.0000,0.1492,0.3537),
    RGB(0.0000,0.1519,0.3606),
    RGB(0.0000,0.1546,0.3676),
    RGB(0.0000,0.1574,0.3746),
    RGB(0.0000,0.1601,0.3817),
    RGB(0.0000,0.1629,0.3888),
    RGB(0.0000,0.1657,0.3960),
    RGB(0.0000,0.1685,0.4031),
    RGB(0.0000,0.1714,0.4102),
    RGB(0.0000,0.1743,0.4172),
    RGB(0.0000,0.1773,0.4241),
    RGB(0.0000,0.1798,0.4307),
    RGB(0.0000,0.1817,0.4347),
    RGB(0.0000,0.1834,0.4363),
    RGB(0.0000,0.1852,0.4368),
    RGB(0.0000,0.1872,0.4368),
    RGB(0.0000,0.1901,0.4365),
    RGB(0.0000,0.1930,0.4361),
    RGB(0.0000,0.1958,0.4356),
    RGB(0.0000,0.1987,0.4349),
    RGB(0.0000,0.2015,0.4343),
    RGB(0.0000,0.2044,0.4336),
    RGB(0.0000,0.2073,0.4329),
    RGB(0.0055,0.2101,0.4322),
    RGB(0.0236,0.2130,0.4314),
    RGB(0.0416,0.2158,0.4308),
    RGB(0.0576,0.2187,0.4301),
    RGB(0.0710,0.2215,0.4293),
    RGB(0.0827,0.2244,0.4287),
    RGB(0.0932,0.2272,0.4280),
    RGB(0.1030,0.2300,0.4274),
    RGB(0.1120,0.2329,0.4268),
    RGB(0.1204,0.2357,0.4262),
    RGB(0.1283,0.2385,0.4256),
    RGB(0.1359,0.2414,0.4251),
    RGB(0.1431,0.2442,0.4245),
    RGB(0.1500,0.2470,0.4241),
    RGB(0.1566,0.2498,0.4236),
    RGB(0.1630,0.2526,0.4232),
    RGB(0.1692,0.2555,0.4228),
    RGB(0.1752,0.2583,0.4224),
    RGB(0.1811,0.2611,0.4220),
    RGB(0.1868,0.2639,0.4217),
    RGB(0.1923,0.2667,0.4214),
    RGB(0.1977,0.2695,0.4212),
    RGB(0.2030,0.2723,0.4209),
    RGB(0.2082,0.2751,0.4207),
    RGB(0.2133,0.2780,0.4205),
    RGB(0.2183,0.2808,0.4204),
    RGB(0.2232,0.2836,0.4203),
    RGB(0.2281,0.2864,0.4202),
    RGB(0.2328,0.2892,0.4201),
    RGB(0.2375,0.2920,0.4200),
    RGB(0.2421,0.2948,0.4200),
    RGB(0.2466,0.2976,0.4200),
    RGB(0.2511,0.3004,0.4201),
    RGB(0.2556,0.3032,0.4201),
    RGB(0.2599,0.3060,0.4202),
    RGB(0.2643,0.3088,0.4203),
    RGB(0.2686,0.3116,0.4205),
    RGB(0.2728,0.3144,0.4206),
    RGB(0.2770,0.3172,0.4208),
    RGB(0.2811,0.3200,0.4210),
    RGB(0.2853,0.3228,0.4212),
    RGB(0.2894,0.3256,0.4215),
    RGB(0.2934,0.3284,0.4218),
    RGB(0.2974,0.3312,0.4221),
    RGB(0.3014,0.3340,0.4224),
    RGB(0.3054,0.3368,0.4227),
    RGB(0.3093,0.3396,0.4231),
    RGB(0.3132,0.3424,0.4236),
    RGB(0.3170,0.3453,0.4240),
    RGB(0.3209,0.3481,0.4244),
    RGB(0.3247,0.3509,0.4249),
    RGB(0.3285,0.3537,0.4254),
    RGB(0.3323,0.3565,0.4259),
    RGB(0.3361,0.3593,0.4264),
    RGB(0.3398,0.3622,0.4270),
    RGB(0.3435,0.3650,0.4276),
    RGB(0.3472,0.3678,0.4282),
    RGB(0.3509,0.3706,0.4288),
    RGB(0.3546,0.3734,0.4294),
    RGB(0.3582,0.3763,0.4302),
    RGB(0.3619,0.3791,0.4308),
    RGB(0.3655,0.3819,0.4316),
    RGB(0.3691,0.3848,0.4322),
    RGB(0.3727,0.3876,0.4331),
    RGB(0.3763,0.3904,0.4338),
    RGB(0.3798,0.3933,0.4346),
    RGB(0.3834,0.3961,0.4355),
    RGB(0.3869,0.3990,0.4364),
    RGB(0.3905,0.4018,0.4372),
    RGB(0.3940,0.4047,0.4381),
    RGB(0.3975,0.4075,0.4390),
    RGB(0.4010,0.4104,0.4400),
    RGB(0.4045,0.4132,0.4409),
    RGB(0.4080,0.4161,0.4419),
    RGB(0.4114,0.4189,0.4430),
    RGB(0.4149,0.4218,0.4440),
    RGB(0.4183,0.4247,0.4450),
    RGB(0.4218,0.4275,0.4462),
    RGB(0.4252,0.4304,0.4473),
    RGB(0.4286,0.4333,0.4485),
    RGB(0.4320,0.4362,0.4496),
    RGB(0.4354,0.4390,0.4508),
    RGB(0.4388,0.4419,0.4521),
    RGB(0.4422,0.4448,0.4534),
    RGB(0.4456,0.4477,0.4547),
    RGB(0.4489,0.4506,0.4561),
    RGB(0.4523,0.4535,0.4575),
    RGB(0.4556,0.4564,0.4589),
    RGB(0.4589,0.4593,0.4604),
    RGB(0.4622,0.4622,0.4620),
    RGB(0.4656,0.4651,0.4635),
    RGB(0.4689,0.4680,0.4650),
    RGB(0.4722,0.4709,0.4665),
    RGB(0.4756,0.4738,0.4679),
    RGB(0.4790,0.4767,0.4691),
    RGB(0.4825,0.4797,0.4701),
    RGB(0.4861,0.4826,0.4707),
    RGB(0.4897,0.4856,0.4714),
    RGB(0.4934,0.4886,0.4719),
    RGB(0.4971,0.4915,0.4723),
    RGB(0.5008,0.4945,0.4727),
    RGB(0.5045,0.4975,0.4730),
    RGB(0.5083,0.5005,0.4732),
    RGB(0.5121,0.5035,0.4734),
    RGB(0.5158,0.5065,0.4736),
    RGB(0.5196,0.5095,0.4737),
    RGB(0.5234,0.5125,0.4738),
    RGB(0.5272,0.5155,0.4739),
    RGB(0.5310,0.5186,0.4739),
    RGB(0.5349,0.5216,0.4738),
    RGB(0.5387,0.5246,0.4739),
    RGB(0.5425,0.5277,0.4738),
    RGB(0.5464,0.5307,0.4736),
    RGB(0.5502,0.5338,0.4735),
    RGB(0.5541,0.5368,0.4733),
    RGB(0.5579,0.5399,0.4732),
    RGB(0.5618,0.5430,0.4729),
    RGB(0.5657,0.5461,0.4727),
    RGB(0.5696,0.5491,0.4723),
    RGB(0.5735,0.5522,0.4720),
    RGB(0.5774,0.5553,0.4717),
    RGB(0.5813,0.5584,0.4714),
    RGB(0.5852,0.5615,0.4709),
    RGB(0.5892,0.5646,0.4705),
    RGB(0.5931,0.5678,0.4701),
    RGB(0.5970,0.5709,0.4696),
    RGB(0.6010,0.5740,0.4691),
    RGB(0.6050,0.5772,0.4685),
    RGB(0.6089,0.5803,0.4680),
    RGB(0.6129,0.5835,0.4673),
    RGB(0.6168,0.5866,0.4668),
    RGB(0.6208,0.5898,0.4662),
    RGB(0.6248,0.5929,0.4655),
    RGB(0.6288,0.5961,0.4649),
    RGB(0.6328,0.5993,0.4641),
    RGB(0.6368,0.6025,0.4632),
    RGB(0.6408,0.6057,0.4625),
    RGB(0.6449,0.6089,0.4617),
    RGB(0.6489,0.6121,0.4609),
    RGB(0.6529,0.6153,0.4600),
    RGB(0.6570,0.6185,0.4591),
    RGB(0.6610,0.6217,0.4583),
    RGB(0.6651,0.6250,0.4573),
    RGB(0.6691,0.6282,0.4562),
    RGB(0.6732,0.6315,0.4553),
    RGB(0.6773,0.6347,0.4543),
    RGB(0.6813,0.6380,0.4532),
    RGB(0.6854,0.6412,0.4521),
    RGB(0.6895,0.6445,0.4511),
    RGB(0.6936,0.6478,0.4499),
    RGB(0.6977,0.6511,0.4487),
    RGB(0.7018,0.6544,0.4475),
    RGB(0.7060,0.6577,0.4463),
    RGB(0.7101,0.6610,0.4450),
    RGB(0.7142,0.6643,0.4437),
    RGB(0.7184,0.6676,0.4424),
    RGB(0.7225,0.6710,0.4409),
    RGB(0.7267,0.6743,0.4396),
    RGB(0.7308,0.6776,0.4382),
    RGB(0.7350,0.6810,0.4368),
    RGB(0.7392,0.6844,0.4352),
    RGB(0.7434,0.6877,0.4338),
    RGB(0.7476,0.6911,0.4322),
    RGB(0.7518,0.6945,0.4307),
    RGB(0.7560,0.6979,0.4290),
    RGB(0.7602,0.7013,0.4273),
    RGB(0.7644,0.7047,0.4258),
    RGB(0.7686,0.7081,0.4241),
    RGB(0.7729,0.7115,0.4223),
    RGB(0.7771,0.7150,0.4205),
    RGB(0.7814,0.7184,0.4188),
    RGB(0.7856,0.7218,0.4168),
    RGB(0.7899,0.7253,0.4150),
    RGB(0.7942,0.7288,0.4129),
    RGB(0.7985,0.7322,0.4111),
    RGB(0.8027,0.7357,0.4090),
    RGB(0.8070,0.7392,0.4070),
    RGB(0.8114,0.7427,0.4049),
    RGB(0.8157,0.7462,0.4028),
    RGB(0.8200,0.7497,0.4007),
    RGB(0.8243,0.7532,0.3984),
    RGB(0.8287,0.7568,0.3961),
    RGB(0.8330,0.7603,0.3938),
    RGB(0.8374,0.7639,0.3915),
    RGB(0.8417,0.7674,0.3892),
    RGB(0.8461,0.7710,0.3869),
    RGB(0.8505,0.7745,0.3843),
    RGB(0.8548,0.7781,0.3818),
    RGB(0.8592,0.7817,0.3793),
    RGB(0.8636,0.7853,0.3766),
    RGB(0.8681,0.7889,0.3739),
    RGB(0.8725,0.7926,0.3712),
    RGB(0.8769,0.7962,0.3684),
    RGB(0.8813,0.7998,0.3657),
    RGB(0.8858,0.8035,0.3627),
    RGB(0.8902,0.8071,0.3599),
    RGB(0.8947,0.8108,0.3569),
    RGB(0.8992,0.8145,0.3538),
    RGB(0.9037,0.8182,0.3507),
    RGB(0.9082,0.8219,0.3474),
    RGB(0.9127,0.8256,0.3442),
    RGB(0.9172,0.8293,0.3409),
    RGB(0.9217,0.8330,0.3374),
    RGB(0.9262,0.8367,0.3340),
    RGB(0.9308,0.8405,0.3306),
    RGB(0.9353,0.8442,0.3268),
    RGB(0.9399,0.8480,0.3232),
    RGB(0.9444,0.8518,0.3195),
    RGB(0.9490,0.8556,0.3155),
    RGB(0.9536,0.8593,0.3116),
    RGB(0.9582,0.8632,0.3076),
    RGB(0.9628,0.8670,0.3034),
    RGB(0.9674,0.8708,0.2990),
    RGB(0.9721,0.8746,0.2947),
    RGB(0.9767,0.8785,0.2901),
    RGB(0.9814,0.8823,0.2856),
    RGB(0.9860,0.8862,0.2807),
    RGB(0.9907,0.8901,0.2759),
    RGB(0.9954,0.8940,0.2708),
    RGB(1.0000,0.8979,0.2655),
    RGB(1.0000,0.9018,0.2600),
    RGB(1.0000,0.9057,0.2593),
    RGB(1.0000,0.9094,0.2634),
    RGB(1.0000,0.9131,0.2680),
    RGB(1.0000,0.9169,0.2731)
], 64), :misc)

ext: jl


text:
"""
		mlp_block(inplanes::Integer, hidden_planes::Integer, outplanes::Integer = inplanes;
							dropout = 0., activation = gelu)

Feedforward block used in many MLPMixer-like and vision-transformer models.

# Arguments
- `inplanes`: Number of dimensions in the input.
- `hidden_planes`: Number of dimensions in the intermediate layer.
- `outplanes`: Number of dimensions in the output - by default it is the same as `inplanes`.
- `dropout`: Dropout rate.
- `activation`: Activation function to use.
"""
function mlp_block(inplanes::Integer, hidden_planes::Integer, outplanes::Integer = inplanes;
									 dropout = 0., activation = gelu)
  Chain(Dense(inplanes, hidden_planes, activation), Dropout(dropout),
        Dense(hidden_planes, outplanes), Dropout(dropout))
end

"""
    gated_mlp(gate_layer, inplanes::Integer, hidden_planes::Integer,
							outplanes::Integer = inplanes; dropout = 0., activation = gelu)

Feedforward block based on the implementation in the paper "Pay Attention to MLPs".
([reference](https://arxiv.org/abs/2105.08050))

# Arguments
- `gate_layer`: Layer to use for the gating.
- `inplanes`: Number of dimensions in the input.
- `hidden_planes`: Number of dimensions in the intermediate layer.
- `outplanes`: Number of dimensions in the output - by default it is the same as `inplanes`.
- `dropout`: Dropout rate.
- `activation`: Activation function to use.
"""
function gated_mlp_block(gate_layer, inplanes::Integer, hidden_planes::Integer,
												 outplanes::Integer = inplanes; dropout = 0., activation = gelu)
	@assert hidden_planes % 2 == 0 "`hidden_planes` must be even for gated MLP"
	return Chain(Dense(inplanes, hidden_planes, activation),
							 Dropout(dropout),
							 gate_layer(hidden_planes),
							 Dense(hidden_planes ÷ 2, outplanes),
							 Dropout(dropout))
end
gated_mlp_block(::typeof(identity), args...; kwargs...) = mlp_block(args...; kwargs...)


text:
module particle

export PositionBasis, MomentumBasis,
        gaussianstate,
        spacing, samplepoints,
        position, momentum, potentialoperator,
        transform

import Base: ==, position
import ..operators

using ..bases, ..states, ..operators, ..operators_dense, ..operators_sparse

"""
    PositionBasis(xmin, xmax, Npoints)
    PositionBasis(b::MomentumBasis)

Basis for a particle in real space.

For simplicity periodic boundaries are assumed which means that
the rightmost point defined by `xmax` is not included in the basis
but is defined to be the same as `xmin`.

When a [`MomentumBasis`](@ref) is given as argument the exact values
of ``x_{min}`` and ``x_{max}`` are due to the periodic boundary conditions
more or less arbitrary and are chosen to be
``-\\pi/dp`` and ``\\pi/dp`` with ``dp=(p_{max}-p_{min})/N``.
"""
mutable struct PositionBasis <: Basis
    shape::Vector{Int}
    xmin::Float64
    xmax::Float64
    N::Int
    PositionBasis(xmin::Real, xmax::Real, N::Int) = new([N], xmin, xmax, N)
end

"""
    MomentumBasis(pmin, pmax, Npoints)
    MomentumBasis(b::PositionBasis)

Basis for a particle in momentum space.

For simplicity periodic boundaries are assumed which means that
`pmax` is not included in the basis but is defined to be the same as `pmin`.

When a [`PositionBasis`](@ref) is given as argument the exact values
of ``p_{min}`` and ``p_{max}`` are due to the periodic boundary conditions
more or less arbitrary and are chosen to be
``-\\pi/dx`` and ``\\pi/dx`` with ``dx=(x_{max}-x_{min})/N``.
"""
mutable struct MomentumBasis <: Basis
    shape::Vector{Int}
    pmin::Float64
    pmax::Float64
    N::Int
    MomentumBasis(pmin::Real, pmax::Real, N::Int) = new([N], pmin, pmax, N)
end

PositionBasis(b::MomentumBasis) = (dp = (b.pmax - b.pmin)/b.N; PositionBasis(-pi/dp, pi/dp, b.N))
MomentumBasis(b::PositionBasis) = (dx = (b.xmax - b.xmin)/b.N; MomentumBasis(-pi/dx, pi/dx, b.N))

==(b1::PositionBasis, b2::PositionBasis) = b1.xmin==b2.xmin && b1.xmax==b2.xmax && b1.N==b2.N
==(b1::MomentumBasis, b2::MomentumBasis) = b1.pmin==b2.pmin && b1.pmax==b2.pmax && b1.N==b2.N

"""
    gaussianstate(b::PositionBasis, x0, p0, sigma)
    gaussianstate(b::MomentumBasis, x0, p0, sigma)

Create a Gaussian state around `x0` and` p0` with width `sigma`.

In real space the gaussian state is defined as

```math
\\Psi(x) = \\frac{\\sqrt{\\Delta x}}{\\pi^{1/4}\\sqrt{\\sigma}}
            e^{i p_0 (x-\\frac{x_0}{2}) - \\frac{(x-x_0)^2}{2 \\sigma^2}}
```

and is connected to the momentum space definition

```math
\\Psi(p) = \\frac{\\sqrt{\\sigma} \\sqrt{\\Delta x}}{\\pi^{1/4}}
            e^{-i x_0 (p-\\frac{p_0}{2}) - \\frac{1}{2}(p-p_0)^2 \\sigma^2}
```

via a Fourier-transformation

```math
\\Psi(p) = \\frac{1}{\\sqrt{2\\pi}}
            \\int_{-\\infty}^{\\infty} e^{-ipx}\\Psi(x) \\mathrm{d}x
```

The state has the properties

* ``⟨p⟩ = p_0``
* ``⟨x⟩ = x_0``
* ``\\mathrm{Var}(x) = \\frac{σ^2}{2}``
* ``\\mathrm{Var}(p) = \\frac{1}{2 σ^2}``

Due to the numerically necessary discretization additional scaling
factora ``\\sqrt{Δx}`` and ``\\sqrt{Δp}`` are used so that
``Ψx_i = \\sqrt{Δ x} Ψ(x_i)`` and ``Ψp_i = \\sqrt{Δ x} Ψ(p_i)`` so
that the resulting Ket state is normalized.
"""
function gaussianstate(b::PositionBasis, x0::Real, p0::Real, sigma::Real)
    psi = Ket(b)
    dx = spacing(b)
    alpha = 1./(pi^(1/4)*sqrt(sigma))*sqrt(dx)
    x = b.xmin
    for i=1:b.N
        psi.data[i] = alpha*exp(1im*p0*(x-x0/2) - (x-x0)^2/(2*sigma^2))
        x += dx
    end
    return psi
end

function gaussianstate(b::MomentumBasis, x0::Real, p0::Real, sigma::Real)
    psi = Ket(b)
    dp = spacing(b)
    alpha = sqrt(sigma)/pi^(1/4)*sqrt(dp)
    p = b.pmin
    for i=1:b.N
        psi.data[i] = alpha*exp(-1im*x0*(p-p0/2) - (p-p0)^2/2*sigma^2)
        p += dp
    end
    return psi
end

"""
    spacing(b::PositionBasis)

Difference between two adjacent points of the real space basis.
"""
spacing(b::PositionBasis) = (b.xmax - b.xmin)/b.N
"""
    spacing(b::MomentumBasis)

Momentum difference between two adjacent points of the momentum basis.
"""
spacing(b::MomentumBasis) = (b.pmax - b.pmin)/b.N

"""
    samplepoints(b::PositionBasis)

x values of the real space basis.
"""
samplepoints(b::PositionBasis) = (dx = spacing(b); Float64[b.xmin + i*dx for i=0:b.N-1])
"""
    samplepoints(b::MomentumBasis)

p values of the momentum basis.
"""
samplepoints(b::MomentumBasis) = (dp = spacing(b); Float64[b.pmin + i*dp for i=0:b.N-1])

"""
    position(b::PositionBasis)

Position operator in real space.
"""
position(b::PositionBasis) = SparseOperator(b, spdiagm(complex(samplepoints(b)), 0, length(b), length(b)))

"""
    position(b:MomentumBasis)

Position operator in momentum space.
"""
function position(b::MomentumBasis)
    b_pos = PositionBasis(b)
    transform(b, b_pos)*full(position(b_pos))*transform(b_pos, b)
end

"""
    momentum(b:MomentumBasis)

Momentum operator in momentum space.
"""
momentum(b::MomentumBasis) = SparseOperator(b, spdiagm(complex(samplepoints(b)), 0, length(b), length(b)))

"""
    momentum(b::PositionBasis)

Momentum operator in real space.
"""
function momentum(b::PositionBasis)
    b_mom = MomentumBasis(b)
    transform(b, b_mom)*full(momentum(b_mom))*transform(b_mom, b)
end

"""
    potentialoperator(b::PositionBasis, V(x))

Operator representing a potential ``V(x)`` in real space.
"""
function potentialoperator(b::PositionBasis, V::Function)
    x = samplepoints(b)
    diagonaloperator(b, V.(x))
end

"""
    potentialoperator(b::MomentumBasis, V(x))

Operator representing a potential ``V(x)`` in momentum space.
"""
function potentialoperator(b::MomentumBasis, V::Function)
    b_pos = PositionBasis(b)
    transform(b, b_pos)*full(potentialoperator(b_pos, V))*transform(b_pos, b)
end

"""
    FFTOperator

Abstract type for all implementations of FFT operators.
"""
abstract type FFTOperator <: Operator end

PlanFFT = Base.DFT.FFTW.cFFTWPlan

"""
    FFTOperators

Operator performing a fast fourier transformation when multiplied with a state
that is a Ket or an Operator.
"""
mutable struct FFTOperators <: FFTOperator
    basis_l::Basis
    basis_r::Basis
    fft_l!::PlanFFT
    fft_r!::PlanFFT
    fft_l2!::PlanFFT
    fft_r2!::PlanFFT
    mul_before::Array{Complex128}
    mul_after::Array{Complex128}
end

"""
    FFTKets

Operator that can only perform fast fourier transformations on Kets.
This is much more memory efficient when only working with Kets.
"""
mutable struct FFTKets <: FFTOperator
    basis_l::Basis
    basis_r::Basis
    fft_l!::PlanFFT
    fft_r!::PlanFFT
    mul_before::Array{Complex128}
    mul_after::Array{Complex128}
end

"""
    transform(b1::MomentumBasis, b2::PositionBasis)
    transform(b1::PositionBasis, b2::MomentumBasis)

Transformation operator between position basis and momentum basis.
"""
function transform(basis_l::MomentumBasis, basis_r::PositionBasis; ket_only::Bool=false)
    Lx = (basis_r.xmax - basis_r.xmin)
    dp = spacing(basis_l)
    dx = spacing(basis_r)
    if basis_l.N != basis_r.N || abs(2*pi/dp - Lx)/Lx > 1e-12
        throw(IncompatibleBases())
    end
    mul_before = exp.(-1im*basis_l.pmin*(samplepoints(basis_r)-basis_r.xmin))
    mul_after = exp.(-1im*basis_r.xmin*samplepoints(basis_l))/sqrt(basis_r.N)
    x = Vector{Complex128}(length(basis_r))
    if ket_only
        FFTKets(basis_l, basis_r, plan_bfft!(x), plan_fft!(x), mul_before, mul_after)
    else
        A = Matrix{Complex128}(length(basis_r), length(basis_r))
        FFTOperators(basis_l, basis_r, plan_bfft!(x), plan_fft!(x), plan_bfft!(A, 2), plan_fft!(A, 1), mul_before, mul_after)
    end
end

"""
    transform(b1::CompositeBasis, b2::CompositeBasis)

Transformation operator between two composite bases. Each of the bases
has to contain bases of type PositionBasis and the other one a corresponding
MomentumBasis.
"""
function transform(basis_l::PositionBasis, basis_r::MomentumBasis; ket_only::Bool=false)
    Lx = (basis_l.xmax - basis_l.xmin)
    dp = spacing(basis_r)
    dx = spacing(basis_l)
    if basis_l.N != basis_r.N || abs(2*pi/dp - Lx)/Lx > 1e-12
        throw(IncompatibleBases())
    end
    mul_before = exp.(1im*basis_l.xmin*(samplepoints(basis_r)-basis_r.pmin))
    mul_after = exp.(1im*basis_r.pmin*samplepoints(basis_l))/sqrt(basis_r.N)
    x = Vector{Complex128}(length(basis_r))
    if ket_only
        FFTKets(basis_l, basis_r, plan_fft!(x), plan_bfft!(x), mul_before, mul_after)
    else
        A = Matrix{Complex128}(length(basis_r), length(basis_r))
        FFTOperators(basis_l, basis_r, plan_fft!(x), plan_bfft!(x), plan_fft!(A, 2), plan_bfft!(A, 1), mul_before, mul_after)
    end
end

function transform(basis_l::CompositeBasis, basis_r::CompositeBasis; ket_only::Bool=false, index::Vector{Int}=Int[])
    @assert length(basis_l.bases) == length(basis_r.bases)
    if length(index) == 0
        check_pos = typeof.(basis_l.bases) .== PositionBasis
        check_mom = typeof.(basis_l.bases) .== MomentumBasis
        if any(check_pos) && !any(check_mom)
            index = [1:length(basis_l.bases);][check_pos]
        elseif any(check_mom) && !any(check_pos)
            index = [1:length(basis_l.bases);][check_mom]
        else
            throw(IncompatibleBases())
        end
    end
    if all(typeof.(basis_l.bases[index]) .== PositionBasis)
        @assert all(typeof.(basis_r.bases[index]) .== MomentumBasis)
        transform_xp(basis_l, basis_r, index; ket_only=ket_only)
    elseif all(typeof.(basis_l.bases[index]) .== MomentumBasis)
        @assert all(typeof.(basis_r.bases[index]) .== PositionBasis)
        transform_px(basis_l, basis_r, index; ket_only=ket_only)
    else
        throw(IncompatibleBases())
    end
end

function transform_xp(basis_l::CompositeBasis, basis_r::CompositeBasis, index::Vector{Int}; ket_only::Bool=false)
    n = length(basis_l.bases)
    Lx = [(b.xmax - b.xmin) for b=basis_l.bases[index]]
    dp = [spacing(b) for b=basis_r.bases[index]]
    dx = [spacing(b) for b=basis_l.bases[index]]
    N = [length(b) for b=basis_l.bases]
    for i=1:n
        if N[i] != length(basis_r.bases[i])
            throw(IncompatibleBases())
        end
    end
    for i=1:length(index)
        if abs(2*pi/dp[i] - Lx[i])/Lx[i] > 1e-12
            throw(IncompatibleBases())
        end
    end

    if index[1] == 1
        mul_before = exp.(1im*basis_l.bases[1].xmin*(samplepoints(basis_r.bases[1])-basis_r.bases[1].pmin))
        mul_after = exp.(1im*basis_r.bases[1].pmin*samplepoints(basis_l.bases[1]))/sqrt(basis_r.bases[1].N)
    else
        mul_before = ones(N[1])
        mul_after = ones(N[1])
    end
    for i=2:n
        if any(i .== index)
            mul_before = kron(exp.(1im*basis_l.bases[i].xmin*(samplepoints(basis_r.bases[i])-basis_r.bases[i].pmin)), mul_before)
            mul_after = kron(exp.(1im*basis_r.bases[i].pmin*samplepoints(basis_l.bases[i]))/sqrt(basis_r.bases[i].N), mul_after)
        else
            mul_before = kron(ones(N[i]), mul_before)
            mul_after = kron(ones(N[i]), mul_after)
        end
    end
    mul_before = reshape(mul_before, (N...))
    mul_after = reshape(mul_after, (N...))

    x = Array{Complex128}(N...)
    if ket_only
        FFTKets(basis_l, basis_r, plan_fft!(x, index), plan_bfft!(x, index), mul_before, mul_after)
    else
        A = Array{Complex128}([N, N;]...)
        FFTOperators(basis_l, basis_r, plan_fft!(x, index), plan_bfft!(x, index), plan_fft!(A, [n + 1:2n;][index]), plan_bfft!(A, [1:n;][index]), mul_before, mul_after)
    end
end

function transform_px(basis_l::CompositeBasis, basis_r::CompositeBasis, index::Vector{Int}; ket_only::Bool=false)
    n = length(basis_l.bases)
    Lx = [(b.xmax - b.xmin) for b=basis_r.bases[index]]
    dp = [spacing(b) for b=basis_l.bases[index]]
    dx = [spacing(b) for b=basis_r.bases[index]]
    N = [length(b) for b=basis_l.bases]
    for i=1:n
        if N[i] != length(basis_r.bases[i])
            throw(IncompatibleBases())
        end
    end
    for i=1:length(index)
        if abs(2*pi/dp[i] - Lx[i])/Lx[i] > 1e-12
            throw(IncompatibleBases())
        end
    end

    if index[1] == 1
        mul_before = exp.(-1im*basis_l.bases[1].pmin*(samplepoints(basis_r.bases[1])-basis_r.bases[1].xmin))
        mul_after = exp.(-1im*basis_r.bases[1].xmin*samplepoints(basis_l.bases[1]))/sqrt(N[1])
    else
        mul_before = ones(N[1])
        mul_after = ones(N[1])
    end
    for i=2:n
        if i in index
            mul_before = kron(exp.(-1im*basis_l.bases[i].pmin*(samplepoints(basis_r.bases[i])-basis_r.bases[i].xmin)), mul_before)
            mul_after = kron(exp.(-1im*basis_r.bases[i].xmin*samplepoints(basis_l.bases[i]))/sqrt(N[i]), mul_after)
        else
            mul_before = kron(ones(N[i]), mul_before)
            mul_after = kron(ones(N[i]), mul_after)
        end
    end
    mul_before = reshape(mul_before, (N...))
    mul_after = reshape(mul_after, (N...))

    x = Array{Complex128}(N...)
    if ket_only
        FFTKets(basis_l, basis_r, plan_bfft!(x, index), plan_fft!(x, index), mul_before, mul_after)
    else
        A = Array{Complex128}([N, N;]...)
        FFTOperators(basis_l, basis_r, plan_bfft!(x, index), plan_fft!(x, index), plan_bfft!(A, [n + 1:2n;][index]), plan_fft!(A, [1:n;][index]), mul_before, mul_after)
    end
end

operators.full(op::FFTOperators) = op*identityoperator(DenseOperator, op.basis_r)

operators.dagger(op::FFTOperators) = transform(op.basis_r, op.basis_l)
operators.dagger(op::FFTKets) = transform(op.basis_r, op.basis_l; ket_only=true)

operators.tensor(A::FFTOperators, B::FFTOperators) = transform(tensor(A.basis_l, B.basis_l), tensor(A.basis_r, B.basis_r))
operators.tensor(A::FFTKets, B::FFTKets) = transform(tensor(A.basis_l, B.basis_l), tensor(A.basis_r, B.basis_r); ket_only=true)

function operators.gemv!(alpha_, M::FFTOperator, b::Ket, beta_, result::Ket)
    alpha = convert(Complex128, alpha_)
    beta = convert(Complex128, beta_)
    N::Int = length(M.basis_r)
    if beta==Complex(0.)
        @inbounds for i=1:N
            result.data[i] = M.mul_before[i] * b.data[i]
        end
        M.fft_r! * reshape(result.data, size(M.mul_before))
        @inbounds for i=1:N
            result.data[i] *= M.mul_after[i] * alpha
        end
    else
        psi_ = Ket(M.basis_l, copy(b.data))
        @inbounds for i=1:N
            psi_.data[i] *= M.mul_before[i]
        end
        M.fft_r! * reshape(psi_.data, size(M.mul_before))
        @inbounds for i=1:N
            result.data[i] = beta*result.data[i] + alpha * psi_.data[i] * M.mul_after[i]
        end
    end
    nothing
end

function operators.gemv!(alpha_, b::Bra, M::FFTOperator, beta_, result::Bra)
    alpha = convert(Complex128, alpha_)
    beta = convert(Complex128, beta_)
    N::Int = length(M.basis_l)
    if beta==Complex(0.)
        @inbounds for i=1:N
            result.data[i] = conj(M.mul_after[i]) * conj(b.data[i])
        end
        M.fft_l! * reshape(result.data, size(M.mul_after))
        @inbounds for i=1:N
            result.data[i] = conj(result.data[i]) * M.mul_before[i] * alpha
        end
    else
        psi_ = Bra(M.basis_r, conj(b.data))
        @inbounds for i=1:N
            psi_.data[i] *= conj(M.mul_after[i])
        end
        M.fft_l! * reshape(psi_.data, size(M.mul_after))
        @inbounds for i=1:N
            result.data[i] = beta*result.data[i] + alpha * conj(psi_.data[i]) * M.mul_before[i]
        end
    end
    nothing
end

function operators.gemm!(alpha_, A::DenseOperator, B::FFTOperators, beta_, result::DenseOperator)
    alpha = convert(Complex128, alpha_)
    beta = convert(Complex128, beta_)
    if beta != Complex(0.)
        data = Matrix{Complex128}(size(result.data, 1), size(result.data, 2))
    else
        data = result.data
    end
    copy!(data, A.data)
    scale!(data, B.mul_after[:])
    conj!(data)
    B.fft_l2! * reshape(data, [size(B.mul_after)..., size(B.mul_after)...;]...)
    conj!(data)
    scale!(data, B.mul_before[:])
    if alpha != Complex(1.)
        scale!(alpha, data)
    end
    if beta != Complex(0.)
        scale!(result.data, beta)
        result.data += data
    end
    nothing
end

function operators.gemm!(alpha_, A::FFTOperators, B::DenseOperator, beta_, result::DenseOperator)
    alpha = convert(Complex128, alpha_)
    beta = convert(Complex128, beta_)
    if beta != Complex(0.)
        data = Matrix{Complex128}(size(result.data, 1), size(result.data, 2))
    else
        data = result.data
    end
    copy!(data, B.data)
    scale!(A.mul_before[:], data)
    A.fft_r2! * reshape(data, [size(A.mul_before)..., size(A.mul_before)...;]...)
    scale!(A.mul_after[:], data)
    if alpha != Complex(1.)
        scale!(alpha, data)
    end
    if beta != Complex(0.)
        scale!(result.data, beta)
        result.data += data
    end
    nothing
end

end # module