text:
# Electrical circuit gate sizing example adapted from Boyd, Kim, Vandenberghe, and Hassibi, "A Tutorial on Geometric Programming", p28
# CVX code adapted from simple_dig_ckt_sizing.m by Almir Mutapcic
# http://cvxr.com/cvx/examples/gp_tutorial/html/simple_dig_ckt_sizing.html
#
# Choose gate scale factors y_i to minimize ckt delay, subject to limits on the total area and power
#   minimize     D
#   subject to   P <= Pmax,
#                A <= Amax
#                y >= 1
# Where variables are scale factors y, which are constrained to equal discrete values
# Data is hard-coded from the circuit shown in figure 4 of "A Tutorial on Geometric Programming"

# source: https://github.com/JuliaOpt/Pajarito.jl/blob/master/examples/gatesizing.jl

using Convex, Pajarito

# Set up Convex.jl model, solve, print solution
function gatesizing(yUB)
    fe = [1, 0.8, 1, 0.7, 0.7, 0.5, 0.5] .* [1, 2, 1, 1.5, 1.5, 1, 2]
    Cout6 = 10
    Cout7 = 10

    y = Variable(7, Positive())
    z = Variable(yUB, 7, Positive(), :Bin)

    D1 = exp(-y[1]) + exp(-y[1] + y[4])
    D2 = 2 * exp(-y[2]) + exp(-y[2] + y[4]) + exp(-y[2] + y[5])
    D3 = 2 * exp(-y[3]) + exp(-y[3] + y[5]) + exp(-y[3] + y[7])
    D4 = 2 * exp(-y[4]) + exp(-y[4] + y[6]) + exp(-y[4] + y[7])
    D5 = exp(-y[5]) + exp(-y[5] + y[7])
    D6 = Cout6 * exp(-y[6])
    D7 = Cout7 * exp(-y[7])

    P = minimize(
        maximum([(D1+D4+D6), (D1+D4+D7), (D2+D4+D6), (D2+D4+D7), (D2+D5+D7), (D3+D5+D6), (D3+D7)]),
        sum(fe .* exp(y)) <= 20,
        sum(exp(y)) <= 100)

    for i in 1:7
        P.constraints += (sum(z[:,i]) == 1)
        P.constraints += (y[i] == sum([log(j) * z[j,i] for j=1:yUB]))
    end

    solve!(P, solver)

    println("\nCircuit delay (obj) = $(P.optval)")
    println("Scale factors (exp(y)):\n$(exp.(y.value))")
    println("Value indicators (z):\n$(round.(z.value))")
end

#=========================================================
Choose solvers and options
=========================================================#

mip_solver_drives = false
log_level = 3
rel_gap = 1e-5

# using Cbc
# mip_solver = CbcSolver()

using CPLEX
mip_solver = CplexSolver(
    CPX_PARAM_SCRIND=(mip_solver_drives ? 1 : 0),
    # CPX_PARAM_SCRIND=1,
    CPX_PARAM_EPINT=1e-8,
    CPX_PARAM_EPRHS=1e-7,
    CPX_PARAM_EPGAP=(mip_solver_drives ? 1e-5 : 1e-9)
)

using SCS
cont_solver = SCSSolver(eps=1e-5, max_iters=1000000, verbose=1)

# using ECOS
# cont_solver = ECOSSolver(verbose=false)

solver = PajaritoSolver(
    mip_solver_drives=mip_solver_drives,
    log_level=log_level,
    rel_gap=rel_gap,
	mip_solver=mip_solver,
	cont_solver=cont_solver,
    init_exp=true,
)

#=========================================================
Specify data
=========================================================#

# Integer upper bound on integer gate size factors
# Optimal solution for this instances has all exp(y) <= 3,
# so yUB larger than 2 is not constraining, but impacts number of variables
yUB = 3

#=========================================================
Solve Convex.jl model
=========================================================#

# Optimal solution is [2,3,3,3,2,3,3] with value 8.3333
gatesizing(yUB)


text:
using Printf
using NPZ
include("util.jl")
include("kmeans_filters.jl")
include("quantum_filters.jl")

log_file = open("run_filters.log", "a")

# for name in ARGS
    # target_label = parse(Int, split(name, "-")[3][end:end])
name = "test"
reps = npzread("/content/gdrive/My Drive/Resnet/Cheat/ones.np")'
n = size(reps)[2]
polluted_mark = npzread("/content/gdrive/My Drive/Resnet/Cheat/mark.np")
eps = sum(polluted_mark)
# removed = round(Int, 1.5*eps)
removed = round(Int, 1.5*3424)

    # @printf("%s: Running PCA filter\n", name)
    # reps_pca, U = pca(reps, 1)
    # pca_poison_ind = k_lowest_ind(-abs.(mean(reps_pca[1, :]) .- reps_pca[1, :]), round(Int, 1.5*eps))
    # poison_removed = sum(pca_poison_ind[end-eps+1:end])
    # clean_removed = removed - poison_removed
    # @show poison_removed, clean_removed
    # @printf(log_file, "%s-pca: %d, %d\n", name, poison_removed, clean_removed)
    # npzwrite("output/$(name)/mask-pca-target.npy", pca_poison_ind)

    # @printf("%s: Running kmeans filter\n", name)
    # kmeans_poison_ind = .! kmeans_filter2(reps, eps)
    # poison_removed = sum(kmeans_poison_ind[end-eps+1:end])
    # clean_removed = removed - poison_removed
    # @show poison_removed, clean_removed
    # @printf(log_file, "%s-kmeans: %d, %d\n", name, poison_removed, clean_removed)
    # npzwrite("output/$(name)/mask-kmeans-target.npy", kmeans_poison_ind)

@printf("Running quantum filter\n")
quantum_poison_ind = .! rcov_auto_quantum_filter(polluted_mark, reps, eps)
# quantum_poison_ind = .! rcov_quantum_filter(reps, eps, 3)
pr = 0
for i in 1:n
    if polluted_mark[i] && quantum_poison_ind[i]
        global pr += 1
    end
    quantum_poison_ind[i] = !quantum_poison_ind[i]
end

# npzwrite("/content/gdrive/My Drive/Resnet/Cheat/Retrain/select(v2).np", quantum_poison_ind)
clean_removed = removed - pr
@show pr, clean_removed, eps
# end

ext: jl

max_stars_repo_name: Lteins/qfilter

max_issues_repo_name: Lteins/qfilter

max_forks_repo_name: Lteins/qfilter


text:
if !m._obbt_performed_flag
    if m._nonlinear_evaluator_created
        set_node!(m._working_problem._relaxed_evaluator, m._current_node)
        set_node_flag!(m)
        set_reference_point!(m)
        fill!(m._working_problem._relaxed_evaluator.subexpressions_eval, false)
    end
    update_relaxed_problem_box!(m)
    relax_constraints!(m)
end
relax_objective!(m)

# Optimizes the object
MOI.set(m.relaxed_optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)

"""

"""
function intrepret_relaxed_solution!(m::Optimizer, d::T) where T
    valid_flag, feasible_flag = is_globally_optimal(m._lower_termination_status,
                                                    m._lower_result_status)
    if valid_flag && feasible_flag
        set_dual!(m)
        m._cut_add_flag = true
        m._lower_feasibility = true
        m._lower_objective_value = MOI.get(d, MOI.ObjectiveValue())
        for i = 1:_variable_num(_working_problem(m))
             m._lower_solution[i] = MOI.get(d, MOI.VariablePrimal(),
                                               m._relaxed_variable_index[i])
        end
    elseif valid_flag
        m._cut_add_flag = false
        m._lower_feasibility  = false
        m._lower_objective_value = -Inf
    end
    return valid_flag
end

function _modify_existing_relax(m)
    flag = MOIU.supports_default_copy_to(m)
    flag &=
    return flag
end

function MOIU.supports_default_copy_to(opt::T)
end

function lower_problem!(t::ExtensionType, m::Optimizer)
    if _modify_existing_relax(m)
        _modify_relax!(m)
    elseif _add_delete_existing_relax(m)
        _add_delete_relax!(m)
    else
        MOI.empty!(opt)
        if MOIU.supports_default_copy_to(m.relaxed_optimizer)
            _add_relax!(m, m.relaxed_optimizer)
        else
            MOI.empty!(m.cached_opt)
            _add_relax!(m, cached_opt)
            MOIU.copyto!(m.relaxed_optimizer, m.cached_opt)
        end
    end

    MOI.optimize!(m.relaxed_optimizer)
    if !intrepret_relaxed_solution(m, m.relaxed_optimizer)
        fallback_interval_lower_bound!(m, n)
    end
    return
end

if supports_constraint_mod(opt)
    # modify
else
    MOI.empty!(opt)
    if MOIU.supports_default_copy_to(opt::T)
        _relax!(opt)
    else
        _relax!(cached_opt)
        MOIU.copyto!(opt, cached_opt)
    end
end


text:
struct PeriodicArray{T,N} <: AbstractArray{T,N}
    data::Array{T,N}
end
PeriodicArray{T}(initializer, args...) where T =
    PeriodicArray(Array{T}(initializer, args...))
PeriodicArray{T,N}(initializer, args...) where {T,N} =
    PeriodicArray(Array{T,N}(undef, args...))

Base.size(a::PeriodicArray) = size(a.data)
Base.size(a::PeriodicArray,i) = size(a.data,i)
Base.length(a::PeriodicArray) = length(a.data)
Base.getindex(a::PeriodicArray{T,N}, I::Vararg{Int,N}) where {T,N} =
    @inbounds getindex(a.data, map(mod1, I, size(a.data))...)
Base.setindex!(a::PeriodicArray{T,N}, v, I::Vararg{Int,N}) where {T,N} =
    @inbounds setindex!(a.data, v, map(mod1, I, size(a.data))...)

Base.similar(a::PeriodicArray, dims::Union{Integer, AbstractUnitRange}...) =
    PeriodicArray(similar(a.data, dims...))
Base.similar(a::PeriodicArray, T::Type, dims::Union{Integer, AbstractUnitRange}...) =
    PeriodicArray(similar(a.data, T, dims...))

Base.copy(a::PeriodicArray) = PeriodicArray(copy(a.data))
function Base.copyto!(dst::PeriodicArray,src::PeriodicArray)
    copyto!(dst.data,src.data);
    dst
end
# not necessary but maybe more efficient

Base.convert(::Type{PeriodicArray{T}}, a::PeriodicArray) where {T} =
    PeriodicArray(convert(Array{T}, a.data))
Base.convert(::Type{PeriodicArray{T,N}}, a::PeriodicArray) where {T,N} =
    PeriodicArray(convert(Array{T,N}, a.data))

Base.checkbounds(a::PeriodicArray, I...) = true

Base.circshift(t::PeriodicArray{T,N},tup::Tuple{Vararg{Integer,N}}) where{T,N}= PeriodicArray{T,N}(circshift(t.data,tup))
Base.repeat(t::PeriodicArray,args::Vararg{Integer,N} where N) = PeriodicArray(repeat(t.data,args...))


text:
# matrix of constats

export joConstants
"""
    julia> op = joConstants(m,a;[DDT=...,][RDT=...,][name=...])
    julia> op = joConstants(m,n,a;[DDT=...,][RDT=...,][name=...])

Operator equivalent to matrix of same elements

# Signature

    joConstants(m::Integer,a::EDT;
        DDT::DataType=EDT,RDT::DataType=promote_type(EDT,DDT),
        name::String="joConstants") where {EDT<:Number}
    joConstants(m::Integer,n::Integer,a::EDT;
        DDT::DataType=EDT,RDT::DataType=promote_type(EDT,DDT),
        name::String="joConstants") where {EDT<:Number}

# Arguments

- `m`: size
- `a`: element value, default DDT/RDT will have the same type as `a`
- optional
    - `n`: 2nd dimension if not square
- keywords
    - `DDT`: domain data type
    - `RDT`: range data type
    - `name`: custom name

# Examples

    joConstants(m,a)
    joConstants(m,n,a)

examples with DDT/RDT

    joConstants(m,a; DDT=Float32)
    joConstants(m,a; DDT=Float32,RDT=Float64)

"""
joConstants(m::Integer,a::EDT;
    DDT::DataType=EDT,RDT::DataType=promote_type(EDT,DDT),
    name::String="joConstants") where {EDT<:Number} = joConstants(m,m,a;DDT=DDT,RDT=RDT,name=name)
joConstants(m::Integer,n::Integer,a::EDT;
    DDT::DataType=EDT,RDT::DataType=promote_type(EDT,DDT),
    name::String="joConstants") where {EDT<:Number} =
    joMatrix{DDT,RDT}(name,m,n,
        v1->jo_convert(RDT,a*ones(eltype(v1),m,1)*sum(v1,dims=1),false),
        v2->jo_convert(DDT,a*ones(eltype(v2),n,1)*sum(v2,dims=1),false),
        v3->jo_convert(DDT,conj(a)*ones(eltype(v3),n,1)*sum(v3,dims=1),false),
        v4->jo_convert(RDT,conj(a)*ones(eltype(v4),m,1)*sum(v4,dims=1),false),
        @joNF, @joNF, @joNF, @joNF
        )


text:
# Small example showing how to instrument Julia code at IR level
# in order to collect basic block coverage information during execution.

# We want to measure basic block coverage when testing this method:
function f(x, y)
    # Just insert a delay so we can measure overhead more realistically below...
    sleep(0.0001)
    if x < y
        return x + y
    else
        return x - y
    end
end

# We will use a coverage object to count the number of executions of each basic block.
abstract type CodeTracer end
mutable struct BasicBlockCoverage <: CodeTracer
    blockexecutions::Vector{Int}
end

BasicBlockCoverage(n::Int) = BasicBlockCoverage(zeros(Int, n))

function enterfunction!(c::BasicBlockCoverage, numblocks::Int)
    @assert numblocks == length(c.blockexecutions)
end

enterblock!(c::BasicBlockCoverage, block::Int) =
    c.blockexecutions[block] += 1

function coverage(c::BasicBlockCoverage)
    length(c.blockexecutions) == 0 && return 0.0
    sum(c.blockexecutions .> 0) / length(c.blockexecutions)
end

# Now use IRTools so we can instrument functions to insert callbacks
# to a tracer object that we will take in as a new first argument.
using IRTools
using IRTools: @code_ir, blocks, argument!, func, xcall

