text:
module Core

using UnstructuredGrids.Helpers
using UnstructuredGrids.Kernels

import Base: ndims
import Base: show
import Base: ==
import UnstructuredGrids.Kernels: generate_dual_connections
import UnstructuredGrids.Kernels: generate_cell_to_faces
import UnstructuredGrids.Kernels: find_cell_to_faces
import UnstructuredGrids.Kernels: generate_face_to_ftype
import UnstructuredGrids.Kernels: generate_face_to_vertices
import UnstructuredGrids.Kernels: generate_facet_to_isboundary
import UnstructuredGrids.Kernels: generate_face_to_isboundary

export Connections
export RefCell
export UGrid
export VERTEX
export list
export ptrs
export coordinates
export connections
export reffaces
export facetypes
export vtkid
export vtknodes
export celltypes
export refcells
export generate_ftype_to_refface
export generate_grid_graph
export generate_full_grid_graph

struct Connections{L<:AbstractVector{<:Integer},P<:AbstractVector{<:Integer}}
  list::L
  ptrs::P
end

list(c::Connections) = c.list

ptrs(c::Connections) = c.ptrs

function Connections(c::AbstractVector{<:AbstractVector{<:Integer}})
  list, ptrs = generate_data_and_ptrs(c)
  Connections(list,ptrs)
end

function show(io::IO,c::Connections)
  clist = list(c)
  cptrs = ptrs(c)
  ncells = length(cptrs)-1
  for cell in 1:ncells
    a = cptrs[cell]
    b = cptrs[cell+1]-1
    println(io,"$cell -> $(clist[a:b])")
  end
end

function (==)(a::Connections,b::Connections)
  a.list == b.list && a.ptrs == b.ptrs
end

struct RefCell
  ndims::Int
  faces::Vector{Vector{Vector{Int}}}
  facetypes::Vector{Vector{Int}}
  reffaces::Vector{Vector{RefCell}}
  coordinates::Array{Float64,2}
  vtkid::Int
  vtknodes::Vector{Int}
end

ndims(r::RefCell) = r.ndims

connections(r::RefCell, dim::Integer) =r.faces[dim+1]

facetypes(r::RefCell, dim::Integer) = r.facetypes[dim+1]

reffaces(r::RefCell, dim::Integer) =r.reffaces[dim+1]

coordinates(r::RefCell) = r.coordinates

vtkid(r::RefCell) = r.vtkid

vtknodes(r::RefCell) = r.vtknodes

function RefCell(;
  ndims::Int,
  faces::Vector{Vector{Vector{Int}}},
  facetypes::Vector{Vector{Int}} = fill(Int[],ndims),
  reffaces::Vector{Vector{RefCell}} = fill(RefCell[],ndims),
  coordinates::Array{Float64,2} = zeros(ndims,0),
  vtkid::Int = UNSET,
  vtknodes::Vector{Int} = Int[])

  @assert ndims == length(faces)
  @assert ndims == length(facetypes)
  @assert ndims == length(reffaces)
  @assert ndims == size(coordinates,1)

  RefCell( ndims, faces, facetypes, reffaces, coordinates, vtkid, vtknodes)

end

const VERTEX = RefCell(
  ndims = 0, faces = fill([Int[]],0), vtkid = 1, vtknodes = [1] )

struct UGrid{
  C<:Connections,
  T<:AbstractVector{<:Integer},
  X<:AbstractArray{<:Number,2}}

  cells::C
  celltypes::T
  refcells::Vector{RefCell}
  coordinates::X

end

connections(g::UGrid) = g.cells

celltypes(g::UGrid) = g.celltypes

refcells(g::UGrid) = g.refcells

coordinates(g::UGrid) = g.coordinates

function ndims(g::UGrid)
  maximum((ndims(r) for r in refcells(g)))
end

function refconnections(g::UGrid,dim::Integer)
  refconnections(refcells(g),dim)
end

function refconnections(refcells::Vector{RefCell}, dim::Integer)
  [ connections(refcell,dim) for refcell in refcells ]
end

function UGrid( cellsdata, cellsptrs, celltypes, refcells, coordinates)
  cells = Connections(cellsdata, cellsptrs)
  UGrid( cells, celltypes, refcells, coordinates)
end

function UGrid(r::RefCell;dim::Integer)
  cells = connections(r,dim)
  celltypes = facetypes(r,dim)
  refcells = reffaces(r,dim)
  coords = coordinates(r)
  UGrid( Connections(cells), celltypes, refcells, coords)
end

function UGrid(
  grid::UGrid,
  dim::Integer,
  vertex_to_cells=generate_dual_connections(connections(grid)),
  cell_to_faces=generate_cell_to_faces(dim,grid,vertex_to_cells))
  cell_to_vertices = connections(grid)
  cell_to_ctype = celltypes(grid)
  ctype_to_refcell = refcells(grid)
  ftype_to_refface, ctype_to_lface_to_ftype = generate_ftype_to_refface(dim,ctype_to_refcell)
  face_to_ftype = generate_face_to_ftype(
    cell_to_faces, cell_to_ctype, ctype_to_lface_to_ftype)
  face_to_vertices = generate_face_to_vertices(
    dim, cell_to_vertices, cell_to_faces, cell_to_ctype, ctype_to_refcell)
  point_to_coords = coordinates(grid)
  UGrid(face_to_vertices, face_to_ftype, ftype_to_refface, point_to_coords)
end

function generate_dual_connections(cell_to_faces::Connections)
  cell_to_faces_data = list(cell_to_faces)
  cell_to_faces_ptrs = ptrs(cell_to_faces)
  face_to_cells_data, face_to_cells_ptrs = generate_dual_connections(
  cell_to_faces_data, cell_to_faces_ptrs)
  Connections(face_to_cells_data, face_to_cells_ptrs)
end

function generate_cell_to_faces(
  dim::Integer,
  grid::UGrid,
  vertex_to_cells=generate_dual_connections(connections(grid)))
  cell_to_vertices = connections(grid)
  cell_to_ctype = celltypes(grid)
  ctype_to_refcell = refcells(grid)
  cell_to_faces = generate_cell_to_faces(
      cell_to_vertices, vertex_to_cells, cell_to_ctype, ctype_to_refcell,dim)
  cell_to_faces
end

function generate_cell_to_faces(
  cell_to_vertices::Connections,
  vertex_to_cells::Connections,
  cell_to_ctype::AbstractVector{<:Integer},
  ctype_to_refcell::Vector{RefCell},
  dim::Integer)

  cell_to_vertices_data = list(cell_to_vertices)
  cell_to_vertices_ptrs = ptrs(cell_to_vertices)

  vertex_to_cells_data = list(vertex_to_cells)
  vertex_to_cells_ptrs = ptrs(vertex_to_cells)

  ctype_to_lface_to_lvertices = refconnections(ctype_to_refcell, dim)

  cell_to_faces_data, cell_to_faces_ptrs = generate_cell_to_faces(
    cell_to_vertices_data,
    cell_to_vertices_ptrs,
    ctype_to_lface_to_lvertices,
    cell_to_ctype,
    vertex_to_cells_data,
    vertex_to_cells_ptrs)

  Connections(cell_to_faces_data, cell_to_faces_ptrs)

end

function find_cell_to_faces(
  cell_to_vertices::Connections,
  vertex_to_faces::Connections,
  cell_to_ctype::AbstractVector{<:Integer},
  ctype_to_refcell::Vector{RefCell},
  dim::Integer)

  cell_to_vertices_data = list(cell_to_vertices)
  cell_to_vertices_ptrs = ptrs(cell_to_vertices)

  vertex_to_faces_data = list(vertex_to_faces)
  vertex_to_faces_ptrs = ptrs(vertex_to_faces)

  ctype_to_lface_to_lvertices = refconnections( ctype_to_refcell, dim)

  cell_to_faces_data, cell_to_faces_ptrs = find_cell_to_faces(
    cell_to_vertices_data,
    cell_to_vertices_ptrs,
    ctype_to_lface_to_lvertices,
    cell_to_ctype,
    vertex_to_faces_data,
    vertex_to_faces_ptrs)

  Connections(cell_to_faces_data, cell_to_faces_ptrs)

end

function find_cell_to_faces(grid::UGrid, fgrid::UGrid)
  face_to_vertices = connections(fgrid)
  vertex_to_faces = generate_dual_connections(face_to_vertices)
  dim = ndims(fgrid)
  find_cell_to_faces(grid, vertex_to_faces, dim)
end

function find_cell_to_faces(
  grid::UGrid,
  vertex_to_faces::Connections,
  dim::Integer)

  cell_to_vertices = connections(grid)
  cell_to_ctype = celltypes(grid)
  ctype_to_refcell = refcells(grid)
  find_cell_to_faces(
    cell_to_vertices,
    vertex_to_faces,
    cell_to_ctype,
    ctype_to_refcell,
    dim)
end

function generate_face_to_vertices(
  dim::Integer,
  cell_to_vertices::Connections,
  cell_to_faces::Connections,
  cell_to_ctype::AbstractVector{<:Integer},
  ctype_to_refcell::Vector{RefCell})

  ctype_to_lface_to_lvertices = refconnections( ctype_to_refcell, dim)

  l, p = generate_face_to_vertices(
    list(cell_to_vertices),
    ptrs(cell_to_vertices),
    list(cell_to_faces),
    ptrs(cell_to_faces),
    cell_to_ctype,
    ctype_to_lface_to_lvertices)

  Connections(l,p)

end

function generate_face_to_ftype(
  cell_to_faces::Connections,
  cell_to_ctype::AbstractVector{<:Integer},
  ctype_to_lface_to_ftype::AbstractVector{<:AbstractVector{<:Integer}})
  generate_face_to_ftype(
    list(cell_to_faces),
    ptrs(cell_to_faces),
    cell_to_ctype,
    ctype_to_lface_to_ftype)
end

function generate_ftype_to_refface(dim,ctype_to_refcell)

  i_to_refface = Vector{RefCell}(undef,0)
  nctypes = length(ctype_to_refcell)
  ctype_to_lftype_to_i = Vector{Vector{Int}}(undef,nctypes)

  i = 1
  for (ctype,refcell) in enumerate(ctype_to_refcell)
    lftype_to_refface = reffaces(refcell,dim)
    lftype_to_i = Vector{Int}(undef,length(lftype_to_refface))
    for (lftype,refface) in enumerate(lftype_to_refface)
      push!(i_to_refface,refface)
      lftype_to_i[lftype] = i
      i +=1
    end
    ctype_to_lftype_to_i[ctype] = lftype_to_i
  end

  ftype_to_refface = unique(i_to_refface)
  i_to_ftype = indexin(i_to_refface,ftype_to_refface)

  ctype_to_lftype_to_ftype = copy(ctype_to_lftype_to_i)
  for ctype in 1:length(ctype_to_lftype_to_i)
    for lftype in 1:length(ctype_to_lftype_to_i[ctype])
      i = ctype_to_lftype_to_i[ctype][lftype]
      ftype = i_to_ftype[i]
      ctype_to_lftype_to_ftype[ctype][lftype] = ftype
    end
  end

  ctype_to_lface_to_ftype = Vector{Vector{Int}}(undef,nctypes)
  for (ctype,refcell) in enumerate(ctype_to_refcell)
    lface_to_lftype = facetypes(refcell,dim)
    lftype_to_ftype = ctype_to_lftype_to_ftype[ctype]
    lface_to_ftype = lftype_to_ftype[lface_to_lftype]
    ctype_to_lface_to_ftype[ctype] = lface_to_ftype
  end

  (ftype_to_refface, ctype_to_lface_to_ftype)

end

function generate_facet_to_isboundary(face_to_cells::Connections)
  face_to_cells_ptrs = ptrs(face_to_cells)
  generate_facet_to_isboundary(face_to_cells_ptrs)
end

function generate_face_to_isboundary(
  facet_to_isboundary::AbstractVector{Bool},
  face_to_facets::Connections)

  face_to_facets_data = list(face_to_facets)
  face_to_facets_ptrs = ptrs(face_to_facets)

  generate_face_to_isboundary(
    facet_to_isboundary,
    face_to_facets_data,
    face_to_facets_ptrs)

end

function generate_grid_graph(grid::UGrid)

  D = ndims(grid)
  cell_to_vertices = connections(grid)
  vertex_to_cells = generate_dual_connections(cell_to_vertices)

  T = typeof(cell_to_vertices)
  primal = Vector{T}(undef,D+1)
  dual = Vector{T}(undef,D+1)

  d = 0
  primal[d+1] = cell_to_vertices
  dual[d+1] = vertex_to_cells

  for d in 1:(D-1)
    cell_to_dfaces = generate_cell_to_faces( d, grid, vertex_to_cells)
    dfaces_to_cells = generate_dual_connections(cell_to_dfaces)
    primal[d+1] = cell_to_dfaces
    dual[d+1] = dfaces_to_cells
  end

  (primal, dual)

end

function generate_full_grid_graph(grid::UGrid)

  D = ndims(grid)
  data = Array{Connections}(undef,D+1,D+1)

  primal, dual = generate_grid_graph(grid)
  for d in 0:(D-1)
    data[D+1,d+1] = primal[d+1]
    data[d+1,D+1] = dual[d+1]
  end

  fgrids = [ UGrid(grid,d,dual[0+1],primal[d+1]) for d in 1:(D-1)]
  for d in 1:(D-1)
    fgrid = fgrids[d]
    face_to_vertices = connections(fgrid)
    vertex_to_faces = generate_dual_connections(face_to_vertices)
    data[d+1,0+1] = face_to_vertices
    data[0+1,d+1] = vertex_to_faces
  end

  for d in 1:(D-1)
    fgrid = fgrids[d]
    for j in 1:(d-1)
      vertex_to_jfaces = data[0+1,j+1]
      face_to_jfaces = find_cell_to_faces(fgrid, vertex_to_jfaces, j)
      jface_to_faces = generate_dual_connections(face_to_jfaces)
      data[d+1,j+1] = face_to_jfaces
      data[j+1,d+1] = jface_to_faces
    end
  end

  for d in 0:(D-1)
    dfaces_to_cells = data[d+1,D+1]
    data[d+1,d+1] = _identity_connections(dfaces_to_cells)
  end

  cell_to_vertices = data[D+1,0+1]
  data[D+1,D+1] = _identity_connections(cell_to_vertices)

  data

end

function _identity_connections(dfaces_to_cells)
    ndfaces = length(dfaces_to_cells.ptrs)-1
    Connections( [i for i in 1:ndfaces], [i for i in 1:(ndfaces+1)])
end

# Helpers

end # module Core


text:
using TestImages, FixedPointNumbers, Colors, AxisArrays
using Test, ReferenceTests, Suppressor
using ImageContrastAdjustment

using TestImages: shepp_logan # will be exported in future release

# make sure all files in remotefiles are valid image files
foreach(TestImages.remotefiles) do img_name
    img = testimage(img_name)
    @test isa(img, AbstractArray{<:Colorant})
end

include("utils.jl")

img = testimage("cameraman")
@test isa(img, Matrix{Gray{N0f8}})
img = testimage("mri-stack")
@test isa(img, AxisArray)
@test map(step, axisvalues(img)) == (1,1,5)
@test_nowarn testimage("c")

# mismatch handling
# Note: the behavior might change as `remotefiles` changes
err_str = @except_str testimage("nonexistence.png") ArgumentError
@test_reference "references/nonexistence_err.txt" err_str

err_str = @except_str testimage("leans") ArgumentError
@test_reference "references/leans_err.txt" err_str

