text:
#=
Type 'EventClass' for Event History objects
Subtypes:
Surv:      (Right)-censored survival outcome
SurvInt:   Interval-censored outcome (Left,Right,Interval)
SurvTrunc: Right-censoring + left-truncation
CompRisk:  Competing risks model

constructor: Event
Access methods: Time, Status, Cause, ...
=#

const CensNot = 0
const CensLeft = 1
const CensRight = 2
const CensInt = 3

#=
Types/classes: Surv,CompRisk,...
=#

@compat abstract type EventClass end
#abstract EventInt <: EventClass

immutable Surv <: EventClass
    Time::Number  # Exit time
    Status::Bool  # Censoring status
end

immutable SurvTrunc <: EventClass
    Entry::Number # Entry time
    Time::Number  # Exit time
    Status::Bool  # Censoring status
end

immutable CompRisk <: EventClass
    Entry::Number # Entry time
    Time::Number  # Exit time
    Status::Bool  # Censoring status
    Cause::Number
    function CompRisk(Entry,Time,Cause) new(Entry,Time,Cause!=0,Cause) end
end

immutable SurvInt <: EventClass
    Time::Number     # Time 1
    Time2::Number    # Interval [Time;Time2]. Use Inf/-Inf for left/right censoring
    Status::Int      # Censoring (0:none,1:right,2:left,3:interval)
    function SurvInt(Time,Time2)
        if Time>Time2 error("time 1 larger than time 2") end
        if (Time==Time2)
            status = EventHistory.CensNot
        elseif Time2==Inf
            status = EventHistory.CensRight
        elseif Time==-Inf
            status = EventHistory.CensLeft
        else
            status = EventHistory.CensInt
        end
        new(Time,Time2,status)
    end
end

const Events = Vector{EventClass}
const Vec = Union{Vector,Matrix}
const BoolVec = Union{Vector{Bool},Matrix{Bool},BitVector}

#=
show methods
=#

function show(io::IO, obj::EventClass)
    print(io, obj.Time, obj.Status>0 ? "":"+")
end

function show(io::IO, obj::SurvTrunc)
    print(io, "(", obj.Entry, ";", obj.Time, obj.Status>0 ? "":"+","]")
end

function show(io::IO, obj::CompRisk)
    print(io, "(", obj.Entry, ";", obj.Time, ":", obj.Status==0 ? "+" : obj.Cause,"]")
end

function show(io::IO, obj::SurvInt)
    if obj.Status==EventHistory.CensNot val=obj.Time
    elseif obj.Status==EventHistory.CensLeft val=string("(-Inf;",obj.Time2,"]")
    elseif obj.Status==EventHistory.CensRight val=string("[",Time,";Inf)")
    else val = string("[",obj.Time,";",obj.Time2,"]"); end
    print(io,val)
end

#=
Accessors
=#

# Meta-programming definitions of Time,Entry,Status,Cause access methods
for key = (:Time, :Entry, :Status, :Cause)
    for typ = (:Vector, :Matrix) #, :(DataFrames.DataVector))
        @eval function $(key){T<:EventClass}(e::$(typ){T})
            ##$(symbol(string("Event_",key)))(e)
            n = length(e)
            res = Array{typeof(e[1].$(key))}(n)
            for i=1:n
                res[i] = e[i].$(key)
            end
            res
        end
    end
end

export Time,Entry,Status,Cause

#=
Event constructor
=#

function Event(time::Vec,
               status::BoolVec)
    n = length(time)
    sz = size(time)
    if (length(sz)==1)
        E = Array{EventHistory.Surv}(n)
    else
        E = Array{EventHistory.Surv}(sz[1],sz[2])
    end
    for i=1:n
        E[i] = EventHistory.Surv(time[i],status[i])
        ## E[i] = EventHistory.Surv(time[i],convert(Bool,status[i]))
    end
    E
end

# function Event(time::Union{Vector,DataFrames.DataVector},
#                status::Union{Vector,DataFrames.DataVector},
#                method::Symbol=:comprisk)
#     Event(time,status,string(method))
# end

function Event(time::Vec,
               status::Vec,
               method::AbstractString)
    n = length(time)
    sz = size(time)
    ## interval
    if lowercase(method)=="interval"
        if (length(sz)==1)
            E = Array{EventHistory.SurvInt}(n)
        else
            E = Array{EventHistory.SurvInt}(sz[1],sz[2])
        end
        for i=1:n
            E[i] = EventHistory.SurvInt(time[i],status[i])
        end
        return E
    end
    ## comprisk:
    if (length(sz)==1)
        E = Array{EventHistory.CompRisk}(n)
    else
        E = Array{EventHistory.CompRisk}(sz[1],sz[2])
    end
    for i=1:n
        E[i] = EventHistory.CompRisk(0,time[i],status[i])
    end
    E
end

function Event(entry::Vec,
               time::Vec,
               status::BoolVec)
    n = length(time)
    sz = size(time)
    if (length(sz)==1)
        E = Array{EventHistory.SurvTrunc}(n)
    else
        E = Array{EventHistory.SurvTrunc}(sz[1],sz[2])
    end
    for i=1:n
        E[i] = EventHistory.SurvTrunc(entry[i],time[i],status[i])
    end
    E
end

function Event(entry::Vec,
               time::Vec,
               status::Vec)
    n = length(time)
    sz = size(time)
    # issurv = typeof(status[1])==Bool
    # if (issurv)
    #     E = Array{EventHistory.SurvTrunc}(sz[1],sz[2])
    #     for i=1:n
    #         E[i] = EventHistory.SurvTrunc(entry[i],time[i],status[i])
    #     end
    #     return E
    # end
    if (length(sz)==1)
        E = Array{EventHistory.CompRisk}(n)
    else
        E = Array{EventHistory.CompRisk}(sz[1],sz[2])
    end
    for i=1:n
        E[i] = EventHistory.CompRisk(entry[i],time[i],status[i])
    end
    E
end

function Event(var::Vector{Symbol}, data::DataFrame, censdef::Function=x->x.>0)
    p = size(var,1)
    if p==2
        Status = data[var[2]]
    else
        Status = data[var[3]]
    end
    Status = censdef(Status)
    if p==2
        return Event(data[var[1]],Status)
    end
    return Event(data[var[1]],data[var[2]],Status)
end

#=
Misc
=#

# function transpose(x::EventHistory.EventClass)
#     x
# end
# export transpose


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "YasmBuilder"
version = v"1.3.0"

# Collection of sources required to build YasmBuilder
sources = [
    "http://www.tortall.net/projects/yasm/releases/yasm-1.3.0.tar.gz" =>
    "3dce6601b495f5b3d45b59f7d2492a340ee7e84b5beca17e48f862502bd5603f",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd yasm-1.3.0/
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain
make -j${nproc}
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, :glibc),
    Linux(:x86_64, :musl)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libyasm", :libyasm),
    ExecutableProduct(prefix, "", :yasm),
    ExecutableProduct(prefix, "", :ytasm),
    LibraryProduct(prefix, "libyasmstd", :libyasmstd),
    ExecutableProduct(prefix, "", :vsyasm)
]

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

]

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

ext: jl

max_forks_repo_name: jpsamaroo/YasmBuilder


text:
warn("This script requires a large amount of memory to run!")

#BLAS.set_num_threads(6)

using TensorOperations
using JLD
using LinearMaps

include("src/PlotTools.jl")

include("src/EDutils.jl")
import EDutils
ED = EDutils

include("src/LatCFT.jl")
import LatCFT
LCFT = LatCFT

modelname = "ANNNI3";

ham_comp = (λ, hz, hx, Δ, Δd) -> begin
    h1, h2 = ED.ham_ising_comp(λ, hz, hx)
    h2 = h2 + Δd * kron(ED.pauliZ, ED.pauliZ)
    h3 = Δ * kron(ED.pauliX, eye(2), ED.pauliX)
    h1,h2,h3
end

ham = (λ, hz, hx, Δ, Δd) -> begin
    h1, h2, h3 = ham_comp(λ, hz, hx, Δ, Δd)
    h3 + 0.5(kron(h2, eye(2)) + kron(eye(2), h2)) + 1/3 * (kron(h1, eye(4)) + kron(eye(2), h1, eye(2)) + kron(eye(4), h1))
end

λ = 1.0; hz = 1.0; Δ = -0.5; Δd = -0.5 #Self-dual ANNNI model remains at Ising criticality for ferromagnetic interactions.
h = reshape(ham(λ, hz, 0.0, Δ, Δd), (2,2,2,2,2,2))
h_comp = ham_comp(λ, hz, 0.0, Δ, Δd)

nev = 80

eV_data = Dict{Int,Any}()
ed_data = Dict{Int,Any}()
Hs_data = Dict{Int,Any}()

datapath = "data/"

ED.ed!(eV_data, collect(6:1:8), h, nev)

ps_allowed = Float64[-4.0,-3.0,-2.0,-1.0,0.0,1.0,2.0,3.0,4.0]
ED.ed!(eV_data, collect(9:1:22), h, nev, symm_func=ED.get_Pmoms(ps_allowed))

save("$(datapath)$(modelname)_redP_eV.jld", "eV", eV_data)

for (N, eVd) in eV_data
    ed_data[N] = (eVd[1:2]...)
end

Hs_data = Dict{Int,Any}()
LCFT.Hns_ed!(Hs_data, eV_data, h_comp, 2, 2)

LCFT.normalise_Hn!(Hs_data, ed_data)

save("$(datapath)$(modelname)_redP_LCFTdata.jld", "ev", ed_data, "Hs", Hs_data)

max_stars_repo_name: amilsted/spinchain-virasoro-demo

max_issues_repo_name: amilsted/spinchain-virasoro-demo

max_forks_repo_name: amilsted/spinchain-virasoro-demo


text:
import Nemo.isone, Nemo.divexact, Base.copy
export divexact!, gcd_into!, coprime_base, coprime_base_insert

function divexact!(a::fmpz, b::fmpz)
  ccall((:fmpz_divexact, libflint), Nothing,
          (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), a, a, b)
  return a
end

function gcd_into!(a::fmpz, b::fmpz, c::fmpz)
  ccall((:fmpz_gcd, libflint), Nothing,
          (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), a, b, c)
  return a
end

function gcd_into!(a::T, b::T, c::T) where T <: Integer
  return gcd(b, c)::T
end

function mul_into!(a::T, b::T, c::T) where T
  return (b*c)::T
end

function mul_into!(a::fmpz, b::fmpz, c::fmpz)
  mul!(a, b, c)
  return a
end

function copy(a::fmpz)
  return deepcopy(a)
end

function copy_into!(a, b)
  return b
end

function copy_into!(a::MPolyElem, b::MPolyElem)
  return copy(b)
end

function copy_into!(a::fmpz, b::fmpz)
  ccall((:fmpz_set, libflint), Nothing, (Ref{fmpz}, Ref{fmpz}), a, b)
  return a
end

function copy_into!(a::fmpq, b::fmpq)
  ccall((:fmpq_set, libflint), Nothing, (Ref{fmpq}, Ref{fmpq}), a, b)
  return a
end

#for larger lists much better than Bill's (Nemo's) prod function
# the build-in Julia is better than Bill's Nemo function anyway
"""
  Data structure for a product-tree evaluation, need only
    O(nbits(length)) storage in comparison to O(n) for the julia
  or flint version.

  See implementation of `my_prod` for usage.
"""
mutable struct ProdEnv{T}
  level :: Vector{Int}
  val   :: Vector{T}
  last  :: Int

  function ProdEnv{T}(n::Int) where {T}
    r = new{T}()
    m = nbits(n)
    r.level = Array{Int}(undef, m)
    r.val   = Array{T}(undef, m)
    r.last = 0
    for i=1:m
      r.val[i] = T()
    end
    return r
  end
end

function Base.push!(A::ProdEnv{T}, b::T) where T
  if A.last == 0
    A.level[1] = 1
    A.val[1] = copy_into!(A.val[1], b)
    A.last = 1
    return
  end
  if A.level[A.last] > 1
    A.last += 1
    A.level[A.last] = 1
    A.val[A.last] = copy_into!(A.val[A.last], b)
    return
  end

  lev = 1
  while A.last > 0 && A.level[A.last] <= lev
    b = mul_into!(A.val[A.last], A.val[A.last], b)
    lev += A.level[A.last]
    A.last -= 1
  end
  A.last += 1
  A.level[A.last] = lev
  A.val[A.last] = b
  return
end

function finish(A::ProdEnv)
  b = A.val[A.last]
  while A.last >1
    A.last -= 1
    b = mul_into!(b, A.val[A.last], b)
  end
  return b
end

function my_prod(a::AbstractVector{T}) where T
  b = ProdEnv{T}(length(a))
  for x in a
    push!(b, x)
  end
  return finish(b)
end

"""
 A specialized version of the ProdEnv for 2x2 matrices representing
 polynomials

 f = sum a_i x^i

 can be computed via

    x   0      x   0    ...    x   0  =  x^(n+1   0)
    a_0 1      a_1 1           a_n 1     f(x)     1

 clearly, the second col is not neccessary, and (slightly more complicated)
 the only powers if x that are used in the prod tree are x^(2^i)
 for 0<= i <= nbits(n)
"""
mutable struct EvalEnv{T}
  level :: Vector{Int}
  pow   :: Vector{T}
  val   :: Vector{T}
  last  :: Int

  function EvalEnv{T}(x::T, n::Int) where {T}
    r = new{T}()
    m = nbits(n)
    r.level = Array{Int}(undef, m)
    r.val   = Array{T}(undef, m)
    r.pow   = Array{T}(undef, m+1) #as a temp. var.
    r.last = 0
    r.pow[1] = x
    for i=1:m
      r.val[i] = parent(x)()
    end
    for i=2:m
      r.pow[i] = r.pow[i-1]^2
    end
    r.pow[end] = parent(x)()
    return r
  end
end

function Base.push!(A::EvalEnv{T}, b::T) where T
  if A.last == 0
    A.level[1] = 1
    A.val[1] = copy_into!(A.val[1], b)
    A.last = 1
    return
  end
  if A.level[A.last] > 1
    A.last += 1
    A.level[A.last] = 1
    A.val[A.last] = copy_into!(A.val[A.last], b)
    return
  end

  lev = 1
  #= b enters at level 1:
     x
     b 1
  its multiplied by level 1:
     x
     c 1
  to get
     x^2
     bx+c  1
  this is multiplied by level 2
     x^2
     d 1
  to get
     x^4
     bx^3+cx^2+d 1
  ...
  =#
  b = copy_into!(A.pow[end], b)
  while A.last > 0 && A.level[A.last] <= lev
#    b = b*A.pow[lev] + A.val[A.last]
    b = mul!(b, b, A.pow[lev])
    b = add!(b, b, A.val[A.last])
    lev += 1
    A.last -= 1
  end
  A.last += 1
  A.level[A.last] = lev
  A.val[A.last] = copy_into!(A.val[A.last], b)
  return
end

