text:
include("precompile/precompile_Tullio.jl")
_precompile_()

module _Precompile_Core
    include("precompile/precompile_Core.jl")
end

module _Precompile_Base
    include("precompile/precompile_Base.jl")
end

if VERSION >= v"1.6-"
    _Precompile_Base._precompile_()
    _Precompile_Core._precompile_()
end

#=
# To generate these files, following:
# https://timholy.github.io/SnoopCompile.jl/stable/snoopi/
# For the full benefit, it seems you must first disable them here.

VERSION # 1.6

using Pkg
Pkg.activate(mktempdir())
Pkg.add("SnoopCompile")
using SnoopCompile

using Tullio

inf_timing = @snoopi begin
    Tullio._tullio(:( A[i] := (1:10)[i] ))
    Tullio._tullio(:( A[i+_] := (1:10)[i+j] ), :(i in 1:3))
    Tullio._tullio(:( A[i, J[k]] := B[i] * C[j,k] ), :(grad=Dual))
end

pc = SnoopCompile.parcel(inf_timing)

SnoopCompile.write(joinpath(pkgdir(Tullio), "src", "precompile"), pc)

=#

#=
# Simple test:

using Tullio
@time Tullio._tullio(:( A[i] := (1:10)[i+j] + (1:3)[j]) );

# Julia 1.5.2:
# 7.116702 seconds (17.27 M allocations: 873.755 MiB, 4.21% gc time)
# 5.373509 seconds (5.91 M allocations: 307.566 MiB, 1.39% gc time)

# Julia 1.6-
# 7.703672 seconds (19.13 M allocations: 1.073 GiB, 4.19% gc time)
# 5.538456 seconds (6.53 M allocations: 380.879 MiB, 1.88% gc time)

=#


text:
# # [How to port rules generated by SparseADRules to Zygote](@id power_method)
#
# In this demo we'll show how to insert SparseADRules's sparse gradient implementation to boost Zygote's gradient.
using Zygote
using SparseArrays, LinearAlgebra
using BenchmarkTools