err_str = @except_str testimage("abcd.png") ArgumentError
@test_reference "references/abcd_err.txt" err_str

err_str = @capture_err testimage("camereman")
@test_reference "references/camereman_warning.txt" split(err_str, "\n")[1]

# Shepp-Logan phantom
@testset "shepp_logan" begin
    @test eltype(shepp_logan(64)) == Gray{Float64}
    @test size(shepp_logan(64)) == (64, 64)
    @test size(shepp_logan(64, 128)) == (64, 128)
    @test shepp_logan(64; high_contrast=true) == shepp_logan(64)

    phantom_ct = shepp_logan(128; high_contrast=false)
    phantom_mri = shepp_logan(128; high_contrast=true)

    # CT value exceeds [0, 1] range, so we need to stretch its value to make references
    # FIXME: ReferenceTests is broken for *.png references
    @test_reference "references/shepp_logan_CT_128.txt" adjust_histogram(phantom_ct, LinearStretching())
    @test_reference "references/shepp_logan_MRI_128.txt" phantom_mri

    _norm(x) = sqrt(sum(abs2, x))
    P = [ 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0;
          0.0  0.0  1.0  0.2  0.2  1.0  0.0  0.0;
          0.0  0.0  0.2  0.3  0.3  0.2  0.0  0.0;
          0.0  0.0  0.2  0.0  0.2  0.2  0.0  0.0;
          0.0  0.0  0.2  0.0  0.0  0.2  0.0  0.0;
          0.0  0.0  0.2  0.2  0.2  0.2  0.0  0.0;
          0.0  0.0  1.0  0.2  0.2  1.0  0.0  0.0;
          0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0]
    Q = shepp_logan(8)
    @test _norm(P - Q) < 1e-10
    P = [ 0.0  0.0  0.0   0.0   0.0   0.0   0.0  0.0;
          0.0  0.0  2.0   1.02  1.02  2.0   0.0  0.0;
          0.0  0.0  1.02  1.03  1.03  1.02  0.0  0.0;
          0.0  0.0  1.02  1.0   1.02  1.02  0.0  0.0;
          0.0  0.0  1.02  1.0   1.0   1.02  0.0  0.0;
          0.0  0.0  1.02  1.02  1.02  1.02  0.0  0.0;
          0.0  0.0  2.0   1.02  1.02  2.0   0.0  0.0;
          0.0  0.0  0.0   0.0   0.0   0.0   0.0  0.0]
    Q = shepp_logan(8, high_contrast=false)
    @test _norm(P - Q) < 1e-10
end

ext: jl


text:
function peak_similarity(connectivity, pmap, tspan, init_cond)
  model = create_model(connectivity)
  @nonamespace u₀map = [model.a[1] => init_cond[1], model.a[2] => init_cond[2],
                        model.a[3] => init_cond[3]]
  prob = ODEProblem(model, u₀map, tspan, pmap)
  sol = solve(prob, lsoda())
  pks, ~ = calc_main_freq(sol, Int(5e3))
  return [pks[1] == pks[2], pks[2] == pks[3]]
end

# Goodwin Oscillator
connectivity = [0 0 1; 1 0 0; 0 -1 0]

# Test Goodwin's oscillatory solution
@nonamespace pmap_osci = (model.α[1] => 0.00004564, model.α[2] => 0.16342066,
                     model.α[3] => 0.70351496, model.β[1] => 56.67356362,
                     model.β[2] => 50.87885942, model.β[3] => 30.31904319,
                     model.γ[1] => 1657.39796806, model.γ[2] => 6988.78292809,
                     model.γ[3] => 7461.31128709, model.κ[1] => 0.11196140,
                     model.κ[2] => 0.50912113, model.κ[3] => 0.33210760,
                     model.η[1] => 3.67709117, model.η[2] => 7.17307667,
                     model.η[3] => 6.80130072)
init_cond_osci = [0.21, 0.23, 0.039]
# We expect all periodogram peaks to have the same frequency value for an
# oscillatory solution
@test all(peak_similarity(connectivity, pmap_osci, (0.0, 1.0), init_cond_osci))

# Test Goodwin's steady state solution
@nonamespace pmap_ss = (model.α[1] => 0.00004564, model.α[2] => 0.16342066,
                     model.α[3] => 0.70351496, model.β[1] => 6.67356362,
                     model.β[2] => 50.87885942, model.β[3] => 300.31904319,
                     model.γ[1] => 1657.39796806, model.γ[2] => 6988.78292809,
                     model.γ[3] => 7461.31128709, model.κ[1] => 0.11196140,
                     model.κ[2] => 0.50912113, model.κ[3] => 0.33210760,
                     model.η[1] => 3.67709117, model.η[2] => 7.17307667,
                     model.η[3] => 6.80130072)
init_cond_ss = [0.577, 0.99, 0.027]
# In a steady state solution we expect to see different frequencies for the
# peaks in the periodogram
@test !all(peak_similarity(connectivity, pmap_ss, (0.0, 1.0), init_cond_ss))


text:
"""
    isdense(M) -> Bool

Return true if a matrix is dense.

Note:
It is not exactly same as !isparse, e.g. Diagonal, IMatrix and PermMatrix are both not isdense and not isparse.
"""
function isdense end

isdense(M)::Bool = !issparse(M)
isdense(::Diagonal) = false
isdense(x::Transpose) = isdense(parent(x))
isdense(x::Adjoint) = isdense(parent(x))

"""faster invperm"""
function fast_invperm(order)
    v = similar(order)
    @inbounds @simd for i = 1:length(order)
        v[order[i]] = i
    end
    v
end

dropzeros!(A::Diagonal; trim::Bool = false) = A
nnz(M::Diagonal) = size(M, 1)
nonzeros(M::Diagonal) = M.diag
findnz(M::Diagonal) = (collect(1:size(M, 1)), collect(1:size(M, 1)), M.diag)


text:
"""
    function recallfull(ranks::AbstractArray, rtrue::AbstractArray)

Computes the relative number of matches in the input ranking with respect to the input ground truth
"""
function recallfull(ranks::AbstractArray, rtrue::AbstractArray)
    recall = 0
    nn, m = size(rtrue)
    @assert size(ranks)[1] == size(rtrue)[1]
    for i in 1:m
        recall += size(findall(in(ranks[:,i]),rtrue[:,i]))[1]
    end
    return recall / (nn*m)
end

"""
    recallN(ranks::AbstractArray, rtrue::AbstractArray, N)

Takes a parameter `N` and computes the N@N recall between the input ranking and the input ground truth.
"""
function recallN(ranks::AbstractArray, rtrue::AbstractArray, N)
    recall = 0
    nn, m = size(rtrue)
    for i in 1:m
        recall += size(findall(in(ranks[1:N,i]),rtrue[1:N,i]))[1]
    end
    return recall / (N*m)
end

"""
    recallN(ranks::AbstractArray, rtrue::AbstractArray, N)

Takes a parameter `N` and computes the 1@N recall between the input ranking and the input ground truth.
"""
function recall1atN(ranking::AbstractArray, rtrue::AbstractArray,N)
    score = 0
    n_queries = length(rtrue[1,:])
    for i in 1:n_queries
        if rtrue[1,i] in @view(ranking[1:N,i])
            score+=1
        end
    end
    return score / n_queries
end