numblocks(ir) = length(blocks(ir))

function instrument!(ir)
    bs = blocks(ir)
    tracerarg = argument!(ir; at = 2) # add tracer argument to function
    for bi in eachindex(bs)
        pushfirst!(bs[bi], xcall(Main, :enterblock!, tracerarg, bi))
    end
    # This should come after the loop above so it is really added first in the whole
    # function:
    pushfirst!(ir, xcall(Main, :enterfunction!, tracerarg, numblocks(ir)))
    ir
end

# Get the IR, instrument it, and turn into a new lambda/function:
old_ir = @code_ir f(2, 3)
instrumented_ir = instrument!(old_ir)
f2 = func(instrumented_ir); # Make a new function from the instrumented ir

# Ok, let's test this!

# Ensure same output:
origres = f(2, 3)
c0 = BasicBlockCoverage(numblocks(old_ir))
res = f2(nothing, c0, 2, 3)
@assert origres == res

# If we only run it once we should cover only 2 of 3 basic blocks:
c = BasicBlockCoverage(numblocks(old_ir))
f2(nothing, c, 2, 3)
@assert round(coverage(c), digits = 2) == 0.67

# But if we run 99 more times with random inputs we are very likely
# to get 100% coverage:
for _ in 1:99
    x = rand(0:10)
    y = rand(0:10)
    f2(nothing, c, x, y)
end
@assert coverage(c) == 1.0

# Overhead is not too bad:
using BenchmarkTools
torig = @belapsed f(2, 3)
t     = @belapsed f2(nothing, c, 2, 3)
overhead = round(100.0 * (t/torig - 1.0), digits = 2) # between 3-17% on my machine...

# This will not work for more complex methods with optional args, keyword args etc
# but it is a good start:
macro instrument(expr)
    fname = first(expr.args)
    args = Symbol[gensym() for _ in 1:(length(expr.args)-1)]
    ir, newf = gensym(), gensym()
    quote
        oldir = IRTools.@code_ir $expr
        newfunc = IRTools.func(instrument!(oldir))
        $(esc(fname))(c::CodeTracer, $(args...)) = newfunc(nothing, c, $(args...))
        newfunc # Return the func if they want to do other stuff with it...
    end
end

#@macroexpand @instrument f(2, 3)
@instrument f(2, 3)
c = BasicBlockCoverage(3)
f(c, 2, 3)


text:
# Introduction to DataFrames
# https://www.juliabox.com/notebook/notebooks/tutorials/intro-to-julia-DataFrames/01_constructors.ipynb
# Let's get started by loading the DataFrames package

using InteractiveUtils
import DataFrames: DataFrame
import Random: randstring

versioninfo()

#= Constructors
In this section, you'll see many ways to create a DataFrame using
the DataFrame() constructor.
=#
df = DataFrame(A=1:3, B=rand(3), C=randstring.([3,3,3]))
@show df
#=
df = 3×3 DataFrame
│ Row │ A │ B        │ C   │
├─────┼───┼──────────┼─────┤
│ 1   │ 1 │ 0.643894 │ wpd │
│ 2   │ 2 │ 0.611762 │ 4He │
│ 3   │ 3 │ 0.665899 │ u1Y │
=#

# We can create a DataFrame from a dictionary, in which case

dict = Dict("C"=> [1,2], "B"=>[true,false], "A"=>['a','b'])
df = DataFrame(dict)
@show df

#=
Rather than explicitly creating a dictionary first, as above,
we coult pass DataFrame arguments with the syntax of dictionary
key-value pairs.
=#

df = DataFrame(:A => [1,2], :B => [true,false], :C => ['a','b'])
@show df