# find max eigenvalue by power method
function power_max_eigen(A, x, target; niter=100)
    for i=1:niter
        x = A * x
        x /= norm(x)
    end
    return abs(x' * target)
end

# find min eigenvalue by shift power method
function power_min_eigen(A, x, target, λ; niter=100)
    A = λ*I - A
    for i=1:niter
        x = A * x
        x /= norm(x)
    end
    return abs(x' * target)
end

# generate test data
A = sprand(5000, 5000, 0.1)
x = randn(5000)
target = randn(5000)
λ = 2

# Zygote is able to generate correct gradient for sparse matrix and dense vectors
# in computation form of dense arrays
@btime max_ga_z, max_gx_z, max_gt_z = Zygote.gradient(power_max_eigen, $A, $x, $target) #12.954 s (4272 allocations: 27.91 GiB)
@btime min_ga_z, min_gx_z, min_gt_z = Zygote.gradient(power_min_eigen, $A, $x, $target, λ) # 12.343 s (4349 allocations: 28.00 GiB)
max_ga_z, max_gx_z, max_gt_z = Zygote.gradient(power_max_eigen, A, x, target)
min_ga_z, min_gx_z, min_gt_z = Zygote.gradient(power_min_eigen, A, x, target, λ)

using SparseADRules
# The gradient generated by SparseADRules is much faster
# since it keeps the orginal type in computation process
@btime max_ga, max_gx, max_gt = Zygote.gradient(power_max_eigen, $A, $x, $target) # 6.180 s (5072 allocations: 16.75 GiB)
@btime min_ga, min_gx, min_gt = Zygote.gradient(power_min_eigen, $A, $x, $target, λ) #  7.289 s (5149 allocations: 16.86 GiB)
max_ga, max_gx, max_gt = Zygote.gradient(power_max_eigen, A, x, target)
min_ga, min_gx, min_gt = Zygote.gradient(power_min_eigen, A, x, target, λ)

using Test

# check the results generated by Zygote and SparseADRules are consistent
@testset begin
    @testset "power method for max eigenvalue" begin
        @test max_ga_z ≈ max_ga
        @test max_gx_z ≈ max_gx
        @test max_gt_z ≈ max_gt
    end
    @testset "power method for min eigenvalue" begin
        @test min_ga_z ≈ min_ga
        @test min_gx_z ≈ min_gx
        @test min_gt_z ≈ min_gt
    end
end

ext: jl


text:
"""
    Either{L, R} = Union{Const{L}, Identity{R}}
    Either{Int, Bool}(true) == Identity(true)
    Either{Int, Bool}(1) == Const(1)
    Either{String}("left") == Const("left")
    Either{String}(:anythingelse) == Identity(:anythingelse)

A very common type to represent Result or Failure.
We reuse [`Identity`](@ref) as representing "Success" and [`Const`](@ref) for "Failure".
"""
const Either{L, R} = Union{Const{L}, Identity{R}}
Either{L, R}(x::L) where {L, R} = Const(x)
Either{L, R}(x::R) where {L, R} = Identity(x)
Either{L}(x::R) where {L, R} = Identity(x)
Either{L}(x::L) where {L} = Const(x)

Base.eltype(::Type{Either{L, R}}) where {L, R} = R
Base.eltype(::Type{Either{<:Any, R}}) where {R} = R
Base.eltype(::Type{Either}) = Any

# Helpers for Either
# ------------------

"""
    either(:left, bool_condition, "right")

If `bool_condition` is true, it returns the right value, wrapped into [Identity](@ref).
Else returns left side, wrapped into [Const](@ref).

Example
-------
```jldoctest
either(:left, 1 < 2, "right") == Identity("right")
```
"""
function either(left_false, comparison::Bool, right_true)
  comparison ? Identity(right_true) : Const(left_false)
end

"""
    @either true ? "right" : Symbol("left")
    @either if false
      "right"
    else
      :left
    end

Simple macro to reuse ? operator and simple if-else for constructing Either.
"""
macro either(expr::Expr)
  @assert expr.head == :if && length(expr.args) == 3 "@either macro only works on ? or simple if-else."
  if isa(expr.args[3], Expr) && (expr.args[3].head == :elseif)
    error("Found elseif, however can only deal with simple if-else.")
  end
  esc(quote
    DataTypesBasic.either($(expr.args[3]), $(expr.args[1]), $(expr.args[2]))
  end)
end

"""
    flip_left_right(Const(1)) == Identity(1)
    flip_left_right(Identity(:whatever)) == Const(:whatever)

exchanges left and right, i.e. what was Const, becomes an Identity and the other way around.
"""
flip_left_right(x::Const) = Identity(x.value)
flip_left_right(x::Identity) = Const(x.value)

"""
    iseither(::Const) = true
    iseither(::Identity) = true
    iseither(other) = false

check whether something is an [`Either`](@ref)
"""
iseither(::Const) = true
iseither(::Identity) = true
iseither(other) = false

"""
    isleft(::Const) = true
    isleft(::Identity) = false

Identical to [`isconst`](@ref), but might be easier to read when working with `Either`.
"""
isleft(::Const) = true
isleft(::Identity) = false

"""
    isright(::Const) = false
    isright(::Identity) = true

Identical to [`isidentity`](@ref), but might be easier to read when working with `Either`.
"""
isright(::Const) = false
isright(::Identity) = true

"""
    getleft(Const(:something)) == :something
    getleft(Identity(23))  # throws MethodError

Extract a value from a "left" `Const` value. Will result in loud error when used on anything else.
"""
getleft(e::Const) = e.value

"""
    getright(Identity(23)) == 23
    getright(Const(:something))  # throws MethodError

Extract a value from a "right" `Identity` value. Will result in loud error when used on anything else.
Identical to `Base.get` but explicit about the site (and not defined for other things)
"""
getright(e::Identity) = e.value

"""
    getleftOption(Identity(23)) == Option()
    getleftOption(Const(:something)) == Option(:something)

Convert to option, assuming you want to have the left value be preserved.
"""
getleftOption(e::Const) = Option(e.value)
getleftOption(e::Identity) = Option()

"""
    getrightOption(Identity(23)) == Option(23)
    getrightOption(Const(:something)) == Option()

Convert to option, assuming you want to have the right value be preserved.
Identical to `getOption`, just explicit about the site.
"""
getrightOption(e::Const) = Option()
getrightOption(e::Identity) = e

"""
    getOption(Identity(23)) == Option(23)
    getOption(Const(:something)) == Option()

Convert to option, assuming you want to have the right value be preserved, and the left value represented as
`Option()`.
"""
getOption(e::Const) = Option()
getOption(e::Identity) = e


text:
module Routing
using Pivot
export redirect_to

# TODO: make this DRY
function redirect_to(ctx, route)
  path = split(route, "/", keep=false)
  ep, params = fetch(ctx.data[:router], path)
  method = Pivot.STI[ctx.data[:request][:method]]
  ctx.data[:params] = params
  ctx |> Pivot.proper_method(method, ep.handlermap)
end

"""
# make_partial_route

build an incomplete `handle!` method call
"""
function make_partial_route(method, path, handler)
  :(handle!($(handler), $(method), $(path)))
end

"""
# get

"/" => HomeController.index
"""
macro get(phandler)
  :(make_partial_route(GET, $(phandler.args[1]), $(esc(phandler.args[2]))))
end

macro post(phandler)
  :(make_partial_route(POST, $(phandler.args[1]), $(esc(phandler.args[2]))))
end

macro put(phandler)
  :(make_partial_route(PUT, $(phandler.args[1]), $(esc(phandler.args[2]))))
end

macro delete(phandler)
  :(make_partial_route(DELETE, $(phandler.args[1]), $(esc(phandler.args[2]))))
end

macro resources(phandler)
  base_url = phandler.args[1]
  controller = phandler.args[2]
  if endswith(base_url, "/")
    base_url = base_url[1:end-1]
  end

  quote
    [
    @get    $(base_url)               => $(esc(Expr(:., controller, :(:index   ))))
    @get    $(base_url * "/new")      => $(esc(Expr(:., controller, :(:new     ))))
    @post   $(base_url)               => $(esc(Expr(:., controller, :(:create  ))))
    @get    $(base_url * "/:id")      => $(esc(Expr(:., controller, :(:show    ))))
    @get    $(base_url * "/:id/edit") => $(esc(Expr(:., controller, :(:edit    ))))
    @put    $(base_url * "/:id")      => $(esc(Expr(:., controller, :(:update  ))))
    @delete $(base_url * "/:id")      => $(esc(Expr(:., controller, :(:destroy ))))
    ]
  end
end

end


text:
"""
    NotImplementedError{M}(m)

`Exception` thrown when a method from a NetworkEpidemics interface
is not implemented by a given type. Borrowed from the `LightGraphs` source code.
"""
struct NotImplementedError{M} <: Exception
    m::M
    NotImplementedError(m::M) where {M} = new{M}(m)
end

Base.showerror(io::IO, ie::NotImplementedError) = print(io, "method $(ie.m) not implemented.")

_NI(m) = throw(NotImplementedError(m))


text:
macro juliadots(str::String)
	"Print Julia dots before bolded string output."
	"For vanilla strings."

	expr_out = :(
				print(Crayon(foreground=:red, bold=true), " ●");
				print(Crayon(foreground=:green, bold=true), "●");
				print(Crayon(foreground=:blue, bold=true), "● ");
				print(Crayon(foreground=:white, bold=true), $str);
				)

	return esc(expr_out)
end

macro juliadots(expr::Expr)
	"Print Julia dots before bolded string output."
	"For interpolated strings."

	expr_out = :(
				print(Crayon(foreground=:red, bold=true), " ●");
				print(Crayon(foreground=:green, bold=true), "●");
				print(Crayon(foreground=:blue, bold=true), "● ");
				print(Crayon(foreground=:white, bold=true), $expr)
				)

	return esc(expr_out)
end

macro boink(expr::Expr)
	"Add EPSILON to a numerical variable or array."

	expr_out = :($expr .+ EPSILON)
	return esc(expr_out)
end

macro positive(expr::Expr)
	"Add EPSILON to a numerical variable or array during variable assignment."

	if (expr.head == :.) || (expr.head == :ref)
		expr_out = :($expr .+= EPSILON)

	elseif expr.head == :(=)
		expr_out = :($(expr.args[1]) = EPSILON .+ $(expr.args[2]))
	end

	return esc(expr_out)
end

macro finite(expr::Expr)
	"Prevent overflow of floating point to Inf by returning floatmax() of value."

	if (expr.head == :.) || (expr.head == :ref)
		expr_out = :($expr = sign.($expr) .* min.(abs.($expr), floatmax.($expr)))

	elseif expr.head == :(=)
		expr_out = :($(expr.args[1]) = sign.($(expr.args[2])) .* min.(abs.($(expr.args[2])), floatmax.($(expr.args[2]))))

	elseif expr.head == :(-=)
		expr_out = :($(expr.args[1]) = sign.($(expr.args[1])) .* min.(abs.($(expr.args[1]) - $(expr.args[2])), floatmax.($(expr.args[1]) - $(expr.args[2]))))
	end

	return esc(expr_out)
end

macro buffer(expr::Expr)
	"Load individual variable into buffer memory."

	model = expr.args[1]

	if expr.args[2] == :(:alpha)
		expr_out = :($(esc(model)).alpha_buffer = cl.Buffer(Float32, $(esc(model)).context, (:rw, :copy), hostbuf=$(esc(model)).alpha))

	elseif expr.args[2] == :(:invsigma)
		expr_out = :($(esc(model)).invsigma_buffer = cl.Buffer(Float32, $(esc(model)).context, (:r, :copy), hostbuf=Matrix($(esc(model)).invsigma)))
	end

	return expr_out
end

macro host(expr::Expr)
	"Load individual variable into host memory."

	model = expr.args[1]

	if expr.args[2] == :(:Elogtheta_sum_buffer)
		expr_out = :($(esc(model)).Elogtheta_sum = cl.read($(esc(model)).queue, $(esc(model)).Elogtheta_sum_buffer))

	elseif expr.args[2] == :(:Elogtheta_dist_buffer)
		expr_out = :($(esc(model)).Elogtheta_dist = cl.read($(esc(model)).queue, $(esc(model)).Elogtheta_dist_buffer)[1:$(esc(model)).M])

	elseif expr.args[2] == :(:sigma_buffer)
		expr_out = :($(esc(model)).sigma = Symmetric(reshape(cl.read($(esc(model)).queue, $(esc(model)).sigma_buffer), $(esc(model)).K, $(esc(model)).K)))

	elseif expr.args[2] == :(:lambda_dist_buffer)
		expr_out = :($(esc(model)).lambda_dist = cl.read($(esc(model)).queue, $(esc(model)).lambda_dist_buffer)[1:$(esc(model)).M])

	elseif expr.args[2] == :(:gimel_buffer)
		expr_out =
		quote
		gimel_host = reshape(cl.read($(esc(model)).queue, $(esc(model)).gimel_buffer), $(esc(model)).K, $(esc(model)).M)
		$(esc(model)).gimel = [gimel_host[:,d] for d in 1:$(esc(model)).M]
		end
	end

	return expr_out
end

macro gpu(expr::Expr)
	"Train model on GPU."

	expr.args[1] == :train! || throw(ArgumentError("GPU acceleration only applies to the train! function."))

	quote
		local model = $(esc(expr.args[2]))
		local kwargs = [(kw.args[1], eval(kw.args[2])) for kw in $(esc(expr.args[3:end]))]

		if isa(model, LDA)
			gpumodel = gpuLDA(Corpus(), 1)
			gpumodel.corp = model.corp
			gpumodel.K = model.K
			gpumodel.M = model.M
			gpumodel.V = model.V
			gpumodel.N = model.N
			gpumodel.C = model.C
			gpumodel.topics = model.topics
			gpumodel.alpha = model.alpha
			gpumodel.beta = model.beta
			gpumodel.Elogtheta = model.Elogtheta
			gpumodel.gamma = model.gamma
			gpumodel.phi = [Array{Float32}(undef, 0, 0) for d in 1:model.M]

			for d in 1:model.M
				terms = model.corp[d].terms
				@positive gpumodel.phi[d] = model.beta_old[:,terms] .* exp.(model.Elogtheta_old[d])
				gpumodel.phi[d] ./= sum(gpumodel.phi[d], dims=1)
			end

			gpumodel.elbo = model.elbo

			train!(gpumodel; kwargs...)

			model.topics = gpumodel.topics
			model.alpha = gpumodel.alpha
			model.beta = gpumodel.beta
			model.Elogtheta = gpumodel.Elogtheta
			model.Elogtheta_old = deepcopy(model.Elogtheta)
			model.gamma = gpumodel.gamma
			model.phi = gpumodel.phi[1:min(gpumodel.M, 1)]
			model.elbo = gpumodel.elbo

			model.beta ./= sum(model.beta, dims=2)
			model.beta_old = copy(model.beta)
			model.phi ./= sum(model.phi, dims=1)
			nothing

		elseif isa(model, CTM)
			gpumodel = gpuCTM(Corpus(), 1)
			gpumodel.corp = model.corp
			gpumodel.K = model.K
			gpumodel.M = model.M
			gpumodel.V = model.V
			gpumodel.N = model.N
			gpumodel.C = model.C
			gpumodel.topics = model.topics
			gpumodel.mu = model.mu
			gpumodel.sigma = convert(Symmetric{Float32,Array{Float32,2}}, model.sigma)
			gpumodel.invsigma = convert(Symmetric{Float32,Array{Float32,2}}, model.invsigma)
			gpumodel.beta = model.beta
			gpumodel.lambda = model.lambda
			gpumodel.vsq = model.vsq
			gpumodel.logzeta = model.logzeta
			gpumodel.phi = [Array{Float32}(undef, 0, 0) for d in 1:model.M]

			for d in 1:model.M
				terms = model.corp[d].terms
				@positive gpumodel.phi[d] = additive_logistic(log.(model.beta_old[:,terms]) .+ model.lambda_old[d], dims=1)
				gpumodel.phi[d] ./= sum(gpumodel.phi[d], dims=1)
			end

			gpumodel.elbo = model.elbo

			train!(gpumodel; kwargs...)

			model.topics = gpumodel.topics
			model.mu = gpumodel.mu
			model.sigma = convert(Symmetric{Float64,Array{Float64,2}}, model.sigma)
			model.invsigma = convert(Symmetric{Float64,Array{Float64,2}}, model.invsigma)
			model.beta = gpumodel.beta
			model.lambda = gpumodel.lambda
			model.lambda_old = deepcopy(model.lambda)
			model.vsq = gpumodel.vsq
			model.logzeta = gpumodel.logzeta
			model.phi = gpumodel.phi[1:min(gpumodel.M, 1)]
			model.elbo = gpumodel.elbo

			model.beta ./= sum(model.beta, dims=2)
			model.beta_old = copy(model.beta)
			model.phi ./= sum(model.phi, dims=1)
			nothing

		elseif isa(model, CTPF)
			gpumodel = gpuCTPF(Corpus(), 1)
			gpumodel.corp = model.corp
			gpumodel.topics = model.topics
			gpumodel.scores = model.scores
			gpumodel.libs = model.libs
			gpumodel.drecs = model.drecs
			gpumodel.urecs = model.urecs
			gpumodel.K = model.K
			gpumodel.M = model.M
			gpumodel.V = model.V
			gpumodel.U = model.U
			gpumodel.N = model.N
			gpumodel.C = model.C
			gpumodel.R = model.R
			gpumodel.a = model.a
			gpumodel.b = model.b
			gpumodel.c = model.c
			gpumodel.d = model.d
			gpumodel.e = model.e
			gpumodel.f = model.f
			gpumodel.g = model.g
			gpumodel.h = model.h
			gpumodel.alef = model.alef
			gpumodel.he = model.he
			gpumodel.bet = model.bet
			gpumodel.vav = model.vav
			gpumodel.gimel = model.gimel
			gpumodel.zayin = model.zayin
			gpumodel.dalet = model.dalet
			gpumodel.het = model.het
			gpumodel.phi = [Array{Float32}(undef, 0, 0) for d in 1:model.M]
			gpumodel.xi = [Array{Float32}(undef, 0, 0) for d in 1:model.M]

			for d in 1:model.M
				terms = model.corp[d].terms
				readers = model.corp[d].readers

				gpumodel.phi[d] = exp.(digamma.(model.gimel_old[d]) - log.(model.dalet_old) - log.(model.bet_old) .+ digamma.(model.alef_old[:,terms]))
				gpumodel.phi[d] ./= sum(gpumodel.phi[d], dims=1)

				gpumodel.xi[d] = vcat(exp.(digamma.(model.gimel_old[d]) - log.(model.dalet_old) - log.(model.vav_old) .+ digamma.(model.he_old[:,readers])), exp.(digamma.(model.zayin_old[d]) - log.(model.het_old) - log.(model.vav_old) .+ digamma.(model.he_old[:,readers])))
				gpumodel.xi[d] ./= sum(gpumodel.xi[d], dims=1)
			end

			gpumodel.elbo = model.elbo

			train!(gpumodel; kwargs...)

			model.topics = gpumodel.topics
			model.scores = gpumodel.scores
			model.drecs = gpumodel.drecs
			model.urecs = gpumodel.urecs
			model.alef = gpumodel.alef
			model.alef_old = copy(model.alef)
			model.he = gpumodel.he
			model.he_old = copy(model.he)
			model.bet = gpumodel.bet
			model.bet_old = copy(model.bet)
			model.vav = gpumodel.vav
			model.vav_old = copy(model.vav)
			model.gimel = gpumodel.gimel
			model.gimel_old = deepcopy(model.gimel)
			model.zayin = gpumodel.zayin
			model.zayin_old = deepcopy(model.zayin)
			model.dalet = gpumodel.dalet
			model.dalet_old = copy(model.dalet)
			model.het = gpumodel.het
			model.het_old = copy(model.het)
			model.phi = gpumodel.phi[1:min(gpumodel.M, 1)]
			model.xi = gpumodel.xi[1:min(gpumodel.M, 1)]
			model.elbo = gpumodel.elbo

			model.phi ./= sum(model.phi, dims=1)
			model.xi ./= sum(model.xi, dims=1)
			nothing

		elseif isa(model, fLDA)
			nothing

		elseif isa(model, fCTM)
			nothing

		else
			train!(model; kwargs...)
		end
	end
end


text:
module EvaluateEquationDerivativeModule

using LinearAlgebra
import ..CoreModule: Node, Options
import ..UtilsModule: @return_on_false2, is_bad_array
import ..EquationUtilsModule: count_constants, index_constants, NodeIndex
import ..EvaluateEquationModule: deg0_eval

"""
    eval_diff_tree_array(tree::Node, cX::AbstractMatrix{T}, options::Options, direction::Int)

Compute the forward derivative of an expression, using a similar
structure and optimization to eval_tree_array. `direction` is the index of a particular
variable in the expression. e.g., `direction=1` would indicate derivative with
respect to `x1`.

# Returns

- `(evaluation, derivative, complete)::Tuple{AbstractVector{T}, AbstractVector{T}, Bool}`: the normal evaluation,
    the derivative, and whether the evaluation completed as normal (or encountered a nan or inf).
"""
function eval_diff_tree_array(
    tree::Node, cX::AbstractMatrix{T}, options::Options, direction::Int
)::Tuple{AbstractVector{T},AbstractVector{T},Bool} where {T<:Real}
    # TODO: Implement quick check for whether the variable is actually used
    # in this tree. Otherwise, return zero.
    evaluation, derivative, complete = _eval_diff_tree_array(tree, cX, options, direction)
    @return_on_false2 complete evaluation derivative
    return evaluation, derivative, !(is_bad_array(evaluation) || is_bad_array(derivative))
end

function _eval_diff_tree_array(
    tree::Node, cX::AbstractMatrix{T}, options::Options, direction::Int
)::Tuple{AbstractVector{T},AbstractVector{T},Bool} where {T<:Real}
    if tree.degree == 0
        diff_deg0_eval(tree, cX, options, direction)
    elseif tree.degree == 1
        diff_deg1_eval(tree, cX, Val(tree.op), options, direction)
    else
        diff_deg2_eval(tree, cX, Val(tree.op), options, direction)
    end
end

function diff_deg0_eval(
    tree::Node, cX::AbstractMatrix{T}, options::Options, direction::Int
)::Tuple{AbstractVector{T},AbstractVector{T},Bool} where {T<:Real}
    n = size(cX, 2)
    const_part = deg0_eval(tree, cX, options)[1]
    derivative_part = (tree.feature == direction) ? ones(T, n) : zeros(T, n)
    return (const_part, derivative_part, true)
end

function diff_deg1_eval(
    tree::Node, cX::AbstractMatrix{T}, ::Val{op_idx}, options::Options, direction::Int
)::Tuple{AbstractVector{T},AbstractVector{T},Bool} where {T<:Real,op_idx}
    n = size(cX, 2)
    (cumulator, dcumulator, complete) = eval_diff_tree_array(tree.l, cX, options, direction)
    @return_on_false2 complete cumulator dcumulator

    op = options.unaops[op_idx]
    diff_op = options.diff_unaops[op_idx]

    @inbounds @simd for j in 1:n
        x = op(cumulator[j])::T
        dx = diff_op(cumulator[j])::T * dcumulator[j]

        cumulator[j] = x
        dcumulator[j] = dx
    end
    return (cumulator, dcumulator, true)
end

function diff_deg2_eval(
    tree::Node, cX::AbstractMatrix{T}, ::Val{op_idx}, options::Options, direction::Int
)::Tuple{AbstractVector{T},AbstractVector{T},Bool} where {T<:Real,op_idx}
    n = size(cX, 2)
    (cumulator, dcumulator, complete) = eval_diff_tree_array(tree.l, cX, options, direction)
    @return_on_false2 complete cumulator dcumulator
    (array2, dcumulator2, complete2) = eval_diff_tree_array(tree.r, cX, options, direction)
    @return_on_false2 complete2 array2 dcumulator2

    op = options.binops[op_idx]
    diff_op = options.diff_binops[op_idx]

    @inbounds @simd for j in 1:n
        x = op(cumulator[j], array2[j])

        dx = dot(diff_op(cumulator[j], array2[j]), [dcumulator[j], dcumulator2[j]])

        cumulator[j] = x
        dcumulator[j] = dx
    end
    return (cumulator, dcumulator, true)
end

"""
    eval_grad_tree_array(tree::Node, cX::AbstractMatrix{T}, options::Options; variable::Bool=false)

Compute the forward-mode derivative of an expression, using a similar
structure and optimization to eval_tree_array. `variable` specifies whether
we should take derivatives with respect to features (i.e., cX), or with respect
to every constant in the expression.

# Returns

- `(evaluation, gradient, complete)::Tuple{AbstractVector{T}, AbstractMatrix{T}, Bool}`: the normal evaluation,
    the gradient, and whether the evaluation completed as normal (or encountered a nan or inf).
"""
function eval_grad_tree_array(
    tree::Node, cX::AbstractMatrix{T}, options::Options; variable::Bool=false
)::Tuple{AbstractVector{T},AbstractMatrix{T},Bool} where {T<:Real}
    n = size(cX, 2)
    if variable
        n_gradients = size(cX, 1)
    else
        n_gradients = count_constants(tree)
    end
    index_tree = index_constants(tree, 0)
    return eval_grad_tree_array(
        tree, n, n_gradients, index_tree, cX, options, Val(variable)
    )
end

function eval_grad_tree_array(
    tree::Node,
    n::Int,
    n_gradients::Int,
    index_tree::NodeIndex,
    cX::AbstractMatrix{T},
    options::Options,
    ::Val{variable},
)::Tuple{AbstractVector{T},AbstractMatrix{T},Bool} where {T<:Real,variable}
    evaluation, gradient, complete = _eval_grad_tree_array(
        tree, n, n_gradients, index_tree, cX, options, Val(variable)
    )
    @return_on_false2 complete evaluation gradient
    return evaluation, gradient, !(is_bad_array(evaluation) || is_bad_array(gradient))
end

function _eval_grad_tree_array(
    tree::Node,
    n::Int,
    n_gradients::Int,
    index_tree::NodeIndex,
    cX::AbstractMatrix{T},
    options::Options,
    ::Val{variable},
)::Tuple{AbstractVector{T},AbstractMatrix{T},Bool} where {T<:Real,variable}
    if tree.degree == 0
        grad_deg0_eval(tree, n, n_gradients, index_tree, cX, options, Val(variable))
    elseif tree.degree == 1
        grad_deg1_eval(
            tree, n, n_gradients, index_tree, cX, Val(tree.op), options, Val(variable)
        )
    else
        grad_deg2_eval(
            tree, n, n_gradients, index_tree, cX, Val(tree.op), options, Val(variable)
        )
    end
end

function grad_deg0_eval(
    tree::Node,
    n::Int,
    n_gradients::Int,
    index_tree::NodeIndex,
    cX::AbstractMatrix{T},
    options::Options,
    ::Val{variable},
)::Tuple{AbstractVector{T},AbstractMatrix{T},Bool} where {T<:Real,variable}
    const_part = deg0_eval(tree, cX, options)[1]

    if variable == tree.constant
        return (const_part, zeros(T, n_gradients, n), true)
    end

    index = variable ? tree.feature : index_tree.constant_index
    derivative_part = zeros(T, n_gradients, n)
    derivative_part[index, :] .= T(1)
    return (const_part, derivative_part, true)
end

function grad_deg1_eval(
    tree::Node,
    n::Int,
    n_gradients::Int,
    index_tree::NodeIndex,
    cX::AbstractMatrix{T},
    ::Val{op_idx},
    options::Options,
    ::Val{variable},
)::Tuple{AbstractVector{T},AbstractMatrix{T},Bool} where {T<:Real,op_idx,variable}
    (cumulator, dcumulator, complete) = eval_grad_tree_array(
        tree.l, n, n_gradients, index_tree.l, cX, options, Val(variable)
    )
    @return_on_false2 complete cumulator dcumulator

    op = options.unaops[op_idx]
    diff_op = options.diff_unaops[op_idx]

    @inbounds @simd for j in 1:n
        x = op(cumulator[j])::T
        dx = diff_op(cumulator[j])

        cumulator[j] = x
        for k in 1:n_gradients
            dcumulator[k, j] = dx * dcumulator[k, j]
        end
    end
    return (cumulator, dcumulator, true)
end

function grad_deg2_eval(
    tree::Node,
    n::Int,
    n_gradients::Int,
    index_tree::NodeIndex,
    cX::AbstractMatrix{T},
    ::Val{op_idx},
    options::Options,
    ::Val{variable},
)::Tuple{AbstractVector{T},AbstractMatrix{T},Bool} where {T<:Real,op_idx,variable}
    derivative_part = Array{T,2}(undef, n_gradients, n)
    (cumulator1, dcumulator1, complete) = eval_grad_tree_array(
        tree.l, n, n_gradients, index_tree.l, cX, options, Val(variable)
    )
    @return_on_false2 complete cumulator1 dcumulator1
    (cumulator2, dcumulator2, complete2) = eval_grad_tree_array(
        tree.r, n, n_gradients, index_tree.r, cX, options, Val(variable)
    )
    @return_on_false2 complete2 cumulator1 dcumulator1

    op = options.binops[op_idx]
    diff_op = options.diff_binops[op_idx]

    @inbounds @simd for j in 1:n
        x = op(cumulator1[j], cumulator2[j])
        dx = diff_op(cumulator1[j], cumulator2[j])
        cumulator1[j] = x
        for k in 1:n_gradients
            derivative_part[k, j] = dx[1] * dcumulator1[k, j] + dx[2] * dcumulator2[k, j]
        end
    end

    return (cumulator1, derivative_part, true)
end

end


text:
using Pkg
 if pwd() == pkgdir(RvSpectML) Pkg.activate("examples") end

verbose = true
  make_plots = false
  if verbose && !isdefined(Main,:RvSpectML)  println("# Loading RvSpecML")    end
  using RvSpectML
  if verbose && make_plots && !isdefined(Main,:RvSpectMLPlots)  println("# Loading RvSpecMLPlots")    end
  if make_plots    using RvSpectMLPlots    end
  if verbose   println("# Loading other packages")    end
  using Statistics

# Load data
all_spectra = include(joinpath(pkgdir(EchelleInstruments),"examples/read_expres_data_101501.jl"))
#all_spectra = include(joinpath(pkgdir(EchelleInstruments),"examples/read_neid_solar_data_20190918.jl"))
if typeof(get_inst(all_spectra)) <: AnyEXPRES   continuum_normalize_spectra!(all_spectra)   end
order_list_timeseries = extract_orders(all_spectra, pipeline_plan, orders_to_use=orders_to_use_default(first(all_spectra).inst), recalc=true )

# Preliminary CCF calculation
linelist_for_ccf_fn_w_path = joinpath(pkgdir(EchelleCCFs),"data","masks",linelist_for_ccf_filename)
line_list_df = prepare_line_list(linelist_for_ccf_fn_w_path, all_spectra, pipeline_plan,  orders_to_use=orders_to_use_default(first(all_spectra).inst), v_center_to_avoid_tellurics=ccf_mid_velocity, Δv_to_avoid_tellurics = RvSpectMLBase.max_bc,recalc=true)
(ccfs, v_grid) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan,  mask_type=:tophat, mask_scale_factor=2.0, ccf_mid_velocity=ccf_mid_velocity, recalc=true, calc_ccf_var=false)
rvs_ccf = calc_rvs_from_ccf_total(ccfs, pipeline_plan, v_grid=v_grid, times = order_list_timeseries.times, recalc=true)
if abs(ccf_mid_velocity-mean(rvs_ccf)) > 1e3     ccf_mid_velocity=mean(rvs_ccf)   end

# Measure how much space needed to ensure no telluric contamination from from one high-SNR spectra and recompute line list
ref_obs_idx = RvSpectMLBase.choose_obs_idx_for_init_guess(df_files_use,first(all_spectra).inst)
line_width_05 = RvSpectMLBase.calc_line_width(v_grid,view(ccfs,:,ref_obs_idx),frac_depth=0.05)
line_list_df = prepare_line_list(linelist_for_ccf_fn_w_path, all_spectra, pipeline_plan,  orders_to_use=orders_to_use_default(first(all_spectra).inst), v_center_to_avoid_tellurics=ccf_mid_velocity, Δv_to_avoid_tellurics=RvSpectMLBase.max_bc+2.5*line_width_50, recalc=true)

# Recompute CCF using wider mask (check results from expres_compare_mask_widths for mask suggestions)
((ccfs, ccf_vars), v_grid) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan,  mask_type=:halfcos, mask_scale_factor=10.0, ccf_mid_velocity=ccf_mid_velocity, calc_ccf_var=true, recalc=true )
alg_fit_rv = EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=1.0, measure_width_at_frac_depth=0.5, init_guess_ccf_σ=ccf_mid_velocity)
 rvs_ccf = calc_rvs_from_ccf_total(ccfs, ccf_vars, pipeline_plan, v_grid=v_grid, times = order_list_timeseries.times, alg_fit_rv=alg_fit_rv, recalc=true)