"""
    approx_error(qd::QuantizerData, data::AbstractMatrix, rtrue::AbstractArray, queries::AbstractArray)

Takes the input approximation data generated by a constructor, the exact data set, the true nearest neighbors and a set of test queries.
Takes the nearest neighbor to each query, computes the approximated inner product with the query and the true inner product with the query
Returns the average normalized absolute between the two: `|<q, ̃x-x>-<q,x>|`
"""
function approx_error(qd::QuantizerData, data::AbstractMatrix, rtrue::AbstractArray, queries::AbstractArray)
    dist = zeros(size(queries)[2])
    for i in 1:size(queries)[2]
        qTx = dot(queries[:,i], data[:,rtrue[1,i]])
        qTx_approx = dot(queries[:,i], qd.I.B[rtrue[1,i]]'qd.C)
        dist[i] = abs((qTx - qTx_approx)/qTx)
    end
    return mean(dist)
end

"""
    function get1atNscores(yhat::AbstractMatrix, groundtruth::AbstractMatrix, n_neighbors::Int)

Computes the recall1atN score for every value between 1 and the input `n_neighbors` for the estimated
ranking and the groundtruth and outputs as a vector. Used for creating comparative 1@N plots.
"""
function get1atNscores(yhat::AbstractMatrix, groundtruth::AbstractMatrix, n_neighbors::Int)
    scores = zeros(n_neighbors)
    for i in 1:n_neighbors
        scores[i] = recall1atN(yhat, groundtruth, i)
    end
    return scores
end

export recallfull
export recallN
export recall1atN
export approx_error
export get1atNscores


text:
using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using IterativeSolvers
using LinearAlgebra

# The user can create his own orthogonalization function to use in IAR
function doubleGS_function!(VV, vv, h)
    h[:]=VV'*vv; vv[:]=vv-VV*h; g=VV'*vv; vv[:]=vv-VV*g;
    h[:] = h[:]+g[:]; β=norm(vv); vv[:]=vv/β; return β
end
# Then it is needed to create a type to access to this function
abstract type DoubleGS <: IterativeSolvers.OrthogonalizationMethod end
# And then introduce a function dispatch for this new type in order to use
# the defined orthogonalization function
import IterativeSolvers.orthogonalize_and_normalize!
function orthogonalize_and_normalize!(V,v,h,::Type{DoubleGS})
    doubleGS_function!(V, v, h) end

@testset "IAR" begin
    dep=nep_gallery("dep0");
    n=size(dep,1);

    @bench @testset "accuracy eigenpairs" begin
        (λ,Q)=iar(dep,σ=1.1,neigs=5,v=ones(n),
                  maxit=100,tol=eps()*100,errmeasure=ResidualErrmeasure(dep));
        verify_lambdas(5, dep, λ, Q, eps()*100)
    end

    @bench @testset "Solve by projection" begin
        (λ,Q)=iar(dep,σ=1.1,neigs=5,v=ones(n),
                  maxit=100,tol=eps()*100,errmeasure=ResidualErrmeasure(dep), proj_solve=true);
        verify_lambdas(5, dep, λ, Q, eps()*100)
    end

    @testset "Compute as many eigenpairs as possible (neigs=Inf)" begin
        (λ,Q)=iar(dep,σ=1.1,neigs=Inf,v=ones(n),
                  maxit=38,tol=eps()*100);
        verify_lambdas(6, dep, λ, Q, eps()*200)
    end

    @testset "orthogonalization" begin
        # NOW TEST DIFFERENT ORTHOGONALIZATION METHODS

        @bench @testset "DGKS" begin
            (λ,Q,V)=iar(dep,orthmethod=DGKS,σ=1.1,neigs=5,v=ones(n),maxit=100,tol=eps()*100)
            @test opnorm(V'*V - I) < 1e-6
        end

        @bench @testset "User provided doubleGS" begin
            (λ,Q,V)=iar(dep,orthmethod=DoubleGS,σ=1.1,neigs=5,v=ones(n),maxit=100,tol=eps()*100)
            @test opnorm(V'*V - I) < 1e-6
        end

        @bench @testset "ModifiedGramSchmidt" begin
            (λ,Q,V)=iar(dep,orthmethod=ModifiedGramSchmidt,σ=1.1,neigs=5,v=ones(n),maxit=100,tol=eps()*100)
            @test opnorm(V'*V - I) < 1e-6
        end

        @bench @testset "ClassicalGramSchmidt" begin
            (λ,Q,V)=iar(dep,orthmethod=ClassicalGramSchmidt,σ=1.1,neigs=5,v=ones(n),maxit=100,tol=eps()*100)
            @test opnorm(V'*V - I) < 1e-6
        end
    end

    @testset "Errors thrown" begin
        np=100;
        dep=nep_gallery("dep0",np);
        @test_throws NEPCore.NoConvergenceException (λ,Q)=iar(dep,σ=1.1,neigs=6,v=ones(np),
                  maxit=7,tol=eps()*100);
    end

end


text:
#=

approxChol Laplacian solver by Daniel A. Spielman, 2017.
This algorithm is an implementation of an approximate edge-by-edge elimination
algorithm inspired by the Approximate Gaussian Elimination algorithm of
Kyng and Sachdeva.

For usage exaples, see http://danspielman.github.io/Laplacians.jl/latest/usingSolvers/index.html

There are two versions of this solver:
one that fixes the order of elimination beforehand,
and one that adapts the order to eliminate verties of low degree.
These use different data structures.
LLOrdMat is for the fixed order, and LLmatp is for the adaptive order.

These coes produce a structure we call LDLinv that is then used in the solve.
The structure of this code is as follows:

The data structures appear in approxCholTypes.jl
We then have the outline:

* constructors for LLmatp and LLMatOrd
* get_ll_col and compress_ll_col : used inside the elimination
* approxChol : the main routine
* LDLsolver, and its forward and backward solve the apply LDLinv
* approxchol_lap: the main solver, which calls approxchol_lap1 on connected
    components.
    This then calls one of approxchol_lapWdeg, approxchol_lapGiven or approxchol_lapGreedy,
    depending on the parameters.

* approxchol_lapChol - for producing a Cholesky factor instead of an LDLinv.
  might be useful if optimized.
* data structures that are used for the adaptive low-degree version to
  choose the next vertex.

=#

"""
    params = ApproxCholParams(order, output)
order can be one of
* :deg (by degree, adaptive),
* :wdeg (by original wted degree, nonadaptive),
* :given
"""
mutable struct ApproxCholParams
    order::Symbol
    stag_test::Integer
end

ApproxCholParams() = ApproxCholParams(:deg, 5)
ApproxCholParams(sym::Symbol) = ApproxCholParams(sym, 5)

LDLinv(a::SparseMatrixCSC{Tval,Tind}) where {Tind,Tval} =
  LDLinv(zeros(Tind,a.n-1), zeros(Tind,a.n),Tind[],Tval[],zeros(Tval,a.n))

LDLinv(a::LLMatOrd{Tind,Tval}) where {Tind,Tval} =
  LDLinv(zeros(Tind,a.n-1), zeros(Tind,a.n),Tind[],Tval[],zeros(Tval,a.n))

LDLinv(a::LLmatp{Tind,Tval}) where {Tind,Tval} =
  LDLinv(zeros(Tind,a.n-1), zeros(Tind,a.n),Tind[],Tval[],zeros(Tval,a.n))

function LLmatp(a::SparseMatrixCSC{Tval,Tind}) where {Tind,Tval}
    n = size(a,1)
    m = nnz(a)

    degs = zeros(Tind,n)

    flips = flipIndex(a)

    cols = Array{LLp{Tind,Tval}}(undef, n)
    llelems = Array{LLp{Tind,Tval}}(undef, m)

    @inbounds for i in 1:n
        degs[i] = a.colptr[i+1] - a.colptr[i]

        ind = a.colptr[i]
        j = a.rowval[ind]
        v = a.nzval[ind]
        llpend = LLp{Tind,Tval}(j,v)
        next = llelems[ind] = llpend
        for ind in (a.colptr[i]+one(Tind)):(a.colptr[i+1]-one(Tind))
            j = a.rowval[ind]
            v = a.nzval[ind]
            next = llelems[ind] = LLp{Tind,Tval}(j,v,next)
        end
        cols[i] = next
    end

    @inbounds for i in 1:n
        for ind in a.colptr[i]:(a.colptr[i+1]-one(Tind))
            llelems[ind].reverse = llelems[flips[ind]]
        end
    end

    return LLmatp{Tind,Tval}(n, degs, cols, llelems)
end

"""
  Print a column in an LLmatp matrix.
  This is here for diagnostics.
"""
function print_ll_col(llmat::LLmatp, i::Int)
    ll = llmat.cols[i]
    println("col $i, row $(ll.row) : $(ll.val)")

    while ll.next != ll
        ll = ll.next
        println("col $i, row $(ll.row) : $(ll.val)")
    end
end

function LLMatOrd(a::SparseMatrixCSC{Tval,Tind}) where {Tind,Tval}
    n = size(a,1)
    m = nnz(a)

    cols = zeros(Tind, n)
    llelems = Array{LLord{Tind,Tval}}(undef, m)

    ptr = one(Tind)

    @inbounds for i in Tind(1):Tind(n-1)
        next = zero(Tind)

        for ind in (a.colptr[i]):(a.colptr[i+1]-one(Tind))
            j = a.rowval[ind]
            if (i < j)

              v = a.nzval[ind]
              llelems[ptr] = LLord{Tind,Tval}(j, next, v)
              next = ptr
              ptr += one(Tind)

            end
        end
        cols[i] = next
    end

    return LLMatOrd{Tind,Tval}(n, cols, llelems)
end

function LLMatOrd(a::SparseMatrixCSC{Tval,Tind}, perm::Array) where {Tind,Tval}
    n = size(a,1)
    m = nnz(a)

    invp = invperm(perm)

    cols = zeros(Tind, n)
    llelems = Array{LLord{Tind,Tval}}(undef, m)

    ptr = one(Tind)

    @inbounds for i0 in Tind(1):Tind(n)
        i = invp[i0]
        next = zero(Tind)

        for ind in (a.colptr[i0]):(a.colptr[i0+1]-one(Tind))
            j = invp[a.rowval[ind]]
            if (i < j)

              v = a.nzval[ind]
              llelems[ptr] = LLord{Tind,Tval}(j, next, v)
              next = ptr
              ptr += one(ptr)

            end
        end
        cols[i] = next
    end

    return LLMatOrd{Tind,Tval}(n, cols, llelems)
end

"""
  Print a column in an LLMatOrd matrix.
  This is here for diagnostics.
"""
function print_ll_col(llmat::LLMatOrd, i::Int)
    ptr = llmat.cols[i]
    while ptr != 0
      ll = llmat.lles[ptr]
      println("col $i, row $(ll.row) : $(ll.val)")

      ptr = ll.next
    end
end

#=============================================================

The approximate factorization

=============================================================#

function get_ll_col(llmat::LLmatp{Tind,Tval},
  i,
  colspace::Vector{LLp{Tind,Tval}}) where {Tind,Tval}

    ll = llmat.cols[i]
    len = 0
    @inbounds while ll.next != ll

        if ll.val > zero(Tval)
            len = len+1
            if (len > length(colspace))
                push!(colspace,ll)
            else
                colspace[len] = ll
            end
        end

        ll = ll.next
    end

    if ll.val > zero(Tval)
        len = len+1
        if (len > length(colspace))
            push!(colspace,ll)
        else
            colspace[len] = ll
        end
    end

    return len
end

function get_ll_col(llmat::LLMatOrd{Tind,Tval},
  i,
  colspace::Vector{LLcol{Tind,Tval}}) where {Tind,Tval}

    ptr = llmat.cols[i]
    len = 0
    @inbounds while ptr != 0

        #if ll.val > 0
            len = len+1

            # should not be an lles - is an abuse
            item = LLcol(llmat.lles[ptr].row, ptr, llmat.lles[ptr].val)
            if (len > length(colspace))
                push!(colspace,item)
            else
                colspace[len] = item
            end
        #end

        ptr = llmat.lles[ptr].next
    end

    return len
end

function compressCol!(a::LLmatp{Tind,Tval},
  colspace::Vector{LLp{Tind,Tval}},
  len::Int,
  pq::ApproxCholPQ{Tind}) where {Tind,Tval}

    o = Base.Order.ord(isless, x->x.row, false, Base.Order.Forward)

    sort!(colspace, 1, len, QuickSort, o)

    ptr = 0
    currow::Tind = 0

    c = colspace

    @inbounds for i in 1:len

        if c[i].row != currow
            currow = c[i].row
            ptr = ptr+1
            c[ptr] = c[i]

        else
            c[ptr].val = c[ptr].val + c[i].val
            c[i].reverse.val = zero(Tval)

            approxCholPQDec!(pq, currow)
        end
    end

    o = Base.Order.ord(isless, x->x.val, false, Base.Order.Forward)
    sort!(colspace, 1, ptr, QuickSort, o)

    return ptr
end

function compressCol!(
  colspace::Vector{LLcol{Tind,Tval}},
  len::Int
  ) where {Tind,Tval}

    o = Base.Order.ord(isless, x->x.row, false, Base.Order.Forward)

    sort!(colspace, one(len), len, QuickSort, o)

    c = colspace

    ptr = 0
    currow = c[1].row
    curval = c[1].val
    curptr = c[1].ptr

    @inbounds for i in 2:len

        if c[i].row != currow

            ptr = ptr+1
            c[ptr] = LLcol(currow, curptr, curval)  # next is abuse here: reall keep where it came from.

            currow = c[i].row
            curval = c[i].val
            curptr = c[i].ptr

        else

            curval = curval + c[i].val

        end

    end

    # emit the last row

    ptr = ptr+1
    c[ptr] = LLcol(currow, curptr, curval)

    o = Base.Order.ord(isless, x->x.val, false, Base.Order.Forward)
    sort!(colspace, one(ptr), ptr, QuickSort, o)

    return ptr
end

function approxChol(a::LLMatOrd{Tind,Tval}) where {Tind,Tval}
    n = a.n

    # need to make custom one without col info later.
    ldli = LDLinv(a)
    ldli_row_ptr = one(Tind)

    d = zeros(Tval,n)

    colspace = Array{LLcol{Tind,Tval}}(undef, n)
    cumspace = Array{Tval}(undef, n)
    #vals = Array(Tval,n) # will be able to delete this

    o = Base.Order.ord(isless, identity, false, Base.Order.Forward)

    for i in Tind(1):Tind(n-1)

        ldli.col[i] = i  # will get rid of this with new data type
        ldli.colptr[i] = ldli_row_ptr

        len = get_ll_col(a, i, colspace)

        len = compressCol!(colspace, len)

        csum = zero(Tval)
        for ii in 1:len
            #vals[ii] = colspace[ii].val    # if immut, no need for vals
            csum = csum + colspace[ii].val
            cumspace[ii] = csum
        end
        wdeg = csum

        colScale = one(Tval)

        for joffset in 1:(len-1)

            llcol = colspace[joffset]
            w = llcol.val * colScale
            j = llcol.row

            f = w/(wdeg)

            #vals[joffset] = zero(Tval)

            # kind = Laplacians.blockSample(vals,k=1)[1]
            r = rand() * (csum - cumspace[joffset]) + cumspace[joffset]
            koff = searchsortedfirst(cumspace,r,one(len),len,o)

            k = colspace[koff].row

            newEdgeVal = w*(one(Tval)-f)

            # create edge (j,k) with newEdgeVal
            # do it by reassigning ll
            if j < k # put it in col j
                jhead = a.cols[j]
                a.lles[llcol.ptr] = LLord(k, jhead, newEdgeVal)
                #ll.next = jhead
                #ll.val = newEdgeVal
                #ll.row = k
                a.cols[j] = llcol.ptr
            else # put it in col k
              khead = a.cols[k]
              a.lles[llcol.ptr] = LLord(j, khead, newEdgeVal)
              #ll.next = khead
              #ll.val = newEdgeVal
              #ll.row = j
              a.cols[k] = llcol.ptr
            end

            colScale = colScale*(one(Tval)-f)
            #wdeg = wdeg*(1.0-f)^2
            wdeg = wdeg - 2w + w^2/wdeg

            push!(ldli.rowval,j)
            push!(ldli.fval, f)
            ldli_row_ptr = ldli_row_ptr + one(Tind)

            # push!(ops, IJop(i,j,1-f,f))  # another time suck

        end # for joffset

        llcol = colspace[len]
        w = llcol.val * colScale
        j = llcol.row

        push!(ldli.rowval,j)
        push!(ldli.fval, one(Tval))
        ldli_row_ptr = ldli_row_ptr + one(Tind)

        d[i] = w

    end # for i

    ldli.colptr[n] = ldli_row_ptr

    ldli.d = d

    return ldli
end

# this one is greedy on the degree - also a big win
function approxChol(a::LLmatp{Tind,Tval}) where {Tind,Tval}
    n = a.n

    ldli = LDLinv(a)
    ldli_row_ptr = one(Tind)

    d = zeros(n)

    pq = ApproxCholPQ(a.degs)

    it = 1

    colspace = Array{LLp{Tind,Tval}}(undef, n)
    cumspace = Array{Tval}(undef, n)
    vals = Array{Tval}(undef, n) # will be able to delete this

    o = Base.Order.ord(isless, identity, false, Base.Order.Forward)

    @inbounds while it < n

        i = approxCholPQPop!(pq)

        ldli.col[it] = i # conversion!
        ldli.colptr[it] = ldli_row_ptr

        it = it + 1

        len = get_ll_col(a, i, colspace)

        len = compressCol!(a, colspace, len, pq)  #3hog

        csum = zero(Tval)
        for ii in 1:len
            vals[ii] = colspace[ii].val
            csum = csum + colspace[ii].val
            cumspace[ii] = csum
        end
        wdeg = csum

        colScale = one(Tval)

        for joffset in 1:(len-1)

            ll = colspace[joffset]
            w = vals[joffset] * colScale
            j = ll.row
            revj = ll.reverse

            f = w/(wdeg)

            vals[joffset] = zero(Tval)

            # kind = Laplacians.blockSample(vals,k=1)[1]
            r = rand() * (csum - cumspace[joffset]) + cumspace[joffset]
            koff = searchsortedfirst(cumspace,r,one(len),len,o)

            k = colspace[koff].row

            approxCholPQInc!(pq, k)

            newEdgeVal = f*(one(Tval)-f)*wdeg

            # fix row k in col j
            revj.row = k   # dense time hog: presumably becaus of cache
            revj.val = newEdgeVal
            revj.reverse = ll

            # fix row j in col k
            khead = a.cols[k]
            a.cols[k] = ll
            ll.next = khead
            ll.reverse = revj
            ll.val = newEdgeVal
            ll.row = j

            colScale = colScale*(one(Tval)-f)
            wdeg = wdeg*(one(Tval)-f)^2

            push!(ldli.rowval,j)
            push!(ldli.fval, f)
            ldli_row_ptr = ldli_row_ptr + one(Tind)

            # push!(ops, IJop(i,j,1-f,f))  # another time suck

        end # for

        ll = colspace[len]
        w = vals[len] * colScale
        j = ll.row
        revj = ll.reverse

        if it < n
            approxCholPQDec!(pq, j)
        end

        revj.val = zero(Tval)

        push!(ldli.rowval,j)
        push!(ldli.fval, one(Tval))
        ldli_row_ptr = ldli_row_ptr + one(Tind)

        d[i] = w

    end

    ldli.colptr[it] = ldli_row_ptr

    ldli.d = d

    return ldli
end

#=============================================================

The routines that do the solve.

=============================================================#

function LDLsolver(ldli::LDLinv, b::Vector)
    y = copy(b)

    forward!(ldli, y)

    @inbounds for i in 1:(length(ldli.d))
        if ldli.d[i] != 0
            y[i] /= ldli.d[i]
        end
    end

    backward!(ldli, y)

    mu = mean(y)
    @inbounds for i in eachindex(y)
        y[i] = y[i] - mu
    end

    return y
end

# Must return y because of Krylov/LinearOperators usage
function LDLsolver!(y::Vector, ldli::LDLinv, b::Vector)
    y .= b

    forward!(ldli, y)

    @inbounds for i in 1:(length(ldli.d))
        if ldli.d[i] != 0
            y[i] /= ldli.d[i]
        end
    end

    backward!(ldli, y)

    mu = mean(y)
    @inbounds for i in eachindex(y)
        y[i] = y[i] - mu
    end

    return y
end

function forward!(ldli::LDLinv{Tind,Tval}, y::Vector) where {Tind,Tval}

    @inbounds for ii in 1:length(ldli.col)
        i = ldli.col[ii]

        j0 = ldli.colptr[ii]
        j1 = ldli.colptr[ii+1]-one(Tind)

        yi = y[i]

        for jj in j0:(j1-1)
            j = ldli.rowval[jj]
            y[j] += ldli.fval[jj] * yi
            yi *= (one(Tval)-ldli.fval[jj])
        end
        j = ldli.rowval[j1]
        y[j] += yi
        y[i] = yi
    end
end

function backward!(ldli::LDLinv{Tind,Tval}, y::Vector) where {Tind,Tval}
    o = one(Tind)
    @inbounds for ii in length(ldli.col):-1:1
        i = ldli.col[ii]

        j0 = ldli.colptr[ii]
        j1 = ldli.colptr[ii+1]-o

        j = ldli.rowval[j1]
        yi = y[i]
        yi = yi + y[j]

        for jj in (j1-o):-o:j0
            j = ldli.rowval[jj]
            yi = (one(Tval)-ldli.fval[jj])*yi + ldli.fval[jj]*y[j]
        end
        y[i] = yi
    end
end

#=
  An attempt at an efficient solver for the case when y is a matrix.
  Have not yet found a meaningful speedup

function LDLsolver(ldli::LDLinv, b::Matrix)
    y = copy(b)

    (d, n) = size(y)
    @assert n == length(ldli.col)+1

    forward!(ldli, y)

    @inbounds for i in 1:(length(ldli.d))
        if ldli.d[i] != 0
          @simd for j in 1:d
            y[j,i] = y[j,i] / ldli.d[i]
          end
        end
    end

    backward!(ldli, y)

    @inbounds for j in 1:size(y,1)
        mu = mean(y[j,:])

        for i in 1:size(y,2)
            y[j,i] = y[j,i] - mu
        end
    end

    return y
end

function forward!{Tind,Tval}(ldli::LDLinv{Tind,Tval}, y::Matrix)

    (d, n) = size(y)
    @assert n == length(ldli.col)+1

    #yi = zeros(y[:,1])

    @inbounds for ii in 1:length(ldli.col)
        i = ldli.col[ii]

        j0 = ldli.colptr[ii]
        j1 = ldli.colptr[ii+1]-one(Tind)

        for jj in j0:(j1-1)
            j = ldli.rowval[jj]
            @simd for k in 1:d
              y[k,j] = y[k,j] + (ldli.fval[jj] * y[k,i])
              y[k,i] = y[k,i] * (one(Tval)-ldli.fval[jj])
          end
        end
        j = ldli.rowval[j1]

        @simd for k in 1:d
          y[k,j] = y[k,j] + y[k,i]
        end
    end
end

function backward!{Tind,Tval}(ldli::LDLinv{Tind,Tval}, y::Matrix)
    o = one(Tind)

    (d, n) = size(y)
    @assert n == length(ldli.col)+1

    yi = zeros(y[:,1])

    @inbounds for ii in length(ldli.col):-1:1
        i = ldli.col[ii]

        j0 = ldli.colptr[ii]
        j1 = ldli.colptr[ii+1]-o

        j = ldli.rowval[j1]
        #copy!(yi, y[:,i])

        @simd for k in 1:d
          y[k,i] = y[k,i] + y[k,j]
        end

        for jj in (j1-o):-o:j0
            j = ldli.rowval[jj]
            @simd for k in 1:d
              y[k,i] = (one(Tval)-ldli.fval[jj])*y[k,i] + ldli.fval[jj].*y[k,j]
            end
        end
        #y[:,i] = yi
    end
end

=#

"""
    solver = approxchol_lap(a); x = solver(b);
    solver = approxchol_lap(a; tol::Real=1e-6, maxits=1000, maxtime=Inf, verbose=false, pcgIts=Int[], params=ApproxCholParams())

A heuristic by Daniel Spielman inspired by the linear system solver in https://arxiv.org/abs/1605.02353 by Rasmus Kyng and Sushant Sachdeva.  Whereas that paper eliminates vertices one at a time, this eliminates edges one at a time.  It is probably possible to analyze it.
The `ApproxCholParams` let you choose one of three orderings to perform the elimination.

* ApproxCholParams(:given) - in the order given.
    This is the fastest for construction the preconditioner, but the slowest solve.
* ApproxCholParams(:deg) - always eliminate the node of lowest degree.
    This is the slowest build, but the fastest solve.
* ApproxCholParams(:wdeg) - go by a perturbed order of wted degree.

For more info, see http://danspielman.github.io/Laplacians.jl/latest/usingSolvers/index.html
"""
function approxchol_lap(a::SparseMatrixCSC{Tv,Ti};
  tol::Real=1e-6,
  maxits=1000,
  maxtime=Inf,
  verbose=false,
  pcgIts=Int[],
  params=ApproxCholParams()) where {Tv,Ti}

  if minimum(a.nzval) < 0
      error("Adjacency matrix can not have negative edge weights")
  end

    return Laplacians.lapWrapComponents(approxchol_lap1, a,
    verbose=verbose,
    tol=tol,
    maxits=maxits,
    maxtime=maxtime,
    pcgIts=pcgIts,
    params=params)

end

function approxchol_lapGreedy(a::SparseMatrixCSC;
  tol::Real=1e-6,
  maxits=1000,
  maxtime=Inf,
  verbose=false,
  pcgIts=Int[],
  params=ApproxCholParams())

  tol_ =tol
  maxits_ =maxits
  maxtime_ =maxtime
  verbose_ =verbose
  pcgIts_ =pcgIts

  t1 = time()

  la = lap(a) # a hit !?

  llmat = LLmatp(a)
  ldli = approxChol(llmat)
  F(b) = LDLsolver(ldli, b)

  if verbose
    println("Using greedy degree ordering. Factorization time: ", time()-t1)
    println("Ratio of operator edges to original edges: ", 2 * length(ldli.fval) / nnz(a))
  end

  if verbose
      println("ratio of max to min diagonal of laplacian : ", maximum(diag(la))/minimum(diag(la)))
  end

  f(b;tol=tol_,maxits=maxits_, maxtime=maxtime_, verbose=verbose_, pcgIts=pcgIts_) = pcg(la, b .- mean(b), F, tol=tol, maxits=maxits, maxtime=maxtime, pcgIts=pcgIts, verbose=verbose, stag_test = params.stag_test)

end

function approxchol_lapGiven(a::SparseMatrixCSC;
  tol::Real=1e-6,
  maxits=1000,
  maxtime=Inf,
  verbose=false,
  pcgIts=Int[],
  params=ApproxCholParams())

  tol_ =tol
  maxits_ =maxits
  maxtime_ =maxtime
  verbose_ =verbose
  pcgIts_ =pcgIts

  t1 = time()

  la = lap(a)

  llmat = LLMatOrd(a)
  ldli = approxChol(llmat)
  F(b) = LDLsolver(ldli, b)

  if verbose
    println("Using given ordering. Factorization time: ", time()-t1)
    println("Ratio of operator edges to original edges: ", 2 * length(ldli.fval) / nnz(a))
  end

  if verbose
      println("ratio of max to min diagonal of laplacian : ", maximum(diag(la))/minimum(diag(la)))
  end

  f(b;tol=tol_,maxits=maxits_, maxtime=maxtime_, verbose=verbose_, pcgIts=pcgIts_) = pcg(la, b .- mean(b), F, tol=tol, maxits=maxits, maxtime=maxtime, pcgIts=pcgIts, verbose=verbose, stag_test = params.stag_test)

end

function approxchol_lapWdeg(a::SparseMatrixCSC;
  tol::Real=1e-6,
  maxits=1000,
  maxtime=Inf,
  verbose=false,
  pcgIts=Int[],
  params=ApproxCholParams())

  tol_ =tol
  maxits_ =maxits
  maxtime_ =maxtime
  verbose_ =verbose
  pcgIts_ =pcgIts

  t1 = time()

  la = lap(a)

  v = vec(sum(a,dims=1))
  v = v .* (1 .+ rand(length(v)))
  p = sortperm(v)

  llmat = LLMatOrd(a,p)
  ldli = approxChol(llmat)

  ip = invperm(p)
  ldlip = LDLinv(p[ldli.col], ldli.colptr, p[ldli.rowval], ldli.fval, ldli.d[ip]);

  F = function(b)
    x = zeros(size(b))
    x = LDLsolver(ldlip, b)
    #x[p] = LDLsolver(ldli, b[p])
    return x
  end

  if verbose
    println("Using wted degree ordering. Factorization time: ", time()-t1)
    println("Ratio of operator edges to original edges: ", 2 * length(ldli.fval) / nnz(a))
  end

  if verbose
      println("ratio of max to min diagonal of laplacian : ", maximum(diag(la))/minimum(diag(la)))
  end

  f(b;tol=tol_,maxits=maxits_, maxtime=maxtime_, verbose=verbose_, pcgIts=pcgIts_) = pcg(la, b .- mean(b), F, tol=tol, maxits=maxits, maxtime=maxtime, pcgIts=pcgIts, verbose=verbose, stag_test = params.stag_test)

end

function approxchol_lap1(a::SparseMatrixCSC{Tv,Ti};
  tol::Real=1e-6,
  maxits=1000,
  maxtime=Inf,
  verbose=false,
  pcgIts=Int[],
  params=ApproxCholParams()) where {Tv,Ti}

    tol_ =tol
    maxits_ =maxits
    maxtime_ =maxtime
    verbose_ =verbose
    pcgIts_ =pcgIts

    if params.order == :deg

      return approxchol_lapGreedy(a,
        verbose=verbose,
        tol=tol,
        maxits=maxits,
        maxtime=maxtime,
        pcgIts=pcgIts,
        params=params)

    elseif params.order == :wdeg

      return approxchol_lapWdeg(a,
        verbose=verbose,
        tol=tol,
        maxits=maxits,
        maxtime=maxtime,
        pcgIts=pcgIts,
        params=params)

    else
      return approxchol_lapGiven(a,
        verbose=verbose,
        tol=tol,
        maxits=maxits,
        maxtime=maxtime,
        pcgIts=pcgIts,
        params=params)

    end

end

# TODO: remove later
#=
"""
    solver = approxchol_sddm(sddm); x = solver(b);
    solver = approxchol_sddm(sddm; tol=1e-6, maxits=1000, maxtime=Inf, verbose=false, pcgIts=Int[], params=ApproxCholParams())

Solves sddm systems by wrapping approxchol_lap.
Not yet optimized directly for sddm.

For more info, see http://danspielman.github.io/Laplacians.jl/latest/usingSolvers/index.html
"""
approxchol_sddm = sddmWrapLap(approxchol_lap)
=#
#===============================

 Compute only the LDLinv preconditioner

 ==============================#

 function approxchol_lap_pc(a::SparseMatrixCSC{Tv,Ti};
   tol::Real=1e-6,
   maxits=1000,
   maxtime=Inf,
   verbose=false,
   pcgIts=Int[],
   params=ApproxCholParams()) where {Tv,Ti}

   if minimum(a.nzval) < 0
       error("Adjacency matrix can not have negative edge weights")
   end

     return approxchol_lap1_pc(a,
     verbose=verbose,
     tol=tol,
     maxits=maxits,
     maxtime=maxtime,
     pcgIts=pcgIts,
     params=params)

 end

 function approxchol_lapGreedy_pc(a::SparseMatrixCSC;
   tol::Real=1e-6,
   maxits=1000,
   maxtime=Inf,
   verbose=false,
   pcgIts=Int[],
   params=ApproxCholParams())

   tol_ =tol
   maxits_ =maxits
   maxtime_ =maxtime
   verbose_ =verbose
   pcgIts_ =pcgIts

   t1 = time()

   la = lap(a) # a hit !?

   llmat = LLmatp(a)
   ldli = approxChol(llmat)
   return ldli

 end

 function approxchol_lapGiven_pc(a::SparseMatrixCSC;
   tol::Real=1e-6,
   maxits=1000,
   maxtime=Inf,
   verbose=false,
   pcgIts=Int[],
   params=ApproxCholParams())

   tol_ =tol
   maxits_ =maxits
   maxtime_ =maxtime
   verbose_ =verbose
   pcgIts_ =pcgIts

   t1 = time()

   la = lap(a)

   llmat = LLMatOrd(a)
   ldli = approxChol(llmat)

 end

 function approxchol_lapWdeg_pc(a::SparseMatrixCSC;
   tol::Real=1e-6,
   maxits=1000,
   maxtime=Inf,
   verbose=false,
   pcgIts=Int[],
   params=ApproxCholParams())

   tol_ =tol
   maxits_ =maxits
   maxtime_ =maxtime
   verbose_ =verbose
   pcgIts_ =pcgIts

   t1 = time()

   la = lap(a)

   v = vec(sum(a,dims=1))
   v = v .* (1 .+ rand(length(v)))
   p = sortperm(v)

   llmat = LLMatOrd(a,p)
   ldli = approxChol(llmat)

   ip = invperm(p)
   ldlip = LDLinv(p[ldli.col], ldli.colptr, p[ldli.rowval], ldli.fval, ldli.d[ip]);
   return ldlip

 end

 function approxchol_lap1_pc(a::SparseMatrixCSC{Tv,Ti};
   tol::Real=1e-6,
   maxits=1000,
   maxtime=Inf,
   verbose=false,
   pcgIts=Int[],
   params=ApproxCholParams()) where {Tv,Ti}

     tol_ =tol
     maxits_ =maxits
     maxtime_ =maxtime
     verbose_ =verbose
     pcgIts_ =pcgIts

     if params.order == :deg

       return approxchol_lapGreedy_pc(a,
         verbose=verbose,
         tol=tol,
         maxits=maxits,
         maxtime=maxtime,
         pcgIts=pcgIts,
         params=params)

     elseif params.order == :wdeg

       return approxchol_lapWdeg_pc(a,
         verbose=verbose,
         tol=tol,
         maxits=maxits,
         maxtime=maxtime,
         pcgIts=pcgIts,
         params=params)

     else
       return approxchol_lapGiven_pc(a,
         verbose=verbose,
         tol=tol,
         maxits=maxits,
         maxtime=maxtime,
         pcgIts=pcgIts,
         params=params)

     end

 end

#===============================

  Checking the condition number

=================================#

"""
    cn = condNumber(a, ldli; verbose=false)

Given an adjacency matrix a and an ldli computed by approxChol,
this computes the condition number.
"""
function condNumber(a, ldli; verbose=false)
  la = lap(a)

  # construct the square operator
  g = function(b)

    y = copy(b)

    #=
    mu = mean(y)
    @inbounds for i in eachindex(y)
        y[i] = y[i] - mu
    end
      =#

    @inbounds for i in 1:(length(ldli.d))
        if ldli.d[i] != 0
            y[i] /= (ldli.d[i])^(1/2)
        else
            y[i] = 0
        end
    end

    backward!(ldli, y)

    y = la * y

    forward!(ldli, y)

    @inbounds for i in 1:(length(ldli.d))
        if ldli.d[i] != 0
            y[i] /= (ldli.d[i])^(1/2)
        else
            y[i] = 0
        end
    end

    #=
    mu = mean(y)
    @inbounds for i in eachindex(y)
        y[i] = y[i] - mu
    end
    =#

    return y
  end

  gOp = SqLinOp(true,1.0,size(a,1),g)
  upper = eigs(gOp;nev=1,which=:LM,tol=1e-2)[1][1]

  g2(b) = upper*b - g(b)
  g2Op = SqLinOp(true,1.0,size(a,1),g2)
  lower = upper - eigs(g2Op;nev=2,which=:LM,tol=1e-2)[1][2]

  if verbose
      println("lower: ", lower, ", upper: ", upper);
  end

  return upper/lower

end

#===========================================

  Alternate solver approach

===========================================#

"""
    L = ldli2Chol(ldli)
This produces a matrix L so that L L^T approximate the original Laplacians.
It is not quite a Cholesky factor, because it is off by a perm
(and the all-1s vector orthogonality.
"""
function ldli2Chol(ldli)
    n = length(ldli.colptr)
    m = n + length(ldli.fval)
    li = zeros(Int,m)
    lj = zeros(Int,m)
    lv = zeros(Float64,m)
    lptr = 0

    dhi = zeros(n)
    for i in 1:n
        if ldli.d[i] == 0
            dhi[i] = 1.0
        else
            dhi[i] = sqrt(ldli.d[i])
        end
    end

    scales = ones(n)
    for ii in 1:(n-1)
        i = ldli.col[ii]
        j0 = ldli.colptr[ii]
        j1 = ldli.colptr[ii+1]-1
        scales[i] = prod(1.0 .- ldli.fval[j0:(j1-1)])
    end

    for ii in 1:(n-1)
        i = ldli.col[ii]
        j0 = ldli.colptr[ii]
        j1 = ldli.colptr[ii+1]-1
        scale = scales[i] / dhi[i]

        scj = 1
        for jj in j0:(j1-1)
            j = ldli.rowval[jj]
            f = ldli.fval[jj]

            lptr += 1
            li[lptr] = i
            lj[lptr] = j
            lv[lptr] = -f*scj/scale

            scj = scj*(1-f)
        end
        j = ldli.rowval[j1]

        lptr += 1
        li[lptr] = i
        lj[lptr] = j
        lv[lptr] = -dhi[i]

        lptr += 1
        li[lptr] = i
        lj[lptr] = i
        lv[lptr] = 1/scale

    end

    for i in 1:n
        if ldli.d[i] == 0
            lptr += 1
            li[lptr] = i
            lj[lptr] = i
            lv[lptr] = 1.0
        end
    end

    return sparse(li,lj,lv,n,n)
    #return li, lj, lv
end

function LDLsolver(L::SparseMatrixCSC, b::Array)
    y = x6 = L \ (L' \ b)
    return y .- mean(y)
end

"""
This variation of approxChol creates a cholesky factor to do the elimination.
It has not yet been optimized, and does not yet make the cholesky factor lower triangular
"""
function approxchol_lapChol(a::SparseMatrixCSC{Tv,Ti}; tol::Real=1e-6, maxits=1000, maxtime=Inf, verbose=false, pcgIts=Int[]) where {Tv,Ti}

    tol_ =tol
    maxits_ =maxits
    maxtime_ =maxtime
    verbose_ =verbose
    pcgIts_ =pcgIts

    t1 = time()
    llmat = LLmatp(a)

    ldli = approxChol(llmat)

    chL = ldli2Chol(ldli)

    if verbose
      println("Factorization time: ", time()-t1)
      println("Ratio of operator edges to original edges: ", 2 * length(ldli.fval) / nnz(a))
    end

    F(b) = LDLsolver(chL, b)

    la = lap(a)

    if verbose
        println("ratio of max to min diagonal of laplacian : ", maximum(diag(la))/minimum(diag(la)))
    end

    f(b;tol=tol_,maxits=maxits_, maxtime=maxtime_, verbose=verbose_, pcgIts=pcgIts_) = pcg(la, b .- mean(b), F, tol=tol, maxits=maxits, maxtime=maxtime, pcgIts=pcgIts, verbose=verbose)

    return f
end

#=============================================================

ApproxCholPQ
It only implements pop, increment key, and decrement key.
All nodes with degrees 1 through n appear in their own doubly-linked lists.
Nodes of higher degrees are bundled together.

=============================================================#

function keyMap(x, n)
    return x <= n ? x : n + div(x,n)
end

function ApproxCholPQ(a::Vector{Tind}) where Tind

    n = length(a)
    elems = Array{ApproxCholPQElem{Tind}}(undef, n)
    lists = zeros(Tind, 2*n+1)
    minlist = one(n)

    for i in 1:length(a)
        key = a[i]
        head = lists[key]

        if head > zero(Tind)
            elems[i] = ApproxCholPQElem{Tind}(zero(Tind), head, key)

            elems[head] = ApproxCholPQElem{Tind}(i, elems[head].next, elems[head].key)
        else
            elems[i] = ApproxCholPQElem{Tind}(zero(Tind), zero(Tind), key)

        end

        lists[key] = i
    end

    return ApproxCholPQ(elems, lists, minlist, n, n)
end

function approxCholPQPop!(pq::ApproxCholPQ{Tind}) where Tind
    if pq.nitems == 0
        error("ApproxPQ is empty")
    end
    while pq.lists[pq.minlist] == 0
        pq.minlist = pq.minlist + 1
    end
    i = pq.lists[pq.minlist]
    next = pq.elems[i].next

    pq.lists[pq.minlist] = next
    if next > 0
        pq.elems[next] = ApproxCholPQElem(zero(Tind), pq.elems[next].next, pq.elems[next].key)
    end

    pq.nitems -= 1

    return i
end

function approxCholPQMove!(pq::ApproxCholPQ{Tind}, i, newkey, oldlist, newlist) where Tind

    prev = pq.elems[i].prev
    next = pq.elems[i].next

    # remove i from its old list
    if next > zero(Tind)
        pq.elems[next] = ApproxCholPQElem{Tind}(prev, pq.elems[next].next, pq.elems[next].key)
    end
    if prev > zero(Tind)
        pq.elems[prev] = ApproxCholPQElem{Tind}(pq.elems[prev].prev, next, pq.elems[prev].key)

    else
        pq.lists[oldlist] = next
    end

    # insert i into its new list
    head = pq.lists[newlist]
    if head > 0
        pq.elems[head] = ApproxCholPQElem{Tind}(i, pq.elems[head].next, pq.elems[head].key)
    end
    pq.lists[newlist] = i

    pq.elems[i] = ApproxCholPQElem{Tind}(zero(Tind), head, newkey)

    return nothing
end

"""
    Decrement the key of element i
    This could crash if i exceeds the maxkey
"""
function approxCholPQDec!(pq::ApproxCholPQ{Tind}, i) where Tind

    oldlist = keyMap(pq.elems[i].key, pq.n)
    newlist = keyMap(pq.elems[i].key - one(Tind), pq.n)

    if newlist != oldlist

        approxCholPQMove!(pq, i, pq.elems[i].key - one(Tind), oldlist, newlist)

        if newlist < pq.minlist
            pq.minlist = newlist
        end

    else
        pq.elems[i] = ApproxCholPQElem{Tind}(pq.elems[i].prev, pq.elems[i].next, pq.elems[i].key - one(Tind))
    end

    return nothing
end

"""
    Increment the key of element i
    This could crash if i exceeds the maxkey
"""
function approxCholPQInc!(pq::ApproxCholPQ{Tind}, i) where Tind

    oldlist = keyMap(pq.elems[i].key, pq.n)
    newlist = keyMap(pq.elems[i].key + one(Tind), pq.n)

    if newlist != oldlist

        approxCholPQMove!(pq, i, pq.elems[i].key + one(Tind), oldlist, newlist)

    else
        pq.elems[i] = ApproxCholPQElem{Tind}(pq.elems[i].prev, pq.elems[i].next, pq.elems[i].key + one(Tind))
    end

    return nothing
end


text:
############################################################
# joAbstractFosterLinearOperatorInplace subtypes ###########
############################################################

############################################################
# joMatrixInplace ##########################################

export joMatrixInplace, joMatrixInplaceException

# type definition
"""
joMatrixInplace type

# TYPE PARAMETERS
- DDT::DataType : domain DataType
- RDT::DataType : range DataType

# FIELDS
- name::String : given name
- m::Integer : # of rows
- n::Integer : # of columns
- fop::Function : forward matrix
- fop_T::Function : transpose matrix
- fop_A::Function : adjoint matrix
- fop_C::Function : conj matrix
- iop::Nullable{Function} : inverse for fop
- iop_T::Nullable{Function} : inverse for fop_T
- iop_A::Nullable{Function} : inverse for fop_A
- iop_C::Nullable{Function} : inverse for fop_C

"""
struct joMatrixInplace{DDT<:Number,RDT<:Number} <: joAbstractLinearOperatorInplace{DDT,RDT}
    name::String
    m::Integer
    n::Integer
    fop::Function    # forward
    fop_T::Function  # transpose
    fop_A::Function  # adjoint
    fop_C::Function  # conj
    iop::Nullable{Function}
    iop_T::Nullable{Function}
    iop_A::Nullable{Function}
    iop_C::Nullable{Function}
end

# type exception
struct joMatrixInplaceException <: Exception
    msg :: String
end

############################################################
# joLooseMatrixInplace #####################################

export joLooseMatrixInplace, joLooseMatrixInplaceException

# type definition
"""
joLooseMatrixInplace type

# TYPE PARAMETERS
- DDT::DataType : domain DataType
- RDT::DataType : range DataType

# FIELDS
- name::String : given name
- m::Integer : # of rows
- n::Integer : # of columns
- fop::Function : forward matrix
- fop_T::Function : transpose matrix
- fop_A::Function : adjoint matrix
- iop::Nullable{Function} : inverse for fop
- iop_T::Nullable{Function} : inverse for fop_T
- iop_A::Nullable{Function} : inverse for fop_A

"""
struct joLooseMatrixInplace{DDT<:Number,RDT<:Number} <: joAbstractLinearOperatorInplace{DDT,RDT}
    name::String
    m::Integer
    n::Integer
    fop::Function    # forward
    fop_T::Function  # transpose
    fop_A::Function  # adjoint
    fop_C::Function  # conj
    iop::Nullable{Function}
    iop_T::Nullable{Function}
    iop_A::Nullable{Function}
    iop_C::Nullable{Function}
end

# type exception
struct joLooseMatrixInplaceException <: Exception
    msg :: String
end

############################################################
# joLinearFunctionInplace ##################################

export joLinearFunctionInplace, joLinearFunctionInplaceException

# type definition
"""
joLinearFunctionInplace type

# TYPE PARAMETERS
- DDT::DataType : domain DataType
- RDT::DataType : range DataType

# FIELDS
- name::String : given name
- m::Integer : # of rows
- n::Integer : # of columns
- fop::Function : forward function
- fop_T::Nullable{Function} : transpose function
- fop_A::Nullable{Function} : adjoint function
- iop::Nullable{Function} : inverse for fop
- iop_T::Nullable{Function} : inverse for fop_T
- iop_A::Nullable{Function} : inverse for fop_A

"""
struct joLinearFunctionInplace{DDT<:Number,RDT<:Number} <: joAbstractLinearOperatorInplace{DDT,RDT}
    name::String
    m::Integer
    n::Integer
    fop::Function              # forward
    fop_T::Nullable{Function}  # transpose
    fop_A::Nullable{Function}  # adjoint
    fop_C::Nullable{Function}  # conj
    iop::Nullable{Function}
    iop_T::Nullable{Function}
    iop_A::Nullable{Function}
    iop_C::Nullable{Function}
end

# type exceptions
struct joLinearFunctionInplaceException <: Exception
    msg :: String
end

############################################################
# joLooseLinearFunctionInplace #############################

export joLooseLinearFunctionInplace, joLooseLinearFunctionInplaceException

# type definition
"""
joLooseLinearFunctionInplace type

# TYPE PARAMETERS
- DDT::DataType : domain DataType
- RDT::DataType : range DataType

# FIELDS
- name::String : given name
- m::Integer : # of rows
- n::Integer : # of columns
- fop::Function : forward function
- fop_T::Nullable{Function} : transpose function
- fop_A::Nullable{Function} : adjoint function
- iop::Nullable{Function} : inverse for fop
- iop_T::Nullable{Function} : inverse for fop_T
- iop_A::Nullable{Function} : inverse for fop_A

"""
struct joLooseLinearFunctionInplace{DDT<:Number,RDT<:Number} <: joAbstractLinearOperatorInplace{DDT,RDT}
    name::String
    m::Integer
    n::Integer
    fop::Function              # forward
    fop_T::Nullable{Function}  # transpose
    fop_A::Nullable{Function}  # adjoint
    fop_C::Nullable{Function}  # conj
    iop::Nullable{Function}
    iop_T::Nullable{Function}
    iop_A::Nullable{Function}
    iop_C::Nullable{Function}
end

# type exceptions
struct joLooseLinearFunctionInplaceException <: Exception
    msg :: String
end


text:
GenCoNumbers = 5
genCoHorizon = 20
totalHorizon = 50
day2hour = 24
year2day = 6
milestoneYear = 5
interestRate = 0.05
segmentNumbers = 1
downpaymentRatio = 0.0

priceCap = 3000

###Technology parameters
technologyType = [:Base,:Mid,:Peak]
emission_ton_perMWh = [0,0,0]
lifeTime = [20,20,20]
domesticLearningRates = [0,0,0]
nonDomesticLearningRates = [0,0,0]

availability = [1,1,1]
constructionTimeStamp = [0,0,0]
constructionTime = [0,0,0]

deconstructionTime = [1,1,1]
efficiency = [0.4,0.48,0.6] #Conversion fuel -> electricity 0.0-1.0
fuelPrice = [3,15,25]

unitCapacity = [100,100,100] #MW
minCapacity = [0,0,0]

#EconomicProperty
technologyVOM = [5,4,4] #EUR/MWh
technologyFOM_perMWa = [80,40,17] * 1000 #EUR/MWa
constructionCost_perMW = [3000,1200,800] * 1000 #EUR/MW
deconstructionCost_perMW = constructionCost_perMW .* [0.1,0.1,0.1] #EUR/MW

#Technological cost of each unit
unitVOMCost = technologyVOM
unitFOMCost = technologyFOM_perMWa .* unitCapacity
unitConstructionCost = constructionCost_perMW .* unitCapacity
unitDeconstructionCost = deconstructionCost_perMW .* unitCapacity

emissionPrice = 0.0


text:
using SpikingNetworks

parseParameters"""
N=10000
Ne=N÷5*4;                 Ni=N÷5;
re=rand(Ne);            ri=rand(Ni);
a=[0.02*ones(Ne);       0.02+0.08*ri];
b=[0.2*ones(Ne);        0.25-0.05*ri];"""

parseEquations"""
dv/dt=0.04*v^2+5*v+140-u+I(t)
du/dt=a*(b*v-u)
I(t)=[5*randn(Ne);2*randn(Ni)]"""

net=init"""
N=N
v=-65*ones(N)
u=b.*v
c=[-65+15*re.^2;        -65*ones(Ni)]
d=[8-6*re.^2;           2*ones(Ni)]
S=[0.5*sprand(N,Ne,0.02) sprand(N,Ni,0.02)]"""

spike"v.>30"
reset"""
v[i]=c[i]
u[i]+=d[i]"""
synapse"v+=S"

net.dt=1
@time run!(net,1000/net.dt)
using MatlabPlot
t,neuron=rasterPlot(net)
mplot(t,neuron,".k")


text:
# test peakassign
known_energies = [3691, 4090, 4510, 4952, 5414, 5898, 6403, 6929, 7477, 8040, 8637] # kalpha energies Ca-Zn

filt_values = known_energies.^0.9 # caluclate filt values with some amount of non-linearity

tested = 0
peakassignwrong=0

nc = 7
tested+=length(combinations(1:length(known_energies),nc))
for combo in combinations(1:length(known_energies),nc)
	locations_arb = filt_values
	locations_true = known_energies[combo]
	best_combo, best_fom = peakassign(locations_arb, locations_true)
	if !all(best_combo .== filt_values[combo])
		peakassignwrong+=1
	end
end

@test peakassignwrong == 0

# test findpeaks
energy_value_distribution = MixtureModel([Normal(mu,15) for mu in known_energies])

be = 0:1:20000 # bin edges
events = rand(energy_value_distribution, 10000)
be, c = hist(events, be)

# plotpeaks(be, c, 10,10) # dont plot in regular testing
peakinds = findpeaks(midpoints(be), c, fwhm=50,n=length(known_energies))
@test length(peakinds)==length(known_energies)
peakx = sort(midpoints(be)[peakinds])
@test all(abs(peakx-known_energies).<5);

# test findpeaks on Histogram
bin_edges_v = 0:1:20000
h = Histogram(bin_edges_v)
addcounts!(h, events)
peakinds = findpeaks(h,fwhm=50,n=length(known_energies))
@test length(peakinds)==length(known_energies)
peakx = sort(midpoints(be)[peakinds])
@test all(abs(peakx-known_energies).<5);

calibrate_nofit(h, known_energies)


text:
###################
# DISPLAY METHODS #
###################

# mathjax = "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.7/MathJax.js?config=TeX-AMS-MML_SVG"
#COLORS
function theme_chooser(darkmode)
	if darkmode
		bgcolor = "rgba(0,0,0,0)"#"rgb(13,16,17)"
		text_color = "gray"
		plot_bgcolor = "rgb(22,26,29)" #"rgb(33,37,41)"
		grid_color = "rgb(40,52,66)" #rgb(40,40,40)
		sep_color = "white"
	else
		bgcolor = "rgba(0,0,0,0)"#"white"
		text_color = "gray"#"rgb(49,70,101)"
		plot_bgcolor = "rgb(229,236,246)"
		grid_color = "white"
		sep_color = "white"
	end
	bgcolor, text_color, plot_bgcolor, grid_color, sep_color
end
function interp_map(c_map, t_interp)
	idx = [c[1] for c = c_map]
	R = [parse.(Int, split(c[2][5:end-1],", "))[1] for c = c_map]
	G = [parse.(Int, split(c[2][5:end-1],", "))[2] for c = c_map]
	B = [parse.(Int, split(c[2][5:end-1],", "))[3] for c = c_map]
	r = LinearInterpolation(idx,R)(t_interp)
	g = LinearInterpolation(idx,G)(t_interp)
	b = LinearInterpolation(idx,B)(t_interp)
	c_map_interp = ["hsv($r, $g, $b)" for (r,g,b)=zip(r,g,b)]
	c_map_interp
end

plot_seq(seq::Sequence; width=nothing, height=nothing, slider=true, show_seq_blocks=false, darkmode=false, max_rf_samples=100) = begin
	bgcolor, text_color, plot_bgcolor, grid_color, sep_color = theme_chooser(darkmode)
	idx = ["Gx" "Gy" "Gz"]
	N = length(seq)
	O = size(seq.RF,1)
	ΔT = durs(seq)
	T0 = cumsum([0; ΔT],dims=1)
	off_val = Inf #This removes the unnecessary points in the plot
	#GRADS
	t1x = vcat([get_theo_t(seq.GR[1,i]) .+ T0[i] for i=1:N]...)
	t1y = vcat([get_theo_t(seq.GR[2,i]) .+ T0[i] for i=1:N]...)
	t1z = vcat([get_theo_t(seq.GR[3,i]) .+ T0[i] for i=1:N]...)
	Gx =  vcat([get_theo_A(seq.GR[1,i];off_val) for i=1:N]...)
	Gy =  vcat([get_theo_A(seq.GR[2,i];off_val) for i=1:N]...)
	Gz =  vcat([get_theo_A(seq.GR[3,i];off_val) for i=1:N]...)
	#RFS
	t2 =  vcat([get_theo_t(seq.RF[1,i];max_rf_samples) .+ T0[i] for i=1:N]...)
	R =   vcat([get_theo_A(r;off_val,max_rf_samples) for r = seq.RF]...)
	#ADC
	t3 =  vcat([get_theo_t(seq.ADC[i])  .+ T0[i] for i=1:N]...)
	D =   vcat([get_theo_A(d;off_val) for d = seq.ADC]...)
	#Shapes
	shapes = []
	if show_seq_blocks
		aux = [line(
			xref="x", yref="paper",
			x0=T0[i]*1e3, y0=0,
			x1=T0[i]*1e3, y1=1,
			line=attr(color=sep_color, width=1),
			) for i = 1:N+1]
		append!(shapes, aux)
	end
	# The commented code here is to visually check the simulation blocks
	# if show_sim_blocks
	# 	aux = [line(
	# 		xref="x", yref="paper",
	# 		x0=t_sim_parts[i]*1e3, y0=0,
	# 		x1=t_sim_parts[i]*1e3, y1=1,
	# 		line=attr(color="Red", width=1),
	# 		) for i = 1:length(t_sim_parts)]
	# 		append!(shapes, aux)
	# end
	l = PlotlyJS.Layout(; hovermode="closest",
			xaxis_title="",
			modebar=attr(orientation="h",yanchor="bottom",xanchor="right",y=1,x=0,bgcolor=bgcolor,color=text_color,activecolor=plot_bgcolor),
			legend=attr(orientation="h",yanchor="bottom",xanchor="left",y=1,x=0),
			plot_bgcolor=plot_bgcolor,
			paper_bgcolor=bgcolor,
			xaxis_gridcolor=grid_color,
			yaxis_gridcolor=grid_color,
			xaxis_zerolinecolor=grid_color,
			yaxis_zerolinecolor=grid_color,
			font_color=text_color,
			yaxis_fixedrange = false,
			xaxis=attr(
				ticksuffix=" ms",
				range=[0.,min(20,1e3*dur(seq))],
				rangeslider=attr(visible=slider),
				rangeselector=attr(
					buttons=[
						attr(count=1,
						label="1m",
						step=10,
						stepmode="backward"),
						attr(step="all")
						]
					),
				),
			shapes = shapes,
			margin=attr(t=0,l=0,r=0,b=0),
			height=height, width=width
			)
	plotter = PlotlyJS.scatter #using scattergl speeds up the plotting but does not show the sequence in the slider below
	p = [plotter() for j=1:(3+2O+1)]
	#GR
	p[1] = plotter(x=t1x*1e3, y=Gx*1e3,name=idx[1],hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m)")
	p[2] = plotter(x=t1y*1e3, y=Gy*1e3,name=idx[2],hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m)")
	p[3] = plotter(x=t1z*1e3, y=Gz*1e3,name=idx[3],hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m)")
	#RF
	for j=1:O
		phase = angle.(R[:,j])
		phase[R[:,j] .== Inf] .= Inf
		p[2j-1+3] = plotter(x=t2*1e3, y=abs.(R[:,j])*1e6,name="|RF_$j|",hovertemplate="(%{x:.4f} ms, %{y:.2f} μT)")
		p[2j+3] =   plotter(x=t2*1e3, y=phase, text=ones(size(t2)), name="<RF_$j",hovertemplate="(%{x:.4f} ms, ∠B1: %{y:.4f} rad)", visible="legendonly")
	end
	#ADC
	p[2O+3+1] = plotter(x=t3*1e3, y=D*1., name="ADC",hovertemplate="(%{x:.4f} ms, %{y:i})")
	config = PlotConfig(
		displaylogo=false,
		toImageButtonOptions=attr(
			format="svg", # one of png, svg, jpeg, webp
		).fields,
		modeBarButtonsToRemove=["zoom", "autoScale", "resetScale2d", "pan", "tableRotation", "resetCameraLastSave", "zoomIn", "zoomOut"]
	)
	PlotlyJS.plot(p, l; config) #, options=Dict(:displayModeBar => false))
end

function plot_image(image; height=750, width=nothing, zmin=minimum(abs.(image[:])), zmax=maximum(abs.(image[:])), darkmode=false, title="")
	#Layout
	bgcolor, text_color, plot_bgcolor, grid_color, sep_color = theme_chooser(darkmode)
	l = PlotlyJS.Layout(;title=title,yaxis_title="y",
    xaxis_title="x",height=height,wifth=width,
    yaxis=attr(scaleanchor="x"),
	font_color=text_color,
    modebar=attr(orientation="v",bgcolor=bgcolor,color=text_color,activecolor=plot_bgcolor),xaxis=attr(constrain="domain"),hovermode="closest",
	paper_bgcolor=bgcolor,
	plot_bgcolor=plot_bgcolor,
	xaxis_gridcolor=grid_color,
	yaxis_gridcolor=grid_color,
	xaxis_zerolinecolor=grid_color,
	yaxis_zerolinecolor=grid_color)
	#Plot
	p = PlotlyJS.heatmap(z=image,transpose=false,zmin=zmin,zmax=zmax,colorscale="Greys")
	config = PlotConfig(
		displaylogo=false,
		toImageButtonOptions=attr(
			format="svg", # one of png, svg, jpeg, webp
		).fields,
		modeBarButtonsToRemove=["zoom", "autoScale", "resetScale2d", "pan", "tableRotation", "resetCameraLastSave", "zoomIn", "zoomOut"]
	)
	PlotlyJS.plot(p,l;config)
end

function plot_kspace(seq; width=nothing, height=nothing, darkmode=false)
	bgcolor, text_color, plot_bgcolor, grid_color, sep_color = theme_chooser(darkmode)
	#Calculations of theoretical k-space
	kspace, kspace_adc = get_kspace(seq; Δt=1) #simParams["Δt"])
	t_adc = get_sample_times(seq)
	#Colormap
	c_map = [[t, "hsv($(floor(Int,(1-t)*255)), 100, 50)"] for t=range(0,1,10)]
	c = "gray"
	c2_idx = []
	counter = 0
	for s in seq
		if is_ADC_on(s)
			N = s.ADC.N[1]
			append!(c2_idx, counter:N+counter-1)
			counter += N
		end
	end
	c2 = interp_map(c_map, c2_idx ./ maximum(c2_idx))
	#Layout
	mink = minimum(kspace_adc,dims=1)
	maxk = maximum(kspace_adc,dims=1)
	dW = maximum(maxk .- mink, dims=2) * .3
	mink .-= dW
	maxk .+= dW
	#Layout
	l = PlotlyJS.Layout(;
		paper_bgcolor=bgcolor,
		scene=attr(xaxis=attr(title="kx [m⁻¹]",range=[mink[1],maxk[1]],backgroundcolor=plot_bgcolor,gridcolor=grid_color,zerolinecolor=grid_color),
				   yaxis=attr(title="ky [m⁻¹]",range=[mink[2],maxk[2]],backgroundcolor=plot_bgcolor,gridcolor=grid_color,zerolinecolor=grid_color),
				   zaxis=attr(title="kz [m⁻¹]",range=[mink[3],maxk[3]],backgroundcolor=plot_bgcolor,gridcolor=grid_color,zerolinecolor=grid_color)),
		modebar=attr(orientation="h",yanchor="bottom",xanchor="right",y=1,x=0,bgcolor=bgcolor,color=text_color,activecolor=plot_bgcolor),
		legend=attr(orientation="h",yanchor="bottom",xanchor="left",y=1,x=0),
		height=height, width=width,font_color=text_color,
		scene_camera_eye=attr(x=0, y=0, z=1.7),
		scene_camera_up=attr(x=0, y=1., z=0),
		scene_aspectmode="cube",
		margin=attr(t=0,l=0,r=0))
	#Plot
	p = [PlotlyJS.scatter() for j=1:3]
	p[1] = PlotlyJS.scatter3d(x=kspace[:,1],y=kspace[:,2],z=kspace[:,3],mode="lines",
			line=attr(color=c),name="Trajectory",hoverinfo="skip")
	p[2] = PlotlyJS.scatter3d(x=kspace_adc[:,1],y=kspace_adc[:,2],z=kspace_adc[:,3],text=round.(t_adc*1e3,digits=3),mode="markers",
			line=attr(color=c2),marker=attr(size=2),name="ADC",hovertemplate="kx: %{x:.1f} m⁻¹
ky: %{y:.1f} m⁻¹
kz: %{z:.1f} m⁻¹
t_acq: %{text} ms<extra></extra>")
	p[3] = PlotlyJS.scatter3d(x=[0],y=[0],z=[0],name="k=0",marker=attr(symbol="cross",size=10,color="red"))
	config = PlotConfig(
		displaylogo=false,
		toImageButtonOptions=attr(
			format="svg", # one of png, svg, jpeg, webp
		).fields,
		modeBarButtonsToRemove=["zoom", "pan", "tableRotation", "resetCameraLastSave3d", "orbitRotation", "resetCameraDefault3d"]
	)
	PlotlyJS.plot(p,l; config)
end

function plot_M0(seq; height=nothing, width=nothing, slider=true, darkmode=false)
	#Times
	dt = 1
	t, Δt = Koma.get_uniform_times(seq, dt)
	#kx,ky
	ts = t .+ Δt
	rf_idx, rf_type = Koma.get_RF_types(seq, t)
	k, _ =  Koma.get_kspace(seq; Δt=dt)

	#plots k(t)
	bgcolor, text_color, plot_bgcolor, grid_color, sep_color = theme_chooser(darkmode)
	l = PlotlyJS.Layout(;yaxis_title="M0", hovermode="closest",
			xaxis_title="",
			modebar=attr(orientation="h",yanchor="bottom",xanchor="right",y=1,x=0,bgcolor=bgcolor,color=text_color,activecolor=plot_bgcolor),
			legend=attr(orientation="h",yanchor="bottom",xanchor="left",y=1,x=0),
			plot_bgcolor=plot_bgcolor,
			paper_bgcolor=bgcolor,
			xaxis_gridcolor=grid_color,
			yaxis_gridcolor=grid_color,
			xaxis_zerolinecolor=grid_color,
			yaxis_zerolinecolor=grid_color,
			font_color=text_color,
			yaxis_fixedrange = false,
			xaxis=attr(
				ticksuffix=" ms",
				range=[0.,min(20,1e3*dur(seq))],
				rangeslider=attr(visible=slider),
				rangeselector=attr(
					buttons=[
						attr(count=1,
						label="1m",
						step=10,
						stepmode="backward"),
						attr(step="all")
						]
					),
				),
			margin=attr(t=0,l=0,r=0),
			height=height, width=width
			)
	plotter = PlotlyJS.scatter #using scattergl speeds up the plotting but does not show the sequence in the slider below
	p = [plotter() for j=1:4]
	p[1] = plotter(x=ts*1e3, y=k[:,1], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms)", name="x")
	p[2] = plotter(x=ts*1e3, y=k[:,2], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms)", name="y")
	p[3] = plotter(x=ts*1e3, y=k[:,3], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms)", name="z")
	p[4] = plotter(x=t[rf_idx]*1e3,y=rf_type,name="RFs",marker=attr(symbol="cross",size=8,color="orange"),mode="markers")
	config = PlotConfig(
		displaylogo=false,
		toImageButtonOptions=attr(
			format="svg", # one of png, svg, jpeg, webp
		).fields,
		modeBarButtonsToRemove=["zoom", "select", "select2d","lasso2d", "autoScale", "resetScale2d", "pan", "tableRotation", "resetCameraLastSave", "zoomIn", "zoomOut"]
	)
	PlotlyJS.plot(p, l; config)
end

"""Plot phantom map, the keys could be {:ρ, :T1, :T2, :T2s, :x, :y, :z}. Optional parameter: t0 [ms] to see displacement.

Colormaps from https://github.com/markgriswold/MRFColormaps
Towards Unified Colormaps for Quantitative MRF Data, Mark Griswold, et al. (2018)."""
function plot_phantom_map(ph::Phantom, key::Symbol; t0=0, height=700, width=nothing, darkmode=false)
	cmin_key = minimum(getproperty(ph,key))
	cmax_key = maximum(getproperty(ph,key))
	if key == :T1 || key == :T2 || key == :T2s
		cmin_key = 0
		factor = 1e3
		unit = " ms"
		if key  == :T1
			cmax_key = 2500/factor
			colors = MAT.matread("src/ui/assets/T1cm.mat")["T1colormap"]
			N, _ = size(colors)
			idx = range(0,1,N)
			colormap = [[idx[n], "rgb($(floor(Int,colors[n,1]*255)),$(floor(Int,colors[n,2]*255)),$(floor(Int,colors[n,3]*255)))"] for n=1:N]
		elseif key == :T2 || key == :T2s
			if key == :T2
				cmax_key = 250/factor
			end
    		colors = MAT.matread("src/ui/assets/T2cm.mat")["T2colormap"]
			N, _ = size(colors)
			idx = range(0,1,N)
			colormap = [[idx[n], "rgb($(floor(Int,colors[n,1]*255)),$(floor(Int,colors[n,2]*255)),$(floor(Int,colors[n,3]*255)))"] for n=1:N]
		end
	elseif key == :x || key == :y || key == :z
		factor = 1e2
		unit = " cm"
		colormap="Greys"
	elseif key == :Δw
		factor = 1/(2π)
		unit = " Hz"
		colormap="Greys"
	else
		factor = 1
		cmin_key = 0
		unit=""
		colormap="Greys"
	end
	cmin_key *= factor
	cmax_key *= factor
	x0 = minimum([ph.x;ph.y;ph.z])*1e2
    xf = maximum([ph.x;ph.y;ph.z])*1e2
	#Layout
	bgcolor, text_color, plot_bgcolor, grid_color, sep_color = theme_chooser(darkmode)
	l = PlotlyJS.Layout(;title=ph.name, height=height, width=width,
		paper_bgcolor=bgcolor,font_color=text_color,
		scene=attr(
			xaxis=attr(title="x",range=[x0,xf],ticksuffix=" cm",backgroundcolor=plot_bgcolor,gridcolor=grid_color,zerolinecolor=grid_color),
			yaxis=attr(title="y",range=[x0,xf],ticksuffix=" cm",backgroundcolor=plot_bgcolor,gridcolor=grid_color,zerolinecolor=grid_color),
			zaxis=attr(title="z",range=[x0,xf],ticksuffix=" cm",backgroundcolor=plot_bgcolor,gridcolor=grid_color,zerolinecolor=grid_color),
		),
		modebar=attr(orientation="h",bgcolor=bgcolor,color=text_color,activecolor=plot_bgcolor),xaxis=attr(constrain="domain"),
		hovermode="closest",aspectmode="cube", margin=attr(t=25,l=0,r=0))
	h = PlotlyJS.scatter3d( x=(ph.x .+ ph.ux(ph.x,ph.y,ph.z,t0*1e-3))*1e2,
							y=(ph.y .+ ph.uy(ph.x,ph.y,ph.z,t0*1e-3))*1e2,
							z=(ph.z .+ ph.uz(ph.x,ph.y,ph.z,t0*1e-3))*1e2,
							mode="markers",
							marker=attr(color=getproperty(ph,key)*factor,
										showscale=true,
										colorscale=colormap,
										colorbar=attr(ticksuffix=unit, title=string(key)),
										cmin=cmin_key,
										cmax=cmax_key,
										size=2
										),
							text=round.(getproperty(ph,key)*factor,digits=4),
							hovertemplate="x: %{x:.1f} cm
y: %{y:.1f} cm
z: %{z:.1f} cm
$(string(key)): %{text}$unit<extra></extra>")
	config = PlotConfig(
		displaylogo=false,
		toImageButtonOptions=attr(
			format="svg", # one of png, svg, jpeg, webp
		).fields,
		modeBarButtonsToRemove=["zoom", "pan", "tableRotation", "resetCameraLastSave3d", "orbitRotation", "resetCameraDefault3d"]
	)
	p = PlotlyJS.plot(h,l;config)
end

function plot_signal(raw::RawAcquisitionData; height=nothing, width=nothing, darkmode=false)
	not_Koma = raw.params["systemVendor"] != "Koma.jl"
	t = []
	signal = []
	current_t0 = 0
	for p in raw.profiles
		dt = p.head.sample_time_us != 0 ? p.head.sample_time_us * 1e-3 : 1
		t0 = p.head.acquisition_time_stamp * 1e-3 #This parameter is used in Koma to store the time offset
		N =  p.head.number_of_samples != 0 ? p.head.number_of_samples : 1
		if not_Koma
			t0 = current_t0 * dt
			current_t0 += N
		end
		append!(t, t0.+(0:dt:dt*(N-1)))
		append!(signal, p.data[:,1]) #Just one coil
		#To generate gap
		append!(t, t[end])
		append!(signal, [Inf + Inf*1im])
	end
	#Time window
	Tmin = minimum(t)
	Tmax = maximum(t)
	Tacq = Tmax - Tmin
	#PLOT
	bgcolor, text_color, plot_bgcolor, grid_color, sep_color = theme_chooser(darkmode)
	l = PlotlyJS.Layout(;title="",
		font_color=text_color,
		xaxis=attr(ticksuffix=" ms",range=[Tmin+Tacq/2-.025*Tacq,Tmin+Tacq/2+.025*Tacq],
					rangeslider=attr(visible=true),
					gridcolor=grid_color,zerolinecolor=grid_color),
		yaxis=attr(gridcolor=grid_color,zerolinecolor=grid_color),
		paper_bgcolor=bgcolor,
		plot_bgcolor=plot_bgcolor, height=height, width=width,
		yaxis_fixedrange = false,
		legend=attr(orientation="h",yanchor="bottom",xanchor="left",y=1.02,x=0),
		modebar=attr(orientation="h",yanchor="bottom",xanchor="right",y=1.02,x=1,bgcolor=bgcolor,color=text_color,activecolor=plot_bgcolor))
	absS = PlotlyJS.scatter(x=t,y=abs.(signal), name="|S(t)|")
	reS =  PlotlyJS.scatter(x=t,y=real.(signal),name="Re{S(t)}")
	imS =  PlotlyJS.scatter(x=t,y=imag.(signal),name="Im{S(t)}")
	config = PlotConfig(
		displaylogo=false,
		toImageButtonOptions=attr(
			format="svg", # one of png, svg, jpeg, webp
		).fields,
		modeBarButtonsToRemove=["zoom", "select2d", "lasso2d", "autoScale", "resetScale2d", "pan", "tableRotation", "resetCameraLastSave", "zoomIn", "zoomOut"]
	)
	p = PlotlyJS.plot([absS,reS,imS],l;config)
end

function plot_dict(dict::Dict)
	html = """

			#
			Name
			Value

	"""
	i = 1
	for (key,val) = dict
		html *= """

				$i
				$(string(key))
				$(string(val))

		"""
		i += 1
	end
	html *= "

"
end

# """Plots gradient moments M0, M1 and M2 for specified axis. Does NOT consider RF pulses."""
# plot_grads_moments(seq::Sequence, idx::Int=1; title="", mode="quick") = begin
# 	if mode == "quick"
# 		plotter = PlotlyJS.scattergl
# 	else
# 		plotter = PlotlyJS.scatter
# 	end
# 	names = ["Gx", "Gy", "Gz"]
# 	M, N = size(seq.GR)
# 	G = [seq.GR[j,floor(Int,i/2)+1].A for i=0:2*N-1, j=1:M]
# 	T = [seq.GR[1,i].T for i=1:N]
# 	t = [sum(T[1:i]) for i=1:N]
# 	t = [t[floor(Int,i/2)+1] for i=0:2*N-1]
# 	t = [0; t[1:end-1]]

# 	M0, M1, M2 = get_M0_M1_M2(seq)
# 	M0t = [M0(t)[idx] for t=t][:]; M0max = maximum(abs.(M0t))*1.1e6
# 	M1t = [M1(t)[idx] for t=t][:]; M1max = maximum(abs.(M1t))*1.1e9
# 	M2t = [M2(t)[idx] for t=t][:]; M2max = maximum(abs.(M2t))*1.1e12
# 	Gmax = maximum(abs.(G))*1.1e3

# 	p = [plotter() for k = 1:4]

# 	l = PlotlyJS.Layout(;
# 	title=title,
# 	xaxis=attr(domain = [0, .75]),
# 	yaxis=attr(title=attr(text="G [mT/m]", standoff=0),
# 		side="left", tickfont=attr(color="#636efa"), titlefont=attr(color="#636efa"),range=[-Gmax, Gmax]),
# 	yaxis2=attr(title=attr(text="M0 [mT/m⋅ms]", standoff=1, position=1), showgrid=false,
# 		anchor="free", overlaying="y",position=0.75, side="right",
# 		tickfont=attr(color="#ef553b"), titlefont=attr(color="#ef553b"), range=[-M0max, M0max]),
# 	yaxis3=attr(title=attr(text="M1 [mT/m⋅ms²]", standoff=1), showgrid=false,
# 		anchor="free",overlaying="y",position=0.75+.25/3,side="right",
# 		tickfont=attr(color="#45d9b2"), titlefont=attr(color="#45d9b2"), range=[-M1max, M1max]),
# 	yaxis4=attr(title=attr(text="M2 [mT/m⋅ms³]", standoff=1), showgrid=false,
# 		anchor="free",overlaying="y",position=0.75+.5/3,side="right",
# 		tickfont=attr(color="#b373fa"), titlefont=attr(color="#b373fa"), range=[-M2max, M2max]),
# 	xaxis_title="t [ms]", height=300)
# 	p[1] = plotter(x=t*1e3, y=G[:,idx]*1e3, name=names[idx] , line_shape="hv")
# 	p[2] = plotter(x=t*1e3, y=M0t*1e6, name="M0", yaxis="y2",line=attr(dash="dash"))
# 	p[3] = plotter(x=t*1e3, y=M1t*1e9, name="M1", yaxis="y3",line=attr(dash="dash"))
# 	p[4] = plotter(x=t*1e3, y=M2t*1e12, name="M2", yaxis="y4",line=attr(dash="dash"))
# 	PlotlyJS.plot(p, l)
# end
# plot_Phantom(obj::Phantom,filename::String) = begin
# 	# Phantom
# 	p1 = heatmap(obj.x*1e2,obj.y*1e2,obj.ρ,aspect_ratio=:equal)#,clims=(0,maximum(T2[:])))
# 	xaxis!("x [cm]"); yaxis!("y [cm]")
# 	title!("Proton density")
# 	p2 = heatmap(obj.x*1e2,obj.y*1e2,obj.T2*1e3,aspect_ratio=:equal)#,clims=(0,maximum(T2[:])))
# 	xaxis!("x [cm]"); yaxis!("y [cm]")
# 	title!("T2 [ms]")
# 	p3 = heatmap(obj.x*1e2,obj.y*1e2,obj.Δw/2π,aspect_ratio=:equal)#,clims=(0,maximum(T2[:])))
# 	xaxis!("x [cm]"); yaxis!("y [cm]")
# 	title!("df [Hz]")
# 	P(i,j) = rotz(Dθ[i,j])[1:2,1:2]; D(i,j) = [obj.Dλ1[i,j] 0;0 obj.Dλ2[i,j]]
# 	nx = [1;0]; ny = [0;1]
# 	Dx = [nx'*P(i,j)'*D(i,j)*P(i,j)*nx for i=1:size(obj.Dλ1,1),j=1:size(obj.Dλ1,2)]
# 	Dy = [ny'*P(i,j)'*D(i,j)*P(i,j)*ny for i=1:size(obj.Dλ1,1),j=1:size(obj.Dλ1,2)]
# 	p4 = heatmap(obj.x*1e2,obj.y*1e2,Dx*1e12,aspect_ratio=:equal)#,clims=(0,maximum(T2[:])))
# 	xaxis!("x [cm]"); yaxis!("y [cm]")
# 	title!("Dx [um2/s]")
# 	p5 = heatmap(obj.x*1e2,obj.y*1e2,Dy*1e12,aspect_ratio=:equal)#,clims=(0,maximum(T2[:])))
# 	xaxis!("x [cm]"); yaxis!("y [cm]")
# 	title!("Dy [um2/s]")
# 	p = plot(p1,p2,p3,p4,p5,size=(1300,500),layout=@layout [a b c; d e])
# 	savefig(p,filename)
# end
# plot_sim_res(obj::Phantom,SEQ::Array{Sequence},S::Array{ComplexF64},
# 	t::Array{Float64,2},filename::String,t_k0::Float64,Ga::Float64) = begin
# 	T2m = sum(obj.T2.*obj.ρ)/sum(obj.ρ)
# 	b, n = get_bvalue(SEQ[1])
# 	println("Sequence with b-value: "*string(round(b*1e-6,digits=2))*" s/mm2") # s/mm2
# 	P(i) = rotz(obj.Dθ[i])[1:2,1:2]; D(i) = [obj.Dλ1[i] 0;0 obj.Dλ2[i]]
# 	Deq = [n'*P(i)'*D(i)*P(i)*n for i=1:prod(size(obj.Dλ1))]
# 	Dm = sum(Deq.*obj.ρ)/sum(obj.ρ) #MODIFYY

# 	p = plot_grads(SEQ,t,t_k0,Ga)
# 	p3 = plot([t_k0; t_k0]*1e3,[0; 1.1],linewidth=0.25,color=:black,label="k=0")
# 	plot!(t[:]*1e3,abs.(S),label="|S|")
# 	plot!(t[:]*1e3,exp.(-t[:]/T2m),linestyle=:dash,label="exp(-t/T2)")
# 	# exp(-b*D) <-> exp(-4*π^2*(Δ-δ/3)*q'*D*q)
# 	plot!(t[:]*1e3,exp.(-t[:]/T2m.-b*Dm).*(t[:].>=(Δ+δ)),linestyle=:dash,label="exp(-t/T2-bD)")

# 	xlabel!("Time [ms]")
# 	ylabel!("Signal [a.u]")
# 	xlims!((minimum(t),minimum(t).+dur(sum(SEQ))).*1e3)
# 	p = plot(p[1],p[2],p3,size=(800,600),layout=@layout [a ; b; c])
# 	savefig(p,filename)
# end
# plot_ksapce_trajectory(ACQ::Sequence,t::Array{Float64,2},filename::String) = begin
# 	k = get_designed_kspace(ACQ)
# 	p = plot(k[:,1],k[:,2],legend=:none,aspect_ratio=:equal)
# 	k = get_actual_kspace(ACQ,t)
# 	scatter!(k[:,1],k[:,2],legend=:none,markersize=1)
# 	xlabel!("kx [1/m]"); ylabel!("ky [1/m]")
# 	savefig(p,filename)
# end
# plot_recon(kdata::Array{ComplexF64},rec::Array{ComplexF64},
# 	Δx_pix::Float64,Δy_pix::Float64,filename::String,title::String) = begin
# 	Nx, Ny = size(rec)
# 	xr = -Δx_pix*(Nx-1)/2:Δx_pix:Δx_pix*(Nx-1)/2
# 	yr = -Δy_pix*(Ny-1)/2:Δy_pix:Δy_pix*(Ny-1)/2
# 	Wx, Wy = 1/Δx_pix, 1/Δy_pix
# 	kx = range(-Wx/2,stop=Wx/2,length=Nx)
# 	ky = range(-Wy/2,stop=Wy/2,length=Ny)
# 	p1 = heatmap(kx,ky,log.(abs.(kdata).+1),aspect_ratio=:equal,
# 				legend=:none,size=(400,400))
# 	xaxis!("kx [1/m]"); yaxis!("ky [1/m]")
# 	title!("k-space")
# 	p2 = heatmap(xr*1e2,yr*1e2,abs.(rec),aspect_ratio=:equal,
# 				legend=:none,size=(400,400))
# 	xaxis!("x [cm]"); yaxis!("y [cm]")
# 	title!(title)
# 	#savefig(p1,filename*"_ksp.pdf")
# 	savefig(p2,filename*".pdf")
# end
# plot_Eq(vx,vy,S,S0,SEQ) = begin
# 	Nq, Nθ = size(SEQ)
# 	rs = range(0,stop=1,length=Nq)
# 	θs = range(0,stop=π,length=Nθ); θs = [θs[1:end-1]; π.+θs[1:end-1]];
#     Eq(vx,vy) = begin
#         Eq = [S[i,j][vx,vy] for i=1:Nq,j=1:Nθ]
#         Eq = [Eq[:,1:end-1] Eq[:,1:end-1]]
#     end
#     c = RGB{Float64}(1,1,1)
#     α = real.(abs.(S0[vx,vy]))/maximum(real.(abs.(S0)))
#     pyplot()
#     hm = heatmap(θs,rs,Eq(vx,vy)*(α<.2 ? 0 : 1),proj=:polar,aspect_ratio=:equal,
#     	legend=:none,grid=false,xticks=:none,yticks=:none,
#     	background_color_subplot=α*c)
# end
# plot_Pr(vx,vy,S,S0,SEQ) = begin
# 	Nq, Nθ = size(SEQ)
# 	rs = range(0,stop=1,length=10)*25e-6
# 	θs = range(0,stop=2π,length=10)
#     Pr(vx,vy) = begin
#         Eq = [S[i,j][vx,vy] for i=1:Nq,j=1:Nθ]
# 		b = [get_bvalue(SEQ[i,j])[1] for i=1:Nq,j=1:Nθ]
# 		n = [get_bvalue(SEQ[i,j])[2] for i=1:Nq,j=1:Nθ]
# 		B = [-b[i,j]*[n[i,j][1]^2
# 					2*n[i,j][1]*n[i,j][2]
# 					  n[i,j][2]^2] for i=1:Nq,j=1:Nθ]
# 		B = [B[i][j] for i=1:Nq*Nθ,j=1:3] #reshape
# 		y = log.(abs.(Eq))[:]
# 		# exp(-b*n'D n) <-> exp(-4*π^2*(Δ-δ/3)*q'*D*q)
# 		# min D || -b*n'D*n - log(S) ||_2 + λ ||D||_2
# 		# min D || B D - y ||_2 + λ ||D||_2 => D = (B'B+λI)^-1 B' y
# 		λ = 0
# 		Id = Matrix{Float64}(I,3,3); D = [0;0;0]
# 		for n = 0:2 #Tikhonov regularized and iterativly weighted
# 			W = n==0 ? diagm(0=>Eq[:].^2) : diagm(0=>exp.(2*B*D))
# 			D = (B'*W*B + λ*Id)^-1*B'*W*y
# 		end
# 		# Diffusion propagator
# 		D_inv = [D[1] D[2];
# 				 D[2] D[3]]^-1
# 		pr = exp.(-[([r*cos(θ) r*sin(θ)]*D_inv*[r*cos(θ);r*sin(θ)])[1] for r=rs, θ=θs])
#     end
#     c = RGB{Float64}(1,1,1)
#     α = real.(abs.(S0[vx,vy]))/maximum(real.(abs.(S0)))
#     pyplot()
#     hm = heatmap(θs,rs,Pr(vx,vy)*(α<.2 ? 0 : 1),proj=:polar,aspect_ratio=:equal,
#     	legend=:none,grid=false,xticks=:none,yticks=:none,
#     	background_color_subplot=α*c)
# end


text:
using TensorLabXD
using MPSLabXD
using Test
using TestExtras

const TX = TensorLabXD

time_initial = time()
include("mps_test.jl")
time_final = time()
printstyled("Finished all tests in ",
            string(round((time_final-time_initial)/60; sigdigits=3)),
            " minutes."; bold = true, color = Base.info_color())
println()


text:
using FighterJets
using Test

@testset "FighterJets.jl" begin
    # Write your tests here.
end


text:
HELP_TEXT = "show this help message and exit"

struct ArgumentParser
    prog::Optional{String}
    usage::Optional{String}
    description::Optional{String}
    epilog::Optional{String}
    positional_args::Vector{Argument}
    optional_args::Vector{Argument}
    flag_args::LittleDict{String,Argument}
    has_numeric_flags::Ref{Bool}
end

function ArgumentParser(;
    prog = nothing,
    usage = nothing,
    description = nothing,
    epilog = nothing,
    add_help = true)

    positional_args = Argument[]
    optional_args = Argument[]
    flag_args = LittleDict{String,Argument}()
    has_numeric_flags = Ref{Bool}(false)

    parser = ArgumentParser(prog, usage, description, epilog, positional_args, optional_args, flag_args, has_numeric_flags)
    add_help && add_argument!(parser, "-h", "--help", action = "help", help = HELP_TEXT)

    return parser
end

function argument_adder(parser::ArgumentParser)
    function add_argument!(name_or_flags::String...; kwargs...)
        ArgParse2.add_argument!(parser, name_or_flags...; kwargs...)
    end
end

function add_argument(args...; kwargs...)
    @warn """add_argument(...) is deprecated and will be removed before ArgParse2 is registered.
             Please use add_argument!(...)""" maxlog = 1
    add_argument!(args...; kwargs...)
end

function add_argument!(parser::ArgumentParser, name_or_flags::String...; kwargs...)
    arg = Argument(name_or_flags...; kwargs...)

    if is_positional(arg)
        push!(parser.positional_args, arg)
    else
        for flag in arg.flags
            flag in keys(parser.flag_args) && throw(ArgumentError("flag `$flag` defined multiple times"))
            parser.flag_args[flag] = arg
            if is_numeric(flag)
                parser.has_numeric_flags[] = true
            end
        end
        push!(parser.optional_args, arg)
    end

    nothing
end

is_positional(arg::Argument) = isempty(arg.flags)


text:
macro mustimplement(sig)
    :($(esc(sig)) = error("must implement ", $(string(sig))))
end

@inline propertytype(x, name::Symbol) = typeof(getproperty(x, name))

argerror(s::AbstractString) = throw(ArgumentError(s))

function genpath(path::AbstractString; force::Bool = false, sep = '_')
    if force
        rm(path, force = true, recursive = true)
        return path
    elseif ispath(path)
        parts = splitpath(path)
        root = parts[1:length(parts)-1]
        base = parts[end]
        if isfile(path)
            filename, ext = splitext(base)
            f = i -> joinpath(root..., "$(filename)$(sep)$(i)$(ext)")
        elseif isdir(path)
            f = i -> joinpath(root..., "$(base)$(sep)$(i)")
        else
            error("Not a file or directory: $path")
        end
        i = 1
        while ispath(f(i))
            i += 1
        end
        return f(i)
    else
        return path
    end
end


text:
"""
    cxeltype(a)

Returns Complex{eltype(a)} if eltype(a) <: Real; eltype(a) otherwise.
"""
cxeltype(a) = eltype(a) <: Real ? Complex{eltype(a)} : eltype(a)

const steppers = [
  "ForwardEuler",
  "FilteredForwardEuler",
  "AB3",
  "FilteredAB3",
  "RK4",
  "FilteredRK4",
  "ETDRK4",
  "FilteredETDRK4",
]

const filteredsteppers = [
  "FilteredForwardEuler",
  "FilteredAB3",
  "FilteredRK4",
  "FilteredETDRK4",
]

"Returns a time-stepper of type `steppernameTimeStepper'."
TimeStepper(steppername, args...) = autoconstructtimestepper(steppername, args...)

"""
    autoconstructtimestepper(stepper, dt, sol, g=ZeroDGrid())

Returns a time-stepper type defined by the prefix `stepper`, timestep `dt`
solution `sol` (used to construct variables
with identical type and size as the solution vector), and grid `g`.
"""
function autoconstructtimestepper(stepper, dt, sol, g::AbstractGrid=ZeroDGrid())
  fullsteppername = Symbol(stepper, :TimeStepper)
  tsexpr = occursin("Filtered", stepper) ?
      Expr(:call, fullsteppername, dt, sol, g) : Expr(:call, fullsteppername, dt, sol)
  eval(tsexpr)
end

"""
    autoconstructtimestepper(stepper, dt, solc, solr, g=ZeroDGrid())

Returns a time-stepper type defined by the prefix `stepper`, timestep `dt`
complex solution `solc` and real solution `solr` (used to construct variables
with identical type and size as the solution vector), and grid `g`.

"""
function autoconstructtimestepper(stepper, dt, solc, solr, g::AbstractGrid=ZeroDGrid())
  fullsteppername = Symbol(stepper, :TimeStepper)
  tsexpr = occursin("Filtered", stepper) ?
    Expr(:call, fullsteppername, dt, solc, solr, g) : Expr(:call, fullsteppername, dt, solc, solr)
  eval(tsexpr)
end

"""
    @createarrays T dims a b c...

Create arrays of all zeros with element type `T`, size `dims`, and global names
`a`, `b`, `c` (for example). An arbitrary number of arrays may be created.
"""
macro createarrays(T, dims, vars...)
  expr = Expr(:block)
  append!(expr.args, [:($(esc(var)) = zeros($(esc(T)), $(esc(dims))); ) for var in vars])
  expr
end

"Returns an expression that defines a Composite Type of the AbstractVars variety."
function structvarsexpr(name, physfields, transfields; vardims=2, parent=:AbstractVars, T=Float64, arraytype=:Array)
  physexprs = [:($fld::$arraytype{T,$vardims}) for fld in physfields]
  transexprs = [:($fld::$arraytype{Complex{T},$vardims}) for fld in transfields]
  expr = :(struct $name{T} <: $parent; $(physexprs...); $(transexprs...); end)
end

"""
    structvarsexpr(name, fieldspecs; parent=nothing)

Returns an expression that defines a composite type whose fields are given by
the name::type pairs specifed by the tuples in fieldspecs. The convention is
name = fieldspecs[i][1] and type = fieldspecs[i][2] for the ith element of
fieldspecs.
"""
function structvarsexpr(name, fieldspecs; parent=nothing)
  # name = spec[1]; type = spec[2]
  # example: fieldspecs[1] = (:u, Array{Float64,2})
  fieldexprs = [ :( $(spec[1])::$(spec[2]) ) for spec in fieldspecs ]
  if parent == nothing
    expr = :(struct $name{T}; $(fieldexprs...); end)
  else
    expr = :(struct $name{T} <: $parent; $(fieldexprs...); end)
  end
  expr
end

"""
    getfieldspecs(fieldnames, fieldtype)

Returns an array of (fieldname[i], fieldtype) tuples that can be given to the
function getstructexpr. This function makes it convenient to construct
fieldspecs for lists of variables of the same type.
"""
getfieldspecs(fieldnames, fieldtype) = collect(zip(fieldnames, [ fieldtype for name in fieldnames ]))

"""
    fftwavenums(n; L=1)

Return the fftwavenumber vector with length n and domain size L.
"""
fftwavenums(n::Int; L=1) = 2π/L*cat(0:n/2, -n/2+1:-1, dims=1)

"""
    rms(q)

Return the root-mean-square of an array.
"""
rms(q) = sqrt(mean(q.^2))

"""
    peakedisotropicspectrum(g, kpeak, E0; mask=mask, allones=false)

Generate a real and random two-dimensional vorticity field q(x, y) with
a Fourier spectrum peaked around a central non-dimensional wavenumber kpeak and
normalized so that its total energy is E0.
"""
function peakedisotropicspectrum(g::TwoDGrid, kpeak::Real, E0::Real; mask=ones(size(g.Kr)), allones=false)
  if g.Lx !== g.Ly
      error("the domain is not square")
  else
    k0 = kpeak*2π/g.Lx
    modk = sqrt.(g.KKrsq)
    psik = zeros(g.nk, g.nl)
    psik =  (modk.^2 .* (1 .+ (modk/k0).^4)).^(-0.5)
    psik[1, 1] = 0.0
    psih = (randn(g.nkr, g.nl)+im*randn(g.nkr, g.nl)).*psik
    if allones; psih = psik; end
    psih = psih.*mask
    Ein = real(sum(g.KKrsq.*abs2.(psih)/(g.nx*g.ny)^2))
    psih = psih*sqrt(E0/Ein)
    q = -irfft(g.KKrsq.*psih, g.nx)
  end
end

"""
    parsevalsum2(uh, g)

Calculate ∫|u|² = Σ|uh|² on a 2D grid, where uh is the Fourier transform of u.
Accounts for DFT normalization, grid resolution, and whether or not uh
is the product of fft or rfft.
"""
function parsevalsum2(uh, g::TwoDGrid)
  norm = g.Lx*g.Ly/(g.nx^2*g.ny^2)    # weird normalization for dft

  if size(uh)[1] == g.nkr             # uh is conjugate symmetric
    U = sum(abs2, uh[1, :])           # k=0 modes
    U += 2*sum(abs2, uh[2:end, :])    # sum k>0 modes twice
  else                                # count every mode once
    U = sum(abs2, uh)
  end

  norm*U
end

"""
    parsevalsum(uh, g)

Calculate real(Σ uh) on a 2D grid.  Accounts for DFT normalization,
grid resolution, and whether or not uh is in a conjugate-symmetric form to
save memory.
"""
function parsevalsum(uh, g::TwoDGrid)
  norm = g.Lx*g.Ly/(g.nx^2*g.ny^2) # weird normalization for dft

  if size(uh)[1] == g.nkr       # uh is conjugate symmetric
    U = sum(uh[1, :])           # k=0 modes
    U += 2*sum(uh[2:end, :])    # sum k>0 modes twice
  else # count every mode once
    U = sum(uh)
  end

  norm*real(U)
end

"""
    jacobianh(a, b, g)

Returns the transform of the Jacobian of two fields a, b on the grid g.
"""
function jacobianh(a, b, g::TwoDGrid)
  if eltype(a) <: Real
    bh = rfft(b)
    bx = irfft(im*g.kr.*bh, g.nx)
    by = irfft(im*g.l.*bh, g.nx)
    return im*g.kr.*rfft(a.*by)-im*g.l.*rfft(a.*bx)
  else
    # J(a, b) = dx(a b_y) - dy(a b_x)
    bh = fft(b)
    bx = ifft(im*g.k.*bh)
    by = ifft(im*g.l.*bh)
    return im*g.k.*fft(a.*by).-im*g.l.*fft(a.*bx)
  end
end

"""
    jacobian(a, b, g)

Returns the Jacobian of a and b.
"""
function jacobian(a, b, g::TwoDGrid)
  if eltype(a) <: Real
   return irfft(jacobianh(a, b, g), g.nx)
  else
   return ifft(jacobianh(a, b, g))
  end
end

"""
    radialspectrum(ah, g; n=nothing, m=nothing, refinement=2)

Returns `aρ = ∫ ah(ρ,θ) ρ dρ dθ`, the radial spectrum of `ah` known on the
Cartesian wavenumber grid (k,l).

`aρ` is found by intepolating `ah` onto a polar wavenumber grid (ρ,θ), and
then integrating over `θ` to find `aρ`. The default resolution (n,m) for the
polar wave number grid is `n=refinement*maximum(nk, nl),
m=refinement*maximum(nk, nl)`, where `refinement=2` by default. If
`ah` is in conjugate symmetric form only the upper half plane in `θ` is
represented on the polar grid.
"""
function radialspectrum(ah, g::TwoDGrid; n=nothing, m=nothing, refinement=2)

  n = n == nothing ? refinement*maximum([g.nk, g.nl]) : n
  m = m == nothing ? refinement*maximum([g.nk, g.nl]) : m

  # Calcualte shifted k and l
  lshift = range(-g.nl/2+1, stop=g.nl/2, length=g.nl)*2π/g.Ly

  if size(ah)[1] == g.nkr # conjugate symmetric form
    m = Int(m/2)                         # => half resolution in θ
    θ = range(-π/2, stop=π/2, length=m)  # θ-grid from k=0 to max(kr)
    ahshift = fftshift(ah, 2)            # shifted ah
    kshift = range(0, stop=g.nkr-1, length=g.nkr)*2π/g.Lx
  else # ordinary form
    θ = range(0, stop=2π, length=m)      # θ grid
    ahshift = fftshift(ah, [1, 2])       # shifted ah
    kshift = range(-g.nk/2+1, stop=g.nk/2, length=g.nk)*2π/g.Lx
  end

  # Interpolator for ah
  itp = scale(interpolate(ahshift, BSpline(Linear())), kshift, lshift)

  # Get radial wavenumber vector
  ρmax = minimum([(g.nk/2-1)*2π/g.Lx, (g.nl/2-1)*2π/g.Ly])
  ρ = range(0, stop=ρmax, length=n)

  # Interpolate ah onto fine grid in (ρ,θ).
  ahρθ = zeros(eltype(ahshift), (n, m))

  for i=2:n, j=1:m # ignore zeroth mode
    kk = ρ[i]*cos(θ[j])
    ll = ρ[i]*sin(θ[j])
    ahρθ[i, j] = itp(kk, ll)
  end

  # ahρ = ρ ∫ ah(ρ,θ) dθ  =>  Ah = ∫ ahρ dρ = ∫∫ ah dk dl
  dθ = θ[2]-θ[1]
  if size(ah)[1] == g.nkr
    ahρ = 2ρ.*sum(ahρθ, dims=2)*dθ # multiply by 2 for conjugate symmetry
  else
    ahρ = ρ.*sum(ahρθ, dims=2)*dθ
  end

  ahρ[1] = ah[1, 1] # zeroth mode

  ρ, ahρ
end

# Moments and cumulants
"Compute the average of `c` on the grid `g`."
domainaverage(c, g) = g.dx*g.dy*sum(c)/(g.Lx*g.Ly)

"Compute the `n`th x-moment of `c` on the grid `g`."
xmoment(c, g, n=1) = sum(g.X.^n.*c)/sum(c)

"Compute the `n`th y-moment of `c` on the grid `g`."
ymoment(c, g, n=1) = sum(g.Y.^n.*c)/sum(c)