text:
using Elefridge
using ZfpCompression
using PyPlot
using PyCall
using Blosc
using ColorSchemes
using StatsBase, Statistics
using TranscodingStreams, CodecZstd
xr = pyimport("xarray")
ccrs = pyimport("cartopy.crs")

ZstdCompressorL3 = ZstdCompressor(level=3)
TranscodingStreams.initialize(ZstdCompressorL3)

ZstdCompressorL10 = ZstdCompressor(level=3)
TranscodingStreams.initialize(ZstdCompressorL10)

ZstdCompressorL22 = ZstdCompressor(level=3)
TranscodingStreams.initialize(ZstdCompressorL22)

path = "/Users/milan/cams/gridded/"
filelist = filter(x->endswith(x,"_go3.grib"),readdir(path))
grib = xr.open_dataarray(joinpath(path,filelist[end]),engine="cfgrib")
X = grib.data
lat = grib.latitude.data
lon = grib.longitude.data

level = 85
X = permutedims(X,[3,2,1])      # for lossless longitude first
o3 = X[:,:,level]

## compression round + LZ4HC/Zstd
rbits_ll = [23,7,5,3,1,0]
cfs_ll = fill(0.0,length(rbits_ll))     # compression factors
decerr_ll = fill(0.0,length(rbits_ll))
# Blosc.set_compressor("lz4hc")

for (i,r) in enumerate(rbits_ll)
    Xr = round(X[1:450,:,:],r)      # use only half the longitudes for speed-up

    Xr8 = unsafe_wrap(Array, Ptr{UInt8}(pointer(Xr)), sizeof(Xr))
    Xc = transcode(ZstdCompressorL22,Xr8)

    # Xc = compress(Xr,level=5)
    cfs_ll[i] = sizeof(Xr)/sizeof(Xc)*2/(900*451)*(348528)    # relative to Float64 and unstructured grid
    # decerr_ll[i] = maximum(vec(abs.(log2.(abs.(X./o3r)))))
end

## compression zfp
X = grib.data
# X = permutedims(X,[3,2,1])      # for zfp vert x lat x lon
rbits_zfp = [23,13,10,8,5,4]               # corresponds to the sigbits from above
rbits_zfp = [23,1e-1,5e-1,4,20,50]
cfs_zfp = fill(0.0,length(rbits_zfp))     # compression factors
# decerr_zfp = fill(0.0,length(rbits_zfp))   # decimal error

O3plot = fill(0f0,length(rbits_zfp)-1,size(o3)...)

# zfp lossless, relative to Float64 and unstructured grid
cfs_zfp[1] = sizeof(X)/sizeof(zfp_compress(X))*2/(900*451)*(348528)

for (i,r) in enumerate(rbits_zfp[2:end])
    Xc = zfp_compress(log.(X),tol=r)
    cfs_zfp[i+1] = sizeof(X)/sizeof(Xc)*2/(900*451)*(348528)
    O3plot[i,:,:] = exp.(zfp_decompress(Xc)[level,:,:]')
    # local o3r = zfp_decompress(o3c)
    # decerr_zfp[i+1] = median(vec(abs.(log2.(abs.(X./o3r)))))
end

lat_div(n::Integer) = Array(-90:180/(n-1):90)
lon_div(n::Integer) = Array(0:360/n:360)

## PLOT
# pygui(true)
# ion()
ioff()
cmap = ColorMap(ColorSchemes.turku.colors)
vmin,vmax = extrema(o3)
fig,(ax0,ax1) = plt.subplots(2,1,figsize=(7,8),sharex=true,sharey=true,
                    subplot_kw=Dict("projection"=>ccrs.Robinson(central_longitude=180)))

pos = ax1.get_position()
cax = fig.add_axes([pos.x0,0.07,pos.x1-pos.x0,0.02])

lond = lon_div(length(rbits_ll))

# Float32 / lossless
scale = 1e7
pcm = ax0.pcolormesh(lon[lon .<= lond[2]],lat,scale*o3[lon .<= lond[2],:]',
    vmin=scale*vmin,vmax=scale*vmax,transform=ccrs.PlateCarree(),cmap=cmap)
ax1.pcolormesh(lon[lon .<= lond[2]],lat,o3[lon .<= lond[2],:]',
    vmin=vmin,vmax=vmax,transform=ccrs.PlateCarree(),cmap=cmap)

cbar = colorbar(pcm,cax=cax,orientation="horizontal")
cbar.set_label(L"mixing ratio [10$^{-7}$ kg/kg]")
cbar.set_ticks([1,2,3])

# round + lossless
for i in 3:length(lond)
    lonwhere = (lon .<= lond[i]) .& (lon .> lond[i-1]-0.5)
    ax0.pcolormesh(lon[lonwhere],lat,round(o3,rbits_ll[i-1])[lonwhere,:]',
        vmin=vmin,vmax=vmax,transform=ccrs.PlateCarree(),cmap=cmap)
end

# zfp
for i in 3:length(lond)
    lonwhere = (lon .<= lond[i]) .& (lon .> lond[i-1]-0.5)
    ax1.pcolormesh(lon[lonwhere],lat,O3plot[i-2,lonwhere,:]',
        vmin=vmin,vmax=vmax,transform=ccrs.PlateCarree(),cmap=cmap)
end

ax0.coastlines(color="k",lw=0.5)
ax1.coastlines(color="k",lw=0.5)

for longi in lond[2:end-1]
    ax0.plot(repeat([longi,],31),lat_div(31),"w",transform=ccrs.PlateCarree(),lw=0.25)
    ax1.plot(repeat([longi,],31),lat_div(31),"w",transform=ccrs.PlateCarree(),lw=0.25)
end

ax0.text(0.02,0.02,"Compression\n               factor",transform=ax0.transAxes,fontsize=8)
ax0.text(0.02,0.93,"         Significant\nbits retained",transform=ax0.transAxes,fontsize=8)
ax0.text(0.88,0.87,"Real information\n       preserved",transform=ax0.transAxes,fontsize=8)

info_preserved = ["100%","99.9%","99%","95%","82%","71%"]

for (i,x) in enumerate([0.16,0.28,0.41,0.54,0.65,0.77])
    ax0.text(x,0.89,info_preserved[i],transform=ax0.transAxes,fontweight="bold",
                fontsize=10,color="white")
    ax1.text(x,0.89,info_preserved[i],transform=ax1.transAxes,fontweight="bold",
                fontsize=10,color="white")
end

for (i,x) in enumerate([0.23,0.33,0.43,0.53,0.63,0.74])
    c = Int(round(cfs_ll[i]))
    r = rbits_ll[i]
    ax0.text(x,0.03,"$(c)x",transform=ax0.transAxes,fontweight="bold",fontsize=12,
                            color="white")
    ax0.text(x,0.98,"$(r)",transform=ax0.transAxes,fontweight="bold",fontsize=10,
                            color="white",va="top")
    ax1.text(x,0.98,"$(r)",transform=ax1.transAxes,fontweight="bold",fontsize=10,
                            color="white",va="top")
    c = Int(round(cfs_zfp[i]))
    ax1.text(x,0.03,"$(c)x",transform=ax1.transAxes,fontweight="bold",fontsize=12,
                            color="white")
end

ax0.set_title(L"O$_3$ compression, round+lossless")
ax0.set_title("a",fontweight="bold",loc="right")

ax1.set_title(L"O$_3$ zfp compression")
ax1.set_title("b",fontweight="bold",loc="right")

tight_layout(rect=[0.02,0.08,1,1])
savefig("/Users/milan/git/Elefridge.jl/maps/zfp_lossless_logo3_$level.png",dpi=200)
close(fig)


text:
# Delta is not following the PreMetric rules since d(x, x) == 1
struct Delta <: Distances.UnionPreMetric end

(dist::Delta)(a::Number, b::Number) = a == b
Base.@propagate_inbounds function (dist::Delta)(
    a::AbstractArray{<:Number}, b::AbstractArray{<:Number}
)
    @boundscheck if length(a) != length(b)
        throw(
            DimensionMismatch(
                "first array has length $(length(a)) which does not match the length of the second, $(length(b)).",
            ),
        )
    end
    return a == b
end

Distances.result_type(::Delta, Ta::Type, Tb::Type) = Bool

ext: jl


text:
const MISCELLANEOUS_JULIA_RELATED_REPOS = MirrorUpdater.Types.SrcDestPair[
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/libuv/libuv",
        destination_repo_name = "libuv-libuv",),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/JeffBezanson/femtolisp",
        destination_repo_name = "JeffBezanson-femtolisp",),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/JuliaMath/openlibm",
        destination_repo_name = "JuliaMath-openlibm",),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/JuliaMath/openspecfun",
        destination_repo_name = "JuliaMath-openspecfun",),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/JuliaStrings/utf8proc",
        destination_repo_name = "JuliaStrings-utf8proc",),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/djsegal/julia_observer",
        destination_repo_name = "djsegal-julia_observer",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/JuliaInterop/libcxxwrap-julia",
        destination_repo_name = "JuliaInterop-libcxxwrap-julia",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/JuliaStats/RmathBuilder",
        destination_repo_name = "JuliaStats-RmathBuilder",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/staticfloat/homebrew-julia",
        destination_repo_name = "staticfloat-homebrew-julia",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/staticfloat/homebrew-juliadeps",
        destination_repo_name = "staticfloat-homebrew-juliadeps",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/staticfloat/julia-docker",
        destination_repo_name = "staticfloat-julia-docker",
        ),
    ]


text:
using ExodusII
using Base.Test

# write your own tests here

@testset "ExodusII: $(num_dim)D" for num_dim in 1:3
  # Start by writing a new file
  filename = "TEMP" # tempname()
  dbtitle = "ExodusII.jl test title"
  num_nodes = 6
  num_elem = 5
  num_elem_blk = 2
  num_node_sets = 2
  num_side_sets = 0
  f = ExodusII.ex_create(filename,
                         dbtitle,
                         num_dim,
                         num_nodes,
                         num_elem,
                         num_elem_blk,
                         num_node_sets,
                         num_side_sets)
  # Set Coords and int maps
  coord_names = ["First C","Second C","Third C"][1:num_dim]
  ExodusII.ex_put_coord_names(f,coord_names)
  coord = [rand() for _1 in 1:num_dim, _2 in 1:num_nodes]
  ExodusII.ex_put_coord(f,coord)

  node_map = [2,5,3,4,9,10]
  ExodusII.ex_put_node_num_map(f,node_map)
  elem_map = [19,15,13,14,19]
  ExodusII.ex_put_elem_num_map(f,elem_map)

  # Set elem blocks
  block_id1 = 1234
  elem_type1 = "EDGE2"
  blk_num_elem1 = 4
  nodes_per_elem1 = 2
  el_attr1 = 0
  ExodusII.ex_put_elem_block(f,block_id1,elem_type1,blk_num_elem1,nodes_per_elem1,el_attr1)
  conn1 = [[2,4] [2,9] [3,10] [4,5]]
  ExodusII.ex_put_elem_connections(f,block_id1,conn1)

  block_id2 = 41923
  elem_type2 = "EDGE2"
  blk_num_elem2 = 3
  nodes_per_elem2 = 2
  el_attr2 = 0
  ExodusII.ex_put_elem_block(f,block_id2,elem_type2,blk_num_elem2,nodes_per_elem2,el_attr2)
  conn2 = [[5,2] [2,3] [5,9]]
  ExodusII.ex_put_elem_connections(f,block_id2,conn2)

  # Set node sets
  set_id1 = 4123
  set_name1 = "Node Set 1"
  set_num_node1 = 2
  nd_dist1 = 0
  ExodusII.ex_put_node_set_param(f,set_id1,set_num_node1,nd_dist1)
  node_set1 = [5,9]
  ExodusII.ex_put_node_set(f,set_id1,node_set1)
  ExodusII.ex_put_name(f,ExodusII.EX_NODE_SET,set_id1,set_name1)

  set_id2 = 914
  set_name2 = "Node Set 2"
  set_num_node2 = 1
  nd_dist2 = 0
  ExodusII.ex_put_node_set_param(f,set_id2,set_num_node2)
  node_set2 = [10]
  ExodusII.ex_put_node_set(f,set_id2,node_set2)
  ExodusII.ex_put_name(f,ExodusII.EX_NODE_SET,set_id2,set_name2)

  # Set variables
  nvar_names = ["Nvar1","Nvar2","Nvar3"]
  ExodusII.ex_put_num_node_vars(f,length(nvar_names))
  ExodusII.ex_put_nodal_var_names(f,nvar_names)

  node_valsv3t1 = [1.7,2.3,3.9,4.5,5.7,6.8]
  ExodusII.ex_put_node_var(f,1,3,node_valsv3t1)

  evar_names = ["Evar1"]
  ExodusII.ex_put_num_elem_vars(f,length(evar_names))
  ExodusII.ex_put_elem_var_names(f,evar_names)

  elem_tab = [1,1]'
  ExodusII.ex_put_elem_var_tab(f,elem_tab)
  elem_valsv1b1t2 = [3.14,9.82,4.45,7.23]
  ExodusII.ex_put_elem_var(f,2,1,block_id1,elem_valsv1b1t2)
  elem_valsv1b2t3 = [1.0,2.0,3.0]
  ExodusII.ex_put_elem_var(f,3,1,block_id2,elem_valsv1b2t3)

  # Set times
  times = [1.0,5.0,6.0,8.0]
  for (i,t) in enumerate(times)
    ExodusII.ex_put_time(f,i,t)
  end

  # Close the file
  ExodusII.ex_close(f)

  # Read the file and check it all matches
  f = ExodusII.ex_open(filename)
  @test f.title == dbtitle
  @test f.num_dim == num_dim
  @test f.num_nodes == num_nodes
  @test f.num_elem == num_elem
  @test f.num_elem_blk == num_elem_blk
  @test f.num_node_sets == num_node_sets
  @test f.num_side_sets == num_side_sets

  # Check Coords and int maps
  read_coord = ExodusII.ex_get_coord(f)
  @test read_coord == coord
  read_coord_names = ExodusII.ex_get_coord_names(f)
  @test read_coord_names == coord_names

  read_node_map = ExodusII.ex_get_node_num_map(f)
  @test read_node_map == node_map
  read_elem_map = ExodusII.ex_get_elem_num_map(f)
  @test read_elem_map == elem_map

  # Check elem blocks
  read_elem_blocks_ids = ExodusII.ex_get_elem_block_ids(f)
  @test read_elem_blocks_ids == [block_id1,block_id2]

  read_type,read_num_elem,read_nodes_per_elem,read_attr = ExodusII.ex_get_elem_block(f,block_id1)
  @test read_type == elem_type1
  @test read_num_elem == blk_num_elem1
  @test read_nodes_per_elem == nodes_per_elem1
  @test read_attr == el_attr1
  read_conn = ExodusII.ex_get_elem_connections(f,block_id1)
  @test read_conn == conn1

  read_type,read_num_elem,read_nodes_per_elem,read_attr = ExodusII.ex_get_elem_block(f,block_id2)
  @test read_type == elem_type2
  @test read_num_elem == blk_num_elem2
  @test read_nodes_per_elem == nodes_per_elem2
  @test read_attr == el_attr2
  read_conn = ExodusII.ex_get_elem_connections(f,block_id2)
  @test read_conn == conn2

  # Check node sets
  read_node_set_ids = ExodusII.ex_get_node_set_ids(f)
  @test read_node_set_ids == [set_id1,set_id2]

  read_num_node,read_num_dist = ExodusII.ex_get_node_set_param(f,set_id1)
  @test read_num_node == set_num_node1
  @test read_num_dist == nd_dist1
  read_ns = ExodusII.ex_get_node_set(f,set_id1)
  @test read_ns == node_set1
  read_ns_name = ExodusII.ex_get_name(f,ExodusII.EX_NODE_SET,set_id1)
  @test read_ns_name == set_name1

  read_num_node,read_num_dist = ExodusII.ex_get_node_set_param(f,set_id2)
  @test read_num_node == set_num_node2
  @test read_num_dist == nd_dist2
  read_ns = ExodusII.ex_get_node_set(f,set_id2)
  @test read_ns == node_set2
  read_ns_name = ExodusII.ex_get_name(f,ExodusII.EX_NODE_SET,set_id2)
  @test read_ns_name == set_name2

  # Check variables
  read_nvar_names = ExodusII.ex_get_nodal_var_names(f)
  @test read_nvar_names == nvar_names

  read_vals = ExodusII.ex_get_node_var_vals(f,nvar_names[3],1)
  @test read_vals == node_valsv3t1

  read_evar_names = ExodusII.ex_get_elem_var_names(f)
  @test read_evar_names == evar_names

  read_elem_tab = ExodusII.ex_get_elem_var_tab(f)
  @test read_elem_tab == elem_tab

  read_var = ExodusII.ex_get_elem_var_vals(f,2,1,block_id1)
  @test read_var == elem_valsv1b1t2
  read_var = ExodusII.ex_get_elem_var_vals(f,3,1,block_id2)
  @test read_var == elem_valsv1b2t3

  # Check times
  read_times = ExodusII.ex_get_all_times(f)
  @test read_times == times
  read_t = ExodusII.ex_get_time(f,3)
  @test read_t == times[3]

  ExodusII.ex_close(f)