#create a DataFrame from a matrix
@show DataFrame(transpose([1 2 3]))
df = DataFrame([1 2 3; 4 5 6; 7 8 9])
@show df
m = convert(Matrix, df)
@show m
# pass a second argument to give the columns names
@show DataFrame([1:3,4:6,7:9], [:A,:B,:C])

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# test index notation using @tensor macro
#-----------------------------------------
withblas = TensorContractionsXD.use_blas() ? "with" : "without"
withcache = TensorContractionsXD.use_cache() ? "with" : "without"
@testset "Index Notation $withblas BLAS and $withcache cache" begin
    t0 = time()
    A = randn(Float64, (3,5,4,6))
    p = (4,1,3,2)
    C1 = permutedims(A, p)
    @tensor C2[4,1,3,2] := A[1,2,3,4]
    @test C1 ≈ C2
    @test_throws TensorContractionsXD.IndexError begin
        @tensor C[1,2,3,4] := A[1,2,3]
    end
    @test_throws TensorContractionsXD.IndexError begin
        @tensor C[1,2,3,4] := A[1,2,2,4]
    end
    println("tensorcopy: $(time()-t0) seconds")
    t0 = time()

    B=randn(Float64, (5,6,3,4))
    p=[3,1,4,2]
    @tensor C1[3,1,4,2] := A[3,1,4,2] + B[1,2,3,4]
    C2=A+permutedims(B, p)
    @test C1 ≈ C2
    @test_throws DimensionMismatch begin
        @tensor C[1,2,3,4] := A[1,2,3,4] + B[1,2,3,4]
    end
    println("tensoradd: $(time()-t0) seconds")
    t0 = time()

    A=randn(Float64, (50,100,100))
    @tensor C1[a] := A[a, b', b']
    C2=zeros(50)
    for i=1:50
        for j=1:100
            C2[i]+=A[i, j, j]
        end
    end
    @test C1 ≈ C2
    A=randn(Float64, (3,20,5,3,20,4,5))
    @tensor C1[e, a, d] := A[a, b, c, d, b, e, c]
    C2=zeros(4,3,3)
    for i1=1:4, i2=1:3, i3=1:3
        for j1=1:20, j2=1:5
            C2[i1, i2, i3]+=A[i2, j1, j2, i3, j1, i1, j2]
        end
    end
    @test C1 ≈ C2
    println("tensortrace: $(time()-t0) seconds")
    t0 = time()

    A=randn(Float64, (3,20,5,3,4))
    B=randn(Float64, (5,6,20,3))
    @tensor C1[a, g, e, d, f] := A[a, b, c, d, e] * B[c, f, b, g]
    C2=zeros(3,3,4,3,6)
    for a=1:3, b=1:20, c=1:5, d=1:3, e=1:4, f=1:6, g=1:3
        C2[a, g, e, d, f] += A[a, b, c, d, e] * B[c, f, b, g]
    end
    @test C1 ≈ C2
    @test_throws TensorContractionsXD.IndexError begin
        @tensor A[a, b, c, d] * B[c, f, b, g]
    end
    println("tensorcontract 1: $(time()-t0) seconds")
    t0 = time()

    A=randn(Float64, (5,5,5,5))
    B=rand(ComplexF64, (5,5,5,5))
    @tensor C1[1,2,5,6,3,4,7,8] := A[1,2,3,4] * B[5,6,7,8]
    C2=reshape(kron(reshape(B, (25,25)), reshape(A, (25,25))), (5,5,5,5,5,5,5,5))
    @test C1 ≈ C2
    @test_throws TensorContractionsXD.IndexError begin
        @tensor C[a, b, c, d, e, f, g, i] := A[a, b, c, d] * B[e, f, g, h]
    end
    println("tensorcontract 2: $(time()-t0) seconds")
    t0 = time()

    Da=10
    Db=15
    Dc=4
    Dd=8
    De=6
    Df=7
    Dg=3
    Dh=2
    A=rand(ComplexF64, (Da, Dc, Df, Da, De, Db, Db, Dg))
    B=rand(ComplexF64, (Dc, Dh, Dg, De, Dd))
    C=rand(ComplexF64, (Dd, Dh, Df))
    @tensor D1[d, f, h] := A[a, c, f, a, e, b, b, g] * B[c, h, g, e, d]+0.5 * C[d, h, f]
    D2=zeros(ComplexF64, (Dd, Df, Dh))
    for d=1:Dd, f=1:Df, h=1:Dh
        D2[d, f, h] += 0.5 * C[d, h, f]
        for a=1:Da, b=1:Db, c=1:Dc, e=1:De, g=1:Dg
            D2[d, f, h] += A[a, c, f, a, e, b, b, g] * B[c, h, g, e, d]
        end
    end
    @test D1 ≈ D2
    @test norm(vec(D1)) ≈ sqrt(abs((@tensor scalar(D1[d, f, h] * conj(D1[d, f, h])))))
    println("tensorcontract 3: $(time()-t0) seconds")
    t0 = time()

    Abig=randn(Float64, (30,30,30,30))
    A=view(Abig,1 .+ 3 .* (0:9),2 .+ 2 .* (0:6),5 .+ 4 .* (0:6),4 .+ 3 .* (0:8))
    p=[3,1,4,2]
    Cbig=zeros(ComplexF64, (50,50,50,50))
    C=view(Cbig,13 .+ (0:6),11 .+ 4 .* (0:9),15 .+ 4 .* (0:8),4 .+ 3 .* (0:6))
    Acopy = copy(A)
    Ccopy = copy(C)
    @tensor C[3,1,4,2] = A[1,2,3,4]
    @tensor Ccopy[3,1,4,2] = Acopy[1,2,3,4]
    @test C ≈ Ccopy
    @test_throws TensorContractionsXD.IndexError begin
        @tensor C[3,1,4,2] = A[1,2,3]
    end
    @test_throws DimensionMismatch begin
        @tensor C[3,1,4,2] = A[3,1,4,2]
    end
    @test_throws TensorContractionsXD.IndexError begin
        @tensor C[1,1,2,3] = A[1,2,3,4]
    end
    println("views: $(time()-t0) seconds")
    t0 = time()

    Cbig=zeros(ComplexF64, (50,50,50,50))
    C=view(Cbig, 13 .+ (0:6), 11 .+ 4 .* (0:9), 15 .+ 4 .* (0:8), 4 .+ 3 .* (0:6))
    Acopy=permutedims(copy(A), p)
    Ccopy=copy(C)
    α=randn(Float64)
    β=randn(Float64)
    @tensor C[3,1,4,2] = β * C[3,1,4,2] + α * A[1,2,3,4]
    Ccopy=β * Ccopy+α * Acopy
    @test C ≈ Ccopy
    @test_throws TensorContractionsXD.IndexError begin
        @tensor C[3,1,4,2] = 0.5 * C[3,1,4,2] + 1.2 * A[1,2,3]
    end
    @test_throws DimensionMismatch  begin
        @tensor C[3,1,4,2] = 0.5 * C[3,1,4,2] + 1.2 * A[3,1,2,4]
    end
    @test_throws TensorContractionsXD.IndexError  begin
        @tensor C[1,1,2,3] = 0.5 * C[1,1,2,3] + 1.2 * A[1,2,3,4]
    end
    println("more views: $(time()-t0) seconds")
    t0 = time()

    Abig=rand(Float64, (30,30,30,30))
    A=view(Abig,1 .+ 3 .* (0:8),2 .+ 2 .* (0:14),5 .+ 4 .* (0:6),7 .+ 2 .* (0:8))
    Bbig=rand(ComplexF64, (50,50))
    B=view(Bbig,13 .+ (0:14),3 .+ 5 .* (0:6))
    Acopy=copy(A)
    Bcopy=copy(B)
    α=randn(Float64)
    @tensor B[b, c] += α * A[a, b, c, a]
    for i=1 .+ (0:8)
        Bcopy += α * view(A, i, :, :, i)
    end
    @test B ≈ Bcopy
    @test_throws TensorContractionsXD.IndexError begin
        @tensor B[b, c] += α * A[a, b, c]
    end
    @test_throws DimensionMismatch begin
        @tensor B[c, b] += α * A[a, b, c, a]
    end
    @test_throws TensorContractionsXD.IndexError begin
        @tensor B[c, b] += α * A[a, b, a, a]
    end
    @test_throws DimensionMismatch begin
        @tensor B[c, b] += α * A[a, b, a, c]
    end
    println("even more views: $(time()-t0) seconds")
    t0 = time()

    Abig=rand(Float64, (30,30,30,30))
    A=view(Abig, 1 .+ 3 .* (0:8), 2 .+ 2 .* (0:14), 5 .+ 4 .* (0:6), 7 .+ 2 .* (0:8))
    Bbig=rand(ComplexF64, (50,50,50))
    B=view(Bbig, 3 .+ 5 .* (0:6), 7 .+ 2 .* (0:7), 13 .+ (0:14))
    Cbig=rand(ComplexF32, (40,40,40))
    C=view(Cbig,3 .+ 2 .* (0:8),13 .+ (0:8),7 .+ 3 .* (0:7))
    Acopy=copy(A)
    Bcopy=copy(B)
    Ccopy=copy(C)
    α=randn(Float64)
    for d=1 .+ (0:8), a=1 .+ (0:8), e=1 .+ (0:7)
        for b=1 .+ (0:14), c=1 .+ (0:6)
            Ccopy[d, a, e] -=α * A[a, b, c, d] * conj(B[c, e, b])
        end
    end
    @tensor C[d, a, e] -= α * A[a, b, c, d] * conj(B[c, e, b])
    @test C ≈ Ccopy
    println("and some more views: $(time()-t0) seconds")
    t0 = time()

    Cbig=rand(ComplexF32, (40,40,40))
    C=view(Cbig,3 .+ 2 .* (0:8),13 .+ (0:8),7 .+ 3 .* (0:7))
    Ccopy=copy(C)
    for d=1 .+ (0:8), a=1 .+ (0:8), e=1 .+ (0:7)
        for b=1 .+ (0:14), c=1 .+ (0:6)
            Ccopy[d, a, e] += α * A[a, b, c, d] * conj(B[c, e, b])
        end
    end
    @tensor C[d, a, e] += α * A[a, b, c, d] * conj(B[c, e, b])
    @test C ≈ Ccopy
    @test_throws TensorContractionsXD.IndexError begin
        @tensor C[d, a, e] += α * A[a, b, c, a] * B[c, e, b]
    end
    @test_throws TensorContractionsXD.IndexError begin
        @tensor C[d, a, e] += α * A[a, b, c, d] * B[c, b]
    end
    @test_throws TensorContractionsXD.IndexError begin
        @tensor C[d, e] += α * A[a, b, c, d] * B[c, e, b]
    end
    @test_throws DimensionMismatch begin
        @tensor C[d, e, a] += α * A[a, b, c, d] * B[c, e, b]
    end
    println("Float32 views: $(time()-t0) seconds")
    t0 = time()

    # Simple function example
    @tensor function f(A, b)
        w[x] := (1//2)*A[x,y]*b[y]
        return w
    end
    for T in (Float32, Float64, ComplexF32, ComplexF64, BigFloat)
        A = rand(T, 10, 10)
        b = rand(T, 10)
        @test f(A,b) ≈ (1//2)*A*b
    end

    # Example from README.md
    α=randn()
    A=randn(5,5,5,5,5,5)
    B=randn(5,5,5)
    C=randn(5,5,5)
    D=zeros(5,5,5)
    @tensor begin
        D[a, b, c] = A[a, e, f, c, f, g] * B[g, b, e] + α * C[c, a, b]
        E[a, b, c] := A[a, e, f, c, f, g] * B[g, b, e] + α * C[c, a, b]
    end
    @test D == E
    println("readme example: $(time()-t0) seconds")
    t0 = time()

    # Some tensor network examples
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64, BigFloat)
        D1, D2, D3 = 30, 40, 20
        d1, d2 = 2, 3
        A1 = rand(T, D1, d1, D2) .- 1//2
        A2 = rand(T, D2, d2, D3) .- 1//2
        rhoL = rand(T, D1, D1) .- 1//2
        rhoR = rand(T, D3, D3) .- 1//2
        H = rand(T, d1, d2, d1, d2) .- 1//2
        A12 = reshape(reshape(A1, D1 * d1, D2) * reshape(A2, D2, d2 * D3), (D1, d1, d2, D3))
        rA12 = reshape(reshape(rhoL * reshape(A12, (D1, d1*d2*D3)), (D1*d1*d2, D3)) * rhoR, (D1, d1, d2, D3))
        HrA12 = permutedims(reshape(reshape(H, (d1 * d2, d1*d2)) * reshape(permutedims(rA12, (2,3,1,4)), (d1 * d2, D1 * D3)), (d1, d2, D1, D3)), (3,1,2,4))
        E = dot(A12, HrA12)
        @tensor HrA12′[a, s1, s2, c] := rhoL[a, a'] * A1[a', t1, b] * A2[b, t2, c'] * rhoR[c', c] * H[s1, s2, t1, t2]
        @tensor HrA12′′[:] := rhoL[-1, 1] * H[-2, -3, 4, 5] * A2[2, 5, 3] * rhoR[3, -4] * A1[1, 4, 2] # should be contracted in exactly same order
        @tensor HrA12′′′[a, s1, s2, c] := rhoL[a, a'] * H[s1, s2, t1, t2] * A2[b, t2, c'] * rhoR[c', c] * A1[a', t1, b] order=(a',b,c',t1,t2)# should be contracted in exactly same order
        @tensoropt HrA12′′′′[:] := rhoL[-1, 1] * H[-2, -3, 4, 5] * A2[2, 5, 3] * rhoR[3, -4] * A1[1, 4, 2]

        @test HrA12′ == HrA12′′ == HrA12′′′ # should be exactly equal
        @test HrA12 ≈ HrA12′
        @test HrA12 ≈ HrA12′′′′
        @test HrA12′′ == ncon([rhoL, H, A2, rhoR, A1],
                                [[-1,1],[-2,-3,4,5],[2,5,3],[3,-4],[1,4,2]])
        @test HrA12′′ == @ncon([rhoL, H, A2, rhoR, A1],
                                [[-1,1],[-2,-3,4,5],[2,5,3],[3,-4],[1,4,2]];
                                order = [1,2,3,4,5], output=[-1,-2,-3,-4])
        @test E ≈ @tensor scalar(rhoL[a', a] * A1[a, s, b] * A2[b, s', c] * rhoR[c, c'] * H[t, t', s, s'] * conj(A1[a', t, b']) * conj(A2[b', t', c']))
    end
    println("tensor network examples: $(time()-t0) seconds")
    t0 = time()

    @testset "tensoropt" begin
        A=randn(5,5,5,5)
        B=randn(5,5,5)
        C=randn(5,5,5)
        @tensoropt (a=>χ,b=>χ^2,c=>2*χ,d=>χ,e=>5,f=>2*χ) D1[a,b,c,d] := A[a,e,c,f]*B[g,d,e]*C[g,f,b]
        @tensoropt (a=χ,b=χ^2,c=2*χ,d=χ,e=5,f=2*χ) D2[a,b,c,d] := A[a,e,c,f]*B[g,d,e]*C[g,f,b]
        @tensoropt ((a,d)=>χ,b=>χ^2,(c,f)=>2*χ,e=>5) D3[a,b,c,d] := A[a,e,c,f]*B[g,d,e]*C[g,f,b]
        @tensoropt ((a,d)=χ,b=χ^2,(c,f)=2*χ,e=5) D4[a,b,c,d] := A[a,e,c,f]*B[g,d,e]*C[g,f,b]
        @test D1 == D2 == D3 == D4
        _optdata = optex -> TensorContractionsXD.optdata(optex, :(D1[a,b,c,d] := A[a,e,c,f]*B[g,d,e]*C[g,f,b]))
        optex1 = :((a=>χ,b=>χ^2,c=>2*χ,d=>χ,e=>5,f=>2*χ))
        optex2 = :((a=χ,b=χ^2,c=2*χ,d=χ,e=5,f=2*χ))
        optex3 = :(((a,d)=>χ,b=>χ^2,(c,f)=>2*χ,e=>5))
        optex4 = :(((a,d)=χ,b=χ^2,(c,f)=2*χ,e=5))
        optex5 = :(((a,)=>χ,b=>χ^2,(c,)=>2χ,d=>χ,e=>5,f=>χ*2,()=>12345))
        @test _optdata(optex1) == _optdata(optex2) == _optdata(optex3) == _optdata(optex4) == _optdata(optex5)
        optex6 = :(((a,b,c)=χ,))
        optex7 = :((a,b,c))
        @test _optdata(optex6) == _optdata(optex7)
        optex8 = :(((a,b,c)=1,(d,e,f,g)=χ))
        optex9 = :(!(a,b,c))
        @test _optdata(optex8) == _optdata(optex9)
        optex10 = :((a=>χ,b=>χ^2,c=2*χ,d=>χ,e=>5,f=2*χ))
        optex11 = :((a=χ,b=χ^2,c=2*χ,d,e=5,f))
        @test_throws ErrorException _optdata(optex10)
        @test_throws ErrorException _optdata(optex11)
    end

    @testset "Issue 83" begin
        op1 = randn(2, 2)
        op2 = randn(2, 2)
        op3 = randn(2, 2)

        f83(op,op3) = @ncon((op, op3), ([-1 3], [3 -3]))

        b = f83(op1,op3)
        bcopy = deepcopy(b)
        c = f83(op2,op3)
        @test b == bcopy
        @test b != c
    end

    # diagonal matrices
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        A = randn(T, 10, 10, 10, 10)
        @tensor C[3,1,4,2] := A[1,2,3,4]
        U,S,V = svd(reshape(C, (100, 100)))
        U2 = reshape(view(U, 1:2:100, :), (5,10,100))
        S2 = Diagonal(S)
        V2 = reshape(view(V, 1:2:100, :), (5, 10, 100))
        @tensor A2[1,2,3,4] := U2[3,1,a]*S2[a,a']*conj(V2[4,2,a'])
        @test A2 ≈ view(A, :, :, 1:2:10, 1:2:10)

        B = randn(T, 10)
        @tensor C1[3,1,4,2] := A[a,2,3,4]*Diagonal(B)[a,1]
        @tensor C2[3,1,4,2] := A[1,a,3,4]*conj(Diagonal(B)[a,2])
        @tensor C3[3,1,4,2] := conj(A[1,2,a,4])*Diagonal(B)[a,3]
        @tensor C4[3,1,4,2] := conj(A[1,2,3,a])*conj(Diagonal(B)[a,4])
        @tensor C1′[3,1,4,2] := Diagonal(B)[a,1]*A[a,2,3,4]
        @tensor C2′[3,1,4,2] := conj(Diagonal(B)[a,2])*A[1,a,3,4]
        @tensor C3′[3,1,4,2] := Diagonal(B)[a,3]*conj(A[1,2,a,4])
        @tensor C4′[3,1,4,2] := conj(Diagonal(B)[a,4])*conj(A[1,2,3,a])
        @test C1 == C1′
        @test C2 == C2′
        @test C3 == C3′
        @test C4 == C4′
        for i = 1:10
            @test C1[:,i,:,:] ≈ permutedims(A[i,:,:,:],(2,3,1))*B[i]
            @test C2[:,:,:,i] ≈ permutedims(A[:,i,:,:],(2,1,3))*conj(B[i])
            @test C3[i,:,:,:] ≈ conj(permutedims(A[:,:,i,:],(1,3,2)))*B[i]
            @test C4[:,:,i,:] ≈ conj(permutedims(A[:,:,:,i],(3,1,2)))*conj(B[i])
        end

        @tensor D1[1,2] := A[1,2,a,b]*Diagonal(B)[a,b]
        @tensor D2[1,2] := A[1,b,2,a]*conj(Diagonal(B)[a,b])
        @tensor D3[1,2] := conj(A[a,2,1,b])*Diagonal(B)[a,b]
        @tensor D4[1,2] := conj(A[a,1,b,2])*conj(Diagonal(B)[a,b])
        @tensor D1′[1,2] := Diagonal(B)[a,b]*A[1,2,a,b]
        @tensor D2′[1,2] := conj(Diagonal(B)[a,b])*A[1,b,2,a]
        @tensor D3′[1,2] := Diagonal(B)[a,b]*conj(A[a,2,1,b])
        @tensor D4′[1,2] := conj(Diagonal(B)[a,b])*conj(A[a,1,b,2])
        @test D1 == D1′
        @test D2 == D2′
        @test D3 == D3′
        @test D4 == D4′

        E1 = zero(D1)
        E2 = zero(D2)
        E3 = zero(D3)
        E4 = zero(D4)
        for i = 1:10
            E1[:,:] += A[:,:,i,i]*B[i]
            E2[:,:] += A[:,i,:,i]*conj(B[i])
            E3[:,:] += A[i,:,:,i]'*B[i]
            E4[:,:] += conj(A[i,:,i,:])*conj(B[i])
        end
        @test D1 ≈ E1
        @test D2 ≈ E2
        @test D3 ≈ E3
        @test D4 ≈ E4

        F = randn(T,(10,10))
        @tensor G[a,c,b,d] := F[a,b]*Diagonal(B)[c,d]
        @test reshape(G,(100,100)) ≈ kron(Diagonal(B), F)
    end
    println("diagonal examples: $(time()-t0) seconds")
    t0 = time()

end


text:
using HorizonSideRobots
inverse(side::HorizonSide) = HorizonSide(mod(Int(side) + 2, 4))

function find_prohod(r::Robot, side::HorizonSide)
    counter = 1
    while (isborder(r, Nord))==true
        for i in 1:counter
            move!(r, side)
            if !(isborder(r, Nord))
                break
            end
        end

        if !(isborder(r, Nord))
            break
        end

        counter += 1

        for i in 1:counter
            move!(r, inverse(side))
            if !(isborder(r, Nord))
                break
            end
        end
        counter += 1
    end
end

find_prohod(r, Ost)


text:
"""
This creates an enum which details how extrapolation from the interpolation domain should be done.
The possible enum values are:

* `Curve` - Curve extrapolation extends out the quadratic form at the edges. This can lead to a nonmonotonic result (as the curve can eventually change direction)
* `Linear` - Linear extrapolation extends out the gradient from at the edges. This will always lead to a monotonic result.
* `Constant` - Linear extrapolation extends out the value from at the edges. This leads to flat values being extended out.
"""
@enum Schumaker_ExtrapolationSchemes begin
    Curve = 0
    Linear = 1
    Constant = 2
end

"""
    Schumaker(x::Array{T,1},y::Array{T,1} ; gradients::Union{Missing,Array{T,1}} = missing, extrapolation::Schumaker_ExtrapolationSchemes = Curve,
                  left_gradient::Union{Missing,T} = missing, right_gradient::Union{Missing,T} = missing)
    Schumaker(x::Array{Int,1},y::Array{T,1} ; gradients::Union{Missing,Array{T,1}} = missing, extrapolation::Schumaker_ExtrapolationSchemes = Curve,
                  left_gradient::Union{Missing,T} = missing, right_gradient::Union{Missing,T} = missing)
    Schumaker(x::Array{Date,1},y::Array{T,1} ; gradients::Union{Missing,Array{T,1}} = missing, extrapolation::Schumaker_ExtrapolationSchemes = Curve,
                  left_gradient::Union{Missing,T} = missing, right_gradient::Union{Missing,T} = missing)
Creates a Schumaker spline.
### Inputs
* `x` - A vector of x coordinates.
* `y` - A vector of y coordinates.
* `extrapolation` - This should be `Curve`, `Linear` or `Constant` specifying how to interpolate outside of the sample domain.
* `gradients` - A vector of gradients at each point. If not supplied these are imputed from x and y.
* `left_gradient` - The gradient at the lowest value of x in the domain. This will override the gradient imputed or submitted in the gradients optional argument (if it is submitted there)
* `right_gradient` - The gradient at the highest value of x in the domain. This will override the gradient imputed or submitted in the gradients optional argument (if it is submitted there)

### Returns
* A `Schumaker` object which contains the spline. This object can then be evaluated with evaluate or evaluate_integral.
 """
