text:
cd(@__DIR__); include("setups/simpletree.jl")
gr(dpi = 200)

## build frames
SGWT = sgwt_frame(Matrix(W); nf = 6)
SGWT = reshape(SGWT, (N, :))
SGWT_dual = (SGWT * SGWT') \ SGWT

distROT = natural_eigdist(𝚽, 𝛌, Q; α = 1.0, input_format = :pmf1, distance = :ROT)
rNGWF, dic_l2x = rngwf_all_vectors(distROT, 𝚽; σ = 0.1 * maximum(distROT), thres = 0.15)
rNGWF_dual = (rNGWF * rNGWF') \ rNGWF
Γ = rngwf_lx(dic_l2x)

## (a)
plt = plot(size = (600, 500))
    gplot!(W, X, width = 1)
    scatter_gplot!(X; marker = f, ms = 4)
    plot!(frame = :none, cbar = true, xlim = [-10, 14])
savefig(plt, "../figs/simpletree_f.png")

## (b)
rel_approx_sgwt, f_approx_sgwt = frame_approx(f, SGWT, SGWT_dual; num_kept = 3 * N)
rel_approx_rngwf, f_approx_rngwf = frame_approx(f, rNGWF, rNGWF_dual; num_kept = 3 * N)
plt = plot(0:length(rel_approx_rngwf)-1, [rel_approx_sgwt rel_approx_rngwf],
    grid = false, lw = 2, c = [:blue :green],
    xlab = "Number of Coefficients Retained",
    ylab = "Relative Approximation Error", yaxis = :log,
    lab = ["SGWT" "rNGWF"])
    plot!(xguidefontsize = 14, yguidefontsize = 14, legendfontsize = 11, size = (600, 500))
savefig(plt, "../figs/simpletree_approx_f.png")


text:
export View, isviewable, view!, next!, trytoview

import Base: size, getindex, setindex!

mutable struct View{T,N} <: AbstractArray{T,N}
    parent::AbstractArray
    v::SubArray{T,N}
    ind::Int
end

size(a::View, args...) = size(a.v, args...)
getindex(a::View, args...) = getindex(a.v, args...)
setindex!(a::View, value::AbstractArray, args...) = a.v[args...] = value
setindex!(a::View, value, args::Integer...) where N = a.v[args...] = value
setindex!(a::View, value, args...) = a.v[args...] .= value

isviewable(a::Union{DenseArray,SharedArray}{T}) where {T<:Number} = true
isviewable(a) = false

function view(a::AbstractArray{T,1}, i::Int = 1) where T
    View(a, Base.view(a, i:i), i)
end

function view(a::AbstractArray{T,2}, i::Int = 1) where T
    View(a, selectdim(a, ndims(a), i:i), i)
end

function view(a::AbstractArray{T,N}, i::Int = 1) where T where N
    View(a, selectdim(a, ndims(a), i), i)
end

function view!(a::AbstractArray{T,N}, i::Int, v::View{T,N}) where T where N
    @assert v.parent == a
    v.v = view(a, i)
    v
end

function view(a::AbstractArray, ind::UnitRange)
    selectdim(a, ndims(a), ind)
end

function next!(v::View)
    v.ind += v.ind < len(v.parent) ? 1 : 0
    if ndims(v.parent) <= 2
        v.v = selectdim(v.parent, ndims(v.parent), v.ind:v.ind)
    else
        v.v = selectdim(v.parent, ndims(v.parent), v.ind)
    end
    v
end

trytoview(a::Vector{T}, i = 1) where {T<:Number} = at(a,i)
trytoview(a::Vector{T}, i, v::View) where {T<:Number} = at(a,i)
trytoview(a::DenseArray{T}, i = 1) where {T<:Number} = view(a, i)
trytoview(a::DenseArray{T}, i, v::View) where {T<:Number} = view!(a, i, v)
trytoview(a, i) = at(a, i)
trytoview(a, i, v) = at(a, i)
# trytoview(a, i) = view(a, i)
# trytoview(a, i, v) = view(a, i)

import Base.read!
export read!
read!(io, v::View) = read!(io, v.v)
function read!(io, v::SubArray)
    if v.stride1 <= 1
        GC.@preserve v unsafe_read(io, pointer(v), sizeof(v))
    else
        error("FunctionalData: There is no suitable read! for SubArrays with v.stride != 1")
    end
end

ext: jl


text:
#=
Rejection sampler based on algorithm from Robert (1995)

  - Available at http://arxiv.org/abs/0907.4010

Implementation from Distributions.jl, with few modifications to make it more generic.
Also here we avoid the overhead of computing the CDF, which `truncated(Normal(), a, b)`
always does when it constructs the truncated distribution (`tp`).
We specialize to the case where `b = Inf`.
=#

"""
    randnt([rng], a)

Random standard normal lower truncated at `a` (that is, Z ≥ a).
"""
randnt(rng::AbstractRNG, a::Real) = randnt(rng, float(a))
randnt(rng::AbstractRNG, a::BigFloat) = randnt(rng, Float64(a))
randnt(a::Real) = randnt(GLOBAL_RNG, a)

function randnt(rng::AbstractRNG, a::Base.IEEEFloat)
    if a ≤ 0
        while true
            r = randn(rng, typeof(a))
            r ≥ a && return r
        end
    else
        t = sqrt1half(a)
        !(t < Inf) && return a
        while true
            r = a + randexp(rng, typeof(a)) / t
            u = rand(rng, typeof(a))
            if u < exp(-(r - t)^2 / 2)
                return r
            end
        end
    end
end

"""
    sqrt1half(x)

Accurate computation of sqrt(1 + (x/2)^2) + |x|/2.
"""
sqrt1half(x::Real) = _sqrt1half(float(abs(x)))

function _sqrt1half(x::Real)
    if x > 2/sqrt(eps(x))
        return x
    else
        return sqrt(one(x) + (x/2)^2) + x/2
    end
end

"""
    randnt_half([rng], μ, σ)

Samples the normal distribution with mean `μ` and standard deviation `σ`
truncated to positive values.
"""
function randnt_half(rng::AbstractRNG, μ::Real, σ::Real)
    z = randnt(rng, -μ / σ)
    return μ + σ * z
end

randnt_half(μ::Real, σ::Real) = randnt_half(GLOBAL_RNG, μ, σ)

"""
    tnmean(a)

Mean of the standard normal distribution,
truncated to the interval (a, +∞).
"""
tnmean(a::Real) = sqrt(two(a)/π) / erfcx(a/√two(a))

"""
    tnvar(a)

Variance of the standard normal distribution,
truncated to the interval (a, +∞).
WARNING: Fails for very very large values of `a`.
"""
function tnvar(a::Real)
    μ = tnmean(a)
    return one(μ) - (μ - a) * μ
end

"""
    tnmeanvar(a)

Mean and variance of the standard normal distribution truncated to the interval (a, +∞).
Equivalent to `tnmean(a), tnvar(a)` but saves some common computations.
WARNING: `tnvar(a) can fail for very very large values of `a`.
"""
function tnmeanvar(a::Real)
    μ = tnmean(a)
    ν = one(μ) - (μ - a) * μ
    return μ, ν
end

#= tnmeanvar.(A) produces an array of tuples, which is usually not what you want.
So I have this method which broadcasts and returns the tuple of arrays instead. =#
function tnmeanvar(a::AbstractArray)
    μ = tnmean.(a)
    ν = one.(μ) - (μ - a) .* μ
    return μ, ν
end


text:
using Documenter
using CombinatorialSpaces

makedocs(
  sitename = "CombinatorialSpaces.jl",
  format = Documenter.HTML(),
  modules = [CombinatorialSpaces],
  checkdocs = :exports,
  pages = [
    "simplicial_sets.md",
    "discrete_exterior_calculus.md",
    "combinatorial_maps.md",
    "meshes.md"
  ]
)

deploydocs(
  target = "build",
  repo = "github.com/AlgebraicJulia/CombinatorialSpaces.jl.git",
  branch = "gh-pages",
  devbranch = "main"
)


text:
# Julia wrapper for header: /Volumes/BOOTCAMP/programs/compa_libs/libim_/include/im.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@static Sys.iswindows() ? (const libim_ = "im") : (const libim_ = "libim")  # Name of im shared lib.

const IM_INTERRUPTED = -1
const IM_ALL = -1
const IM_COUNTER_CB = 0
const IM_RESOLUTION_CB = 1
const IM_GIF_TRANSPARENT_COLOR_CB = 0
const IM_TIF_IMAGE_DESCRIPTION_CB = 0

# begin enum imDataType
const imDataType = UInt32
const IM_BYTE = 0
const IM_SHORT = 1
const IM_USHORT = 2
const IM_INT = 3
const IM_FLOAT = 4
const IM_CFLOAT = 5
# end enum imDataType

# begin enum imColorSpace
const imColorSpace = UInt32
const IM_RGB = 0
const IM_MAP = 1
const IM_GRAY = 2
const IM_BINARY = 3
const IM_CMYK = 4
const IM_YCBCR = 5
const IM_LAB = 6
const IM_LUV = 7
const IM_XYZ = 8
# end enum imColorSpace

# begin enum imColorModeConfig
const imColorModeConfig = UInt32
const IM_ALPHA = 256
const IM_PACKED = 512
const IM_TOPDOWN = 1024
# end enum imColorModeConfig

# begin enum imErrorCodes
const imErrorCodes = UInt32
const IM_ERR_NONE = 0
const IM_ERR_OPEN = 1
const IM_ERR_ACCESS = 2
const IM_ERR_FORMAT = 3
const IM_ERR_DATA = 4
const IM_ERR_COMPRESS = 5
const IM_ERR_MEM = 6
const IM_ERR_COUNTER = 7
const IM_ERR_READ = IM_ERR_ACCESS
const IM_ERR_WRITE = IM_ERR_ACCESS
const IM_ERR_TYPE = IM_ERR_DATA
const IM_ERR_COMP = IM_ERR_COMPRESS
# end enum imErrorCodes
const imFile = Cvoid

# begin enum ANONYMOUS_1
const ANONYMOUS_1 = UInt32
const IM_BMP = 0
const IM_PCX = 1
const IM_GIF = 2
const IM_TIF = 3
const IM_RAS = 4
const IM_SGI = 5
const IM_JPG = 6
const IM_LED = 7
const IM_TGA = 8
# end enum ANONYMOUS_1

# begin enum ANONYMOUS_2
const ANONYMOUS_2 = UInt32
const IM_NONE = 0
const IM_DEFAULT = 256
const IM_COMPRESSED = 512
# end enum ANONYMOUS_2
const imCallback = Ptr{Cvoid}
const imFileCounterCallback = Ptr{Cvoid}
const imResolutionCallback = Ptr{Cvoid}

# begin enum ANONYMOUS_3
const ANONYMOUS_3 = UInt32
const IM_RES_NONE = 0
const IM_RES_DPI = 1
const IM_RES_DPC = 2
# end enum ANONYMOUS_3
const imGifTranspIndex = Ptr{Cvoid}
const imTiffImageDesc = Ptr{Cvoid}

function imFileOpen(file_name::String, error::Ptr{Cint})
  ccall( (:imFileOpen, libim_), Ptr{imFile}, (Ptr{UInt8}, Ptr{Cint}), file_name, error)
end

function imFileOpenAs(file_name::String, format::String, error::Ptr{Cint})
  ccall( (:imFileOpenAs, libim_), Ptr{imFile}, (Ptr{UInt8}, Ptr{UInt8}, Ptr{Cint}), file_name, format, error)
end

function imFileNew(file_name::String, format::String, error::Ptr{Cint})
  ccall( (:imFileNew, libim_), Ptr{imFile}, (Ptr{UInt8}, Ptr{UInt8}, Ptr{Cint}), file_name, format, error)
end

function imFileClose(ifile::Ptr{imFile})
  ccall( (:imFileClose, libim_), Cvoid, (Ptr{imFile},), ifile)
end

function imFileHandle(ifile::Ptr{imFile}, index::Integer)
  ccall( (:imFileHandle, libim_), Ptr{Cvoid}, (Ptr{imFile}, Cint), ifile, index)
end

function imFileGetInfo(ifile::Ptr{imFile}, format::String, compression::String, image_count::Ptr{Cint})
  ccall( (:imFileGetInfo, libim_), Cvoid, (Ptr{imFile}, Ptr{UInt8}, Ptr{UInt8}, Ptr{Cint}), ifile, format, compression, image_count)
end

function imFileSetInfo(ifile::Ptr{imFile}, compression::String)
  ccall( (:imFileSetInfo, libim_), Cvoid, (Ptr{imFile}, Ptr{UInt8}), ifile, compression)
end

function imFileSetAttribute(ifile::Ptr{imFile}, attrib::String, data_type::Integer, count::Integer, data::Ptr{Cvoid})
  ccall( (:imFileSetAttribute, libim_), Cvoid, (Ptr{imFile}, Ptr{UInt8}, Cint, Cint, Ptr{Cvoid}), ifile, attrib, data_type, count, data)
end

function imFileGetAttribute(ifile::Ptr{imFile}, attrib::String, data_type::Ptr{Cint}, count::Ptr{Cint})
  ccall( (:imFileGetAttribute, libim_), Ptr{Cvoid}, (Ptr{imFile}, Ptr{UInt8}, Ptr{Cint}, Ptr{Cint}), ifile, attrib, data_type, count)
end

function imFileGetAttributeList(ifile::Ptr{imFile}, attrib::Ptr{Ptr{UInt8}}, attrib_count::Ptr{Cint})
  ccall( (:imFileGetAttributeList, libim_), Cvoid, (Ptr{imFile}, Ptr{Ptr{UInt8}}, Ptr{Cint}), ifile, attrib, attrib_count)
end

function imFileGetPalette(ifile::Ptr{imFile}, palette::Ptr{Clong}, palette_count::Ptr{Cint})
  ccall( (:imFileGetPalette, libim_), Cvoid, (Ptr{imFile}, Ptr{Clong}, Ptr{Cint}), ifile, palette, palette_count)
end

function imFileSetPalette(ifile::Ptr{imFile}, palette::Ptr{Clong}, palette_count::Integer)
  ccall( (:imFileSetPalette, libim_), Cvoid, (Ptr{imFile}, Ptr{Clong}, Cint), ifile, palette, palette_count)
end

function imFileReadImageInfo(ifile::Ptr{imFile}, index::Integer, width::Ptr{Cint}, height::Ptr{Cint}, file_color_mode::Ptr{Cint}, file_data_type::Ptr{Cint})
  ccall( (:imFileReadImageInfo, libim_), Cint, (Ptr{imFile}, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), ifile, index, width, height, file_color_mode, file_data_type)
end

function imFileWriteImageInfo(ifile::Ptr{imFile}, width::Integer, height::Integer, user_color_mode::Integer, user_data_type::Integer)
  ccall( (:imFileWriteImageInfo, libim_), Cint, (Ptr{imFile}, Cint, Cint, Cint, Cint), ifile, width, height, user_color_mode, user_data_type)
end

function imFileReadImageData(ifile::Ptr{imFile}, data::Ptr{Cvoid}, convert2bitmap::Integer, color_mode_flags::Integer)
  ccall( (:imFileReadImageData, libim_), Cint, (Ptr{imFile}, Ptr{Cvoid}, Cint, Cint), ifile, data, convert2bitmap, color_mode_flags)
end

function imFileWriteImageData(ifile::Ptr{imFile}, data::Ptr{Cvoid})
  ccall( (:imFileWriteImageData, libim_), Cint, (Ptr{imFile}, Ptr{Cvoid}), ifile, data)
end

function imFormatRegisterInternal()
  ccall( (:imFormatRegisterInternal, libim_), Cvoid, (), )
end

function imFormatRemoveAll()
  ccall( (:imFormatRemoveAll, libim_), Cvoid, (), )
end

function imFormatList(format_list::Ptr{Ptr{UInt8}}, format_count::Ptr{Cint})
  ccall( (:imFormatList, libim_), Cvoid, (Ptr{Ptr{UInt8}}, Ptr{Cint}), format_list, format_count)
end

function imFormatInfo(format::String, desc::String, ext::String, can_sequence::Ptr{Cint})
  ccall( (:imFormatInfo, libim_), Cint, (Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{Cint}), format, desc, ext, can_sequence)
end

function imFormatInfoExtra(format::String, extra::String)
  ccall( (:imFormatInfoExtra, libim_), Cint, (Ptr{UInt8}, Ptr{UInt8}), format, extra)
end

function imFormatCompressions(format::String, comp::Ptr{Ptr{UInt8}}, comp_count::Ptr{Cint}, color_mode::Integer, data_type::Integer)
  ccall( (:imFormatCompressions, libim_), Cint, (Ptr{UInt8}, Ptr{Ptr{UInt8}}, Ptr{Cint}, Cint, Cint), format, comp, comp_count, color_mode, data_type)
end

function imFormatCanWriteImage(format::String, compression::String, color_mode::Integer, data_type::Integer)
  ccall( (:imFormatCanWriteImage, libim_), Cint, (Ptr{UInt8}, Ptr{UInt8}, Cint, Cint), format, compression, color_mode, data_type)
end

function imEncodeColor(red::Cuchar, green::Cuchar, blue::Cuchar)
  ccall( (:imEncodeColor, libim_), Clong, (Cuchar, Cuchar, Cuchar), red, green, blue)
end

function imDecodeColor(red::Ptr{Cuchar}, green::Ptr{Cuchar}, blue::Ptr{Cuchar}, palette::Clong)
  ccall( (:imDecodeColor, libim_), Cvoid, (Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Clong), red, green, blue, palette)
end

function imFileFormat(filename::String, format::Ptr{Cint})
  ccall( (:imFileFormat, libim_), Cint, (Ptr{UInt8}, Ptr{Cint}), filename, format)
end

function imImageInfo(filename::String, width::Ptr{Cint}, height::Ptr{Cint}, _type::Ptr{Cint}, palette_count::Ptr{Cint})
  ccall( (:imImageInfo, libim_), Cint, (Ptr{UInt8}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), filename, width, height, _type, palette_count)
end

function imLoadRGB(filename::String, red::Ptr{Cuchar}, green::Ptr{Cuchar}, blue::Ptr{Cuchar})
  ccall( (:imLoadRGB, libim_), Cint, (Ptr{UInt8}, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}), filename, red, green, blue)
end

function imSaveRGB(width::Integer, height::Integer, format::Integer, red::Ptr{Cuchar}, green::Ptr{Cuchar}, blue::Ptr{Cuchar}, filename::String)
  ccall( (:imSaveRGB, libim_), Cint, (Cint, Cint, Cint, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{UInt8}), width, height, format, red, green, blue, filename)
end

function imLoadMap(filename::String, map::Ptr{Cuchar}, palette::Ptr{Clong})
  ccall( (:imLoadMap, libim_), Cint, (Ptr{UInt8}, Ptr{Cuchar}, Ptr{Clong}), filename, map, palette)
end

function imSaveMap(width::Integer, height::Integer, format::Integer, map::Ptr{Cuchar}, palette_count::Integer, palette::Ptr{Clong}, filename::String)
  ccall( (:imSaveMap, libim_), Cint, (Cint, Cint, Cint, Ptr{Cuchar}, Cint, Ptr{Clong}, Ptr{UInt8}), width, height, format, map, palette_count, palette, filename)
end

function imRGB2Map(width::Integer, height::Integer, red::Ptr{Cuchar}, green::Ptr{Cuchar}, blue::Ptr{Cuchar}, map::Ptr{Cuchar}, palette_count::Integer, palette::Ptr{Clong})
  ccall( (:imRGB2Map, libim_), Cvoid, (Cint, Cint, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Cint, Ptr{Clong}), width, height, red, green, blue, map, palette_count, palette)
end

function imMap2RGB(width::Integer, height::Integer, map::Ptr{Cuchar}, palette_count::Integer, colors::Ptr{Clong}, red::Ptr{Cuchar}, green::Ptr{Cuchar}, blue::Ptr{Cuchar})
  ccall( (:imMap2RGB, libim_), Cvoid, (Cint, Cint, Ptr{Cuchar}, Cint, Ptr{Clong}, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}), width, height, map, palette_count, colors, red, green, blue)