if !make_plots     make_no_plots!(pipeline_plan)    end
if make_plot(pipeline_plan, :ccf_total)
   using Plots
   using RvSpectMLPlots
   plt = RvSpectMLPlots.make_plot_ccf_vs_time(ccfs,v_grid = v_grid, output_path=output_dir, save_fig = save_plot(pipeline_plan,:ccf_total) )
   display(plt)
end

if make_plot(pipeline_plan, :ccf_total)
   plt = RvSpectMLPlots.make_heatmap_ccf_vs_time(ccfs,v_grid = v_grid, times=order_list_timeseries.times, output_path=output_dir, save_fig = save_plot(pipeline_plan,:ccf_total) )
   display(plt)
end

need_to!(pipeline_plan,:scalpels)
if need_to(pipeline_plan,:scalpels)
   rvs_scalpels = map(n->Scalpels.clean_rvs_scalpels(rvs_ccf, ccfs, num_basis=n), 1:5)
   println("RMS RVs cleaned by Scalpels: ",std.(rvs_scalpels) )
   dont_need_to!(pipeline_plan,:scalpels)
end

make_plot!(pipeline_plan,:scalpels)
if make_plot(pipeline_plan, :scalpels)
   @assert !need_to(pipeline_plan, :rvs_ccf_total)
   @assert !need_to(pipeline_plan, :ccf_total)
   plt = make_plots_scalpels(rvs_ccf, ccfs, max_num_basis=2, v_grid=v_grid, times=order_list_timeseries.times, output_path=output_dir, save_fig = save_plot(pipeline_plan, :scalpels) )
   display(plt)
end