struct Schumaker{T<:AbstractFloat}
    IntStarts_::Array{T,1}
    coefficient_matrix_::Array{T,2}
    function Schumaker(x::Array{T,1},y::Array{<:Real,1} ; gradients::Union{Missing,Array{<:Real,1}} = missing, left_gradient::Union{Missing,<:Real} = missing, right_gradient::Union{Missing,<:Real} = missing,
                       extrapolation::Tuple{Schumaker_ExtrapolationSchemes,Schumaker_ExtrapolationSchemes} = (Curve,Curve)) where T<:Real
        if length(x) == 0
            error("Zero length x vector is insufficient to create Schumaker Spline.")
        elseif length(x) == 1
            IntStarts = Array{T,1}(x)
            @inbounds SpCoefs = [0 0 y[1]]
            # Note that this hardcodes in constant extrapolation. This is only
            # feasible one as we do not have derivative or curve information.
            return new{T}(IntStarts, SpCoefs)
        elseif length(x) == 2
            @inbounds IntStarts = Array{T,1}([x[1]])
            @inbounds IntEnds   = Array{T,1}([x[2]])
            @inbounds linear_coefficient = (y[2]- y[1]) / (x[2]-x[1])
            @inbounds SpCoefs = [0 linear_coefficient y[1]]
            # In this case it defaults to curve extrapolation (which is same as linear here)
            # So we just alter in case constant is specified.
            if extrapolation[1] == Constant || extrapolation[2] == Constant
                matrix_without_extrapolation = hcat(IntStarts , SpCoefs)
                @inbounds matrix_with_extrapolation    = extrapolate(matrix_without_extrapolation, extrapolation, x[2], y)
                return @inbounds new{T}(matrix_with_extrapolation[:,1], matrix_with_extrapolation[:,2:4])
            else
                return new{T}(IntStarts, SpCoefs)
            end
        end
        if ismissing(gradients)
           gradients = imputeGradients(x,y)
        end
        if !ismissing(left_gradient)
            @inbounds gradients[1] = left_gradient
        end
        if !ismissing(right_gradient)
            @inbounds gradients[length(gradients)] = right_gradient
        end
        IntStarts, SpCoefs = getCoefficientMatrix(x,y,gradients, extrapolation)
        if T<:AbstractFloat
            G = T
        else
            G = Float64
        end
        return new{G}(G.(IntStarts), G.(SpCoefs))
     end
    function Schumaker(IntStarts_::Array{T,1}, coefficient_matrix_::Array{T,2}) where {T<:Real}
        return new{T}(IntStarts_, coefficient_matrix_)
    end
    function Schumaker{Q}(IntStarts_::Array{<:Real,1}, coefficient_matrix_::Array{<:Real,2}) where Q<:Real
        return new{Q}(Q.(IntStarts_), Q.(coefficient_matrix_))
    end
    function Schumaker(x::Array{Date,1},y::Array{<:Real,1} ; gradients::Union{Missing,Array{<:Real,1}} = missing, left_gradient::Union{Missing,<:Real} = missing, right_gradient::Union{Missing,<:Real} = missing,
                       extrapolation::Tuple{Schumaker_ExtrapolationSchemes,Schumaker_ExtrapolationSchemes} = (Curve,Curve))
        days_as_ints = Dates.days.(x)
        T = promote_type(eltype(y), eltype(days_as_ints))
        if T<:AbstractFloat
            G = T
        else
            G = Float64
        end
        return Schumaker{G}(days_as_ints , y; gradients = gradients , extrapolation = extrapolation, left_gradient = left_gradient, right_gradient = right_gradient)
    end
    function Schumaker{Q}(x::AbstractArray, y::AbstractArray ; gradients::Union{Missing,AbstractArray} = missing, left_gradient::Union{Missing,Real} = missing, right_gradient::Union{Missing,Real} = missing,
                          extrapolation::Tuple{Schumaker_ExtrapolationSchemes,Schumaker_ExtrapolationSchemes} = (Curve,Curve)) where Q<:Real
        got_both = (!).(ismissing.(x) .| ismissing.(y))
        @inbounds new_x = convert.(Q, x[got_both])
        @inbounds new_y = convert.(Q, y[got_both])
        @inbounds new_gradients = ismissing(gradients)        ? missing : convert.(Q, gradients[got_both])
        converted_left  = ismissing(left_gradient)  ? missing : convert(Q, left_gradient)
        converted_right = ismissing(right_gradient) ? missing : convert(Q, right_gradient)
        new_left  = got_both[1]                     ? converted_left  : missing
        new_right = got_both[length(got_both)]      ? converted_right : missing
        if length(new_x) == 0 error("After removing missing elements there are no points left to estimate schumaker spline") end
        return Schumaker(new_x , new_y; gradients = new_gradients , extrapolation = extrapolation, left_gradient = new_left, right_gradient = new_right)
    end
    function Schumaker(x::Union{AbstractArray{T,1},AbstractArray{Union{Missing,T},1}},y::Union{AbstractArray{R,1},AbstractArray{Union{Missing,R},1}} ; gradients::Union{Missing,AbstractArray{<:Real,1}} = missing, left_gradient::Union{Missing,Real} = missing, right_gradient::Union{Missing,Real} = missing,
                       extrapolation::Tuple{Schumaker_ExtrapolationSchemes,Schumaker_ExtrapolationSchemes} = (Curve,Curve)) where T<:Real where R<:Real
       promo_type = promote_type(T,R)
       return Schumaker{promo_type}(x, y; gradients = gradients, extrapolation = extrapolation, left_gradient = left_gradient, right_gradient = right_gradient)
   end
end
Base.broadcastable(e::Schumaker) = Ref(e)

"""
Evaluates the spline at a point. The point can be specified as a Real number (Int, Float, etc) or a Date.
Derivatives can also be taken.
### Inputs
 * `spline` - A `Schumaker` type spline
 * `PointToExamine` - The point at which to evaluate the integral
 * `derivative` - The derivative being sought. This should be 0 to just evaluate the spline, 1 for the first derivative or 2 for a second derivative.
 Higher derivatives are all zero (because it is a quadratic spline). Negative values do not give integrals. Use evaluate_integral instead.
### Returns
 * A value of the spline or appropriate derivative in the same format as specified in the spline.
"""
function evaluate(spline::Schumaker, PointToExamine::T,  derivative::Integer = 0) where T<:Real
    # Derivative of 0 means normal spline evaluation.
    # Derivative of 1, 2 are first and second derivatives respectively.
    IntervalNum = searchsortedlast(spline.IntStarts_, PointToExamine)
    IntervalNum = max(IntervalNum, 1)
    @inbounds xmt = AbstractFloat(PointToExamine - spline.IntStarts_[IntervalNum])
    Coefs = @inbounds spline.coefficient_matrix_[IntervalNum,:]
    if derivative == 0
        return sum(@. Coefs .* [xmt^2, xmt, 1.0])
    elseif derivative == 1
        return sum(@. Coefs .* [2*xmt, 1.0, 0.0])
    elseif derivative == 2
        return sum(@. Coefs .* [2.0, 0.0, 0.0])
    elseif derivative < 0
        error("This function cannot do integrals. Use evaluate_integral instead")
    else
        return 0.0
    end
end
function evaluate(spline::Schumaker, PointToExamine::Date,  derivative::Int = 0)
    days_as_int = Dates.days.(PointToExamine)
    return evaluate(spline,days_as_int,  derivative)
end
function (s::Schumaker)(x::Union{Real,Date}, deriv::Integer = 0)
    return evaluate(s, x, deriv)
end

"""
Estimates the integral of the spline between lhs and rhs. These end points can be input
as Reals or Dates.
### Inputs
 * `spline` - A Schumaker type spline
 * `lhs` - The left hand limit of the integral
 * `rhs` - The right hand limit of the integral
### Returns
 * A `Float64` value of the integral.
"""
function evaluate_integral(spline::Schumaker, lhs::Real, rhs::Real)
    first_interval = searchsortedlast(spline.IntStarts_, lhs)
    last_interval = searchsortedlast(spline.IntStarts_, rhs)
    number_of_intervals = last_interval - first_interval
    if number_of_intervals == 0
        return section_integral(spline , lhs, rhs)
    elseif number_of_intervals == 1
        @inbounds first = section_integral(spline , lhs , spline.IntStarts_[first_interval + 1])
        @inbounds lst  = section_integral(spline , spline.IntStarts_[last_interval]  , rhs)
        return first + lst
    else
        interior_areas = 0.0
        @inbounds first = section_integral(spline , lhs , spline.IntStarts_[first_interval + 1])
        @simd for i in 1:(number_of_intervals-1)
            @inbounds sec_int =  section_integral(spline ,  spline.IntStarts_[first_interval + i] , spline.IntStarts_[first_interval + i+1] )
            interior_areas += sec_int
        end
        @inbounds last  = section_integral(spline , spline.IntStarts_[last_interval]  , rhs)
        return first + interior_areas + last
    end
end
function evaluate_integral(spline::Schumaker, lhs::Date, rhs::Date)
    return evaluate_integral(spline, Dates.days.(lhs) , Dates.days.(rhs))
end
function section_integral(spline::Schumaker, lhs::Real,  rhs::Real)
    # Note that the lhs is used to infer the interval.
    IntervalNum = searchsortedlast(spline.IntStarts_, lhs)
    IntervalNum = max(IntervalNum, 1)
    @inbounds Coefs = spline.coefficient_matrix_[ IntervalNum , :]
    @inbounds r_xmt = rhs - spline.IntStarts_[IntervalNum]
    @inbounds l_xmt = lhs - spline.IntStarts_[IntervalNum]

    Lint_array =  [(1/3)*l_xmt^3, 0.5*l_xmt^2, l_xmt]
    Rint_array =  [(1/3)*r_xmt^3, 0.5*r_xmt^2, r_xmt]
    return sum(@. Coefs .* Rint_array) - sum(@. Coefs .* Lint_array)
end

"""
find_derivative_spline(spline::Schumaker)
    Returns a SchumakerSpline that is the derivative of the input spline
"""
function find_derivative_spline(spline::Schumaker)
    coefficient_matrix = Array{Float64,2}(undef, size(spline.coefficient_matrix_)...   )
    coefficient_matrix[:,3] .= spline.coefficient_matrix_[:,2]
    coefficient_matrix[:,2] .= 2 .* spline.coefficient_matrix_[:,1]
    coefficient_matrix[:,1] .= 0.0
    return Schumaker(spline.IntStarts_, coefficient_matrix)
end

"""
    imputeGradients(x::Vector{T}, y::Vector{T})
Imputes gradients based on a vector of x and y coordinates.
"""
function imputeGradients(x::Vector{<:Real}, y::Vector{<:Real})
     n = length(x)
     # Judd (1998), page 233, second last equation
     @inbounds L = sqrt.( (x[2:n]-x[1:(n-1)]).^2 + (y[2:n]-y[1:(n-1)]).^2)
     # Judd (1998), page 233, last equation
     @inbounds d = (y[2:n]-y[1:(n-1)])./(x[2:n]-x[1:(n-1)])
     # Judd (1998), page 234, Eqn 6.11.6
     @inbounds Conditionsi = d[1:(n-2)].*d[2:(n-1)] .> 0
     @inbounds MiddleSiwithoutApplyingCondition = (L[1:(n-2)].*d[1:(n-2)]+L[2:(n-1)].* d[2:(n-1)]) ./ (L[1:(n-2)]+L[2:(n-1)])
     sb = Conditionsi .* MiddleSiwithoutApplyingCondition
     # Judd (1998), page 234, Second Equation line plus 6.11.6 gives this array of slopes.
     @inbounds ff = [((-sb[1]+3*d[1])/2);  sb ;  ((3*d[n-1]-sb[n-2])/2)]
     return ff
 end

"""
Splits an interval into 2 subintervals and creates the quadratic coefficients
### Inputs
 * `gradients` - A 2 entry vector with gradients at either end of the interval
 * `y` - A 2 entry vector with y values at either end of the interval
 * `x` - A 2 entry vector with x values at either end of the interval
### Returns
 * A 2 x 4 matrix. The first column is the x values of start of the two subintervals. The last 3 columns are quadratic coefficients in two subintervals.
"""
function schumakerIndInterval(gradients::Vector{<:Real}, y::Vector{<:Real}, x::Vector{<:Real}) where T<:Real
   # The SchumakerIndInterval function takes in each interval individually
   # and returns the location of the knot as well as the quadratic coefficients in each subinterval.

   # Judd (1998), page 232, Lemma 6.11.1 provides this if condition:
   if (sum(gradients)*(x[2]-x[1]) == 2*(y[2]-y[1]))
     tsi = x[2]
   else
     # Judd (1998), page 233, Algorithm 6.3 along with equations 6.11.4 and 6.11.5 provide this whole section
     @inbounds delta = (y[2] -y[1])/(x[2]-x[1])
     @inbounds Condition = ((gradients[1]-delta)*(gradients[2]-delta) >= 0)
     @inbounds Condition2 = abs(gradients[2]-delta) < abs(gradients[1]-delta)
     if (Condition)
       tsi = sum(x)/2
     elseif (Condition2)
       @inbounds tsi = (x[1] + (x[2]-x[1])*(gradients[2]-delta)/(gradients[2]-gradients[1]))
     else
       @inbounds  tsi = (x[2] + (x[2]-x[1])*(gradients[1]-delta)/(gradients[2]-gradients[1]))
     end
   end

   # Judd (1998), page 232, 3rd last equation of page.
   alpha = tsi-x[1]
   beta = x[2]-tsi
   # Judd (1998), page 232, 4th last equation of page.
    @inbounds sbar = (2*(y[2]-y[1])-(alpha*gradients[1]+beta*gradients[2]))/(x[2]-x[1])
   # Judd (1998), page 232, 3rd equation of page. (C1, B1, A1)
    @inbounds Coeffs1 = [ (sbar-gradients[1])/(2*alpha)  gradients[1]  y[1] ]
   if (beta == 0)
     Coeffs2 = Coeffs1
   else
     # Judd (1998), page 232, 4th equation of page. (C2, B2, A2)
     @inbounds Coeffs2 = [ (gradients[2]-sbar)/(2*beta)  sbar  Coeffs1 * [alpha^2, alpha, 1] ]
   end
   Machine4Epsilon = 4*eps()
     if (tsi  <  x[1] + Machine4Epsilon )
         return [x[1] Coeffs2]
     elseif (tsi + Machine4Epsilon > x[2] )
         return [x[1]  Coeffs1]
     else
         return [x[1] Coeffs1 ; tsi Coeffs2]
     end
 end

 """
 Calls `SchumakerIndInterval` many times to get full set of spline intervals and coefficients. Then calls extrapolation for out of sample behaviour
### Inputs
 * `gradients` - A vector of gradients at each point
 * `x` - A vector of `x` coordinates
 * `y` - A vector of `y` coordinates
 * `extrapolation` - A string in (Curve, Linear or Constant) that gives behaviour outside of interpolation range.
### Returns
 * A vector of interval starts
 * A vector of interval ends
 * A matrix of all coefficients
  """
 function getCoefficientMatrix(x::Array{<:Real,1}, y::Array{<:Real,1}, gradients::Array{<:Real,1},
                               extrapolation::Tuple{Schumaker_ExtrapolationSchemes,Schumaker_ExtrapolationSchemes})
   n = length(x)
   fullMatrix = schumakerIndInterval([gradients[1], gradients[2]], [y[1], y[2]], [x[1], x[2]] )
    for intrval = 2:(n-1)
      xs = [ x[intrval] , x[intrval + 1] ]
      ys = [ y[intrval], y[intrval + 1] ]
      grads = [ gradients[intrval], gradients[intrval + 1] ]
      intMatrix = schumakerIndInterval(grads,ys,xs)
      fullMatrix = vcat(fullMatrix,intMatrix)
    end
    fullMatrix = extrapolate(fullMatrix, extrapolation, x[n], y)
   return fullMatrix[:,1], fullMatrix[:,2:4]
 end