end

function imRGB2Gray(width::Integer, height::Integer, red::Ptr{Cuchar}, green::Ptr{Cuchar}, blue::Ptr{Cuchar}, map::Ptr{Cuchar}, grays::Ptr{Clong})
  ccall( (:imRGB2Gray, libim_), Cvoid, (Cint, Cint, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Clong}), width, height, red, green, blue, map, grays)
end

function imMap2Gray(width::Integer, height::Integer, map::Ptr{Cuchar}, palette_count::Integer, colors::Ptr{Clong}, grey_map::Ptr{Cuchar}, grays::Ptr{Clong})
  ccall( (:imMap2Gray, libim_), Cvoid, (Cint, Cint, Ptr{Cuchar}, Cint, Ptr{Clong}, Ptr{Cuchar}, Ptr{Clong}), width, height, map, palette_count, colors, grey_map, grays)
end

function imResize(src_width::Integer, src_height::Integer, src_map::Ptr{Cuchar}, dst_width::Integer, dst_height::Integer, dst_map::Ptr{Cuchar})
  ccall( (:imResize, libim_), Cvoid, (Cint, Cint, Ptr{Cuchar}, Cint, Cint, Ptr{Cuchar}), src_width, src_height, src_map, dst_width, dst_height, dst_map)
end

function imStretch(src_width::Integer, src_height::Integer, src_map::Ptr{Cuchar}, dst_width::Integer, dst_height::Integer, dst_map::Ptr{Cuchar})
  ccall( (:imStretch, libim_), Cvoid, (Cint, Cint, Ptr{Cuchar}, Cint, Cint, Ptr{Cuchar}), src_width, src_height, src_map, dst_width, dst_height, dst_map)
end

function imRegisterCallback(cb::imCallback, cb_id::Integer, format::Integer)
  ccall( (:imRegisterCallback, libim_), Cint, (imCallback, Cint, Cint), cb, cb_id, format)
end


text:
using ParserCombinator.Parsers.GML

example1 = "
graph [
  comment \"This is a sample graph\"
  directed 1
  id 42
  label \"Hello, I am a graph\"
  node [
    id 1
    label \"node 1\"
    thisIsASampleAttribute 42
  ]
  node [
    id 2
    label \"node 2\"
    thisIsASampleAttribute 43
  ]
  node [
    id 3
    label \"node 3\"
    thisIsASampleAttribute 44
  ]
  edge [
    source 1
    target 2
    label \"Edge from node 1 to node 2\"
  ]
  edge [
    source 2
    target 3
    label \"Edge from node 2 to node 3\"
  ]
  edge [
    source 3
    target 1
    label \"Edge from node 3 to node 1\"
  ]
]
"

root = parse_dict(example1)

for graph in root[:graph]
    println("graph $(graph[:label])")
    for node in graph[:node]
        println(" node $(node[:id])")
    end
    for edge in graph[:edge]
        println(" edge $(edge[:label]): $(edge[:source]) - $(edge[:target])")
    end
end


text:
# functions for applying SBP operators in Cartesian directions (rather than
# parametric directions).  These functions should eventually be migrated into
# SummationByParts.jl

#------------------------------------------------------------------------------
# Functions for computing [Dx * w, Dy * w, Dz * w]

"""
  Computes Qx^T * w in all Cartesian directions.

  **Inputs**

   * sbp: the SBP operator
   * w: `numDofPerNode` x `numNodesPerElement` array of values
   * dxidx: `dim` x `dim` x `numNodesPerElement` array containing the metrics
   * op: a UnaryFunctor that determines if values are added or subtracted to
         `wx`

  **Inputs/Outputs**

   * wxi: temporary array, `numDofPerNode` x `numNodesPerElement` x `dim`
   * wx: array updates with the result, same size as `wxi`.
"""
function applyQxTransposed(sbp, w::AbstractMatrix, dxidx::Abstract3DArray,
                             wxi::Abstract3DArray, wx::Abstract3DArray,
                             op::SummationByParts.UnaryFunctor=SummationByParts.Add())

  # compute dw/dxi
  numDofPerNode, numNodesPerElement, dim = size(wx)

  wxi_d = sview(wxi, :, :, 1)
  for d1=1:dim  # cartesian direction

    for d2=1:dim  # parametric dimension (summed)
      for j=1:numNodesPerElement
        for k=1:numDofPerNode
          wxi_d[k, j] = op(dxidx[d2, d1, j]*w[k, j])
        end
      end

      Q_d = sview(sbp.Q, :, :, d2)
      wx_d = sview(wx, :, :, d1)
      smallmatmat_kernel!(wxi_d, Q_d, wx_d, 1, 1)
    end
  end

  return nothing
end

"""
  Computes Dx * w in all Cartian directions.  See [`applyQxTranspose`](@ref)
  for most of the arguments.  Note that the user is responsible for zeroing out
  `wx` if required.

  **Inputs**

   * sbp
   * w
   * dxidx
   * jac: determinant of the mapping Jacobian, length `numNodesPerElement`.
          This is required because `dxidx` contains a factor of `1/|J|`
   * op

  **Inputs/Outputs**

   * wxi
   * wx
"""
function applyDx(sbp, w::AbstractMatrix, dxidx::Abstract3DArray,
                 jac::AbstractVector,
                 wxi::Abstract3DArray, wx::Abstract3DArray,
                 op::SummationByParts.UnaryFunctor=SummationByParts.Add())

  numDofPerNode, numNodesPerElement, dim = size(wx)
  fill!(wxi, 0)
  for d=1:dim
    differentiateElement!(sbp, d, w, sview(wxi, :, :, d))
  end

  @simd for d1=1:dim
    @simd for j=1:numNodesPerElement
      @simd for d2=1:dim
        fac = jac[j]*dxidx[d2, d1, j]
        @simd for k=1:numDofPerNode
          wx[k, j, d1] += op(wxi[k, j, d2] * fac)
        end
      end
    end
  end

  return nothing
end

function applyDxTransposed(sbp, w::AbstractMatrix, dxidx::Abstract3DArray,
                 jac::AbstractVector,
                 wxi::Abstract3DArray, wx::Abstract3DArray,
                 op::SummationByParts.UnaryFunctor=SummationByParts.Add())

  numDofPerNode, numNodesPerElement, dim = size(wx)

  # Dx^T * w = (Dxi^T * |dxi/dx| * u + Deta^T * |deta/dx|* w)
  # so rotate w using mapping jacobian first, then apply D
  for d1=1:dim

    for d2=1:dim
      for j=1:numNodesPerElement
        for k=1:numDofPerNode
          # we only need temporary storage in the shape of w, but to be
          # consistent with the other functions we use the wxi array
          wxi[k, j, d2] = w[k, j]*(dxidx[d2, d1, j]*jac[j])
        end
      end

      differentiateElement!(sbp, d2, ro_sview(wxi, :, :, d2), sview(wx, :, :, d1), op, true)
    end  # end d2
  end  # end d1

  return nothing
end