function finish(A::EvalEnv)
  lst = A.last
  b = A.val[lst]
  while lst >1
    lst -= 1
    lev = A.level[lst]
    b = mul!(b, b, A.pow[lev])
    b = add!(b, b, A.val[lst])
  end
  return b
end

function _my_eval(a, x::T) where T
  b = EvalEnv{T}(x, length(a))
  for x in a
    push!(b, x)
  end
  return finish(b)
end
my_eval(a::AbstractVector{T}, x::T) where {T} = _my_eval(a, x)
function my_eval(f::PolyElem{T}, x::T) where T
  return _my_eval(coefficients(f), x)
end

#coprime base Bach/ Schallit/ ???

function pair_bach(a::E, b::E) where E
  if is_unit(a)
    if is_unit(b)
      return Vector{E}()
    else
      return [b]
    end
  end
  if is_unit(b)
    return [a]
  end

  n = [a, b]
  i = 1
  while i < length(n)
    g = gcd(n[i], n[i+1])
    if isone(g)
      i += 1
    else
      n[i] = divexact!(n[i], g)
      n[i+1] = divexact!(n[i+1], g)
      insert!(n, i+1, g)
      if is_unit(n[i+2])
        deleteat!(n, i+2)
      end
      if is_unit(n[i])
        deleteat!(n, i)
      end
    end
  end

  return n
end

function augment_bach(S::Vector{E}, m::E) where E
  T = Vector{E}()
  i = 1
  while i <= length(S) && !is_unit(m)
    if !is_unit(S[i])
      Ts = pair_bach(m, S[i])
      T = vcat(T, view(Ts, 2:length(Ts)))
      m = Ts[1]
    end
    i += 1
  end
  if i <= length(S)
    T = vcat(T, view(S, i:length(S)))
  end
  if !is_unit(m)
    push!(T, m)
  end
  return T
end

function coprime_base_bach(a::Vector{E}) where E #T need to support GCDs
  if length(a) < 2
    return a
  end

  T = pair_bach(a[1], a[2])
  j = 3
  while j <= length(a)
    T = augment_bach(T, a[j])
    j += 1
  end
  return T
end

# Bernstein: coprime bases
#
#Note: Bernstein needs bigints, either Integer or fmpz
#      well, polys are also OK, small integers are not.

# could/ should be optimsed using divexact! and gcd_into!
# probably should also be combined with ppio somewhere

function ppgle(a::E, b::E) where E
  n = gcd(a,b)
  r = divexact(a, n)
  m = n
  g = gcd(r, n)
  while !isone(g)
    r *= g
    n = divexact(n, g)
    g = gcd(r, n)
  end
  return m, r, n
end

function pair_bernstein(a::E, b::E) where E
  T = Vector{E}()
  if is_unit(b)
    if is_unit(a)
      return T
    else
      return push!(T, a)
    end
  end
  if is_unit(a)
    return push!(T, b)
  end

  a,r = ppio(a,b)
  if !is_unit(r)
    push!(T, r)
  end
  g,h,c = ppgle(a, b)
  c0 = c
  r = c0
  k = 1
  while !isone(h)
    g,h,c = ppgle(h, g^2)
    d = gcd(c, b)
    if isone(d)
      push!(T, c)
      continue
    end
    r *= d
    n = d^(2^(k-1))
    T = vcat(T, pair_bernstein(divexact(c, n), d))
    k += 1
  end
  T = vcat(T, pair_bernstein(divexact(b, r), c0))
  return T
end

function split_bernstein(a::T, P::Hecke.node{T}) where T
  b = Hecke.ppio(a, P.content)[1]
  if !isdefined(P, :left)
    if !isdefined(P, :right)
      return [(P.content, b)]
    else
      return split_bernstein(b, P.right)
    end
  else
    if !isdefined(P, :right)
      return split_bernstein(b, P.left)
    else
      return vcat(split_bernstein(b, P.left), split_bernstein(b, P.right))
    end
  end
end

function split_bernstein(a::T, P::Vector{T}) where T
  if length(P) == 0
    return P
  end
  F = FactorBase(P, check = false)
  b = Hecke.ppio(a, F.prod)[1]
  if length(P)==1
    return [(P[1], b)]
  end
  return vcat(split_bernstein(b, F.ptree.left), split_bernstein(b, F.ptree.right))
end

function augment_bernstein(P::Vector{E}, b::E) where E
  T = Vector{E}()
  if length(P) == 0
    if is_unit(b)
      return T
    else
      return push!(T, b)
    end
  end
  F = FactorBase(P, check = false)
  a,r = Hecke.ppio(b, F.prod)
  if ! is_unit(r)
    push!(T, r)
  end
  S = split_bernstein(a, F.ptree)
  for X in S
    T = vcat(T, pair_bach(X[1], X[2]))
  end
  return T
end

function merge_bernstein(P::Vector{E}, Q::Vector{E}) where E
  m = length(Q)
  b = nbits(m)
  S = P
  i = 0
  while i<=b
    R = prod(view(Q, findall(x -> x & (2^i) ==0, 1:length(Q))))
    T = augment_bernstein(S, R)
    R = prod(view(Q, findall(x -> x & (2^i) !=0, 1:length(Q))))
    S = augment_bernstein(T, R)
    i += 1
  end
  return S
end

function coprime_base_bernstein(S::Vector{E}) where E
  if length(S)<2
    return S
  end
  P1 = coprime_base_bernstein([S[i] for i=1:div(length(S), 2)])
  P2 = coprime_base_bernstein([S[i] for i=(div(length(S), 2)+1):length(S)])
  return merge_bernstein(P1, P2)
end

function augment_steel(S::Vector{E}, a::E, start::Int = 1) where E
  i = start
  if is_unit(a)
    return S
  end

  g = a
  new = true

  while i<=length(S) && !isone(a)
    if new
      g = gcd(S[i], a)
      new = false
    else
      g = gcd_into!(g, S[i], a)
    end
    if is_unit(g)
      i += 1
      continue
    end
    si = divexact(S[i], g)
    a = divexact(a, g)
    if is_unit(si) # g = S[i] and S[i] | a
      continue
    end
    S[i] = si
    if is_unit(a) # g = a and a | S[i]
      a = copy(g)
      continue
    end
    augment_steel(S, copy(g), i)
    continue
  end
  if !is_unit(a)
    push!(S, a)
  end

  return S;
end

function coprime_base_steel(S::Vector{E}) where E
  @assert !isempty(S)
  T = Array{E}(undef, 1)
  T[1] = S[1]
  for i=2:length(S)
    augment_steel(T, S[i])
  end
  return T
end

##implemented
# Bernstein: asymptotically fastest, linear in the total input size
#   pointless for small ints as it requires intermediate numbers to be
#   huge
# Bach/Shallit/???: not too bad, source Epure's Masters thesis
#   can operate on Int types as no intermediate is larger than the input
# Steel: a misnomer: similar to Magma, basically implements a memory
#   optimised version of Bach
#   faster than Magma on
# > I := [Random(1, 10000) * Random(1, 10000) : x in [1..10000]];
# > time c := CoprimeBasis(I);
# julia> I = [fmpz(rand(1:10000))*rand(1:10000) for i in 1:10000];
#
# experimentally, unless the input is enormous, Steel wins
# on smallish input Bach is better than Bernstein, on larger this
# changes
#
# needs
# isone, gcd_into!, divexact!, copy
# (some more for Bernstein: FactorBase, gcd, divexact)

@doc Markdown.doc"""
    coprime_base{E}(S::Vector{E}) -> Vector{E}

Returns a coprime base for $S$, i.e. the resulting array contains pairwise coprime objects that multiplicatively generate the same set as the input array.
"""
coprime_base(x) = coprime_base_steel(x)

@doc Markdown.doc"""
    coprime_base_insert{E}(S::Vector{E}, a::E) -> Vector{E}

Given a coprime array $S$, insert a new element, i.e. find a coprime base for `push(S, a)`.
"""
coprime_base_insert(S, a) = augment_steel(S, a)


text:
using LinearAlgebra
using YAML
using Plots
using Random

parameters = YAML.load(open("Julia/Modules/cfg.yml"))

mutable struct Weights

    e::Array{Array{Float64,2}}
    ie::Array{Any}
    ei::Array{Any}
    de::Array{Array{Float64,2}}

    function Weights(arch::Array{Int64}, param::Dict = parameters)
        # exc
        #Random.seed!(0)
        e = Array{Float64,2}[]
        de = Array{Float64,2}[]
        for l in 1:length(arch)-1
            push!(e, clamp.(randn(arch[l],arch[l+1]) .* param["std"] .+ param["m"], 0, param["wmax"]))
            push!(de, zeros(arch[l],arch[l+1]))
        end
        # inh
        ie = Array{Any}[]
        ei = Array{Any}[]
        for l in 1:length(arch)
            if l>1 && l!=length(arch)
                push!(ie, param["wie"] .* (ones(arch[l],arch[l])-I(arch[l])))
                push!(ei, param["wei"] .* I(arch[l]))
            else
                push!(ie, [[] for _ in 1:arch[l]])
                push!(ei, [[] for _ in 1:arch[l]])
            end
        end
        new(e,ie,ei,de)
    end
end

mutable struct Ind

    exc::Array{Int64}
    inh::Array{Int64}
    el::Array{Array{Int64}}
    il::Array{Array{Int64}}

    function Ind(arch)
        n_exc = sum(arch)
        n_inh = sum(arch[2:end-1])
        n_neurons = n_exc + n_inh
        exc = collect(1:n_exc)
        inh = collect(n_exc+1:n_neurons)

        architr = [0; arch]
        el = Array{Int64}[]
        il = Array{Int64}[]

        for l in 2:length(architr)
            push!(el, collect(sum(architr[1:l-1])+1:sum(architr[1:l])))
            if l > 2 && l!=length(architr)
                push!(il, collect(sum(architr[2:l-1])+1:sum(architr[2:l])) .+ n_exc .- arch[1])
            else
                push!(il, [])
            end
        end
        new(exc,inh,el,il)
    end
end

mutable struct Network

    n_exc::Int64
    n_inh::Int64
    n_neurons::Int64

    i::Ind

    v::Array{Float64}
    u::Array{Float64}
    d::Array{Float64}
    a::Array{Float64}

    s::Weights
    trace_e::Array{Float64}
    ho::Array{Float64}

    I::Array{Float64}

    #RL
    da::Float64
    tderror::Float64
    arch::Array{Int64}
    arch_inh::Array{Int64}
    param::Dict

    function Network(arch::Array{Int64}, param::Dict = parameters)
        n_exc = sum(arch)
        n_inh = sum(arch[2:end-1])
        n_neurons = n_exc + n_inh

        i = Ind(arch)

        v = param["c"]*ones(n_neurons)
        u = 0.2 * v
        d = [(i<=n_exc) ? param["de"] : param["di"] for i in 1:n_neurons]
        a = [(i<=n_exc) ? param["ae"] : param["ai"] for i in 1:n_neurons]

        s = Weights(arch,param)
        trace_e = zeros(n_exc)
        ho = zeros(n_neurons)

        I = zeros(n_neurons)
        da = 0.0
        tderror = 0.0
        arch = arch
        arch_inh = Int64[]
        for l in 1:length(arch)
            if l == 1 || l == length(arch)
                push!(arch_inh, 0)
            else
                push!(arch_inh, div(arch[l],4))
            end
        end

        new(n_exc,n_inh,n_neurons,i,v,u,d,a,s,trace_e,ho,I,da,tderror,arch, arch_inh,param)
    end
end

# %% functions

# direct spike of the neuron
function input!(net::Network, input_spikes::Array{Int64})
    net.v[input_spikes] .= 40.0
end

######## LIF #########

# superthreshold current processed by the neuron
function input_LIF!(net::Network, input_spikes::Array{Tuple{Int64,Float64}})

    net.I = zeros(net.n_neurons) #.+ [(rand()-0.5)*0.2 for k in 1:net.n_neurons]
    max_current = net.param["imin"] + 1.0
    min_current = net.param["imin"]
    indices = [input_spikes[i][1] for i in eachindex(input_spikes)]
    intensity = [input_spikes[i][2] for i in eachindex(input_spikes)]
    net.I[indices] .+= (max_current-min_current) .* intensity .+ min_current
    l1 = collect(1:net.arch[1])
    net.v[l1] .= @. net.v[l1] + (1/net.param["taulif"])*(-net.v[l1]+net.I[l1]+net.param["c"])
end

function spike_LIF!(net::Network)
    # spikes
    net.I = zeros(net.n_neurons)# .+ [(rand()-0.5)*0.2 for k in 1:net.n_neurons]
    spiked = findall(x->x>net.param["vthresh"],net.v .- net.ho)
    spiked_e = filter(x->(x in net.i.exc), spiked)
    spiked_ho = filter(x-> x > net.arch[1], spiked_e)
    net.ho[spiked_e] .+= net.param["theta"] #homeostasis on every layer
    #net.ho[spiked_ho] .+= net.param["theta"] #homeostasis on all but first layer
    net.v[spiked] .= net.param["c"]
    # increment current
    for neuron in spiked
        for l in 1:length(net.arch)-1
            if neuron in net.i.el[l]
                neur = neuron - (net.i.el[l][1] - 1)
                net.I[net.i.el[l+1]] .+= net.s.e[l][neur,:]
                net.I[net.i.il[l]] .+= net.s.ei[l][neur,:]
            end
            if neuron in net.i.il[l]
                neur = neuron - (net.i.il[l][1] - 1)
                net.I[net.i.el[l]] .+= net.s.ie[l][neur,:]
            end
        end
    end

    lr = collect(net.arch[1]+1:net.n_exc) # indexes of all excitatory neuron outside of first layer ( layer rest )
    net.v[lr] .+= net.I[lr]
    net.v[lr] .= @. net.v[lr] + (1/net.param["taulif"])*(-net.v[lr]+net.param["c"])
    net.v[lr] .= max.(net.v[lr], 0.)

    #time step (LIF model)
    net.trace_e[spiked_e] .= 1.0 # not += allows comparison of only last spikes emitted by a neuron

    # LTP & LTD (STDP)
    for neuron in spiked_e
        for l in 1:length(net.arch)
            if neuron in net.i.el[l]
                neur = neuron - (net.i.el[l][1] - 1) # get column number
                if l > 1
                    net.s.de[l-1][:,neur] .+= net.param["apost"] .* net.trace_e[net.i.el[l-1]]
                end
                if l < length(net.arch)
                    net.s.de[l][neur,:] .-= net.param["apre"] .* net.trace_e[net.i.el[l+1]]
                end
            end
        end
    end

    # time decay
    net.ho .*= net.param["ttheta"]
    net.trace_e .*= net.param["ttrace"] # tau_+ = 25 ms #  a mettre en parametre
    net.da *= net.param["tde"] # should be remembering as much as eligibility (de)
    for l in 1:length(net.arch)-1
        net.s.de[l] .*= net.param["tde"]
    end
    spiked
end