end


text:
function binarysearch(lst::Vector{T}, value::T, low=1, high=length(lst)) where T
    if isempty(lst) return 0 end
    if low ≥ high
        if low > high || lst[low] != value
            return 0
        else
            return low
        end
    end
    mid = (low + high) ÷ 2
    if lst[mid] > value
        return binarysearch(lst, value, low, mid-1)
    elseif lst[mid] < value
        return binarysearch(lst, value, mid+1, high)
    else
        return mid
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@inline function block_mat_mat_mul!(::Multithreaded, C, A, B, sz)
    @inbounds @views begin
        C11 = C[1:sz,     1:sz]; C12 = C[1:sz,     sz+1:end]
        C21 = C[sz+1:end, 1:sz]; C22 = C[sz+1:end, sz+1:end]

        A11 = A[1:sz,     1:sz]; A12 = A[1:sz,     sz+1:end]
        A21 = A[sz+1:end, 1:sz]; A22 = A[sz+1:end, sz+1:end]

        B11 = B[1:sz,     1:sz]; B12 = B[1:sz,     sz+1:end]
        B21 = B[sz+1:end, 1:sz]; B22 = B[sz+1:end, sz+1:end]
    end
    @sync begin
        Threads.@spawn begin
            #_mul!(multithreaded,     C11, A11, B11, sz)
            gemm_kernel!(C11, A11, B11)
            _mul_add!(multithreaded, C11, A12, B21, sz)
        end
        Threads.@spawn begin
            _mul!(multithreaded,     C12, A11, B12, sz)
            _mul_add!(multithreaded, C12, A12, B22, sz)
        end
        Threads.@spawn begin
            _mul!(multithreaded,     C21, A21, B11, sz)
            _mul_add!(multithreaded, C21, A22, B21, sz)
        end
        _mul!(multithreaded,     C22, A21, B12, sz)
        _mul_add!(multithreaded, C22, A22, B22, sz)
    end
end

@inline function block_mat_mat_mul!(::Singlethreaded, C, A, B, sz)
    @inbounds @views begin
        C11 = C[1:sz,     1:sz]; C12 = C[1:sz,     sz+1:end]
        C21 = C[sz+1:end, 1:sz]; C22 = C[sz+1:end, sz+1:end]

        A11 = A[1:sz,     1:sz]; A12 = A[1:sz,     sz+1:end]
        A21 = A[sz+1:end, 1:sz]; A22 = A[sz+1:end, sz+1:end]

        B11 = B[1:sz,     1:sz]; B12 = B[1:sz,     sz+1:end]
        B21 = B[sz+1:end, 1:sz]; B22 = B[sz+1:end, sz+1:end]
    end
    #_mul!(singlethreaded,     C11, A11, B11, sz)
    gemm_kernel!(C11, A11, B11)
    _mul_add!(singlethreaded, C11, A12, B21, sz)
    _mul!(singlethreaded,     C12, A11, B12, sz)
    _mul_add!(singlethreaded, C12, A12, B22, sz)
    _mul!(singlethreaded,     C21, A21, B11, sz)
    _mul_add!(singlethreaded, C21, A22, B21, sz)
    _mul!(singlethreaded,     C22, A21, B12, sz)
    _mul_add!(singlethreaded, C22, A22, B22, sz)
end

function block_mat_vec_mul!(::Multithreaded, C, A, B, sz)
    @inbounds @views begin
        C11 = C[1:sz,     1:end];
        C21 = C[sz+1:end, 1:end];

        A11 = A[1:sz,     1:sz]; A12 = A[1:sz,     sz+1:end]
        A21 = A[sz+1:end, 1:sz]; A22 = A[sz+1:end, sz+1:end]

        B11 = B[1:sz,     1:end];
        B21 = B[sz+1:end, 1:end];
    end
    @sync begin
        Threads.@spawn begin
            #_mul!(multithreaded,     C11, A11, B11, sz)
            gemm_kernel!(C11, A11, B11)
            _mul_add!(multithreaded, C11, A12, B21, sz)
        end
        _mul!(multithreaded,     C21, A21, B11, sz)
        _mul_add!(multithreaded, C21, A22, B21, sz)
    end
end

function block_mat_vec_mul!(::Singlethreaded, C, A, B, sz)
    @inbounds @views begin
        C11 = C[1:sz,     1:end];
        C21 = C[sz+1:end, 1:end];

        A11 = A[1:sz,     1:sz]; A12 = A[1:sz,     sz+1:end]
        A21 = A[sz+1:end, 1:sz]; A22 = A[sz+1:end, sz+1:end]

        B11 = B[1:sz,     1:end];
        B21 = B[sz+1:end, 1:end];
    end
    #_mul!(singlethreaded,     C11, A11, B11, sz)
    gemm_kernel!(C11, A11, B11)
    _mul_add!(singlethreaded, C11, A12, B21, sz)
    _mul!(singlethreaded,     C21, A21, B11, sz)
    _mul_add!(singlethreaded, C21, A22, B21, sz)
end

function block_mat_vec_mul!(::Multithreaded, C::VecTypes, A, B::VecTypes, sz)
    @inbounds @views begin
        C11 = C[1:sz    ];
        C21 = C[sz+1:end];

        A11 = A[1:sz,     1:sz]; A12 = A[1:sz,     sz+1:end]
        A21 = A[sz+1:end, 1:sz]; A22 = A[sz+1:end, sz+1:end]

        B11 = B[1:sz    ];
        B21 = B[sz+1:end];
    end
    @sync begin
        Threads.@spawn begin
            gemm_kernel!(C11, A11, B11)
            _mul_add!(multithreaded, C11, A12, B21, sz)
        end
        _mul!(multithreaded,     C21, A21, B11, sz)
        _mul_add!(multithreaded, C21, A22, B21, sz)
    end
end

function block_mat_vec_mul!(::Singlethreaded, C::VecTypes, A, B::VecTypes, sz)
    @inbounds @views begin
        C11 = C[1:sz    ];
        C21 = C[sz+1:end];

        A11 = A[1:sz,     1:sz]; A12 = A[1:sz,     sz+1:end]
        A21 = A[sz+1:end, 1:sz]; A22 = A[sz+1:end, sz+1:end]

        B11 = B[1:sz    ];
        B21 = B[sz+1:end];
    end
    gemm_kernel!(C11, A11, B11)
    _mul_add!(singlethreaded, C11, A12, B21, sz)
    _mul!(singlethreaded,     C21, A21, B11, sz)
    _mul_add!(singlethreaded, C21, A22, B21, sz)
end

function block_covec_mat_mul!(::Multithreaded, C, A, B, sz)
    @inbounds @views begin
        C11 = C[1:end,    1:sz]; C12 = C[1:end,    sz+1:end]

        A11 = A[1:end,    1:sz]; A12 = A[1:end,    sz+1:end]

        B11 = B[1:sz,     1:sz]; B12 = B[1:sz,     sz+1:end]
        B21 = B[sz+1:end, 1:sz]; B22 = B[sz+1:end, sz+1:end]
    end
    @sync begin
        Threads.@spawn begin
            #_mul!(multithreaded,     C11, A11, B11, sz)
            gemm_kernel!(C11, A11, B11)
            _mul_add!(multithreaded, C11, A12, B21, sz)
        end
        _mul!(multithreaded,     C12, A11, B12, sz)
        _mul_add!(multithreaded, C12, A12, B22, sz)
    end
end

function block_covec_mat_mul!(::Singlethreaded, C, A, B, sz)
    @inbounds @views begin
        C11 = C[1:end,    1:sz]; C12 = C[1:end,    sz+1:end]

        A11 = A[1:end,    1:sz]; A12 = A[1:end,    sz+1:end]

        B11 = B[1:sz,     1:sz]; B12 = B[1:sz,     sz+1:end]
        B21 = B[sz+1:end, 1:sz]; B22 = B[sz+1:end, sz+1:end]
    end
    #_mul!(threading,     C11, A11, B11, sz)
    gemm_kernel!(C11, A11, B11)
    _mul_add!(multithreaded, C11, A12, B21, sz)
    _mul!(multithreaded,     C12, A11, B12, sz)
    _mul_add!(multithreaded, C12, A12, B22, sz)
end

function block_vec_covec_mul!(::Multithreaded, C, A, B, sz)
    @inbounds @views begin
        C11 = C[1:sz,     1:sz]; C12 = C[1:sz,     sz+1:end]
        C21 = C[sz+1:end, 1:sz]; C22 = C[sz+1:end, sz+1:end]

        A11 = A[1:sz,     1:end];
        A21 = A[sz+1:end, 1:end];

        B11 = B[1:end,     1:sz]; B12 = B[1:end,     sz+1:end]
    end
    @sync begin
        Threads.@spawn begin
            #_mul!(multithreaded,     C11, A11, B11, sz)
            gemm_kernel!(C11, A11, B11)
        end
        Threads.@spawn begin
            _mul!(multithreaded, C12, A11, B12, sz)
        end
        Threads.@spawn begin
            _mul!(multithreaded, C21, A21, B11, sz)
        end
        _mul!(multithreaded, C22, A21, B12, sz)
    end
end

function block_vec_covec_mul!(::Singlethreaded, C, A, B, sz)
    @inbounds @views begin
        C11 = C[1:sz,     1:sz]; C12 = C[1:sz,     sz+1:end]
        C21 = C[sz+1:end, 1:sz]; C22 = C[sz+1:end, sz+1:end]

        A11 = A[1:sz,     1:end];
        A21 = A[sz+1:end, 1:end];

        B11 = B[1:end,     1:sz]; B12 = B[1:end,     sz+1:end]
    end
    #_mul!(singlethreaded,     C11, A11, B11, sz)
    gemm_kernel!(C11, A11, B11)
    _mul!(singlethreaded, C12, A11, B12, sz)
    _mul!(singlethreaded, C21, A21, B11, sz)
    _mul!(singlethreaded, C22, A21, B12, sz)
end