"""
 Adds a row on top and bottom of coefficient matrix to give out of sample prediction.
### Inputs
 * `fullMatrix` - output from `GetCoefficientMatrix` first few lines
 * `extrapolation` - A tuple with two enums in (Curve, Linear, Constant) that gives behaviour outside of interpolation range.
 * `x` - A vector of x coordinates
 * `y` - A vector of y coordinates
### Returns
  * A new version of fullMatrix with out of sample prediction built into it.
"""
function extrapolate(fullMatrix::Array{<:Real,2}, extrapolation::Tuple{Schumaker_ExtrapolationSchemes,Schumaker_ExtrapolationSchemes}, Topx::Real, y::Array{<:Real,1})
    # In the fullMatrix the first column are the x values and then the next 3 columns are a, b, c in the expression a(x-start)^2 + b(x-start) + c
  if (extrapolation[1] == Curve) && (extrapolation[2] == Curve)
    return fullMatrix
  end
  # Preliminaries used throughout
  dim = size(fullMatrix)[1]
  Botx   = fullMatrix[1,1]
  Boty   = y[1]
  Topy = y[length(y)]
  # Initialising variable so their domain is not restricted to the if statement spaces.
  BotB = 0.0
  BotC = 0.0
  TopA = 0.0
  TopB = 0.0
  TopC = 0.0
  # Now creating the extrapolation to the left.
  if extrapolation[1] == Linear
    BotB = fullMatrix[1 , 3]
    BotC   = Boty - BotB
  elseif extrapolation[1] == Constant
    BotB = 0.0
    BotC = Boty
  end # Note for the curve case we will simply not append the new block.
  BotRow = [Botx-1e-10, 0.0, BotB, BotC]

  # Now doing the extrapolation to the right.
  if extrapolation[2] == Linear # This is a bit more complicated than before because the
                                # coefficients by themselves give the gradients at the left
                                # of the interval. Here we want the gradient at the right.
    last_interval_width = Topx - fullMatrix[dim , 1]
    @inbounds grad_at_right = 2 * fullMatrix[dim , 2] * last_interval_width + fullMatrix[dim , 3]
    TopB = grad_at_right
    TopC = Topy
  elseif extrapolation[2] == Constant
    TopB = 0.0
    TopC = Topy
  elseif extrapolation[2] == Curve # We are just going to add this one on regardless as otherwise end of data interval information is lost.
    Gap  =  Topx - fullMatrix[dim ,1]
    TopA = fullMatrix[dim ,2]
    @inbounds TopB = fullMatrix[dim ,3] + 2 * TopA * Gap
    @inbounds TopC = fullMatrix[dim ,4] + TopB * Gap - TopA * Gap*Gap
  end
  TopRow = [ Topx, TopA ,TopB ,TopC]
  # Appending blocks and returning.
  if extrapolation[1] != Curve fullMatrix = vcat(BotRow' , fullMatrix) end
  fullMatrix = vcat(fullMatrix,  TopRow')
  return fullMatrix
end


text:
using GeoStats
using SpectralGaussianSimulation
using SparseArrays
using Plots
using Random
using BlockArrays
using Distances
using StatsFuns
using JLD

include("./SortSparse.jl")
include("./KoLesky.jl")
include("./CovFuncs.jl")
include("./Utils.jl")

Random.seed!(123)

n = 1000
N = n^2
NTest = 20000
NTrain = N - NTest
NSamples = 1000

ρList = [1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8,0, 8.5, 9.0, 9.5, 10.0]
λ = 1.0
l = 0.5
covfunc!(r) = matern32!(r, l)
vario = MaternVariogram(range=l, order=1.5)
padFact = 5

# Prepare output Arrays
RMSE = zeros(length(ρList))
RMSE_int = zeros(length(ρList))
RMSE_exp = zeros(length(ρList))
coverage = zeros(length(ρList))
coverage_int = zeros(length(ρList))
coverage_exp = zeros(length(ρList))

# Sampling the process
if (isfile("./out/jld/FFTData_n_$n.jld"))
  println("Loading FFT data from file")
  ld = load("./out/jld/FFTData_n_$n.jld")
  x = ld["x"]
  Y = ld["Y"]
else
  println("FFT data not found. Computing it from scratch")
  @time x, Y = specSampleGrid2D(n, NSamples, padFact, vario)
  save("./out/jld/FFTData_n_$n.jld", "x", x, "Y", Y)
end

# random shuffle of indices
shuffleIDs = shuffle(1:N)
# Selecting groups of nodes as test data
innerInds = [i for i in 1:N if (norm(x[:, i] .- 0.3) < 0.05 || norm((x[:,i] .- 0.75) .* [1.0, 0.1]) < 0.01)]
shuffleIDs = vcat(shuffleIDs, innerInds)
shuffleIDs = unique(shuffleIDs[end:-1:1])[end:-1:1]

# Creating data to be used by each run of for loop
#Reordering data
xBase = x[:, shuffleIDs]
YBase = Y[shuffleIDs, :]
#Split data into train and testset
xTrainBase = xBase[:, 1 : NTrain]
xTestBase = xBase[:, (NTrain + 1) : end]

n_extrapolation = min(NTest, length(innerInds))
interpolation_indices_base = (NTrain + 1) : (NTrain + NTest - n_extrapolation)
extrapolation_indices_base = (NTrain + NTest - n_extrapolation + 1) : N

for (ρInd, ρ) in enumerate(ρList)
  x = copy(xBase)
  Y = copy(YBase)
  xTrain = copy(xTrainBase)
  xTest = copy(xTestBase)
  global RMSE
  global RMSE_int
  global RMSE_exp
  global coverage
  global coverage_int
  global coverage_exp

  colptr, rowval, P, revP, distances = sortSparseRev(xTrain, xTest, ρ, 1 )
  x = x[:, P]

  Y = Y[P, :]
  Y = PseudoBlockArray(Y, [NTest, NTrain], [NSamples])

  xTest = x[:, 1 : NTest]
  xTrain = x[:, (NTest + 1) : end]

  interpolation_indices = revP[interpolation_indices_base]
  extrapolation_indices = revP[extrapolation_indices_base]

  skeletons = construct_skeletons(colptr, rowval, distances, λ)
  @time L = assembleL(skeletons, x, covfunc!)
  L = PseudoBlockArray(L, [NTest, NTrain], [NTest, NTrain])

  # Compute posterior covariance
  PΘApprox = Matrix(Hermitian(inv(Matrix(L[Block(1), Block(1)] * L[Block(1), Block(1)]'))))
  # Compute posterior mean
  μApprox = -L[Block(1), Block(1)]' \ (L[Block(2), Block(1)]' * Y[Block(2), :])

  # Computing coverage of the 0.9 confidence interval
  α = 0.1
  coverageMatrix = Float64.(abs.(μApprox - Y[Block(1), :]) ./ sqrt.(diag(PΘApprox)) .<= norminvcdf(1.0 - α/2))

  RMSE[ρInd] = sqrt(mean((μApprox .- Y[Block(1), :]).^2))
  RMSE_int[ρInd] = sqrt(mean(((μApprox .- Y[Block(1), :]).^2)[interpolation_indices, :]))
  RMSE_exp[ρInd] = sqrt(mean(((μApprox .- Y[Block(1), :]).^2)[extrapolation_indices, :]))

  coverage[ρInd] = mean(coverageMatrix)
  coverage_int[ρInd] = mean(coverageMatrix[interpolation_indices, :])
  coverage_exp[ρInd] = mean(coverageMatrix[extrapolation_indices, :])
end

save("./out/jld/FFTVaryRho_lambda_$(Int(10 * λ)).jld",
     "ρList", ρList,
     "RMSE", RMSE,
     "RMSE_int", RMSE_int,
     "RMSE_exp", RMSE_exp,
     "coverage", coverage,
     "coverage_int", coverage_int,
     "coverage_exp", coverage_exp,
     "xTest", xTestBase)


text:
module DeutschlandDaten
using DataDeps
using CSV
using DataFrames
using StatsBase
using StatsPlots

export tot_cgrad
tot_cgrad = cgrad([:white, :black],[-.1,1]) # :reds

export alterstring, rangestart
alterstring(a) = (iszero(a.start) ? "" : "$(a.start)") * (a.stop==DeutschlandDaten.MAX ? "+" : "-$(a.stop)")
rangestart(x) = x.start

export label_geschlecht
label_geschlecht = Dict("Männlich" => "♂", "Weiblich" => "♀")

export endofweek
function endofweek(y, w)
    d = Dates.firstdayofweek(Date(y)) + Dates.Day(2)
    d + Dates.Week(w - (year(d)<y ? 0 : 1)) + Dates.Day(4)
end

export gliding_mean
"""
    gliding_mean(x, left=6, right=left)

compute gliding mean over x with a window left and right (if data available).
"""
function gliding_mean(x, left::Int=6, right::Int=left)
    y = Vector{Union{Missing,Float64}}(missing, length(x))
    for i in 1:length(x)
        y[i] = mean(skipmissing(x[max(1,i-left):min(end,i+right)]))
    end
    y
end

# data sources
include("metadata.jl")
include("deaths.jl")
include("vaccines.jl")
include("population.jl")

using Plots
using StatsPlots
using Plots.PlotMeasures
using Printf
export plotyears

using Printf
function plotyears(mordat, from_date = minimum(mordat[:,:date]); kw...)
    p1 = plotyears(mordat,"Männlich", from_date; title = "männlich ♂", kw...)
    p2 = plotyears(mordat,"Weiblich", from_date; title ="weiblich ♀", ylabel="", yaxis=(formatter = p -> ""), kw...)
    plot(plot(title = "Sterbekurve Deutschland $(year(from_date))-2021",framestyle=nothing,showaxis=false,xticks=false,yticks=false,margin=0mm), p1,p2, layout=@layout [a{0.01h}; grid(1,2)])
end

function plotyears(mordat, geschlecht::String, from_date = minimum(mordat[:,:date]);
                   title = geschlecht,
                   ydata = :PD,
                   ylabel = "Sterberate pro 10'000 / Monat (log)",
                   yticks = 10000 .* [.00001, .00002, .00005, .0001, .0002,.0005,.001,.002,.005,.01,.02],
                   ylim =  10000 .* (0.00001,0.023),
                   yaxis=(formatter = p -> @sprintf("%.2f",p)),
                   yearticks = (2021-year(from_date) > 10) ? ((year(from_date)):5:2021) : ((year(from_date)):1:2021))
    d = mordat[( mordat[:,:geschlecht] .== geschlecht ) .& (mordat[:,:date] .>= from_date ),:]
    date,maxdate=minimum(d[:,:date]),maximum(d[:,:date])
    @info "background"
    plot([Date(2020,06,1), maxdate],
         color=:yellow, opacity = .25,
         seriestype=:vspan,
         xlim = (date,maxdate), ylim =ylim,
         xticks = yearticks, yticks = yticks,
         yaxis=:log10, title=title, legend=nothing)

    plot!([Date(2020,12,27), maxdate],
         color=:gray, opacity = .25,
         seriestype=:vspan)
    plot!([Date(year,1,1) for year in (year(from_date)+1):2021];
          seriestype=:vline, color=:lightgray)
    @info "PD"
    @df d plot!(:date, :PD; series=:alter, c=:gray, legend=nothing)
    @info "PD_12months"
    p1 = @df d plot!(:date, :PD_12months; group=:alter, width=3,
                     xticks=([Date(y,01,01) for y in yearticks], [y for y in yearticks]),
                     xrotation = -90,
                     # xlabel=nothing,
                     yticks=yticks,
                     ylabel=ylabel,
                     yaxis=yaxis,
                     ylim=ylim,
                     )
    g = first(unique(d.alter))
    annotate!([(maximum(d[:,:date]),
                d[map(x-> x == g,d.alter),:PD_12months][end],
                text("$(g.start)", 6, halign=:left, valign=:vcenter))
               for g in unique(d.alter)])
    p1
end

export plotquelle!
"""
    plotquelle!(a...;
                data = "Datenquellen: Statistisches Bundesamt, Robert-Koch-Institut",
                codepage = "Sterberate.md",
                code = "Quelle: https://github.com/gkappler/DeutschlandDaten/blob/main/\$codepage",
                inset = bbox(.5, 400px, 50px, 50px),
                titlefontsize = 6, margin=0mm,
                title = "\$data\n\$code",
                kw...)

Add a reference for plot updates.
"""
function plotquelle!(a...;
                     data = "Datenquellen: Statistisches Bundesamt, Robert-Koch-Institut",
                     codepage = "Sterberate.md",
                     code = "Quelle: https://github.com/gkappler/DeutschlandDaten/blob/main/$codepage",
                     inset = bbox(.5, 400px, 50px, 50px),
                     titlefontsize = 6, margin=0mm,
                     title = "$data\n$code",
                     kw...)
    plot!(a...; title=title,
          inset = inset,
          titlefontsize = titlefontsize,
          bg_inside = nothing, framestyle=nothing, showaxis=false, xticks=false, yticks=false,
          margin=margin,
          kw...)
end

function savefigs(p, file, quelle=true; kw...)
    quelle && plotquelle!(p; subplot = length(p.subplots)+1, kw...)
    savefig(p,"images/$file.svg")
    savefig(p,"images/$file.png")
end

export plot_unstacked_heatmap
plot_unstacked_heatmap(d; kw...) =
    plot_unstacked_heatmap((p,m)-> (p,m), d; kw...)
function plot_unstacked_heatmap(f::Function, d; labelf, kw...)
    us = unstack(d[:,Symbol[:row,:col,:value]], :col, :value)
    yt = map(labelf, us[:,1])
    p = heatmap(Matrix(us[:,2:end]);
               yticks = (1:length(yt), yt),
               kw...)
    f(p, us)
end

end

max_stars_repo_name: gkappler/DeutschlandDaten

max_issues_repo_name: gkappler/DeutschlandDaten

max_forks_repo_name: gkappler/DeutschlandDaten


text:
module ClimateModels

using UUIDs, Suppressor, OrderedCollections
using Pkg, Git, TOML
import Pkg.Artifacts

p=dirname(pathof(ClimateModels))
artifact_toml = joinpath(p, "../Artifacts.toml")
IPCC_SPM_hash = Artifacts.artifact_hash("IPCC_SPM", artifact_toml)
IPCC_SPM_path = Artifacts.artifact_path(IPCC_SPM_hash)

export AbstractModelConfig, ModelConfig
export clean, build, compile, setup, launch
export put!, take!
export pathof, readdir, log
#export git_log_init, git_log_msg, git_log_fil
#export git_log_prm, git_log_show
#export monitor, help, pause
#export train, compare, analyze
export cmip
export IPCC_SPM_path, IPCC_hexagons #should be commented out?
export OrderedDict, UUID, uuid4, @suppress #should be commented out?

include("interface.jl")
include("access.jl")
include("toy_models.jl")

end


text:
using FromFile
@from "Core.jl" import Node, copyNode, Options, Dataset
@from "EquationUtils.jl" import countNodes, countConstants, countDepth
@from "LossFunctions.jl" import scoreFunc, scoreFuncBatch
@from "CheckConstraints.jl" import check_constraints
@from "PopMember.jl" import PopMember
@from "MutationFunctions.jl" import genRandomTree, mutateConstant, mutateOperator, appendRandomOp, prependRandomOp, insertRandomOp, deleteRandomOp
@from "SimplifyEquation.jl" import simplifyTree, combineOperators, simplifyWithSymbolicUtils

# Go through one simulated options.annealing mutation cycle
#  exp(-delta/T) defines probability of accepting a change
function nextGeneration(dataset::Dataset{T},
                        baseline::T, member::PopMember, temperature::T,
                        curmaxsize::Int, frequencyComplexity::AbstractVector{T},
                        options::Options)::PopMember where {T<:Real}

    prev = member.tree
    tree = prev
    #TODO - reconsider this
    if options.batching
        beforeLoss = scoreFuncBatch(dataset, baseline, prev, options)
    else
        beforeLoss = member.score
    end

    nfeatures = dataset.nfeatures

    mutationChoice = rand()
    #More constants => more likely to do constant mutation
    weightAdjustmentMutateConstant = min(8, countConstants(prev))/8.0
    cur_weights = copy(options.mutationWeights) .* 1.0
    cur_weights[1] *= weightAdjustmentMutateConstant
    n = countNodes(prev)
    depth = countDepth(prev)

    # If equation too big, don't add new operators
    if n >= curmaxsize || depth >= options.maxdepth
        cur_weights[3] = 0.0
        cur_weights[4] = 0.0
    end
    cur_weights /= sum(cur_weights)
    cweights = cumsum(cur_weights)

    successful_mutation = false
    #TODO: Currently we dont take this \/ into account
    is_success_always_possible = true
    attempts = 0
    max_attempts = 10

    #############################################
    # Mutations
    #############################################
    while (!successful_mutation) && attempts < max_attempts
        tree = copyNode(prev)
        successful_mutation = true
        if mutationChoice < cweights[1]
            tree = mutateConstant(tree, temperature, options)

            is_success_always_possible = true
            # Mutating a constant shouldn't invalidate an already-valid function

        elseif mutationChoice < cweights[2]
            tree = mutateOperator(tree, options)

            is_success_always_possible = true
            # Can always mutate to the same operator

        elseif mutationChoice < cweights[3]
            if rand() < 0.5
                tree = appendRandomOp(tree, options, nfeatures)
            else
                tree = prependRandomOp(tree, options, nfeatures)
            end
            is_success_always_possible = false
            # Can potentially have a situation without success
        elseif mutationChoice < cweights[4]
            tree = insertRandomOp(tree, options, nfeatures)
            is_success_always_possible = false
        elseif mutationChoice < cweights[5]
            tree = deleteRandomOp(tree, options, nfeatures)
            is_success_always_possible = true
        elseif mutationChoice < cweights[6]
            tree = simplifyTree(tree, options) # Sometimes we simplify tree
            tree = combineOperators(tree, options) # See if repeated constants at outer levels
            # SymbolicUtils is quite slow, so only rarely
            #  do we use it for simplification.
            if rand() < 0.01
                tree = simplifyWithSymbolicUtils(tree, options, curmaxsize)
            end
            return PopMember(tree, beforeLoss)

            is_success_always_possible = true
            # Simplification shouldn't hurt complexity; unless some non-symmetric constraint
            # to commutative operator...

        elseif mutationChoice < cweights[7]
            tree = genRandomTree(5, options, nfeatures) # Sometimes we generate a new tree completely tree

            is_success_always_possible = true
        else # no mutation applied
            return PopMember(tree, beforeLoss)
        end

        successful_mutation = successful_mutation && check_constraints(tree, options, curmaxsize)

        attempts += 1
    end
    #############################################

    if !successful_mutation
        return PopMember(copyNode(prev), beforeLoss)
    end

    if options.batching
        afterLoss = scoreFuncBatch(dataset, baseline, tree, options)
    else
        afterLoss = scoreFunc(dataset, baseline, tree, options)
    end

    if isnan(afterLoss)
        return PopMember(copyNode(prev), beforeLoss)
    end

    probChange = 1.0
    if options.annealing
        delta = afterLoss - beforeLoss
        probChange *= exp(-delta/(temperature*options.alpha))
    end
    if options.useFrequency
        oldSize = countNodes(prev)
        newSize = countNodes(tree)
        probChange *= frequencyComplexity[oldSize] / frequencyComplexity[newSize]
    end

    if probChange < rand()
        return PopMember(copyNode(prev), beforeLoss)
    else
        return PopMember(tree, afterLoss)
    end
end


text:
# There is very litte written against the abstract type.
# So this is of limited use.

module Z4Groups

import Random

export Z4Group, z4group, ConvertableNumber, AbstractZ4Group

abstract type AbstractZ4Group end

Base.one(n::AbstractZ4Group) = one(typeof(n))

const ConvertableNumber = Union{Integer, Complex{<:Integer}, Float64, Complex{<:Float64}}

"""
    struct Z4Group

Type representing the group `{i, -1, -i, 1}`
"""
struct Z4Group <: AbstractZ4Group
    imag::Bool
    minus::Bool
end

const Z4GROUP_i = Z4Group(1, 0)
const Z4GROUP_minus_one = Z4Group(0, 1)
const Z4GROUP_minus_i = Z4Group(1, 1)
const Z4GROUP_one = Z4Group(0, 0)

#####
##### Constructors
#####

"""
    z4group(n::Integer)::Z4Group

Perform one-based index `n` into `(i, -1, -i, 1)`.
"""
z4group(n::Integer) = (Z4GROUP_i, Z4GROUP_minus_one, Z4GROUP_minus_i, Z4GROUP_one)[n]

function Random.rand(rng::Random.AbstractRNG, ::Random.SamplerType{Z4Group})
    return z4group(rand(rng, 1:4))
end

#####
##### Conversion
#####

Base.promote_type(::Type{Z4Group}, ::Type{T}) where {T <: Real} = Complex{T}
Base.promote_type(::Type{T}, ::Type{Z4Group}) where {T <: Real} = Complex{T}
Base.promote_type(::Type{Complex{T}}, ::Type{Z4Group}) where {T <: Real} = Complex{T}

Z4Group(x::ConvertableNumber) = convert(Z4Group, x)
function Base.convert(::Type{Z4Group}, x::ConvertableNumber)
    if isreal(x)
        isone(x) && return one(Z4Group)
        isone(-x) && return -one(Z4Group)
    else iszero(x.re)
        isone(x.im) && return Z4Group(true, false)
        isone(-x.im) && return Z4Group(true, true)
    end
    throw(InexactError(:convert, Z4Group, x))
end

function Base.convert(::Type{T}, x::Z4Group) where {T<:Real}
    if ~x.imag
        if x.minus
            return -one(T)
        else
            return one(T)
        end
    else
        throw(InexactError(:convert, T, x))
    end
end

function Base.convert(::Type{Complex{T}}, x::Z4Group) where {T<:Real}
    if isreal(x)
        return Complex(convert(T, x), zero(T))
    end
    if x.minus
        return Complex(zero(T), -one(T))
    else
        return Complex(zero(T), one(T))
    end
end

####
#### IO
####

function Base.show(io::IO, ip::Z4Group)
    if ip.minus
        print(io, "-")
    else
        print(io, "+")
    end
    if ip.imag
        print(io, "i")
    else
        print(io, "1")
    end
end

Base.show(::MIME"text/input", ip::Z4Group) = show(stdout, MIME"text/input"(), ip)
function Base.show(io::IO, ::MIME"text/input", ip::Z4Group)
    print(io, "Z4Group(")
    if ip.minus
        print(io, "-")
    else
        print(io, "+")
    end
    if ip.imag
        print(io, "im")
    else
        print(io, "1")
    end
    print(io, ")")
end

####
#### Container interface
####

Base.length(::Z4Group) = 1
Base.iterate(n::Z4Group) = (n, nothing)
Base.iterate(n::Z4Group, ::Nothing) = nothing

####
#### Compare / predicates
####

Base.one(::Type{Z4Group}) = Z4Group(false, false)
Base.isreal(x::Z4Group) = ~x.imag

function Base.real(x::Z4Group)
    if x.imag
        return 0
    else
        return x.minus ? -1 : 1
    end
end

function Base.imag(x::Z4Group)
    if ~x.imag
        return 0
    else
        return x.minus ? -1 : 1
    end
end

## Convert n in (0, 1, 2, 3) to binary as a Tuple of Bools
_int_to_bin(n::Integer) = (Bool((n & 2) >> 1), Bool(n & 1))

####
#### Algebra / mathematical operations
####

Base.:*(x::ConvertableNumber, y::Z4Group) = convert(Z4Group, x) * y
Base.:*(y::Z4Group, x::ConvertableNumber) = x * y

function Base.:*(x::Z4Group, y::Z4Group)
    minus = ~(~(x.minus ⊻ y.minus) ⊻ (x.imag & y.imag))
    imag = x.imag ⊻ y.imag
    return Z4Group(imag, minus)
end

Base.:-(x::Z4Group) = Z4Group(x.imag, ~x.minus)

end # module Z4Groups

# TODO: finish this
# function Base.:^(x::Z4Group, n::Integer)
#     nr = mod(n, 5)
#     if nr == 0
#         ip = Z4Group(1)
#     elseif nr == 1
#         ip = Z4Group(im)
#     elseif nr == 2
#         ip = Z4Group(-im)
# end


text:
"""
`module  JAC.LandeZeeman`
    ... a submodel of JAC that contains all methods for computing Lande factors and Zeeman properties for some level(s).
"""
module LandeZeeman

    using Printf, ..AngularMomentum, ..Basics, ..Defaults, ..InteractionStrength, ..ManyElectron, ..Nuclear, ..Radial, ..TableStrings

    """
    `struct  LandeZeeman.SublevelJ`  ... defines a type to specify a magnetic sublevel with well-defined J.

        + M                      ::AngularM64        ... M_J-value
        + energy                 ::Float64           ... energy of this sublevel
    """
    struct SublevelJ
        M                        ::AngularM64
        energy                   ::Float64
    end

    # `Base.show(io::IO, Jsublevel::LandeZeeman.SublevelJ)`  ... prepares a proper printout of the variable Jsublevel::LandeZeeman.SublevelJ.
    function Base.show(io::IO, Jsublevel::LandeZeeman.SublevelJ)
        println(io, "Sublevel [M=$(Jsublevel.M); energy = $(Jsublevel.energy)]")
    end

    # SublevelF defines a magnetic hyperfine sublevel
    """
    `struct  LandeZeeman.SublevelF`  ... defines a type to specify a magnetic hyperfine sublevel with well-defined F.

        + F                      ::AngularJ64        ... F-value
        + M                      ::AngularM64        ... M_F-value
        + energy                 ::Float64           ... energy of this sublevel
    """
    struct SublevelF
        F                        ::AngularJ64
        M                        ::AngularM64
        energy                   ::Float64
    end

    # `Base.show(io::IO, Fsublevel::LandeZeeman.SublevelF)`  ... prepares a proper printout of the variable Fsublevel::LandeZeeman.SublevelF.
    function Base.show(io::IO, Fsublevel::LandeZeeman.SublevelF)
        println(io, "Sublevel [F=$(Fsublevel.F); M=$(Fsublevel.M); energy = $(Fsublevel.energy)]")
    end

    """
    `struct  LandeZeeman.Outcome`
        ... defines a type to keep the Lande Factors and Zeeman spittling parameters of a fine-structure level.

        + Jlevel                 ::Level             ... Fine-structure levels to which the results refer to.
        + LandeJ                 ::Float64           ... Lande-J factor of this fine-structure level
        + amplitudeN1            ::Complex{Float64}  ... N1 amplitude
        + amplitudeDeltaN1       ::Complex{Float64}  ... Delta-N1 amplitude
        + nuclearI               ::AngularJ64        ... nuclear spin
        + hasJsublevels          ::Bool              ... True, if information about the magnetic sublevels are provided and false otherwise.
        + Jsublevels             ::Array{LandeZeeman.SublevelJ,1}   ... List of the magnetic fine-structure sublevels and data with
                                                                        well-defined J-value
        + hasFsublevels          ::Bool              ... True, if information about the magnetic hyperfine sublevels are provided and
                                                         false otherwise.
        + Fsublevels             ::Array{LandeZeeman.SublevelF,1}   ... List of the magnetic hyperfine sublevels and data with
                                                                        well-defined F-value
    """
    struct Outcome
        Jlevel                   ::Level
        LandeJ                   ::Float64
        amplitudeN1              ::Complex{Float64}
        amplitudeDeltaN1         ::Complex{Float64}
        nuclearI                 ::AngularJ64
        hasJsublevels            ::Bool
        Jsublevels               ::Array{LandeZeeman.SublevelJ,1}
        hasFsublevels            ::Bool
        Fsublevels               ::Array{LandeZeeman.SublevelF,1}
    end

    """
    `LandeZeeman.Outcome()`  ... constructor for an `empty` instance of LandeZeeman.Outcome.
    """
    function Outcome()
        Outcome(Level(), 0., 0., 0., AngularJ64(0), false, SublevelJ[],  false, SublevelF[])
    end

    # `Base.show(io::IO, outcome::LandeZeeman.Outcome)`  ... prepares a proper printout of the variable LandeZeeman.Outcome.
    function Base.show(io::IO, outcome::LandeZeeman.Outcome)
        println(io, "Jlevel:           $(outcome.Jlevel)  ")
        println(io, "LandeJ:           $(outcome.LandeJ)  ")
        println(io, "amplitudeN1:      $(outcome.amplitudeN1)  ")
        println(io, "amplitudeDeltaN1: $(outcome.amplitudeDeltaN1)  ")
        println(io, "nuclearI:         $(outcome.nuclearI)  ")
        println(io, "hasJsublevels:    $(outcome.hasJsublevels)  ")
        println(io, "Fsublevels:       $(outcome.Fsublevels)  ")
        println(io, "hasFsublevels:    $(outcome.hasFsublevels)  ")
        println(io, "Fsublevels:       $(outcome.Fsublevels)  ")
    end

    """
    `struct  LandeZeeman.Settings`
        ... defines a type for the details and parameters of computing the Lande Factors and Zeeman spittling
            of fine-structure levels.

        + calcLandeJ             ::Bool              ... True if Lande-J factors are to be calculated.
        + calcLandeF             ::Bool              ... True if Lande-F factors are to be calculated.
        + calcZeeman             ::Bool              ... True if Zeeman splittings are to be calculated.
        + includeSchwinger       ::Bool              ... True if Schwinger's QED correction ``\\Delta N^(1)`` is to be included,
                                                         and false otherwise.
        + BField                 ::Float64           ... Strength of the magnetic field in [Tesla]
        + printBefore            ::Bool              ... True if a list of selected levels is printed before the actual computations start.
        + levelSelection         ::LevelSelection    ... Specifies the selected levels, if any.
    """
    struct Settings
        calcLandeJ               ::Bool
        calcLandeF               ::Bool
        calcZeeman               ::Bool
        includeSchwinger         ::Bool
        BField                   ::Float64
        printBefore              ::Bool
        levelSelection           ::LevelSelection
     end

    """
    `LandeZeeman.Settings()`
        ... constructor for an `empty` instance of ZeemanSettings for the computation of isotope M and F parameters.
     """
     function Settings()
         Settings(false, false, false, false, 0., false, LevelSelection() )
    end

    # `Base.show(io::IO, settings::LandeZeeman.Settings)`  ... prepares a proper printout of the variable settings::LandeZeeman.Settings.
    function Base.show(io::IO, settings::LandeZeeman.Settings)
        println(io, "calcLandeJ:               $(settings.calcLandeJ)  ")
        println(io, "calcLandeF:               $(settings.calcLandeF)  ")
        println(io, "calcZeeman:               $(settings.calcZeeman)  ")
        println(io, "includeSchwinger:         $(settings.includeSchwinger)  ")
        println(io, "BField:                   $(settings.BField)  ")
        println(io, "printBefore:              $(settings.printBefore)  ")
        println(io, "levelSelection:           $(settings.levelSelection)  ")
    end

    """
    `LandeZeeman.amplitude(kind::String, rLevel::Level, sLevel::Level, grid::Radial.Grid)`
        ... to compute either the  N^(1) or Delta N^(1) Zeeman amplitude <alpha_r J_r || N^(1)) || alpha_s J_s>
            for a given pair of levels. A value::ComplexF64 is returned.
    """
    function  amplitude(kind::String, rLevel::Level, sLevel::Level, grid::Radial.Grid)
        nr = length(rLevel.basis.csfs);    ns = length(sLevel.basis.csfs);    matrix = zeros(ComplexF64, nr, ns)
        printstyled("Compute Zeeman $(kind[1:5]) matrix of dimension $nr x $ns ... ", color=:light_green)
        #
        if  rLevel.parity != sLevel.parity   return( ComplexF64(0.) )   end
        #
        for  r = 1:nr
            for  s = 1:ns
                me = 0.
                if  rLevel.basis.csfs[r].parity  != rLevel.parity    ||  sLevel.basis.csfs[s].parity  != sLevel.parity  ||
                    rLevel.parity != sLevel.parity    continue
                end
                #
                if      kind == "N^(1) amplitude"
                #--------------------------------
                    wa = compute("angular coefficients: 1-p, Grasp92", 0, 1, rLevel.basis.csfs[r], sLevel.basis.csfs[s])
                    for  coeff in wa
                        ja = Basics.subshell_2j(rLevel.basis.orbitals[coeff.a].subshell)
                        jb = Basics.subshell_2j(sLevel.basis.orbitals[coeff.b].subshell)
                        tamp  = InteractionStrength.zeeman_n1(rLevel.basis.orbitals[coeff.a], sLevel.basis.orbitals[coeff.b], grid)
                        me = me + coeff.T * tamp
                    end
                #
                elseif  kind == "Delta N^(1) amplitude"
                #--------------------------------------
                    wa = compute("angular coefficients: 1-p, Grasp92", 0, 1, rLevel.basis.csfs[r], sLevel.basis.csfs[s])
                    for  coeff in wa
                        ja = Basics.subshell_2j(rLevel.basis.orbitals[coeff.a].subshell)
                        jb = Basics.subshell_2j(sLevel.basis.orbitals[coeff.b].subshell)
                        tamp  = InteractionStrength.zeeman_Delta_n1(rLevel.basis.orbitals[coeff.a], sLevel.basis.orbitals[coeff.b], grid)
                        me = me + coeff.T * tamp
                    end
                #
                else    error("stop a")
                end
                #
                matrix[r,s] = me
            end
        end
        printstyled("done.\n", color=:light_green)
        amplitude = transpose(rLevel.mc) * matrix * sLevel.mc
        #
        return( amplitude )
    end

    """
    `LandeZeeman.computeAmplitudesProperties(outcome::LandeZeeman.Outcome, grid::Radial.Grid, settings::LandeZeeman.Settings)`
        ... to compute all amplitudes and properties of for a given level; an outcome::LandeZeeman.Outcome is returned for
            which the amplitudes and all requested properties are now evaluated explicitly.
    """
    function  computeAmplitudesProperties(outcome::LandeZeeman.Outcome, grid::Radial.Grid, settings::LandeZeeman.Settings)
        LandeJ = 0.0;   amplitudeN1 = 0.0;   amplitudeDeltaN1 = 0.0;    J = AngularMomentum.oneJ(outcome.Jlevel.J)
        amplitudeN1 = LandeZeeman.amplitude("N^(1) amplitude", outcome.Jlevel, outcome.Jlevel, grid)
        #
        if  settings.includeSchwinger
            amplitudeDeltaN1 = LandeZeeman.amplitude("Delta N^(1) amplitude", outcome.Jlevel, outcome.Jlevel, grid)
        end
        #
        if  settings.calcLandeJ    LandeJ = 2*(amplitudeN1 + amplitudeDeltaN1) / sqrt(J*(J+1))    end

        newOutcome = LandeZeeman.Outcome( outcome.Jlevel, LandeJ, amplitudeN1, amplitudeDeltaN1, outcome.nuclearI,
                                          outcome.hasJsublevels, outcome.Jsublevels, outcome.hasFsublevels, outcome.Fsublevels)
        return( newOutcome )
    end

    """
    `LandeZeeman.computeOutcomes(multiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid, settings::LandeZeeman.Settings; output=true)`
        ... to compute (as selected) the Lande J or F factors for the levels of the given multiplet and as specified by the given settings. The
            results are printed in neat tables to screen but nothing is returned otherwise.
    """
    function computeOutcomes(multiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid, settings::LandeZeeman.Settings; output=true)
        println("")
        printstyled("LandeZeeman.computeOutcomes(): The computation of the Zeeman amplitudes and Lande factors starts now ... \n", color=:light_green)
        printstyled("-------------------------------------------------------------------------------------------------------- \n", color=:light_green)
        println("")
        outcomes = LandeZeeman.determineOutcomes(multiplet, settings)
        # Display all selected levels before the computations start
        if  settings.printBefore    LandeZeeman.displayOutcomes(outcomes)    end
        # Calculate all amplitudes and requested properties
        newOutcomes = LandeZeeman.Outcome[]
        for  outcome in outcomes
            newOutcome = LandeZeeman.computeAmplitudesProperties(outcome, grid, settings)
            push!( newOutcomes, newOutcome)
        end
        # Print all results to screen
        LandeZeeman.displayResults(stdout, newOutcomes, nm, settings)
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary    LandeZeeman.displayResults(iostream, newOutcomes, nm, settings)   end
        #
        if    output    return( newOutcomes )
        else            return( nothing )
        end
    end

    """
    `LandeZeeman.determineOutcomes(multiplet::Multiplet, settings::LandeZeeman.Settings)`
        ... to determine a list of Outcomes's for the computation of Lande factors and/or Zeeman splittings for
            levels from the given multiplet. It takes into account the particular selections and settings.
            An Array{LandeZeeman.Outcome,1} is returned. Apart from the level specification, all physical properties are
            still set to zero during the initialization process.
    """
    function  determineOutcomes(multiplet::Multiplet, settings::LandeZeeman.Settings)
        # Define values that depend on the requested computations
        nuclearI = AngularJ64(0);    hasJsublevels = false;    Jsublevels = LandeZeeman.SublevelJ[]
                                     hasFsublevels = false;    Fsublevels = LandeZeeman.SublevelF[]

        outcomes = LandeZeeman.Outcome[]
        for  level  in  multiplet.levels
            if  Basics.selectLevel(level, settings.levelSelection)
                push!( outcomes, LandeZeeman.Outcome(level, 0., 0., 0., nuclearI, hasJsublevels, Jsublevels, hasFsublevels, Fsublevels) )
            end
        end
        return( outcomes )
    end

    """
    `LandeZeeman.displayOutcomes(outcomes::Array{LandeZeeman.Outcome,1})`
        ... to display a list of levels that have been selected for the computations. A small neat table of all selected
            levels and their energies is printed but nothing is returned otherwise.
    """
    function  displayOutcomes(outcomes::Array{LandeZeeman.Outcome,1})
        nx = 43
        println(" ")
        println("  Selected Lande-Zeeman levels:")
        println(" ")
        println("  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(10, "Level"; na=2);                             sb = sb * TableStrings.hBlank(12)
        sa = sa * TableStrings.center(10, "J^P";   na=4);                             sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(14, "Energy"; na=4);
        sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=4)
        println(sa);    println(sb);    println("  ", TableStrings.hLine(nx))
        #
        for  outcome in outcomes
            sa  = "  ";    sym = LevelSymmetry( outcome.Jlevel.J, outcome.Jlevel.parity)
            sa = sa * TableStrings.center(10, TableStrings.level(outcome.Jlevel.index); na=2)
            sa = sa * TableStrings.center(10, string(sym); na=4)
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", outcome.Jlevel.energy)) * "    "
            println( sa )
        end
        println("  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `LandeZeeman.displayResults(stream::IO, outcomes::Array{LandeZeeman.Outcome,1}, nm::Nuclear.Model, settings::LandeZeeman.Settings)`
        ... to display the energies, Lande factors, Zeeman amplitudes etc. for the selected levels. A neat table is printed but nothing is
            returned otherwise.
    """
    function  displayResults(stream::IO, outcomes::Array{LandeZeeman.Outcome,1}, nm::Nuclear.Model, settings::LandeZeeman.Settings)
        #
        if  settings.calcLandeJ
            nx = 135
            println(stream, " ")
            println(stream, "  Lande g_J factors and Zeeman amplitudes:")
            println(stream, " ")
            println(stream, "  ", TableStrings.hLine(nx))
            sa = "  ";   sb = "  "
            sa = sa * TableStrings.center(10, "Level"; na=2);                             sb = sb * TableStrings.hBlank(12)
            sa = sa * TableStrings.center(10, "J^P";   na=4);                             sb = sb * TableStrings.hBlank(14)
            sa = sa * TableStrings.center(14, "Energy"; na=5)
            sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=5)
            sa = sa * TableStrings.center(12, "Lande-J"; na=8);                           sb = sb * TableStrings.hBlank(20)
            sa = sa * TableStrings.center(62, "N1   -- Zeeman Amplitudes --   Delta N1"    ; na=4)
            sb = sb * TableStrings.center(62, "   Re              Im                 Re              Im"; na=5)
            println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
            #
            for  outcome in outcomes
                sa  = "  ";    sym = LevelSymmetry( outcome.Jlevel.J, outcome.Jlevel.parity)
                sa = sa * TableStrings.center(10, TableStrings.level(outcome.Jlevel.index); na=2)
                sa = sa * TableStrings.center(10, string(sym); na=4)
                energy = 1.0
                sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", energy))                    * "    "
                sa = sa * TableStrings.flushright(15, @sprintf("%.8e", outcome.LandeJ) )              * "    "
                sa = sa * TableStrings.flushright(15, @sprintf("%.8e", outcome.amplitudeN1.re) )      * " "
                sa = sa * TableStrings.flushright(15, @sprintf("%.8e", outcome.amplitudeN1.im) )      * "    "
                sa = sa * TableStrings.flushright(15, @sprintf("%.8e", outcome.amplitudeDeltaN1.re) ) * " "
                sa = sa * TableStrings.flushright(15, @sprintf("%.8e", outcome.amplitudeDeltaN1.im) ) * "    "
                println(stream, sa )
            end
            println(stream, "  ", TableStrings.hLine(nx))
        end
        #
        if  settings.calcLandeF
            nx = 80
            println(stream, " ")
            println(stream, "  Hyperfine levels and Lande g_F factors:")
            println(stream, " ")
            println(stream, "  ", TableStrings.hLine(nx))
            sa = "  ";   sb = "  "
            sa = sa * TableStrings.center(10, "Level"; na=2);                             sb = sb * TableStrings.hBlank(12)
            sa = sa * TableStrings.center(10, "J^P";   na=4);                             sb = sb * TableStrings.hBlank(14)
            sa = sa * TableStrings.center(14, "Energy"; na=5)
            sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=5)
            sa = sa * TableStrings.center(10, "F^P";   na=4);                             sb = sb * TableStrings.hBlank(14)
            sa = sa * TableStrings.center(12, "Lande-F"; na=8);                           sb = sb * TableStrings.hBlank(20)
            println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
            #
            for  outcome in outcomes
                sa  = "  ";    sym = LevelSymmetry( outcome.Jlevel.J, outcome.Jlevel.parity)
                sa = sa * TableStrings.center(10, TableStrings.level(outcome.Jlevel.index); na=2)
                sa = sa * TableStrings.center(10, string(sym); na=4)
                energy = 1.0
                sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", energy))                    * "    "
                println(stream, sa )
                #
                sc = TableStrings.hBlank( length(sa) + 1 )
                Flist = oplus(nm.spinI, outcome.Jlevel.J)
                for  F in Flist
                    symf   = LevelSymmetry( F, outcome.Jlevel.parity);      Fx = AngularMomentum.oneJ(F)
                    Jx     = AngularMomentum.oneJ(outcome.Jlevel.J);    Ix = AngularMomentum.oneJ(nm.spinI)
                    LandeF = (Fx*(Fx+1) + Jx*(Jx+1) - Ix*(Ix+1)) / (2*Fx*(Fx+1)) * outcome.LandeJ
                    sa = sc * TableStrings.center(10, string(symf); na=4)
                    sa = sa * TableStrings.flushright(15, @sprintf("%.8e", LandeF) )
                    println(stream, sa )
                end
            end
            println(stream, "  ", TableStrings.hLine(nx))
        end
        #
        if  settings.calcZeeman
            nx = 135
            println(stream, " ")
            println(stream, "  Zeeman splittings of (hyper-) fine-structure levels:")
            println(stream, " ")
            println(stream, "  ", TableStrings.hLine(nx))
            println(stream, "  ", TableStrings.hLine(nx))
        end
        #
        return( nothing )
    end

end # module


text:
export short_stringify!

# Original code courtesy of Nils Gudat
using ShortStrings: ShortString
using Missings: passmissing
using PooledArrays

# Functions to turn String columns into ShortStrings
function short_stringify(x::AbstractVector)
    y = ShortString("a"^maximum(length.(skipmissing(x))))
    return passmissing(typeof(y)).(x)
end

function short_stringify!(df::DataFrame)
    cols = unique([names(df, String); names(df, Union{Missing, String})])
    for c ∈ cols
        df[!, c] = PooledArray(short_stringify(df[!, c]))
    end
    return df
end

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


text:
# Lag, Spelade, Vunna, Oavgjorda, Förlorade, Gjorda Mål, Insläppta Mål, Måldifferens, Poäng, Poängsnitt
# T, GP, W, D, L, GF, GA, GD, PTS, PA
# Team, GamesPlayed, Wins, Draws, Losses, GoalsFor, GoalsAgainst, GoalDifference, Points, PointAverage

function outcome(homeScore,awayScore)
	if homeScore>awayScore
		:HomeVictory
	elseif homeScore<awayScore
		:AwayVictory
	else
		:Draw
	end
end

function sorttable!(table::DataFrame)
	sort!(table, [order(:Points, rev=true),
		          order(:GoalDifference, rev=true),
		          order(:GoalsFor, rev=true)])
end

# construct Home/Away tables by using incHome/incAway
function addgame!(T::DataFrame, games::DataFrame, g::Int; winPoints=3, incHome=true, incAway=true)
	teams = T[:Team]
	homeRow = findfirst(isequal(games[g,:HomeTeam]),teams)
	awayRow = findfirst(isequal(games[g,:AwayTeam]),teams)

	homeScore = games[g,:HomeScore]
	awayScore = games[g,:AwayScore]

	incHome && (T[homeRow,:GamesPlayed] += 1)
	incAway && (T[awayRow,:GamesPlayed] += 1)

	incHome && (T[homeRow,:GoalsFor] += homeScore)
	incAway && (T[awayRow,:GoalsFor] += awayScore)

	incHome && (T[homeRow,:GoalsAgainst] += awayScore)
	incAway && (T[awayRow,:GoalsAgainst] += homeScore)

	incHome && (T[homeRow,:GoalDifference] += homeScore-awayScore)
	incAway && (T[awayRow,:GoalDifference] += awayScore-homeScore)

	oc = outcome(homeScore,awayScore)

	# possibly correct outcome
	if haskey(games,:Outcome) && !ismissing(games[g,:Outcome])
		oc = Symbol(games[g,:Outcome])
	end

	incHome && (T[homeRow,:Wins] += oc==:HomeVictory)
	incAway && (T[awayRow,:Wins] += oc==:AwayVictory)

	incHome && (T[homeRow,:Draws] += oc==:Draw)
	incAway && (T[awayRow,:Draws] += oc==:Draw)

	incHome && (T[homeRow,:Losses] += oc==:AwayVictory)
	incAway && (T[awayRow,:Losses] += oc==:HomeVictory)

	homePoints = (oc==:HomeVictory)*winPoints + (oc==:Draw)
	awayPoints = (oc==:AwayVictory)*winPoints + (oc==:Draw)

	# possibly correct points
	if haskey(games,:PointsDeductedHome) && !ismissing(games[g,:PointsDeductedHome])
		homePoints -= games[g,:PointsDeductedHome]
	end
	if haskey(games,:PointsDeductedAway) && !ismissing(games[g,:PointsDeductedAway])
		awayPoints -= games[g,:PointsDeductedAway]
	end

	incHome && (T[homeRow,:Points] += homePoints)
	incAway && (T[awayRow,:Points] += awayPoints)
	T
end

function table(games::DataFrame; kwargs...)
	teams = sort(unique(vcat(games[:HomeTeam], games[:AwayTeam])))

	# unsorted table
	T = DataFrame(Team=teams)
	T[[:GamesPlayed, :Wins, :Draws, :Losses, :GoalsFor, :GoalsAgainst, :GoalDifference, :Points]] = 0

	# Compute everything for this game
	for g=1:size(games,1)
		addgame!(T, games, g; kwargs...)
	end

	T[:PointAverage] = T[:Points] ./ T[:GamesPlayed]
	sorttable!(T)
end

# assumes dest and table have exactly the same columns
function _addtable(dest, table)
	@assert names(dest)==names(table)

	for i=1:size(table,1)
		team = table[i,:Team]
		destRow = findfirst(isequal(team), dest[:Team])

		if destRow!=nothing # add to existing row
			for j=2:size(table,2)
				dest[destRow,j] += table[i,j]
			end
		else # create new row
			append!(dest, table[i,:])
		end
	end
end

# create marathon table - adds Seasons column
function mergetables(tables...)
	merged = DataFrame(Team=String[], GamesPlayed=Int[], Wins=Int[], Draws=Int[], Losses=Int[], GoalsFor=Int[], GoalsAgainst=Int[], GoalDifference=Int[], Points=Int[], Seasons=Int[])

	for table in tables
		table = copy(table) # copy before changing
		table[:Seasons] = 1 # all teams get one season from this table
		deletecols!(table, :PointAverage) # remove point average and recompute it later

		@assert names(merged) == names(table)

		_addtable(merged, table)
	end
	merged[:PointAverage] = merged[:Points] ./ merged[:GamesPlayed]
	sorttable!(merged)
end

defaultresultsfolder() = joinpath(dirname(pathof(Allsvenskan)),"..","data","results")

function allgamesbyseason(folder::String=defaultresultsfolder())
	files = readdir(folder)
	filter!(x->lowercase(splitext(x)[2])==".csv", files)
	# exit.(joinpath.([folder],files))
	# CSV.read.(joinpath.([folder],files), missingstring="NA")
	CSV.read.(joinpath.([folder],files), categorical=false,missingstring="NA",strings=:raw,rows_for_type_detect=10000)
end

allgames(G::AbstractArray{DataFrame}) = _vcat(G...)
allgames(folder::String=defaultresultsfolder()) = allgames(allgamesbyseason(folder))

alltables(G::AbstractArray{DataFrame}) = table.(G)
alltables(folder::String=defaultresultsfolder()) = alltables(allgamesbyseason(folder))

allteams(g::DataFrame) = unique(convert(Vector{String},vcat(g[:HomeTeam],g[:AwayTeam])))
allteams(folder::String=defaultresultsfolder()) = allteams(allgames(folder))

function filterbyteam!(G, allTeams, teams, cols=[:HomeTeam,:AwayTeam])
    if !all(x->x in allTeams, teams)
        i = findfirst(x->!(x in allTeams), teams)
        error("Unknown team: \"", teams[i], "\". Did you mean: \"",
              didyoumean(allTeams, teams[i]), "\"?")
    end

    for col in cols
        for i=1:size(G,1)
            mask = map!(x-> x in teams, BitVector(undef,size(G[i],1)), G[i][col])
            G[i] = G[i][mask,:]
        end
    end
    G
end

function marathon(; folder=defaultresultsfolder(), seasons=1:100000,
                  teams=nothing, turf=:both, kwargs...)
    @assert turf in [:both, :home, :away]

    G = allgamesbyseason(folder)
    allTeams = unique(vcat( map(x->x[:HomeTeam],G)..., map(x->x[:AwayTeam],G)... ))

    filter!(x-> Dates.year(Date(x[1,:Date])) in seasons, G)

    teams==nothing || filterbyteam!(G, allTeams, teams)

    incHome = turf==:both || turf==:home
    incAway = turf==:both || turf==:away

    mergetables(table.(G;incHome=incHome,incAway=incAway,kwargs...)...)
end


text:
## --- Load required packages
    using StatGeochem
    using Plots; gr()
    using Statistics, DelimitedFiles, SpecialFunctions

## --- # # # # # # # # # # # pMelts equil. batch melting # # # # # # # # # # # #
    # Note: melts_configure requires a working alphamelts installation.
    # If you don't have alphamelts, you'll need to download, install, and
    # configure it. The absolute path to the alphamelts perl script must be input here:
    meltspath = "/usr/local/bin/run_alphamelts.command"

    # The absolute or relative path to a direcory where you want to run MELTS -
    # A number of configuration and output files will be written here.
    scratchdir = "scratch/"

    # Conditions
    P_range = [20000,20000]
    T_range = [1700,800]
    # Starting composition
    elements = ["SiO2",  "TiO2","Al2O3","Fe2O3","Cr2O3",  "FeO",  "MnO",  "MgO",   "NiO",  "CoO",  "CaO",  "Na2O", "K2O", "P2O5", "H2O",]
    composition=[44.8030, 0.1991, 4.4305, 0.9778, 0.3823, 7.1350, 0.1344, 37.6345, 0.2489, 0.0129, 3.5345, 0.3584, 0.0289, 0.0209, 0.15,] #mcdbse (McDonough Pyrolite)
    # Run simulation
    melts_configure(meltspath, scratchdir, composition, elements, T_range, P_range,
        batchstring="1\nsc.melts\n10\n1\n3\n1\nliquid\n1\n1.0\n0\n10\n0\n4\n0\n",
        dT=-10, dP=0, index=1, version="pMELTS",mode="isobaric",fo2path="FMQ")

    # Read results
    melt_comp = melts_query_liquid(scratchdir, index=1)
    solid_comp = melts_query_solid(scratchdir, index=1)
    modes = melts_query_modes(scratchdir, index=1)

## --- Plot melt composition
    h = plot(xlabel="Percent melt",ylabel="Abudance (wt. %) in melt")
    for e in ["SiO2","Al2O3","CaO","MgO","FeO","Na2O","K2O"]
        plot!(h,melt_comp["mass"],melt_comp[e],label=e)
    end
    plot!(h,xlims=(0,100),framestyle=:box,fg_color_legend=:white,legend=:right)
    display(h)

## --- Plot solid composition
    h = plot(xlabel="Percent melt",ylabel="Abudance (wt. %) in solid")
    for e in ["SiO2","Al2O3","CaO","MgO","FeO","Na2O","K2O"]
        plot!(h,100-solid_comp["mass"],solid_comp[e],label=e)
    end
    plot!(h,xlims=(0,100),framestyle=:box,fg_color_legend=:white,legend=:right)
    display(h)

## --- Plot phase modes
    h = plot(xlabel="Temperature (C)",ylabel="Abudance (wt. %)")
    for m in modes["elements"][4:end]
        plot!(h,modes["Temperature"],modes[m],label=m)
    end
    plot!(h,ylims=(0,100),framestyle=:box,fg_color_legend=:white,legend=:topleft)
    display(h)

## --- End of File


text:
"""
    flat_repeat(x, n)

Shortcut for `Iterators.repeat([x], n)`
"""
flat_repeat(x, n) = Iterators.repeat([x], n)

function subscript_number(n::Int)
    sub_dict = Dict(
        '0' => '₀',
        '1' => '₁',
        '2' => '₂',
        '3' => '₃',
        '4' => '₄',
        '5' => '₅',
        '6' => '₆',
        '7' => '₇',
        '8' => '₈',
        '9' => '₉'
    )
    str_n = string(n)
    join([sub_dict[c] for c ∈ str_n])
end

truncate(arr::Vector, maxlen::Int) = arr[1:min(maxlen, length(arr))]
"""
    FakeVector{X, T}(creator::X, fieldname::Symbol, vec::Vector{T})

A utility structure which calls a function `setvectorindex!` when the `Base.setindex!` function is called on it. This
allows for some smart updating of expression trees when mutated.
"""
struct FakeVector{X, T}
    creator::X
    fieldname::Symbol
    vec::Vector{T}
end

function Base.getindex(fv::FakeVector{X, T}, i) where {X, T}
    throw(ErrorException("attempted to get index of `$(typeof(fv))`"))
end
function Base.setindex!(fv::FakeVector{X, T}, val, inds...) where {X, T}
    if length(inds) > 1
        @warn "indexes greater than length 1 (provided is length $(length(inds)) will be truncated to length 1"
    end
    setvectorindex!(fv.creator, fv.fieldname, val, first(inds))
end

struct Empty end
Base.show(::IO, ::Empty) = nothing