function learn!(net::Network) # ( learn with DA-STDP )
    #increase the synaptic weigths
    for l in 1:length(net.arch)-1
        #net.da = clamp(net.da, -param["damax"], param["damax"])
        net.s.e[l] .+= net.da .* net.s.de[l]
        net.s.e[l] .= clamp.(net.s.e[l], 0, net.param["wmax"])
    end
end

function learn_critic!(net::Network)
    #increase the synaptic weigths
    for l in 1:length(net.arch)-1
        #net.da = clamp(net.da, -param["damax"], param["damax"])
        net.s.e[l] .+= net.param["eta"] .* net.da .* net.s.de[l]
        net.s.e[l] .= clamp.(net.s.e[l], 0, net.param["wmax"])
    end
    for l in 1:length(net.arch_inh)-1
        n_inh = net.arch_inh[l]
        for i in 1:n_inh
            net.s.e[l][net.arch[l]+1-i, :] .= net.param["winh"]
        end
    end
end

# %% main functions

function step!(net::Network, train::Bool=false)
    spiked = spike!(net)
    if train
        learn!(net)
    end
    spiked
end

function step!(net::Network, input_spikes::Array{Int64}, train::Bool=false)
    input!(net,input_spikes)
    spiked = spike!(net)
    if train
        learn!(net)
    end
    spiked
end

function step_LIF!(net::Network, input_spikes::Array{Tuple{Int64,Float64}}, train::Bool=false)
    input_LIF!(net,input_spikes)
    spiked = spike_LIF!(net)
    if train
        learn!(net)
    end
    spiked
end

function step_LIF_critic!(net::Network, input_spikes::Array{Tuple{Int64,Float64}}, train::Bool=false)
    input_LIF!(net,input_spikes)
    spiked = spike_LIF!(net)
    if train
        learn_critic!(net)
    end
    spiked
end

function step!(net::Network, input_spikes::Array{Tuple{Int64,Float64}}, train::Bool=false)
    input!(net,input_spikes)
    spiked = spike!(net)
    if train
        learn!(net)
    end
    spiked
end

# %% visualization function

function weight_viz(net::Network, method::String = "histogram")
    weights = net.s.e
    if method == "histogram"
        h = []
        for i in 1:length(net.arch)-1
            vec = filter(x-> x>=0 ,weights[i])
            h_tmp = histogram(vec,
                          legend = false,
                          xlabel = "Synaptic strength layer $(i)",
                          ylabel = "Synapse number",
                          bins = 10)
            push!(h, h_tmp)
        end
        display(plot(h..., layout = length(net.arch)-1))

    elseif method == "heatmap"
        h = []
        for i in 1:length(net.arch)-1
            h_tmp = heatmap(weights[i],
                            xlabel = "Connection matrix layer $(i) - $(i+1)")
            push!(h, h_tmp)
        end
        display(plot(h..., layout = length(net.arch)-1))

    end
end

########## Izhikevic #############

# superthreshold current processed by the neuron
function input!(net::Network, input_spikes::Array{Tuple{Int64,Float64}})
    #net.I = 13 .* (rand(net.n_neurons) .- 0.5)
    net.I = zeros(net.n_neurons)
    max_current = net.param["imax"]
    min_current = net.param["imin"]
    indices = [input_spikes[i][1] for i in eachindex(input_spikes)]
    intensity = [input_spikes[i][2] for i in eachindex(input_spikes)]
    net.I[indices] .+= (max_current-min_current) .* intensity .+ min_current
    net.v[indices] .= @. net.v[indices]+0.5*((0.04*net.v[indices]+5)*net.v[indices]+140-net.u[indices]+net.I[indices])
    net.v[indices] .= @. net.v[indices]+0.5*((0.04*net.v[indices]+5)*net.v[indices]+140-net.u[indices]+net.I[indices])
end

function spike!(net::Network)
    # spikes
    net.I = zeros(net.n_neurons)
    spiked = findall(x->x>net.param["vthresh"],net.v .- net.ho)
    spiked_e = filter(x->(x in net.i.exc), spiked)
    spiked_ho = filter(x-> x > net.arch[1], spiked_e)
    net.ho[spiked_ho] .+= net.param["theta"]
    net.v[spiked] .= net.param["c"]
    net.u[spiked] .+= net.d[spiked]
    # increment current
    for neuron in spiked
        for l in 1:length(net.arch)-1
            if neuron in net.i.el[l]
                neur = neuron - (net.i.el[l][1] - 1)
                net.I[net.i.el[l+1]] .+= net.s.e[l][neur,:]
                net.I[net.i.il[l]] .+= net.s.ei[l][neur,:]
            end
            if neuron in net.i.il[l]
                neur = neuron - (net.i.il[l][1] - 1)
                net.I[net.i.el[l]] .+= net.s.ie[l][neur,:]
            end
        end
    end

    #time step (Izhikevic model)
    net.v .= @. net.v+0.5*((0.04*net.v+5)*net.v+140-net.u+net.I)
    net.v .= @. net.v+0.5*((0.04*net.v+5)*net.v+140-net.u+net.I)
    net.u .= @. net.u+net.a*(0.2*net.v-net.u)
    net.trace_e[spiked_e] .= 0.1

    # LTP & LTD (STDP)
    for neuron in spiked_e
        for l in 1:length(net.arch)
            if neuron in net.i.el[l]
                neur = neuron - (net.i.el[l][1] - 1)
                if l > 1
                    net.s.de[l-1][:,neur] .+= net.param["apost"] .* net.trace_e[net.i.el[l-1]]
                end
                if l < length(net.arch)
                    net.s.de[l][neur,:] .-= net.param["apre"] .* net.trace_e[net.i.el[l+1]]
                end
            end
        end
    end

    # time decay
    net.ho .*= net.param["ttheta"]
    net.trace_e .*= net.param["ttrace"]
    net.da *= net.param["tda"]
    spiked
end

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
# Undefined column fail test
@test_throws MethodError DataView()

# Set up data for data view call overloading
test00id = rand(UInt64)
test00index = rand(Int64)
test00name = rand([collect('A':'Z'); collect('a':'z')])
test00symbol = Symbol(test00name)

# Build the singleton types
_DataView(test00id, test00symbol, test00index)
@test DataView(Val{test00id}, Val{test00symbol}) == test00index
@test DataView(Val{test00id}, Val{test00index}) == test00symbol

# Set up test data for courtesy functions
test01id = rand(UInt64)
test01index = rand(Int64)
test01name = rand([collect('A':'Z'); collect('a':'z')])
test01symbol = Symbol(test01name)

# Build the singleton types
setcolumnname(test01id, test01symbol, test01index)
@test DataView(Val{test01id}, Val{test01symbol}) == test01index
@test DataView(Val{test01id}, Val{test01index}) == test01symbol

# Set up test data for courtesy functions
test02id = rand(UInt64)
test02index = rand(Int64)
test02name = rand([collect('A':'Z'); collect('a':'z')])
test02symbol = Symbol(test02name)

# Build the singleton types
setcolumnname(test02id, test02name, test02index)
@test DataView(Val{test02id}, Val{test02symbol}) == test02index
@test DataView(Val{test02id}, Val{test02index}) == test02symbol

# Make sure unique check works
@test_throws ErrorException setcolumnnames(rand(UInt64), fill(Symbol(rand([collect('A':'Z'); collect('a':'z')])), rand(1:10)))

# Set up test data for courtesy functions
test03id = rand(UInt64)
test03indexes = collect(1:52)
test03names = [collect('A':'Z'); collect('a':'z')][randperm(52)]
test03symbols = [Symbol(t) for t in test03names]

# Build the singleton types
setcolumnnames(test03id, test03symbols)
@test all([DataView(Val{test03id}, Val{test03symbols[i]}) == test03indexes[i] for i = 1:52])
@test all([DataView(Val{test03id}, Val{test03indexes[i]}) == test03symbols[i] for i = 1:52])

# Set up test data for courtesy functions
test04id = rand(UInt64)
test04indexes = collect(1:52)
test04names = [collect('A':'Z'); collect('a':'z')][randperm(52)]
test04symbols = [Symbol(t) for t in test04names]

# Build the singleton types
setcolumnnames(test04id, test04names)
@test all([(DataView(Val{test04id}, Val{test04symbols[i]}) == test04indexes[i])::Bool for i = 1:52])
@test all([(DataView(Val{test04id}, Val{test04indexes[i]}) == test04symbols[i])::Bool for i = 1:52])

# Set up test data for courtesy functions
test05id = rand(UInt64)
test05index = rand(Int64)
test05name = rand([collect('A':'Z'); collect('a':'z')])
test05symbol = Symbol(test05name)

# Build the singleton types
setcolumnname(test05id, test05name, test05index)
@test getcolumnindex(test05id, test05symbol) == test05index
@test getcolumnindex(test05id, test05name) == test05index
@test getcolumnname(test05id, test05index) == test05symbol


text:
"""
abstract type Model

Models hold the parameters of interest and the associated likelihood/cost
functions used in optimization. Individual models should be derived from
this base object.
"""
abstract type Model end

"""
Results

Created upon optimizing a model, contains the model and all of
its parameters as well as the optimized results found by
minimizing the cost function.
"""
mutable struct Results
  opt
  model::Model
  min_objective::Float64
  min_parameters::Array{Float64}
  iterations::Int64
  lower_bounds::Array{Float64}
  upper_bounds::Array{Float64}
end

"""
    optimize_model!(m::Model, nll::NLogLikelihood;
                    lower_bounds=nothing, upper_bounds=nothing)

Currently NLopt is used for minimizing the objective function.
optimize_model! produces the objective function from the given
NLogLikelihood which is then passed to NLopt to be minimized.
optimize_model! returns a set of Results.

TODO: Move lower_bounds and upper_bounds into options dict
"""
function optimize_model!(m::Model, nll::NLogLikelihood;
                         lower_bounds=nothing, upper_bounds=nothing,
                         options=Dict() )
  # Selection of optimizers.
  optimizer_selector = Dict(
                            "sbplx"  => :LN_SBPLX,
                            "cobyla" => :LN_COBYLA,
                            "bobyqa" => :LN_BOBYQA,
                            "praxis" => :LN_PRAXIS,
                            "nelder" => :LN_NELDERMEAD
                           )
  opt_string = get(options, "optimizer", "sbplx")
  @debug "Using optimizer" opt_string
  opt_use = optimizer_selector[opt_string]
  function objective(x::Vector, grad::Vector)
    if length(grad)>0
      grad = 2x
    end
    nll.objective(x)
  end
  #opt_use = get(options, "optimizer", :LN_SBPLX)
  opt = Opt(opt_use, nll.numparams)
  # A note here on tolerence. With nlopt one can specify four
  # types of tolerence for a stopping condition. Either on the
  # value of the objective function, or on the values of the parameters.
  # For both of these it can be relative to the scale of each, or
  # absolute. For simplicity I have forced in tolerance on the objective
  # function at an absolute scale. This is because I am assuming that
  # the function is a log-likelihood, and thus has statistical information
  # in its derivative, but may be shifting by a constant NOT a scale factor.
  p0 = get(options, "p0", [v.init for v in nll.parameters])
  trials = get(options, "trials", 1)

  opt.ftol_abs     = get(options, "ftol_abs", 1e-6)
  opt.ftol_rel     = get(options, "ftol_rel", 0)
  opt.xtol_rel     = get(options, "xtol_rel", 0)
  opt.xtol_abs     = get(options, "xtol_abs", 0)
  opt.maxeval      = get(options, "maxeval", 0)
  opt.maxtime      = get(options, "maxtime", 0)
  initial_step     = get(options, "initial_step", p0./10.0)
  initial_step[(initial_step .== 0)] .= 0.1
  opt.initial_step = initial_step
  opt.stopval      = get(options, "stopval", -Inf)
  opt.lower_bounds = nll.lower_bounds
  opt.upper_bounds = nll.upper_bounds

  opt.min_objective = objective

  # Global Opt

  #= Skipping the global optimization
  gopt = Opt(:GN_CRS2_LM, nll.numparams)
  global_upper = p0 * 100
  #global_upper[(global_upper .< 100)] .= 100
  gopt.lower_bounds = -global_upper
  gopt.upper_bounds = global_upper
  #gopt.ftol_rel = opt.ftol_rel
  gopt.xtol_abs = opt.xtol_abs .* 1_000.0
  gopt.min_objective = objective
  (minf, minx, ret) = optimize!(gopt, p0)
  p0 = minx
  =#

  # FIXME: This always fails
  #try
  #  nll.objective(p0)
  #catch e
  #  @error "Problems contructing objective function"
  #  println(e)
  #end
  # We will loop over trials, if trials > 1 then we will random sample
  best_minf = 1e9
  if trials > 1
    # We would like to sample a few fits then start near the best
    minf_set = Array{Float64}(undef, 0)
    p0_set = Array{Array{Float64}}(undef, 0)
    for tr in collect(1:trials)
      p0 = Array{Float64}(undef, 0)
      for v in nll.parameters
        if v.samplewidth != Inf
          push!(p0, randn()*v.samplewidth + v.init)
        else
          push!(p0, v.init)
        end
      end
      (minf, minx, ret) = optimize!(opt, p0)
      push!(p0_set, p0)
      push!(minf_set, minf)
    end
    p0 = p0_set[ (minf_set .== minimum(minf_set)) ][1]
    best_minf = minimum(minf_set)
  end

  # Test before optimizing
  # (minf, minx, ret) = optimize!(opt, p0)
  ## If we fit twice we can go coarse the first time and then re-evaluate
  # initial_step = minx ./ 10.0
  initial_step[(initial_step .== 0)] .= 0.1
  opt.initial_step = initial_step
  (minf, minx, ret) = optimize!(opt, p0)

  if minf > best_minf
    @show best_min minf
  end

  numevals = opt.numevals
  for (i,p) in enumerate(m.params)
    p.fit = minx[i]
  end
  Results( opt, m, copy(minf), copy(minx), numevals,
           copy(opt.lower_bounds), copy(opt.upper_bounds) )
end

function getparam(m::Model, name::String)
  for (i, p) in enumerate(m.params)
    if p.name == name
      return p
    end
  end
  return nothing
end

function add_likelihood!( m::Model, nll::NLogLikelihood )
  m.nll = nll
  m.params = nll.parameters
end


text:
function useproc(proc)
    global PROC
    if (proc=="CPU") | (proc=="GPU") | (proc=="GPU32")
        PROC=proc
    end
    include(ADDIR*"/src/gpumacros.jl")

    f=open(ADDIR*"/src/initfile.jl","w")
    if (proc=="GPU")
        write(f,"global PROC; PROC=\"GPU\"\n")
    end
    if (proc=="GPU32")
        write(f,"global PROC; PROC=\"GPU32\"\n")
    end
    if (proc=="CPU")
        write(f,"global PROC; PROC=\"CPU\"\n")
    end
    close(f)