# Continue, if want to compute and plot order CCFs (e.g., to figure out which to include/reject)
#(order_ccfs_novars, v_grid_order_ccfs_novars) = ccf_orders(order_list_timeseries, line_list_df, pipeline_plan,  ccf_mid_velocity=ccf_mid_velocity, recalc=true)
mask_shape = EchelleCCFs.TopHatCCFMask(order_list_timeseries.inst, scale_factor=2)
line_list = EchelleCCFs.BasicLineList2D(line_list_df)
v_step = 250
ccf_plan = EchelleCCFs.BasicCCFPlan(mask_shape = mask_shape, line_list=line_list, midpoint=ccf_mid_velocity, range_no_mask_change=RvSpectMLBase.max_bc, #= step=v_step,=#  max=RvSpectMLBase.max_bc  )
v_grid_orders_ccfs = EchelleCCFs.calc_ccf_v_grid(ccf_plan)
(order_ccfs, order_ccf_vars ) = EchelleCCFs.calc_order_ccf_and_var_chunklist_timeseries(order_list_timeseries, ccf_plan)
#all(order_ccfs.-order_ccfs_novars .== 0)

plt = plot()
 obsid = 1
 for chid in 1:size(order_ccfs,2)
   plot!(plt,v_grid_orders_ccfs,order_ccfs[:,chid,1],label=:none)
 end
 display(plt)

 plt = plot()
  obsid = 1
  for chid in 1:size(order_ccfs,2)
    plot!(plt,v_grid_orders_ccfs,order_ccfs[:,chid,1]./maximum(order_ccfs[:,chid,1]),label=:none)
  end
  display(plt)

plt = plot()
  chid = 1
  for obsid in 1:size(order_ccfs,3)
    plot!(plt,v_grid_orders_ccfs,(order_ccfs[:,chid,obsid])./maximum(order_ccfs[:,chid,obsid]),label=:none)
    if obsid ==1
    #  scatter!(plt,v_grid_orders_ccfs,(order_ccfs_novars[:,chid,obsid])./maximum(order_ccfs_novars[:,chid,obsid]),markersize=1.3,label=:none)
    end
  end
  display(plt)

plt = plot()
  chid = 10
  for obsid in 1:size(order_ccfs,3)
    plot!(plt,v_grid_orders_ccfs,(order_ccfs[:,chid,obsid])./maximum(order_ccfs[:,chid,obsid]) .-
                   mean(order_ccfs[:,chid,:]./maximum(order_ccfs[:,chid,:]),dims=2) ,label=:none)
  end
  begin
    obsid = 1
    scatter!(plt,v_grid_orders_ccfs,(order_ccfs[:,chid,obsid])./maximum(order_ccfs[:,chid,obsid]) .- mean(order_ccfs[:,chid,:]./maximum(order_ccfs[:,chid,:]),dims=2), markersize=1.3, label=:none)
  end
  display(plt)

if make_plot(pipeline_plan, :ccf_orders)
   plt = make_heatmap_ccf_vs_order(order_ccfs, v_grid=v_grid_orders_ccfs, order_labels=map(c->order_list_timeseries.chunk_list[1].data[c].λ.indices[2], 1:size(order_ccfs,2)),
                      output_path=output_dir * target_subdir, save_fig=save_plot(pipeline_plan,:ccf_orders)   )
   display(plt)
end


text:
struct GatesCuthillMcKeePermuter
    rev::Bool
    warnunconnected::Bool
    assumesymmetry::Bool
    checksymmetry::Bool
end

GatesCuthillMcKeePermuter(; rev=true, warnunconnected=false, assumesymmetry=false, checksymmetry=true) = GatesCuthillMcKeePermuter(rev, warnunconnected, assumesymmetry, checksymmetry)

function permute_stripe(A::Symmetric{Tv, SparseMatrixCSC{Tv, Ti}}, method::GatesCuthillMcKeePermuter; kwargs...) where {Tv, Ti}
    prm = CuthillMcKee.symrcm(SparseMatrixCSC(A), method.rev, false, method.warnunconnected)
    return DomainPermutation(prm)
end

function permute_stripe(A::SparseMatrixCSC, method::GatesCuthillMcKeePermuter; kwargs...)
    return permute_plaid(A, method; kwargs...)[end]
end

function permute_plaid(A::Symmetric{Tv, SparseMatrixCSC{Tv, Ti}}, method::GatesCuthillMcKeePermuter; kwargs...) where {Tv, Ti}
    prm = CuthillMcKee.symrcm(SparseMatrixCSC(A), method.rev, false, method.warnunconnected)
    return (DomainPermutation(prm), DomainPermutation(prm))
end

function permute_plaid(A::SparseMatrixCSC, method::GatesCuthillMcKeePermuter; adj_A = nothing, kwargs...)
    if method.assumesymmetry || (method.checksymmetry && issymmetric(pattern(A)))
        prm = CuthillMcKee.symrcm(A, method.rev, false, method.warnunconnected)
        return (DomainPermutation(prm), DomainPermutation(prm))
    else
        if adj_A === nothing
            adj_A = adjointpattern(A)
        end
        m, n = size(A)
        B = hvcat((2, 2), spzeros(m, m), A, adj_A, spzeros(n, n))
        prm = CuthillMcKee.symrcm(B, method.rev, false, method.warnunconnected)
        iprm = Int[]
        jprm = Int[]
        for v in prm
            if v > m
                push!(jprm, v - m)
            else
                push!(iprm, v)
            end
        end
        return (DomainPermutation(iprm), DomainPermutation(jprm))
    end
end


text:
using Attenuations
using Test
import Unitful: cm, g, keV

@testset "Elements - Table 1" begin
    @test Elements.Hydrogen.Z == 1
    @test μᵨ(Elements.Aluminum, 70keV) ≈ 0.2301 * cm^2 / g
    @test μ(Elements.Aluminum, 70keV) ≈ 0.6210399 * cm^-1
end

@testset "Compounds" begin
    @test μᵨ(Compound("H2O"), 70keV) ≈ 0.1929 * cm^2 / g
end

@testset "Mixtures" begin
    @test μᵨ(Mixture(Dict("H2O" => 0.9, "NaCl" => 0.1)), 70keV) ≈
          0.2015 * cm^2 / g
end

@testset "Materials - Table 2" begin
    @test μᵨ(Materials.water, 70keV) ≈ 0.1929 * cm^2 / g
    @test μ(Materials.water, 70keV) ≈ 0.1929 * cm^-1

    energies = 20keV:1keV:120keV
    @test length(μ(Elements.Gold, energies)) == length(energies)
end

@testset "Materials - composition sum to one" begin
    for m ∈ Materials
        @info m.name
        @test (m.composition |> values |> sum) ≈ 1 atol = 1e-4
    end
end

@testset "Values from data" begin
    @test typeof(val(μ(Elements.Gold, 70keV))) <: AbstractFloat
end


text:
#####
##### Chain
#####

# Note: Knet serialization requires everything to be mutable

mutable struct Chain
  layers
  Chain(xs...) = new(xs)
end

children(c::Chain) = c.layers

mapchildren(f, c::Chain) = Chain((f(l) for l in c.layers)...)

function (c::Chain)(x)
  for l in c.layers
    x = l(x)
  end
  return x
end

#####
##### Dense
#####

mutable struct Dense
  W
  b
  σ
end

function Dense(in::Int, out::Int, σ=identity)
  W = Knet.param(Float32, out, in, atype=Array)
  b = Knet.param0(Float32, out, atype=Array)
  return Dense(W, b, σ)
end

children(c::Dense) = (c.W, c.b, c.σ)

mapchildren(f, c::Dense) = Dense(f(c.W), f(c.b), f(c.σ))

(d::Dense)(x) = d.σ.(d.W * x .+ d.b)

#####
##### Conv
#####

mutable struct Conv
  W
  pad
end

function Conv(ksize::Tuple, nchans::Pair; pad=0)
  W = Knet.param(Float32, ksize..., nchans..., atype=Array)
  return Conv(W, pad)
end

(c::Conv)(x) = Knet.conv4(c.W, x, padding=c.pad)

children(c::Conv) = (c.W,)

mapchildren(f, c::Conv) = Conv(f(c.W), c.pad)

#####
##### BatchNorm
#####

mutable struct BatchNorm
  moments
  params
  train :: Bool
  activation
end

function BatchNorm(nchans::Int, activation=identity; momentum)
  moments = Knet.bnmoments(momentum=momentum)
  params = Knet.bnparams(Float32, nchans) |> Knet.param
  BatchNorm(moments, params, true, activation)
end

function (l::BatchNorm)(x)
  y = Knet.batchnorm(x, l.moments, l.params, training=l.train)
  return l.activation.(y)
end

#####
##### Skip connection
#####

mutable struct SkipConnection
  block
  connection
end

children(c::SkipConnection) = (c.block,)

function mapchildren(f, c::SkipConnection)
  SkipConnection(f(c.block), c.connection)
end

function (skip::SkipConnection)(input)
  skip.connection(skip.block(input), input)
end

#####
##### Uilities
#####

using Knet: relu

flatten(x) = reshape(x, :, size(x)[end])

softmax(x) = Knet.softmax(x, dims=1)


text:
@testset "Replace" begin
    @testset "Default" begin
        # Tester that replaces with 0.0
        tester = ImputorTester(Replace; values=0.0)

        # Defining our own equality because an empty constructor isn't supported
        @testset "Equality" begin
            @test tester.imp(; tester.kwargs...) == tester.imp(; tester.kwargs...)
        end

        test_vector(tester)
        test_matrix(tester)
        test_dataframe(tester)
        # groupby test also fail because it tries to call an empty constructor
        # test_groupby(tester)
        test_axisarray(tester)
        test_nameddimsarray(tester)
        test_keyedarray(tester)
        test_columntable(tester)
        test_rowtable(tester)
    end

    @testset "Multiple values over tables" begin
        imp = Replace(; values=(DateTime(0), -9999, NaN, ""))
        df_table = DataFrame(
            :time => [missing, DateTime(2020, 02, 02), DateTime(2121, 12, 12)],
            :loc => [12, -5, missing],
            :val => [1.5, missing, 3.0],
            :desc => ["foo", "bar", missing],
        )
        df_expected = DataFrame(
            :time => [DateTime(0), DateTime(2020, 02, 02), DateTime(2121, 12, 12)],
            :loc => [12, -5, -9999],
            :val => [1.5, NaN, 3.0],
            :desc => ["foo", "bar", ""],
        )

        @testset "DataFrame" begin
            table = copy(df_table)
            expected = copy(df_expected)

            result = impute(table, imp)
            @test isequal(result, expected)

            table2 = deepcopy(table)
            impute!(table2, imp)
            @test isequal(table2, expected)
        end

        @testset "Column Table" begin
            table = Tables.columntable(df_table)
            expected = Tables.columntable(df_expected)

            result = impute(table, imp)
            @test isequal(result, expected)

            table2 = deepcopy(table)
            impute!(table2, imp)
            @test isequal(table2, expected)
        end

        @testset "Row Table" begin
            table = Tables.rowtable(df_table)
            expected = Tables.rowtable(df_expected)

            result = impute(table, imp)
            @test isequal(result, expected)

            table2 = deepcopy(table)
            impute!(table2, imp)
            @test !isequal(table2, expected)
            @test isequal(table2, table)
        end
    end
end


text:
import Base.run
import Base.read
using HDF5

"""
A named collection of jobs, a hdf5 source/sink specification, and properties
which can apply to any job
"""
type H5Task
    file
    source
    sink
    jobs
    props
end

"""
A job which can be run parallel to other jobs with source arguments stored in
the hdf5 file
"""
type H5Job
    task
    script
    func
    args
    props
    id
end

"""
Create a job which consists of running 'func'(map(read,'args')...) after loading
'script'.
This job in placed in 'task's collection of jobs
"""
function H5Job(task, script, func, args)
    hh = H5Job(task, script, func, args, Dict{Any,Any}(), length(task.jobs)+1)
    push!(task.jobs, hh)
    k = [keys(task.props)...]
    for j=1:length(k)
        hh[k[j]] = task.props[k[j]]
    end
    hh
end

"""
Provide an argument which can be accessed by any job
"""
function setindex!(task::H5Task, val, name::ASCIIString)
    task.props[name] = val
end

"""
Provide an argument for the specific job.
NOTE: this argument overrides any specified for the task
"""
function setindex!(job::H5Job, val, name::ASCIIString)
    job.props[name] = val
end

"""
Read the results of a particular job
"""
function read(job::H5Job)
    h5open(job.task.file) do file
        read(file[string("/results/",job.id,"/result")])
    end
end

"""
Run the task via N different child processes
"""
function run(task::H5Task, N::Int)
    gc()
    taskfile = task.file
    lockfile = string(taskfile, ".lock")
    println("touch $lockfile")
    run(`touch $lockfile`)

    file = h5open(task.file, false, true, true, false, false)
    if(exists(file, task.source))
        o_delete(file[task.source])
    end
    if(exists(file, task.sink))
        o_delete(file[task.sink])
    end
    for i=1:length(task.jobs)
        k = [keys(task.jobs[i].props)...]
        for j=1:length(k)
            file[string(task.source,i,"/",k[j])] = task.jobs[i].props[k[j]]
        end
    end

    if(exists(file, "/sched/"))
        o_delete(file["/sched/"])
    end
    for i=1:length(task.jobs)
        file[string("/sched/",i,"/taken")]  = 0
        file[string("/sched/",i,"/script")] = task.jobs[i].script
        file[string("/sched/",i,"/func")]   = task.jobs[i].func
        file[string("/sched/",i,"/args")]   = task.jobs[i].args
    end

    close(file)

    running = Any[]

    println("Starting Jobs...")
    for i=1:N
        #println(`julia -e \"reload(\"child.jl\");main(\"$taskfile\",\"$lockfile\")\"`)
        #println(readall(`julia -e reload(\"child.jl\");main(\"$taskfile\",\"$lockfile\")`))
        push!(running, spawn(`julia --depwarn=yes -e
        "using LocalJobs;child_main(\"$taskfile\",\"$lockfile\")"`,
        (open("/dev/null","r"), open("stdout.$i.log","w"), open("stderr.$i.log","w"))))
    end

    println("Merging Jobs...")
    for i=1:N
        wait(running[i])
    end

    println("Job Output...")
    for i=1:N
        out = readall("stdout.$i.log")
        err = readall("stderr.$i.log")
        if(length(err) != 0)
            println("Worker $i Errors:")
            println(err)
        end
        if(length(out) != 0)
            println("Worker $i Output:")
            println(out)
        end
    end

    println("Cleanup Temporary Job Files...")
    #for i=1:N
    #    rm("stdout.$i.log")
    #    rm("stderr.$i.log")
    #end
    rm(lockfile)

    println("Done...")

end


text:
using ProgressMeter, Plots
import KitBase

begin
    # case
    matter = "gas"
    case = "relaxation"
    space = "1d3f1v"
    flux = "kfvs"
    collision = "rykov"
    nSpecies = 1
    interpOrder = 1
    limiter = "vanleer"
    boundary = "fix"
    cfl = 0.8
    maxTime = 5.0

    # phase space
    x0 = -1.0
    x1 = 1.0
    nx = 1
    pMeshType = "uniform"
    nxg = 0

    # velocity space
    umin = -5.0
    umax = 5.0
    nu = 32
    vMeshType = "rectangle"
    nug = 0

    # gas
    knudsen = 0.1
    mach = 0.0
    prandtl = 0.72
    inK = 2.0
    inKr = 2.0
    omega = 0.81
    alphaRef = 1.0
    omegaRef = 0.5
    Tr0 = 91.5 / 273
    Z0 = 18.1
    sigma = 1 / 1.55
    omega1 = 0.2354
    omega2 = 0.3049
end

γ = KitBase.heat_capacity_ratio(inK, inKr, 1)
set = KitBase.Setup(
    matter,
    case,
    space,
    flux,
    collision,
    nSpecies,
    interpOrder,
    limiter,
    boundary,
    cfl,
    maxTime,
)
pSpace = KitBase.PSpace1D(x0, x1, nx)
vSpace = KitBase.VSpace1D(umin, umax, nu)
gas = KitBase.DiatomicGas(
    knudsen,
    mach,
    prandtl,
    inK,
    inKr,
    γ,
    omega,
    alphaRef,
    omegaRef,
    Tr0,
    Z0,
    sigma,
    omega1,
    omega2,
)

prim0 = [1.0, 0.0, 1.6556, 1.0, 100.0]
w0 = KitBase.prim_conserve(prim0, γ, inKr)

h0 = zeros(vSpace.nu)
b0 = similar(h0)
r0 = similar(h0)

Ht = similar(h0)
Bt = similar(b0)
Rt = similar(r0)
Hr = similar(h0)
Br = similar(b0)
Rr = similar(r0)

KitBase.maxwellian!(Ht, Bt, Rt, Hr, Br, Rr, vSpace.u, prim0, inK, inKr)

Zr = KitBase.rykov_zr(1 / prim0[4], Tr0, Z0)

@. h0 = (1.0 - 1.0 / Zr) * Ht + 1.0 / Zr * Hr
@. b0 = (1.0 - 1.0 / Zr) * Bt + 1.0 / Zr * Br
@. r0 = (1.0 - 1.0 / Zr) * Rt + 1.0 / Zr * Rr

w0 .= KitBase.diatomic_moments_conserve(h0, b0, r0, vSpace.u, vSpace.weights)
prim0 .= KitBase.conserve_prim(w0, inK, inKr)

ib = KitBase.IB3F(
    w0,
    prim0,
    h0,
    b0,
    r0,
    prim0,
    nothing,
    nothing,
    nothing,
    w0,
    prim0,
    h0,
    b0,
    r0,
    prim0,
    nothing,
    nothing,
    nothing,
)

cd(@__DIR__)
folder = pwd()

ks = KitBase.SolverSet(set, pSpace, vSpace, gas, ib, folder)

w = deepcopy(w0)
prim = deepcopy(prim0)
h = deepcopy(h0)
b = deepcopy(b0)
r = deepcopy(r0)

function step(KS, w, prim, h, b, r, dt)

    w_old = deepcopy(w)
    q = KitBase.heat_flux(h, b, r, prim, KS.vSpace.u, KS.vSpace.weights)

    MHT = similar(h)
    MBT = similar(b)
    MRT = similar(r)
    MHR = similar(h)
    MBR = similar(b)
    MRR = similar(r)

    KitBase.maxwellian!(
        MHT,
        MBT,
        MRT,
        MHR,
        MBR,
        MRR,
        KS.vSpace.u,
        prim,
        KS.gas.K,
        KS.gas.Kr,
    )
    τ_old = KitBase.vhs_collision_time(prim[1:end-1], KS.gas.μᵣ, KS.gas.ω)
    Zr = KitBase.rykov_zr(1.0 / prim[4], KS.gas.T₀, KS.gas.Z₀)
    Er0_old = 0.5 * sum(@. KS.vSpace.weights * ((1.0 / Zr) * MRR + (1.0 - 1.0 / Zr) * MRT))

    w[4] += dt * (Er0_old - w_old[4]) / τ_old
    prim .= KitBase.conserve_prim(w, KS.gas.K, KS.gas.Kr)

    #--- calculate M^{n+1} and tau^{n+1} ---#
    KitBase.maxwellian!(
        MHT,
        MBT,
        MRT,
        MHR,
        MBR,
        MRR,
        KS.vSpace.u,
        prim,
        KS.gas.K,
        KS.gas.Kr,
    )

    SHT = similar(h)
    SBT = similar(b)
    SRT = similar(r)
    SHR = similar(h)
    SBR = similar(b)
    SRR = similar(r)
    KitBase.rykov!(
        SHT,
        SBT,
        SRT,
        SHR,
        SBR,
        SRR,
        KS.vSpace.u,
        MHT,
        MBT,
        MRT,
        MHR,
        MBR,
        MRR,
        q,
        prim,
        KS.gas.Pr,
        KS.gas.K,
        KS.gas.σ,
        KS.gas.ω₁,
        KS.gas.ω₂,
    )

    MH = (1.0 - 1.0 / Zr) * (MHT + SHT) + 1.0 / Zr * (MHR + SHR)
    MB = (1.0 - 1.0 / Zr) * (MBT + SBT) + 1.0 / Zr * (MBR + SBR)
    MR = (1.0 - 1.0 / Zr) * (MRT + SRT) + 1.0 / Zr * (MRR + SRR)

    τ = KitBase.vhs_collision_time(prim[1:end-1], KS.gas.μᵣ, KS.gas.ω)

    #--- update distribution function ---#
    for i in eachindex(h)
        h[i] = (h[i] + dt / τ * MH[i]) / (1.0 + dt / τ)
        b[i] = (b[i] + dt / τ * MB[i]) / (1.0 + dt / τ)
        r[i] = (r[i] + dt / τ * MR[i]) / (1.0 + dt / τ)
    end

end

dt = 0.001
nt = 100

function solve(KS, w, prim, h, b, r, dt, nt)
    whis = zeros(4, nt)

    @showprogress for iter = 1:nt
        step(KS, w, prim, h, b, r, dt)
        whis[:, iter] .= w
    end

    return whis
end

w_his = solve(ks, w, prim, h, b, r, dt, nt)
prim_his = zeros(5, nt)
for i = 1:nt
    prim_his[:, i] .= KitBase.conserve_prim(w_his[:, i], ks.gas.K, ks.gas.Kr)
end

plot(1 ./ prim_his[4, :])
plot!(1 ./ prim_his[5, :])


text:
# epsxc is always of type Array{Float64,1}
# Vxc is always of type Array{Float64,2}

# Rhoe can be spinpol or not
function calc_epsxc_VWN( xc_calc::XCCalculator, Rhoe )
    Npoints = size(Rhoe,1)
    epsxc = zeros(Float64,Npoints)
    calc_epsxc_VWN!( xc_calc, Rhoe, epsxc )
    return epsxc
end

function calc_Vxc_VWN( xc_calc::XCCalculator, Rhoe )
    Npoints = size(Rhoe,1)
    Nspin = size(Rhoe,2)
    #Vxc = zeros(Float64, Npoints, Nspin)
    Vxc = similar(Rhoe)
    if Nspin == 1
        calc_Vxc_VWN!( xc_calc, Rhoe[:,1], Vxc )
        return Vxc
    else
        calc_Vxc_VWN!( xc_calc, Rhoe, Vxc )
        return Vxc
    end

end

#
# Inplace version
#

# Spin unpolarized version
function calc_epsxc_VWN!( xc_calc::XCCalculator, Rhoe::Array{Float64,1}, epsxc )
    Npoints = size(Rhoe,1)
    for ip in 1:Npoints
        ss_x, _ = XC_x_slater( Rhoe[ip] )
        ss_c, _ = XC_c_vwn( Rhoe[ip] )
        epsxc[ip] = ss_x + ss_c
    end
    return
end

# Spin unpolarized version
function calc_Vxc_VWN!( xc_calc::XCCalculator, Rhoe::Array{Float64,1}, Vxc )
    Npoints = size(Rhoe,1)
    for ip in 1:Npoints
        _, vx = XC_x_slater( Rhoe[ip] )
        _, vc = XC_c_vwn( Rhoe[ip] )
        Vxc[ip] = vx + vc
    end
    return
end

# Spin-polarized version
function calc_epsxc_VWN!( xc_calc::XCCalculator, Rhoe::Array{Float64,2}, epsxc::Array{Float64,1} )
    Npoints = size(Rhoe,1)
    Nspin = size(Rhoe,2)
    if Nspin == 1
        calc_epsxc_VWN!( xc_calc, Rhoe[:,1], epsxc )
        return
    end

    for ip in 1:Npoints
        ρ = Rhoe[ip,1] + Rhoe[ip,2]
        ζ = (Rhoe[ip,1] - Rhoe[ip,2])/ρ

        ss_x, _, _ = XC_x_slater_spin( ρ, ζ )
        ss_c, _, _ = XC_c_vwn_spin( ρ, ζ )

        epsxc[ip] = ss_x + ss_c
    end
    return
end

function calc_Vxc_VWN!( xc_calc::XCCalculator, Rhoe::Array{Float64,2}, Vxc::Array{Float64,2} )
    Npoints = size(Rhoe,1)
    Nspin = size(Rhoe,2)
    if Nspin == 1
        calc_Vxc_VWN!( xc_calc, Rhoe[:,1], Vxc )
        return
    end
    for ip in 1:Npoints
        ρ = Rhoe[ip,1] + Rhoe[ip,2]
        ζ = (Rhoe[ip,1] - Rhoe[ip,2])/ρ

        _, vxup, vxdn = XC_x_slater_spin( ρ, ζ )
        _, vcup, vcdn = XC_c_vwn_spin( ρ, ζ )

        Vxc[ip,1] = vxup + vcup
        Vxc[ip,2] = vxdn + vcdn
    end
    return
end


text:
function create_splittable_from_params(params::cluster_parameters, α::Float32)
    params_l = deepcopy(params)
    params_l.distribution = sample_distribution(params.posterior_hyperparams)
    params_r = deepcopy(params)
    params_r.distribution = sample_distribution(params.posterior_hyperparams)
    lr_weights = rand(Dirichlet(Float64.([α / 2, α / 2])))
    return splittable_cluster_params(params,params_l,params_r,lr_weights, false, [-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf])
end

function merge_clusters_to_splittable(cpl::cluster_parameters,cpr::cluster_parameters, α::Float32)
    suff_stats = aggregate_suff_stats(cpl.suff_statistics,cpr.suff_statistics)
    posterior_hyperparams = calc_posterior(cpl.hyperparams, suff_stats)
    lr_weights = rand(Dirichlet(Float64.([cpl.suff_statistics.N + (α / 2), cpr.suff_statistics.N + (α / 2)])))
    cp = cluster_parameters(cpl.hyperparams, cpl.distribution, suff_stats, posterior_hyperparams)
    return splittable_cluster_params(cp,cpl,cpr,lr_weights, false, [-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf])
end

function should_merge!(should_merge::AbstractArray{Float32,1}, cpl::cluster_parameters,cpr::cluster_parameters, α::Float32, final::Bool)
    new_suff = aggregate_suff_stats(cpl.suff_statistics, cpr.suff_statistics)
    cp = cluster_parameters(cpl.hyperparams, cpl.distribution, new_suff,cpl.posterior_hyperparams)
    cp.posterior_hyperparams = calc_posterior(cp.hyperparams, cp.suff_statistics)
    log_likihood_l = log_marginal_likelihood(cpl.hyperparams, cpl.posterior_hyperparams, cpl.suff_statistics)
    log_likihood_r = log_marginal_likelihood(cpr.hyperparams, cpr.posterior_hyperparams, cpr.suff_statistics)
    log_likihood = log_marginal_likelihood(cp.hyperparams, cp.posterior_hyperparams, cp.suff_statistics)
    log_HR = (-log(α) + lgamma(α) -2*lgamma(0.5*α) + lgamma(cp.suff_statistics.N) -lgamma(cp.suff_statistics.N + α) +
        lgamma(cpl.suff_statistics.N + 0.5*α)-lgamma(cpl.suff_statistics.N)  - lgamma(cpr.suff_statistics.N) +
        lgamma(cpr.suff_statistics.N + 0.5*α)+ log_likihood- log_likihood_l- log_likihood_r)
    # log_HR = -(log(α) +
    #     lgamma(cpl.suff_statistics.N) + log_likihood_l +
    #     lgamma(cpr.suff_statistics.N) + log_likihood_r -
    #     (lgamma(cp.suff_statistics.N) + log_likihood))
    if (log_HR > log(rand())) || (final && log_HR > log(0.1))
        should_merge .= 1
    end
end

function sample_cluster_params(params::splittable_cluster_params, α::Float32, first::Bool)
    points_count = Vector{Float32}()
    params.cluster_params.distribution = sample_distribution(first ? params.cluster_params.hyperparams : params.cluster_params.posterior_hyperparams)
    params.cluster_params_l.distribution = sample_distribution(first ? params.cluster_params_l.hyperparams : params.cluster_params_l.posterior_hyperparams)
    params.cluster_params_r.distribution = sample_distribution(first ? params.cluster_params_r.hyperparams : params.cluster_params_r.posterior_hyperparams)
    push!(points_count, params.cluster_params_l.suff_statistics.N)
    push!(points_count, params.cluster_params_r.suff_statistics.N)
    points_count .+= α / 2
    params.lr_weights = rand(Dirichlet(Float64.(points_count)))

    log_likihood_l = log_marginal_likelihood(params.cluster_params_l.hyperparams,params.cluster_params_l.posterior_hyperparams, params.cluster_params_l.suff_statistics)
    log_likihood_r = log_marginal_likelihood(params.cluster_params_r.hyperparams,params.cluster_params_r.posterior_hyperparams, params.cluster_params_r.suff_statistics)

    params.logsublikelihood_hist[1:burnout_period-1] = params.logsublikelihood_hist[2:burnout_period]
    params.logsublikelihood_hist[burnout_period] = log_likihood_l + log_likihood_r
    logsublikelihood_now = 0.0
    for i=1:burnout_period
        logsublikelihood_now += params.logsublikelihood_hist[i] *(1/(burnout_period-0.1))
    end
    if logsublikelihood_now != -Inf && logsublikelihood_now - params.logsublikelihood_hist[burnout_period] < 1e-2 # propogate abs change to other versions?
        # println(params.logsublikelihood_hist)
        params.splittable = true
    end

    return params
end


text:
module Solvers

    using LinearAlgebra

    using ..Config
    using ..Utils

    export LinearSolver, LUSolver, LUSolverLAPACK,
           factorize!, solve!

    include("solvers/linear/linear_solvers.jl")
    include("solvers/linear/lu_solver.jl")
    include("solvers/linear/lu_solver_lapack.jl")

    export JacobianParameters, JacobianParametersAD, JacobianParametersFD,
           JacobianParametersUser, getJacobianParameters,
           computeJacobian, computeJacobianAD, computeJacobianFD

    export computeJacobian, check_jacobian, print_jacobian

    export NonlinearSolver, AbstractNewtonSolver, NewtonSolver, QuasiNewtonSolver,
           residual_initial!, residual_absolute!, residual_relative!,
           print_solver_status, check_solver_converged, check_solver_status,
           solve!

    include("solvers/nonlinear/nonlinear_solvers.jl")
    include("solvers/nonlinear/jacobian.jl")
    include("solvers/nonlinear/abstract_fixed_point_solver.jl")
    include("solvers/nonlinear/abstract_newton_solver.jl")
    include("solvers/nonlinear/fixed_point_solver.jl")
    include("solvers/nonlinear/newton_solver.jl")
    include("solvers/nonlinear/quasi_newton_solver.jl")

    function __init__()
        default_params = (
            (:ls_solver, :lapack),
            (:nls_atol,  2eps()),
            (:nls_rtol,  2eps()),
            (:nls_stol,  2eps()),
            (:nls_atol_break,  1E-0),
            (:nls_rtol_break,  1E-0),
            (:nls_stol_break,  1E-0),
            (:nls_nmax,  10000),
            (:nls_nmin,  0),
            (:nls_nwarn, 100),
            (:nls_solver, QuasiNewtonSolver),
            (:jacobian_autodiff, true),
            (:jacobian_fd_ϵ, 8sqrt(eps())),
            (:quasi_newton_refactorize, 5),
            (:linesearch_nmax, 50),
            (:linesearch_armijo_λ₀, 1.0),
            (:linesearch_armijo_σ₀, 0.1),
            (:linesearch_armijo_σ₁, 0.5),
            (:linesearch_armijo_ϵ,  0.5),
        )

        for param in default_params
            add_config(param...)
        end
    end

end


text:
using Printf, Statistics
pp(x;kwargs...) = scatter(x[1,:], x[2,:]; kwargs...)

function dists(points)
    dims,n = size(points)
    dists = fill(Inf, n)
    for i = 1:n
        for j = 1:n
            i == j && continue
            d = norm(points[:,i] - points[:,j])
            dists[i] = min(d, dists[i])
        end
    end
    dists
end

function bluenoise(;
    dims = 2,
    nsamples = 100,
    iters = 10,
    points = rand(dims, nsamples)
    )
    dd = dists(points)
    @printf("Initial dist dist mean: %8.3g std: %8.3g\n", mean(dd), std(dd))
    hist1 = histogram(dd, bins=20, title= "Initial distance to closest")
    n = nsamples
    rdmax = [1/(2n), sqrt(1/(2sqrt(3)*n)), (1/(4sqrt(2)*n))^(1/3), (1/(8n))^(1/4)]
    w = 1 ./ rdmax
    w ./= sum(w)
    fadelength = iters ÷ 2
    fadein(i) = min(i,fadelength)*w[1]/fadelength
    fig1 = pp(points, title="Uniform random")
    scatter!(zeros(nsamples), points[2,:], m=(2,))
    scatter!(points[1,:], zeros(nsamples), m=(2,)) |> display
    points2 = copy(points)
    α = 0.005
    grad = zeros(dims)
    d = zeros(dims)
    @progress for i = 1:iters
        i == iters ÷ 2 && (α *= 0.2)
        perms = [sortperm(points[dim,:]) for dim = 1:dims]
        for p1 = 1:nsamples
            grad .*= 0
            for dim = 1:dims
                optimal = (-0.5 + perms[dim][p1])/n
                grad[dim] += fadein(i)*sign(optimal - points[dim,p1])
            end
            for p2 = 1:nsamples
                d .= points[:,p1] .- points[:,p2]
                nd = norm(d)
                # nd > 0.2 && continue
                grad .+= w[dims]*d./n/(0.001 + nd)^3
            end
            points2[:,p1] .+= α*grad
        end
        clamp!(points2,0.,1.)
        copyto!(points, points2)
        # if i % 2 == 0
        #     pp(points, xlims=(0,1), ylims=(0,1))
        #     scatter!(zeros(nsamples), points[2,:], m=(2,))
        #     scatter!(points[1,:], zeros(nsamples), m=(2,)) |> display
        # end
    end
    fig2 = pp(points, title="Optimized sample")
    scatter!(zeros(nsamples), points[2,:], m=(2,))
    scatter!(points[1,:], zeros(nsamples), m=(2,)) |> display
    dd = dists(points)
    @printf("Final dist dist mean: %8.3g std: %8.3g\n", mean(dd), std(dd))
    hist2 = histogram(dd, bins=20, title= "Optimized distance to closest")
    plot(fig1, fig2, hist1, hist2) |> display
    points
end

default(legend=false)
points = bluenoise(iters=100, nsamples=30)
# points = bluenoise(points=points, iters=10000, nsamples=20)

using StatPlots
histogram([std(dists(rand(1,20))) for i = 1:10000])

histogram([mean(dists(rand(1,20))) for i = 1:10000])
vline!([1/21 mean(diff((-0.5 .+ 1:20)./20)) mean(dists(points[1:1,:]))])


text:
export minres

"""
 x,flag,err,iter,resvec = minres(A,b,;x=[],maxIter=20,tol=1e-10,out=1)

 Minimal Residual Method (MINRES) for solving

    A x = b,

 where A is symmetric and possibly indefinite.

 Implementation follows page 26 of

 Choi, S.-C. T. (2006).
 Iterative Methods for Singular Linear Equations and Least-squares Problems.
 Phd thesis, Stanford University.

 See also:

 Paige, C. C., & Saunders, M. A. (1975).
 Solution of sparse indefinite systems of linear equations.
 SIAM Journal on Numerical Analysis.

 Required input:

   A       - function that performs matrix-vector product, e.g., x -> A*x
   b       - right hand side vector

 Optional input:

   x        - starting guess (optional)
   maxIter  - maximum number of iterations
   btol     - error tolerance for Lanczos step
   rtol     - error tolerance for estimated relative residual
   gtol     - error tolerance for estimated gradient norm
   condlim  - limit on condition number esitmate of A
   out      - flag for output (-1: no output, 0: only errors, 1: final status, 2: residual norm at each iteration)

 Output:

  x       - approximate solution
  flag    - exit flag (0: converged, -1: maxIter, -2: Lanczos, -3: condition)
  relres  - estimated relative residual at final iteration
  nrmG    - estimated gradient norm at final iteration
  nrmA    - estimated Frobenius norm of A
  conA    - condition number estimate
  phi     - relres for each iteration
"""
function minres{T1,T2}(A::SparseMatrixCSC{T1,Int},b::Array{T2,1}; kwargs...)
	x = zeros(promote_type(T1,T2),size(A,2)) # pre-allocate
	return minres(v -> At_mul_B!(1.0,A,v,0.0,x),b;kwargs...) # multiply with transpose of A for efficiency
end

minres(A,b::Vector;kwargs...) = minres(x -> A*x,b::Vector;kwargs...)

function minres(A::Function,b::Vector;x=[],sigma=0.0,btol=1e-10,rtol=1e-10,gtol=1e-10,condlim=1e10,maxIter=10,out=1,storeInterm::Bool=false)

    n      = length(b)
    nres   = norm(b)
	if nres==0; return zeros(eltype(b),n),-9,0.0,0.0,0.0,[0.0]; end

    if !isempty(x)
		b -= A(x)
	else
    	x = zeros(eltype(b),n)
	end
	if storeInterm
        X = zeros(eltype(b),n,maxIter)	# allocate space for intermediates
    end

    # initialize scalars (or vectors of size maxIter)
    alpha    = zeros(maxIter+1)
    beta     = zeros(maxIter+2)
    beta[2]  = norm(b)
    gamma    = zeros(maxIter+2)
    gamma[1] = beta[2]
    epsil    = zeros(maxIter+2)
    delta    = zeros(maxIter+2)
    delta[2] = 0.0
    phi      = zeros(maxIter+1)
    phi[1]   = beta[2]
    psi      = zeros(maxIter+1)
    tau      = zeros(maxIter+1)
    tau[1]   = beta[2]

    conA     =  1.0
    gammaMin =  Inf
    nrmG     =  0.0
    c        = -1.0
    s        =  0.0

    # allocate vectors of size n
    vk   = b/beta[2]
    vkm1 = zeros(eltype(b),n)
    dkm1 = zeros(eltype(b),n)
    dkm2 = zeros(eltype(b),n)
    tt   = zeros(eltype(b),n)

    if out>1
		println("=== minres ===")
        @printf "iter\t|r_k|\t\t|A*r_k|\t\t|A|\t\t\tcond(A)\n"
    end
    flag = -1

    k = 2
    nrmA = alpha[1]*beta[2]
    for k=2:maxIter+1
        alpha[k],beta[k+1],vk,vkm1 = LanczosStep!(A,vk,vkm1,beta[k],sigma=sigma,tol=btol)

        if beta[k+1] < btol
            flag = -2
            break
        end

        deltab     = c*delta[k] + s*alpha[k]
        gamma[k]   = s*delta[k] - c*alpha[k]
        epsil[k+1] = s*beta[k+1]
        delta[k+1] =-c*beta[k+1]

        c,s,gammab = symOrtho(gamma[k],beta[k+1])

        tau[k]   = c*phi[k-1]
        phi[k]   = s*phi[k-1]
        nrmG     = phi[k]*sqrt(gamma[k+1]^2+delta[k+1]^2)

        psi[k-1] = phi[k-1] * sqrt(gamma[k]^2+delta[k+1]^2)

        nrmA = max(nrmA,sqrt(beta[k]^2 + alpha[k]^2 + beta[k+1]^2))

        if gammab != 0
            # the following line is equivalent to: tt = copy(dkm1)
            tt = BLAS.blascopy!(n,dkm1,1,tt,1)
            # the following two lines are equivalent to dkm1 = vkm1 - deltab*dkm1
            dkm1 = BLAS.scal!(n,-deltab,dkm1,1)
            dkm1 = BLAS.axpy!(n,1.0,vkm1,1,dkm1,1)
            dkm1 = BLAS.axpy!(n,-epsil[k],dkm2,1,dkm1,1)
            dkm1 = BLAS.scal!(n,1.0/gammab,dkm1,1)
            dkm2 = BLAS.blascopy!(n,tt,1,dkm2,1)

            # the following line is equivalent to x    += tau[k]*dkm1
            x    = BLAS.axpy!(n,tau[k],dkm1,1,x,1)
            gammaMin = min(gammaMin,gammab)
            conA     = nrmA/gammaMin
        end
        if storeInterm; X[:,k-1] = x; end
        if out>1
            @printf "%d\t%1.2e\t%1.2e\t%1.2e\t\t%1.2e\n" k-1 phi[k]/nres nrmG nrmA conA
        end
        if (phi[k]/nres) < rtol && (nrmG < gtol)
            flag = 0
            break
        end
        if conA > condlim
            flag = -3
            break
        end
    end
    if out>=0
        if flag==-1
            println(@sprintf("minres iterated maxIter (=%d) times but reached only residual norm %1.2e instead of tol=%1.2e.",
            maxIter,phi[k]/nres,rtol))
        end
        if flag==-2
            println(@sprintf("minres stopped because beta (=%1.2e) became smaller than the tolerance (btol=%1.2e) at iteration %d.",
            beta[k+1],btol,k-1))
        end
        if flag==-3
            println(@sprintf("The estimated condition number (=%1.2e) in minres is bigger than the tolerance (=%1.2e) at iteration %d.",
            conA,condlim,k-1))
        end
        if (out>0) && flag==0
            println(@sprintf("minres converged at iteration %d. |A r_k|=%1.2e and |r_k|=%1.2e.",k-1,nrmG,phi[k]/nres))
        end
    end

    x = (storeInterm) ? X[:,1:k-1] : x
    return x,flag,phi[k]/nres,nrmG, nrmA,phi
end


text:
using DirectX
using Base.Test

include("dx9adl_tests.jl")


text:
"""
    SeisDecimate(in; <keyword arguments>)

Apply random or regular decimation to a multidimensional array of data. Input and output have the same dimensions.

# Arguments
- `in`: input data as 2D, or 3D,4D,5D tensors. The first dimension is time.

# Keyword arguments
- `mode="random"`: decimation mode. Random is default, else decimates uniformly.
- `perc=50`: percentage of traces equal to zero (decimated). Only for random mode
- `incx1=1`: consecutive traces zeroed in first dimension. Only for regular decimation
- `incx2=1`: consecutive traces zeroed in second dimension.
- `incx3=1`: consecutive traces zeroed in third dimension.
- `incx4=1`: consecutive traces zeroed in fourth dimension.

# Example
```julia
julia> d = SeisLinearEvents(); deci = SeisDecimate(d);
```

*Credits: Aaron Stanton,2017*

"""
function SeisDecimate(in;mode="random",perc=50,incx1=1,incx2=1,incx3=1,incx4=1)

    if (mode=="random") # decimate data randomly
        out = copy(in)
        out = reshape(out,size(in,1),:)
	mask = rand(1,size(in,2)*size(in,3)*size(in,4)*size(in,5));
	mask[(LinearIndices(mask .< perc/100))[findall(mask .< perc/100)]] .= 0;
	mask[(LinearIndices(mask .>= perc/100))[findall(mask .>= perc/100)] ] .= 1;
	for it = 1 : size(in,1)
	    out[it,:] = out[it:it,:].*mask;
	end
    out = reshape(out,size(in));
    else # decimate data regularly with respect to 4 spatial dimensions
        out = zero(in)
	if (size(in,5)) > 1
	    out[:,1:incx1:end,1:incx2:end,1:incx3:end,1:incx4:end] = in[:,1:incx1:end,1:incx2:end,1:incx3:end,1:incx4:end]
	elseif (size(in,4)) > 1
            out[:,1:incx1:end,1:incx2:end,1:incx3:end] = in[:,1:incx1:end,1:incx2:end,1:incx3:end]
        elseif (size(in,3)) > 1
            out[:,1:incx1:end,1:incx2:end] = in[:,1:incx1:end,1:incx2:end]
        elseif (size(in,2)) > 1
            out[:,1:incx1:end] = in[:,1:incx1:end]
        end

    end

    return out
end


text:
using LinearAlgebra, SparseArrays, Logging, Plots, ProgressMeter

"""
    Bad eigen value distribution, they are distributed from 0.001 to 1 like a
    geometric series.

    * A lot of small eigenvalues are clustered close to zero, few larger ones are
    far away from the other Eigenvalues.

    Parameters:
        rho: Number:
            An number between 0 and 1 that paramaterize the distribution of the
            eigenvalues for the PSD matrix.
        N=20:
            The size of the PSD matrix.

    returns:
        a diagonal matrix.
"""
function GetNastyPSDMatrix(
    N=20,
    rho::Number=0.9,
    eigen_min=1e-4,
    eigen_max=1,
    inverted=false
)
    @assert rho <= 1 && rho >= 0
    EigenValues = zeros(N)
    EigenValues[1] = eigen_min
    for IdexI in 2:N
        EigenValues[IdexI] = eigen_min +
            ((IdexI - 1)/(N - 1))*(eigen_max - eigen_min)*rho^(N - IdexI)  # formulas
    end
    if inverted
        return Diagonal(1 + eigen_min - EigenValues)
    end
    return Diagonal(EigenValues)
end

"""
    Get a diagonal PSD matrix whose eigenvalues are squared of a standard
    normal distributions.
"""
function GetNormalPSDMatrix(N=20)
    v = randn(N).^2
    return Diagonal(v)
end

"""
    Return a diagonal matrix whose entries are draw from the uniform distribution on the
    [0, 1] interval.
"""
function GetUniformPSDMatrix(N)
    v = rand(N)
    return Diagonal(v)
end

"""
    Given the opeartr, RHS vector, and the list of guesses, compute the
    relative error in energy norm induced by linear operator A.
"""
function ResRelEnergyNorm(A, b, xs)
    x̄ = A\b
    x0 = xs[1]
    return xs .|> x -> dot(x - x̄, A, x - x̄)/dot(x0 - x̄, A, x0 - x̄)
end

"""
    Compute the Theoretical Error Bound for a given Linear Opeartor and the
    iterations count that you are looking at.
"""
function TheoreticalErrorBound(A::AbstractMatrix, k)
    κ = cond(A)
    y = (sqrt(κ) - 1)/(sqrt(κ) + 1)
    return 2*(y^k + y^(-k))^(-1)
end

"""
    Given an instance of a plot, and the file name, it
    stores it to a folder named: 'plots' under the
    directory of the running script.
"""
function SaveFigToCurrentScriptDir(fig::Plots.Plot, file_name::String)
    scriptDir = @__DIR__
    directory = scriptDir*"/plots"
    path = mkpath(directory)
    savefig(fig, directory*"/"*file_name)
    return path
end


text:
using Comonicon
using Comonicon.Types
using Comonicon.BuildTools: precompile_script, install
using Test

module Dummy
using Comonicon
using Test

"""
foo foo

# Arguments

- `x`: an argument

# Options

- `--foo <foo>`: foo foo
"""
@cast function foo(x::Int, y; foo::Int = 1, hulu::Float64 = 2.0, flag::Bool = false) where {T}
    @test x == 1
    @test y == "2.0"
    @test foo == 2
    @test hulu == 3.0
    @test flag == true
end

"""
goo goog gooasd dasdas

goo asdas dasd assadas

# Arguments

- `ala`: ala ahjsd asd wvxzj
- `gaga`: djknawd sddasd kw

# Options

- `-g,--giao <name>`: huhuhuhuhuuhu

# Flags

- `-f,--flag`: dadsa fasf gas
"""
@cast function goo(ala::Int, gaga; giao = "Bob", flag::Bool = false)
    @test ala == 1
    @test gaga == "2.0"
    @test giao == "Sam"
    @test flag == true
end

"""
tick tick.

# Arguments

- `xx`: xxxxxxxxxxxx
- `yy`: yyyyyyyyyyyy
"""
@cast function tick(xx::Int, yy::Float64 = 1.0)
    @test xx == 1
    @test yy in [1.0, 2.0]
end

@main name = "dummy" doc = """
    dummy command. dasdas dsadasdnaskdas dsadasdnaskdas
    sdasdasdasdasdasd adsdasdas dsadasdas dasdasd dasda
    """
end

@test Dummy.command_main(String["foo", "1.0", "2.0", "--foo", "2", "--hulu=3.0", "--flag"]) == 0
@test Dummy.command_main(String["foo", "1.0", "2.0", "--foo", "2", "--hulu=3.0", "-f"]) == 0
@test Dummy.command_main(String["goo", "1.0", "2.0", "-gSam", "-f"]) == 0
@test_throws ErrorException LeafCommand(
    identity;
    name = "foo",
    options = [Option("huhu"; short = true)],
)
@test Dummy.command_main(String["tick", "1.0", "2.0"]) == 0
@test Dummy.command_main(String["tick", "1.0"]) == 0

@test precompile_script(Dummy) == """
using Main.Dummy;
Main.Dummy.command_main(["-h"]);
Main.Dummy.command_main(["goo", "-h"]);
Main.Dummy.command_main(["tick", "-h"]);
Main.Dummy.command_main(["foo", "-h"]);
"""

empty!(ARGS)
append!(ARGS, ["2", "--opt1", "3"])

"""
ArgParse example implemented in Comonicon.

# Arguments

- `x`: an argument

# Options

- `--opt1 <arg>`: an option
- `-o, --opt2 <arg>`: another option

# Flags

- `-f, --flag`: a flag
"""
@main function main(x; opt1 = 1, opt2::Int = 2, flag = false)
    @test flag == false
    @test x == "2"
    @test opt1 == "3"
    @test opt2 == 2
end

Comonicon.install(
    Dummy;
    bin = Comonicon.PATH.project("test", "bin"),
    completion = false,
    quiet = false,
)

@test isfile(Comonicon.PATH.project("test", "bin", "dummy"))
@test isfile(Comonicon.PATH.project("test", "bin", "dummy.jl"))

@testset "default_name" begin
    @test Comonicon.Parse.default_name("Foo.jl") == "foo"
    @test Comonicon.Parse.default_name(sin) == "sin"
end

cmd = @cast(f_issue_47(xs::Int...) = xs)

@testset "issue/#47" begin
    @test cmd.args[1].type == Int
    @test cmd.args[1].vararg == true
end


text:
# This file is a part of ElasticArrays.jl, licensed under the MIT License (MIT).

using ElasticArrays
using Test

using Random, LinearAlgebra

@testset "elasticarray" begin
    test_dims = (3, 2, 4)
    test_kernel_dims = Base.front(test_dims)

    function test_A(test_code)
        A = rand!(Array{Int}(undef, test_dims...), 0:99)
        test_code(A)
    end

    function test_E(test_code)
        E = rand!(ElasticArray{Int}(undef, test_dims...), 0:99)
        test_code(E)
    end

    function test_E_A(test_code)
        E = rand!(ElasticArray{Int}(undef, test_dims...), 0:99)
        A = rand!(Array{Int}(undef, test_dims...), 0:99)
        test_code(E, A)
    end

    function test_E1_E2(test_code)
        E1 = rand!(ElasticArray{Int}(undef, test_dims...), 0:99)
        E2 = rand!(ElasticArray{Int}(undef, test_dims...), 0:99)
        test_code(E1, E2)
    end

    function test_E_V(test_code)
        E = ElasticArray{Float64}(undef, test_kernel_dims..., 0)
        V = Vector{Array{Float64,length(test_kernel_dims)}}()
        test_code(E, V)
    end

    lastdim_slice_idxs(A::AbstractArray{T,N}, i::Integer) where {T,N} = (ntuple(_ -> :, Val(N - 1))..., i)

    test_comp(E::ElasticArray, V::Vector{<:Array}) =
        all(i -> @view(E[lastdim_slice_idxs(E, i)...]) == V[i], eachindex(V))

    @testset "ctors" begin
        @test (@inferred ElasticArray{Int}(undef, 2, 3, 4)).kernel_size == (2, 3)
        @test ElasticArray{Int}(undef, 2, 3, 4).kernel_length == Base.MultiplicativeInverses.SignedMultiplicativeInverse(2 * 3)
        @test size(ElasticArray{Int}(undef, 2, 3, 4).data) == (2 * 3 * 4,)

        @test fill!((@inferred ElasticArray{Int}(undef, 2, 3, 4)), 42) == fill!(ElasticArray{Int}(undef, 2, 3, 4), 42)

        @test size(@inferred(ElasticArray{Float64}(undef, 0, 2, 0))) == (0, 2, 0)
        @test_throws ArgumentError ElasticArray{Float64}(undef, 0, 2, 3)
    end

    @testset "size, length and index style" begin
        let A = ElasticArray{Float32}(undef, 2, 5)
            @test vec(A) === A.data
            @test vec(A) == vec(Array(A))
        end

        @test (4,) == @inferred size(@inferred ElasticArray{Int}(undef, 4))
        @test (2,3,4) == @inferred size(@inferred ElasticArray{Int}(undef, 2,3,4))

        test_E() do E
            @test length(E) == prod(size(E))
            @test IndexStyle(E) == IndexLinear()
            @test eachindex(E) == eachindex(vec(E))
            @test sizeof(E) == sizeof(E.data)
            @test Base.eltype(E) == Base.eltype(E.data)
        end
    end

    @testset "getindex and setindex!" begin
        test_E_A() do E, A
            for i in eachindex(E, A)
                E[i] = A[i]
            end
            @test all(i -> E[i] == A[i], eachindex(E, A))
            @test all(i -> E[i] == A[i], CartesianIndices(size(A)))
        end

        @test all(x -> x == 42, @inferred fill!(ElasticArray{Int}(undef, 2,3,4), 42))
    end

    @testset "equality" begin
        test_E() do E
            @test E == @inferred deepcopy(E)
        end

        test_E1_E2() do E1, E2
            @test E1 != E2
        end

        @test fill!(ElasticArray{Int}(undef, 2, 3, 4), 0) == fill!(ElasticArray{Int}(undef, 2, 3, 4), 0)
        @test fill!(ElasticArray{Int}(undef, 2, 3, 4), 0) != fill!(ElasticArray{Int}(undef, 2, 4, 3), 0)
    end

    @testset "mightalias and dataids" begin
        E1 = ElasticArray{Int}(undef, 10, 5)
        E2 = ElasticArray{Int}(undef, 10, 5)
        @test Base.dataids(vec(E1)) == @inferred Base.dataids(E1)
        @test @inferred !Base.mightalias(E1, E2)
        @test @inferred !Base.mightalias(view(E1, 2:3, 1:2), view(E1, 4:5, 1:2))
        @test @inferred Base.mightalias(view(E1, 2:4, 1:2), view(E1, 3:5, 1:2))
        @test @inferred !Base.mightalias(view(E1, 2:4, 1:2), view(E2, 3:5, 1:2))
    end

    @testset "copyto!, conversion ctor, convert and similar" begin
        test_E_A() do E, A
            @test E === @inferred copyto!(E, A)
            @test E == A
        end

        test_E_A() do E, A
            A2 = Array(deepcopy(E))
            @test E === @inferred copyto!(E, 3, A, 5, 7)
            copyto!(A2, 3, A, 5, 7)
            @test E == A2
        end

        test_E_A() do E, A
            @test A === @inferred copyto!(A, E)
            @test A == E
        end

        test_E_A() do E, A
            A2 = deepcopy(A)
            @test A === @inferred copyto!(A, 3, E, 5, 7)
            copyto!(A2, 3, Array(deepcopy(E)), 5, 7)
            @test A == A2
        end

        test_E1_E2() do E1, E2
            @test E1 === @inferred copyto!(E1, E2)
            @test E1 == E2
        end

        test_E1_E2() do E1, E2
            A1 = Array(deepcopy(E1))
            A2 = Array(deepcopy(E2))
            @test E1 === @inferred copyto!(E1, 3, E2, 5, 7)
            copyto!(A1, 3, A2, 5, 7)
            @test E1 == A1
        end

        test_A() do A
            E = @inferred ElasticArray{Float64}(A)
            @test E isa ElasticArray
            @test E == A
            @test eltype(E) == Float64
        end

        test_A() do A
            E = @inferred ElasticArray{Float32}(A)
            @test E isa ElasticArray
            @test E == A
            @test eltype(E) == Float32
        end

        test_A() do A
            E = @inferred ElasticArray{Float64,3}(A)
            @test E isa ElasticArray
            @test E == A
            @test eltype(E) == Float64
        end

        test_A() do A
            E = @inferred ElasticArray{Float32,3}(A)
            @test E isa ElasticArray
            @test E == A
            @test eltype(E) == Float32
        end

        test_A() do A
            E = @inferred ElasticArray(A)
            @test E isa ElasticArray
            @test E == A
            @test eltype(E) == eltype(A)
        end

        test_A() do A
            E = @inferred convert(ElasticArray{Float64,3}, A)
            @test E isa ElasticArray{Float64,3,2}
            @test E == A
            @test eltype(E) == Float64

            @test convert(ElasticArray{Float64,3,2}, E) === E
            @test convert(ElasticArray{Float64,3}, E) === E
            @test convert(ElasticArray{Float64}, E) === E
            @test convert(ElasticArray, E) === E

            @test E == @inferred convert(ElasticArray{Float64,3,2}, A)
            @test_throws ArgumentError convert(ElasticArray{Float64,3,1}, A)
        end

        test_A() do A
            E = @inferred convert(ElasticArray{Float32,3}, A)
            @test E isa ElasticArray{Float32,3,2}
            @test E == A
            @test eltype(E) == Float32

            @test convert(ElasticArray{Float32,3,2}, E) === E
            @test convert(ElasticArray{Float32,3}, E) === E
            @test convert(ElasticArray{Float32}, E) === E
            @test convert(ElasticArray, E) === E

            @test E == @inferred convert(ElasticArray{Float32,3,2}, A)
            @test_throws ArgumentError convert(ElasticArray{Float32,3,4}, A)
        end

        test_A() do A
            E = @inferred convert(ElasticArray{Float64}, A)
            @test E isa ElasticArray
            @test E == A
            @test eltype(E) == Float64
        end

        test_A() do A
            E = @inferred convert(ElasticArray{Float32}, A)
            @test E isa ElasticArray
            @test E == A
            @test eltype(E) == Float32
        end

        test_A() do A
            E = @inferred convert(ElasticArray, A)
            @test E isa ElasticArray
            @test E == A
            @test eltype(E) == eltype(A)
        end

        @test typeof(@inferred similar(ElasticArray{Int}, (2,3,4))) == ElasticArray{Int,3,2,Vector{Int}}
        @test size(@inferred similar(ElasticArray{Int}, (2,3,4))) == (2,3,4)

        @test typeof(@inferred similar(ElasticArray{Int}(undef,4,5), (2,3,4))) == ElasticArray{Int,3,2,Vector{Int}}
        @test size(@inferred similar(ElasticArray{Int}(undef,4,5), (2,3,4))) == (2,3,4)

        @test typeof(@inferred similar(ElasticArray{Int}(undef,4,5), ())) == Array{Int,0}
        @test size(@inferred similar(ElasticArray{Int}(undef,4,5), ())) == ()

        @test size(@inferred(similar(ElasticArray{Int}(undef, 3, 4), Float32, (0, 2)))) == (0, 0)
    end

    @testset "pointer and unsafe_convert" begin
        test_E() do E
            @test pointer(E) == pointer(vec(E))
            @test pointer(E, 4) == pointer(vec(E), 4)
        end

        test_E() do E
            @test Base.unsafe_convert(Ptr{eltype(E)}, E) == Base.unsafe_convert(Ptr{eltype(E)}, vec(E))
        end
    end

    @testset "resize!" begin
        function resize_test(delta::Integer)
            test_E() do E
                A = Array(deepcopy(E))
                new_size = (Base.front(size(E))..., size(E, ndims(E)) + delta)
                cmp_idxs = (Base.front(axes(E))..., 1:(last(size(E)) + min(0, delta)))
                @test E === @inferred sizehint!(E, new_size...)
                @test E === @inferred resize!(E, new_size...)
                @test size(E) == new_size
                @test E[cmp_idxs...] == A[cmp_idxs...]
            end
        end

        resize_test(0)
        resize_test(2)
        resize_test(-2)

        function resize_lastdim_test(delta::Integer)
            test_E() do E
                A = Array(deepcopy(E))
                new_size = (Base.front(size(E))..., size(E, ndims(E)) + delta)
                cmp_idxs = (Base.front(axes(E))..., 1:(last(size(E)) + min(0, delta)))
                @test E === @inferred ElasticArrays.sizehint_lastdim!(E, size(E, ndims(E)) + delta)
                @test E === @inferred ElasticArrays.resize_lastdim!(E, size(E, ndims(E)) + delta)
                @test size(E) == new_size
                @test E[cmp_idxs...] == A[cmp_idxs...]
            end
        end

        resize_lastdim_test(0)
        resize_lastdim_test(2)
        resize_lastdim_test(-2)
    end

    @testset "append! and prepend!" begin
       test_A() do A
            E = @inferred convert(ElasticArray{Float64}, A)
            @test E isa ElasticArray
            @test E == A
            @test eltype(E) == Float64
        end

        test_E_V() do E, V
            dims = Base.front(size(E))
            for i in 1:4
                push!(V, rand(dims...))
                @inferred append!(E, last(V))
            end
            @test size(E) == (dims..., length(V))
            @test test_comp(E, V)
        end

        test_E_V() do E, V
            dims = Base.front(size(E))
            for i in 1:4
                pushfirst!(V, rand(dims...))
                @inferred prepend!(E, first(V))
            end
            @test size(E) == (dims..., length(V))
            @test test_comp(E, V)
        end

        test_E() do E
            kernel_size = Base.front(size(E))
            kernel_length = last(size(E))
            V = rand(Int, prod(kernel_size))
            @inferred append!(E, (el for el in V))
            @test size(E) == (kernel_size..., kernel_length + 1)
            @test vec(E[:, :, kernel_length + 1]) == V
        end

        test_E() do E
            kernel_size = Base.front(size(E))
            kernel_length = last(size(E))
            V = rand(Int, prod(kernel_size))
            @inferred prepend!(E, (el for el in V))
            @test size(E) == (kernel_size..., kernel_length + 1)
            @test vec(E[:, :, 1]) == V
        end
    end

    @testset "basic math" begin
        T = Float64
        E1 = rand!(ElasticArray{T}(undef, 9, 9))
        E2 = rand!(ElasticArray{T}(undef, 9, 9))
        E3 = rand!(ElasticArray{T}(undef, 9, 7))

        A1 = Array(E1)
        A2 = Array(E2)
        A3 = Array(E3)

        @test @inferred(2 * E1) isa ElasticArray{T,2}
        @test 2 * E1 == 2 * A1

        @test @inferred(E1 .+ 2) isa ElasticArray{T,2}
        @test E1 .+ 2 == A1 .+ 2

        @test @inferred(E1 + E2) isa ElasticArray{T,2}
        @test E1 + E2 == A1 + A2

        @test @inferred(E1 * E2) isa ElasticArray{T,2}
        @test E1 * E2 == A1 * A2
        @test E1 * E3 == A1 * A3

        @test E1^3 == A1^3
        @test inv(E1) == inv(A1)
        @test inv(E1) isa ElasticArray{T,2}
    end

    @testset "linear algebra" begin
        lmul!(cholesky(Array{Float32}([1 0; 0 1])).L, view(ElasticArray(Float32[1 1; 1 1]), :, 1)) == [1, 1]
    end

    @testset "ElasticVector" begin
        @test ElasticVector{Int} === ElasticArray{Int,1,0}
        @test ElasticArray(rand(Float32, 2)) isa ElasticVector{Float32}
        @test size(ElasticVector{Float64}(undef, 2)) == (2,)
        @test size(ElasticVector([-1, 1, 0])) == (3,)
        @test ElasticVector([1, 2, 3]) == [1, 2, 3]
        @test eltype(ElasticVector([1.0, 2.1])) === Float64
        @test eltype(ElasticVector{Float32}([1, 2, 3])) === Float32
        @test_throws Exception ElasticVector([1 2 3])
        @test_throws Exception ElasticVector(zeros(2, 3, 4))
    end

    @testset "ElasticMatrix" begin
        @test ElasticMatrix{Int} === ElasticArray{Int,2,1}
        @test ElasticArray(rand(Float32, 2, 3)) isa ElasticMatrix{Float32}
        @test size(ElasticMatrix{Float64}(undef, 2, 3)) == (2, 3)
        @test size(ElasticMatrix([1 2 3; 4 5 6])) == (2, 3)
        @test ElasticMatrix([1 2; 3 4]) == [1 2; 3 4]
        @test eltype(ElasticMatrix([1.0 2; 3 4.1])) === Float64
        @test eltype(ElasticMatrix{Float32}([1 2])) === Float32
        @test_throws Exception ElasticMatrix([1, 2, 3])
        @test_throws Exception ElasticMatrix(zeros(2, 3, 4))
    end

    @testset "Regression tests" begin
        let A = ElasticArray([3.2 0.5; 0.5 2.0])
            @test eigvals(A) == eigvals(Array(A))
        end
    end
end


text:
using PackageCompiler
create_sysimage([:CSV, :DataFrames, :IJulia, :Pluto, :PlutoUI,
                :MeshArrays, :ClimateModels, :MITgcmTools,
                :UnicodePlots, :CairoMakie, :Plots];
                precompile_execution_file = "sysimage/warmup.jl",
                sysimage_path="ExampleSysimage.so",
                cpu_target = PackageCompiler.default_app_cpu_target())

max_stars_repo_name: gaelforget/Notebooks

max_issues_repo_name: gaelforget/Notebooks

max_forks_repo_name: gaelforget/Notebooks


text:
abstract type Transform end

Base.broadcastable(tr::Transform) = Ref(tr)

struct ComposedTransform{TR1<:Transform,TR2<:Transform} <: Transform
  """
  Outer transform.
  """
  t1::TR1
  """
  Inner transform.
  """
  t2::TR2
end
(t::ComposedTransform)(p::AbstractVector) = t.t1(t.t2(p))

AbstractTrees.children(tr::ComposedTransform) = (tr.t1, tr.t2)

transforms(tr::ComposedTransform) = Leaves(tr)

Base.:∘(t1::Transform, t2::Transform) = ComposedTransform(t1, t2)

inv(t::ComposedTransform) = inv(t.t2) ∘ inv(t.t1)

struct Transformed{O,TR<:Transform}
  obj::O
  transf::TR
  Transformed{O,TR}(obj::O, args...) where {O,TR<:Transform} = Transformed{O,TR}(obj, TR(args...))
  Transformed{O,TR}(obj::O, transf::TR) where {O,TR<:Transform} = new{O,TR}(obj, transf)
  Transformed{TR1,TR2}(obj::TR1, transf::TR2) where {TR1<:Transformed,TR2<:Transform} = Transformed(obj.obj, transf ∘ obj.transf)
  Transformed(obj::O, transf::TR) where {O,TR<:Transform} = Transformed{O,TR}(obj, transf)
end
(tr::Transformed)(p) = tr.obj(inv(tr.transf)(p))

(≈)(x::Transformed, y::Transformed) = typeof(x) == typeof(y) && x.obj ≈ y.obj && x.transf ≈ y.transf

include("transforms/scaling.jl")
include("transforms/translation.jl")
include("transforms/rotation.jl")


text:
FRGRealTime.ppfunps_delta(10.0, 10.0, 20.0, m2fun(20.0), 145.0)

plot(
    k -> FRGRealTime.dkF1TildeintqsAll_delta2(
        50.0,
        1.0,
        10.0,
        k,
        m2fun(k),
        145.0,
    ),
    10,
    50,
)

dkVReintqs_delta2(p0, ps, q0, qsmax, k, m, T, Npi, lam4pik)
plot(
    k -> FRGRealTime.dkVReintqs_delta2(
        50.0,
        1.0,
        5.0,
        10.0,
        k,
        m2fun(k),
        145.0,
        4.0,
        lamfun(k),
    ),10,50
)

plot(
    k -> FRGRealTime.dkF1Tildeps_delta2(50.0, 10.0, k, m2fun(k), 145.0),
    10,
    50,
)

FRGRealTime.fastpropReintqs_delta2(
    50.0,
    1.0,
    145.0,
    1.0,
    800.0,
    4.0,
    m2fun,
    lamfun;maxevals=2000000)


text:
function fieldline_images(bx::Function, by::Function, bz::Function, ds::Float64, grd::Grid)

    deltaPhi = 2*pi / grd.Nz

    chnl = RemoteChannel(()->Channel{Bool}(), 1)
    p = Progress(grd.Nk, desc="Mapping fieldlines... ")
    @async while take!(chnl)
        next!(p)
    end

    img = DArray((6, grd.Nk), workers(), [1, length(workers())]) do inds

        Nk = length(inds[2])
        points0 = grd.points[:, inds[2]]
        imgpoints = zeros(6, Nk)

        for k=1:Nk

            try
                imgpoints[1:3,k] = trace_fieldline(points0[:,k]..., bx, by, bz, ds, deltaPhi)
                imgpoints[4:6,k] = trace_fieldline(points0[:,k]..., bx, by, bz, -ds, deltaPhi)
            catch
                imgpoints[1:3,k] = [points0[:,k]..., deltaPhi]
                imgpoints[4:6,k] = [points0[:,k]..., deltaPhi]
            end
            put!(chnl, true)

        end
        imgpoints
    end

    put!(chnl, false)
    return Matrix(img)
end

function fieldline_derivatives(img::Matrix{Float64}, mx::Int64, my::Int64, MinvT::Matrix{Float64}, grd::Grid)

    dS1 = img[3,:]
    dS2 = img[6,:]
    fwd_pts = img[1:2,:]
    bck_pts = img[4:5,:]

    # make diagonal matrices for distance coefficients
    _A = Diagonal(-dS1 ./ (dS2.^2 + dS1.*dS2))
    _B = Diagonal(1 ./ dS2 - 1 ./ dS1)
    _C = Diagonal(dS2 ./ (dS1.^2 + dS1.*dS2))
    _D = Diagonal(1 ./ (dS2.^2 + dS1.*dS2))
    _E = Diagonal(-1 ./ (dS1.*dS2))
    _F = Diagonal(1 ./ (dS1.^2 + dS1.*dS2))

    # make off diagonal matrices for kronecker product
    Ia = spdiagm(-1 => ones(Float64, grd.Nz-1))
    Ia[1,end] = 1
    Ib = sparse(I, grd.Nz, grd.Nz)
    Ic = spdiagm(1 => ones(Float64, grd.Nz-1))
    Ic[end,1] = 1

    # make fieldline map matrices

    chnl = RemoteChannel(()->Channel{Bool}(), 1)
    p = Progress(grd.Nk, desc="Generating fieldline derivatives... ")
    @async while take!(chnl)
        next!(p)
    end

    is = vcat([k*ones(mx*my) for k=1:grd.Nk]...)
    c = pmap(1:grd.Nk) do k

        fwd_row_dat, fwd_row_j = interpolation_row(fwd_pts[:,k]..., mx, my, MinvT, grd)
        bck_row_dat, bck_row_j = interpolation_row(bck_pts[:,k]..., mx, my, MinvT, grd)
        if any(x -> !(0 < x < grd._Nx*grd._Ny), fwd_row_j)
            fwd_row_dat, fwd_row_j = interpolation_row(grd.points[:,k]..., mx, my, MinvT, grd)
        end
        if any(x -> !(0 < x < grd._Nx*grd._Ny), bck_row_j)
            bck_row_dat, bck_row_j = interpolation_row(grd.points[:,k]..., mx, my, MinvT, grd)
        end

        put!(chnl, true)
        hcat(fwd_row_dat, fwd_row_j, bck_row_dat, bck_row_j)
    end

    M = vcat(c...)
    fwd_dat = M[:,1]
    fwd_j = Int64[M[:,2]...]
    bck_dat = M[:,3]
    bck_j = Int64[M[:,4]...]

    # fwd_j = Int64[]
    # fwd_dat = Float64[]
    # bck_j = Int64[]
    # bck_dat = Float64[]
    # for k=1:grd.Nk
    #     fwd_row_dat, fwd_row_j = interpolation_row(fwd_pts[:,k]..., mx, my, MinvT, grd)
    #     bck_row_dat, bck_row_j = interpolation_row(bck_pts[:,k]..., mx, my, MinvT, grd)
    #     fwd_j = [fwd_j; fwd_row_j]
    #     fwd_dat = [fwd_dat; fwd_row_dat]
    #     bck_j = [bck_j; bck_row_j]
    #     bck_dat = [bck_dat; bck_row_dat]
    #
    #     @info "$k/$(grd.Nk)"
    # end

    FWD = sparse(is, fwd_j, fwd_dat, grd.Nk, grd._Nx*grd._Ny) * transpose(grd.Proj)
    BCK = sparse(is, bck_j, bck_dat, grd.Nk, grd._Nx*grd._Ny) * transpose(grd.Proj)

    # compute fieldline derivatives
    Ds = kron(Ia, _A*BCK) + kron(Ib, _B) + kron(Ic, _C*FWD)
    Dss = kron(Ia, _D*BCK) + kron(Ib, _E) + kron(Ic, _F*FWD)

    put!(chnl, false)
    return Ds, Dss
end

#
#
# # NOTE: tracing parallel or anti-parallel to b can be achieved with positive or negative ds
# function fieldline_map_matrix(bx::Function, by::Function, bz::Function, ds::Float64, mx::Int64, my::Int64, MinvT::Matrix{Float64}, Nz::Int64, grd::Grid)
#
#     deltaPhi = 2*pi / Nz
#
#     # define function we want to map onto grd
#     function f(x0, y0, grd)
#
#         x, y, deltaS = trace_fieldline(x0, y0, bx, by, bz, ds, deltaPhi)
#         row_dat, row_js = interpolation_row(x, y, mx, my, MinvT, grd)
#         dS_vec = zeros(Float64, mx*my)
#         dS_vec[1] = deltaS
#
#         return hcat(row_dat, row_js, dS_vec)
#     end
#
#     M = grid_map(f, (mx*my,3), grd)
#     dS = M[1:mx*my:end,3]
#     dat = M[:,1]
#     js = Int64[M[:,2]...]
#     is = vcat([k*ones(Int64, mx*my) for k=1:grd.Nk]...)
#
#     matrix = sparse(is, js, dat, grd.Nk, grd._Nx*grd._Ny) * transpose(grd.Proj)
#     return matrix, dS
# end
#
#
# function fieldline_derivatives(bx::Function, by::Function, bz::Function, ds::Float64, mx::Int64, my::Int64, MinvT::Matrix{Float64}, Nz::Int64, grd::Grid)
#
#     # compute forward and backward maps
#     FM, dS1 = fieldline_map_matrix(bx, by, bz, ds, mx, my, MinvT, Nz, grd)
#     @info "Calculated forward map."
#     BM, dS2 = fieldline_map_matrix(bx, by, bz, -ds, mx, my, MinvT, Nz, grd)
#     @info "Calculated backward map."
#
#     # make diagonal matrices for distance coefficients
#     _A = Diagonal(-dS1 ./ (dS2.^2 + dS1.*dS2))
#     _B = Diagonal(1 ./ dS2 - 1 ./ dS1)
#     _C = Diagonal(dS2 ./ (dS1.^2 + dS1.*dS2))
#     _D = Diagonal(1 ./ (dS2.^2 + dS1.*dS2))
#     _E = Diagonal(-1 ./ (dS1.*dS2))
#     _F = Diagonal(1 ./ (dS1.^2 + dS1.*dS2))
#
#     # make off diagonal matrices for kronecker product
#     Ia = spdiagm(-1 => ones(Float64, Nz-1))
#     Ia[1,end] = 1
#     Ib = sparse(I, Nz, Nz)
#     Ic = spdiagm(1 => ones(Float64, Nz-1))
#     Ic[end,1] = 1
#
#     # construct derivative matrices
#     Ds = kron(Ia, _A*BM) + kron(Ib, _B) + kron(Ic, _C*FM)
#     Dss = kron(Ia, _D*BM) + kron(Ib, _E) + kron(Ic, _F*FM)
#
#     return Ds, Dss
# end

function flux_surface_average(psi::Function, psi_val, Ntheta, r0::Vector{Float64}, mx::Int64, my::Int64, MinvT::Matrix{Float64}, grd::Grid; dr=0.01)

    points = zeros(2, Ntheta)
    thetas = LinRange(0, 2*pi, Ntheta+1)
    for i=1:Ntheta

        et = Float64[cos(thetas[i]), sin(thetas[i])]
        r = r0 + dr*et
        sgn = sign(psi(r...) - psi_val)

        while true
            y = r + dr*et
            if sgn*(psi(y...) - psi_val) < 0
                points[:,i] = 0.5*(y + r)
                break
            end
            r = y
        end

    end

    return line_average(points, mx, my, MinvT, grd)
end


text:
function randargmin2(a)
    indices = eachindex(a)
    y = iterate(indices)
    y === nothing && throw(ArgumentError("collection must be non-empty"))
    (idx, state) = y
    minval = a[idx]
    bestidx = idx
    bestcount = 1
    y = iterate(indices, state)
    while y !== nothing
        (idx, state) = y
        curval = a[idx]
        if isless(curval, minval)
            minval = curval
            bestidx = idx
            bestcount = 1
        elseif isequal(curval, minval)
            bestcount += 1
            rand() * bestcount < 1 && (bestidx = idx)
        end
        y = iterate(indices, state)
    end
    bestidx
end

max_stars_repo_path: Chapter02/01. Finding the index of a random minimum element in an array/randargmin2.jl


text:
using Pixell
using Test
import Pixell: degree, arcminute

using DelimitedFiles

include("test_geometry.jl")  # creating geometries and sky ↔ pix
include("test_enmap.jl")     # enmap features and manipulation
include("test_transforms.jl")
include("test_distance_transform.jl")
include("test_io.jl")
include("test_plot.jl")


text:
# test mlogloss
let pd = [0.0 0.0 0.0 1.0;1.0 0.0 0.0 0.0], gt = [0.0 0.0 0.0 1.0;1.0 0.0 0.0 0.0]

    @test AutoML.mlogloss(pd,gt) == 1

end


text:
struct TwoDLinearSolverSplineMass
    nx1::Int
    nx2::Int
    factor::Float64
    eigvals1::Vector{Float64}
    eigvals2::Vector{Float64}
    wk::Array{ComplexF64,2}

    function TwoDLinearSolverSplineMass(nx1, nx2, eigvals1, eigvals2)
        factor = 1.0
        wk = zeros(ComplexF64, nx1, nx2)
        return new(nx1, nx2, factor, eigvals1, eigvals2, wk)
    end
end

function solve(solver::TwoDLinearSolverSplineMass, rhs::Vector{Float64})::Vector{Float64}
    nx1, nx2 = solver.nx1, solver.nx2

    solver.wk .= reshape(rhs, nx1, nx2)

    fft!(solver.wk)

    for j in 1:nx2, i in 1:nx1
        @inbounds solver.wk[i, j] /= (
            solver.eigvals1[i] * solver.eigvals2[j] * solver.factor
        )
    end

    ifft!(solver.wk)

    return vec(real(solver.wk))
end