"""
  Computes Qx * w in all Cartesian directions.  See [`applyQxTranspose`](@ref)
  for arguments.  Note that the user must zero out `wx` if required.
"""
function applyQx(sbp, w::AbstractMatrix, dxidx::Abstract3DArray,
                 wxi::Abstract3DArray, wx::Abstract3DArray,
                 op::SummationByParts.UnaryFunctor=SummationByParts.Add())

  numDofPerNode, numNodesPerElement, dim = size(wx)
  fill!(wxi, 0)
  for d=1:dim
    weakDifferentiateElement!(sbp, d, w, sview(wxi, :, :, d))
  end

  @simd for d1=1:dim
    @simd for j=1:numNodesPerElement
      @simd for d2=1:dim
        @simd for k=1:numDofPerNode
          wx[k, j, d1] += op(wxi[k, j, d2] * dxidx[d2, d1, j])
        end
      end
    end
  end

  return nothing
end

#------------------------------------------------------------------------------
# Functions that apply the operators to different vectors in each direction and
# sums the result.
# Eg. given wx, wy, wz, compute w = Dx * wx + Dy * wy + Dz * wz

"""
  Performs  res += Qx^T * w[:, :, 1] + Qy^T * w[:, :, 2] (and the z contribution
  in 3D.

  Unlike the other method, `w` is a 3D array and `wx` is a a 2D array

  **Inputs**

   * sbp
   * w: `numDofPerNode` x `numNodesPerElement` x `dim` containing the values
        for each dimension.
   * dxidx: the metric terms, same as other method
   * op: a UnaryFunctor that determines if values are added or subtracted to
         `wx`

  **Inputs/Outputs**

   * wxi: work array, same as other method
   * wx: array to be updated with the result, `numDofPerNode` x
         `numNodesPerElement`.

"""
function applyQxTransposed(sbp, w::Abstract3DArray, dxidx::Abstract3DArray,
                             wxi::Abstract3DArray, wx::AbstractMatrix,
                             op::SummationByParts.UnaryFunctor=SummationByParts.Add())

  numDofPerNode, numNodesPerElement, dim = size(wxi)

  wxi_d = sview(wxi, :, :, 1)
  for d1=1:dim  # cartesian direction
    for d2=1:dim  # parametric dimension (summed)

      for j=1:numNodesPerElement
        for k=1:numDofPerNode
          wxi_d[k, j] = op(dxidx[d2, d1, j]*w[k, j, d1])
        end
      end

      Q_d = sview(sbp.Q, :, :, d2)
      smallmatmat_kernel!(wxi_d, Q_d, wx, 1, 1)
    end
  end

  return nothing
end

"""
  Like [`applyQxTransposed`](@ref), but applies D instead of Q^T.  See
  that function for details.  Note that the user must zero out `wx` if
  required.

  **Inputs**

   * sbp
   * w
   * dxidx
   * jac: determinant of the mapping Jacobian, length `numNodesPerElement`.
   * op

  **Inputs/Outputs**

   * wxi
   * wx
"""
function applyDx(sbp, w::Abstract3DArray, dxidx::Abstract3DArray,
                 jac::AbstractVector,
                 wxi::Abstract3DArray, wx::AbstractMatrix,
                 op::SummationByParts.UnaryFunctor=SummationByParts.Add())

  numDofPerNode, numNodesPerElement, dim = size(w)

  @simd for d1=1:dim
    fill!(wxi, 0)
    for d2=1:dim
      differentiateElement!(sbp, d2, sview(w, :, :, d1), sview(wxi, :, :, d2))
    end

    @simd for j=1:numNodesPerElement
      @simd for d2=1:dim
        fac = jac[j]*dxidx[d2, d1, j]
        @simd for k=1:numDofPerNode
          wx[k, j] += op(wxi[k, j, d2] * fac)
        end
      end
    end
  end

  return nothing
end

function applyDxTransposed(sbp, w::Abstract3DArray, dxidx::Abstract3DArray,
                 jac::AbstractVector,
                 wxi::Abstract3DArray, wx::AbstractMatrix,
                 op::SummationByParts.UnaryFunctor=SummationByParts.Add())

  numDofPerNode, numNodesPerElement, dim = size(w)

  for d1=1:dim

    for d2=1:dim
      for j=1:numNodesPerElement
        for k=1:numDofPerNode
          # we only need temporary storage in the shape of w, but to be
          # consistent with the other functions we use the wxi array
          wxi[k, j, d2] = w[k, j, d1]*(dxidx[d2, d1, j]*jac[j])
        end
      end

      differentiateElement!(sbp, d2, ro_sview(wxi, :, :, d2), wx, op, true)
    end  # end d2
  end  # end d1

  return nothing
end

"""
  Like [`applyQxTransposed`](@ref), but applies Q instead of Q^T.  See
  that function for details.  Note that the user must zero out `wx` if
  required.
"""
function applyQx(sbp, w::Abstract3DArray, dxidx::Abstract3DArray,
                 wxi::Abstract3DArray, wx::AbstractMatrix,
                 op::SummationByParts.UnaryFunctor=SummationByParts.Add())

  numDofPerNode, numNodesPerElement, dim = size(w)

  @simd for d1=1:dim
    fill!(wxi, 0)
    for d2=1:dim
      weakDifferentiateElement!(sbp, d2, sview(w, :, :, d1), sview(wxi, :, :, d2))
    end

    @simd for d2=1:dim
      @simd for j=1:numNodesPerElement
        @simd for k=1:numDofPerNode
          wx[k, j] += op(wxi[k, j, d2] * dxidx[d2, d1, j])
        end
      end
    end
  end

  return nothing
end


text:
using BinaryBuilder

# Collection of pre-build pandoc binaries
name = "pandoc"
version = v"2.11.4"
pandoc_ver = "2.11.4"