end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
@doc raw"""
Lie-Trotter Splitting A

For a vector field $\dot{x} = f_1 (t,x) +  f_2 (t,x)$, the splitting reads
```math
\Phi_{\Delta t} = \varphi^{2}_{\Delta t} \circ \varphi^{1}_{\Delta t}
```

Reference:

    H. F. Trotter.
    On the product of semi-groups of operators.
    Proceedings of the American Mathematical Society, Volume 10, Pages 545-551, 1959.
    doi: 10.1090/S0002-9939-1959-0108732-6.

"""
function TableauLieA(::Type{T}=Float64) where {T}
    a = Array{T}([ 1 ])
    b = Array{T}([ 0 ])
    TableauSplittingNS(:LieTrotterSplittingA, 1, a, b)
end

@doc raw"""
Lie-Trotter Splitting B

For a vector field $\dot{x} = f_1 (t,x) +  f_2 (t,x)$, the splitting reads
```math
\Phi_{\Delta t} = \varphi^{1}_{\Delta t} \circ \varphi^{2}_{\Delta t}
```

Reference:

    H. F. Trotter.
    On the product of semi-groups of operators.
    Proceedings of the American Mathematical Society, Volume 10, Pages 545-551, 1959.
    doi: 10.1090/S0002-9939-1959-0108732-6.

"""
function TableauLieB(::Type{T}=Float64) where {T}
    a = Array{T}([ 0 ])
    b = Array{T}([ 1 ])
    TableauSplittingNS(:LieTrotterSplittingB, 1, a, b)
end