function block_covec_vec_mul!(threading::Threading, C, A, B, sz)
    @inbounds @views begin
        A11 = A[1:end,    1:sz]; A12 = A[1:end,     sz+1:end]

        B11 = B[1:sz,     1:end];
        B21 = B[sz+1:end, 1:end];
    end
    gemm_kernel!(C, A11, B11)
    #_mul!(threading,     C, A11, B11, sz)
    _mul_add!(threading, C, A12, B21, sz)
end

function block_covec_vec_mul!(threading::Threading, C::VecTypes, A, B::VecTypes, sz)
    @inbounds @views begin
        A11 = A[1:end,    1:sz]; A12 = A[1:end,     sz+1:end]

        B11 = B[1:sz    ];
        B21 = B[sz+1:end];
    end
    gemm_kernel!(C, A11, B11)
    _mul_add!(threading, C, A12, B21, sz)
end

@inline function block_mat_mat_mul_add!(::Multithreaded, C, A, B, sz, ::Val{factor} = Val(1)) where {factor}
    @inbounds @views begin
        C11 = C[1:sz,     1:sz]; C12 = C[1:sz,     sz+1:end]
        C21 = C[sz+1:end, 1:sz]; C22 = C[sz+1:end, sz+1:end]

        A11 = A[1:sz,     1:sz]; A12 = A[1:sz,     sz+1:end]
        A21 = A[sz+1:end, 1:sz]; A22 = A[sz+1:end, sz+1:end]

        B11 = B[1:sz,     1:sz]; B12 = B[1:sz,     sz+1:end]
        B21 = B[sz+1:end, 1:sz]; B22 = B[sz+1:end, sz+1:end]
    end
    @sync begin
        Threads.@spawn begin
            add_gemm_kernel!(C11, A11, B11, Val(factor))
            _mul_add!(multithreaded, C11, A12, B21, sz, Val(factor))
        end
        Threads.@spawn begin
            _mul_add!(multithreaded, C12, A11, B12, sz, Val(factor))
            _mul_add!(multithreaded, C12, A12, B22, sz, Val(factor))
        end
        Threads.@spawn begin
            _mul_add!(multithreaded, C21, A21, B11, sz, Val(factor))
            _mul_add!(multithreaded, C21, A22, B21, sz, Val(factor))
        end
        _mul_add!(multithreaded, C22, A21, B12, sz, Val(factor))
        _mul_add!(multithreaded, C22, A22, B22, sz, Val(factor))
    end
end

@inline function block_mat_mat_mul_add!(::Singlethreaded, C, A, B, sz, ::Val{factor} = Val(1)) where {factor}
    @inbounds @views begin
        C11 = C[1:sz,     1:sz]; C12 = C[1:sz,     sz+1:end]
        C21 = C[sz+1:end, 1:sz]; C22 = C[sz+1:end, sz+1:end]

        A11 = A[1:sz,     1:sz]; A12 = A[1:sz,     sz+1:end]
        A21 = A[sz+1:end, 1:sz]; A22 = A[sz+1:end, sz+1:end]

        B11 = B[1:sz,     1:sz]; B12 = B[1:sz,     sz+1:end]
        B21 = B[sz+1:end, 1:sz]; B22 = B[sz+1:end, sz+1:end]
    end
    add_gemm_kernel!(C11, A11, B11, Val(factor))
    _mul_add!(singlethreaded, C11, A12, B21, sz, Val(factor))
    _mul_add!(singlethreaded, C12, A11, B12, sz, Val(factor))
    _mul_add!(singlethreaded, C12, A12, B22, sz, Val(factor))
    _mul_add!(singlethreaded, C21, A21, B11, sz, Val(factor))
    _mul_add!(singlethreaded, C21, A22, B21, sz, Val(factor))
    _mul_add!(singlethreaded, C22, A21, B12, sz, Val(factor))
    _mul_add!(singlethreaded, C22, A22, B22, sz, Val(factor))
end

function block_mat_vec_mul_add!(::Multithreaded, C, A, B, sz, ::Val{factor} = Val(1)) where {factor}
    @inbounds @views begin
        C11 = C[1:sz,     1:end];
        C21 = C[sz+1:end, 1:end];

        A11 = A[1:sz,     1:sz]; A12 = A[1:sz,     sz+1:end]
        A21 = A[sz+1:end, 1:sz]; A22 = A[sz+1:end, sz+1:end]

        B11 = B[1:sz,     1:end];
        B21 = B[sz+1:end, 1:end];
    end
    @sync begin
        Threads.@spawn begin
            add_gemm_kernel!(C11, A11, B11, Val(factor))
            _mul_add!(multithreaded, C11, A12, B21, sz, Val(factor))
        end
        _mul_add!(multithreaded, C21, A21, B11, sz, Val(factor))
        _mul_add!(multithreaded, C21, A22, B21, sz, Val(factor))
    end
end

function block_mat_vec_mul_add!(::Singlethreaded, C, A, B, sz, ::Val{factor} = Val(1)) where {factor}
    @inbounds @views begin
        C11 = C[1:sz,     1:end];
        C21 = C[sz+1:end, 1:end];

        A11 = A[1:sz,     1:sz]; A12 = A[1:sz,     sz+1:end]
        A21 = A[sz+1:end, 1:sz]; A22 = A[sz+1:end, sz+1:end]

        B11 = B[1:sz,     1:end];
        B21 = B[sz+1:end, 1:end];
    end
    add_gemm_kernel!(C11, A11, B11, Val(factor))
    _mul_add!(singlethreaded, C11, A12, B21, sz, Val(factor))
    _mul_add!(singlethreaded, C21, A21, B11, sz, Val(factor))
    _mul_add!(singlethreaded, C21, A22, B21, sz, Val(factor))
end

function block_mat_vec_mul_add!(::Multithreaded, C::VecTypes, A, B::VecTypes, sz, ::Val{factor} = Val(1)) where {factor}
    @inbounds @views begin
        C11 = C[1:sz    ];
        C21 = C[sz+1:end];

        A11 = A[1:sz,     1:sz]; A12 = A[1:sz,     sz+1:end]
        A21 = A[sz+1:end, 1:sz]; A22 = A[sz+1:end, sz+1:end]

        B11 = B[1:sz    ];
        B21 = B[sz+1:end];
    end
    @sync begin
        Threads.@spawn begin
            add_gemm_kernel!(C11, A11, B11, Val(factor))
            _mul_add!(multithreaded, C11, A12, B21, sz, Val(factor))
        end
        _mul_add!(multithreaded, C21, A21, B11, sz, Val(factor))
        _mul_add!(multithreaded, C21, A22, B21, sz, Val(factor))
    end
end

function block_mat_vec_mul_add!(::Singlethreaded, C::VecTypes, A, B::VecTypes, sz, ::Val{factor} = Val(1)) where {factor}
    @inbounds @views begin
        C11 = C[1:sz    ];
        C21 = C[sz+1:end];

        A11 = A[1:sz,     1:sz]; A12 = A[1:sz,     sz+1:end]
        A21 = A[sz+1:end, 1:sz]; A22 = A[sz+1:end, sz+1:end]

        B11 = B[1:sz    ];
        B21 = B[sz+1:end];
    end
    add_gemm_kernel!(C11, A11, B11, Val(factor))
    _mul_add!(multithreaded, C11, A12, B21, sz, Val(factor))
    _mul_add!(multithreaded, C21, A21, B11, sz, Val(factor))
    _mul_add!(multithreaded, C21, A22, B21, sz, Val(factor))
end

function block_covec_mat_mul_add!(::Multithreaded, C, A, B, sz, ::Val{factor} = Val(1)) where {factor}
    @inbounds @views begin
        C11 = C[1:end,    1:sz]; C12 = C[1:end,     sz+1:end]

        A11 = A[1:end,    1:sz]; A12 = A[1:end,    sz+1:end]

        B11 = B[1:sz,     1:sz]; B12 = B[1:sz,     sz+1:end]
        B21 = B[sz+1:end, 1:sz]; B22 = B[sz+1:end, sz+1:end]
    end
    @sync begin
        Threads.@spawn begin
            add_gemm_kernel!(C11, A11, B11, Val(factor))
            _mul_add!(multithreaded, C11, A12, B21, sz, Val(factor))
        end
        _mul_add!(multithreaded, C12, A11, B12, sz, Val(factor))
        _mul_add!(multithreaded, C12, A12, B22, sz, Val(factor))
    end
end

function block_covec_mat_mul_add!(::Singlethreaded, C, A, B, sz, ::Val{factor} = Val(1)) where {factor}
    @inbounds @views begin
        C11 = C[1:end,    1:sz]; C12 = C[1:end,     sz+1:end]

        A11 = A[1:end,    1:sz]; A12 = A[1:end,    sz+1:end]

        B11 = B[1:sz,     1:sz]; B12 = B[1:sz,     sz+1:end]
        B21 = B[sz+1:end, 1:sz]; B22 = B[sz+1:end, sz+1:end]
    end
    add_gemm_kernel!(C11, A11, B11, Val(factor))
    _mul_add!(singlethreaded, C11, A12, B21, sz, Val(factor))
    _mul_add!(singlethreaded, C12, A11, B12, sz, Val(factor))
    _mul_add!(singlethreaded, C12, A12, B22, sz, Val(factor))
end

function block_vec_covec_mul_add!(::Multithreaded, C, A, B, sz, ::Val{factor} = Val(1)) where {factor}
    @inbounds @views begin
        C11 = C[1:sz,     1:sz]; C12 = C[1:sz,     sz+1:end]
        C21 = C[sz+1:end, 1:sz]; C22 = C[sz+1:end, sz+1:end]

        A11 = A[1:sz,     1:end];
        A21 = A[sz+1:end, 1:end];

        B11 = B[1:end,    1:sz]; B12 = B[1:end,    sz+1:end]
    end
    @sync begin
        Threads.@spawn begin
            add_gemm_kernel!(C11, A11, B11, Val(factor))
        end
        Threads.@spawn begin
            _mul_add!(multithreaded, C12, A11, B12, sz, Val(factor))
        end
        Threads.@spawn begin
            _mul_add!(multithreaded, C21, A21, B11, sz, Val(factor))
        end
        _mul_add!(multithreaded, C22, A21, B12, sz, Val(factor))
    end
end

function block_vec_covec_mul_add!(::Singlethreaded, C, A, B, sz, ::Val{factor} = Val(1)) where {factor}
    @inbounds @views begin
        C11 = C[1:sz,     1:sz]; C12 = C[1:sz,     sz+1:end]
        C21 = C[sz+1:end, 1:sz]; C22 = C[sz+1:end, sz+1:end]

        A11 = A[1:sz,     1:end];
        A21 = A[sz+1:end, 1:end];

        B11 = B[1:end,    1:sz]; B12 = B[1:end,    sz+1:end]
    end
    add_gemm_kernel!(C11, A11, B11, Val(factor))
    _mul_add!(singlethreaded, C12, A11, B12, sz, Val(factor))
    _mul_add!(singlethreaded, C21, A21, B11, sz, Val(factor))
    _mul_add!(singlethreaded, C22, A21, B12, sz, Val(factor))
end

function block_covec_vec_mul_add!(threading::Threading, C, A, B, sz, ::Val{factor} = Val(1)) where {factor}
    @inbounds @views begin
        A11 = A[1:end,    1:sz]; A12 = A[1:end,    sz+1:end]

        B11 = B[1:sz,     1:end];
        B21 = B[sz+1:end, 1:end];
    end
    add_gemm_kernel!(C, A11, B11, Val(factor))
    _mul_add!(threading, C, A12, B21, sz, Val(factor))
end

function block_covec_vec_mul_add!(threading::Threading, C::VecTypes, A, B::VecTypes, sz, ::Val{factor} = Val(1)) where {factor}
    @inbounds @views begin
        A11 = A[1:end,    1:sz]; A12 = A[1:end,    sz+1:end]

        B11 = B[1:sz    ];
        B21 = B[sz+1:end];
    end
    add_gemm_kernel!(C, A11, B11, Val(factor))
    _mul_add!(threading, C, A12, B21, sz, Val(factor))
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct Dimension <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    displayName::Any # spec type: Union{ Nothing, String } # spec name: displayName

    function Dimension(;name=nothing, displayName=nothing)
        o = new()
        validate_property(Dimension, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(Dimension, Symbol("displayName"), displayName)
        setfield!(o, Symbol("displayName"), displayName)
        o
    end
end # type Dimension

const _property_map_Dimension = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("displayName")=>Symbol("displayName"))
const _property_types_Dimension = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("displayName")=>"String")
Base.propertynames(::Type{ Dimension }) = collect(keys(_property_map_Dimension))
Swagger.property_type(::Type{ Dimension }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_Dimension[name]))}
Swagger.field_name(::Type{ Dimension }, property_name::Symbol) =  _property_map_Dimension[property_name]

function check_required(o::Dimension)
    true
end

function validate_property(::Type{ Dimension }, name::Symbol, val)
end


text:
using HyperbolicPlane, Plots

function circle_example()
    P = HPoint(1,pi/6)
    C = HCircle(P,1.5)
    X,Y,Z = points_on_circle(C)
    T = X+Y+Z
    set_color(T,:red)
    set_radius(P,2)
    L1 = bisector(X,Y)
    L2 = bisector(X,Z)
    L3 = bisector(Y,Z)
    set_color(L1,:green)
    set_color(L2,:green)
    set_color(L3,:green)
    plot()
    draw(C,T,L1,L2,L3,P,HPlane())
    finish()
end

function rand_circle_example()
    X,Y,Z = [ RandomHPoint() for j=1:3 ]
    C = HCircle()
    try
        C = HCircle(X,Y,Z)
    catch
        println("Sorry: The three points we generated are not contained in a common circle")
        println("Please try again")
        return
    end
    set_thickness(C,2)
    P = get_center(C)

    T = X+Y+Z
    set_color(T,:red)
    set_radius(P,2)
    L1 = bisector(X,Y)
    L2 = bisector(X,Z)
    L3 = bisector(Y,Z)
    set_color(L1,:green)
    set_color(L2,:green)
    set_color(L3,:green)

    plot()
    draw(C,T,L1,L2,L3,P,HPlane())
    finish()