url_prefix = "https://github.com/jgm/pandoc/releases/download/$(pandoc_ver)/pandoc-$(pandoc_ver)"
sources = [
    ArchiveSource("$(url_prefix)-linux-amd64.tar.gz", "b15ce6009ab833fb51fc472bf8bb9683cd2bd7f8ac948f3ddeb6b8f9a366d69a"; unpack_target = "x86_64-linux-gnu"),
    ArchiveSource("$(url_prefix)-macOS.zip", "13b8597860afa6ab802993a684b340be3f31f4d2a06c50b6601f9e726cf76f71"; unpack_target = "x86_64-apple-darwin14"),
    ArchiveSource("$(url_prefix)-windows-x86_64.zip", "ee1b0c4d0f539ee8316d6cebb29f6aa709aa3a72be2b2b4ab6e9e4a77a01a50b"; unpack_target = "x86_64-w64-mingw32"),
    FileSource("https://raw.githubusercontent.com/jgm/pandoc/54d8c6959cc53e78cfea2093b33504672f81ed74/COPYING.md", "2b0d4dda4bf8034e1506507a67f80f982131137afe62bf144d248f9faea31da4"),
    FileSource("https://raw.githubusercontent.com/jgm/pandoc/54d8c6959cc53e78cfea2093b33504672f81ed74/COPYRIGHT", "7087ebeaa4c3131b7a989ea4310da1442610e9e2874c3bc9d3e8a52d618fd97a"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/
mkdir -p "${bindir}"
if [[ "${target}" != *-mingw* ]]; then
    subdir="bin/"
fi
cp ${target}/pandoc-*/${subdir}pandoc${exeext} ${bindir}
chmod +x ${bindir}/*
install_license COPYRIGHT
install_license COPYING.md
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"),
    Platform("x86_64", "macos"),
    Platform("x86_64", "windows"),
]

# The products that we will ensure are always built
products = [
    ExecutableProduct("pandoc", :pandoc),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
function is_degenerate(sides)
    any(v->(v == 0), sides) ||
    sides[1] + sides[2] < sides[3] ||
    sides[2] + sides[3] < sides[1] ||
    sides[3] + sides[1] < sides[2]
end

function num_unique(sides)
    length(unique(sides))
end

function is_equilateral(sides)
    !is_degenerate(sides) &&
    num_unique(sides) == 1
end

function is_isosceles(sides)
    !is_degenerate(sides) &&
    num_unique(sides) <= 2
end

function is_scalene(sides)
    !is_degenerate(sides) &&
    num_unique(sides) == 3
end

max_stars_repo_name: cjtallman/exercism

max_issues_repo_name: cjtallman/exercism

max_forks_repo_name: cjtallman/exercism


text:
using WAV
using PlotlyJS
using DSP
using DSP.Filters

function center_clip(xs; thresh=0.3)
    cl = thresh * maximum(abs.(xs))
    [abs(x) > cl ? (x - (sign(x) * cl)) : 0.0 for x in xs]
end

function clip(xs)
    max.(min.(xs, 1.), -1.)
end

function freq_trace(ys, fs; name="spectrum")
    nbins = length(ys)
    xs = round.(Int, fs / 2 * (0:nbins-1) ./ (nbins))
    spec = 20.*log10.(ys)
    scatter(x=xs, y=spec, name=name)
end

function fft_trace(x, fs; name="spectrum")
    ys = abs.(fftshift(fft(x)))
    nbins = div(length(ys), 2)
    xs = round.(Int, fs / 2 * (0:nbins-1) ./ (nbins))
    spec = 20.*log10.(ys[nbins+1:end])
    scatter(x=xs, y=spec, name=name)
end

function lpc_trace(x, fs, lpcOrder; name="LP order $lpcOrder")
    A, erra = lpc(x, lpcOrder)
    G = sqrt(erra)
    ws = linspace(0, π, length(x))
    Aw = zeros(Complex, size(ws))
    for (i, w) in enumerate(ws)
        Aw[i] = 1.
        for k = 1:length(A)
            Aw[i] += A[k] * exp(-im * w * k)
        end
    end
    freq_trace(abs.(G./Aw), fs; name=name)
end

function plot_spectrum(xs, fs)
    if eltype(xs) <: Vector
        traces = [fft_trace(x, fs) for x in xs]
    else
        traces = [fft_trace(xs, fs)]
    end
    plot(traces,
        Layout(
            yaxis_title="magnitude (dB)",
            xaxis_type="log",
            xaxis_title="frequency (Hz)"
        )
    )
end

function plot_lpc_spectrum(xs, fs; lpcOrders=[10,20,30], title="LP spectral envelope")
    plot([fft_trace(xs, fs; name="actual"); [lpc_trace(x, fs, o) for o in lpcOrders]],
        Layout(
            title=title,
            yaxis_title="magnitude (dB)",
            xaxis_type="log",
            xaxis_title="frequency (Hz)"
        )
    )
end

# http://www.emptyloop.com/technotes/a%20tutorial%20on%20linear%20prediction%20and%20levinson-durbin.pdf
# https://gist.github.com/jameslyons/8779402
function levinson{T<:Real}(r::Vector{T}, m::Int)
    A = zeros(m, m)
    err = zeros(1, m)

    # for k = 1
    A[1, 1] = - (r[2] / r[1])
    err[1] = r[1] * (1 - A[1, 1] ^ 2)

    # for k = 2,3,4,..m
    for k = 2:m
        λ = -(sum(A[k-1, 1:k-1] .* r[k:-1:2]) + r[k+1])/err[k-1]
        A[k, k] = λ
        A[k, 1:k-1] .= A[k-1, 1:k-1] + λ * A[k-1, k-1:-1:1]
        err[k] = err[k-1] * (1 - λ^2)
    end
    A[m, :], err[m]
end

function pitch_period_samples(xs; min_period_samples=50)
    findmax(xcorr(xs, xs)[length(xs)+min_period_samples:end])[2] + min_period_samples
end

function impulse_train(n; period=n)
    nrepeats = div(n + period - 1, period)
    repeat([1.; [0. for _ = 1:period-1]], outer=nrepeats)[1:n]
end

function lpc(xs, m)
    r = xcorr(xs, xs)[length(xs):end]
    A, err = levinson(r, m)
    return A, err
end

function zero_crossing_rate(x)
    sum(abs.(diff(0.5 .* sign.(x)))) / length(x)
end

@enum FrameType VoicedFrame UnvoicedFrame

function dump_voicing_labels(x, fs, block_size_samples, labels, outfile)
    markers = Dict{UInt32, WAVMarker}()
    j = 0
    last_label = labels[1]
    i = 1
    for l in labels[2:end]
        if l != last_label
            markers[length(markers)] = WAVMarker(
                last_label == VoicedFrame ? "voiced" : "unvoiced",
                floor(Int, j * block_size_samples),
                floor(Int, (i - j) * block_size_samples)
            )
            j = i
            last_label = l
        end
        i += 1
    end

    l = labels[end]
    markers[length(markers)] = WAVMarker(
        last_label == VoicedFrame ? "voiced" : "unvoiced",
        floor(Int, j * block_size_samples),
        floor(Int, (length(labels) - j) * block_size_samples)
    )

    out_chunks = wav_cue_write(markers)
    wavwrite(x, outfile, Fs=fs, compression=WAV.WAVE_FORMAT_PCM, chunks=out_chunks)
end

function lpc_encode_decode(xs, fs; excitationFile=nothing, lpcOrder=30, blockSizeMs=15, unvoicedZCR=0.12, unvoicedHz=300.)
    out = zeros(xs)

    if excitationFile != nothing
        excitationSamples, excitationFs = wavread(excitationFile)
        assert(excitationFs == fs)
    end
    n = ceil(Int, blockSizeMs / 1000. * fs)
    noverlap = div(n, 2)
    nblocks = div(length(xs), n)
    frameTypes = []

    index = 1
    lastOutput = zeros(n)
    while index + n < length(xs)
        block = xs[index:(index+n-1)]
        block = hamming(n) .* block
        z = zero_crossing_rate(block)
        A, erra = lpc(block, lpcOrder)
        G = sqrt(erra)
        p = pitch_period_samples(center_clip(block))
        freq_hz = fs / p
        ft = (z > unvoicedZCR || freq_hz > unvoicedHz) ? UnvoicedFrame : VoicedFrame

        if excitationFile == nothing
            if ft == VoicedFrame
                excitation = impulse_train(n; period=(p*4))
            else
                excitation = (rand(n) - 0.5) * 0.1
            end
        else
            excitation = excitationSamples[index:(index+n-1)]
        end
        push!(frameTypes, ft)
        f = PolynomialRatio([G], [1; A])
        predicted = hamming(n) .* filt(f, excitation)
        out[index:(index+noverlap-1)] = lastOutput[(noverlap+1):n]
        out[index:(index+n-1)] .+= predicted
        lastOutput = predicted
        index += noverlap
    end

    return out, noverlap, frameTypes
end


text:
##########################################################################
#
# Benchmark optimization algorithms by tracking:
#
# * Number of iterations
# * Number of f_calls
# * Number of g_calls
# * Memory requirements (TODO)
#
####
# unconstrained
options = Optim.Options()
cutest_problems = CUTEst.select(max_var=10,contype = :unc, custom_filter=x->x["derivative_order"]>=2)
n = length(default_solvers)
m = length(cutest_problems)
f = open(join([version_dir, "cutest_benchmark.csv"], "/"), "w")
write(f, join(["Problem", "Optimizer", "Converged", "Time", "Minimum", "Iterations", "f_calls", "g_calls", "h_calls", "f_hat", "f_error", "x_error"], ","))
write(f, "\n")
@showprogress 1 "Benchmarking..." for prob in cutest_problems
    @show prob
    output = []

    nlp = CUTEstModel(prob)
    x0 = initial_x(nlp)

    x_hat = copy(x0)
    f_hat = obj(nlp, x_hat)
    xs = []
    times = []
        for i = 1:n
        	try
                d, constraints = optim_problem(nlp)
        		result = optimize(d, default_solvers[i], options)
                println(result)
                mintime = @elapsed optimize(d, x0, default_solvers[i], options)
                if f_hat > Optim.minimum(result)
                    f_hat = Optim.minimum(result)
                    x_hat[:] = Optim.minimizer(result)
                end
                push!(output, [prob,
                               default_names[i],
                               Optim.converged(result),
                               mintime,
                               Optim.minimum(result),
                               Optim.iterations(result),
                               Optim.f_calls(result),
                               Optim.g_calls(result),
                               Optim.h_calls(result)])
                push!(xs, Optim.minimizer(result))
            catch
                push!(output, ([prob,
                               default_names[i],
                               false,
                               Inf,
                               Inf,
                               Inf,
                               Inf,
                               Inf,
                               Inf]))
                push!(xs, fill(Inf, length(x0)))
        	end
        end
        for i = 1:n
            write(f, join(map(x->"$x",output[i]),",")*",")
            write(f, join(map(x->"$x",[f_hat, output[i][5]-f_hat, norm(xs[i]-x_hat, Inf)]),","))
            write(f, "\n")
        end
    finalize(nlp)
end
close(f)


text:
struct EquivalentCircuit
    Circuit::String
    Parameters::NamedTuple
end

Base.display(circuit::EquivalentCircuit) = println(circuit.Circuit)

function parametertuple(circuit,parameters)
    elements = foldl(replace,["["=>"","]"=>"","-"=>"",","=>""],init = denumber_circuit(circuit))
    name_strings = []
    number = 1
    for e in elements
        if e == 'P'
            push!(name_strings,"P$(number)w")
            push!(name_strings,"P$(number)n")
        else
            push!(name_strings,"$(e)$(number)")
        end
        number += 1
    end
    name_symbols = Tuple(map(Symbol, name_strings))
    return NamedTuple{name_symbols}(Tuple(parameters))
end

function initializecircuit(head=8,terminals="RCLP")
    karva = generatekarva(head,terminals)
    parameters = karva_parameters(karva)
    return Circuit(karva,parameters,nothing)
end

function initializepopulation(size=20,head=8,terminals="RCLP")
    return [initializecircuit(head,terminals) for i in 1:size]
end

function initializevariedpopulation(size=30,head=8)
    subpopulation_size = floor(Int(size/3))
    RCs = initializepopulation(subpopulation_size,head,"RC")
    RCLs = initializepopulation(subpopulation_size,head,"RCL")
    RCLPs = initializepopulation(subpopulation_size,head,"RCLP")
    return vcat(RCs,RCLs,RCLPs)
end

function simplifypopulation!(population,terminals="RCPL")
    for circuit in population
        simplifycircuit!(circuit,terminals)
    end
end

function circuit_offspring(circuit_1,circuit_2,terminals = "RCLP")
    offspring = ""
    rand_ = rand()
    if rand_ > 0.5
        offspring = one_point_crossover(circuit_1,circuit_2)
    elseif rand_ > 0.1
        offspring = two_point_crossover(circuit_1,circuit_2)
    else
        offspring = rand([circuit_1,circuit_2])
    end
    return mutate(offspring,terminals)
end

function circuitfitness(circuit,measurements,frequencies)
    tree = get_circuit_tree(circuit)
    circfunc,params,upper,param_inds = func_and_params_for_optim(tree)
    objective = objectivefunction(circfunc,measurements,frequencies)
    optparams,fitness = optimizeparameters(objective,params,upper)
    optparams = max.(min.(optparams,upper),0) #workaround of optim.jl's bug
    return deflatten_parameters(optparams,tree,param_inds), fitness, param_inds
end

function generate_offspring(population,terminals="RCLP")
    population_size = length(population)
    elite_size = Int(ceil(0.1*population_size))
    mating_pool_size = population_size-elite_size
    progenitors = tournamentselection(population,mating_pool_size,elite_size)
    elites = population[1:elite_size]
    offspring = Array{Circuit}(undef, mating_pool_size)
    for e in 1:mating_pool_size
        offspring[e] = circuit_offspring(progenitors[e],progenitors[mating_pool_size-e+1],terminals)
    end
    vcat(elites,offspring)
end

function evaluate_fitness!(population,measurements,frequencies)
    params = []
    param_inds = []
    for circuit in population
        params,circuit.fitness,param_inds = circuitfitness(circuit,measurements,frequencies)
        circuit.parameters[param_inds] = params
    end
end

"""
circuitevolution(measurements::Array{Complex{Float64},1},frequencies::Array{Float64,1}; <keyword arguments>)

Identify an equivalent electrical circuit that fits a given set of electrochemical impedance spectroscopy measurements.

The inputs are a circuit (e.g. "R1-[C2,R3]-P4") an array of complex-valued impedance measurements and their corresponding
frequencies. The output is a EquivalentCircuit object containing a field Circuit, which is a string denoting the identified circuit
and a field Parameters, which is a NamedTuple of the circuit's components with their corresponding parameter values.

# Arguments
- `generations::Integer=10`: the maximum number of iterations of the evolutionary algorithm.
- `population_size::Integer=30`: the number of individuals in the population during each iteration.
- `terminals::String="RCLP"`: the circuit components that are to be included in the circuit identification.
- `head::Integer=8`: a hyperparameter than controls the maximum considered complexity of the circuits.
- `cutoff::Float64=0.8`: a hyperparameter that controls the circuit complexity by removing redundant components.
- `initial_population::Array{Circuit,1}=nothing`:the option to provide an initial population of circuits
(obtained by using the loadpopulation function) with which the algorithm starts.

# Example
```julia

julia> using EquivalentCircuits, Random

julia> Random.seed!(25);

julia> measurements = [5919.9 - 15.7, 5918.1 - 67.5im, 5887.1 - 285.7im, 5428.9 - 997.1im, 3871.8 - 978.9im,
3442.9 - 315.4im, 3405.5 - 242.5im, 3249.6 - 742.0im, 1779.4 - 1698.9im,  208.2 - 777.6im, 65.9 - 392.5im];

julia> frequencies = [0.10, 0.43, 1.83, 7.85, 33.60, 143.84, 615.85,  2636.65, 11288.38, 48329.30, 100000.00];

julia> circuitevolution(measurements, frequencies , generations= 15, terminals = "RC")
[R1,C2]-[C3,R4]-R5
```
"""
function circuitevolution(measurements,frequencies;generations::Real=10,population_size=30,terminals = "RCLP",head=8,cutoff=0.8,initial_population=nothing)
    # Either initialize a new population, or work with a provided initial population.
    if isnothing(initial_population)
        population = initializepopulation(population_size,head,terminals) #initializevariedpopulation(population_size,head)
    else
        population = initial_population
    end
    # Theoretical simplification of the initial population.
        simplifypopulation!(population,terminals)
        evaluate_fitness!(population,measurements,frequencies)
        sort!(population)
        generation = 0
        convergence_threshold = 5e-4
        # Keep track of the fittest individual circuit.
        elite = minimum(population)
        min_fitness = elite.fitness
        while (min_fitness > convergence_threshold) && (generation<=generations)
            population = generate_offspring(population,terminals)
            simplifypopulation!(population,terminals)
            evaluate_fitness!(population,measurements,frequencies)
            sort!(population)
            elite = minimum(population).fitness < elite.fitness ? minimum(population) : elite
            min_fitness = elite.fitness
            population[1] = redundacy_testing(population[1],measurements,frequencies,terminals,cutoff)
            generation += 1
        end
        replace_redundant_cpes!(population)
        population = removeduplicates(sort!(vcat(population,elite)))
        for i in 1:3
            population[i] = removeredundancy(population[i],measurements,frequencies,terminals,cutoff)
        end
        # extract the converged circuits.
    population = filter(p -> p.fitness ≤ convergence_threshold, population)
    # in case of no converged circuits => alternate output print statement "Algorithm did not converge"
    if length(population) == 0
        println("Algorithm did not converge")
    else
        best_circuit = readablecircuit(population[1])
        return EquivalentCircuit(best_circuit,parameteroptimisation(best_circuit,measurements,frequencies)) #readablecircuit.(population[1:top_n])
    end
end

"""
    circuitevolution(filepath::String; <keyword arguments>)

 Identify an equivalent electrical circuit that fits a given set of electrochemical impedance spectroscopy measurements.

 The inputs are a circuit (e.g. "R1-[C2,R3]-P4") and a filepath to a CSV file containing the three following columns:
 the real part of the impedance, the imaginary part of the impedance, and the frequencies corresponding to the measurements.
 The output is NamedTuple of the circuit's components with their corresponding parameter values.

 # Arguments
- `filepath::String`: filepath of the CSV file containing the impedance measurements.
- `generations::Integer=10`: the maximum number of iterations of the evolutionary algorithm.
- `population_size::Integer=30`: the number of individuals in the population during each iteration.
- `terminals::String="RCLP"`: the circuit components that are to be included in the circuit identification.
- `head::Integer=8`: a hyperparameter than controls the maximum considered complexity of the circuits.
- `cutoff::Float64=0.8`: a hyperparameter that controls the circuit complexity by removing redundant components.
- `initial_population::Array{Circuit,1}=nothing`:the option to provide an initial population of circuits
(obtained by using the loadpopulation function) with which the algorithm starts.

"""

function circuitevolution(filepath::String;generations::Real=10,population_size=30,terminals = "RCLP",head=8,cutoff=0.8,initial_population = nothing)
    # Read the measurement file.
    meansurement_file = readdlm(filepath,',')
    # convert the measurement data into usable format.
    reals = meansurement_file[:,1]
    imags = meansurement_file[:,2]
    frequencies = meansurement_file[:,3]
    measurements = reals + imags*im
    # Either initialize a new population, or work with a provided initial population.
    if isnothing(initial_population)
        population = initializepopulation(population_size,head,terminals) #initializevariedpopulation(population_size,head)
    else
        population = initial_population
    end
    # Theoretical simplification of the initial population.
    simplifypopulation!(population)
    # Calculate the fitness for each individual in the population and sort the population by fitness.
    evaluate_fitness!(population,measurements,frequencies)
    sort!(population)
    generation = 0
    convergence_threshold = 5e-4
    # Keep track of the fittest individual circuit.
    elite = minimum(population)
    min_fitness = elite.fitness
    # apply genetic operators and tournament selection to obtain new generations of circuits
    # as long as none of the termination criteria are satisfied.
    while (min_fitness > convergence_threshold) && (generation<=generations)
        population = generate_offspring(population,terminals)
        simplifypopulation!(population)
        evaluate_fitness!(population,measurements,frequencies)
        sort!(population)
        elite = minimum(population).fitness < elite.fitness ? minimum(population) : elite
        min_fitness = elite.fitness
        population[1] = redundacy_testing(population[1],measurements,frequencies,terminals,cutoff)
        generation += 1
    end
    # Convert CPEs with second parameter equal to 0 or 1 to capacitor or resistor.
    replace_redundant_cpes!(population)
    population = removeduplicates(sort!(vcat(population,elite)))
    for i in 1:3
        population[i] = removeredundancy(population[i],measurements,frequencies,terminals,cutoff)
    end
    # extract the converged circuits.
    population = filter(p -> p.fitness ≤ convergence_threshold, population)
    best_circuit = readablecircuit(population[1])
    # in case of no converged circuits => alternate output print statement "Algorithm did not converge"
    if length(population) == 0
        println("Algorithm did not converge")
    else
        return EquivalentCircuit(best_circuit,parameteroptimisation(best_circuit,measurements,frequencies)) #readablecircuit.(population[1:top_n])
    end
end


text:
using Arrowhead
using LinearAlgebra
using Polynomials
using Test

# There are four tests, a random matrix test and three tests from the arrowhead paper (see README for details)
println("There are four tests for arrowhead matrices, a random matrix test and three")
println("tests from the arrowhead paper [1] (see README for details)")

tols = [1e2,1e2,1e2,1e2,1e2] # set the tolerances for eig
import Random
Random.seed!(421)

println("   1st test - Random SymArrow matrix")
A = GenSymArrow( 10, 10 )
E, info = eigen( A, tols ) # arrowhead solver
Residual = A * E.vectors - E.vectors * Diagonal(E.values) # compute the residual
normR = norm( Residual )
@test normR < 100.0 * eps() # norm of residual should be small

println("   2nd test - Example 1, p. 15 from [1]")
A = SymArrow( [ 2e-3, 1e-7, 0.0, -1e-7, -2e-3 ], [ 1e7, 1e7, 1.0, 1e7, 1e7 ], 1e20, 6 )
E, info = eigen( A, tols )
Λ=E.values
Λₜ=[ 1.0e20, 0.001999001249000113, 4.987562099722814e-9, -9.99999999998e-21, -2.0049855621017174e-6, -0.002001001251000111 ] # accurate eigenvalues from [1]
@test maximum( abs.( Λ - Λₜ ) ./ abs.( Λₜ ) ) < 3.0 * eps() # relative errors in eigenvalues should be small
v4 = [ 4.999999999985e-11, 9.999999999969e-7, 0.9999999999989999, -9.999999999970999e-7, -4.999999999985e-11, -9.99999999997e-21] # v_4 from [1]
@test maximum( abs.( E.vectors[:,4] - v4 ) ./ abs.( v4 ) ) < 3.0 * eps() # component-wise relative errors in eigenvector(s) should be small

println("   3rd test - Example 2, p. 16 from [1]")
A = SymArrow( [ 1+4eps(), 1+3eps(), 1+2eps(), 1+eps() ], [ 1.0, 2.0, 3.0, 4.0 ], 0.0, 5 )
E, info = eigen( A, tols )
Λ=E.values
Λₜ = [ 6.000000000000001, 1+4eps(), 1+3eps(), 1+2eps(), -4.999999999999999 ] # accurate eigenvalues from [1]
@test maximum( abs.( Λ - Λₜ ) ./ abs.( Λₜ ) ) < 3.0 * eps() # relative errors in eigenvalues should be small

println("   4th test - Example 3, p. 17 from [1]")
A = SymArrow( [ 1e10+1.0/3.0, 4.0, 3.0, 2.0, 1.0 ], [ 1e10 - 1.0/3.0, 1.0, 1.0, 1.0, 1.0 ], 1e10, 6 )
E, info = eigen( A, tols )
Λ=E.values
Λₜ = [ 1.999999999983333e10, 4.174722501468362, 3.188318635336404 , 2.223251566590035, 1.261850509234366, -0.3481422590562395 ] # accurate eigenvalues from [1]
@test maximum( abs.( Λ - Λₜ ) ./ abs.( Λₜ ) ) < 3.0 * eps() # relative errors in eigenvalues should be small

# There are four tests for DPR1, a random matrix test and three tests from the arrowhead paper (see README for details)
println("\n","There are four tests for DPR1 matrices, a random matrix test and three")
println("tests from the DPR1 paper [2] (see README for details)")

tols = [1e3,1e3,1e3,1e3,1e3] # set the tolerances for eig

println("   1st test - Random SymDPR1 matrix")
A = GenSymDPR1( 10 )
E, info = eigen( A, tols ) # DPR1 solver
Residual = A * E.vectors - E.vectors * Diagonal(E.values) # compute the residual
normR = norm( Residual )
@test normR < 300.0 * eps() # norm of residual should be small

println("   2nd test - Example 1, p. 9 from [2]")
A = SymDPR1( [ 1e10, 5.0, 4e-3, 0.0, -4e-3,-5.0 ], [ 1e10, 1.0, 1.0, 1e-7, 1.0,1.0 ], 1.0 )
E, info = eigen( A, tols ) # DPR1 solver
Λ=E.values
U=E.vectors
Λₜ=[ 1.000000000100000e20, 5.000000000100000, 4.000000100000001e-3, 9.999999998999997e-25, -3.999999900000001e-3, -4.999999999900000] # accurate eigenvalues from [1]
@test maximum( abs.( Λₜ - Λ ) ./ abs.( Λ ) ) < 3.0 * eps() # relative errors in eigenvalues should be small
v4 = [ 9.999999998999996e-18, 1.999999999800000e-18, 2.499999999749999e-15, -1.000000000000000, -2.499999999749999e-15, -1.999999999800000e-18] # v_4 from [1]
@test maximum( abs.( U[:,4] - v4 ) ./ abs.( v4 ) ) < 3.0 * eps() # component-wise relative errors in eigenvector(s) should be small

println("   3rd test - Example 2, p. 10 from [2]")
A = SymDPR1( [ 1+40eps(), 1+30eps(), 1+20eps(), 1+10eps() ], [ 1.0, 2.0, 2.0, 1.0 ], 1.0 )
E, info = eigen( A, tols ) # DPR1 solver
Λ=E.values
Λₜ= [ 11.0+32eps(), 1.0+39eps(), 1.0+25eps(), 1.0+11eps()] # accurate eigenvalues from [1]
@test maximum( abs.( Λ - Λₜ ) ./ abs.( Λₜ ) ) < 3.0 * eps() # relative errors in eigenvalues should be small

println("   4th test - Example 3, p. 11 from [2], see also Ming Gu's paper")
A = SymDPR1( [ 10.0/3.0, 2.0+1e-7, 2.0-1e-7, 1.0 ], [ 2.0, 1e-7, 1e-7, 2.0], 1.0 )
E, info = eigen( A, tols ) # DPR1 solver
U=E.vectors
@test norm(U'*U - I) < 3.0 * eps() # the eigenvectors are orthogonal (AND componentwise accurate)

println("\n","There is one test for SVD of HalfArrow matrices")
tols = [1e2,1e2,1e2,1e2,1e2]
println("   HalfArrow with entries varying in magnitude")
A=HalfArrow(sort(exp.(20*(rand(8).-0.5)),rev=true),(exp.(20*(rand(8).-0.5))))
S, info = svd( A, tols )
N1=norm(S.U'*S.U-I)
N2=norm(S.V*S.Vt-I)
@test N1 < 300.0*eps() && N2 < 300.0*eps() # Test the orthogonality of the eigenvectors
# @test norm(A*V-U*Diagonal(Σ))< 300.0*eps() # Test residual

println("\n","We test tridiagonal divide and conquer on the Wlikinson's matrix W21")
W=SymTridiagonal(abs.(collect(-10.0:10)),ones(20))
E=tdc(W)
Λₜ= [10.746194182903393,
 10.746194182903322,
 9.210678647361332,
 9.210678647304919,
 8.038941122829023,
 8.038941115814273,
 7.003952209528675,
 7.003951798616375,
 6.000234031584167,
 6.000217522257098,
 5.000244425001913,
 4.999782477742902,
 4.004354023440857,
 3.9960482013836245,
 3.0430992925788236,
 2.9610588841857264,
 2.1302092193625057,
 1.7893213526950813,
 0.9475343675292932,
 0.25380581709667804,
 -1.1254415221199847]
Λ=E.values
@test norm(sort(Λ)-sort(Λₜ))<5.0*eps()

println("\n","There are three tests for roots of polynomials")
println("   Example 1 from [3] - the Wilkinson's polynomial p18")
p18=Poly([   6402373705728000,
 -22376988058521600,
  34012249593822720,
 -30321254007719424,
  17950712280921504,
  -7551527592063024,
   2353125040549984,
   -557921681547048,
    102417740732658,
    -14710753408923,
      1661573386473,
      -147560703732,
        10246937272,
         -549789282,
           22323822,
            -662796,
              13566,
               -171,
                1])
D=roots(polyder(p18))
R,Qout=rootsah(p18,D)
@test norm(R-collect(18.0:-1:1.0))<5*eps()

println("   Example 2 from [3]")
e5=[-618970019642690000010608640,
4181389724724490601097907890741292883247104,
-6277101735386680066937501969125693243111159424202737451008,
713623846352979940529142984724747568191373312,
-20282409603651670423947251286016, 1]
p5=Poly(map(Float64,e5))
pd=polyder(p5)
pd=pd/pd[end]
companion = diagm(-1=>ones(3))
an = pd[end]
companion[1,:] = -pd[3:-1:0] / an
r = eigvals(companion)
D=rootsWDK(pd,r,1)
R,Qout=rootsah(p5,D)
Rtrue=  [2.028240960365167e31,
 1.759218623050247e13,
 1.7592185858329531e13,
 4.440892098500623e-16,
2.2204460492503136e-16]
@test norm(sort(R)-sort(Rtrue))<5.0*eps()

println("   Example 3 from [3] - the Chebyshev polynomial T[374]")
# Generate first 374 Chebyshev polynomials in BigFloat
n=374
T=Array{Any}(undef,n+1)
T[1]=Poly([BigFloat(1)])
T[2]=Poly([BigFloat(0),1])
for i=3:n+1
    T[i]=2*T[2]*T[i-1]-T[i-2]
end
# True zeros
Rtrue=map(Float64,[cos((2*k-1)*pi/(2*n)) for k=1:n])
# Zeros of T[n]
D=map(Float64,[cos((2*k-1)*pi/(2*(n-1))) for k=1:n-1])
# Computed zeros
R,qout=rootsah(T[n+1],D);
@test norm(sort(R)-sort(Rtrue))<10.0*eps()


text:
# real(conj(x) * y) avoiding computing the imaginary part if possible
@inline _realconjtimes(x, y) = real(conj(x) * y)
@inline _realconjtimes(x::Complex, y::Complex) = muladd(real(x), real(y), imag(x) * imag(y))
@inline _realconjtimes(x::Real, y::Complex) = x * real(y)
@inline _realconjtimes(x::Complex, y::Real) = real(x) * y
@inline _realconjtimes(x::Real, y::Real) = x * y

# imag(conj(x) * y) avoiding computing the real part if possible
@inline _imagconjtimes(x, y) = imag(conj(x) * y)
@inline function _imagconjtimes(x::Complex, y::Complex)
    return muladd(-imag(x), real(y), real(x) * imag(y))
end
@inline _imagconjtimes(x::Real, y::Complex) = x * imag(y)
@inline _imagconjtimes(x::Complex, y::Real) = -imag(x) * y
@inline _imagconjtimes(x::Real, y::Real) = Zero()


text:
import Base.Threads.@spawn

mutable struct Tasks
    d::Dict{Int, Any}
    c::Int # counter
end

const tasks = Tasks(Dict(), 1)

function getticket(tasks)
    return tasks.c += 1
end

function asyncall(f; kwargs...)
    global tasks
    ticket = getticket(tasks)
    id = @spawn f(; kwargs...)
    push!(tasks.d, ticket=>id)
    return (;ticket)
end

function checkjob(; ticket)
    global tasks
    id = tasks.d[ticket]
    done = istaskdone(id)
    if !done
        return (; done)
    else
        result = fetch(id)
        pop!(tasks.d, ticket)
        return (; done, result)
    end
end

asyncfunctions = (; checkjob)


text:
const UP = "👍"
const DOWN = "👎"
const client = D.BotClient(get(ENV, "DISCORD_TOKEN", ""))

@testset "A bunch of HTTP interactions" begin
    playback("responses.bson") do
        # Create a guild to play around in.
        guild = D.create_guild(client; name="MyGuild")
        user = D.get_user(client)
        try
            # Check that we can use either the object or the ID.
            @test D.get_guild(client, guild.id).id == guild.id
            @test D.get_guild(client, guild).id == guild.id

            # Loading channels.
            channels = D.get_guild_channels(client, guild)
            @test length(channels) == 4
            @test map(ch -> (ch.type, ch.name), channels) == [
                (D.ChannelType.GUILD_CATEGORY, "Text Channels"),
                (D.ChannelType.GUILD_CATEGORY, "Voice Channels"),
                (D.ChannelType.GUILD_TEXT, "general"),
                (D.ChannelType.GUILD_VOICE, "General"),
            ]

            # Reordering channels (a good test for ArrayBody).
            text_general = channels[3]
            channel = D.create_guild_channel(client, guild; name="test")
            D.update_guild_channel_positions(client, guild; positions=[
                (; id=text_general.id, position=channel.position),
                (; id=channel.id, position=text_general.position),
            ])
            @test D.get_channel(client, text_general.id).position == channel.position
            @test D.get_channel(client, channel.id).position == text_general.position

            # Update + delete channel.
            channel = D.update_channel(client, channel; name="test2")
            @test channel.name == "test2"
            D.delete_channel(client, channel)

            # Create + get messages.
            channel = text_general
            message = D.create_message(client, channel; content="hello")
            @test message.content == "hello"
            @test message.timestamp isa DateTime
            @test D.get_channel_message(client, channel, message).id == message.id
            messages = D.get_channel_messages(client, channel)
            @test length(messages) == 1 && messages[1].id == message.id
            @test messages[1].timestamp isa DateTime

            # Creating, listing, deleting reactions.
            D.create_reaction(client, channel, message, UP)
            users = D.get_reactions(client, channel, message, UP)
            @test length(users) == 1 && users[1].id == user.id
            D.delete_reaction(client, channel, message, UP)
            D.create_reaction(client, channel, message, UP)
            D.create_reaction(client, channel, message, DOWN)
            D.delete_all_reactions(client, channel, message, UP)
            @test isempty(D.get_reactions(client, channel, message, UP))
            @test !isempty(D.get_reactions(client, channel, message, DOWN))
            D.delete_all_reactions(client, channel, message)
            @test isempty(D.get_reactions(client, channel, message, DOWN))

            # Updating + deleting messages.
            @test D.update_message(client, channel, message; content="a").content == "a"
            D.delete_message(client, channel, message)
            @test isempty(D.get_channel_messages(client, channel))
            ids = map(i -> D.create_message(client, channel; content="hi").id, 1:5)
            @test length(D.get_channel_messages(client, channel)) == 5
            D.delete_messages(client, channel; messages=ids)
            @test isempty(D.get_channel_messages(client, channel))

            # Listing, creating, deleting invites.
            @test isempty(D.get_channel_invites(client, channel))
            @test isempty(D.get_guild_invites(client, guild))
            invite = D.create_channel_invite(client, channel)
            invites = D.get_channel_invites(client, channel)
            @test length(invites) == 1 && invites[1].code == invite.code
            invites = D.get_guild_invites(client, guild)
            @test length(invites) == 1 && invites[1].code == invite.code
            @test D.get_invite(client, invite).code == invite.code
            @test D.delete_invite(client, invite).code == invite.code
            @test isempty(D.get_channel_invites(client, channel))
            @test isempty(D.get_guild_invites(client, guild))

            # Pins + typing.
            D.create_typing_indicator(client, channel)
            @test isempty(D.get_pinned_messages(client, channel))
            message = D.create_message(client, channel; content="pin")
            D.create_pinned_channel_message(client, channel, message)
            pins = D.get_pinned_messages(client, channel)
            @test length(pins) == 1 && pins[1].id == message.id
            D.delete_pinned_channel_message(client, channel, message)
            @test isempty(D.get_pinned_messages(client, channel))

            # Emojis.
            @test isempty(D.get_guild_emojis(client, guild))
            # TODO: Figure out how to properly send the payload.
            # png = read(joinpath(@__DIR__, "emoji.png"), String)
            # emoji = D.create_guild_emoji(client, guild; name="test", image=png)
            # emojis = D.get_guild_emojis(client, guild)
            # @test length(emojis) == 1 && emojis[1].id == emoji.id
            # @test D.update_guild_emoji(client, guild, emoji; name="test2").name == "test2"
            # @test D.get_guild_emoji(client, guild, emoji).name == "test2"
            # D.delete_guild_emoji(client, guild, emoji)
            # @test isempty(D.get_guild_emojis(client, guild))

            # Guilds/members.
            @test D.get_guild(client, guild).id == guild.id
            @test D.update_guild(client, guild; name="test2").name == "test2"
            @test D.get_guild_member(client, guild, user).user.id == user.id
            members = D.get_guild_members(client, guild)
            @test length(members) == 1 && members[1].user.id == user.id
            @test D.modify_user_nick(client, guild; nick="test3").nick == "test3"
            @test D.get_guild_member(client, guild, user).nick == "test3"

            # Roles.
            roles = D.get_guild_roles(client, guild)
            @test length(roles) == 1 && roles[1].name == "@everyone"
            role = D.create_guild_role(client, guild; name="test")
            @test role.name == "test"
            roles = D.get_guild_roles(client, guild)
            @test length(roles) == 2 && roles[2].id == role.id
            @test D.update_guild_role(client, guild, role; name="test2").name == "test2"
            D.create_guild_member_role(client, guild, user, role)
            role_ids = D.get_guild_member(client, guild, user).roles
            @test length(role_ids) == 1 && role_ids[1] == role.id
            D.delete_guild_member_role(client, guild, user, role)
            @test isempty(D.get_guild_member(client, guild, user).roles)
            D.delete_guild_role(client, guild, role)
            roles = D.get_guild_roles(client, guild)
            @test length(roles) == 1 && roles[1].name == "@everyone"

            # Webhooks.
            @test isempty(D.get_channel_webhooks(client, channel))
            @test isempty(D.get_guild_webhooks(client, guild))
            webhook = D.create_webhook(client, channel; name="test")
            @test webhook.name == "test"
            @test D.get_webhook(client, webhook).id == webhook.id
            webhooks = D.get_channel_webhooks(client, channel)
            @test length(webhooks) == 1 && webhooks[1].id == webhook.id
            webhooks = D.get_guild_webhooks(client, guild)
            @test length(webhooks) == 1 && webhooks[1].id == webhook.id
            @test D.update_webhook(client, webhook; name="test2").name == "test2"
            message = D.execute_webhook(client, webhook, webhook.token; content="hi")
            @test message.content == "hi"
            D.delete_webhook(client, webhook)
            @test isempty(D.get_channel_webhooks(client, channel))
            @test isempty(D.get_guild_webhooks(client, guild))

            # Misc.
            @test eltype(D.get_guild_voice_regions(client, guild)) === D.VoiceRegion
            @test eltype(D.get_voice_regions(client)) === D.VoiceRegion
            @test isempty(D.get_guild_bans(client, guild))
            @test isempty(D.get_guild_integrations(client, guild))
            @test D.get_guild_prune_count(client, guild).pruned == 0
            @test !D.get_guild_widget(client, guild).enabled
        finally
            D.delete_guild(client, guild)
        end
    end
end

@testset "Sending files" begin
    playback("files.bson") do
        guild = D.create_guild(client; name="MyGuild")
        try
            channel = D.create_guild_channel(client, guild; name="test")
            msg = mktempdir() do dir
                file = joinpath(dir, "foo.json")
                write(file, "{}")
                open(file) do f
                    D.create_message(client, channel; file=f, __boundary__="abcdef")
                end
            end
            attachments = msg.attachments
            @test length(attachments) == 1
            attachment = attachments[1]
            @test attachment.filename == "foo.json"
            @test attachment.size == 2
        finally
            D.delete_guild(client, guild)
        end
    end
end


text:
# Load Data
path = "X://scJulia/data/10X_PBMC/filtered_feature_bc_matrix/"
TestMat = read_10X("S://PissaroData/10X_PBMC/filtered_feature_bc_matrix/");
B = deepcopy(TestMat["Gene Expression"]);

@time NormalizeData(B.array)

@time pissaro = init_pissaro(B, min_features = 0, min_cells = 5);
@time NormalizeData!(pissaro); #0.17s
@time lul = ScaleData(pissaro, assay = "RNA") #0.6
@time ScaleData!(pissaro, assay = "RNA") #0.6
dimnames(pissaro.assays["RNA"].norm_count)
intersect(names(pissaro.assays["RNA"].norm_count)[1], ["BPREB3", "CD79A", "CD79B", "DNTT"])
peepee = names(lul.assays["RNA"].scaled_count)
pissaro.assays["RNA"].norm_count[["CD79A", "CD79B", "DNTT"],:]
lul.assays["RNA"].scaled_count
pissaro.assays["RNA"].scaled_count
scale_factor = 10000
@time lul, dick, shit = tsvd(pissaro.assays["RNA"].scaled_count, 50)
ScaleData(pissaro.assays["RNA"].norm_count.array)
Base.summarysize(flunch)/(1024^2)
GetNormCount(flunch)
@time scaled_count(flunch)
FindVariableFeatures(flunch)
GetVariableFeatures(flunch)
GetFeatureData(flunch)
PercentageFeatureSet(pissaro)
pissaro.assays["RNA"].scaled_count["CD79A",:]
flunch.assays["RNA"].norm_count.array
A = convert(SparseMatrixCSC,B)
k = 10
m, n = size(A)
F = Arpack.svds(A, nsv = k)[1]
U::Array{Float32, 2}, Z::Array{Float32, 1}, V::Array{Float32, 2} = F.U, F.S, F.V
U, Z, V = convert(Array{Float32, 2}, F.U), convert(Array{Float32, 1}, F.S), convert(Array{Float32, 2}, F.V)
@assert size(U, 2) == length(Z) == size(V, 2) == k
return U, Z, V

B.array

lol = tsvd(B, 6)
lol[1]
using Plots
plot(lol[1][:,1], lol[1][:,2], seriestype = :scatter)

using MultivariateStats
MultivariateStats.ica
ok = fit(ICA, A, 5)

ScaledData = (C .- mean(C,dims=2)) ./ std(C,dims=2)

A1 = [[1]]
A2 = copy(A1)
A1[1] === A2[1]                 # true
A3 = deepcopy(A1)
A1[1] === A3[1]

C = deepcopy(flunch.assays["RNA"].norm_count.array)
StatsBase.standardize(ZScoreTransform, C.array[1:100, :], dims = 2)

GetCount(flunch, assay = "DICK")
in("DICK",keys(flunch.assays))

function test1(x::Int; y)
    x+y
end

Matrix(scaled_count(flunch))
scale_factor = 10000
GetCount(X, assay = assay)
GetNormCount(X, assay = assay)
NormalizeData(GetCount(flunch, assay = "RNA"), scale_factor = 2000)
X = flunch.assays["RNA"].count
log1p.(scale_factor * (X ./ sum(X, dims = 1)))
GetNormCount(flunch, assay = "RNA", named = false)

@time for i in ProgressBar(1:10000000) #wrap any iterator
    sqrt(i)
 end
@time sqrt.(1:(10^9))

function copy_cols(x::Vector{T}) where T
    inds = axes(x, 1)
    out = similar(Array{T}, inds, inds)
    for i = inds
        out[:, i] = x
    end
    return out
end

function copy_rows(x::Vector{T}) where T
    inds = axes(x, 1)
    out = similar(Array{T}, inds, inds)
    for i = inds
        out[i, :] = x
    end
    return out
end

function copy_col_row(x::Vector{T}) where T
    inds = axes(x, 1)
    out = similar(Array{Float64}, inds, inds)
    for col = inds, row = inds
        out[row, col] = x[row]
    end
    return out
end

function copy_row_col(x::Vector{T}) where T
    inds = axes(x, 1)
    out = similar(Array{T}, inds, inds)
    for row = inds, col = inds
        out[row, col] = x[col]
    end
    return out
end

x = randn(10000);

@time copy_cols(x)';
@time copy_col_row(x);
@time copy_rows(x);

fmt(f) = println(rpad(string(f)*": ", 14, ' '), @elapsed f(x))

map(fmt, [copy_cols, copy_rows, copy_col_row, copy_row_col]);

df = DataFrame(val=rand(8), injury = (repeat(["shoulder", "leg"],4)))

injuries["shoulder"]
A = feat
B = unique(A)
for i in B
    A[A .== i] = A[A .== i] .* append!([""],(string.(range(1,stop =length(A[A .== i])-1, step = 1)).* "."))
end

dict = Dict(i for i in A)
A = ["dog", "dog", "cat", "dog", "cat"]

function make_unique(x::Array{String,1})::Array{String,1}
    dict = Dict{String, Int}()
    for j in unique(x)
        dict[j] = 0
    end
    y = Array{Int,1}()
    for i in x
        push!(y, dict[i])
        dict[i] = dict[i]+1
    end
    x .*  replace("."  .* string.(y), ".0" => "")
end

mmfile = GZip.open(rmat_path)
@time s = deepcopy(readlines(fh));

open(path * "bite.mtx", "w") do io
    writedlm(io, readlines(fh))
end

open(path * "bite.mtx","w") do io
    write(io,s)
end

writedlm(path * "bite.mtx", readlines(fh))

@time lol = mmread(rmat_path);
@time lol2 = mmread2(rmat_path .* ".gz")
mmfile = GZip.open(rmat_path, "r")

function make_unique(x::Array{String,1}; ram::Bool = true)
    if(ram)
	    r=similar(x)
	    count=1
	    sp=sortperm(x)
	    r[sp[1]]=x[sp[1]]
	    for index in 2:length(x)
	    	if x[sp[index]]==x[sp[index-1]]
	    		r[sp[index]]=x[sp[index]]*"."*string(count)
	    		count+=1
	    	else
	    		r[sp[index]]=x[sp[index]]
	    		count=1
	    	end
	    end
	    return r
    else
        counts = Dict{String, Int}()
    result = Vector{String}()
    for e in x
        if e in keys(counts)
            push!(result, string(e, ".", counts[e]))
            counts[e] += 1
        else
            push!(result, e)
            counts[e] = 1
        end
    end
    return result
    end
end

pissaro.assays["RNA"].count.array =
assays["RNA"].scaled_count

span = 0.3

function FindVariableFeatures(mat::AbstractArray; selection_method::String = "vst", loess_span::Float64 = 0.3,
    clip_max::Union{String, Float64} = "auto", mean_function::Function = FastExpMean, dispersion_function::Function = FastLogVMR,
    num_bin::Int64 = 20, binning_method::String = "equal_width")::DataFrame
    mat = pissaro.assays["RNA"].count.array
    if clip_max == "auto"
            clip_max = sqrt(mat.n)
    end
    hvf_mean = mean(mat, dims = 2)[:,1]
    hvf_var  = var(mat, dims = 2)[:,1]
    hvf_expected = Vector{Float64}(undef, length(hvf_var))
    hvf_not_const = hvf_var .> 0
    res = loess(log10.(hvf_mean[hvf_not_const]), log10.(hvf_var[hvf_not_const]), span = 0.3)
    hvf_expected[hvf_not_const] .= 10 .^(Loess.predict(res, log10.(hvf_mean)))
    a, b, c = findnz(mat)
    nonzero_var =  min.((c .- hvf_mean[a]) ./ sqrt.(hvf_expected)[a], clip_max) .^2
    df = DataFrame(A = a, B = b, C = c, D = nonzero_var)
    df = sort(combine(DataFrames.groupby(df, :A), :D => sum), :A)
    hvf_std = (df.D_sum .+ zero_var)/(mat.n - 1)
    DataFrame(vst_mean = hvf_mean, vst_variance = hvf_var, vst_variance_expected = hvf_expected, vst_variance_standardized = hvf_std)
end

@time lul =  FindVariableFeatures(pissaro)
@time FindVariableFeatures!(pissaro)
FindVariableFeatures(B.assays["RNA"].count)
fieldnames(typeof(lul.assays["RNA"]))
(sort(pissaro.assays["RNA"].feature_data, :vst_variance_standardized, rev = true).features)[1:2000]
B = deepcopy(pissaro)
DataFrames.select(B.assays["RNA"].feature_data, DataFrames.Not(["vst_mean", "vst_variance", "vst_variance_expected", "vst_variance_standardized"]))
var_col = ["vst_mean", "vst_variance", "vst_variance_expected", "vst_variance_standardized"]

function our_type(number)
    add5() = number += 5
    double() = number *= 2
    get_number() = number
    (var) -> (add5;double;number; get_number)
end

bite = our_type(10)
bite.add5()
bite.get_number()

max_stars_repo_name: Cortalak/Pissaro

max_issues_repo_name: Cortalak/Pissaro

max_forks_repo_name: Cortalak/Pissaro


text:
module Baselines

using Random
using DataFrames
using StatsBase: mode

using AutoMLPipeline.Utils
using AutoMLPipeline.AbsTypes: Machine, Transformer, Learner, Workflow, Computer

import AutoMLPipeline.AbsTypes: fit!, transform!

export fit!,transform!
export Baseline, Identity

"""
    Baseline(
       default_args = Dict(
	       :name => "baseline",
          :output => :class,
          :strat => mode,
          :impl_args => Dict()
       )
    )

Baseline model that returns the mode during classification.
"""
mutable struct Baseline <: Learner
    name::String
    model::Dict
    args::Dict

    function Baseline(args=Dict())
        default_args = Dict(
            :name      => "baseline",
            :output    => :class,
            :strat     => mode,
            :impl_args => Dict()
			  )
		  cargs = nested_dict_merge(default_args, args)
		  cargs[:name] = cargs[:name]*"_"*randstring(3)
		  new(cargs[:name],Dict(),cargs)
    end
end

"""
    fit!(bsl::Baseline,x::DataFrame,y::Vector)

Get the mode of the training data.
"""
function fit!(bsl::Baseline,x::DataFrame,y::Vector)
  bsl.model = Dict(:choice => bsl.args[:strat](y))
end

"""
    transform!(bsl::Baseline,x::DataFrame)

Return the mode in classification.
"""
function transform!(bsl::Baseline,x::DataFrame)
  fill(bsl.model[:choice],size(x,1))
end

"""
    Identity(args=Dict())

Returns the input as output.
"""
mutable struct Identity <: Transformer
  name::String
  model::Dict
  args::Dict

  function Identity(args=Dict())
	 default_args = Dict{Symbol,Any}(
				:name => "identity",
				:impl_args => Dict()
			)
	 cargs = nested_dict_merge(default_args, args)
    cargs[:name] = cargs[:name]*"_"*randstring(3)
	 new(cargs[:name],Dict(),cargs)
  end
end

"""
    fit!(idy::Identity,x::DataFrame,y::Vector)

Does nothing.
"""
function fit!(idy::Identity,x::DataFrame,y::Vector)
    nothing
end

"""
    transform!(idy::Identity,x::DataFrame)

Return the input as output.
"""
function transform!(idy::Identity,x::DataFrame)
    return x
end

end


text:
################################################################################
#
#   DEFINITION OF GENERIC INTERFACE
#
################################################################################

# generic interface function
function getUnitcell(
            identifier :: Symbol,
            args...
        )

    # parse the identifier to a Val-type
    return getUnitcell(Val(identifier), args...)
end

# generic fallback
function getUnitcell(
            identifier :: Val{V},
            args...
        ) where {V}

    # throw an error
    error("Unitcell " * string(V) * " not implemented yet")
end

# EXPORT OF THE FUNCTION
export getUnitcell


text:
using SigmoidNumbers

if length(ARGS) < 2
  print("""
  USE:
  julia fft.jl [vector_length] [posit_precision] [posit_es] [distribution = {-1:0:1}]

  executes a fast fourier transform/inverse transform of a random vector and
  measures the deviation from the original in standard floats or posits.

  julia fft.jl 1024 32
  julia fft.jl 256  16 randn
  """)

  exit()
end

vector_size = parse(ARGS[1])
posit_precision = parse(ARGS[2])
posit_es = parse(ARGS[3])

F = posit_precision <= 16 ? Float16 : Float32
P = Posit{posit_precision, posit_es}

if length(ARGS) < 4
  randomizer() = rand(-1:1, vector_size)
elseif ARGS[4] == "randn"
  randomizer() = randn(vector_size)
elseif ARGS[4] == "rand"
  randomizer() = rand(vector_size)
end

function ffttest()

  sv = randomizer()

  v = F.(sv)
  t = ifft(fft(v))
  d = sum(abs.(Complex{Float64}.(t .- v)))
  println("$F deviation: ", d)

  vp = P.(sv)
  tp = ifft(fft(vp))
  dp = sum(abs.(Complex{Float64}.(tp .- vp)))
  println("$P deviation: ", dp)

  tp = ifft(fft(vp, normalize=:always), normalize=:always)
  dp = sum(abs.(Complex{Float64}.(tp .- vp)))
  println("$P deviation: ", dp)

  tp = ifft(fft(vp, normalize=:continuous), normalize=:continuous)
  dp = sum(abs.(Complex{Float64}.(tp .- vp)))
  println("$P deviation: ", dp)

end

ffttest()

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
#
# Dealing with negatives
#
# Handles warnings, errors, or ignoring negative values

immutable NegativeMode{T} end

const NegativeWarn = NegativeMode{:warn}()
const NegativeError = NegativeMode{:error}()
const NegativeIgnore = NegativeMode{:ignore}()

# I know, this is a really bad idea
global _hneg = NegativeIgnore

@inline function set_negative_mode(h::NegativeMode)
    global _hneg

    old_h, _hneg = _hneg, h

    return old_h
end

_check_negatives(ϕ::Factor) = _check_negatives(ϕ.potential)
_check_negatives(ϕ::Factor, h::NegativeMode) = _check_negatives(ϕ.potential, h)

_check_negatives{T<:Real}(p::Array{T}) = _check_negatives(p, _hneg)

_check_negatives(::Array, ::NegativeMode{:ignore}) = nothing

_check_negatives{T<:Real}(p::Array{T}, ::NegativeMode{:warn}) =
    any(p .< 0) && warn("potential has negative values")

# TODO find a better error type for this
_check_negatives{T<:Real}(p::Array{T}, ::NegativeMode{:error}) =
    any(p .< 0) && throw(ArgumentError("potential has negative values"))

_check_negatives(::Array, h::NegativeMode) = throw(ArgumentError("unkown handle type, " * repr(h)))


text:
module TestVectors

include("preamble.jl")

@testset "homogeneous size" begin
    xs = UnionVector(Any[1, 2.0])
    @test xs[1] === 1
    @test xs[2] === 2.0
    @test (xs[1] = 3.0) isa Any
    @test xs[1] === 3.0
end

@testset "heterogeneous size" begin
    xs = UnionVector(Any[UInt8(1), 2.0, (a=1, b=2)])
    @test xs[1] === UInt8(1)
    @test xs[2] === 2.0
    @test xs[3] === (a=1, b=2)
    @test (xs[1] = 3.0) isa Any
    @test xs[1] === 3.0
end

end  # module


text:
"""
    call_tetgen(surfaces, refinement)

Call Tetgen on surface geometry. A refinement is applied if `refinement < Inf`.
"""
function call_tetgen(surfaces, refinement)
    (; points, facets, facetmarkers, regions) = surfaces

    nregion = size(regions, 2)

    input = RawTetGenIO{Cdouble}()
    input.pointlist = points

    facetlist!(input, facets)

    input.facetmarkerlist = facetmarkers

    input.regionlist = [
        regions
        collect(1:nregion)'
        fill(0.1, 1, nregion)
    ]

    if isinf(refinement)
        tetgen = tetrahedralize(input, "pqA")
    else
        tetgen = tetrahedralize(input, "pqAa$(refinement)")
    end

    points = tetgen.pointlist
    facets = Int.(tetgen.trifacelist)
    facetmarkers = Int.(tetgen.trifacemarkerlist)
    elements = Int.(tetgen.tetrahedronlist)
    elementmarkers = Int.(tetgen.tetrahedronattributelist[:])

    (; points, facets, facetmarkers, elements, elementmarkers)
end


text:
function score(x, y)
    distance= sqrt(x^2+y^2)
    revenue=0
    if distance<=1
        revenue=10
    elseif distance<=5
        revenue=5
    elseif distance<=10
        revenue=1
    end
    return revenue
end

max_stars_repo_name: stepinski/julia-in

max_issues_repo_name: stepinski/julia-init

max_forks_repo_name: stepinski/julia-init


text:
function what(
    a::Int64,
    asdf::Float64)

end

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
using StructuralEquationModels, Distributions, Random, Optim, LineSearches
import StructuralEquationModels as SEM
include(
     joinpath(chop(dirname(pathof(StructuralEquationModels)), tail = 3),
     "test/examples/helper.jl")
     )

x = Symbol.("x", 1:13)

S = [:x1 0 0 0 0 0 0 0
     0 :x2 0 0 0 0 0 0
     0 0 :x3 0 0 0 0 0
     0 0 0 :x4 0 0 0 0
     0 0 0 0 :x5 0 0 0
     0 0 0 0 0 :x6 0 0
     0 0 0 0 0 0 :x7 :x9
     0 0 0 0 0 0 :x9 :x8]

F = [1.0 0 0 0 0 0 0 0
     0 1.0 0 0 0 0 0 0
     0 0 1.0 0 0 0 0 0
     0 0 0 1.0 0 0 0 0
     0 0 0 0 1.0 0 0 0
     0 0 0 0 0 1.0 0 0]

A = [0 0 0 0 0 0 1.0   0
     0 0 0 0 0 0 :x10 0
     0 0 0 0 0 0 :x11 0
     0 0 0 0 0 0 0    1.0
     0 0 0 0 0 0 0    :x12
     0 0 0 0 0 0 0    :x13
     0 0 0 0 0 0 0    0
     0 0 0 0 0 0 0    0]

ram_matrices = RAMMatrices(;A = A, S = S, F = F, parameters = x, colnames = nothing)

true_val = [repeat([1], 8)
            0.4
            repeat([0.8], 4)]

start = [repeat([1], 9)
             repeat([0.5], 4)]

imply_ml = RAMSymbolic(;specification = ram_matrices, start_val = start)

imply_ml.Σ_function(imply_ml.Σ, true_val)

true_dist = MultivariateNormal(imply_ml.Σ)

Random.seed!(1234)
x = transpose(rand(true_dist, 100000))
semobserved = SemObservedData(data = x, specification = nothing)

loss_ml = SemLoss(SEM.SemML(;observed = semobserved, n_par = length(start)))

optimizer =
    SemOptimizerOptim(
        BFGS(;linesearch = BackTracking(order=3), alphaguess = InitialHagerZhang()),# m = 100),
        Optim.Options(
            ;f_tol = 1e-10,
            x_tol = 1.5e-8))

model_ml = Sem(semobserved, imply_ml, loss_ml, optimizer)
objective!(model_ml, true_val)
solution_ml = sem_fit(model_ml)

@test isapprox(true_val, solution(solution_ml); atol = .05)


text:
# Using IterationControl.jl to control a particle tracker

using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using Plots
pyplot(size = (600, 300*(sqrt(5) -1)))

using IterationControl

# Consider an object that tracks a particle in one dimension, moving
# in discrete time at a speed proportional to the distance away from
# some moving target. We are viewing this particle as a toy example of
# online learning and call the constant of proportionality `η` a
# *learning rate*.

# The particle is initially at rest at the origin, which is also where
# the target begins.

# Calling `train!` on the particle moves it along for the specified
# number of time steps, while calling `ingest!` updates the target
# position ("data ingestion").

mutable struct Particle
    position::Float64
    target::Float64
    velocity::Float64
    η::Float64 # learning rate
    Particle(η) = new(0.0, 0.0, 0.0, η)
end

Particle(; η=0.1) = Particle(η)

loss(model::Particle) =  abs(model.target - model.position)

# A train method for one time step:

function train!(model::Particle)
    model.velocity = model.η*(model.target - model.position)
    model.position = model.position + model.velocity
end

# And for multiple time steps:

function train!(model, n)
    foreach(1:n) do _
        train!(model)
    end
    return nothing
end

# A method for updating the target:

function ingest!(model::Particle, target)
    model.target = target
    return nothing
end

# Lifting the `train!`, `loss` and `ingest!` methods:

IterationControl.train!(model::Particle, n) = train!(model, n)
IterationControl.loss(model::Particle) = loss(model)
IterationControl.ingest!(model::Particle, datum) = ingest!(model, datum)

# Lets `train` our particle by alternatively setting the target
# position to -1 or +1 every 5 time steps, for a total of 40 time
# steps (four -1/+1 cycles). We will then then allow our particle to "cruise"
# with the target fixed thereafter (the provided `data` iterator
# having been exhausted). This cruising is allowed to continue until
# the particle is within 0.01 of the target.

# We use a learning rate of 0.1.

model = Particle(0.1)
data = repeat([-1, 1], outer=4)

#-

target_positions = Float64[]
particle_positions = Float64[]

function callback!(model)
    push!(target_positions, model.target)
    push!(particle_positions, model.position)
end

IterationControl.train!(model,
                        Data(data),
                        Step(5),
                        Callback(callback!),
                        Threshold(0.01),
                        TimeLimit(0.005));

scatter(target_positions, label="target", ms=6)
scatter!(particle_positions, label="particle", ms=6)

#-

loss(model)

# Switching the target just once more to the -1 position and
# continuing training until the particle and target are withing 0.01:

data = [-1,]

IterationControl.train!(model,
                        Data(data),
                        Step(5),
                        Callback(callback!),
                        Threshold(0.01),
                        TimeLimit(0.005));

scatter(target_positions, label="target", ms=6)
scatter!(particle_positions, label="particle", ms=6)

using Literate #src
Literate.markdown(@__FILE__, @__DIR__, execute=false) #src
Literate.notebook(@__FILE__, @__DIR__, execute=true) #src


text:
# Proportional selection using the stohcastic acceptance algorithm of Lipowski and Lipowska

@doc """ propsel_funct()
Proportional selection applied to a population of attribute values.
Fitness is given by: function  fitness( x::Real64 )
"""
function propsel_funct( attributes::Vector{Float64}, fitness::Function, initial_value::Float64, fit_slope::Float64 )
  #println("propsel funct: initial value: ",initial_value,"  fit_slope: ",fit_slope)
  #println("  fit: ",fitness(attributes[1],initial_value,fit_slope))
  #println("propsel funct: pop: ",attributes)
  n = length(attributes)
  fitvec = map( a->fitness(a,initial_value,fit_slope), attributes )
  #println("fitvec: ",fitvec)
  fmax = maximum(fitvec)
  if fmax == 0.0
    # all elements have fitness zero
    println("Warning: all elements have fitness zero")
    return deepcopy(attributes)  # do not modify attributes
  end
  #new_attributes = zeros(Float64,n)
  fitvec /= fmax
  #println("fitvec: ",fitvec)
  selected = zeros(Int64, n)
  k = 0
  while k < n
    i = rand(1:n)
    w = fitvec[i]
    if rand() < w
      selected[k + 1] = i
      k += 1
    end
  end
  attributes[selected]
end

@doc """ test_propsel_funct()
  Does a simple test of the propsel_funct function.
  Will fail with an assertion error if the test fails.
"""
function test_propsel_funct()
  n = 10
  # create an attribute array [1.0, 2.0, . . . ,10.0]
  attr = [ convert(Float64,i) for i=1:n]
  # x->(x==5.0 ? :2.0 : 0.0)  is an anonymous function that gives only 5.0 fitness 2.0
  #   and all other numbers fitness 0.0
  new_attr = propsel_funct( attr, x->(x==5.0 ? :2.0 : 0.0) )
  @assert( new_attr == fill(5.0,n) )  # new_attr should be [5.0, 5.0, . . . ,5.0]
end


text:
using Sigma
# Approximately, what’s the probability that in a room filled with 23 people at least one pair of people have the same birthday?
birthdays = iid(Int, i->discreteuniform(1,366),30)
function same_birthday(birthdays)
  bday = false
  for i = 1:length(birthdays), j = 1:length(birthdays)
    if i != j
      bday |= (birthdays[i] == birthdays[j])
    end
  end
  bday
end

rand(same_birthday(birthdays))
prob(same_birthday(birthdays))


text:
mutable struct ShowSettings
    plain::Bool
    tableformat::TextFormat
    floatformat::String
    showfixed::Bool
    showevals::Bool
    border::Crayon
    header::Crayon
    subheader::Crayon
    fixed::Crayon
    error::Crayon
    highlighted::Crayon
    section::Crayon
    ShowSettings() = new(false, tf_unicode_rounded, "%9.4g", true, true,
                         crayon"light_blue", crayon"light_blue negative bold",
                         crayon"dark_gray bold", crayon"dark_gray",
                         crayon"light_red blink", crayon"negative", crayon"green bold")
end

const showsettings = ShowSettings()

function printtable(args...; formatters=(), hlines=:none, kw...)
    if showsettings.plain
        c = crayon"default"
        pretty_table(args...; formatters=formatters, alignment=:l, crop=:none, tf=tf_compact, hlines=hlines,
                     border_crayon=c, header_crayon=c, subheader_crayon=c,
                     highlighters=())
    else
        pretty_table(args...; formatters=formatters, alignment=:l, crop=:none, tf=showsettings.tableformat, hlines=hlines,
                     border_crayon=showsettings.border, header_crayon=showsettings.header, subheader_crayon=showsettings.subheader,
                     kw...)
    end
end

function section(io, args...)
    if showsettings.plain
        println(io, args...,)
    else
        println(io, showsettings.section, args..., crayon"default")
    end
end

function left(s::String, maxlen::Int)
    (length(s) <= maxlen)  &&  (return s)
    return s[1:maxlen]
end

show(io::IO, mime::MIME"text/plain", domain::AbstractDomain) = show(io, domain)
function show(io::IO, dom::AbstractDomain)
    section(io, string(typeof(dom)) * " (ndims: ", ndims(dom), ", length: ", length(dom), ")")
    hrule = Vector{Int}()
    push!(hrule, 0, 1, ndims(dom)+1)
    table = Matrix{Union{Int,Float64}}(undef, ndims(dom), 6)
    for i in 1:ndims(dom)
        if isa(dom, Domain{1})
            a = convert(Vector{Float64}, dom)
        else
            a = isa(dom, Domain)  ?  dom[i]  :  axis(dom, i)
        end
        b = 0
        if length(a) > 1
            b = a .- circshift(a, 1)
            b = b[2:end]
        end
        table[i, 1] = i
        table[i, 2] = length(a)
        table[i, 3:6] = [minimum(a), maximum(a), minimum(b), maximum(b)]
    end
    if isa(dom, Domain)  &&  (ndims(dom) >= 2) # treat Domain{1} as a Cartesian one, despite it is a linear one.
        printtable(io, table[:, 1:4], ["Dim", "Size", "Min val", "Max val"],
                   hlines=hrule, formatters=ft_printf(showsettings.floatformat, 3:4))
    else
        printtable(io, table, ["Dim", "Size", "Min val", "Max val", "Min step", "Max step"],
                   hlines=hrule, formatters=ft_printf(showsettings.floatformat, 3:6))
    end
end

#=
The following is needed since AbstractData <: AbstractArray:
without it the show(bstractArray) would be invoked.
=#
show(io::IO, mime::MIME"text/plain", data::AbstractData) = show(io, data)

function show(io::IO, data::AbstractData)
    section(io, typeof(data), ": (length: ", (length(data.val)), ")")
    table = Matrix{Union{String,Float64}}(undef, 0, 7)
    hrule = Vector{Int}()
    push!(hrule, 0, 1)

    names = fieldnames(typeof(data))
    error = Vector{Bool}()
    for name in names
        a = getfield(data, name)
        nan = length(findall(isnan.(a))) + length(findall(isinf.(a)))
        a = a[findall(isfinite.(a))]
        push!(error, nan > 0)
        table = vcat(table, [string(name) minimum(a) maximum(a) mean(a) median(a) std(a) (nan > 0  ?  string(nan)  :  "") ])
    end
    push!(hrule, 0, size(table)[1]+1)
    printtable(io, table, ["", "Min", "Max", "Mean", "Median", "Std. dev.", "Nan/Inf"],
               hlines=hrule, formatters=ft_printf(showsettings.floatformat, 2:6),
               highlighters=(Highlighter((data,i,j) -> error[i], showsettings.error)))
end

function show(io::IO, par::Parameter)
    if par.fixed
        println(io, "Value : ", par.val, "   (FIXED)")
    elseif isnan(par.unc)
        println(io, "Value : ", par.val, "  [", par.low , " : ", par.high, "]")
    else
        if par.val == par.patched
            println(io, "Value : ", par.val, " ± ", par.unc,  "  [", par.low , " : ", par.high, "] ")
        else
            println(io, "Value : ", par.val, " ± ", par.unc,  "  [", par.low , " : ", par.high, "], patched: " , par.patched)
        end
    end
end

function preparetable(comp::AbstractComponent; cname::String="?", cfixed=false)
    table = Matrix{Union{String,Float64}}(undef, 0, 7)
    fixed = Vector{Bool}()
    error = Vector{Bool}()
    watch = Vector{Bool}()

    ctype = split(string(typeof(comp)), ".")
    (ctype[1] == "GFit")  &&   (ctype = ctype[2:end])
    ctype = join(ctype, ".")

    for (pid, param) in getparams(comp)
        parname = string(pid.name)
        if pid.index >= 1
            parname *= "[" * string(pid.index) * "]"
        end
        parname *= (param.fixed  ?  " (FIXED)"  :  "")
        (!showsettings.showfixed)  &&  param.fixed  &&  continue
        range = strip(@sprintf("%7.2g:%-7.2g", param.low, param.high))
        (range == "-Inf:Inf")  &&  (range = "")
        table = vcat(table,
                     [cname * (cfixed  ?  " (FIXED)"  :  "") ctype parname range param.val param.unc param.patched])
        push!(fixed, param.fixed)
        push!(error, !(param.low <= param.val <= param.high))
        push!(watch, param.val != param.patched)
        if !showsettings.plain
            cname = ""
            ctype = ""
        end
    end
    if length(table) == 0
        table = vcat(table, [cname ctype "" NaN ""])
        push!(fixed, false)
        push!(error, false)
    end
    return (table, fixed, error, watch)
end

function show(io::IO, comp::AbstractComponent)
    (table, fixed, error, watch) = preparetable(comp)
    printtable(io, table, ["Component", "Type", "Param.", "Range", "Value", "Uncert.", "Patched"],
               formatters=ft_printf(showsettings.floatformat, 5:7),
               highlighters=(Highlighter((data,i,j) -> (watch[i]  &&  (j==7)), showsettings.highlighted),
                             Highlighter((data,i,j) -> (error[i]  && (j in (3,4,5))), showsettings.error),
                             Highlighter((data,i,j) -> fixed[i], showsettings.fixed)))
end

function show(io::IO, red::ExprReducer)
    println(io, red.ef.expr)
end

function show(io::IO, red::SumReducer)
    println(io, join(string.(red.list), " + "))
end

function show(io::IO, model::Model)
    println(io)
    section(io, "Components:")
    (length(model.cevals) == 0)  &&  (return nothing)

    table = Matrix{Union{String,Float64}}(undef, 0, 7)
    fixed = Vector{Bool}()
    error = Vector{Bool}()
    watch = Vector{Bool}()
    hrule = Vector{Int}()
    push!(hrule, 0, 1)
    for (cname, ceval) in model.cevals
        comp = ceval.comp
        (t, f, e, w) = preparetable(comp, cname=string(cname), cfixed=(ceval.cfixed >= 1))
        table = vcat(table, t)
        append!(fixed, f .| (ceval.cfixed >= 1))
        append!(error, e)
        append!(watch, w)
        push!(hrule, length(error)+1)
    end
    printtable(io, table, ["Component", "Type", "Param.", "Range", "Value", "Uncert.", "Patched"],
               hlines=hrule, formatters=ft_printf(showsettings.floatformat, 5:7),
               highlighters=(Highlighter((data,i,j) -> (watch[i]  &&  (j==7)), showsettings.highlighted),
                             Highlighter((data,i,j) -> (error[i] &&  (j in (3,4,5))), showsettings.error),
                             Highlighter((data,i,j) -> fixed[i], showsettings.fixed)))

    i = 1
    error = Vector{Bool}()
    table = Matrix{Union{String,Int,Float64}}(undef,
                                              length(model.cevals) + length(model.revals), 6)
    for (cname, ceval) in model.cevals
        result = ceval.buffer
        v = view(result, findall(isfinite.(result)))
        (length(v) == 0)  &&  (v = [NaN])
        nan = length(findall(isnan.(result)))
        inf = length(findall(isinf.(result)))
        table[i, 1] = string(cname)
        table[i, 2] = ceval.counter
        table[i, 3:5] = [minimum(v), maximum(v), mean(v)]
        table[i, 6] = (nan+inf > 0  ?  string(nan)  :  "")
        push!(error, (nan+inf > 0))
        i += 1
    end

    for (rname, reval) in model.revals
        result = reval.buffer
        v = view(result, findall(isfinite.(result)))
        (length(v) == 0)  &&  (v = [NaN])
        nan = length(findall(isnan.(result)))
        inf = length(findall(isinf.(result)))
        table[i, 1] = string(rname)
        table[i, 2] = reval.counter
        table[i, 3:5] = [minimum(v), maximum(v), mean(v)]
        table[i, 6] = (nan+inf > 0  ?  string(nan)  :  "")
        push!(error, (nan+inf > 0))
        i += 1
    end

    if length(model.patchfuncts) > 0
        section(io, "Patch expressions:")
        for pf in model.patchfuncts
            println(io, string(pf.expr))
        end
    end

    println(io)
    section(io, "Reducers:")
    for (rname, reval) in model.revals
        print(rname, ": ")
        show(io, reval.red)
        println(io)
    end

    if showsettings.showevals
        section(io, "Evaluations:")
        printtable(io, table, ["Component", "Eval. count", "Min", "Max", "Mean", "NaN/Inf"],
                   hlines=[0,1, length(model.cevals)+1,  length(model.cevals)+length(model.revals)+1],
                   formatters=ft_printf(showsettings.floatformat, 3:5),
                   highlighters=(Highlighter((data,i,j) -> (error[i] && j==5), showsettings.error)))
    end
end

function show(io::IO, multi::MultiModel)
    for id in 1:length(multi.models)
        section(io, "\n=====================================================================")
        section(io, "Model $id:")
        show(io, multi.models[id])
    end
    println(io)

    if length(multi.patchfuncts) > 0
        section(io, "Multi model patch expressions:")
        for pf in multi.patchfuncts
            println(io, string(pf.expr))
        end
    end
end

function show(io::IO, res::FitResult)
    section(io, "Fit results:")

    println(io, @sprintf("    #Data  : %8d              Fit-stat: %-10.5g", res.ndata, res.fitstat))
    println(io, @sprintf("    #Param : %8d              Red. GOF: %-10.4g", res.nfree, res.gofstat / res.dof))
    println(io, @sprintf("    DOF    : %8d              Prob.   : %-10.4g", res.dof, 10^res.log10testprob))

    print(io,            "    Status : ")
    (crayon, status, message) = as_string(res.mzer)
    if showsettings.plain
        print(io, @sprintf("%8s", status))
    else
        print(io, crayon, @sprintf("%8s", status), crayon"default")
    end

    println(io, @sprintf("              Elapsed : %-10.4g s", res.elapsed))

    if message != ""
        println(io, crayon, message, crayon"default")
    end
end