@doc raw"""
Strang Splitting

For a vector field $\dot{x} = f_1 (t,x) + f_2 (t,x)$, the splitting reads
```math
\Phi_{\Delta t} = \varphi^{1}_{\Delta t / 2} \circ \varphi^{2}_{\Delta t / 2} \circ \varphi^{2}_{\Delta t / 2} \circ \varphi^{1}_{\Delta t / 2}
```

For vector fields with two components, this is not the most efficient implementation.
For such cases [`TableauStrangA`](@ref) or [`TableauStrangB`](@ref) should be used instead.

References:

    Gilbert Strang.
    On the construction and comparison of difference schemes.
    SIAM Journal on Numerical Analysis, Volume 5, Pages 506-517, 1968.
    doi: 10.1137/0705041.

    Gurij Ivanovich Marchuk.
    Some applications of splitting-up methods to the solution of mathematical physics problems.
    Aplikace Matematiky, Volume 13, Pages 103-132, 1968.
    doi: 10.21136/AM.1968.103142.

"""
function TableauStrang(::Type{T}=Float64) where {T}
    a = Array{T}([ 1//2 ])
    b = Array{T}([ 1//2 ])
    TableauSplittingNS(:StrangSplitting, 2, a, b)
end

"Alias for [`TableauStrang`](@ref)"
TableauMarchuk = TableauStrang

@doc raw"""
Strang Splitting A for a vector field $\dot{x} = f_1 (t,x) + f_2 (t,x)$.

The splitting reads
```math
\Phi_{\Delta t} = \varphi^{1}_{\Delta t / 2} \circ \varphi^{2}_{\Delta t} \circ \varphi^{1}_{\Delta t / 2}
```

References:

    Gilbert Strang.
    On the construction and comparison of difference schemes.
    SIAM Journal on Numerical Analysis, Volume 5, Pages 506-517, 1968.
    doi: 10.1137/0705041.

    Gurij Ivanovich Marchuk.
    Some applications of splitting-up methods to the solution of mathematical physics problems.
    Aplikace Matematiky, Volume 13, Pages 103-132, 1968.
    doi: 10.21136/AM.1968.103142.

"""
function TableauStrangA(::Type{T}=Float64) where {T}
    a = Array{T}([ 1//2, 1//2 ])
    b = Array{T}([ 1//1, 0//1 ])
    TableauSplitting(:StrangSplittingA, 2, a, b)
end

@doc raw"""
Strang Splitting B for a vector field $\dot{x} = f_1 (t,x) + f_2 (t,x)$

The splitting reads
```math
\Phi_{\Delta t} = \varphi^{2}_{\Delta t / 2} \circ \varphi^{1}_{\Delta t} \circ \varphi^{2}_{\Delta t / 2}
```

References:

    Gilbert Strang.
    On the construction and comparison of difference schemes.
    SIAM Journal on Numerical Analysis, Volume 5, Pages 506-517, 1968.
    doi: 10.1137/0705041.

    Gurij Ivanovich Marchuk.
    Some applications of splitting-up methods to the solution of mathematical physics problems.
    Aplikace Matematiky, Volume 13, Pages 103-132, 1968.
    doi: 10.21136/AM.1968.103142.

"""
function TableauStrangB(::Type{T}=Float64) where {T}
    a = Array{T}([ 0//1, 1//1 ])
    b = Array{T}([ 1//2, 1//2 ])
    TableauSplitting(:StrangSplittingB, 2, a, b)
end

@doc raw"""
McLachlan's 2nd order symmetric, minimum error composition method

The composition reads
```math
\Phi_{\Delta t} = \varphi_{\alpha \Delta t} \circ \varphi^{*}_{(1/2 - \alpha) \Delta t} \circ \varphi_{(1/2 - \alpha) \Delta t} \circ \varphi^{*}_{\alpha \Delta t} ,
```
where the parameter $\alpha$ can be optimized, e.g., to minimize the solution error.
McLachlan arrives at $\alpha  = 0.1932$ as a generally useful value.

Reference:

    Robert I. McLachlan.
    On the Numerical Integration of Ordinary Differential Equations by Symmetric Composition Methods
    SIAM Journal on Scientific Computing, Volume 16, Pages 151-168, 1995.
    doi: 10.1137/0916010.

"""
function TableauMcLachlan2(::Type{T}=Float64; α=0.1932) where {T}
    a = Array{T}([ α, 0.5 - α ])
    b = Array{T}([ 0.5 - α, α ])
    TableauSplittingNS(:McLachlanSplitting, 2, a, b)
end

@doc raw"""
McLachlan's 4th order symmetric, minimum error composition method

The composition reads
```math
\Phi_{\Delta t} = \varphi_{\alpha_5 \Delta t} \circ \varphi^{*}_{\beta_5 \Delta t} \circ \dotsc \circ \varphi_{\alpha_2 \Delta t} \circ \varphi^{*}_{\beta_2 \Delta t} \circ \varphi_{\alpha_1 \Delta t} \circ \varphi^{*}_{\beta_1 \Delta t} ,
```
with
```math
\begin{aligned}
\beta_1 &= \alpha_5 = \frac{14 - \sqrt{19}}{108} , &
\alpha_1 &= \beta_5 = \frac{146 + 5 \sqrt{19}}{540} , & & \\
\beta_2 &= \alpha_4 = \frac{- 23 - 20 \sqrt{19}}{270} , &
\alpha_2 &= \beta_4 = \frac{-2 + 10 \sqrt{19}}{135} , &
\beta_3 &= \alpha_3 = \frac{1}{5} .
\end{aligned}
```
The coefficients are optimised to provide an integrator with minimal solution error.

Reference:

    Robert I. McLachlan.
    On the Numerical Integration of Ordinary Differential Equations by Symmetric Composition Methods
    SIAM Journal on Scientific Computing, Volume 16, Pages 151-168, 1995.
    doi: 10.1137/0916010.

"""
function TableauMcLachlan4(::Type{T}=Float64) where {T}
    a = Array{T}(@big [ (146 +  5*√19) / 540,
                        ( -2 + 10*√19) / 135,
                                     1 / 5,
                        (-23 - 20*√19) / 270,
                        ( 14 -    √19) / 108])
    TableauSplittingNS(:McLachlanSplitting, 4, a, a[end:-1:1])
end

@doc raw"""
4th order "Triple Jump" composition method.

The composition reads
```math
\Phi_{\Delta t} = \varphi_{\gamma_3 \Delta t} \circ \varphi_{\gamma_2 \Delta t} \circ \varphi_{\gamma_1 \Delta t}
```
with
```math
\gamma_1 = \gamma_3 = \frac{1}{2 - 2^{1/(p+1)}} , \qquad
\gamma_2 = - \frac{2^{1/(p+1)}}{2 - 2^{1/(p+1)}} .
```

References:

    Michael Creutz and Andreas Gocksch.
    Higher-order hybrid Monte Carlo algorithms.
    Physical Review Letters, Volume 63, Pages 9-12, 1989.
    doi: 10.1103/PhysRevLett.63.9.

    Etienne Forest.
    Canonical integrators as tracking codes (or how to integrate perturbation theory with tracking).
    AIP Conference Proceedings, Volume 184, Pages 1106-1136, 1989.
    doi: 10.1063/1.38062.

    Masuo Suzuki
    Fractal decomposition of exponential operators with applications to many-body theories and Monte Carlo simulations.
    Physics Letters A, Volume 146, Pages 319-323, 1990.
    doi: 10.1016/0375-9601(90)90962-N

    Haruo Yoshida.
    Construction of higher order symplectic integrators.
    Physics Letters A, Volume 150, Pages 262-268, 1990.
    doi: 10.1016/0375-9601(90)90092-3

"""
function TableauTripleJump(::Type{T}=Float64) where {T}
    fac = @big 2^(1/3)
    den = @big 1/(2-fac)
    a = Array{T}([ den, -fac*den ])
    TableauSplittingSS(:TripleJumpSplitting, 4, a)
end

@doc raw"""
Suzuki's 4th order "fractal" composition method

The composition reads
```math
\Phi_{\Delta t} = \varphi_{\gamma_5 \Delta t} \circ \varphi_{\gamma_4 \Delta t} \circ \varphi_{\gamma_3 \Delta t} \circ \varphi_{\gamma_2 \Delta t} \circ \varphi_{\gamma_1 \Delta t}
```
with
```math
\gamma_1 = \gamma_2 = \gamma_4 = \gamma_5 = \frac{1}{4 - 4^{1/(p+1)}} , \qquad
\gamma_3 = - \frac{4^{1/(p+1)}}{4 - 4^{1/(p+1)}} .
```

Reference:

    Masuo Suzuki
    Fractal decomposition of exponential operators with applications to many-body theories and Monte Carlo simulations.
    Physics Letters A, Volume 146, Pages 319-323, 1990.
    doi: 10.1016/0375-9601(90)90962-N

"""
function TableauSuzukiFractal(::Type{T}=Float64) where {T}
    fac = @big 4^(1/3)
    den = @big 1/(4-fac)
    a = Array{T}([ den, den, -fac*den ])
    TableauSplittingSS(:SuzukiFractalSplitting, 4, a)
end


text:
@testset "Matrix algebras" begin

  @testset "Construction" begin

    M1 = matrix(FlintQQ, [ 1 1; 0 1])
    M2 = matrix(FlintQQ, [ 1 0; 1 1])
    A = matrix_algebra(FlintQQ, [ M1, M2 ])
    @test dim(A) == 4

    QG = group_algebra(FlintQQ, small_group(2, 1))
    M1 = matrix(QG, [ 1 1; 0 1])
    M2 = matrix(QG, [ 1 0; 1 1])
    A = matrix_algebra(FlintQQ, QG, [ M1, M2 ])
    @test dim(A) == 4
  end

end


text:
using EzXML
using DataStructures
using LightGraphs
using Vulkan_Headers_jll:vk_xml

xdoc = readxml(vk_xml)
xroot = xdoc.root

include("utils.jl")
include("handles.jl")
include("graph.jl")

base_types_exceptions = Dict(
    "CAMetalLayer" => "void",
    "ANativeWindow" => "void",
    "AHardwareBuffer" => "void",
)

vk_base_types_mapping = Dict(
    ("uint$(size)_t" => "UInt$size" for size ∈ (8, 16, 32, 64))...,
    ("int$(size)_t" => "Int$size" for size ∈ (8, 16, 32, 64))...,
    "float" => "Float32",
    "double" => "Float64",
    "void" => "Cvoid",
)

function translate_base_type_c(base_type)
    base_type ∉ keys(vk_base_types_mapping) && error("Unknown base type $base_type")
    vk_base_types_mapping[base_type]
end

function fetch_base_types(xroot)
    vk_base_types_nodes = findall("//type[@category='basetype']", xroot)
    names = member_attr.(vk_base_types_nodes, "name")
    println.(vk_base_types_nodes)
    res = Dict()
    for (i, name) ∈ enumerate(names)
        res[name] = translate_base_type_c(name ∈ keys(base_types_exceptions) ? base_types_exceptions[name] : member_attr(vk_base_types_nodes[i], "type"))
    end
    res
end

base_types_vk = fetch_base_types(xroot)

function translate_base_type_vk(base_type)
    base_type ∉ keys(base_types_vk) && error("Unknown vulkan type $base_type")
    base_types_vk[base_type]
end

"""
    translate_type(base_type)

Translates a type from available C-based or Vulkan-based definitions.

# Examples
```
julia> translate_type("uint_32_t")
"UInt32"
julia> translate_type("VkBool32")
"UInt32"
```
"""
function translate_type(base_type)
    base_type ∈ keys(base_types_vk) && return translate_base_type_vk(base_type)
    base_type ∈ keys(vk_base_types_mapping) && return translate_base_type_c(base_type)
    error("Unknown type $base_type")
end


text:
using ArgParse, Logging

function isprinted(level, logger=global_logger())
    return Logging.min_enabled_level(logger) <= level
end

function isdebug(logger=global_logger())
    return isprinted(Logging.Debug, logger)
end

include("types.jl")

function parse_cmdargs()
    s = ArgParseSettings(prog="jkMST", version="1.0", add_version=true)

    @add_arg_table! s begin
        "--directory", "-d"
            help = "A directory containing graph files. All subdirectories are considered. Only files with extension .dat are read."
            nargs = '*'
            arg_type = String
        "--file", "-f"
            help = "A graph represented as file. Any file extension is permitted."
            nargs = '*'
            arg_type = String
        "--mode", "-m"
            help = "The solving method. Every instance is solved with every argument. (Allowed values: $(string([string(v) for v in instances(SolvingMode)])))"
            default = ["mtz"]
            nargs = '+'
            arg_type = String
        "--size", "-k"
            help = "Size of the k-MST. If k ∈ [0,1] interpreted as fraction of |V|. If k > 1 interpreted as absolute value. Every instance is solved with every argument."
            default = ["1"]
            nargs = '+'
            arg_type = String
        "--solver", "-s"
            help = "The solver. Every instance is solved with every argument. (Allowed values: $(string([string(v) for v in instances(Solver)]))"
            default = ["cplex"]
            nargs = '+'
            arg_type = String
        "--timeout", "-t"
            help = "Maximum time in seconds spent on a single instance (3600s = 1h)."
            default = Inf64
            arg_type = Float64
        "--generate-timeout"
            help = "Maximum time in seconds spent on generating a model (600s = 10m, 3600s = 1h)."
            default = Inf64
            arg_type = Float64
        "--no-preprocess-instances"
            help = "Do not run Kruskal and Prim before solving file inorder to obtain better starting bounds."
            action = :store_true
        "--preprocess-timeout"
            help = "Maximum amount of time spent on preprocessing instances."
            default = Inf64
            arg_type = Float64
        "--no-print-solution-graphs", "--no-print-solutions"
            help = "Suppress the printing of the solution graphs after a solution is found."
            action = :store_true
        "--intermediate-table-interval"
            help = "Set the interval in seconds for printing the table in between instances."
            default = 60.0
            arg_type = Float64
        "--no-print-intermediate-table"
            help = "Do not print the table in between instances."
            action = :store_true
        "--no-print-models"
            help = "Do not print the models in debug mode."
            action = :store_true
        "--table-style", "--table-format"
            help = "Choose table style. (Allowed values: [\"text\", \"html\", \"latex\"])"
            default = "text"
            arg_type = String
        "--system-info"
            help = "Display what properties the system detects."
            action = :store_true
        "--verbose", "-v"
            help = "Increase verbosity."
            action = :count_invocations
        "--quiet", "-q"
            help = "Decrease verbosity."
            action = :count_invocations
    end

    return parse_args(s)
end

parsed_args = parse_cmdargs()

system_info = parsed_args["system-info"]

function test_system()
    rng = MersenneTwister()

    ops::Int = 10^7
    flops::Int = ops * 4
    operation_time::Float64 = 0.0
        for j in 1:ops
            s = rand(rng, Float64)
            m = rand(rng, Float64)
            operation_time += @elapsed r = m * s
            operation_time += @elapsed r = r - m
            operation_time += @elapsed r = s / m
            operation_time += @elapsed r = 5.0 + m
    end

    flopss = flops/operation_time
    prefix = ""
    if flopss > 1000
        prefix = "Kilo"
        flopss /= 1000
        if flopss > 1000
            prefix = "Mega"
            flopss /= 1000
            if flopss > 1000
                prefix = "Giga"
                flopss /= 1000
                if flopss > 1000
                    prefix = "Terra"
                    flopss /= 1000
                    if flopss > 1000
                        prefix = "Peta"
                        flopss /= 1000
                    end
                end
            end
        end
    end

    @info "CPUs: $(Base.Sys.CPU_THREADS) Threads: $(Base.Threads.nthreads()), estimated FLOPS: $(round(flopss; digits = 2)) $(prefix) FLOP/s, Total memory: $(round(Int,Base.Sys.total_memory() // 2^20)) MiB, Free memory: $(round(Int, Base.Sys.free_memory() // 2^20)) MiB"

end

if system_info
    using Random, Dates
    test_system()
    exit(0)
end

verbosity = parsed_args["verbose"]
quiet = parsed_args["quiet"]
level = loglevels[max(1, min(4, 2 + (quiet - verbosity)))]
logger = Logging.ConsoleLogger(stdout, level)
global_logger(logger)
@debug "LogLevel: $(Logging.min_enabled_level(logger)), isdebug: $(isdebug()), threads: $(Base.Threads.nthreads())."
files = parsed_args["file"]
directories = parsed_args["directory"]
modestrings = parsed_args["mode"]
kstrings = parsed_args["size"]
solverstrings = parsed_args["solver"]
timeout = parsed_args["timeout"]
generate_timeout = parsed_args["generate-timeout"]
preprocessinstances = !parsed_args["no-preprocess-instances"]
preprocesstimeout = parsed_args["preprocess-timeout"]
printsolutiongraphs = !parsed_args["no-print-solution-graphs"]
intermediatetableinterval = parsed_args["intermediate-table-interval"]
printintermediatetable = !parsed_args["no-print-intermediate-table"]
debugmodels = !parsed_args["no-print-models"]
tablestyle = parsed_args["table-style"]

if isempty(files) && isempty(directories)
    exit()
end

@debug "Compiling program."
includetime = @elapsed include("jkMST.jl")
@debug "Compiled program in $(format_seconds_readable(includetime))."

main(files=files,
    directories=directories,
    modestrings=modestrings,
    kstrings=kstrings,
    solverstrings=solverstrings,
    timeout_sec=timeout,
    generate_timeout_sec=generate_timeout,
    preprocessinstances=preprocessinstances,
    preprocesstimeout=preprocesstimeout,
    printsolutiongraphs=printsolutiongraphs,
    printintermediatetable=printintermediatetable,
    debugmodels=debugmodels,
    intermediatetableinterval=intermediatetableinterval,
    tablestyle=tablestyle
    )

max_stars_repo_name: Entze/jkMST

max_issues_repo_name: Entze/jkMST

max_forks_repo_name: Entze/jkMST


text:
"""
    alt_xy(r::Integer, c::Integer)
    alt_xy(rc::Tuple{T, U}) where {T <: Integer, U <: Integer}
Construct a matrix with `r` rows and `c` columns, with elements
alternating between `+1.0` and `-1.0`.
The first element is +1.
"""
function alt_xy(r::Integer, c::Integer)
    a = repeat([1.0, -1.0], outer = div(r + 1, 2))
    isodd(r) && (a = a[1:end - 1])
    b = repeat(vcat(a, -a), outer = div(c + 1, 2))
    isodd(c) && (b = b[1:r * c])
    reshape(b, r, c)
end

alt_xy(rc::Tuple{T, U}) where {T <: Integer, U <: Integer} = alt_xy(rc[1], rc[2])

"""
    alt_x(r::Integer, c::Integer)
    alt_x(rc::Tuple{T, U}) where {T <: Integer, U <: Integer}
Construct a matrix with `r` rows and `c` columns, with columns
alternating between `+1.0` and `-1.0`.
The first column is +1.
"""
function alt_x(r::Integer, c::Integer)
    a = repeat([1.0, -1.0], inner = r)
    b = repeat(a, outer = div(c + 1, 2))
    isodd(c) && (b = b[1:r * c])
    reshape(b, r, c)
end

alt_x(rc::Tuple{T, U}) where {T <: Integer, U <: Integer} = alt_x(rc[1], rc[2])

"""
    alt_y(r::Integer, c::Integer)
    alt_y(rc::Tuple{T, U}) where {T <: Integer, U <: Integer}
Construct a matrix with `r` rows and `c` columns, with rows
alternating between `+1.0` and `-1.0`.
The first row is +1.
"""
function alt_y(r::Integer, c::Integer)
    a = repeat([1.0, -1.0], outer = div(r + 1, 2))
    isodd(r) && (a = a[1:end-1])
    b = repeat(a, outer = c)
    reshape(b, r, c)
end

alt_y(rc::Tuple{T, U}) where {T <: Integer, U <: Integer} = alt_y(rc[1], rc[2])

"""
    fft(eb[, recentre = true])
Returns the two-dimensional DFT of Eigenbrot `eb`. If
`recentre` is `true` the DFT is stored so that the value
at the coordinate origin is at the central location of
the data.
"""
function FFTW.fft(eb::Eigenbrot, recentre::Bool = true)
    w = width(eb)
    h = height(eb)
    scale = 1.0 / sqrt(h * w)
    if recentre
        alt = alt_xy(h, w)
        data = alt .* eb.vals
        eb.fft ? fft!(data) : bfft!(data)
        Eigenbrot(scale * (alt .* data), !eb.fft)
    else
        data = copy(eb.vals)
        eb.fft ? fft!(data) : bfft!(data)
        Eigenbrot(scale * data, !eb.fft)
    end
end

"""
    fftx(eb[, recentre = true])
Returns the one-dimensional DFT of Eigenbrot `eb` in
the `x`-direction. If `recentre` is `true` the DFT is
stored so that the value at the coordinate origin is
at the central location of the data.
"""
function fftx(eb::Eigenbrot, recentre::Bool = true)
    w = width(eb)
    h = height(eb)
    scale = 1.0 / sqrt(w)
    if recentre
        alt = alt_x(h, w)
        data =  alt .* eb.vals
        eb.fft ? fft!(data, 2) : bfft!(data, 2)
        Eigenbrot(scale * (alt .* data), !eb.fft)
    else
        data = copy(eb.vals)
        eb.fft ? fft!(data, 2) : bfft!(data, 2)
        Eigenbrot(scale * data, !eb.fft)
    end
end

"""
    ffty(eb[, recentre = true])
Returns the one-dimensional DFT of Eigenbrot `eb` in
the `y`-direction. If `recentre` is `true` the DFT is
stored so that the value at the coordinate origin is
at the central location of the data.
"""
function ffty(eb::Eigenbrot, recentre::Bool = true)
    w = width(eb)
    h = height(eb)
    scale = 1.0 / sqrt(h)
    if recentre
        alt = alt_y(h, w)
        data = alt .* eb.vals
        eb.fft ? fft!(data, 1) : bfft!(data, 1)
        Eigenbrot(scale * (alt .* data), !eb.fft)
    else
        data = copy(eb.vals)
        eb.fft ? fft!(data, 1) : bfft!(data, 1)
        Eigenbrot(scale * data, !eb.fft)
    end
end

"""
    removeDC(eb)
Return an Eigenbrot that is the Eigenbrot `eb`
with the zero-frequency component removed.
"""
function removeDC(eb::Eigenbrot)
  r = fft(eb, false)
  r[1, 1] = 0
  return fft(r, false)
end


text:
function circ_get(str::String, index::Integer)
  while index > length(str)
    index = index - length(str)
  end
  return str[index]
end

# Read input
input = readline(STDIN)

## Part 1
# Sum
sum = 0

# Find pairs
for (i, c) in enumerate(input)
  if c == circ_get(input, i + 1)
    sum = sum + parse(Int64, c)
  end
end

println("Part 1: $sum");

## Part 2
sum = 0
for (i, c) in enumerate(input)
  if c == circ_get(input, i + convert(Integer, length(input) / 2))
    sum = sum + parse(Int64, c)
  end
end

println("Part 2: $sum");


text:
# API definitions for builtin types

# Integers
left(i::Integer) = i > 0 ? Singleton(i - one(i)) : ∅
right(i::Integer) = i < 0 ? Singleton(i + one(i)) : ∅
generation(i::Integer) = abs(i)

show_inner(io::IO, ::MIME"text/latex", i::Integer) = print(io, i)


text:
using SeisIO, Printf, PlotlyJS, JLD2, FileIO, Statistics, ORCA

#----------------#
finame = "./dataset/NZ_GEONETquake.jld2"

#eventID = 2 # choose station and components from stationlist in SeisData
channel = "BHZ" # choose station and components from stationlist in SeisData
#----------------#

# load data
c = jldopen(finame)
event=c["info/event"]

for eventID = 1:length(event)

	origin = event[eventID]["origin"]
	lat1 = parse(Float64, origin["latitude"])
	lon1 = parse(Float64, origin["longitude"])
	time = origin["time"][1:end-4]
	mag = parse(Float64,origin["magnitude"])

	# search channel
	S = c["event/$time"]
	Sall = SeisData()
	for i = 1:size(S)[1]
		Stemp = S[i]

		#apply filters
		detrend!(Stemp)
		filtfilt!(Stemp, fl=0.1, fh=2.0)
 		taper!(Stemp)
		if occursin(channel, Stemp.id)
			push!(Sall, Stemp)
		else
			continue;
		end
	end

	# compute angular distance

	# chose plotting data
	numoftrace = size(Sall)[1]

	# for sorting plot order
	plotorder = zeros(numoftrace, 2)
	for i = 1:numoftrace
		lat2 = Sall[i].loc.lat
		lon2 = Sall[i].loc.lon

		Δ = acosd(sind(lat1)*sind(lat2) + cosd(lat1)*cosd(lat2)*cosd(lon1 - lon2))
	 	Sall[i].misc["angulardist"] = Δ
		plotorder[i,:] = [float(i), Δ]
	end

	#plotting

	tvec = collect(0:Sall[1].t[end,1]) ./ Sall[1].fs

	trace1 = Array{GenericTrace{Dict{Symbol,Any}}, 1}(undef, numoftrace)
	traceid = 0

	dist_α = 1.0

	plotidlist = round.(Int64, sortslices(plotorder, dims=1, lt=(x,y) -> isless(x[2], y[2]), rev=true)[:,1])
	for i = plotidlist
		normalized_amp = maximum(Sall[plotidlist[end]].x)
		traceid += 1
		trace1[traceid] = scatter(;x=tvec, y=Sall[i].x ./ normalized_amp .+ dist_α*Sall[i].misc["angulardist"], mode="lines",
	     name=@sprintf("%s", Sall[i].id))
	end

	layout = Layout(width=800, height=600,
					xaxis=attr(title="Time [sec]"),
					yaxis=attr(title="Angular distance [°]"),
	                font =attr(size=18),
	                showlegend=true,
					xaxis_range=[0, 180],
					yaxis_range=[-0.5, 8],
	                title = @sprintf("M:%4.1f Time:%s Lat:%4.2f Lon:%4.2f", mag ,time, lat1, lon1)
					)

	p = plot(trace1, layout)
	display(p)

	mkpath("./fig")
	figname = @sprintf("./fig/%s_M%4.2f.png", time, mag)

	savefig(p, figname)

	println("Press enter to display next event:")
	readline()

end

println("to finish plotting: press any key...")
readline()


text:
# By Lazaro Alonso
using GLMakie, ColorSchemes, Random
GLMakie.activate!() # HIDE
function plotVolColormaps(; rseed = 123)
    Random.seed!(rseed)
    # all colormaps from colorschemes
    cs = collect(keys(colorschemes))
    cmapIdx = Node(1)
    cmap = @lift(cs[$cmapIdx]) # @lift(Reverse(cs[$cmapIdx]))
    ncolors = 101
    colors = @lift(to_colormap($cmap, ncolors))
    # some transparencies in the colors
    g(x) = 0.5 + 0.5*tanh((x+3)/3)
    alphas = g.(LinRange(-10,10, ncolors))
    cmap_alpha = @lift(RGBAf0.($colors, alphas))
    # the actual figure
    n = 15
    x = randn(n)
    y = randn(n)
    z = randn(n)
    A = randn(n, n, n) ./ 3

    fig = Figure(resolution = (800,800))
    ax = Axis3(fig, aspect = (1,1,1), perspectiveness = 0.5, elevation = π/9)
    volume!(ax, x, y ,z, A, colormap = cmap_alpha, transparency =false)
    sl = Slider(fig[1, 2], range = 1:length(cs), startvalue = 400, horizontal = false)
    connect!(cmapIdx, sl.value)
    fig[1,1] = ax
    fig[0,1] = Label(fig, @lift("Colormap: $(cs[$cmapIdx])"), textsize = 20,
        tellheight = true, tellwidth = false)
    fig
end
fig = with_theme(plotVolColormaps, theme_black())
save(joinpath(@__DIR__, "output", "plotVolColormaps.png"), fig, px_per_unit = 2.0) # HIDE
using Pkg # HIDE
Pkg.status(["GLMakie", "ColorSchemes"]) #HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
module TestPipeline

using Random
using Test
using AMLPipelineBase
using AMLPipelineBase.Pipelines
using AMLPipelineBase.DecisionTreeLearners
using AMLPipelineBase.CrossValidators
using AMLPipelineBase.Utils

const data     = getiris()
const features = data[:,1:4]
const X        = data[:,1:5]
const Y               = data[:,5] |> Vector
X[!,5]                = X[!,5] .|> string

const ohe = OneHotEncoder()
const noop = Identity()
const rf  = RandomForest()
const ada = Adaboost()
const pt  = PrunedTree()

function test_pipeline()
  # test initialization of types
  ohe = OneHotEncoder()
  linear1 = Pipeline(Dict(:name=>"lp",:machines => [ohe]))
  linear2 = Pipeline(Dict(:name=>"lp",:machines => [ohe]))
  combo1 = ComboPipeline(Dict(:machines=>[ohe,ohe]))
  combo2 = ComboPipeline(Dict(:machines=>[linear1,linear2]))
  linear1 = Pipeline([ohe])
  linear2 = Pipeline([ohe])
  combo1 = ComboPipeline([ohe,ohe])
  combo2 = ComboPipeline([linear1,linear2])
  # test fit/transform workflow
  fit!(combo1,X)
  res1=transform!(combo1,X)
  res2=fit_transform!(combo1,X)
  @test (res1 .== res2) |> Matrix |> sum == 2100
  fit!(combo2,X)
  res3=transform!(combo2,X)
  res4=fit_transform!(combo2,X)
  @test (res3 .== res4) |> Matrix |> sum == 2100
  pcombo1 = @pipeline ohe + ohe
  pres1 = fit_transform!(pcombo1,X)
  @test (pres1 .== res1) |> Matrix |> sum == 2100
end
@testset "Pipelines" begin
  Random.seed!(123)
  test_pipeline()
end

acc(X,Y) = score(:accuracy,X,Y)

function test_sympipeline()
  pcombo5 = @pipeline :((ohe + noop) |> (ada * rf * pt))
  @test crossvalidate(pcombo5,X,Y,acc,5,false).mean >= 0.90
  expr = :((ohe + noop) |> (ada * rf * pt))
  processexpr!(expr.args)
  @test crossvalidate(eval(expr),X,Y,acc,5,false).mean >= 0.90
  pcombo6 = sympipeline(expr) |> eval
  @test crossvalidate(pcombo6,X,Y,acc,5,false).mean >= 0.90
  pcombo7 = (@pipelinez expr) |> eval
  @test crossvalidate(pcombo7,X,Y,acc,5,false).mean >= 0.90
end
@testset "Symbolic Pipeline: Global Scope" begin
  Random.seed!(123)
  test_sympipeline()
end

function test_pipeline()
  # test symbolic pipeline expression
  pcombo2 = @pipeline ohe + noop
  @test fit_transform!(pcombo2,features) |> Matrix |> size |> collect |> sum == 158
  pcombo2 = @pipeline ohe + noop |> rf
  @test crossvalidate(pcombo2,X,Y,acc,5,false).mean >= 0.90
end
@testset "Symbolic Pipeline: Local Scope" begin
  Random.seed!(123)
  test_pipeline()
end

end


text: # TODO: PaddedString


text:
export rosenbrock, rastrigin

"""
The rastrigin function is notorious for having multiple local minima and only one global optima at (0,0)

"""
function rastrigin(x::Vector{T}) where {T <: AbstractFloat}
    a = one(eltype(x)) * 10
    return size(x, 1) * a + sum(@. x^2 - a * cos(2 * π * x))
end

function rosenbrock(x::Vector{T}) where {T <: AbstractFloat}
    f = zero(eltype(x))
    a = one(eltype(x))
    b = 100 * a
    for i in 1:(length(x) - 1)
       @inbounds f += b * (x[i + 1] - x[i]^2)^2 + (a - x[i])^2
    end
    return f
end


text:
struct HealpyNotImported end
getproperty(::HealpyNotImported, ::Symbol) = error("Healpy needs to be installed to use this function.")
@init try
    global hp = pyimport(:healpy)
catch
    global hp = HealpyNotImported()
end

# use generated function to memoize ring info
ringinfo(Nside) = ringinfo(Val(Nside))
@generated function ringinfo(::Val{Nside}) where {Nside}
    r = hp.ringinfo(Nside, collect(1:4Nside-1))
    ring_lengths = r[2]
    ring_starts = [1; cumsum(ring_lengths)[1:end-1] .+ 1]
    ring_ranges = [range(a,length=b) for (a,b) in tuple.(ring_starts, ring_lengths)]
    (ring_lengths=ring_lengths, ring_starts=ring_starts, ring_ranges=ring_ranges, cosθ=r[3], sinθ=r[4], θ=acos.(r[3]))
end

# GradientCache holds precomputed quantities that let us do fast spherical gradients
struct GradientCache{Nside, T, Nobs, NB, W}
    neighbors        :: NB
    W_covariant      :: W
    W_contravariant  :: W
    Wᵀ_covariant     :: W
    Wᵀ_contravariant :: W

    function GradientCache{Nside,T}(Nobs, order=Val(1)) where {Nside,T}
        N_coeffs    = (order == Val(1)) ? 3 : 6
        N_neighbors = (order == Val(1)) ? 4 : 8

        # the commented out line uses the pixel itself, but that gives us a non-zero trace so we don't use it
        # neighbors_mat = [(0:(Nobs-1))'; hp.get_all_neighbours(Nside,collect(0:(Nobs-1)))[(order == Val(1) ? (1:2:end) : 1:end),:]::Matrix{Int}] .+ 1
        neighbors_mat = hp.get_all_neighbours(Nside,collect(0:(Nobs-1)))[(order == Val(1) ? (2:2:end) : 1:end),:]::Matrix{Int} .+ 1
        # technically we dont need the +1 below, but this makes it so the entire
        # very last ring is stored (for uniformity's sake)
        Ntot = maximum(neighbors_mat)+1
        neighbors = SVector{N_neighbors}.(eachcol(Int32.(neighbors_mat)))

        (θs, ϕs) = convert.(Vector{T}, hp.pix2ang(Nside,collect(0:Ntot))::Tuple{Vector{Float64},Vector{Float64}})

        # derivative
        W_covariant, W_contravariant = [], []
        @showprogress 1 "∇  precomputation: " for (i,(ni,θ,ϕ)) in enumerate(zip(neighbors, θs, ϕs))
            if all(ni .<= Nobs)
                Δθ = @.     θs[ni]-θ
                Δϕ = @. rem(ϕs[ni]-ϕ+T(π), T(2π), RoundDown) - T(π)
                if order == Val(1)
                    P = @. [Δθ Δϕ ones(T,N_neighbors)]
                else
                    P = @. [Δθ Δϕ ones(T,N_neighbors) Δθ^2 Δϕ^2 Δθ*Δϕ]
                end
                Q,R = qr(P)
                W = inv(R)[1:2,:]*Q'
            else
                W = zeros(T, 2, N_neighbors)
                neighbors[i] = map(j->(j>Nobs ? i : j), neighbors[i])
            end
            push!(W_covariant,     SMatrix{2,N_neighbors}(W))
            push!(W_contravariant, SMatrix{2,N_neighbors}(W ./ [1, sin(θ)^2]))
        end
        _W_covariant     = collect(typeof(W_covariant[1]),    W_covariant)
        _W_contravariant = collect(typeof(W_contravariant[1]),W_contravariant)

        # transpose derivative
        Wᵀ_covariant     = fill(T(0), length(neighbors), 2, N_neighbors)
        Wᵀ_contravariant = fill(T(0), length(neighbors), 2, N_neighbors)
        @showprogress 1 "∇' precomputation: " for (i,Ni) in collect(enumerate(neighbors))
            for (j,Nij) in enumerate(Ni)
                if Nij<=length(neighbors)
                    j′ = first(indexin(i,neighbors[Nij]))
                    if j′ != nothing
                        Wᵀ_covariant[i,:,j]     = _W_covariant[Nij][:,j′]
                        Wᵀ_contravariant[i,:,j] = _W_contravariant[Nij][:,j′]
                    end
                end
            end
        end
        _Wᵀ_covariant     = SMatrix{2,N_neighbors,T}.(@views [Wᵀ_covariant[i,:,:]     for i=1:length(neighbors)])
        _Wᵀ_contravariant = SMatrix{2,N_neighbors,T}.(@views [Wᵀ_contravariant[i,:,:] for i=1:length(neighbors)])

        new{Nside,T,Nobs,typeof(neighbors),typeof(_W_covariant)}(neighbors, _W_covariant, _W_contravariant, _Wᵀ_covariant, _Wᵀ_contravariant)
    end

end
function GradientCache(Nside, T; θmax, order=1)
    i = hp.ang2pix(Nside, deg2rad(θmax), 0)
    Nobs = ringinfo(Nside).ring_starts[findfirst(j->j>i, ringinfo(Nside).ring_starts)] - 1
    GradientCache{Nside,T}(Nobs, Val(order))
end

# now define the Healpix Fields (called HealpixCap)

struct HpxPix{Nside} <: Pix end
abstract type HealpixCap{Nside,T,Nobs,B,S,P<:HpxPix} <: Field{B,S,P} end

## Spin-0
struct HealpixS0Cap{Nside, T, Nobs, GC<:Union{Nothing,GradientCache{Nside,T,Nobs}}} <: HealpixCap{Nside, T, Nobs, Map, S0, HpxPix{Nside}}
    Ix::Vector{T}
    gradient_cache::GC
end
struct HealpixS2Cap{Nside, T, Nobs, GC<:Union{Nothing,GradientCache{Nside,T,Nobs}}} <: HealpixCap{Nside, T, Nobs, QUMap, S2, HpxPix{Nside}}
    QUx::Matrix{T}
    gradient_cache::GC
end
function (::Type{H})(x::VecOrMat{T}; gradient_order=nothing, Nside=hp.npix2nside(size(x,1))) where {T, H<:HealpixCap}
    Nobs = maximum(findall(!isnan,x[:,1]))
    gc = gradient_order==nothing ? nothing : GradientCache{Nside,T}(Nobs,Val(gradient_order))
    H{Nside,T,Nobs,typeof(gc)}(x, gc)
end

getproperty(f::HealpixS2Cap, ::Val{:Qx}) = f.QUx[:,1]
getproperty(f::HealpixS2Cap, ::Val{:Ux}) = f.QUx[:,2]

# convertable_fields is screwed up for HealpixCap so say by hand here there's
# none (todo: fix this)
convertable_fields(::Type{F}) where {B,S,P<:HpxPix,F<:Field{B,S,P}} = []

similar(f::F) where {F<:HealpixCap} = F(similar(first(fieldvalues(f))), f.gradient_cache)
copy(f::F)    where {F<:HealpixCap} = F(copy(first(fieldvalues(f))),    f.gradient_cache)

## derivatives

DerivBasis(::Type{<:HealpixS0Cap}) = Map
DerivBasis(::Type{<:HealpixS2Cap}) = QUMap

# picking out the right weight matrix for a given ∇ operator
get_W(∇::AbstractArray, gc::GradientCache) = get_W(first(∇), gc)
get_W(::∇i{<:Any,true},           gc) = gc.W_covariant
get_W(::∇i{<:Any,false},          gc) = gc.W_contravariant
get_W(::AdjOp{<:∇i{<:Any,true}},  gc) = gc.Wᵀ_covariant
get_W(::AdjOp{<:∇i{<:Any,false}}, gc) = gc.Wᵀ_contravariant
get_W(::Any, gc::Nothing) = error("Can't take gradients of this field, gradient cache not precomputed.")

function mul!(∇f::FieldVector{F}, ∇Op::∇Op, f::F) where {Nside,T,Nobs,F<:HealpixS0Cap{Nside,T,Nobs}}
    gc = f.gradient_cache
    W = get_W(∇Op, gc)
    @inbounds for i in eachindex(gc.neighbors)
        Ix = @view f.Ix[gc.neighbors[i]]
        ∇f[1].Ix[i], ∇f[2].Ix[i] = W[i] * Ix
    end
    ∇f
end
function mul!(∇f::FieldVector, ∇Op::Union{∇Op,Adjoint{∇i,∇Op}}, f::HealpixS2Cap)
    gc = f.gradient_cache
    W = get_W(∇op, gc)
    @inbounds for i in eachindex(gc.neighbors)
        Qx = @view f.QUx[gc.neighbors[i],1]
        Ux = @view f.QUx[gc.neighbors[i],2]
        ∇f[1].Qx[i,1], ∇f[2].Qx[i,1] = W[i] * Qx
        ∇f[1].Ux[i,2], ∇f[2].Ux[i,2] = W[i] * Ux
    end
    imax = gc.neighbors[end][1] + 1
    ∇f[1].Qx[imax:end] .= ∇f[2].Qx[imax:end] .= ∇f[1].Ux[imax:end] .= ∇f[2].Ux[imax:end] .= NaN
    ∇f
end
*(∇Op::Union{∇Op,Adjoint{∇i,<:∇Op}}, f::HealpixCap) where {B} =  mul!(allocate_result(∇Op,f),∇Op,f)

*(∇Op::Adjoint{∇i,<:∇Op}, v::FieldVector{<:HealpixS0Cap}) = mul!(similar(v[1]), ∇Op, v)
function mul!(f′::F, ∇Op::Adjoint{∇i,<:∇Op}, v::FieldVector{F}, memf′::F=v[1]) where {Nside,T,Nobs,F<:HealpixS0Cap{Nside,T,Nobs}}
    gc = f′.gradient_cache
    W = get_W(∇Op, gc)
    for i in eachindex(gc.neighbors)
        f′.Ix[i] = -(  (W[i] * @view v[1].Ix[gc.neighbors[i]])[1]
                     + (W[i] * @view v[2].Ix[gc.neighbors[i]])[2])
    end
    f′.Ix[Nobs+1:end] .= NaN
    f′
end

# individual components of the gradient (which is wasteful), but needed for the
# way `negδvelocityᴴ!` is currently written. todo: remove this and change the
# way that that's written instead.
function mul!(f′::F, ∇Op::AdjOp{<:∇i{component}}, f::F) where {component,Nside,T,Nobs,F<:HealpixS0Cap{Nside,T,Nobs}}
    if f′===f; f = copy(f); end # even more wasteful...
    gc = f′.gradient_cache
    W = get_W(∇Op, gc)
    for i in eachindex(gc.neighbors)
        f′.Ix[i] = -(W[i] * @view(f.Ix[gc.neighbors[i]]))[component+1]
    end
    f′
end

dot(a::H, b::H) where {Nside, H<:HealpixS0Cap{Nside}} = dot(nan2zero.(a.Ix), nan2zero.(b.Ix))  * hp.nside2pixarea(Nside)
dot(a::H, b::H) where {Nside, H<:HealpixS2Cap{Nside}} = dot(nan2zero.(a.QUx),nan2zero.(b.QUx)) * hp.nside2pixarea(Nside)

function plot(f::HealpixS0Cap, args...; cmap="RdBu_r", vlim=nothing, plot_type=(PyPlot.isinteractive() ? :mollzoom : :mollview), kwargs...)
    kwargs = Dict(kwargs...)
    cmap = get_cmap(cmap)
    if vlim!=nothing
        kwargs[:min], kwargs[:max] = -vlim, vlim
    end
    getproperty(hp,plot_type)(full(f), args...; cmap=cmap, kwargs...)
end

function plot(f::HealpixS2Cap, args...; kwargs...)
    plot(HealpixS0Cap(f.QUx[:,1], f.gradient_cache), args...; kwargs...)
    plot(HealpixS0Cap(f.QUx[:,2], f.gradient_cache), args...; kwargs...)
end

## conversion to flat sky maps
function azeqproj(f::HealpixS0Cap{<:Any,T}, θpix, Nside, lamb=false) where {T}
    pylab = pyimport("matplotlib.pylab")
    wasinteractive = pylab.isinteractive()
    try
        pylab.ioff()
        Ix = hp.azeqview(full(f), rot=(0,90), reso=θpix, xsize=Nside, lamb=lamb, return_projected_map=true)
        close()
        FlatS0Map{T,Flat{θpix,Nside,fourier∂}}(Ix)
    finally
        wasinteractive && pylab.ion()
    end
end
function azeqproj(f::HealpixS2Cap, θpix, Nside)
    FlatS2QUMap((azeqproj(HealpixS0Cap(getproperty(f,k), f.gradient_cache), θpix, Nside) for k in (:Qx,:Ux))...)
end
azeqproj(f::HealpixCap{Nside,T,Nobs}) where {Nside,T,Nobs} = azeqproj(f, round(600rad2deg(hp.nside2resol(Nside)))/10, Nside)

## broadcasting
broadcast_data(::Type{F}, f::F) where {F<:HealpixS0Cap} = (f.Ix,)
broadcast_data(::Type{F}, f::F) where {F<:HealpixS2Cap} = (f.QUx,)
metadata(::Type{F}, f::F) where {F<:HealpixCap} = (f.gradient_cache,)
metadata_reduce((m1,)::Tuple{GC}, (m2,)::Tuple{GC}) where {GC<:Union{Nothing,GradientCache}} = (m1,)

LenseBasis(::Type{<:HealpixS0Cap}) = Map
LenseBasis(::Type{<:HealpixS2Cap}) = QUMap

adjoint(f::HealpixS0Cap) = f
length(::Type{<:HealpixS0Cap{Nside,T,Nobs}}) where {Nside,T,Nobs} = Nobs

function full(f::HealpixS0Cap{Nside,T}) where {Nside,T}
    Ix = fill(T(NaN),12*Nside^2)
    Ix[1:length(f.Ix)] .= f.Ix
    Ix
end

## Harmonic Operators

struct IsotropicHarmonicCov{Nside, T, Nobs, N, GC<:GradientCache{Nside, T, Nobs}} <: LinOp{Basis, Spin, HpxPix}
    Cℓ :: Array{T,N}
    gc :: GC
end
IsotropicHarmonicCov(Cℓ::Array,  gc::GradientCache{Nside,T}) where {Nside,T} = IsotropicHarmonicCov(T.(Cℓ), gc)

*(L::BandPassOp, f::HealpixCap{Nside, T}) where {Nside, T} = IsotropicHarmonicCov(T.(L.Wℓ), f.gradient_cache) * f
function mul!(f′::F, L::IsotropicHarmonicCov{Nside, T, Nobs}, f::F) where {Nside, T, Nobs, F<:HealpixCap{Nside, T, Nobs}}
    ℓmax = min(size(L.Cℓ,1)-1, 3Nside)
    aℓms = map2alm(f, ℓmax=ℓmax)
    @. aℓms *= L.Cℓ[1:ℓmax+1,:]'
    alm2map!(f′, aℓms)
end
function ldiv!(f′::F, L::IsotropicHarmonicCov{Nside, T, Nobs}, f::F) where {Nside, T, Nobs, F<:HealpixCap{Nside, T, Nobs}}
    ℓmax = min(size(L.Cℓ,1)-1, 3Nside)
    aℓms = map2alm(f, ℓmax=ℓmax)
    @. aℓms = nan2zero(L.Cℓ[1:ℓmax+1,:]' \ aℓms)
    alm2map!(f′, aℓms)
end

# one ring before the ring on which the first unobserved pixel is
get_zbounds(Nside, Nobs) = [ringinfo(Nside).cosθ[hp.pix2ring(Nside,[Nobs])[1]], 1]

function map2alm(f::HealpixCap{Nside,T,Nobs}; ℓmax=2Nside) where {Nside,T,Nobs}
    zbounds = get_zbounds(Nside, Nobs)
    map2alm(first(fieldvalues(f)), Nside=Nside, ℓmax=ℓmax, zbounds=zbounds)
end

function alm2map!(f::HealpixCap{Nside,T,Nobs}, aℓms::Array{Complex{T},3}) where {Nside,T,Nobs}
    zbounds = get_zbounds(Nside, Nobs)
    alm2map!(first(fieldvalues(f)), aℓms, Nside=Nside, zbounds=zbounds)
    first(fieldvalues(f))[Nobs+1:end,:] .= NaN
    f
end

# todo: implement proper broadcasting
+(Σa::I, Σb::I) where {I<:IsotropicHarmonicCov} = IsotropicHarmonicCov(Σa.Cℓ.+Σb.Cℓ, Σa.gc)
*(Σa::I, Σb::I) where {I<:IsotropicHarmonicCov} = IsotropicHarmonicCov(Σa.Cℓ.*Σb.Cℓ, Σa.gc)
*(Σ::IsotropicHarmonicCov, α::Real) = IsotropicHarmonicCov(α*Σ.Cℓ, Σ.gc)
*(α::Real, Σ::IsotropicHarmonicCov) = IsotropicHarmonicCov(α*Σ.Cℓ, Σ.gc)
inv(Σ::IsotropicHarmonicCov) = IsotropicHarmonicCov(nan2zero.(inv.(Σ.Cℓ)), Σ.gc)
sqrt(Σ::IsotropicHarmonicCov) = IsotropicHarmonicCov(sqrt.(Σ.Cℓ), Σ.gc)
simulate(Σ::IsotropicHarmonicCov{Nside,T,Nobs,1}) where {Nside,T,Nobs} =
    sqrt(Σ) * HealpixS0Cap(randn(T,Nobs)/T(hp.nside2resol(Nside)), Σ.gc)
simulate(Σ::IsotropicHarmonicCov{Nside,T,Nobs,2}) where {Nside,T,Nobs} =
    sqrt(Σ) * HealpixS2Cap(randn(T,Nobs,2)/T(hp.nside2resol(Nside)), Σ.gc)
zero(Σ::IsotropicHarmonicCov{Nside,T,Nobs}) where {Nside,T,Nobs} =
    HealpixS0Cap(zeros(T,Nobs), Σ.gc)

##
function HealpixCapMask(θmax, Δθapod, Nside)
    W(θ) = θ<θmax ? 1 : θmax<θ<(θmax+Δθapod) ? (cos(π*(θ-θmax)/Δθapod)+1)/2 : 0
    i,l = hp.ringinfo(Nside, collect(1:4Nside))[1:2]
    m = zeros(12*Nside^2)
    for (r,w) in zip(broadcast(:, i, i .+ l .- 1), W.(rad2deg.(acos.(hp.ringinfo(Nside, collect(1:4nside))[3]))))
        m[r.+1] .= w
    end
    m
end

##

# Healpy doesn't yet have the zbounds option, so call into the Healpix Fortran
# libraries directly.

# needing to manually dlopen this is probably a bug in Healpix or gfortan, not
# sure which...
using Libdl
@init try
    Libdl.dlopen("libgomp",Libdl.RTLD_GLOBAL)
catch
    @warn "Failed to load libgomp. Healpix support may not work."
end

@generated function map2alm(maps::Array{T,N}; Nside=hp.npix2nside(size(maps,1)), ℓmax=2Nside, mmax=ℓmax, zbounds=[-1,1]) where {N,T<:Union{Float32,Float64}}
    (spin, Tspin) = if (N==1)
        (), ()
    elseif (N==2)
        (2,), (Ref{Int32},)
    else
        error("maps should be Npix-×-1 or Npix-×-2")
    end
    fn_name = "__alm_tools_MOD_map2alm_$(N==1 ? "sc" : "spin")_$((T==Float32) ? "s" : "d")"
    quote
        aℓms = Array{Complex{T}}(undef, ($N,ℓmax+1,mmax+1))
        aℓms .= NaN
        ccall(
            ($fn_name, "libhealpix"), Nothing,
            (Ref{Int32}, Ref{Int32}, Ref{Int32}, $(Tspin...), Ref{T}, Ref{Complex{T}}, Ref{Float64}, Ref{Nothing}),
            Nside, ℓmax, mmax, $(spin...), maps, aℓms, Float64.(zbounds), C_NULL
        )
        aℓms
    end
end

@generated function alm2map!(maps::Array{T,N}, aℓms::Array{Complex{T},3}; Nside=hp.npix2nside(size(maps,1)), ℓmax=(size(aℓms,2)-1), mmax=(size(aℓms,3)-1), zbounds=[-1,1]) where {N,T<:Union{Float32,Float64}}
    (spin, Tspin) = if (N==1)
        (), ()
    elseif (N==2)
        (2,), (Ref{Int32},)
    else
        error("maps should be Npix-×-1 or Npix-×-2")
    end
    fn_name = "__alm_tools_MOD_alm2map_$(N==1 ? "sc_wrapper" : "spin")_$((T==Float32) ? "s" : "d")"
    quote
        ccall(
           ($fn_name, "libhealpix"), Nothing,
           (Ref{Int32}, Ref{Int32}, Ref{Int32}, $(Tspin...), Ref{Complex{T}}, Ref{T}, Ref{Float64}, Ref{Nothing}),
           Nside, ℓmax, mmax, $(spin...), aℓms, maps, Float64.(zbounds), C_NULL
        )
        maps
    end
end

function alm2cl(alm::Matrix{Complex{T}}; ℓmax=(size(alm,1)-1)) where {T}
    InterpolatedCℓs(0:ℓmax, [(abs2(alm[ℓ+1,1]) + 2sum(abs2.(alm[ℓ+1, 2:ℓ+1])))/(2ℓ+1) for ℓ=0:ℓmax])
end

function alm2cl(alm1::Matrix{Complex{T}}, alm2::Matrix{Complex{T}}; ℓmax=(size(alm,1)-1)) where {T}
    InterpolatedCℓs(0:ℓmax, real.([(dot(alm1[ℓ+1,1], alm2[ℓ+1,1]) + 2dot(alm1[ℓ+1,2:ℓ+1], alm2[ℓ+1,2:ℓ+1]))/(2ℓ+1) for ℓ=0:ℓmax]))
end

get_Cℓ(mp::HealpixS0Cap{Nside}; ℓmax=2Nside) where {Nside} = alm2cl(map2alm(mp, ℓmax=ℓmax))
get_Cℓ(mp1::HealpixS0Cap{Nside}, mp2::HealpixS0Cap{Nside}; ℓmax=2Nside) where {Nside} = alm2cl(map2alm(mp1, ℓmax=ℓmax), map2alm(mp2, ℓmax=ℓmax), ℓmax=ℓmax)


text:
# represents s*log(ε) + c as ε -> 0
struct LogNumber <: Number
    s::ComplexF64
    c::ComplexF64
end

@inline logpart(z::Number) = zero(z)
@inline finitepart(z::Number) = z

@inline logpart(l::LogNumber) = l.s
@inline finitepart(l::LogNumber) = l.c

Base.promote_rule(::Type{LogNumber}, ::Type{<:Number}) = LogNumber
Base.convert(::Type{LogNumber}, z::LogNumber) = z
Base.convert(::Type{LogNumber}, z::Number) = LogNumber(0, z)

==(a::LogNumber, b::LogNumber) = logpart(a) == logpart(b) && finitepart(a) == finitepart(b)
Base.isapprox(a::LogNumber, b::LogNumber; opts...) = ≈(logpart(a), logpart(b); opts...) && ≈(finitepart(a), finitepart(b); opts...)

(l::LogNumber)(ε) = logpart(l)*log(ε) + finitepart(l)

for op in (:real, :imag, :conj)
    @eval $op(l::LogNumber) = LogNumber($op(logpart(l)), $op(finitepart(l)))
end

for f in (:+, :-)
    @eval begin
        $f(a::LogNumber, b::LogNumber) = LogNumber($f(a.s, b.s), $f(a.c, b.c))
        $f(l::LogNumber, b::Number) = LogNumber(l.s, $f(l.c, b))
        $f(a::Number, l::LogNumber) = LogNumber($f(l.s), $f(a, l.c))
    end
end

-(l::LogNumber) = LogNumber(-l.s, -l.c)

for Typ in (:Bool, :Number)
    @eval begin
        *(l::LogNumber, b::$Typ) = LogNumber(l.s*b, l.c*b)
        *(a::$Typ, l::LogNumber) = LogNumber(a*l.s, a*l.c)
    end
end
/(l::LogNumber, b::Number) = LogNumber(l.s/b, l.c/b)

function exp(l::LogNumber)::ComplexF64
    if real(l.s) > 0
        0.0+0.0im
    elseif real(l.s) < 0
        Inf+Inf*im
    elseif real(l.s) == 0 && imag(l.s) == 0
        log(l.c)
    else
        NaN + NaN*im
    end
end

# This is a relative version of dual number, in the sense that its realpart*(1+epsilon)
struct RiemannDual{T} <: Number
    realpart::T
    epsilon::T
end

RiemannDual(x, y) = RiemannDual(promote(x, y)...)

RiemannDual(x::Dual) = RiemannDual(realpart(x), epsilon(x))
Dual(x::RiemannDual) = Dual(realpart(x), epsilon(x))
dual(x::RiemannDual) = Dual(x)

# the relative perturbation
realpart(r::RiemannDual) = r.realpart
epsilon(r::RiemannDual) = r.epsilon
undirected(r::RiemannDual) = undirected(realpart(r))
isinf(r::RiemannDual) = isinf(realpart(r))

in(x::RiemannDual, d::Domain) = realpart(x) ∈ d
in(x::RiemannDual, d::TypedEndpointsInterval{:closed,:closed}) = realpart(x) ∈ d

for f in (:-,)
    @eval $f(x::RiemannDual) = RiemannDual($f(realpart(x)),$f(epsilon(x)))
end

for f in (:sqrt,)
    @eval $f(x::RiemannDual) = RiemannDual($f(dual(x)))
end

for f in (:+, :-, :*)
    @eval $f(x::RiemannDual, y::RiemannDual) = RiemannDual($f(dual(x),dual(y)))
    for Typ in (:Bool, :Number)
        @eval begin
            $f(x::RiemannDual, p::$Typ) = RiemannDual($f(dual(x),p))
            $f(p::$Typ, x::RiemannDual) = RiemannDual($f(p,dual(x)))
        end
    end
end

for OP in (:*, :+, :-, :/)
    @eval begin
        $OP(a::Directed{s}, b::RiemannDual) where {s} = Directed{s}($OP(a.x,b))
        $OP(a::RiemannDual, b::Directed{s}) where {s} = Directed{s}($OP(a,b.x))
    end
end

function inv(z::RiemannDual)
    tol = 1E-15
    abs(realpart(z)) ≤ tol && return RiemannDual(Inf, inv(epsilon(z)))
    RiemannDual(inv(dual(z)))
end

/(z::RiemannDual, x::RiemannDual) = z*inv(x)
/(z::RiemannDual, x::Number) = z*inv(x)
/(x::Number, z::RiemannDual) = x*inv(z)

# loses sign information
for f in (:real, :imag, :abs)
    @eval $f(z::RiemannDual) = $f(realpart(z))
end

function log(z::RiemannDual)
    tol = 1E-14
    abs(realpart(z)) ≤ tol  || abs(realpart(z)) ≥ 1/tol || throw(ArgumentError("Cannot eval at $z"))
    LogNumber(abs(realpart(z)) ≤ tol ? 1 : -1, log(abs(epsilon(z))) + im*angle(epsilon(z)))
end

function atanh(z::RiemannDual)
    if realpart(z) ≈ 1
        LogNumber(-0.5,log(2)/2  - log(abs(epsilon(z)))/2 - im/2*angle(-epsilon(z)))
    elseif realpart(z) ≈ -1
        LogNumber(0.5,-log(2)/2  + log(abs(epsilon(z)))/2 + im/2*angle(epsilon(z)))
    else
        error("Not implemented")
    end
end

log1p(z::RiemannDual) = log(z+1)

SingularIntegralEquations.HypergeometricFunctions.speciallog(x::RiemannDual) =
    (s = sqrt(x); 3(atanh(s)-realpart(s))/realpart(s)^3)

Base.show(io::IO, x::RiemannDual) = show(io, Dual(x))
Base.show(io::IO, x::LogNumber) = print(io, "($(logpart(x)))log ε + $(finitepart(x))")

# # (s*log(M) + c)*(p*M
# function /(l::LogNumber, b::RiemannDual)
#     @assert isinf(realpart(b))
#     LogNumber(l.s/b, l.c/b)

mxa_₂F₁(a, b, c, z::RiemannDual) = directed_mxa_₂F₁(a,b,c,z)
_₂F₁general(a::Number, b::Number, c::Number, z::RiemannDual) = directed_₂F₁general(a,b,c,z)


text:
# single qubit Cliffords

export localclifford

const pX = [0 1; 1 0]
const pY = [0 -im; im 0]
const pZ = [1 0; 0 -1]

C1 = Vector{Clifford}(24)

# identity
C1[1]  = RI

# pi/2 and pi rotations about X,Y,Z
C1[2]  = expm(-im*pi/4*pX)
C1[3]  = im * pX
C1[4]  = expm(im*pi/4*pX)
C1[5]  = expm(-im*pi/4*pY)
C1[6]  = im * pY
C1[7]  = expm(im*pi/4*pY)
C1[8]  = expm(-im*pi/4*pZ)
C1[9]  = im * pZ
C1[10] = expm(im*pi/4*pZ)

# Hadamard class
C1[11] = expm(-im*pi/2/sqrt(2) * (pX+pY))
C1[12] = expm(-im*pi/2/sqrt(2) * (pX-pY))
C1[13] = expm(-im*pi/2/sqrt(2) * (pX+pZ)) # standard Hadamard
C1[14] = expm(-im*pi/2/sqrt(2) * (pX-pZ))
C1[15] = expm(-im*pi/2/sqrt(2) * (pY+pZ))
C1[16] = expm(-im*pi/2/sqrt(2) * (pY-pZ))

# Axis exchange class
C1[17] = expm(-1im*pi/3/sqrt(3) * (pX+pY+pZ))
C1[18] = expm(-2im*pi/3/sqrt(3) * (pX+pY+pZ))
C1[19] = expm(-1im*pi/3/sqrt(3) * (pX-pY+pZ))
C1[20] = expm(-2im*pi/3/sqrt(3) * (pX-pY+pZ))
C1[21] = expm(-1im*pi/3/sqrt(3) * (pX+pY-pZ))
C1[22] = expm(-2im*pi/3/sqrt(3) * (pX+pY-pZ))
C1[23] = expm(-1im*pi/3/sqrt(3) * (-pX+pY+pZ))
C1[24] = expm(-2im*pi/3/sqrt(3) * (-pX+pY+pZ))

rC1 = Dict{Clifford,UInt}()

for (k,v) in enumerate(C1)
    rC1[v] = k
end

localclifford(i::Int) = C1[i]
localclifford(v::Vector) = kron(C1[v]...)

localcliffordindex(c::Clifford) = rC1[c]


text:
include("preamble.jl")

import GR  # Workaround Plots.jl world age problem
using Plots

function nullshow(plt::Plots.Plot)
    nullshow(MIME("image/png"), plt)
end


text:
using Mangal
using Test

my_tests = [
    "backbone.jl",
    "nodes.jl",
    "cache.jl",
    "attributes.jl",
    "datasets.jl",
    "networks.jl",
    "interactions.jl",
    "knownbugs.jl"
]

global test_n
global anyerrors

# Shhhh...
Mangal.login_message()

test_n = 1
anyerrors = false

for my_test in my_tests
  try
    include(my_test)
    println("[TEST $(lpad(test_n,2))] \033[1m\033[32mPASS\033[0m $(my_test)")
  catch e
    global anyerrors = true
    println("[TEST $(lpad(test_n,2))] \033[1m\033[31mFAIL\033[0m $(my_test)")
    showerror(stdout, e, backtrace())
    println()
    throw("TEST FAILED")
  end
  global test_n += 1
end

if anyerrors
  throw("Tests failed")
end


text:
#
# Date created: 2022-03-06
# Author: aradclif
#
#
############################################################################################
#### vfindmax, vfindmin

for (name, f1, fdims) ∈ zip((:vfindmax, :vfindmin), (:_vfindmax, :_vfindmin), (:_vfindmax_dims!, :_vfindmin_dims!))
    @eval function $name(A; dims=nothing)
        dims === nothing && return $f1(A)
        isone(ndims(A)) && (begin vout, iout = $f1(A) end; return [vout], [iout])
        @assert length(dims) == 1
        axes_arg = axes(A)
        axes_out = Base.setindex(axes_arg, length_one_axis(axes_arg[dims]), dims)
        Iout = similar(A, CartesianIndex{ndims(A)}, axes_out)
        Vout = similar(A, axes_out)
        Base.Cartesian.@nif 5 d -> (d ≤ ndims(A) && dims == d) d -> begin
            Rpre = CartesianIndices(ntuple(i -> axes_arg[i], d - 1))
            Rpost = CartesianIndices(ntuple(i -> axes_arg[i + d], ndims(A) - d))
        end d -> begin
            Rpre = CartesianIndices(axes_arg[1:dims-1])
            Rpost = CartesianIndices(axes_arg[dims+1:end])
        end
        $fdims(Vout, Iout, Rpre, 1:size(A, dims), Rpost, A)
    end
    @eval $name(A, ::Colon) = $name(A)
end
for (op, name, init) ∈ zip((:>, :<), (:_vfindmax_dims!, :_vfindmin_dims!), (:typemin, :typemax))
    @eval function $name(Vout, Iout, Rpre, is, Rpost, A)
        s = $init(eltype(A))
        @inbounds for Ipost ∈ Rpost, Ipre ∈ Rpre
            indm = 0
            mval = s
            @turbo for i ∈ is
                newm = $op(A[Ipre, i, Ipost], mval)
                mval = newm ? A[Ipre, i, Ipost] : mval
                indm = newm ? i : indm
            end
            Vout[Ipre, 1, Ipost] = mval
            Iout[Ipre, 1, Ipost] = CartesianIndex(Ipre, indm, Ipost)
        end
        return Vout, Iout
    end
end
for (op, name, init) ∈ zip((:>, :<), (:_vfindmax, :_vfindmin), (:typemin, :typemax))
    @eval function $name(A::AbstractArray{T, N}) where {T, N}
        indm = 0
        mval = $init(T)
        @turbo for i ∈ eachindex(A)
            newm = $op(A[i], mval)
            mval = newm ? A[i] : mval
            indm = newm ? i : indm
        end
        mval, CartesianIndices(A)[indm]
    end
    @eval function $name(A::AbstractVector{T}) where {T}
        indm = 0
        mval = $init(T)
        @turbo for i ∈ eachindex(A)
            newm = $op(A[i], mval)
            mval = newm ? A[i] : mval
            indm = newm ? i : indm
        end
        mval, indm
    end
end

#### vtfindmax, vtfindmin : threaded version
for (name, f1, fdims) ∈ zip((:vtfindmax, :vtfindmin), (:_vtfindmax, :_vtfindmin), (:_vtfindmax_dims!, :_vtfindmin_dims!))
    @eval function $name(A; dims=nothing)
        dims === nothing && return $f1(A)
        isone(ndims(A)) && (begin vout, iout = $f1(A) end; return [vout], [iout])
        @assert length(dims) == 1
        axes_arg = axes(A)
        axes_out = Base.setindex(axes_arg, length_one_axis(axes_arg[dims]), dims)
        Iout = similar(A, CartesianIndex{ndims(A)}, axes_out)
        Vout = similar(A, axes_out)
        Base.Cartesian.@nif 5 d -> (d ≤ ndims(A) && dims == d) d -> begin
            Rpre = CartesianIndices(ntuple(i -> axes_arg[i], d - 1))
            Rpost = CartesianIndices(ntuple(i -> axes_arg[i + d], ndims(A) - d))
        end d -> begin
            Rpre = CartesianIndices(axes_arg[1:dims-1])
            Rpost = CartesianIndices(axes_arg[dims+1:end])
        end
        $fdims(Vout, Iout, Rpre, 1:size(A, dims), Rpost, A)
    end
    @eval $name(A, ::Colon) = $name(A)
end
for (op, name, init) ∈ zip((:>, :<), (:_vfindmax_prepost!, :_vfindmin_prepost!), (:typemin, :typemax))
    @eval function $name(Vout, Iout, Ipre, is, Ipost, A)
        indm = 0
        mval = $init(eltype(A))
        @turbo for i ∈ is
            newm = $op(A[Ipre, i, Ipost], mval)
            mval = newm ? A[Ipre, i, Ipost] : mval
            indm = newm ? i : indm
        end
        Vout[Ipre, 1, Ipost] = mval
        Iout[Ipre, 1, Ipost] = CartesianIndex(Ipre, indm, Ipost)
        return Vout, Iout
    end
end
for (f, name) ∈ zip((:_vfindmax_prepost!, :_vfindmin_prepost!), (:_vtfindmax_dims!, :_vtfindmin_dims!))
    @eval function $name(Vout, Iout, Rpre, is, Rpost, A)
        if length(Rpre) > length(Rpost)
            @batch for Ipre ∈ Rpre, Ipost ∈ Rpost
                $f(Vout, Iout, Ipre, is, Ipost, A)
            end
        else
            @batch for Ipost ∈ Rpost, Ipre ∈ Rpre
                $f(Vout, Iout, Ipre, is, Ipost, A)
            end
        end
        return Vout, Iout
    end
end
for (op, name, init) ∈ zip((:>, :<), (:_vtfindmax, :_vtfindmin), (:typemin, :typemax))
    @eval function $name(A::AbstractArray{T, N}) where {T, N}
        indm = 0
        mval = $init(T)
        @tturbo for i ∈ eachindex(A)
            newm = $op(A[i], mval)
            mval = newm ? A[i] : mval
            indm = newm ? i : indm
        end
        mval, CartesianIndices(A)[indm]
    end
    @eval function $name(A::AbstractVector{T}) where {T}
        indm = 0
        mval = $init(T)
        @tturbo for i ∈ eachindex(A)
            newm = $op(A[i], mval)
            mval = newm ? A[i] : mval
            indm = newm ? i : indm
        end
        mval, indm
    end
end


text:
@testset "probabilityweightedmoment_gp.jl" begin
    n = 10000
    θ = [0.0 ; .2]

    pd = GeneralizedPareto(exp(θ[1]), θ[2])
    y = rand(pd, n)

    @testset "gpfitpwm(y)" begin
        # stationary model building
        fm = Extremes.gpfitpwm(y)

        cinterval = cint(fm)

        @test [x[1] for x in cinterval] <= θ
        @test θ <= [x[2] for x in cinterval]

    end

    @testset "gpfitpwm(df, datacol)" begin
        # stationary model building
        df = DataFrame(y = y)
        fm = Extremes.gpfitpwm(df, :y)

        cinterval = cint(fm)

        @test [x[1] for x in cinterval] <= θ
        @test θ <= [x[2] for x in cinterval]

    end

    @testset "gpfitpwm(model)" begin
        # non-stationary warn
        model = Extremes.ThresholdExceedance(Variable("y", y), logscalecov = [Variable("t", collect(1:n))])

        @test_logs (:warn, "covariates cannot be included in the model when estimating the
            paramters by the probability weighted moment parameter estimation.
            The estimates for the stationary model is returned.") Extremes.gpfitpwm(model)

        # stationary GP fit by pwm
        model = Extremes.ThresholdExceedance(Variable("y", y))

        fm = Extremes.gpfitpwm(model)

        cinterval = cint(fm)

        @test [x[1] for x in cinterval] <= θ
        @test θ <= [x[2] for x in cinterval]

    end
end


text:
using Symbolics: scalarize
@variables u[1:3]
sisequal(b1, b2) = all(isequal.(scalarize(b1), b2 isa Array ? scalarize.(b2) : scalarize(b2)))
@test sisequal(sin_basis(u, 1), sin.(1 .* u))
@test sisequal(cos_basis(u, 2), vcat(cos.(1 .* u), cos.(2 .*u)))
@test sisequal(chebyshev_basis(u, 1), cos.( 1 .* acos.(u)))
@test sisequal(fourier_basis(u, 1), sin.(1 .* u ./2))
@test sisequal(monomial_basis(u, 1), [1; u.^1])
@test sisequal(polynomial_basis(u, 2), [1; u[1]^1; u[1]^2; u[2]^1; u[1]^1*u[2]^1; u[2]^2; u[3]^1; u[1]^1*u[3]^1; u[2]^1*u[3]^1; u[3]^2])

@test sisequal(sin_basis(u, 1:2), vcat([sin.(i .* u) for i in 1:2]...))
@test sisequal(cos_basis(u, 1:5), vcat([cos.(i .* u) for i in 1:5]...))
@test sisequal(chebyshev_basis(u, [1;2]), vcat([cos.(i .* acos.(u)) for i in 1:2]...))
@test sisequal(fourier_basis(u, 1:2), vcat(sin.(1 .* u ./2), cos.( 2 .* u ./ 2)))