end


text:
#=
mri-sense.jl
Methods for generating phantom spectra based on the 2012 paper
[Guerquin-Kern et al.](http://doi.org/10.1109/TMI.2011.2174158)
that combines analytical k-space values of the phantom
with an analytical model for MRI sensitivity maps.
=#

using LazyGrids: ndgrid
using ImagePhantoms: spectrum, Object2d

export mri_smap_basis, mri_smap_fit, spectrum, mri_spectra

"""
   mri_smap_basis(mask ; kmax, kt, ki)

Construct Fourier basis for representing MRI sensitivity maps
in terms of separable complex exponential signals, products of
of the form `basis = (k,N) -> exp.(2im * π * nfun(N) * kfun(k,N))`.
The default is `nfun(N) = -(N÷2):(N÷2)-1` which is suitable for even `N` only.
The default is `kfun(k,N) = k / 2N`, which is DCT-II like frequencies,
leading to better boundary behavior than the DFT frequencies `k/N`.

# Input
* `mask::AbstractArray{Bool,D}` binary support mask for region to reconstruct

# Option
* `kmax::Int = 5` default frequency index -kmax:kmax in all dimensions
* `kfun::Function = (k,N) -> k / (2N)` # DCT-II frequency
* `deltas::NTuple{D,<:Number} = ones(D)` pixel sizes
(For additional options `kmaxs`, `kt`, `ki`, see code.)

# Output
* `(; B, ν)` where `B` is basis matrix of size `count(mask) × nk`
  where typically `nk = (2*kmax+1)^D`
  and `ν` is `nk` frequency tuples;
  each tuple has form `ν = kfun.(Tuple(k), size(mask)) ./ deltas`.
"""
function mri_smap_basis(
    mask::AbstractArray{Bool,D} ;
    kmax::Int = 5,
    kmaxs::NTuple{D, Int} = ntuple(i -> kmax, D),
    kt::NTuple{D, UnitRange{<:Integer}} = ntuple(i -> -kmaxs[i]:kmaxs[i], D),
    ki::CartesianIndices{D} = CartesianIndices(kt),
    deltas::NTuple{D,<:Number} = ones(D),
    kfun::Function = (k,N) -> k / (2N), # DCT-II frequency
    nfun::Function = (N) -> -(N÷2):(N÷2)-1,
    basis::Function = (k,N) -> exp.(2im * π * nfun(N) * kfun(k,N)), # 1D basis
    T::DataType = ComplexF32,
) where D
    B = zeros(T, count(mask), length(ki))
    for (i, k) in enumerate(ki)
        tmp = ndgrid([basis(k[d], size(mask,d)) for d = 1:D]...)
        tmp = .*(tmp...)
        B[:,i] = tmp[mask]
    end

    ν = [kfun.(Tuple(k), size(mask)) ./ deltas for k in ki]
    return (; B, ν)
end

"""
    mri_smap_fit(smaps, embed ; mask, kwargs...)

Fit MRI sensitivity maps `smaps` using `mri_smap_basis(mask ; kwargs...)`.
Caller provides `ImageGeoms.embed` or equivalent.

Return named tuple `(B, ν, coefs, nrmse, smaps)`:
* `(B, ν)` from `mri_smap_basis`
* `coefs::Vector` : `[ncoil]` each of length `nk`
* `nrmse::Real` : `smaps` vs `smaps_fit`
* `smaps::Vector{Array{D}}` : `smaps_fit`
"""
function mri_smap_fit(
    smaps::Vector{<:AbstractArray{<:Number,D}},
    embed::Function ;
    mask::AbstractArray{Bool,D} = trues(size(smaps[1])),
    kwargs...
) where D
    (B, ν) = mri_smap_basis(mask ; kwargs...)
    coefs = [B \ smap[mask] for smap in smaps] # coefficients for each smap

    smaps_fit = [embed(B * coef, mask) for coef in coefs]
    ss(v) = sum(x -> sum(abs2, x), v)
    nrmse = sqrt( ss(smaps_fit - smaps) / ss(smaps) )
    return (; B, ν, coefs, nrmse, smaps = smaps_fit)
end

"""
    spectrum(ob::Object2d, coefs::AbstractVector, f::)
Version of `spectrum(ob)` suitable for parallel MRI with sensitivity maps
that were fit previously using `mri_smap_fit` for a single coil
with fit coefficients `coefs` and frequencies `f` (array of tuples).
"""
function spectrum(ob::Object2d, coefs::AbstractVector, f::Any)
    return (fx,fy) -> sum(
        i -> coefs[i] * spectrum(ob)(fx - f[i][1], fy - f[i][2]),
        1:length(f),
    )
end

"""
    spectrum(ob::Object2d, fit::NamedTuple, coil::Int)
Version of `spectrum(ob)` suitable for parallel MRI with sensitivity maps
that were fit previously using `mri_smap_fit`.
"""
spectrum(ob::Object2d, fit::NamedTuple, coil::Int) =
    spectrum(ob, fit.coefs[coil], fit.ν)

"""
    mri_spectra(fx, fy, oa::Array{<:Object2d}, fit::NamedTuple)
Version of `spectrum` suitable for parallel MRI with sensitivity maps
that were fit previously using `mri_smap_fit`.
Returns a Vector of `ncoil` kspace data Vectors of dimension `length(fx)`
"""
function mri_spectra(
    fx::AbstractVector,
    fy::AbstractVector,
    oa::Array{<:Object2d},
    fit::NamedTuple,
)
    ncoil = length(fit.smaps)
    return [sum(ob -> spectrum(ob, fit, ic).(fx,fy), oa) for ic=1:ncoil]
end


text:
# Convenience type to allow dispatch on solution objects that were created by Trixi
#
# This is a union of a Trixi-specific DiffEqBase.ODESolution and of Trixi's own
# TimeIntegratorSolution.
#
# Note: This is an experimental feature and may be changed in future releases without notice.
const TrixiODESolution = Union{ODESolution{T, N, uType, uType2, DType, tType, rateType, P} where
    {T, N, uType, uType2, DType, tType, rateType, P<:ODEProblem{uType_, tType_, isinplace, P_, F_} where
     {uType_, tType_, isinplace, P_<:AbstractSemidiscretization, F_<:ODEFunction{true, typeof(rhs!)}}}, TimeIntegratorSolution}

"""
    PlotData2D

Holds all relevant data for creating 2D plots of multiple solution variables and to visualize the
mesh.

!!! warning "Experimental implementation"
    This is an experimental feature and may change in future releases.
"""
struct PlotData2D{Coordinates, Data, VariableNames, Vertices}
  x::Coordinates
  y::Coordinates
  data::Data
  variable_names::VariableNames
  mesh_vertices_x::Vertices
  mesh_vertices_y::Vertices
end

"""
    PlotData2D(u, semi;
               solution_variables=cons2prim,
               grid_lines=true, max_supported_level=11, nvisnodes=nothing,
               slice_axis=:z, slice_axis_intercept=0)

Create a new `PlotData2D` object that can be used for visualizing 2D/3D DGSEM solution data array
`u` with `Plots.jl`. All relevant geometrical information is extracted from the semidiscretization
`semi`. By default, the conservative variables from the solution are used for plotting. This can be
changed by passing an appropriate conversion function to `solution_variables`.

If `grid_lines` is `true`, also extract grid vertices for visualizing the mesh. The output
resolution is indirectly set via `max_supported_level`: all data is interpolated to
`2^max_supported_level` uniformly distributed points in each spatial direction, also setting the
maximum allowed refinement level in the solution. `nvisnodes` specifies the number of visualization
nodes to be used. If it is `nothing`, twice the number of solution DG nodes are used for
visualization, and if set to `0`, exactly the number of nodes in the DG elements are used.

When visualizing data from a three-dimensional simulation, a 2D slice is extracted for plotting.
`slice_axis` specifies the axis orthogonal to the slice plane and may be `:x`, `:y`, or `:z`. The
point on the slice axis where it intersects with the slice plane is given in `slice_axis_intercept`.
Both of these values are ignored when visualizing 2D data.

!!! warning "Experimental implementation"
    This is an experimental feature and may change in future releases.

# Examples
```julia
julia> using Trixi, Plots

julia> trixi_include(default_example())
[...]

julia> pd = PlotData2D(sol)
PlotData2D(...)

julia> plot(pd) # To plot all available variables

julia> plot(pd["scalar"]) # To plot only a single variable

julia> plot!(getmesh(pd)) # To add grid lines to the plot
```
"""
function PlotData2D(u, semi;
                    solution_variables=cons2prim,
                    grid_lines=true, max_supported_level=11, nvisnodes=nothing,
                    slice_axis=:z, slice_axis_intercept=0)
  mesh, equations, solver, _ = mesh_equations_solver_cache(semi)
  @assert ndims(mesh) in (2, 3) "unsupported number of dimensions $ndims (must be 2 or 3)"

  # Extract mesh info
  center_level_0 = mesh.tree.center_level_0
  length_level_0 = mesh.tree.length_level_0
  leaf_cell_ids = leaf_cells(mesh.tree)
  coordinates = mesh.tree.coordinates[:, leaf_cell_ids]
  levels = mesh.tree.levels[leaf_cell_ids]

  unstructured_data = get_unstructured_data(u, semi, solution_variables)
  x, y, data, mesh_vertices_x, mesh_vertices_y = get_data_2d(center_level_0, length_level_0,
                                                             leaf_cell_ids, coordinates, levels,
                                                             ndims(mesh), unstructured_data,
                                                             nnodes(solver), grid_lines,
                                                             max_supported_level, nvisnodes,
                                                             slice_axis, slice_axis_intercept)
  variable_names = SVector(varnames(solution_variables, equations))

  return PlotData2D(x, y, data, variable_names, mesh_vertices_x, mesh_vertices_y)
end

"""
    PlotData2D(u_ode::AbstractVector, semi; kwargs...)

Create a `PlotData2D` object from a one-dimensional ODE solution `u_ode` and the semidiscretization
`semi`.

!!! warning "Experimental implementation"
    This is an experimental feature and may change in future releases.
"""
PlotData2D(u_ode::AbstractVector, semi; kwargs...) = PlotData2D(wrap_array(u_ode, semi), semi; kwargs...)

"""
    PlotData2D(sol::Union{DiffEqBase.ODESolution,TimeIntegratorSolution}; kwargs...)

Create a `PlotData2D` object from a solution object created by either `OrdinaryDiffEq.solve!` (which
returns a `DiffEqBase.ODESolution`) or Trixi's own `solve!` (which returns a
`TimeIntegratorSolution`).

!!! warning "Experimental implementation"
    This is an experimental feature and may change in future releases.
"""
PlotData2D(sol::TrixiODESolution; kwargs...) = PlotData2D(sol.u[end], sol.prob.p; kwargs...)

# Auxiliary data structure for visualizing a single variable
#
# Note: This is an experimental feature and may be changed in future releases without notice.
struct PlotDataSeries2D{PD<:PlotData2D}
  plot_data::PD
  variable_id::Int
end

# Show only a truncated output for convenience (the full data does not make sense)
function Base.show(io::IO, pds::PlotDataSeries2D)
  print(io, "PlotDataSeries2D{", typeof(pds.plot_data), "}(<plot_data::PlotData2D>, ",
        pds.variable_id, ")")
end

"""
    Base.getindex(pd::PlotData2D, variable_name)

Extract a single variable `variable_name` from `pd` for plotting with `Plots.plot`.

!!! warning "Experimental implementation"
    This is an experimental feature and may change in future releases.
"""
function Base.getindex(pd::PlotData2D, variable_name)
  variable_id = findfirst(isequal(variable_name), pd.variable_names)

  if isnothing(variable_id)
    throw(KeyError(variable_name))
  end

  return PlotDataSeries2D(pd, variable_id)
end

# Define additional methods for convenience
Base.firstindex(pd::PlotData2D) = first(pd.variable_names)
Base.lastindex(pd::PlotData2D) = last(pd.variable_names)
Base.length(pd::PlotData2D) = length(pd.variable_names)
Base.size(pd::PlotData2D) = (length(pd),)
Base.keys(pd::PlotData2D) = tuple(pd.variable_names...)
Base.eltype(pd::PlotData2D) = Pair{String, PlotDataSeries2D}
function Base.iterate(pd::PlotData2D, state=1)
  if state > length(pd)
    return nothing
  else
    return (pd.variable_names[state] => pd[pd.variable_names[state]], state + 1)
  end
end

# Show only a truncated output for convenience (the full data does not make sense)
function Base.show(io::IO, pd::PlotData2D)
  print(io, "PlotData2D{",
            typeof(pd.x), ",",
            typeof(pd.data), ",",
            typeof(pd.variable_names), ",",
            typeof(pd.mesh_vertices_x),
            "}(<x>, <y>, <data>, <variable_names>, <mesh_vertices_x>, <mesh_vertices_y>)")
end

# Auxiliary data structure for visualizing the mesh
#
# Note: This is an experimental feature and may be changed in future releases without notice.
struct PlotMesh2D{PD<:PlotData2D}
  plot_data::PD
end

# Show only a truncated output for convenience (the full data does not make sense)
function Base.show(io::IO, pm::PlotMesh2D)
  print(io, "PlotMesh2D{", typeof(pm.plot_data), "}(<plot_data::PlotData2D>)")
end

"""
    getmesh(pd::PlotData2D)

Extract grid lines from `pd` for plotting with `Plots.plot`.

!!! warning "Experimental implementation"
    This is an experimental feature and may change in future releases.
"""
getmesh(pd::PlotData2D) = PlotMesh2D(pd)

# Visualize a single variable in a 2D plot (default: heatmap)
#
# Note: This is an experimental feature and may be changed in future releases without notice.
@recipe function f(pds::PlotDataSeries2D)
  @unpack plot_data, variable_id = pds
  @unpack x, y, data, variable_names = plot_data

  # Set geometric properties
  xlims --> (x[begin], x[end])
  ylims --> (y[begin], y[end])
  aspect_ratio --> :equal

  # Set annotation properties
  legend -->  :none
  title --> variable_names[variable_id]
  colorbar --> :true

  # Set series properties
  seriestype --> :heatmap

  # Return data for plotting
  x, y, data[variable_id]
end

# Visualize the mesh in a 2D plot
#
# Note: This is an experimental feature and may be changed in future releases without notice.
@recipe function f(pm::PlotMesh2D)
  @unpack plot_data = pm
  @unpack x, y, mesh_vertices_x, mesh_vertices_y = plot_data

  # Set geometric and annotation properties
  xlims --> (x[begin], x[end])
  ylims --> (y[begin], y[end])
  aspect_ratio --> :equal
  legend -->  :none

  # Set series properties
  seriestype := :path
  linecolor := :black
  linewidth := 1

  # Return data for plotting
  mesh_vertices_x, mesh_vertices_y
end

# Plot all available variables at once for convenience
#
# Note: This is an experimental feature and may be changed in future releases without notice.
@recipe function f(pd::PlotData2D)
  # Create layout that is as square as possible, with a preference for more columns than rows if not
  cols = ceil(Int, sqrt(length(pd)))
  rows = ceil(Int, length(pd)/cols)
  layout := (rows, cols)

  # Plot all existing variables
  for (i, (variable_name, series)) in enumerate(pd)
    @series begin
      subplot := i
      series
    end
  end

  # Fill remaining subplots with empty plot
  for i in (length(pd)+1):(rows*cols)
    @series begin
      subplot := i
      axis := false
      ticks := false
      legend := false
      [], []
    end
  end
end

# Create a PlotData2D plot directly from a TrixiODESolution for convenience
#
# Note: This is an experimental feature and may be changed in future releases without notice.
#
# Note: If you change the defaults values here, you need to also change them in the PlotData2D
#       constructor.
@recipe function f(sol::TrixiODESolution;
                   solution_variables=cons2prim,
                   grid_lines=true, max_supported_level=11, nvisnodes=nothing, slice_axis=:z,
                   slice_axis_intercept=0)
  return PlotData2D(sol;
                    solution_variables=solution_variables,
                    grid_lines=grid_lines, max_supported_level=max_supported_level,
                    nvisnodes=nvisnodes, slice_axis=slice_axis,
                    slice_axis_intercept=slice_axis_intercept)
end


text:
"""
Pre-computes an inverse nonuniform fast Fourier transform of type `N`.

For best performance, choose the right number of threads by `FFTW.set_num_threads(4)`, for example.
"""
struct iNUFFTPlan{N,T,S,PT} <: Plan{T}
    pt::PT
    TP::Toeplitz{T}
    r::Vector{T}
    p::Vector{T}
    Ap::Vector{T}
    ϵ::S
end

"""
Pre-computes an inverse nonuniform fast Fourier transform of type I.
"""
function plan_inufft1(ω::AbstractVector{T}, ϵ::T) where T<:AbstractFloat
    N = length(ω)
    p = plan_nufft1(ω, ϵ)
    pt = plan_nufft2(ω/N, ϵ)
    c = p*ones(Complex{T}, N)
    r = conj(c)
    avg = (r[1]+c[1])/2
    r[1] = avg
    c[1] = avg
    TP = Toeplitz(c, r)
    r = zero(c)
    p = zero(c)
    Ap = zero(c)

    iNUFFTPlan{1, eltype(TP), typeof(ϵ), typeof(pt)}(pt, TP, r, p, Ap, ϵ)
end

"""
Pre-computes an inverse nonuniform fast Fourier transform of type II.
"""
function plan_inufft2(x::AbstractVector{T}, ϵ::T) where T<:AbstractFloat
    N = length(x)
    pt = plan_nufft1(N*x, ϵ)
    r = pt*ones(Complex{T}, N)
    c = conj(r)
    avg = (r[1]+c[1])/2
    r[1] = avg
    c[1] = avg
    TP = Toeplitz(c, r)
    r = zero(c)
    p = zero(c)
    Ap = zero(c)

    iNUFFTPlan{2, eltype(TP), typeof(ϵ), typeof(pt)}(pt, TP, r, p, Ap, ϵ)
end

function (*)(p::iNUFFTPlan{N,T}, x::AbstractVector{V}) where {N,T,V}
    mul!(zeros(promote_type(T,V), length(x)), p, x)
end

function mul!(c::AbstractVector{T}, P::iNUFFTPlan{1,T}, f::AbstractVector{T}) where T
    pt, TP, r, p, Ap, ϵ = P.pt, P.TP, P.r, P.p, P.Ap, P.ϵ
    cg_for_inufft(TP, c, f, r, p, Ap, 50, 100ϵ)
    conj!(mul!(c, pt, conj!(c)))
end

function mul!(c::AbstractVector{T}, P::iNUFFTPlan{2,T}, f::AbstractVector{T}) where T
    pt, TP, r, p, Ap, ϵ = P.pt, P.TP, P.r, P.p, P.Ap, P.ϵ
    cg_for_inufft(TP, c, conj!(pt*conj!(f)), r, p, Ap, 50, 100ϵ)
    conj!(f)
    c
end

"""
Computes an inverse nonuniform fast Fourier transform of type I.
"""
inufft1(c::AbstractVector, ω::AbstractVector{T}, ϵ::T) where {T<:AbstractFloat} = plan_inufft1(ω, ϵ)*c

"""
Computes an inverse nonuniform fast Fourier transform of type II.
"""
inufft2(c::AbstractVector, x::AbstractVector{T}, ϵ::T) where {T<:AbstractFloat} = plan_inufft2(x, ϵ)*c

function cg_for_inufft(A::ToeplitzMatrices.AbstractToeplitz{T}, x::AbstractVector{T}, b::AbstractVector{T}, r::AbstractVector{T}, p::AbstractVector{T}, Ap::AbstractVector{T}, max_it::Integer, tol::Real) where T
	n = length(b)
	n1, n2 = size(A)
	n == n1 == n2 || throw(DimensionMismatch(""))
    nrmb = norm(b)
    if nrmb == 0 nrmb = one(typeof(nrmb)) end
	copyto!(x, b)
    fill!(r, zero(T))
    fill!(p, zero(T))
    fill!(Ap, zero(T))
    # r = b - A*x
    copyto!(r, b)
    mul!(r, A, x, -one(T), one(T))
	copyto!(p, r)
	nrm2 = r⋅r
    for k = 1:max_it
        # Ap = A*p
        mul!(Ap, A, p)
		α = nrm2/(p⋅Ap)
        @inbounds @simd for l = 1:n
            x[l] += α*p[l]
            r[l] -= α*Ap[l]
        end
		nrm2new = r⋅r
        cst = nrm2new/nrm2
        @inbounds @simd for l = 1:n
            p[l] = muladd(cst, p[l], r[l])
        end
		nrm2 = nrm2new
        if sqrt(abs(nrm2)) ≤ tol*nrmb break end
	end
    return x
end


text:
export UniformMedium, Vacuum

"""
    UniformMedium(ϵ, μ)

Represents a medium with a constant permitivitty and permeability
Units are such that ϵ, μ are relative to the vacuum value(ϵ=μ=1)
"""
struct UniformMedium{ϵ, μ}
    function UniformMedium{ϵ, μ}() where {ϵ, μ}
        @assert ϵ isa Number
        @assert μ isa Number
        new{ϵ, μ}()
    end
end

const Vacuum = UniformMedium{1, 1}()

_get_ϵμ(::UniformMedium{ϵ, μ}) where {ϵ, μ} = (ϵ, μ)

"""
    _get_kz(k⃗::Tuple{Vararg{Number}}, k⃗²)
    _get_kz(k⃗::Tuple{Vararg{AbstractVector}}, k⃗²)
    _get_kz(k⃗::Tuple{Vararg{AbstractVector}}, ω, ϵ, μ)
    _get_kz(::PlaneWaves, ::UniformMedium)

Returns an array of kz from the dispersion relation k⃗² = k⃗⋅k⃗ = ω²ϵμ.
"""
_get_kz(pw::PlaneWaves, um::UniformMedium) = _get_kz(pw.k⃗, pw.ω, _get_ϵμ(um)...)
_get_kz(k⃗::Tuple{Vararg{AbstractVector}}, ω, ϵ, μ) = _get_kz(k⃗, ω^2 * ϵ * μ)
_get_kz(k⃗::Tuple{Vararg{AbstractVector}}, k⃗²) = _get_kz.(Iterators.product(k⃗...), k⃗²)
_get_kz(k⃗::Tuple{Vararg{Number}}, k⃗²) = sqrt(Complex(k⃗² - mapreduce(abs2, +, k⃗, init=zero(k⃗²))))


text:
###############################

# patch extraction

#############################################

function extract_patchs(dqx, dux, x, y, dx, dy, σmod, sizbl)
    lengbl = sizbl[1]*sizbl[2]
    gnside = size(dqx,1)
    x1rng  = repmat(collect(1:gnside).',gnside,1)
    x2rng  = x1rng.'

    lnx1 = x + dx
    lnx2 = y + dy

    dqxblks  = Vector{Float64}[]
    duxblks  = Vector{Float64}[]
    lnx1blks = Vector{Float64}[]
    lnx2blks = Vector{Float64}[]
    x1blks = Vector{Float64}[]
    x2blks = Vector{Float64}[]
    σmodblks = Vector{Float64}[]
    iblocks = Vector{Int64}[] # the linear indices of where the local blocks are located withing (x,y)
    for ib in 1:sizbl[1]:gnside, jb in 1:sizbl[2]:gnside
        iblock = find((ib .<= x1rng .< (ib+sizbl[1])) & (jb .<= x2rng .< (jb+sizbl[2])))
        if !isempty(iblock)
            push!(dqxblks,   dqx[iblock])
            push!(duxblks,   dux[iblock])
            push!(lnx1blks, lnx1[iblock])
            push!(lnx2blks, lnx2[iblock])
            push!(x1blks, x[iblock])
            push!(x2blks, y[iblock])
            push!(σmodblks, σmod[iblock])
            push!(iblocks, iblock)
        end
    end
    return lnx1blks, lnx2blks, dqxblks, duxblks, x1blks, x2blks, σmodblks, iblocks
end

function extract_patch_indices(size_vec, sizbl)
    lengbl = sizbl[1]*sizbl[2]
    gnside = size_vec
    x2blks = Vector{Float64}[]
    order = 1:size_vec
    if rem(size_vec+1,lengbl) == 0
        blk_boundaries = collect(1:lengbl:size_vec+1)
    else
        blk_boundaries = vcat(1:lengbl:size_vec+1, size_vec+1)
    end
    block_indices = Vector{Int64}[]
    for k = 2:length(blk_boundaries)
        iblock = order[blk_boundaries[k-1]:blk_boundaries[k]-1]
        if !isempty(iblock)
            push!(block_indices,  iblock)
        end
    end
    return block_indices
end

#######################################################

# XYcov constructs all the covariance matrices

#######################################################

function scalebesselj(n::Int, u::Array{Float64,1})
    rtn = zeros(u)
    @inbounds @simd for i in eachindex(rtn)
        rtn[i] = besselj(n, u[i]) / (u[i] ^ n)
    end
    rtn[u .== 0.0] = (0.5 ^ n) * (1 / gamma(n+1))
    return rtn
end

# Note: this computes Khat(u^2)
function vecϕKn{dm, T}(n, cXXk, ugrid, g::FFTgrid{dm, T})
    #ugrid   = linspace(0.0, √(g.period), 5000) .^ 3
    #ugrid   = vcat(0.0, ugrid)
    # pre-compute the radial power
    ddkbar  = g.deltk ^ dm / ( (2π)^dm )
    magk = unique(g.r)
    Θint = zeros(magk)
    @inbounds for inx in eachindex(magk, Θint)
        Θint[inx] = ddkbar * (-abs2(magk[inx])/2)^n * sum( cXXk[g.r .== magk[inx]] )
    end
    # now loop over ugrid
    vecϕKn  = zeros(ugrid)
    @inbounds @simd for inx in eachindex(ugrid)
        vecϕKn[inx] = dot(Θint, scalebesselj(n, magk.*ugrid[inx]))
    end
    return vecϕKn
end

"""
`YXcov(x1, x2, radKE, radKB)`. Compute the CE, SE, CB, SE cov matrices. Arguments `radKE, radKB`
should be computed by something like the following...
```
cψEk    = mCls.cEEk ./ (g.r.^4)
cψEk[1] = 0.0
vKE2, ugrid = PrimordialB.vecϕKn(2, cψEk, g)
vKE3, ugrid = PrimordialB.vecϕKn(3, cψEk, g)
vKE4, ugrid = PrimordialB.vecϕKn(4, cψEk, g)
vKE5, ugrid = PrimordialB.vecϕKn(5, cψEk, g)
vKE6, ugrid = PrimordialB.vecϕKn(6, cψEk, g)
radKE = hcat(ugrid, vKE2, vKE3, vKE4, vKE5, vKE6)

cψBk    = mCls.cBB0k ./ (g.r.^4)
cψBk[1] = 0.0
vKB2, ugrid = PrimordialB.vecϕKn(2, cψBk, g)
vKB3, ugrid = PrimordialB.vecϕKn(3, cψBk, g)
vKB4, ugrid = PrimordialB.vecϕKn(4, cψBk, g)
vKB5, ugrid = PrimordialB.vecϕKn(5, cψBk, g)
vKB6, ugrid = PrimordialB.vecϕKn(6, cψBk, g)
radKB = hcat(ugrid, vKB2, vKB3, vKB4, vKB5, vKB6)
```
"""
function YXcov(x1, x2, radKE, radKB)
    sqx1   = x1.^2
    sqx2   = x2.^2
    fx1    = x1.^4
    fx2    = x2.^4
    normx  = √(sqx1 + sqx2)
    normx2  = normx.^2
    dnormx2 = sqx1 - sqx2
    dnormx4  = dnormx2.^2

    sKE2xsq = Dierckx.Spline1D(radKE[:,1], radKE[:,2]; k=4, bc="extrapolate", s=0.0)
    sKE3xsq = Dierckx.Spline1D(radKE[:,1], radKE[:,3]; k=4, bc="extrapolate", s=0.0)
    sKE4xsq = Dierckx.Spline1D(radKE[:,1], radKE[:,4]; k=4, bc="extrapolate", s=0.0)
    sKE5xsq = Dierckx.Spline1D(radKE[:,1], radKE[:,5]; k=4, bc="extrapolate", s=0.0)
    sKE6xsq = Dierckx.Spline1D(radKE[:,1], radKE[:,6]; k=4, bc="extrapolate", s=0.0)
    KE2 = reshape(sKE2xsq(vec(normx)), size(normx))
    KE3 = reshape(sKE3xsq(vec(normx)), size(normx))
    KE4 = reshape(sKE4xsq(vec(normx)), size(normx))
    KE5 = reshape(sKE5xsq(vec(normx)), size(normx))
    KE6 = reshape(sKE6xsq(vec(normx)), size(normx))

    sKB2xsq = Dierckx.Spline1D(radKB[:,1], radKB[:,2]; k=4, bc="extrapolate", s=0.0)
    sKB3xsq = Dierckx.Spline1D(radKB[:,1], radKB[:,3]; k=4, bc="extrapolate", s=0.0)
    sKB4xsq = Dierckx.Spline1D(radKB[:,1], radKB[:,4]; k=4, bc="extrapolate", s=0.0)
    sKB5xsq = Dierckx.Spline1D(radKB[:,1], radKB[:,5]; k=4, bc="extrapolate", s=0.0)
    sKB6xsq = Dierckx.Spline1D(radKB[:,1], radKB[:,6]; k=4, bc="extrapolate", s=0.0)
    KB2 = reshape(sKB2xsq(vec(normx)), size(normx))
    KB3 = reshape(sKB3xsq(vec(normx)), size(normx))
    KB4 = reshape(sKB4xsq(vec(normx)), size(normx))
    KB5 = reshape(sKB5xsq(vec(normx)), size(normx))
    KB6 = reshape(sKB6xsq(vec(normx)), size(normx))

    rtnCE, rtnCB         = zeros(normx), zeros(normx)
    rtnSE, rtnSB         = zeros(normx), zeros(normx)
    rtnSCE, rtnSCB       = zeros(normx), zeros(normx)
    rtn∂11CE, rtn∂11CB   = zeros(normx), zeros(normx)
    rtn∂11SE, rtn∂11SB   = zeros(normx), zeros(normx)
    rtn∂11SCE, rtn∂11SCB = zeros(normx), zeros(normx)
    rtn∂22CE, rtn∂22CB   = zeros(normx), zeros(normx)
    rtn∂22SE, rtn∂22SB   = zeros(normx), zeros(normx)
    rtn∂22SCE, rtn∂22SCB = zeros(normx), zeros(normx)
    rtn∂12CE, rtn∂12CB   = zeros(normx), zeros(normx)
    rtn∂12SE, rtn∂12SB   = zeros(normx), zeros(normx)
    rtn∂12SCE, rtn∂12SCB = zeros(normx), zeros(normx)
    @inbounds @simd for inx in eachindex(rtnCE)
        # CE, CB
        t2, t3, t4 = 16, 32*normx2[inx], 16*dnormx4[inx]
        rtnCE[inx] = KE2[inx]*t2 + KE3[inx]*t3 + KE4[inx]*t4
        rtnCB[inx] = KB2[inx]*t2 + KB3[inx]*t3 + KB4[inx]*t4
        # SE, SB
        t2, t3, t4 = 16, 32*normx2[inx], 64*sqx1[inx]*sqx2[inx]
        rtnSE[inx] = KE2[inx]*t2 + KE3[inx]*t3 + KE4[inx]*t4
        rtnSB[inx] = KB2[inx]*t2 + KB3[inx]*t3 + KB4[inx]*t4
        # SCE, SCB
        t4 = 32 * dnormx2[inx] * x1[inx] * x2[inx]
        rtnSCE[inx] =  KE4[inx] * t4
        rtnSCB[inx] = - KB4[inx] * t4

        # ∂11CE, ∂11CB
        t3, t4, t5, t6 = -32*3, -32*18*sqx1[inx], -32*(13fx1[inx]-6sqx1[inx]*sqx2[inx]+fx2[inx]), - 32*(2sqx1[inx]^3-4fx1[inx]*sqx2[inx]+2sqx1[inx]*fx2[inx])
        rtn∂11CE[inx]  = KE3[inx]*t3 + KE4[inx]*t4 + KE5[inx]*t5 + KE6[inx]*t6
        rtn∂11CB[inx]  = KB3[inx]*t3 + KB4[inx]*t4 + KB5[inx]*t5 + KB6[inx]*t6
        # ∂11SE, ∂11SB
        t3, t4, t5, t6 = -32*3, -32*6*(2sqx1[inx]+sqx2[inx]), -32*4*sqx1[inx]*(sqx1[inx]+6sqx2[inx]), -32*8*fx1[inx]*sqx2[inx]
        rtn∂11SE[inx]  = KE3[inx]*t3 + KE4[inx]*t4 + KE5[inx]*t5 + KE6[inx]*t6
        rtn∂11SB[inx]  = KB3[inx]*t3 + KB4[inx]*t4 + KB5[inx]*t5 + KB6[inx]*t6
        # ∂11SCE, ∂11SCB
        t4, t5, t6      = -64*x1[inx]*x2[inx]*3, -64*x1[inx]*x2[inx]*(7sqx1[inx]-3sqx2[inx]), -64*x1[inx]*x2[inx]*2*sqx1[inx]*(sqx1[inx]-sqx2[inx])
        rtn∂11SCE[inx]  = KE4[inx]*t4 + KE5[inx]*t5 + KE6[inx]*t6
        rtn∂11SCB[inx]  = -KB4[inx]*t4 - KB5[inx]*t5 - KB6[inx]*t6

        # ∂22CE, ∂22CB
        t3, t4, t5, t6 = -32*3, -32*18*sqx2[inx], -32*(13fx2[inx]-6sqx2[inx]*sqx1[inx]+fx1[inx]), - 32*(2sqx2[inx]^3-4fx2[inx]*sqx1[inx]+2sqx2[inx]*fx1[inx])
        rtn∂22CE[inx]  = KE3[inx]*t3 + KE4[inx]*t4 + KE5[inx]*t5 + KE6[inx]*t6
        rtn∂22CB[inx]  = KB3[inx]*t3 + KB4[inx]*t4 + KB5[inx]*t5 + KB6[inx]*t6
        # ∂22SE, ∂22SB
        t3, t4, t5, t6 = -32*3, -32*6*(2sqx2[inx]+sqx1[inx]), -32*4*sqx2[inx]*(sqx2[inx]+6sqx1[inx]), -32*8*fx2[inx]*sqx1[inx]
        rtn∂22SE[inx]  = KE3[inx]*t3 + KE4[inx]*t4 + KE5[inx]*t5 + KE6[inx]*t6
        rtn∂22SB[inx]  = KB3[inx]*t3 + KB4[inx]*t4 + KB5[inx]*t5 + KB6[inx]*t6
        # ∂22SCE, ∂22SCB
        t4, t5, t6     = -64*x1[inx]*x2[inx]*(-3), -64*x1[inx]*x2[inx]*(-7sqx2[inx]+3sqx1[inx]), -64*x1[inx]*x2[inx]*2*sqx2[inx]*(sqx1[inx]-sqx2[inx])
        rtn∂22SCE[inx]  = KE4[inx]*t4 + KE5[inx]*t5 + KE6[inx]*t6
        rtn∂22SCB[inx]  = -KB4[inx]*t4 - KB5[inx]*t5 - KB6[inx]*t6

        # ∂12CE, ∂12CB
        t4, t5, t6 = -64*x1[inx]*x2[inx]*3, -64*x1[inx]*x2[inx]*2*normx2[inx], -64*x1[inx]*x2[inx]*dnormx4[inx]
        rtn∂12CE[inx]  = KE4[inx]*t4 + KE5[inx]*t5 + KE6[inx]*t6
        rtn∂12CB[inx]  = KB4[inx]*t4 + KB5[inx]*t5 + KB6[inx]*t6
        # ∂12SE, ∂12SB
        t4, t5, t6 =  -64*x1[inx]*x2[inx]*9, -64*x1[inx]*x2[inx]*6*normx2[inx], -64*x1[inx]*x2[inx]*4*sqx1[inx]*sqx2[inx]
        rtn∂12SE[inx]  = KE4[inx]*t4 + KE5[inx]*t5 + KE6[inx]*t6
        rtn∂12SB[inx]  = KB4[inx]*t4 + KB5[inx]*t5 + KB6[inx]*t6
        # ∂12SCE, ∂12SCB
        t4, t5, t6 = -32*dnormx2[inx]*3, -32*dnormx2[inx]*2*normx2[inx] , -32*dnormx2[inx]*2*2*sqx1[inx]*sqx2[inx]
        rtn∂12SCE[inx]  = KE4[inx]*t4 + KE5[inx]*t5 + KE6[inx]*t6
        rtn∂12SCB[inx]  = -KB4[inx]*t4 - KB5[inx]*t5 - KB6[inx]*t6

    end
    return  rtnCE, rtnCB,
            rtnSE, rtnSB,
            rtnSCE, rtnSCB,
            rtn∂11CE, rtn∂11CB,
            rtn∂11SE, rtn∂11SB,
            rtn∂11SCE, rtn∂11SCB,
            rtn∂22CE, rtn∂22CB,
            rtn∂22SE, rtn∂22SB,
            rtn∂22SCE, rtn∂22SCB,
            rtn∂12CE, rtn∂12CB,
            rtn∂12SE, rtn∂12SB,
            rtn∂12SCE, rtn∂12SCB
end

# ---- variogram computation
"""
Variogram computation for uncertainty in phi. The  `radKϕ` argument should be constructed by
```
vϕK1, ugrid = PrimordialB.vecϕKn(1, varϕk, g)
vϕK2, ugrid = PrimordialB.vecϕKn(2, varϕk, g)
radKϕ = hcat(ugrid, vϕK1, vϕK2)
```
"""
function varioϕx(x1, x2, radKϕ)
    sx1   = abs2(x1)
    sx2   = abs2(x2)
    normx = √(sx1 + sx2)

    sϕK1 = Dierckx.Spline1D(radKϕ[:,1], radKϕ[:,2]; k=4, bc="extrapolate", s=0.0)
    sϕK2 = Dierckx.Spline1D(radKϕ[:,1], radKϕ[:,3]; k=4, bc="extrapolate", s=0.0)
    Kϕ1 = reshape(sϕK1(vec(normx)), size(normx))
    Kϕ2 = reshape(sϕK2(vec(normx)), size(normx))

    rtn∂11, rtn∂12, rtn∂22 = zeros(normx), zeros(normx), zeros(normx)
    @inbounds for inx in eachindex(rtn∂11)
        rtn∂11[inx]  = Kϕ1[inx]*(-2) + Kϕ2[inx]*(-4*sx1[inx])
        rtn∂22[inx]  = Kϕ1[inx]*(-2) + Kϕ2[inx]*(-4*sx2[inx])
        rtn∂12[inx]  = Kϕ2[inx]*(-4*x1[inx]*x2[inx])
    end
    rtn∂11 -= sϕK1(0.0)*(-2)
    rtn∂22 -= sϕK1(0.0)*(-2)
    return -rtn∂11, -rtn∂22, -rtn∂12
end

#######################################################

#  Generates the covariance matrices, with marginalization over nϕ.

#######################################################

function makeΣen_Σb_withoutnoise(xdiff, ydiff, lnxdiff, lnydiff, radKE, radKB, radKϕ)
    ΣCE, ΣCB, ΣSE, ΣSB, ΣSCE, ΣSCB,
    Σ∂11CE, Σ∂11CB, Σ∂11SE, Σ∂11SB, Σ∂11SCE, Σ∂11SCB,
    Σ∂22CE, Σ∂22CB, Σ∂22SE, Σ∂22SB, Σ∂22SCE, Σ∂22SCB,
    Σ∂12CE, Σ∂12CB, Σ∂12SE, Σ∂12SB, Σ∂12SCE, Σ∂12SCB = PrimordialB.YXcov(lnxdiff, lnydiff, radKE, radKB)

    v11Σϕ, v22Σϕ, v12Σϕ = PrimordialB.varioϕx(xdiff, ydiff, radKϕ)

    Σce, Σcb   = Array{Float64}(size(ΣCE)), Array{Float64}(size(ΣCE))
    Σse, Σsb   = Array{Float64}(size(ΣCE)), Array{Float64}(size(ΣCE))
    Σsce, Σscb = Array{Float64}(size(ΣCE)), Array{Float64}(size(ΣCE))
    @inbounds @simd for inx in eachindex(Σce)
        Σce[inx]   = ΣCE[inx]  - v11Σϕ[inx]*Σ∂11CE[inx] - v22Σϕ[inx]*Σ∂22CE[inx] - 2v12Σϕ[inx]*Σ∂12CE[inx]
        Σcb[inx]   = ΣCB[inx]  - v11Σϕ[inx]*Σ∂11CB[inx] - v22Σϕ[inx]*Σ∂22CB[inx] - 2v12Σϕ[inx]*Σ∂12CB[inx]
        Σse[inx]   = ΣSE[inx]  - v11Σϕ[inx]*Σ∂11SE[inx] - v22Σϕ[inx]*Σ∂22SE[inx] - 2v12Σϕ[inx]*Σ∂12SE[inx]
        Σsb[inx]   = ΣSB[inx]  - v11Σϕ[inx]*Σ∂11SB[inx] - v22Σϕ[inx]*Σ∂22SB[inx] - 2v12Σϕ[inx]*Σ∂12SB[inx]
        Σsce[inx]  = ΣSCE[inx] - v11Σϕ[inx]*Σ∂11SCE[inx] - v22Σϕ[inx]*Σ∂22SCE[inx] - 2v12Σϕ[inx]*Σ∂12SCE[inx]
        Σscb[inx]  = ΣSCB[inx] - v11Σϕ[inx]*Σ∂11SCB[inx] - v22Σϕ[inx]*Σ∂22SCB[inx] - 2v12Σϕ[inx]*Σ∂12SCB[inx]
    end
    Σb    = [Σsb Σscb; Σscb Σcb]
    Σen   = [Σce Σsce; Σsce Σse]

    return Σen, Σb
end

function addΣnoise!(Σ, σEErad, deltx, σmod)
    @assert length(σmod) == size(Σ, 1)
    diag_noise = abs2(σmod .* σEErad / deltx)
    @inbounds for k = 1:size(Σ, 1)
        Σ[k,k] += diag_noise[k]
    end
    return Σ
end

function make_v1Σenv2_v1Σbv2(blk1, blk2, vqb1, vqb2, vub1, vub2, x1b1, x1b2, x2b1, x2b2, lnx1b1, lnx1b2, lnx2b1, lnx2b2, radKE, radKB, radKϕ,  σEErad, deltx, σmod)
    x1diff    = x1b1 .- (x1b2).'
    x2diff    = x2b1 .- (x2b2).'
    lnx1diff  = lnx1b1 .- (lnx1b2).'
    lnx2diff  = lnx2b1 .- (lnx2b2).'
    v1  = vcat(vqb1, vub1)
    v2  = vcat(vqb2, vub2)
    Σen, Σb  = makeΣen_Σb_withoutnoise(x1diff, x2diff, lnx1diff, lnx2diff, radKE, radKB, radKϕ)
    if blk1 == blk2
        addΣnoise!(Σen, σEErad, deltx, σmod)
    end
    vΣen = transpose(v1) * Σen * v2
    vΣb  = transpose(v1) * Σb  * v2
    if blk1 != blk2
        vΣen += transpose(v2) * transpose(Σen) * v1
        vΣb  += transpose(v2) * transpose(Σb)  * v1
    end
    return vΣen, vΣb
end

#######################################################

#  likelihood r profile with marginalizing

#######################################################
function loglike_r(rng, qudata, D, V, Db, Vb)
    VLd    = At_mul_B(V, At_mul_B(Vb, qudata) ./ √(Db))
    rout   = Array{Float64}(length(rng))
    for i=1:length(rng)
        @inbounds rout[i] = - 0.5 * dot(VLd, VLd./(D+rng[i])) - 0.5 * sum(log(D+rng[i])) - 0.5 * length(VLd) * log(2π)
    end
    return rout
end

function loglike_r(rng, qudata, Σen, Σb)
    return loglike_r(rng, qudata, D, V, Db, Vb)
end

function get_D_V_Db_Vb(Σen, Σb)
    Db, Vb = eig(Symmetric(Σb))
    Ipos1  = Db .> 0.0
    Db     = Db[Ipos1]
    Vb     = Vb[:,Ipos1]
    B      = transpose(Vb) * Σen * Vb
    B      = transpose(1./√Db) .* B
    B      = B .* (1./√Db)
    D,V    = eig(Symmetric(B))
    # ------
    #Ipos2  = 1e10 .> D .> 0.0
    #Ipos2  = sort(D, rev=true)[10] .> D .> 0.0
    Ipos2  =  D .> 0.0
    # -----
    D      = D[Ipos2]
    V      = V[:,Ipos2]
    return D, V, Db, Vb
end

#######################################################

#  likelihood r profile with marginalizing and cholesky

#######################################################
# function loglike_r_chol(rng, qudata, Σen, Σb)
#     Ub = chol(Σb)
#     UΣenU = (Ub') \ (Σen / Ub)
#     D, V  = eig(Symmetric(UΣenU))
#     VUd    = At_mul_B(V, At_mul_B(Ub, qudata))
#     rout   = Array{Float64}(length(rng))
#     for i=1:length(rng)
#         @inbounds rout[i] = - 0.5 * dot(VUd, VUd./(D+rng[i])) - 0.5 * sum(log(D+rng[i])) - 0.5 * length(VUd) * log(2π)
#     end
#     return rout
# end
function loglike_r_chol(rng, qudata, Σen, Σb)
    Db, Vb = eig((Σen + Σen')./2)
    Ipos1  = Db .> 0.0
    Db     = Db[Ipos1]
    Vb     = Vb[:,Ipos1]
    B      = transpose(Vb) * Σb * Vb
    B      = transpose(1./√Db) .* B
    B      = B .* (1./√Db)
    D,V    = eig((B + B')./2)
    # ------
    #Ipos2  = 1e10 .> D .> 0.0
    #Ipos2  = sort(D, rev=true)[10] .> D .> 0.0
    Ipos2  =  D .> 0.0
    # -----
    D      = D[Ipos2]
    V      = V[:,Ipos2]
    VLd    = At_mul_B(V, At_mul_B(Vb, qudata) ./ √(Db))
    rout   = Array{Float64}(length(rng))
    for i=1:length(rng)
        @inbounds rout[i] = - 0.5 * dot(VLd, VLd./(1.0+rng[i]*D)) - 0.5 * sum(log(1.0+rng[i]*D)) - 0.5 * length(VLd) * log(2π)
    end
    return rout
end

"""
Exact low ell log likelihood for `r` with marginalized delensing residuals.
```
rout = lowℓ_loglike_marg_profile(rng, dbk, lmax, g, mCls, mat_clnbk)
```
"""
function lminlmax_loglike_marg_profile(rng, dbk, lmin, lmax, g, mCls, mat_clnbk)
    hermitianhalf = (g.k[1] .> 0) | ((g.k[1] .== 0.0) & (g.k[2] .> 0.0))
    #lowellbox     = (-lmax .<= g.k[2] .<= lmax) & (-lmax .<= g.k[1] .<= lmax)
    #highellbox    = !( (-lmin .<= g.k[2] .<= lmin) & (-lmin .< g.k[1] .<= lmin) )
    #lrestriction  = lowellbox & highellbox & hermitianhalf
    lrestriction  = (lmin .<= g.r .<= lmax) & hermitianhalf
    cBB0kmax   = mCls.cBB0k[lrestriction]
    cNNkmax    = mCls.cBBnoisek[lrestriction] + mat_clnbk[lrestriction]
    dbkmax     = dbk[lrestriction]
    δ0         = 1/g.deltk^2
    rout       = Array{Float64}(length(rng))
    for i=1:length(rng)
        cddkmax  = (rng[i]*cBB0kmax + cNNkmax)
        rout[i]  = - sum(abs2(dbkmax) ./ cddkmax) / δ0
        rout[i] += - sum(log(cddkmax))
    end
    return rout
end

######################################

#  periodize coordinates (for computing cov matrices) ... currently not used

######################################

function periodize(x, period)
    cst = 2π / (period)
    return angle(exp(im * cst .* x)) ./ cst
end


text:
# 2D geometry NRL: cross-over is at 6% sparsity!!!
# 3D, NRL: ca. 0.06 as well
for N = 2:6
   at = (N,N,N) * bulk("Si", pbc=false, cubic=true)
   tbm = SKTB.NRLTB.NRLTBModel(:Si, FermiDiracSmearing(1.0))
   H, M = hamiltonian(tbm, at)
   Hs = sparse(H)
   println("N = $N | sparsity: $(nnz(Hs)/length(Hs)) %")
   println("  full LU:")
   @time lufact(H.data)
   @time lufact(H.data)
   @time lufact(H.data)
   println("  sparse LU:")
   @time lufact(Hs)
   @time lufact(Hs)
   @time lufact(Hs)
end


text:
using CCDReduction: parse_name,
                    generate_filename,
                    parse_name_ext,
                    writefits

@testset "fitscollection" begin
    # setting initial data
    dir = joinpath(@__DIR__, "data")

    # the two fits file used for testing have identical header keys
    df = fitscollection(dir)
    relevant_keys = filter(i -> i ∉ ("", "HISTORY"), keys(read_header(M6707HH[1])))

    header_M6707HH = read_header(M6707HH[1])
    relevant_values_M6707HH = map(relevant_keys) do key
        header_M6707HH[key]
    end

    header_M35070V = read_header(M35070V[1])
    relevant_values_M35070V = map(relevant_keys) do key
        header_M35070V[key]
    end
    # testing rows columns
    # an additional 3 columns are present because of varibles path, name and hdu
    @test size(df) == (2, length(relevant_keys) + 3)

    # testing each element in rows and columns
    @test copy(df[1, :]) == (path = joinpath(dir, "M35070V.fits"),
                             name = "M35070V.fits", hdu = 1,
                             zip(Symbol.(relevant_keys), relevant_values_M35070V)...)
    @test copy(df[2, :]) == (path = joinpath(dir, "M6707HH.fits"),
                             name = "M6707HH.fits", hdu = 1,
                             zip(Symbol.(relevant_keys), relevant_values_M6707HH)...)

    # testing special cases
    ## drop extension
    df = fitscollection(dir; keepext = false)
    @test df[1, :name] == "M35070V"
    @test df[2, :name] == "M6707HH"

    ## exclude
    df = fitscollection(dir; exclude = "M35070V.fits")
    @test df[1, :name] == "M6707HH.fits"
    @test size(df) == (1, length(relevant_keys) + 3)

    ## exclude_dir
    df = fitscollection(dir; exclude_dir = "data")
    @test size(df) == (0, 0)
end

@testset "array-generators" begin
    # setting initial data
    dir = joinpath(@__DIR__, "data")
    df = fitscollection(dir)
    arr1 = arrays(df) |> collect
    arr2 = map(eachrow(df)) do row
        getdata(FITS(row.path)[row.hdu])
    end
    @test arr1 == arr2
end

@testset "filename-generators" begin
    # setting initial data
    dir = joinpath(@__DIR__, "data")
    df = fitscollection(dir)
    arr1 = filenames(df) |> collect
    arr2 = map(eachrow(df)) do row
        row.path
    end
    @test arr1 == arr2
end

@testset "ccds-generators" begin
    # setting initial data
    dir = joinpath(@__DIR__, "data")
    df = fitscollection(dir)
    arr1 = ccds(df) |> collect
    arr2 = map(eachrow(df)) do row
        CCDData(FITS(row.path)[row.hdu])
    end

    for (hdu1, hdu2) in zip(arr1, arr2)
        @test hdu1.data == hdu2.data
        test_header(hdu1, hdu2)
    end
end

@testset "saving-ccds" begin
    # copy data to separate temp directory for testing
    source_dir = joinpath(@__DIR__, "data")
    dir = mktempdir()
    cp(source_dir, dir; force=true)

    collection = fitscollection(dir)

    mktempdir() do savedir
        final = ccds(collection; save = true, path = savedir, save_prefix = "test1", save_suffix = "test2") do img
            trim(img, (:, 1040:1059))
        end

        # testing function outputs
        @test final[1].data == trim(CCDData(M35070V[1]), (:, 1040:1059)).data
        @test final[2].data == trim(CCDData(M6707HH[1]), (:, 1040:1059)).data

        collection1 = fitscollection(savedir; recursive = false)

        # generating arrays from collection1
        new_saved_data = map(eachrow(collection1)) do row
            fh = FITS(row.path)
            data = getdata(fh[row.hdu])
            close(fh)
            data
        end

        # testing saved data
        for (ccd, arr2) in zip(final, new_saved_data)
            @test ccd.data == arr2
        end

        # testing saved filenames
        @test collection1[1, :name] == "test1_M35070V_test2.fits"
        @test collection1[2, :name] == "test1_M6707HH_test2.fits"
    end

    # testing save functionality with path = nothing, i.e. at same location as of the input file
    new_saved_data = ccds(collection; save_prefix = "save_without_path") do data
        data
    end
    collection2 = fitscollection(dir)
    @test first(size(collection2)) == 4
end

@testset "saving-filename" begin
    # copy data to separate temp directory for testing
    source_dir = joinpath(@__DIR__, "data")
    dir = mktempdir()
    cp(source_dir, dir; force=true)

    collection = fitscollection(dir)

    mktempdir() do savedir
        final = filenames(collection; save = true, path = savedir, save_prefix = "test1", save_suffix = "test2") do img
            getdata(FITS(img)[1])
        end

        # testing function outputs
        @test final[1] == getdata(M35070V[1])
        @test final[2] == getdata(M6707HH[1])

        collection1 = fitscollection(savedir; recursive = false)

        # generating arrays from collection1
        new_saved_data = map(eachrow(collection1)) do row
            fh = FITS(row.path)
            data = getdata(fh[row.hdu])
            close(fh)
            data
        end

        # testing saved data
        @test final == new_saved_data

        # testing saved filenames
        @test collection1[1, :name] == "test1_M35070V_test2.fits"
        @test collection1[2, :name] == "test1_M6707HH_test2.fits"
    end

    # testing save functionality with path = nothing, i.e. at same location as of the input file
    new_saved_data = filenames(collection; save_prefix = "save_without_path") do filename
        CCDData(filename)
    end
    collection2 = fitscollection(dir)
    @test first(size(collection2)) == 4
end

@testset "saving-arrays" begin
    # copy data to separate temp directory for testing
    source_dir = joinpath(@__DIR__, "data")
    dir = mktempdir()
    cp(source_dir, dir; force=true)

    collection = fitscollection(dir)

    mktempdir() do savedir
        final = arrays(collection; save = true, path = savedir, save_prefix = "test1", save_suffix = "test2") do img
            trim(img, (:, 1040:1059))
        end

        # testing function outputs
        @test final[1] == trim(getdata(M35070V[1]), (:, 1040:1059))
        @test final[2] == trim(getdata(M6707HH[1]), (:, 1040:1059))

        collection1 = fitscollection(savedir; recursive = false)

        # generating arrays from collection1
        new_saved_data = map(eachrow(collection1)) do row
            fh = FITS(row.path)
            data = getdata(fh[row.hdu])
            close(fh)
            data
        end

        # testing saved data
        @test final == new_saved_data

        # testing saved filenames
        @test collection1[1, :name] == "test1_M35070V_test2.fits"
        @test collection1[2, :name] == "test1_M6707HH_test2.fits"
    end

    # testing save functionality with path = nothing, i.e. at same location as of the input file
    new_saved_data = arrays(collection; save_prefix = "save_without_path") do data
        data
    end
    collection2 = fitscollection(dir)
    @test first(size(collection2)) == 4
end

@testset "helper" begin
    # testing parse_name
    @test parse_name("abc.fits", "."*"fits", Val(true)) == "abc.fits"
    @test parse_name("abc.fits.tar.gz", "."*"fits.tar.gz", Val(false)) == "abc"
    @test parse_name("foo.fits.fits", "."*"fits", Val(false)) == "foo.fits"
    @test parse_name("foo.fits.fits", "."*r"fits(\.tar\.gz)?"i, Val(false)) == "foo.fits"
    @test parse_name("foo.fits", "."*r"fits(\.tar\.gz)?"i, Val(false)) == "foo"

    # testing generate_filename
    @test generate_filename("home/abcd.fits", @__DIR__, "test1", "test2", "_", r"fits(\.tar\.gz)?"i) == joinpath(@__DIR__, "test1_abcd_test2.fits")
    @test generate_filename("home/tek/abcd.fits", @__DIR__, nothing, "test2", "_", r"fits(\.tar\.gz)?"i) == joinpath(@__DIR__, "abcd_test2.fits")
    @test generate_filename("home/downloads/abcd.fits", @__DIR__, "test1", nothing, "_", r"fits(\.tar\.gz)?"i) == joinpath(@__DIR__, "test1_abcd.fits")
    @test generate_filename("home/hello/abcd.fits", @__DIR__, nothing, nothing, "_", r"fits(\.tar\.gz)?"i) == joinpath(@__DIR__, "abcd.fits")
    @test generate_filename("~/.julia/abcd.fits", @__DIR__, "test1", nothing, "__", r"fits(\.tar\.gz)?"i) == joinpath(@__DIR__, "test1__abcd.fits")
    @test generate_filename("~/.julia/abcd.abcd", @__DIR__, "test1", nothing, "__", r"fits(\.tar\.gz)?"i) == joinpath(@__DIR__, "test1__abcd.abcd") # case when extension does not match

    # testing parse_filename_ext
    parse_name_ext("11.12.20_HD106754.fits.tar.gz", "." * r"fits(\.tar\.gz)?"i) == ("11.12.20_HD106754", ".fits.tar.gz")
    parse_name_ext("abcd", "." * r"fits(\.tar\.gz)?"i) ==  ("abcd", "")
    parse_name_ext("11.12.20_HD106754.fits", "." * "fits") == ("11.12.20_HD106754", ".fits")
    parse_name_ext("abcd.fits", "." * "fits") == ("abcd", ".fits")
    parse_name_ext("11.12.20_HD106754.Fits.Tar.gz", "." * r"fits(\.tar\.gz)?"i) == ("11.12.20_HD106754", ".Fits.Tar.gz")

    # testing writefits
    mktempdir() do path
        filename = joinpath(path, "test1.fits")
        sample_data = rand(5, 10)
        writefits(filename, sample_data)
        fh = FITS(filename)
        image_array = getdata(fh[1])
        @test image_array == sample_data
        close(fh) # closing handle so that generated file can be deleted

        # writitng CCDData
        filename = joinpath(path, "test2.fits")
        ccd = CCDData(zeros(4, 4))
        writefits(filename, ccd)
        fh = FITS(filename)
        image_array = getdata(fh[1])
        @test image_array == ccd.data
        @test read_header(fh[1])["SIMPLE"] == true
        close(fh)
    end
end


text:
using Docx, Test

if isempty(ARGS) || "all" in ARGS
	all_tests = true
else
	all_tests = false
end

if all_tests || "test.docx" in ARGS
	@testset "test.docx" begin
		docx_path = string(@__DIR__, "/test_files/test.docx")

		doc = Docx.open(docx_path)
		result = Docx.read(doc, String)

		@test result == """
		python-docx was here!
		python-docx was here too!
		"""
	end
end

if all_tests || "style_table_shape.docx" in ARGS
	@testset "style_table_shape.docx" begin
		docx_path = string(@__DIR__, "/test_files/style_table_shape.docx")

		doc = Docx.open(docx_path)
		result = Docx.read(doc, String)

		# TODO: Handle de-duplication of "Shape 1 Text"
		# TODO: Add support for producing a "markdown" table
		expected = """
Docx.jl Sample Document 2

This sample document contains various styles and objects that contain text.

Lorem Ipsum is simply dummy text of the printing and typesetting industry. Lorem Ipsum has been the industry's standard dummy text ever since the 1500s, when an unknown printer took a galley of type and scrambled it to make a type specimen book. It has survived not only five centuries, but also the leap into electronic typesetting, remaining essentially unchanged. It was popularised in the 1960s with the release of Letraset sheets containing Lorem Ipsum passages, and more recently with desktop publishing software like Aldus PageMaker including versions of Lorem Ipsum.

This is an italics paragraph.

This is a BOLD paragraph.

This is an underlined paragraph

Column 1
Column 2
Column 3
R1C1
R1C2
R1C3
R2C1
R2C2
R2C3

Shape 1 Text
Shape 1 Text
"""
		@test length(result) == length(expected)
		for i = 1:length(expected)
			@test result[i] == expected[i]
			@test codepoint(result[i]) ==  codepoint(expected[i])

			if codepoint(result[i]) != codepoint(expected[i])
				println("nth[$i]")
				println(expected[1:i])
				println("---")
				println(result[1:i])
				return
			end
		end
	end
end


text:
rettype(::Union{Cptr{CF}, Cbinding{CF, lib, name}}) where {CF<:Cfunction, lib, name} = rettype(CF)
argtypes(::Union{Cptr{CF}, Cbinding{CF, lib, name}}) where {CF<:Cfunction, lib, name} = argtypes(CF)
convention(::Union{Cptr{CF}, Cbinding{CF, lib, name}}) where {CF<:Cfunction, lib, name} = convention(CF)
rettype(::Type{CF}) where {retT, argTs, conv, CF<:Cfunction{retT, argTs, conv}} = retT
argtypes(::Type{CF}) where {retT, argTs, conv, CF<:Cfunction{retT, argTs, conv}} = argTs
convention(::Type{CF}) where {retT, argTs, conv, CF<:Cfunction{retT, argTs, conv}} = conv

# https://www.gnu.org/software/libc/manual/html_node/How-Variadic.html
# it appears that chars and ints are promoted to 32-bit ints and floats are promoted to doubles when used as varargs
cconvert_vararg(::Val{:cdecl}, ::Type{T}) where {T} = T
cconvert_vararg(::Val{:cdecl}, ::Type{T}) where {T<:Signed} = ifelse(sizeof(T) < sizeof(Cint), Cint, T)
cconvert_vararg(::Val{:cdecl}, ::Type{T}) where {T<:Unsigned} = ifelse(sizeof(T) < sizeof(Cuint), Cuint, T)
cconvert_vararg(::Val{:cdecl}, ::Type{T}) where {T<:Cchar} = ifelse(sizeof(T) < sizeof(Cint), Cint, T)
cconvert_vararg(::Val{:cdecl}, ::Type{T}) where {T<:AbstractFloat} = ifelse(sizeof(T) < sizeof(Cdouble), Cdouble, T)

# these define a type to add to the type tuple when ccalling variadic functions
cconvert_default(::Type{T}) where {T} = error("Sorry, `$(T)` is not yet usable for variadic functions, provide `CBinding.cconvert_default(::Type{$(T)}) = CConvertType` to define how to pass it to a variadic function call")
cconvert_default(::Type{T}) where {T<:AbstractString} = Cstring
cconvert_default(::Type{T}) where {T<:Char} = Cchar
cconvert_default(::Type{T}) where {T<:Union{Integer, AbstractFloat, Ref}} = T   # this covers Cptr too, Cptr <: Ref
cconvert_default(::Type{T}) where {E, T<:AbstractArray{E}} = Cptr{E}

(cb::Cbinding{<:Cfunction})(args...) = funccall(cb, args...)
(func::Cptr{<:Cfunction})(args...) = funccall(func, args...)

@generated function funccall(::Type{retT}, ::Type{argTs}, ::Val{conv}, func::Union{Cptr{<:Cfunction}, Cbinding{<:Cfunction, lib, name}}, args...) where {retT, argTs, conv, lib, name}
	return :(ccall(
		$(func <: Cptr ? :(Core.Intrinsics.bitcast(Ptr{Cvoid}, func)) : lib isa Symbol ? :(($(QuoteNode(name)), $(String(lib)))) : QuoteNode(name)),
		$((conv isa Symbol ? (conv,) : ())...),
		$(retT),
		($(((argT <: Cvariadic ? :(Ptr{Cvoid}...) : argT) for argT in argTs)...),),
		$((:(args[$(i)]) for i in eachindex(args))...)
	))
end

@generated function funccall(func::Union{Cptr{CF}, Cbinding{CF, lib, name}}, args...) where {retT, argTs, conv, CF<:Union{Cfunction{retT, argTs}, Cfunction{retT, argTs, conv}}, lib, name}
	function getbitstypes(Ts, ind)
		ind <= length(args) && Ts <: Tuple{} && return nothing  # more args, but no arg types
		ind > length(args) && return Ts <: Tuple{} ? () : Ts <: Tuple{Cvariadic} ? (Cvariadic,) : nothing  # no more args, are there more arg types left?

		# make up a new type if its variadic, otherwise keep plugging along
		H  = Ts <: Tuple{Cvariadic} ? :(cconvert_vararg(Val(conv), bitstype(cconvert_default($(args[ind]))))) : Base.tuple_type_head(Ts)
		Ts = Ts <: Tuple{Cvariadic} ? Ts : Base.tuple_type_tail(Ts)
		sig = getbitstypes(Ts, ind+1)
		return isnothing(sig) ? nothing : (:(bitstype($(H))), sig...)
	end

	rett  = :(bitstype(retT))
	argts = getbitstypes(argTs, 1)
	convt = conv isa Symbol ? :(conv) : nothing

	return isnothing(argts) ? quote
		throw(MethodError(func, args))
	end : quote
		return funccall($(rett), Tuple{$((argt for argt in argts)...)}, Val($(convt)), func, args...)
	end
end

Base.cconvert(::Type{Cptr{CF}}, func::Function) where {retT, argTs, CF<:Cfunction{retT, argTs}} = Cptr{CF}(func)
Base.cconvert(::Type{Cptr{CF}}, func::Base.CFunction) where {retT, argTs, CF<:Cfunction{retT, argTs}} = Cptr{CF}(Base.unsafe_convert(Ptr{Cvoid}, func))

Cptr{CF}(func::Function) where {CF<:Cfunction} = funccallback(CF, Val(func))
Cptr{CF}(func::Base.CFunction) where {CF<:Cfunction} = Cptr{CF}(func)

@generated function funccallback(::Type{retT}, ::Type{argTs}, ::Val{conv}, ::Val{func}) where {retT, argTs, conv, func}
	return :(Cptr{
		$(conv isa Symbol ? :(Cfunction{retT, argTs, conv}) : :(Cfunction{retT, argTs}))
	}(@cfunction(
		$(func),
		$(retT),
		($((argT for argT in argTs)...),))
	))
end

@generated function funccallback(::Type{CF}, ::Val{func}) where {retT, argTs, conv, CF<:Union{Cfunction{retT, argTs}, Cfunction{retT, argTs, conv}}, func}
	getbitstypes(Ts) = Ts <: Tuple{} ? () : (:(bitstype($(Base.tuple_type_head(Ts)))), getbitstypes(Base.tuple_type_tail(Ts))...)

	rett  = :(bitstype(retT))
	argts = getbitstypes(argTs)
	convt = conv isa Symbol ? :(conv) : nothing

	return quote
		isbitstype(typeof(func)) ?
			funccallback($(rett), Tuple{$((argt for argt in argts)...)}, Val($(convt)), Val($(func))) :
			error("Method `$(func)` is a closure, please use @cfunction to obtain a Base.CFunction first")
	end
end