text:
# QCP example
#    maximize x + y
#
#    s.t.  x, y >= 0
#          x^2 + y^2 <= 1
#
#    solution: (0.71, 0.71) objv = 1.414

@testset "QCQP 01" begin
    env = CPLEX.Env()

    model = CPLEX.Model(env, "qcqp_01")
    CPLEX.set_sense!(model, :Max)

    CPLEX.add_vars!(model, [1., 1.], 0., Inf)

     # add_qpterms!(model, linearindices, linearcoeffs, qrowinds, qcolinds, qcoeffs, sense, rhs)
    CPLEX.add_qconstr!(model, [], [], [1, 2], [1, 2], [1, 1.], '<', 1.0)

    CPLEX.optimize!(model)

    sol = CPLEX.get_solution(model)
    @test isapprox(sol[1], sqrt(2)/2, rtol=1e-4)
    @test isapprox(sol[2], sqrt(2)/2, rtol=1e-4)
    @test isapprox(CPLEX.get_objval(model), sqrt(2), rtol=1e-8)
end


text:
module LightPropagation

using DelimitedFiles
using LsqFit
using DSP
using FFTW
using Parameters

export TPSF_DA_semiinf_refl
export TPSF_DA_slab_refl
export TPSF_DA_paralpip_refl

export getfit

export load_asc_data

export fitDTOF, DTOF_fitparams

include("forwardmodels/TPSF/DiffusionApproximation/refl_TD.jl")
include("data/loadData/loaddata.jl")

include("inversefitting/FitStructures/fitstructures.jl")
include("inversefitting/FitFunctions/fitfunction.jl")

end

ext: jl

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation


text:
using Dates: Period, DatePeriod, TimePeriod

function Base.floor(zdt::ZonedDateTime, p::DatePeriod)
    return ZonedDateTime(floor(localtime(zdt), p), timezone(zdt))
end

function Base.floor(zdt::ZonedDateTime, p::TimePeriod)
    # Rounding is done using the current fixed offset to avoid transitional ambiguities.
    dt = floor(localtime(zdt), p)
    utc_dt = dt - zdt.zone.offset
    return ZonedDateTime(utc_dt, timezone(zdt); from_utc=true)
end

function Base.ceil(zdt::ZonedDateTime, p::DatePeriod)
    return ZonedDateTime(ceil(localtime(zdt), p), timezone(zdt))
end

#function Dates.floorceil(zdt::ZonedDateTime, p::Dates.DatePeriod)
    #return floor(zdt, p), ceil(zdt, p)
#end

"""
    floor(zdt::ZonedDateTime, p::Period) -> ZonedDateTime
    floor(zdt::ZonedDateTime, p::Type{Period}) -> ZonedDateTime

Returns the nearest `ZonedDateTime` less than or equal to `zdt` at resolution `p`. The
result will be in the same time zone as `zdt`.

For convenience, `p` may be a type instead of a value: `floor(zdt, Dates.Hour)` is a
shortcut for `floor(zdt, Dates.Hour(1))`.

`VariableTimeZone` transitions are handled as for `round`.

### Examples

The `America/Winnipeg` time zone transitioned from Central Standard Time (UTC-6:00) to
Central Daylight Time (UTC-5:00) on 2016-03-13, moving directly from 01:59:59 to 03:00:00.

```julia
julia> zdt = ZonedDateTime(2016, 3, 13, 1, 45, TimeZone("America/Winnipeg"))
2016-03-13T01:45:00-06:00

julia> floor(zdt, Dates.Day)
2016-03-13T00:00:00-06:00

julia> floor(zdt, Dates.Hour)
2016-03-13T01:00:00-06:00
```
"""
Base.floor(::TimeZones.ZonedDateTime, ::Union{Period, Type{Period}})

"""
    ceil(zdt::ZonedDateTime, p::Period) -> ZonedDateTime
    ceil(zdt::ZonedDateTime, p::Type{Period}) -> ZonedDateTime

Returns the nearest `ZonedDateTime` greater than or equal to `zdt` at resolution `p`.
The result will be in the same time zone as `zdt`.

For convenience, `p` may be a type instead of a value: `ceil(zdt, Dates.Hour)` is a
shortcut for `ceil(zdt, Dates.Hour(1))`.

`VariableTimeZone` transitions are handled as for `round`.

### Examples

The `America/Winnipeg` time zone transitioned from Central Standard Time (UTC-6:00) to
Central Daylight Time (UTC-5:00) on 2016-03-13, moving directly from 01:59:59 to 03:00:00.

```julia
julia> zdt = ZonedDateTime(2016, 3, 13, 1, 45, TimeZone("America/Winnipeg"))
2016-03-13T01:45:00-06:00

julia> ceil(zdt, Dates.Day)
2016-03-14T00:00:00-05:00

julia> ceil(zdt, Dates.Hour)
2016-03-13T03:00:00-05:00
```
"""
Base.ceil(::TimeZones.ZonedDateTime, ::Union{Period, Type{Period}})

"""
    round(zdt::ZonedDateTime, p::Period, [r::RoundingMode]) -> ZonedDateTime
    round(zdt::ZonedDateTime, p::Type{Period}, [r::RoundingMode]) -> ZonedDateTime

Returns the `ZonedDateTime` nearest to `zdt` at resolution `p`. The result will be in the
same time zone as `zdt`. By default (`RoundNearestTiesUp`), ties (e.g., rounding 9:30 to the
nearest hour) will be rounded up.

For convenience, `p` may be a type instead of a value: `round(zdt, Dates.Hour)` is a
shortcut for `round(zdt, Dates.Hour(1))`.

Valid rounding modes for `round(::TimeType, ::Period, ::RoundingMode)` are
`RoundNearestTiesUp` (default), `RoundDown` (`floor`), and `RoundUp` (`ceil`).

### `VariableTimeZone` Transitions

Instead of performing rounding operations on the `ZonedDateTime`'s internal UTC `DateTime`,
which would be computationally less expensive, rounding is done in the local time zone.
This ensures that rounding behaves as expected and is maximally meaningful.

If rounding were done in UTC, consider how rounding to the nearest day would be resolved for
non-UTC time zones: the result would be 00:00 UTC, which wouldn't be midnight local time.
Similarly, when rounding to the nearest hour in `Australia/Eucla (UTC+08:45)`, the result
wouldn't be on the hour in the local time zone.

When `p` is a `DatePeriod` rounding is done in the local time zone in a straightforward
fashion. When `p` is a `TimePeriod` the likelihood of encountering an ambiguous or
non-existent time (due to daylight saving time transitions) is increased. To resolve this
issue, rounding a `ZonedDateTime` with a `VariableTimeZone` to a `TimePeriod` uses the
`DateTime` value in the appropriate `FixedTimeZone`, then reconverts it to a `ZonedDateTime`
in the appropriate `VariableTimeZone` afterward.

Rounding is not an entirely "safe" operation for `ZonedDateTime`s, as in some cases
historical transitions for some time zones (such as `Asia/Colombo`) occur at midnight. In
such cases rounding to a `DatePeriod` may still result in an `AmbiguousTimeError` or a
`NonExistentTimeError`. (But these events should be relatively rare.)

Regular daylight saving time transitions should be safe.

### Examples

The `America/Winnipeg` time zone transitioned from Central Standard Time (UTC-6:00) to
Central Daylight Time (UTC-5:00) on 2016-03-13, moving directly from 01:59:59 to 03:00:00.

```julia
julia> zdt = ZonedDateTime(2016, 3, 13, 1, 45, TimeZone("America/Winnipeg"))
2016-03-13T01:45:00-06:00

julia> round(zdt, Dates.Hour)
2016-03-13T03:00:00-05:00

julia> round(zdt, Dates.Day)
2016-03-13T00:00:00-06:00
```

The `Asia/Colombo` time zone revised the definition of Lanka Time from UTC+6:30 to UTC+6:00
on 1996-10-26, moving from 00:29:59 back to 00:00:00.

```julia
julia> zdt = ZonedDateTime(1996, 10, 25, 23, 45, TimeZone("Asia/Colombo"))
1996-10-25T23:45:00+06:30

julia> round(zdt, Dates.Hour)
1996-10-26T00:00:00+06:30

julia> round(zdt, Dates.Day)
ERROR: Local DateTime 1996-10-26T00:00:00 is ambiguous
```
"""     # Defined in base/dates/rounding.jl
Base.round(::TimeZones.ZonedDateTime, ::Union{Period, Type{Period}})


text:
#=
Rational Krylov algorithm as described in:
"Rational Krylov approximation of matrix functions: Numerical methods
and optimal pole selection" (Stefan Güttel)
The function is ported from the Rational Krylov Toolbox for Matlab version 1.0
(http://guettel.com/rktoolbox/), with adaptations to make it type-stable.
Later versions of the toolbox add features that are unnecessary for our purposes.
=#
function ratkrylov(A::Mat, b::Vector{N}, p::Vector{Complex128}) where {
	N<:Union{Float64, Complex128}, Mat<:Union{Matrix{N}, SparseMatrixCSC{N}}}

	if size(A, 1) != size(A, 2)
		throw(DimensionMismatch("A must be square"))
	end
	if size(A, 2) != length(b)
		throw(DimensionMismatch("can't do A*b"))
	end
	m, n = length(p), length(b)
	if n == 0
		error("A is empty")
	end
	if m == 0
		error("no poles")
	end
	if m > n
		error("too many poles")
	end

	B = eye(A)
	V = zeros(N, n, m+1)
	K = zeros(N, m+1, m)
	H = zeros(K)
	realopt = N <: Real

	if realopt && !canonical_cplx(p)
		error("Cannot use real arithmetic if the poles are not ordered canonically.")
	end

	p, mu, rho, eta = poles_to_moebius(p)

	# run_krylov:
	bd = false
	bd_tol = eps()
	# Starting vector.
	V[:, 1] = b/norm(b)
	j = 1
	while j <= m
		if realopt && !isreal(p[j])
			# Computing the continuation combination.
			u = ones(Complex128, 1) # U[1:j, j] in the original code
			if j > 1
				Q = qr(K[1:j, 1:j-1]/mu[j] - H[1:j, 1:j-1]/p[j], thin=false)[1]
				u = Q[:, end]
			end

			# Compute new vector.
			w = V[:, 1:j]*u
			w = rho[j]*(A*w) - eta[j]*w
			w = (B/mu[j] - A/p[j]) \ w

			# Orthogonalization.
			V[:, j+1] = real(w)
			V[:, j+2] = imag(w)
			# MGS
			for j = j:j+1
				for reo = 0:1
					for reo_i = 1:j
						vi = view(V, 1:size(V, 1), reo_i)
						vj = view(V, 1:size(V, 1), j+1)
						hh = vi'*vj
						V[:, j+1] = vj - vi*hh
						H[reo_i, j] += hh
					end
				end
				normw = norm(V[:, j+1])
				H[j+1, j] = normw
				V[:, j+1] /= normw
				if normw < bd_tol*norm(H[1:j, j])
					bd = true
					break
				end
			end

			# Setting the decomposition.
			rp, ip = real(1/p[j-1]), imag(1/p[j-1])
			cp = [rp ip; -ip rp]
			u0, h = [real(u) imag(u); 0 0; 0 0], H[1:j+1, j-1:j]
			K[1:j+1, j-1:j] = rho[j-1]*u0 + h*cp
			H[1:j+1, j-1:j] = eta[j-1]*u0 + h/mu[j-1]
		else
			pj = N(p[j])
			# Computing the continuation combination.
			u = ones(N, 1) # U[1:j, j] in the original code
			if j > 1
				Q = qr(K[1:j, 1:j-1]/mu[j] - H[1:j, 1:j-1]/pj, thin=false)[1]
				u = Q[:, end]
			end

			# Compute new vector.
			w = V[:, 1:j]*u
			w = rho[j]*(A*w) - eta[j]*w
			w = (B/mu[j] - A/pj) \ w

			# Orthogonalization, MGS.
			for reo = 0:1
				for reo_i = 1:j
					v = view(V, 1:size(V, 1), reo_i)
					hh = v'*w
					w -= v*hh
					H[reo_i, j] += hh
				end
			end
			normw = norm(w)
			H[j+1, j] = normw
			V[:, j+1] = w/normw
			if normw < bd_tol*norm(H[1:j, j])
				bd = true
				break
			end

			# Setting the decomposition.
			u0, h = [u; 0], H[1:j+1, j]
			K[1:j+1, j] = rho[j]*u0 + h/pj
			H[1:j+1, j] = eta[j]*u0 + h/mu[j]
		end # realopt

		j = j+1
	end # while j <= m

	if bd == true
		println("lucky breakdown")
		V = V[:, 1:j]
		K = K[1:j, 1:j-1]
		H = H[1:j, 1:j-1]
	end
	return V, K, H
end

#=
canonical_cplx checks if the poles p are ordered canonically.
=#
function canonical_cplx(p::Vector{Complex128})
	m = length(p)
	j = 1
	while j <= m
		if isreal(p[j]) || isinf(p[j])
			j = j+1
		else
			if j == m || p[j+1] != conj(p[j])
				return false
			end
			j = j+2
		end
	end
	return true
end

#=
Moebius transformation with poles p.
p != 0 is replaced by (p, mu) := (p,   1) and (rho, eta) := (1, 0),
p  = 0 is replaced by (p, mu) := (1, Inf) and (rho, eta) := (0, 1).
=#
function poles_to_moebius(p::Vector{Complex128})
	p = copy(p)
	mu = ones(length(p))
	rho = ones(length(p))
	eta = zeros(length(p))

	sel = p .== 0

	p[sel] = 1
	mu[sel] = Inf
	rho[sel] = 0
	eta[sel] = 1

	return p, mu, rho, eta
end

#=
Computes f(A)*b using the rational Krylov approximation with poles p.
=#
function ratkrylov(f::Func, A::Mat, b::Vector{N}, p::Vector{Complex128}) where {
	Func, N<:Union{Float64, Complex128}, Mat<:Union{Matrix{N}, SparseMatrixCSC{N}}}

	V, K, H = ratkrylov(A, b, p)

	m = size(V, 2) - 1 # may be < length(p) in case of breakdown
	Am = isinf(p[m]) ? [H/K[1:m,1:m] V'*(A*V[:,end])] : V'*A*V

	return V*(schurparlett(f, Am)*(V'*b))
end

#=
Computes f(A)*b using rational Krylov.
The poles of f are found automatically by AAA, with parameters mmax, tol and Z.
You typically want to manually set at least mmax, to limit the size of the resulting Krylov space.
If the sample set Z is not provided, f will be sampled on the 0-centered disk with radius
min(norm(A, 1), norm(A, Inf), vecnorm(A)).
When A is Real, the algorithm will assume f(conj(x)) == conj(f(x)).
=#
function ratkrylov(f::Func, A::Mat, b::Vector{N},
	mmax::Int64=100, tol::Float64=1e-13, Z::Vector{M}=Vector{Complex128}(0)) where {
	Func, M<:Union{Float64, Complex128}, N<:Union{Float64, Complex128},
	Mat<:Union{Matrix{N}, SparseMatrixCSC{N}}}

	if length(Z) == 0
		rad = min(norm(A, 1), norm(A, Inf), vecnorm(A))
		rad = max(rad, sqrt(eps()))
		nsamples = Mat<:SparseMatrixCSC ? nnz(A)÷2 : prod(size(A))
		nsamples = max(nsamples, 100)
		Z = lpdisk(rad, nsamples)
	end

	pol = Vector{Complex128}(aaa(f, Z, tol, mmax)[2])
	if N <: Real
		# Real ratkrylov wants conjugated and canonically ordered poles.
		pos = pol[imag.(pol) .> 1e-13]
		pol = pol[abs.(imag.(pol)) .<= 1e-13]
		pol[1:end] = real.(pol)
		for i = 1:length(pos)
			pol = [pol; pos[i]; conj(pos[i])]
		end
	end
	# Inf as last pole allows for faster projection of A in the Krylov space:
	pol = [pol; Inf]

	println(length(pol))

	return ratkrylov(f, A, b, pol)
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

import .Serializer: known_object_data, object_number, serialize_cycle, deserialize_cycle, writetag,
                      __deserialized_types__, serialize_typename, deserialize_typename,
                      TYPENAME_TAG, object_numbers

type ClusterSerializer{I<:IO} <: AbstractSerializer
    io::I
    counter::Int
    table::ObjectIdDict

    pid::Int                                     # Worker we are connected to.
    tn_obj_sent::Set{UInt64}                     # TypeName objects sent
    glbs_sent::Dict{UInt64, UInt64}              # (key,value) -> (object_id, hash_value)
    glbs_in_tnobj::Dict{UInt64, Vector{Symbol}}  # Track globals referenced in
                                                 # anonymous functions.
    anonfunc_id::UInt64

    ClusterSerializer(io::I) = new(io, 0, ObjectIdDict(), Base.worker_id_from_socket(io),
                                Set{UInt64}(), Dict{UInt64, UInt64}(), Dict{UInt64, Vector{Symbol}}(), 0)
end
ClusterSerializer(io::IO) = ClusterSerializer{typeof(io)}(io)

function deserialize(s::ClusterSerializer, ::Type{TypeName})
    full_body_sent = deserialize(s)
    number = read(s.io, UInt64)
    if !full_body_sent
        tn = get(known_object_data, number, nothing)::TypeName
        if !haskey(object_numbers, tn)
            # set up reverse mapping for serialize
            object_numbers[tn] = number
        end
        deserialize_cycle(s, tn)
    else
        tn = deserialize_typename(s, number)
    end

    # retrieve arrays of global syms sent if any and deserialize them all.
    foreach(sym->deserialize_global_from_main(s, sym), deserialize(s))
    return tn
end

function serialize(s::ClusterSerializer, t::TypeName)
    serialize_cycle(s, t) && return
    writetag(s.io, TYPENAME_TAG)

    identifier = object_number(t)
    send_whole = !(identifier in s.tn_obj_sent)
    serialize(s, send_whole)
    write(s.io, identifier)
    if send_whole
        # Track globals referenced in this anonymous function.
        # This information is used to resend modified globals when we
        # only send the identifier.
        prev = s.anonfunc_id
        s.anonfunc_id = identifier
        serialize_typename(s, t)
        s.anonfunc_id = prev
        push!(s.tn_obj_sent, identifier)
        finalizer(t, x->cleanup_tname_glbs(s, identifier))
    end

    # Send global refs if required.
    syms = syms_2b_sent(s, identifier)
    serialize(s, syms)
    foreach(sym->serialize_global_from_main(s, sym), syms)
    nothing
end

function serialize(s::ClusterSerializer, g::GlobalRef)
    # Record if required and then invoke the default GlobalRef serializer.
    sym = g.name
    if g.mod === Main && isdefined(g.mod, sym)
        v = getfield(Main, sym)
        if  !isa(v, DataType) && !isa(v, Module) &&
            (binding_module(Main, sym) === Main) && (s.anonfunc_id != 0)
            push!(get!(s.glbs_in_tnobj, s.anonfunc_id, []), sym)
        end
    end

    invoke(serialize, Tuple{AbstractSerializer, GlobalRef}, s, g)
end

# Send/resend a global object if
# a) has not been sent previously, i.e., we are seeing this object_id for the first time, or,
# b) hash value has changed or
# c) is a bitstype
function syms_2b_sent(s::ClusterSerializer, identifier)
    lst = Symbol[]
    check_syms = get(s.glbs_in_tnobj, identifier, [])
    for sym in check_syms
        v = getfield(Main, sym)

        if isbits(v)
            push!(lst, sym)
        else
            oid = object_id(v)
            if haskey(s.glbs_sent, oid)
                # We have sent this object before, see if it has changed.
                s.glbs_sent[oid] != hash(v) && push!(lst, sym)
            else
                push!(lst, sym)
            end
        end
    end
    return unique(lst)
end

function serialize_global_from_main(s::ClusterSerializer, sym)
    v = getfield(Main, sym)

    oid = object_id(v)
    record_v = true
    if isbits(v)
        record_v = false
    elseif !haskey(s.glbs_sent, oid)
        # set up a finalizer the first time this object is sent
        try
            finalizer(v, x->delete_global_tracker(s,x))
        catch ex
            # Do not track objects that cannot be finalized.
            if isa(ex, ErrorException)
                record_v = false
            else
                rethrow(ex)
            end
        end
    end
    record_v && (s.glbs_sent[oid] = hash(v))

    serialize(s, isconst(Main, sym))
    serialize(s, v)
end

function deserialize_global_from_main(s::ClusterSerializer, sym)
    sym_isconst = deserialize(s)
    v = deserialize(s)
    if sym_isconst
        @eval Main const $sym = $v
    else
        @eval Main $sym = $v
    end
end

function delete_global_tracker(s::ClusterSerializer, v)
    oid = object_id(v)
    if haskey(s.glbs_sent, oid)
        delete!(s.glbs_sent, oid)
    end

    # TODO: A global binding is released and gc'ed here but it continues
    # to occupy memory on the remote node. Would be nice to release memory
    # if possible.
end

function cleanup_tname_glbs(s::ClusterSerializer, identifier)
    delete!(s.glbs_in_tnobj, identifier)
end

# TODO: cleanup from s.tn_obj_sent

max_stars_repo_name: oxinabox/julia

max_issues_repo_name: oxinabox/julia

max_forks_repo_name: oxinabox/julia


text:
using Pkg
Pkg.add("LibPQ")
Pkg.add("JSON")
Pkg.add("DataStreams")
Pkg.add("JLD")
Pkg.add("Images")
Pkg.add("FileIO")
Pkg.add("ArgParse")
Pkg.add("YAML")
Pkg.add("StaticArrays")


text:
using ClusteringMetrics
using Test

Z1 = fill(0.0, 2, 10)
Z1[:,1:5] = [0.1 0.2 0.15 0.2 0.22;
             0.15 0.21 0.13 0.21 0.19]
Z1[:,6:end] = Z1[:,1:5] .- [1.0;1.0]
ϵ = [-0.010114902986680585 0.004463792650270601 -0.00153613750975559 -0.0055050913393961485 0.0035611946194614437 -0.0009473083769016829 -0.0011402700190551592 0.00043139456531041 0.007363905500805968 -0.0036852874409562726;
     -0.004970181279034265 -0.008439692559205151 -0.002063148403333957 0.0027397663153164308 -0.009655682877822354 0.002952749805508587 0.010722686475465308 -0.008101446035381798 0.0024392528878517524 -0.007034803358359168]
Z1 .+= ϵ
label1 = [fill(1,5);fill(2,5)]

@testset "Isolation distance" begin
    # large isolation distance
    d = ClusteringMetrics.get_isolation_distance(Z1, label1)
    @test d ≈  [7346.57515900546, 11695.14306811699]
end

@testset "Cluster matching" begin
    Z2 = fill!(similar(Z1), 0.0)
    Z2 .= Z1
    label2 = label1
    dd = ClusteringMetrics.get_isolation_distance(Z1, label1, Z2, label2)
    # large cross distance, indicating high generalisability
    @test dd ≈ [296.9888641641518, 301.87857155988695]
    dd = ClusteringMetrics.get_cluster_distance(Z1, label1, Z2, label2)

    # zero confusion when the sets are the same
    @test dd ≈ [0.0, 0.0]

    n = ClusteringMetrics.get_cluster_isolation(Z1, label1, Z2, label2)
    @test n ≈ [1.0, 1.0]

    Z2[:,1:5] = Z1[:,6:end]
    Z2[:, 6:end] = Z1[:, 1:5]
    dd = ClusteringMetrics.get_isolation_distance(Z1, label1, Z2, label2)
    # low cross distance indicating low generlisability
    @test dd ≈ [0.0033671296154972326, 0.003312590207488839]
    dd = ClusteringMetrics.get_cluster_distance(Z1, label1, Z2, label2)
    @test dd ≈ [64.50554889037353, 64.5055488903735]
    n = ClusteringMetrics.get_cluster_isolation(Z1, label1, Z2, label2)
    @test n ≈ [0.0, 0.0]
end


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2020
#
# This is a Julia package for discrete event simulation
#
# -------------------------------------------
#
# This is the first introductory pet example.
# We simulate the easy life of a pet in the morning with a process-based simulation.
#

using DiscreteEvents, Printf, Random

# again we need some data structures
abstract type PetState end
struct Sleeping  <: PetState end
struct Scuffing  <: PetState end
struct Running   <: PetState end

mutable struct Pet
    clk::Clock
    name::String
    state::PetState
    speak::String
end

# the pet API is now different, see documentation
speak(p, n) = now!(p.clk, fun(println, @sprintf("%5.2f %s: %s", tau(p.clk), p.name, p.speak^n)))

function setstate!(p::Pet, q::PetState)
    p.state = q
    now!(p.clk, fun(println, @sprintf("%5.2f %s: %s", tau(p.clk), p.name, repr(p.state))))
end

# the pet function describes the sequence of states a pet lives through
function pet(clk::Clock, p::Pet)
    setstate!(p, Running());  delay!(clk,  5*rand())
    speak(p, 5);              delay!(clk,    rand())  # get hungry
    setstate!(p, Scuffing()); delay!(clk,  2*rand())
    setstate!(p, Running());  delay!(clk,  2*rand())
    speak(p, 2);              delay!(clk,  2*rand())  # get weary
    setstate!(p, Sleeping()); delay!(clk, 10*rand())
end

# we set it all up and run the simulation
clk = Clock(t0=5)
Random.seed!(123)
snoopy = Pet(clk, "Snoopy", Sleeping(), "huff")
process!(clk, Prc(1, pet, snoopy));
run!(clk, 20)


text:
using FSM

cn = Constants{Float64}()

ebm = EBM{Float64}(
    am=1,
    cm=1,
    dm=1,
    em=1,
    hm=1,
)

@code_warntype soil(ebm)

time_soil(ebm) = @time soil(ebm)

time_soil(ebm)
time_soil(ebm)


text:
import CUDAdrv: CUarray_format

const _type_to_cuarrayformat_dict = Dict{DataType,CUarray_format}(UInt8 => CUDAdrv.CU_AD_FORMAT_UNSIGNED_INT8,
    UInt16 => CUDAdrv.CU_AD_FORMAT_UNSIGNED_INT16,
    UInt32 => CUDAdrv.CU_AD_FORMAT_UNSIGNED_INT32,
    Int8 => CUDAdrv.CU_AD_FORMAT_SIGNED_INT8,
    Int16 => CUDAdrv.CU_AD_FORMAT_SIGNED_INT16,
    Int32 => CUDAdrv.CU_AD_FORMAT_SIGNED_INT32,
    Float16 => CUDAdrv.CU_AD_FORMAT_HALF,
    Float32 => CUDAdrv.CU_AD_FORMAT_FLOAT,
)

@inline _alias_type_to_nchan_and_format(::Type{T}) where {T} = @error "Type `$T` is not a valid alias type for a \"CUDA array\" (texture memory) format"
for (T, cuarrayformat) in _type_to_cuarrayformat_dict
    @eval @inline _alias_type_to_nchan_and_format(::Type{$T})  = 1, $cuarrayformat, $T
    @eval @inline _alias_type_to_nchan_and_format(::Type{NTuple{2,$T}}) = 2, $cuarrayformat, $T
    @eval @inline _alias_type_to_nchan_and_format(::Type{NTuple{4,$T}}) = 4, $cuarrayformat, $T
end

@inline function _assert_alias_size(::Type{T}, ::Type{Ta}) where {T,Ta}
    @assert sizeof(T) == sizeof(Ta) "Error: Julia type `$T` cannot be aliased to the \"CUDA array\" (texture memory) format `$Ta`: sizes in bytes do not match"
end

for (T, _) in _type_to_cuarrayformat_dict
    @eval @inline cuda_texture_alias_type(::Type{$T}) = $T
    for N in (2, 4)
        @eval @inline cuda_texture_alias_type(::Type{NT}) where {NT <: NTuple{$N,$T}} = NT
    end
end

@generated function cuda_texture_alias_type(t::Type{T}) where T
    err = "An alias from the type `$T` to a \"CUDA array\" (texture memory) format could not been inferred."
    isprimitivetype(T) && return :(@error $("Primitive type `$T` does not have a defined alias to a \"CUDA array\" (texture memory) format."))
    isbitstype(T) || return :(@error $(err * " Type in not `isbitstype`."))
    Te = nothing
    N = 0
    datatypes = DataType[T]
    while !isempty(datatypes)
        for Ti in fieldtypes(pop!(datatypes))
            if isprimitivetype(Ti)
                Ti = cuda_texture_alias_type(Ti)
                typeof(Ti) == DataType || return :(@error $(err * " Composed of primitive type with no alias."))
                if !isprimitivetype(Ti)
                    push!(datatypes, Ti)
                    break
                end
                if Te == nothing
                    Te = Ti
                    N = 1
                elseif Te == Ti
                    N += 1
                else
                    return :(@error $(err * " Incompatible elements."))
                end
            else
                push!(datatypes, Ti)
            end
        end
    end
    Ta = NTuple{N,Te}
    sizeof(T) == sizeof(Ta) || return :(@error $(err * " Inferred alias type and original type have different byte sizes."))
    in(N, (1, 2, 4)) || return :(@error $(err * " Incompatible number of elements ($N, but only 1, 2 or 4 supported)."))
    return Ta
end


text:
module DocumentationOverviewTests

include("test_pkgapi.jl")
include("test_find.jl")
include("test_table.jl")
include("test_list.jl")

end  # module DocumentationOverviewTests


text:
using FastLevenshtein
using Documenter

makedocs(;
    modules=[FastLevenshtein],
    authors="Matt Kearns mattkearns@protonmail.com",
    repo="https://github.com/mattkearns/FastLevenshtein.jl/blob/{commit}{path}#L{line}",
    sitename="FastLevenshtein.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://mattkearns.github.io/FastLevenshtein.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/mattkearns/FastLevenshtein.jl",
)


text:
function value_of(::Type{MeasurementType}, arg0::JString)
    return jcall(MeasurementType, "valueOf", MeasurementType, (JString,), arg0)
end

function values(::Type{MeasurementType})
    return jcall(MeasurementType, "values", Vector{MeasurementType}, ())
end


text:
# Code to accelerate gpu_batch but temporarily omitted for simplicity

macro merge_swap_shfl()
    q = quote
    b = shfl_sync(0xFFFFFFFF, sum, batch_floor(idx, 2 * n))
    a = n - b
    d = shfl_sync(0xFFFFFFFF, sum, batch_ceil(idx, 2 * n))
    c = n - d
    swap = value
    sum = d + b
    value = shfl_sync(0xFFFFFFFF, swap, batch_step_swap(idx, 2 * n, a, b, c))
    n *= 2
    end
    return esc(q)
end

macro repeat_inline(N, expr)
    function tail(expr, N)
        N <= 0 ? :() : :($expr; $(tail(expr, N - 1)))
    end
    out = tail(expr, N)
    return esc(out)
end

function shfl_reduce_kernel(src, dest, flags)
    i = threadIdx().x
    j = blockIdx().x
    w = blockDim().x
    floor = w * (j - 1)
    idx = floor + i
    @inbounds value = src[idx]

    @inbounds sum = 1 & flags[idx]
    n = 1
    # five merges => 2 ^ 5 = 32 = size of correct batches
    @repeat_inline 5 @merge_swap_shfl

    @inbounds dest[idx] = value
    return nothing
end

function shfl_reduce_kernel(A, flags)
    i = threadIdx().x
    j = blockIdx().x
    w = blockDim().x
    floor = w * (j - 1)
    idx = floor + i
    value = A[idx]

    sum = 1 & flags[idx]
    n = 1
    # five merges => 2 ^ 5 = 32 = size of correct batches
    @repeat_inline 5 @merge_swap_shfl

    A[idx] = value
    return nothing
end


text:
@testset "test_methods" begin
    @testset "Kalman Filter methods" begin
            x₀ = [1;0]
            F = [1. 1; 0 -1]
            H = [1 0]
            P₀ = 3. *I(2)
            Q = 1. *I(2)
            R = 2. *I(1)
            y = 2
            kf = KalmanFilter(x₀, F, H, P₀, Q, R)
            update!(kf, y)
            P1exp = [26/5 -3; -3 4]
            @test kf.P ≈ P1exp
            K₀ = Matrix{Float64}(undef,2,1)
            K₀[:,:] = [3/5;0]
            xexpected = F*x₀ + K₀*([y]-H*x₀)
            xtrue = predict(kf)
            @test xtrue == xexpected
            pdy = 1 / (2*π) / sqrt(5) * exp(-1 / 2 * 1 / 5)
            @test kf.pdy[] == pdy
    end

    @testset "Minimax MMAE methods" begin
            x₀ = [1;0]
            F = [1. 1; 0 -1]
            H = [1 0]
            P₀ = 3. *I(2)
            Q = 1. *I(2)
            R = 2. *I(1)
            y = 2
            kf1 = KalmanFilter(x₀, F, H, P₀, Q, R)
            kf2 = KalmanFilter(x₀, -F, H, P₀, Q, R)
            γ=10
            K₀ = Matrix{Float64}(undef,2,1)
            K₀[:,:] = [3/5;0]
            x1 = F*x₀ + K₀*([y]-H*x₀)
            x2 = -F*x₀ + -K₀*([y]-H*x₀)
            mini = MinimaxMMAE([kf1, kf2], () ->  Hypatia.Optimizer())
            update!(mini, y)
            @test predict(mini.filterbank[1]) == x1
            @test predict(mini.filterbank[2]) == x2
            @test mini.filterbank[1].c[] == 1/5

            yhatreal, valreal = predict(mini, γ)
            P1exp = [26/5 -3; -3 4]
            valexp = x1'*H'*pinv(I - γ^(-2)*H*P1exp*H')*H*x1 - γ^2*kf1.c[]
            @test yhatreal[1] ≈ 0 atol=1e-12
            @test valexp ≈ valreal
    end

    @testset "Bayesian MMAE methods" begin
        x₀ = [1;0]
        F = [1. 1; 0 -1]
        H = [1 0]
        P₀ = 3. *I(2)
        Q = 1. *I(2)
        R = 2. *I(1)
        y = 2
        kf1 = KalmanFilter(x₀, F, H, P₀, Q, R)
        kf2 = KalmanFilter(x₀, -F, H, P₀, Q, R)
        γ=10
        K₀ = Matrix{Float64}(undef,2,1)
        K₀[:,:] = [3/5;0]
        x1 = F*x₀ + K₀*([y]-H*x₀)
        x2 = -F*x₀ + -K₀*([y]-H*x₀)
        bayesian = BayesianMMAE([kf1, kf2])
        update!(bayesian, y)
        @test predict(bayesian.filterbank[1]) == x1
        @test predict(bayesian.filterbank[2]) == x2
        @test predict(bayesian) == [0]
         # Filters have same covariance and estimation error -> equal probability.
        @test bayesian.filterbank[1].pdy[] == 1/2
    end
end


text:
module FunctionTabulations
using Interpolations
using JLD2
using ProgressMeter
using Unitful
using Unitful: FreeUnits
using Base.Threads
using SHA: sha256

export create_tabulation_1D, create_tabulation_2D, create_tabulation_3D

"""
    scaler(x::Real, scale::Symbol)

Utility function to scale a variable, accepts :linear and :log10 scale
"""
function scaler(x::Real, scale::Symbol)
    if scale == :linear
        return x
    elseif scale == :log10
        return log10(x)
    else
        throw(ArgumentError("Scale <$scale> is not supported"))
    end
end

"""
    un_scaler(x::Real, scale::Symbol)

Utility function to un_scale a variable, accepts :linear and :log10 scale
"""
function un_scaler(x::Real, scale::Symbol)
    if scale == :linear
        return x
    elseif scale == :log10
        return 10 .^ x
    else
        throw(ArgumentError("Scale <$scale> is not supported"))
    end
end

"""
    compute_SHA(func)

Utility function to compute the SHA of a function.
"""
function compute_SHA(func::Function)
    c1 = code_lowered(func)
    h1 = bytes2hex(sha256(string(c1)))
    return h1
end

"""
    test_sha(func, sha[; mode = :warn])

Compares the SHA of the `func` funtion with the provided SHA.
"""
function test_sha(func::Function, sha::String; mode::Symbol = :warn)
    func_sha = compute_SHA(func)
    func_name = nameof(func)
    if sha != func_sha
        if mode == :warn
            @warn "The SHA for `$func_name` did not match the one of the stored tabulated function. Please check if the function definition has changed."
        elseif mode == :throw
            @error "The SHA for `$func_name` did not match the one of the stored tabulated function. Please check if the function definition has changed."
        elseif mode == :none
        else
            @error "Mode not supported. Got $(mode)."
        end
    end
end

"""
    wrap_function_1D_add_units(func, x_units, f_units[, args...; kwargs...])

Add units back to a 1D function.
"""
function wrap_function_1D_add_units(func::Function, x_units, f_units, args...; kwargs...)
    function wrapped_func(x)
        x_v = ustrip(x_units, x)
        return func(x_v, args...; kwargs...) * f_units
    end
    return wrapped_func
end

"""
    wrap_function_2D_add_units(func, x_units, y_units, f_units[, args...; kwargs...])

Add units back to a 2D function.
"""
function wrap_function_2D_add_units(func::Function, x_units, y_units, f_units, args...; kwargs...)
    function wrapped_func(x, y)
        x_v = ustrip(x_units, x)
        y_v = ustrip(y_units, y)
        return func(x_v, y_v, args...; kwargs...) * f_units
    end
    return wrapped_func
end

"""
    wrap_function_3D_add_units(func, x_units, y_units, z_units, f_units[, args...; kwargs...])

Add units back to a 3D function.
"""
function wrap_function_3D_add_units(func::Function, x_units, y_units, z_units, f_units, args...; kwargs...)
    function wrapped_func(x, y, z)
        x_v = ustrip(x_units, x)
        y_v = ustrip(y_units, y)
        z_v = ustrip(z_units, z)
        return func(x_v, y_v, z_v, args...; kwargs...) * f_units
    end
    return wrapped_func
end

"""
    wrap_function_1D_remove_units(func, x_units, f_units[, args...; kwargs...])

Remove units from a 1D function.
"""
function wrap_function_1D_remove_units(func::Function, x_units, f_units, args...; kwargs...)
    function wrapped_function(x::Real)
        x_v = x * x_units
        return ustrip(f_units, func(x_v, args...; kwargs...))
    end
    return wrapped_function
end

"""
    wrap_function_2D_remove_units(func, x_units, y_units, f_units[, args...; kwargs...])

Remove units from a 2D function.
"""
function wrap_function_2D_remove_units(func::Function, x_units, y_units, f_units, args...; kwargs...)
    function wrapped_function(x::Real, y::Real)
        x_v = x * x_units
        y_v = y * y_units
        return ustrip(f_units, func(x_v, y_v, args...; kwargs...))
    end
    return wrapped_function
end

"""
    wrap_function_3D_remove_units(func, x_units, y_units, z_units, f_units[, args...; kwargs...])

Remove units from a 3D function.
"""
function wrap_function_3D_remove_units(func::Function, x_units, y_units, z_units, f_units, args...; kwargs...)
    function wrapped_function(x::Real, y::Real, z::Real)
        x_v = x * x_units
        y_v = y * y_units
        z_v = z * z_units
        return ustrip(f_units, func(x_v, y_v, z_v, args...; kwargs...))
    end
    return wrapped_function
end

# made with <3

"""
    create_tabulation_1D(
        func::Function[,
        args...];
        xmin::T,
        xmax::T,
        npoints::Int[,
        x_scale::Symbol = :linear,
        f_scale = :linear,
        jld_base_path = nothing,
        custom_name = nothing,
        interpolation_type = :linear,
        extrapolation_bc = Throw,
        check_SHA = true,
        check_SHA_mode = :warn,
        kwargs...]
    ) where {T<:Union{Real,Quantity}}

Computes or loads the tabulation for a given function of the kind `f(x)`.

# Arguments
- `func::Function`: the function which should be tabulated. Must be of the kind `f(x)`
- `xmin::Union{Real,Quantity}`: the minimum `x` for the range which will be covered by the tabulation
- `xmax::Union{Real,Quantity}`: the maximum `x` for the range which will be covered by the tabulation
- `npoints::Int`: the number of points for which `f(x)` will be computed
- `x_scale::Symbol`: the scale which will be used for tabulation. If `:linear` is choosen, the `x` points will be evenly spaced. If `:log10` is chosen, the points will be logarithmically spaced. Defaults to `:linear`
- `f_scale::Symbol`: the scale which will be used for tabulation. If `:linear` is choosen, the `f(x)` points will not be alterated. If `:log10` is chosen, log10 will be applied to the `f(x)` points before interpolating. Defaults to `:linear`
- `jld_base_path`: path to the folder where the tabulation should be saved. Defaults to the current folder
- `custom_name`: custom name for the tabulation file, to which `_data.jld2` will be appended. Defaults to the name of the function to be tabulated
- `interpolation_type`: Type of the spline to be used for interpolation. Can either be :linear (1st order spline) or :cubic (3rd order spline)
- `extrapolation_bc`: behaviour of the tabulation outside of the boundaries defined by `[xmin, xmax]`. Possible behaviours are handled by Interpolations.jl and include `Throw`` (throws and eror if a value out of bounds is acessed) or `Line` (extrapolate linearly). Defaults to `Throw`
- `check_SHA::Bool`: Whether to check the SHA of the tabulated function. If the SHA of the tabulated function does not match the SHA stored inside the tabulation file, it might be necessary to recompute the tabulation, since it's likely that the definition of the tabulated function has changed.
- `check_SHA_mode::Symbol`: Describes the behaviour of this function if the SHA of the tabulated function does not match the SHA stored inside the tabulation file. Defaults to :warn (print a warning but load the tabulation all the same). Other options include :throw (throw an error to suggest manual deletion of the old tabulation) and :none (don't do anything).
- `args..., kwargs...`: additional `args` and `kwargs` will be passed to the function which is to be tabulated

# Examples
```jldoctest
julia>  func_1d(x) = sin(x)

julia>  itp_1d_1 = create_tabulation_1D(
            func_1d,
            xmin = 0.0,
            xmax = 3.0,
            npoints = 100,
            x_scale = :linear,
            f_scale = :linear
        )

julia>  isapprox(itp_1d_1(2.0), func_1d(2.0), rtol = 1e-3)
true
```

Measurement units are supported too:

```jldoctest
julia>  using Unitful

julia>  func_1d(x) = x^2

julia>  itp_1d_1 = create_tabulation_1D(
            func_1d,
            xmin = 0.0u"m",
            xmax = 3.0u"m",
            npoints = 100,
            x_scale = :linear,
            f_scale = :linear
        )

julia>  isapprox(itp_1d_1(2.0u"m"), func_1d(2.0u"m"), rtol = 1e-3)
true
```
"""
function create_tabulation_1D(
    func::Function,
    args...;
    jld_base_path = nothing,
    custom_name = nothing,
    xmin::T,
    xmax::T,
    npoints::Int,
    x_scale::Symbol = :linear,
    f_scale = :linear,
    interpolation_type = :linear,
    extrapolation_bc = Throw,
    check_SHA::Bool = true,
    check_SHA_mode::Symbol = :warn,
    kwargs...
) where {T<:Union{Real,Quantity}}

    func_name = nameof(func)

    arg_tmp(x) = func(x, args...; kwargs...)
    x_units = unit(xmin)
    f_units = unit(arg_tmp(xmin))

    xmin_v = ustrip(x_units, xmin)::Real
    xmax_v = ustrip(x_units, xmax)::Real

    arg = wrap_function_1D_remove_units(arg_tmp, x_units, f_units)

    if isnothing(jld_base_path)
        base_path = pwd()
    else
        base_path = jld_base_path
    end

    if !(ispath(base_path))
        mkdir(base_path)
    end

    if isnothing(custom_name)
        filename = "$(func_name)_data.jld2"
    else
        filename = "$(custom_name)_data.jld2"
    end

    filepath = "$(base_path)/$(filename)"

    function load_file()
        if isfile(filepath)
            data = load(filepath)
            if data["xmin"] == xmin && data["xmax"] == xmax && data["npoints"] == npoints
                return true
            else
                @warn "The provided tabulation range or npoints does not match the data stored in the tabulation file. The tabulation will be recomputed."
                return false
            end
        else
            return false
        end
    end

    if load_file()
        data = load(filepath)

        @info "$(filename) loaded!"
        if check_SHA
            func_sha = data["sha"]
            test_sha(func, func_sha, mode = check_SHA_mode)
        end
        x = data["x"]

        data_matrix = data["func"]

    else
        if x_scale == :linear
            x = range(xmin_v, xmax_v, length = npoints)
        elseif x_scale == :log10
            x = range(log10(xmin_v), log10(xmax_v), length = npoints)
        else
            throw(ArgumentError("X scale $x_scale not supported"))
        end

        if isnothing(custom_name)
            @info "Computing $(func_name) Interpolation"
        else
            @info "Computing $(custom_name) Interpolation"
        end

        p = Progress(Int(npoints))
        update!(p, 0)

        data_matrix = zeros(npoints)

        @threads for i = 1:npoints

            data_matrix[i] = arg(un_scaler(x[i], x_scale))

            next!(p)
        end

        data_dict = Dict{
            String,
            Union{
                String,
                Real,
                Quantity,
                Array,
                StepRangeLen
            },
        }()
        sha = compute_SHA(func)
        data_dict["x"] = x
        data_dict["func"] = convert(Array, data_matrix)
        data_dict["sha"] = sha
        data_dict["xmin"] = xmin
        data_dict["xmax"] = xmax
        data_dict["npoints"] = npoints

        save(filepath, data_dict)
        @info "$(filename) created and exported!"
    end

    if f_scale == :log10
        data_matrix[data_matrix.<=1e-300] .= 1e-300
    end

    knots = (x,)
    f_matrix = scaler.(data_matrix, f_scale)

    if interpolation_type == :linear
        itp = LinearInterpolation(knots, f_matrix, extrapolation_bc = extrapolation_bc())
    elseif interpolation_type == :cubic
        itp = CubicSplineInterpolation(knots, f_matrix, extrapolation_bc = extrapolation_bc())
    else
        throw(ArgumentError("$interpolation_type is not a valid interpolation type"))
    end

    function func_interp(x::Real)
        return un_scaler(itp(scaler(x, x_scale)), f_scale)::Float64
    end

    func_interp_units = wrap_function_1D_add_units(func_interp, x_units, f_units)

    return func_interp_units::Function
end

"""
    create_tabulation_1D(
        func::Function,
        x::Vector{T},
        args...;
        jld_base_path = nothing,
        custom_name = nothing,
        x_scale::Symbol = :linear,
        f_scale = :linear,
        interpolation_type = :linear,
        extrapolation_bc = Throw,
        check_SHA = true,
        check_SHA_mode = :warn,
        kwargs...
        ) where {T<:Union{Real,Quantity}}

Provides an interface for `create_tabulation_1D` where the tabulation points are defined by the Vector `x`. Only `LinearInterpolation` is supported.`
"""
function create_tabulation_1D(
    func::Function,
    x::Vector{T},
    args...;
    jld_base_path = nothing,
    custom_name = nothing,
    x_scale::Symbol = :linear,
    f_scale = :linear,
    interpolation_type = :linear,
    extrapolation_bc = Throw,
    check_SHA::Bool = true,
    check_SHA_mode::Symbol = :warn,
    kwargs...
) where {T<:Union{Real,Quantity}}

    func_name = nameof(func)

    arg_tmp(x) = func(x, args...; kwargs...)

    npoints = length(x)
    xmin = minimum(x)
    xmax = maximum(x)

    x_units = unit(xmin)
    f_units = unit(arg_tmp(xmin))

    x_v = ustrip.(x_units, x)
    x_v = scaler.(x_v, x_scale)

    arg = wrap_function_1D_remove_units(arg_tmp, x_units, f_units)

    if isnothing(jld_base_path)
        base_path = pwd()
    else
        base_path = jld_base_path
    end

    if !(ispath(base_path))
        mkdir(base_path)
    end

    if isnothing(custom_name)
        filename = "$(func_name)_data.jld2"
    else
        filename = "$(custom_name)_data.jld2"
    end

    filepath = "$(base_path)/$(filename)"

    function load_file()
        if isfile(filepath)
            data = load(filepath)
            if data["xmin"] == xmin && data["xmax"] == xmax && data["npoints"] == npoints
                return true
            else
                @warn "The provided tabulation range or npoints does not match the data stored in the tabulation file. The tabulation will be recomputed."
                return false
            end
        else
            return false
        end
    end

    if load_file()
        data = load(filepath)
        @info "$(filename) loaded!"
        if check_SHA
            func_sha = data["sha"]
            test_sha(func, func_sha, mode = check_SHA_mode)
        end
        x_v = data["x"]

        data_matrix = data["func"]

    else
        if isnothing(custom_name)
            @info "Computing $(func_name) Interpolation"
        else
            @info "Computing $(custom_name) Interpolation"
        end

        p = Progress(Int(npoints))
        update!(p, 0)

        data_matrix = zeros(npoints)

        @threads for i = 1:npoints

            data_matrix[i] = arg(un_scaler(x_v[i], x_scale))

            next!(p)
        end

        data_dict = Dict{
            String,
            Union{
                String,
                Real,
                Quantity,
                Array,
                StepRangeLen
            },
        }()
        sha = compute_SHA(func)
        data_dict["x"] = x_v
        data_dict["func"] = convert(Array{Float64}, data_matrix)
        data_dict["sha"] = sha
        data_dict["xmin"] = xmin
        data_dict["xmax"] = xmax
        data_dict["npoints"] = npoints

        save(filepath, data_dict)
        @info "$(filename) created and exported!"
    end

    if f_scale == :log10
        data_matrix[data_matrix.<=1e-300] .= 1e-300
    end

    knots = (x_v,)
    f_matrix = scaler.(data_matrix, f_scale)

    if interpolation_type == :linear
        itp = LinearInterpolation(knots, f_matrix, extrapolation_bc = extrapolation_bc())
    elseif interpolation_type == :cubic
        throw(ArgumentError("$interpolation_type is not supported for this tabulation method"))
    else
        throw(ArgumentError("$interpolation_type is not a valid interpolation type"))
    end

    function func_interp(x::Real)
        return un_scaler(itp(scaler(x, x_scale)), f_scale)::Float64
    end

    func_interp_units = wrap_function_1D_add_units(func_interp, x_units, f_units)

    return func_interp_units::Function
end

"""
    create_tabulation_2D(
        func::Function[,
        args...];
        xmin::T,
        xmax::T,
        ymin::V,
        ymax::V,
        npoints_x::Int,
        npoints_y::Int[,
        x_scale::Symbol = :linear,
        y_scale::Symbol = :linear,
        f_scale = :linear,
        jld_base_path = nothing,
        custom_name = nothing,
        interpolation_type = :linear,
        extrapolation_bc = Throw,
        check_SHA = true,
        check_SHA_mode = :warn,
        kwargs...]
    ) where {T<:Union{Real,Quantity},V<:Union{Real,Quantity}}

Computes or loads the tabulation for a given function of the kind `f(x, y)`.

# Arguments
- `func::Function`: the function which should be tabulated. Must be of the kind `f(x, y)`
- `xmin::Union{Real,Quantity}`: the minimum `x` for the range which will be covered by the tabulation
- `xmax::Union{Real,Quantity}`: the maximum `x` for the range which will be covered by the tabulation
- `ymin::Union{Real,Quantity}`: the minimum `y` for the range which will be covered by the tabulation
- `ymax::Union{Real,Quantity}`: the maximum `y` for the range which will be covered by the tabulation
- `npoints_x::Int`: the number of points along the `x` axis for which `f(x,y)` will be computed
- `npoints_y::Int`: the number of points along the `y` axis for which `f(x,y)` will be computed
- `x_scale::Symbol`: the scale which will be used for tabulation. If `:linear` is choosen, the `x` points will be evenly spaced. If `:log10` is chosen, the points will be logarithmically spaced. Defaults to `:linear`
- `y_scale::Symbol`: the scale which will be used for tabulation. If `:linear` is choosen, the `y` points will be evenly spaced. If `:log10` is chosen, the points will be logarithmically spaced. Defaults to `:linear`
- `f_scale::Symbol`: the scale which will be used for tabulation. If `:linear` is choosen, the `f(x)` points will not be alterated. If `:log10` is chosen, log10 will be applied to the `f(x)` points before interpolating. Defaults to `:linear`
- `jld_base_path`: path to the folder where the tabulation should be saved. Defaults to the current folder
- `custom_name`: custom name for the tabulation file, to which `_data.jld2` will be appended. Defaults to the name of the function to be tabulated
- `interpolation_type`: Type of the spline to be used for interpolation. Can either be :linear (1st order spline) or :cubic (3rd order spline)
- `extrapolation_bc`: behaviour of the tabulation outside of the boundaries defined by `[xmin, xmax]`. Possible behaviours are handled by Interpolations.jl and include `Throw`` (throws and eror if a value out of bounds is acessed) or `Line` (extrapolate linearly). Defaults to `Throw`
- `check_SHA::Bool`: Whether to check the SHA of the tabulated function. If the SHA of the tabulated function does not match the SHA stored inside the tabulation file, it might be necessary to recompute the tabulation, since it's likely that the definition of the tabulated function has changed.
- `check_SHA_mode::Symbol`: Describes the behaviour of this function if the SHA of the tabulated function does not match the SHA stored inside the tabulation file. Defaults to :warn (print a warning but load the tabulation all the same). Other options include :throw (throw an error to suggest manual deletion of the old tabulation) and :none (don't do anything).
- `args..., kwargs...`: additional `args` and `kwargs` will be passed to the function which is to be tabulated

# Examples
```jldoctest
julia>  func_2d(x, y) = sin(x) * sin(y)

julia>  itp_2d_1 = create_tabulation_2D(
            func_2d,
            xmin = 0.0,
            xmax = 1.0,
            ymin = 0.0,
            ymax = 2.0,
            npoints_x = 100,
            npoints_y = 100,
        )

julia>  isapprox(itp_2d_1(1.0, 1.3), func_2d(1.0, 1.3), rtol = 1e-3)
true
```

Measurement units are supported too:

```jldoctest
julia>  using Unitful

julia>  func_2d(x, y) = x^2 + y

julia>  itp_2d_1 = create_tabulation_2D(
            func_2d,
            xmin = 0.0u"m",
            xmax = 1.0u"m",
            ymin = 0.0u"m^2",
            ymax = 2.0u"m^2",
            npoints_x = 100,
            npoints_y = 100,
        )

julia>  isapprox(itp_2d_1(1.0u"m", 1.3u"m^2"), func_2d(1.0u"m", 1.3u"m^2"), rtol = 1e-3)
true
```
"""
function create_tabulation_2D(
    func::Function,
    args...;
    jld_base_path = nothing,
    custom_name = nothing,
    xmin::T,
    xmax::T,
    ymin::V,
    ymax::V,
    npoints_x::Int,
    npoints_y::Int,
    x_scale::Symbol = :linear,
    y_scale::Symbol = :linear,
    f_scale = :linear,
    interpolation_type = :linear,
    extrapolation_bc = Throw,
    check_SHA::Bool = true,
    check_SHA_mode::Symbol = :warn,
    kwargs...
) where {T<:Union{Real,Quantity},V<:Union{Real,Quantity}}

    func_name = nameof(func)

    arg_tmp(x, y) = func(x, y, args...; kwargs...)
    x_units = unit(xmin)
    y_units = unit(ymin)
    f_units = unit(arg_tmp(xmin, ymin))

    xmin_v = ustrip(x_units, xmin)::Real
    xmax_v = ustrip(x_units, xmax)::Real
    ymin_v = ustrip(y_units, ymin)::Real
    ymax_v = ustrip(y_units, ymax)::Real

    arg = wrap_function_2D_remove_units(arg_tmp, x_units, y_units, f_units)

    if isnothing(jld_base_path)
        base_path = pwd()
    else
        base_path = jld_base_path
    end

    if !(ispath(base_path))
        mkdir(base_path)
    end

    if isnothing(custom_name)
        filename = "$(func_name)_data.jld2"
    else
        filename = "$(custom_name)_data.jld2"
    end

    filepath = "$(base_path)/$(filename)"

    function load_file()
        if isfile(filepath)
            data = load(filepath)
            if data["xmin"] == xmin && data["xmax"] == xmax && data["ymin"] == ymin && data["ymax"] == ymax && data["npoints_x"] == npoints_x && data["npoints_y"] == npoints_y
                return true
            else
                @warn "The provided tabulation range or npoints does not match the data stored in the tabulation file. The tabulation will be recomputed."

                return false
            end
        else
            return false
        end
    end

    if load_file()
        data = load(filepath)
        @info "$(filename) loaded!"
        if check_SHA
            func_sha = data["sha"]
            test_sha(func, func_sha, mode = check_SHA_mode)
        end
        x = data["x"]
        y = data["y"]

        data_matrix = data["func"]

    else
        if x_scale == :linear
            x = range(xmin_v, xmax_v, length = npoints_x)
        elseif x_scale == :log10
            x = range(log10(xmin_v), log10(xmax_v), length = npoints_x)
        else
            throw(ArgumentError("X scale $x_scale not supported"))
        end

        if y_scale == :linear
            y = range(ymin_v, ymax_v, length = npoints_y)
        elseif y_scale == :log10
            y = range(log10(ymin_v), log10(ymax_v), length = npoints_y)
        else
            throw(ArgumentError("Y scale $y_scale not supported"))
        end

        if isnothing(custom_name)
            @info "Computing $(func_name) Interpolation"
        else
            @info "Computing $(custom_name) Interpolation"
        end

        p = Progress(Int(npoints_y))
        update!(p, 0)

        data_matrix = zeros(npoints_x, npoints_y)

        @threads for j = 1:npoints_y
            for i = 1:npoints_x
                data_matrix[i, j] = arg(un_scaler(x[i], x_scale), un_scaler(y[j], y_scale))
            end
            next!(p)
        end

        data_dict = Dict{
            String,
            Union{
                String,
                Real,
                Quantity,
                Array,
                StepRangeLen
            },
        }()
        sha = compute_SHA(func)
        data_dict["x"] = x
        data_dict["y"] = y
        data_dict["func"] = convert(Array{Float64}, data_matrix)
        data_dict["sha"] = sha
        data_dict["xmin"] = xmin
        data_dict["xmax"] = xmax
        data_dict["ymin"] = ymin
        data_dict["ymax"] = ymax
        data_dict["npoints_x"] = npoints_x
        data_dict["npoints_y"] = npoints_y

        save(filepath, data_dict)
        @info "$(filename) created and exported!"
    end

    if f_scale == :log10
        data_matrix[data_matrix.<=1e-300] .= 1e-300
    end

    knots = (x, y)
    f_matrix = scaler.(data_matrix, f_scale)

    if interpolation_type == :linear
        itp = LinearInterpolation(knots, f_matrix, extrapolation_bc = extrapolation_bc())
    elseif interpolation_type == :cubic
        itp = CubicSplineInterpolation(knots, f_matrix, extrapolation_bc = extrapolation_bc())
    else
        throw(ArgumentError("$interpolation_type is not a valid interpolation type"))
    end

    function func_interp(x::Real, y::Real)
        return un_scaler(itp(scaler(x, x_scale), scaler(y, y_scale)), f_scale)::Float64
    end

    func_interp_units = wrap_function_2D_add_units(func_interp, x_units, y_units, f_units)

    return func_interp_units::Function
end

"""
    function create_tabulation_2D(
        func::Function,
        x::Vector{T},
        y::Vector{V},
        args...;
        jld_base_path = nothing,
        custom_name = nothing,
        x_scale::Symbol = :linear,
        y_scale::Symbol = :linear,
        f_scale = :linear,
        interpolation_type = :linear,
        extrapolation_bc = Throw,
        check_SHA = true,
        check_SHA_mode = :warn,
        kwargs...
    ) where {T<:Union{Real,Quantity},V<:Union{Real,Quantity}}

Provides an interface for `create_tabulation_2D` where the tabulation points are defined by the Vectors `x` and `y`. Only `LinearInterpolation` is supported.`
"""
function create_tabulation_2D(
    func::Function,
    x::Vector{T},
    y::Vector{V},
    args...;
    jld_base_path = nothing,
    custom_name = nothing,
    x_scale::Symbol = :linear,
    y_scale::Symbol = :linear,
    f_scale = :linear,
    interpolation_type = :linear,
    extrapolation_bc = Throw,
    check_SHA::Bool = true,
    check_SHA_mode::Symbol = :warn,
    kwargs...
) where {T<:Union{Real,Quantity},V<:Union{Real,Quantity}}

    func_name = nameof(func)

    arg_tmp(x, y) = func(x, y, args...; kwargs...)

    npoints_x = length(x)
    xmin = minimum(x)
    xmax = maximum(x)

    npoints_y = length(y)
    ymin = minimum(y)
    ymax = maximum(y)

    x_units = unit(xmin)
    y_units = unit(ymin)
    f_units = unit(arg_tmp(xmin, ymin))

    x_v = ustrip.(x_units, x)
    x_v = scaler.(x_v, x_scale)
    y_v = ustrip.(y_units, y)
    y_v = scaler.(y_v, y_scale)

    arg = wrap_function_2D_remove_units(arg_tmp, x_units, y_units, f_units)

    if isnothing(jld_base_path)
        base_path = pwd()
    else
        base_path = jld_base_path
    end

    if !(ispath(base_path))
        mkdir(base_path)
    end

    if isnothing(custom_name)
        filename = "$(func_name)_data.jld2"
    else
        filename = "$(custom_name)_data.jld2"
    end

    filepath = "$(base_path)/$(filename)"

    function load_file()
        if isfile(filepath)
            data = load(filepath)
            if data["xmin"] == xmin && data["xmax"] == xmax && data["ymin"] == ymin && data["ymax"] == ymax && data["npoints_x"] == npoints_x && data["npoints_y"] == npoints_y
                return true
            else
                @warn "The provided tabulation range or npoints does not match the data stored in the tabulation file. The tabulation will be recomputed."

                return false
            end
        else
            return false
        end
    end

    if load_file()
        data = load(filepath)
        @info "$(filename) loaded!"
        if check_SHA
            func_sha = data["sha"]
            test_sha(func, func_sha, mode = check_SHA_mode)
        end
        x_v = data["x"]
        y_v = data["y"]

        data_matrix = data["func"]

    else
        if isnothing(custom_name)
            @info "Computing $(func_name) Interpolation"
        else
            @info "Computing $(custom_name) Interpolation"
        end

        p = Progress(Int(npoints_y))
        update!(p, 0)

        data_matrix = zeros(npoints_x, npoints_y)

        @threads for j = 1:npoints_y
            for i = 1:npoints_x
                data_matrix[i, j] = arg(un_scaler(x_v[i], x_scale), un_scaler(y_v[j], y_scale))
            end
            next!(p)
        end

        data_dict = Dict{
            String,
            Union{
                String,
                Real,
                Quantity,
                Array,
                StepRangeLen
            },
        }()
        sha = compute_SHA(func)
        data_dict["x"] = x_v
        data_dict["y"] = y_v
        data_dict["func"] = convert(Array{Float64}, data_matrix)
        data_dict["sha"] = sha
        data_dict["xmin"] = xmin
        data_dict["xmax"] = xmax
        data_dict["ymin"] = ymin
        data_dict["ymax"] = ymax
        data_dict["npoints_x"] = npoints_x
        data_dict["npoints_y"] = npoints_y

        save(filepath, data_dict)
        @info "$(filename) created and exported!"
    end

    if f_scale == :log10
        data_matrix[data_matrix.<=1e-300] .= 1e-300
    end

    knots = (x_v, y_v)
    f_matrix = scaler.(data_matrix, f_scale)

    if interpolation_type == :linear
        itp = LinearInterpolation(knots, f_matrix, extrapolation_bc = extrapolation_bc())
    elseif interpolation_type == :cubic
        throw(ArgumentError("$interpolation_type is not supported for this tabulation method"))
    else
        throw(ArgumentError("$interpolation_type is not a valid interpolation type"))
    end

    function func_interp(x::Real, y::Real)
        return un_scaler(itp(scaler(x, x_scale), scaler(y, y_scale)), f_scale)::Float64
    end

    func_interp_units = wrap_function_2D_add_units(func_interp, x_units, y_units, f_units)

    return func_interp_units::Function
end

"""
    create_tabulation_3D(
        func::Function[,
        args...];
        xmin::T,
        xmax::T,
        ymin::V,
        ymax::V,
        zmin::W,
        zmax::W,
        npoints_x::Int,
        npoints_y::Int,
        npoints_z::Int[,
        x_scale::Symbol = :linear,
        y_scale::Symbol = :linear,
        z_scale::Symbol = :linear,
        f_scale = :linear,
        jld_base_path = nothing,
        custom_name = nothing,
        interpolation_type = :linear,
        extrapolation_bc = Throw,
        check_SHA = true,
        check_SHA_mode = :warn,
        kwargs...]
    ) where {T<:Union{Real,Quantity},V<:Union{Real,Quantity},W<:Union{Real,Quantity}}

Computes or loads the tabulation for a given function of the kind `f(x, y, z)`.

# Arguments
- `func::Function`: the function which should be tabulated. Must be of the kind `f(x, y)`
- `xmin::Union{Real,Quantity}`: the minimum `x` for the range which will be covered by the tabulation
- `xmax::Union{Real,Quantity}`: the maximum `x` for the range which will be covered by the tabulation
- `ymin::Union{Real,Quantity}`: the minimum `y` for the range which will be covered by the tabulation
- `ymax::Union{Real,Quantity}`: the maximum `y` for the range which will be covered by the tabulation
- `zmin::Union{Real,Quantity}`: the minimum `z` for the range which will be covered by the tabulation
- `zmax::Union{Real,Quantity}`: the maximum `z` for the range which will be covered by the tabulation
- `npoints_x::Int`: the number of points along the `x` axis for which `f(x,y,z)` will be computed
- `npoints_y::Int`: the number of points along the `y` axis for which `f(x,y,z)` will be computed
- `npoints_z::Int`: the number of points along the `z` axis for which `f(x,y,z)` will be computed
- `x_scale::Symbol`: the scale which will be used for tabulation. If `:linear` is choosen, the `x` points will be evenly spaced. If `:log10` is chosen, the points will be logarithmically spaced. Defaults to `:linear`
- `y_scale::Symbol`: the scale which will be used for tabulation. If `:linear` is choosen, the `y` points will be evenly spaced. If `:log10` is chosen, the points will be logarithmically spaced. Defaults to `:linear`
- `z_scale::Symbol`: the scale which will be used for tabulation. If `:linear` is choosen, the `z` points will be evenly spaced. If `:log10` is chosen, the points will be logarithmically spaced. Defaults to `:linear`
- `f_scale::Symbol`: the scale which will be used for tabulation. If `:linear` is choosen, the `f(x)` points will not be alterated. If `:log10` is chosen, log10 will be applied to the `f(x)` points before interpolating. Defaults to `:linear`
- `jld_base_path`: path to the folder where the tabulation should be saved. Defaults to the current folder
- `custom_name`: custom name for the tabulation file, to which `_data.jld2` will be appended. Defaults to the name of the function to be tabulated
- `interpolation_type`: Type of the spline to be used for interpolation. Can either be :linear (1st order spline) or :cubic (3rd order spline)
- `extrapolation_bc`: behaviour of the tabulation outside of the boundaries defined by `[xmin, xmax]`. Possible behaviours are handled by Interpolations.jl and include `Throw`` (throws and eror if a value out of bounds is acessed) or `Line` (extrapolate linearly). Defaults to `Throw`
- `check_SHA::Bool`: Whether to check the SHA of the tabulated function. If the SHA of the tabulated function does not match the SHA stored inside the tabulation file, it might be necessary to recompute the tabulation, since it's likely that the definition of the tabulated function has changed.
- `check_SHA_mode::Symbol`: Describes the behaviour of this function if the SHA of the tabulated function does not match the SHA stored inside the tabulation file. Defaults to :warn (print a warning but load the tabulation all the same). Other options include :throw (throw an error to suggest manual deletion of the old tabulation) and :none (don't do anything).
- `args..., kwargs...`: additional `args` and `kwargs` will be passed to the function which is to be tabulated

# Examples
```jldoctest
julia>  func_3d(x, y, z) = x * y + z

julia>  itp_3d_1 = create_tabulation_3D(
            func_3d,
            xmin = 0.0,
            xmax = 1.0,
            ymin = 0.0,
            ymax = 2.0,
            zmin = 0.0,
            zmax = 3.0,
            npoints_x = 100,
            npoints_y = 100,
            npoints_z = 100,
        )

julia>  isapprox(itp_3d_1(1.0, 1.3, 2.5), func_3d(1.0, 1.3, 2.5), rtol = 1e-3)
true
```

Measurement units are supported too:

```jldoctest
julia>  using Unitful

julia>  func_3d(x, y, z) = x * y + z

julia>  itp_3d_1 = create_tabulation_3D(
            func_3d,
            xmin = 0.0u"m",
            xmax = 1.0u"m",
            ymin = 0.0"s^-1",
            ymax = 2.0"s^-1",
            zmin = 0.0u"m/s",
            zmax = 3.0u"m/s",
            npoints_x = 100,
            npoints_y = 100,
            npoints_z = 100,
        )

julia>  isapprox(itp_3d_1(1.0u"m", 1.3u"s^-1", 2.5u"m/s"), func_3d(1.0u"m", 1.3u"s^-1", 2.5u"m/s"), rtol = 1e-3)
true
```
"""
function create_tabulation_3D(
    func::Function,
    args...;
    jld_base_path = nothing,
    custom_name = nothing,
    xmin::T,
    xmax::T,
    ymin::V,
    ymax::V,
    zmin::W,
    zmax::W,
    npoints_x::Int,
    npoints_y::Int,
    npoints_z::Int,
    x_scale::Symbol = :linear,
    y_scale::Symbol = :linear,
    z_scale::Symbol = :linear,
    f_scale = :linear,
    interpolation_type = :linear,
    extrapolation_bc = Throw,
    check_SHA::Bool = true,
    check_SHA_mode::Symbol = :warn,
    compress::Bool = false,
    kwargs...
) where {T<:Union{Real,Quantity},V<:Union{Real,Quantity},W<:Union{Real,Quantity}}

    func_name = nameof(func)

    arg_tmp(x, y, z) = func(x, y, z, args...; kwargs...)

    x_units = unit(xmin)
    y_units = unit(ymin)
    z_units = unit(zmin)
    f_units = unit(arg_tmp(xmin, ymin, zmin))

    xmin_v = ustrip(x_units, xmin)::Real
    xmax_v = ustrip(x_units, xmax)::Real
    ymin_v = ustrip(y_units, ymin)::Real
    ymax_v = ustrip(y_units, ymax)::Real
    zmin_v = ustrip(z_units, zmin)::Real
    zmax_v = ustrip(z_units, zmax)::Real

    arg = wrap_function_3D_remove_units(arg_tmp, x_units, y_units, z_units, f_units)

    if isnothing(jld_base_path)
        base_path = pwd()
    else
        base_path = jld_base_path
    end

    if !(ispath(base_path))
        mkdir(base_path)
    end

    if isnothing(custom_name)
        filename = "$(func_name)_data.jld2"
    else
        filename = "$(custom_name)_data.jld2"
    end

    filepath = "$(base_path)/$(filename)"

    function load_file()
        if isfile(filepath)
            data = load(filepath)
            if data["xmin"] == xmin && data["xmax"] == xmax && data["ymin"] == ymin && data["ymax"] == ymax && data["zmin"] == zmin && data["zmax"] == zmax && data["npoints_x"] == npoints_x && data["npoints_y"] == npoints_y && data["npoints_z"] == npoints_z
                return true
            else
                @warn "The provided tabulation range or npoints does not match the data stored in the tabulation file. The tabulation will be recomputed."

                return false
            end
        else
            return false
        end
    end

    if load_file()
        data = load(filepath)
        @info "$(filename) loaded!"
        if check_SHA
            func_sha = data["sha"]
            test_sha(func, func_sha, mode = check_SHA_mode)
        end
        x = data["x"]
        y = data["y"]
        z = data["z"]

        data_matrix = data["func"]

    else
        if x_scale == :linear
            x = range(xmin_v, xmax_v, length = npoints_x)
        elseif x_scale == :log10
            x = range(log10(xmin_v), log10(xmax_v), length = npoints_x)
        else
            throw(ArgumentError("X scale $x_scale not supported"))
        end

        if y_scale == :linear
            y = range(ymin_v, ymax_v, length = npoints_y)
        elseif y_scale == :log10
            y = range(log10(ymin_v), log10(ymax_v), length = npoints_y)
        else
            throw(ArgumentError("Y scale $y_scale not supported"))
        end

        if z_scale == :linear
            z = range(zmin_v, zmax_v, length = npoints_z)
        elseif z_scale == :log10
            z = range(log10(zmin_v), log10(zmax_v), length = npoints_z)
        else
            throw(ArgumentError("Z scale $z_scale not supported"))
        end

        if isnothing(custom_name)
            @info "Computing $(func_name) Interpolation"
        else
            @info "Computing $(custom_name) Interpolation"
        end

        p = Progress(Int(npoints_z * npoints_y))
        update!(p, 0)

        data_matrix = zeros(npoints_x, npoints_y, npoints_z)

        @threads for k = 1:npoints_z
            for j = 1:npoints_y
                for i = 1:npoints_x
                    data_matrix[i, j, k] = arg(un_scaler(x[i], x_scale), un_scaler(y[j], y_scale), un_scaler(z[k], z_scale))
                end
                next!(p)
            end
        end

        data_dict = Dict{
            String,
            Union{
                String,
                Real,
                Quantity,
                Array,
                StepRangeLen
            },
        }()
        sha = compute_SHA(func)
        data_dict["x"] = x
        data_dict["y"] = y
        data_dict["z"] = z
        data_dict["func"] = convert(Array{Float64}, data_matrix)
        data_dict["sha"] = sha
        data_dict["xmin"] = xmin
        data_dict["xmax"] = xmax
        data_dict["ymin"] = ymin
        data_dict["ymax"] = ymax
        data_dict["zmin"] = zmin
        data_dict["zmax"] = zmax
        data_dict["npoints_x"] = npoints_x
        data_dict["npoints_y"] = npoints_y
        data_dict["npoints_z"] = npoints_z

        save(filepath, data_dict, compress = compress)
        @info "$(filename) created and exported!"
    end

    if f_scale == :log10
        data_matrix[data_matrix.<=1e-300] .= 1e-300
    end

    knots = (x, y, z)
    f_matrix = scaler.(data_matrix, f_scale)

    if interpolation_type == :linear
        itp = LinearInterpolation(knots, f_matrix, extrapolation_bc = extrapolation_bc())
    elseif interpolation_type == :cubic
        itp = CubicSplineInterpolation(knots, f_matrix, extrapolation_bc = extrapolation_bc())
    else
        throw(ArgumentError("$interpolation_type is not a valid interpolation type"))
    end

    function func_interp(x::Real, y::Real, z::Real)
        return un_scaler(itp(scaler(x, x_scale), scaler(y, y_scale), scaler.(z, z_scale)), f_scale)::Float64
    end

    func_interp_units = wrap_function_3D_add_units(func_interp, x_units, y_units, z_units, f_units)

    return func_interp_units::Function
end

"""
    function create_tabulation_3D(
        func::Function,
        x::Vector{T},
        y::Vector{V},
        z::Vector{W},
        args...;
        jld_base_path = nothing,
        custom_name = nothing,
        x_scale::Symbol = :linear,
        y_scale::Symbol = :linear,
        z_scale::Symbol = :linear,
        f_scale = :linear,
        interpolation_type = :linear,
        extrapolation_bc = Throw,
        check_SHA = true,
        check_SHA_mode = :warn,
        kwargs...
    ) where {T<:Union{Real,Quantity},V<:Union{Real,Quantity},W<:Union{Real,Quantity}}

Provides an interface for `create_tabulation_3D` where the tabulation points are defined by the Vectors `x`, `y` and `z`. Only `LinearInterpolation` is supported.`
"""
function create_tabulation_3D(
    func::Function,
    x::Vector{T},
    y::Vector{V},
    z::Vector{W},
    args...;
    jld_base_path = nothing,
    custom_name = nothing,
    x_scale::Symbol = :linear,
    y_scale::Symbol = :linear,
    z_scale::Symbol = :linear,
    f_scale = :linear,
    interpolation_type = :linear,
    extrapolation_bc = Throw,
    check_SHA::Bool = true,
    check_SHA_mode::Symbol = :warn,
    compress::Bool = false,
    kwargs...
) where {T<:Union{Real,Quantity},V<:Union{Real,Quantity},W<:Union{Real,Quantity}}

    func_name = nameof(func)

    arg_tmp(x, y, z) = func(x, y, z, args...; kwargs...)

    npoints_x = length(x)
    xmin = minimum(x)
    xmax = maximum(x)

    npoints_y = length(y)
    ymin = minimum(y)
    ymax = maximum(y)

    npoints_z = length(z)
    zmin = minimum(z)
    zmax = maximum(z)

    x_units = unit(xmin)
    y_units = unit(ymin)
    z_units = unit(zmin)
    f_units = unit(arg_tmp(xmin, ymin, zmin))

    x_v = ustrip.(x_units, x)
    x_v = scaler.(x_v, x_scale)
    y_v = ustrip.(y_units, y)
    y_v = scaler.(y_v, y_scale)
    z_v = ustrip.(z_units, z)
    z_v = scaler.(z_v, z_scale)

    arg = wrap_function_3D_remove_units(arg_tmp, x_units, y_units, z_units, f_units)

    if isnothing(jld_base_path)
        base_path = pwd()
    else
        base_path = jld_base_path
    end

    if !(ispath(base_path))
        mkdir(base_path)
    end

    if isnothing(custom_name)
        filename = "$(func_name)_data.jld2"
    else
        filename = "$(custom_name)_data.jld2"
    end

    filepath = "$(base_path)/$(filename)"

    function load_file()
        if isfile(filepath)
            data = load(filepath)
            if data["xmin"] == xmin && data["xmax"] == xmax && data["ymin"] == ymin && data["ymax"] == ymax && data["zmin"] == zmin && data["zmax"] == zmax && data["npoints_x"] == npoints_x && data["npoints_y"] == npoints_y && data["npoints_z"] == npoints_z
                return true
            else
                @warn "The provided tabulation range or npoints does not match the data stored in the tabulation file. The tabulation will be recomputed."

                return false
            end
        else
            return false
        end
    end

    if load_file()
        data = load(filepath)
        @info "$(filename) loaded!"
        if check_SHA
            func_sha = data["sha"]
            test_sha(func, func_sha, mode = check_SHA_mode)
        end
        x_v = data["x"]
        y_v = data["y"]
        z_v = data["z"]

        data_matrix = data["func"]

    else

        if isnothing(custom_name)
            @info "Computing $(func_name) Interpolation"
        else
            @info "Computing $(custom_name) Interpolation"
        end

        p = Progress(Int(npoints_z * npoints_y))
        update!(p, 0)

        data_matrix = zeros(npoints_x, npoints_y, npoints_z)

        @threads for k = 1:npoints_z
            for j = 1:npoints_y
                for i = 1:npoints_x
                    data_matrix[i, j, k] = arg(un_scaler(x_v[i], x_scale), un_scaler(y_v[j], y_scale), un_scaler(z_v[k], z_scale))
                end
                next!(p)
            end
        end

        data_dict = Dict{
            String,
            Union{
                String,
                Real,
                Quantity,
                Array,
                StepRangeLen
            },
        }()
        sha = compute_SHA(func)
        data_dict["x"] = x_v
        data_dict["y"] = y_v
        data_dict["z"] = z_v
        data_dict["func"] = convert(Array{Float64}, data_matrix)
        data_dict["sha"] = sha
        data_dict["xmin"] = xmin
        data_dict["xmax"] = xmax
        data_dict["ymin"] = ymin
        data_dict["ymax"] = ymax
        data_dict["zmin"] = zmin
        data_dict["zmax"] = zmax
        data_dict["npoints_x"] = npoints_x
        data_dict["npoints_y"] = npoints_y
        data_dict["npoints_z"] = npoints_z

        save(filepath, data_dict, compress = compress)
        @info "$(filename) created and exported!"
    end

    if f_scale == :log10
        data_matrix[data_matrix.<=1e-300] .= 1e-300
    end

    knots = (x_v, y_v, z_v)
    f_matrix = scaler.(data_matrix, f_scale)

    if interpolation_type == :linear
        itp = LinearInterpolation(knots, f_matrix, extrapolation_bc = extrapolation_bc())
    elseif interpolation_type == :cubic
        throw(ArgumentError("$interpolation_type is not supported for this tabulation method"))
    else
        throw(ArgumentError("$interpolation_type is not a valid interpolation type"))
    end

    function func_interp(x::Real, y::Real, z::Real)
        return un_scaler(itp(scaler(x, x_scale), scaler(y, y_scale), scaler.(z, z_scale)), f_scale)::Float64
    end

    func_interp_units = wrap_function_3D_add_units(func_interp, x_units, y_units, z_units, f_units)

    return func_interp_units::Function
end

end


text:
"""
```
smc(m::AbstractDSGEModel, data::Matrix; verbose::Symbol, old_data::Matrix)
smc(m::AbstractDSGEModel, data::DataFrame)
smc(m::AbstractDSGEModel)
```

### Arguments:

- `m`: A model object, which stores parameter values, prior dists, bounds, and various
    other settings that will be referenced
- `data`: A matrix or dataframe containing the time series of the observables used in
    the calculation of the posterior/likelihood
- `old_data`: A matrix containing the time series of observables of previous data
    (with `data` being the new data) for the purposes of a time tempered estimation
    (that is, using the posterior draws from a previous estimation as the initial set
    of draws for an estimation with new data)

### Keyword Arguments:
- `verbose`: Desired frequency of function progress messages printed to standard out.
	- `:none`: No status updates will be reported.
	- `:low`: Status updates for SMC initialization and recursion will be included.
	- `:high`: Status updates for every iteration of SMC is output, which includes
    the mean and standard deviation of each parameter draw after each iteration,
    as well as calculated acceptance rate, ESS, and number of times resampled.

### Outputs

- `cloud`: The ParticleCloud object containing all of the information about the
    parameter values from the sample, their respective log-likelihoods, the ESS
    schedule, tempering schedule etc., which is saved in the saveroot.

### Overview

Sequential Monte Carlo can be used in lieu of Random Walk Metropolis Hastings to
    generate parameter samples from high-dimensional parameter spaces using
    sequentially constructed proposal densities to be used in iterative importance
    sampling.

The implementation here is based on Edward Herbst and Frank Schorfheide's 2014 paper
    'Sequential Monte Carlo Sampling for DSGE Models' and the code accompanying their
    book 'Bayesian Estimation of DSGE Models'.

SMC is broken up into three main steps:

- `Correction`: Reweight the particles from stage n-1 by defining incremental weights,
    which gradually "temper in" the likelihood function p(Y|θ)^(ϕ_n - ϕ_n-1) into the
    normalized particle weights.
- `Selection`: Resample the particles if the distribution of particles begins to
    degenerate, according to a tolerance level for the ESS.
- `Mutation`: Propagate particles {θ(i), W(n)} via N(MH) steps of a Metropolis
    Hastings algorithm.
"""
function smc(m::AbstractDSGEModel, data::Matrix{Float64};
             verbose::Symbol = :low,
             old_data::Matrix{Float64} = Matrix{Float64}(undef, size(data, 1), 0),
             old_cloud::ParticleCloud  = ParticleCloud(m, 0),
             run_test::Bool = false,
             filestring_addl::Vector{String} = Vector{String}(),
             continue_intermediate::Bool = false, intermediate_stage_start::Int = 0,
             save_intermediate::Bool = false, intermediate_stage_increment::Int = 10)

    ########################################################################################
    ### Setting Parameters
    ########################################################################################
    # Construct closure of mutation function so as to avoid issues with serialization
    # across workers with different Julia system images
    sendto(workers(), m = m)
    sendto(workers(), data = data)
    function mutation_closure(p::Vector{S}, d_μ::Vector{S}, d_Σ::Matrix{S},
                              blocks_free::Vector{Vector{Int64}},
                              blocks_all::Vector{Vector{Int64}},
                              ϕ_n::S, ϕ_n1::S; c::S = 1.0, α::S = 1.0,
                              old_data::T = Matrix{S}(undef, size(data, 1), 0),
                              use_chand_recursion::Bool = false,
                              verbose::Symbol = :low) where {S<:Float64, T<:AbstractMatrix}
        return mutation(m, data, p, d_μ, d_Σ, blocks_free, blocks_all, ϕ_n, ϕ_n1; c = c, α = α,
                        old_data = old_data, use_chand_recursion = use_chand_recursion,
                        verbose = verbose)
    end
    @everywhere function mutation_closure(p::Vector{S}, d_μ::Vector{S}, d_Σ::Matrix{S},
                                          blocks_free::Vector{Vector{Int64}},
                                          blocks_all::Vector{Vector{Int64}},
                                          ϕ_n::S, ϕ_n1::S; c::S = 1.0, α::S = 1.0,
                                          old_data::T = Matrix{S}(undef, size(data, 1), 0),
                                          use_chand_recursion::Bool = false,
                                          verbose::Symbol = :low) where {S<:Float64, T<:Matrix}
        return mutation(m, data, p, d_μ, d_Σ, blocks_free, blocks_all, ϕ_n, ϕ_n1; c = c, α = α,
                        old_data = old_data, use_chand_recursion = use_chand_recursion,
                        verbose = verbose)
    end

    # General
    parallel = get_setting(m, :use_parallel_workers)
    n_parts  = get_setting(m, :n_particles)
    n_blocks = get_setting(m, :n_smc_blocks)
    n_steps  = get_setting(m, :n_mh_steps_smc)
    n_params = n_parameters(m)

    use_chand_recursion = get_setting(m, :use_chand_recursion)
    if any(isnan.(data)) & use_chand_recursion
        error("Cannot use Chandrasekhar recursions with missing data")
    end

    # Time Tempering
    tempered_update = !isempty(old_data)

    # Check that if there's a tempered update, old and current vintages are different
    if tempered_update
        old_vintage = get_setting(m, :previous_data_vintage)
        @assert old_vintage != data_vintage(m)
    end

    # Step 1 (ϕ schedule) settings
    i = 1                         # Index tracking the stage of the algorithm
    j = 2                         # Index tracking the fixed_schedule entry ϕ_prop
    ϕ_n    = 0.                   # Instantiate ϕ_n and ϕ_prop variables for
    ϕ_prop = 0.                   # reference in respective while loop conditions
    resampled_last_period = false # Ensures proper resetting of ESS_bar after resample
    λ      = get_setting(m, :λ)
    n_Φ    = get_setting(m, :n_Φ)
    tempering_target   = get_setting(m, :adaptive_tempering_target_smc)
    use_fixed_schedule = tempering_target == 0.0

    # Step 2 (Correction) settings
    resampling_method = get_setting(m, :resampler_smc)
    threshold_ratio   = get_setting(m, :resampling_threshold)
    threshold         = threshold_ratio * n_parts

    # Step 3 (Mutation) settings
    c      = get_setting(m, :step_size_smc)
    α      = get_setting(m, :mixture_proportion)
    target = accept = get_setting(m, :target_accept)

    para_symbols    = [θ.key for θ in m.parameters]
    fixed_para_inds = findall([θ.fixed for θ in m.parameters])
    free_para_inds  = findall([!θ.fixed for θ in m.parameters])
    n_free_para     = length(free_para_inds)

    # Step 4 (Initialization of Particle Array Cloud)
    cloud = Cloud(m, n_parts)

    #################################################################################
    ### Initialize Algorithm: Draws from prior
    #################################################################################

    if VERBOSITY[verbose] >= VERBOSITY[:low]
        println("\n\n SMC " * (m.testing ? "testing " : "") * "starts ....\n\n")
    end

    if tempered_update
        cloud = if isempty(old_cloud)
            loadpath = rawpath(m, "estimate", "smc_cloud.jld2", filestring_addl)
            loadpath = replace(loadpath, r"vint=[0-9]{6}", "vint=" * old_vintage)
            Cloud(load(loadpath, "cloud"))
        else
            Cloud(old_cloud)
        end
        initialize_cloud_settings!(m, cloud; tempered_update = tempered_update)
        initialize_likelihoods!(m, data, cloud, parallel = parallel, verbose = verbose)
    elseif continue_intermediate
        loadpath = rawpath(m, "estimate", "smc_cloud_stage=$(intermediate_stage_start).jld2",
                           filestring_addl)
        cloud    = Cloud(load(loadpath, "cloud"))
    else
        # Instantiating Cloud object, update draws, loglh, & logpost
        initial_draw!(m, data, cloud, parallel = parallel,
                      use_chand_recursion = use_chand_recursion, verbose = verbose)
        initialize_cloud_settings!(m, cloud; tempered_update = tempered_update)
    end

    # Fixed schedule for construction of ϕ_prop
    if use_fixed_schedule
        cloud.tempering_schedule = ((collect(1:n_Φ) .- 1) / (n_Φ-1)) .^ λ
    else
        proposed_fixed_schedule  = ((collect(1:n_Φ) .- 1) / (n_Φ-1)) .^ λ
    end

    # Instantiate incremental and normalized weight matrices for logMDD calculation
    if continue_intermediate
        z_matrix = load(loadpath, "z")
        w_matrix = load(loadpath, "w")
        W_matrix = load(loadpath, "W")
        j        = load(loadpath, "j")
        i        = cloud.stage_index
        c        = cloud.c
        ϕ_prop   = proposed_fixed_schedule[j]
    else
        z_matrix = ones(1)
        w_matrix = zeros(n_parts, 1)
        W_matrix = tempered_update ? get_weights(cloud) : fill(1/n_parts, (n_parts, 1))
    end

    if VERBOSITY[verbose] >= VERBOSITY[:low]
        init_stage_print(cloud, para_symbols; verbose = verbose,
                         use_fixed_schedule = use_fixed_schedule)
    end

    #################################################################################
    ### Recursion
    #################################################################################
    (VERBOSITY[verbose] >= VERBOSITY[:low]) && println("\n\n SMC recursion starts \n\n")

    while ϕ_n < 1.

        start_time = time_ns()
        cloud.stage_index = i += 1

        #############################################################################
        ### Step 0: Setting ϕ_n (either adaptively or by the fixed schedule)
        #############################################################################
        ϕ_n1 = cloud.tempering_schedule[i-1]

        if use_fixed_schedule
            ϕ_n = cloud.tempering_schedule[i]
        else
            ϕ_n, resampled_last_period, j, ϕ_prop = solve_adaptive_ϕ(cloud,
                                                       proposed_fixed_schedule,
                                                       i, j, ϕ_prop, ϕ_n1,
                                                       tempering_target,
                                                       resampled_last_period)
        end

        #############################################################################
        ### Step 1: Correction
        #############################################################################
        # Calculate incremental weights (if no old data, get_old_loglh(cloud) = 0)
        incremental_weights = exp.((ϕ_n1 - ϕ_n) * get_old_loglh(cloud) +
                                   (ϕ_n - ϕ_n1) * get_loglh(cloud))

        # Update weights
        update_weights!(cloud, incremental_weights)
        mult_weights = get_weights(cloud)

        # Normalize weights
        normalize_weights!(cloud)
        normalized_weights = get_weights(cloud)

        push!(z_matrix, sum(mult_weights))
        w_matrix = hcat(w_matrix, incremental_weights)
        W_matrix = hcat(W_matrix, normalized_weights)

        ##############################################################################
        ### Step 2: Selection
        ##############################################################################

        # Calculate the degeneracy/effective sample size metric
        push!(cloud.ESS, 1 / sum(normalized_weights .^ 2))

        # If this assertion does not hold then there are likely too few particles
        @assert !isnan(cloud.ESS[i]) "no particles have non-zero weight"

        # Resample if degeneracy/ESS metric falls below the accepted threshold
        if (cloud.ESS[i] < threshold)
            # Resample according to particle weights, uniformly reset weights to 1/n_parts
            new_inds = resample(normalized_weights; method = resampling_method)
            cloud.particles = [deepcopy(cloud.particles[i,j]) for i in new_inds,
                               j=1:size(cloud.particles, 2)]
            reset_weights!(cloud)
            cloud.resamples += 1
            resampled_last_period = true
            W_matrix[:, i] = fill(1/n_parts, (n_parts, 1))
        end

        ##############################################################################
        ### Step 3: Mutation
        ##############################################################################

        # Calculate adaptive c-step for use as scaling coefficient in mutation MH step
        c = c * (0.95 + 0.10 * exp(16.0 * (cloud.accept - target)) /
                 (1.0 + exp(16.0 * (cloud.accept - target))))
        cloud.c = c

        θ_bar = weighted_mean(cloud)
        R     = weighted_cov(cloud)

        # Ensures marix is positive semi-definite symmetric
        # (not off due to numerical error) and values haven't changed
        R_fr = (R[free_para_inds, free_para_inds] + R[free_para_inds, free_para_inds]') / 2.

        # MvNormal centered at ̄θ with var-cov ̄Σ, subsetting out the fixed parameters
        θ_bar_fr = θ_bar[free_para_inds]

        # Generate random parameter blocks
        blocks_free = generate_free_blocks(n_free_para, n_blocks)
        blocks_all  = generate_all_blocks(blocks_free, free_para_inds)

        new_particles = if parallel
            @distributed (hcat) for k in 1:n_parts
                mutation_closure(cloud.particles[k, :], θ_bar_fr, R_fr, blocks_free,
                                 blocks_all, ϕ_n, ϕ_n1; c = c, α = α, old_data = old_data,
                                 use_chand_recursion = use_chand_recursion,
                                 verbose = verbose)
            end
        else
            hcat([mutation_closure(cloud.particles[k, :], θ_bar_fr, R_fr,
                                   blocks_free, blocks_all, ϕ_n, ϕ_n1; c = c,
                                   α = α, old_data = old_data,
                                   use_chand_recursion = use_chand_recursion,
                                   verbose = verbose) for k=1:n_parts]...)
        end
        update_cloud!(cloud, new_particles)
        update_acceptance_rate!(cloud)

        ##############################################################################
        ### Timekeeping and Output Generation
        ##############################################################################
        cloud.total_sampling_time += Float64((time_ns() - start_time) * 1e-9)

        if VERBOSITY[verbose] >= VERBOSITY[:low]
            end_stage_print(cloud, para_symbols; verbose = verbose,
                            use_fixed_schedule = use_fixed_schedule)
        end

        if run_test && (i == 3)
            break
        end
        if mod(cloud.stage_index, intermediate_stage_increment) == 0 && save_intermediate
            jldopen(rawpath(m, "estimate", "smc_cloud_stage=$(cloud.stage_index).jld2"),
                    true, true, true, IOStream) do file
                write(file, "cloud", ParticleCloud(cloud, para_symbols))
                write(file, "w", w_matrix)
                write(file, "W", W_matrix)
                write(file, "z", z_matrix)
                write(file, "j", j)
            end
        end
    end

    ##################################################################################
    ### Saving data
    ##################################################################################

    if !m.testing
        simfile = h5open(rawpath(m, "estimate", "smcsave.h5", filestring_addl), "w")
        particle_store = d_create(simfile, "smcparams", datatype(Float64),
                                  dataspace(n_parts, n_params))
        for i in 1:n_parts; particle_store[i,:] = cloud.particles[i, 1:n_params] end
        close(simfile)

        jldopen(rawpath(m, "estimate", "smc_cloud.jld2", filestring_addl),
                true, true, true, IOStream) do file
            write(file, "cloud", ParticleCloud(cloud, para_symbols))
            write(file, "w", w_matrix)
            write(file, "W", W_matrix)
            write(file, "z", z_matrix)
        end
    end
end

function smc(m::AbstractDSGEModel, data::DataFrame; verbose::Symbol = :low,
             save_intermediate::Bool = false, intermediate_stage_increment::Int = 10,
             filestring_addl::Vector{String} = Vector{String}(undef, 0))
    data_mat = df_to_matrix(m, data)
    return smc(m, data_mat, verbose = verbose, save_intermediate = save_intermediate,
               filestring_addl = filestring_addl)
end

function smc(m::AbstractDSGEModel; verbose::Symbol = :low,
             save_intermediate::Bool = false, intermediate_stage_increment::Int = 10,
             filestring_addl::Vector{String} = Vector{String}(undef, 0))
    data = load_data(m)
    data_mat = df_to_matrix(m, data)
    return smc(m, data_mat, verbose=verbose, save_intermediate = save_intermediate,
               filestring_addl = filestring_addl)
end


text:
using AlphaShapes
using Test
using LinearAlgebra
include("../src/utils.jl")

@testset "Test Alpha Shape 2D" begin
   N = 1000
   r1 = 0.75
   r2 = 2.0
   X = randn(N,2)
   ind = [((norm(X[i,:])>r1) + (norm(X[i,:])<r2)) == 2 for i in 1:size(X,1)]
   X = X[ind.==true,:]
   AlphaShape(X;α=nothing,search=(0.0, 10.0),MaxSteps=100)
   t = (@timed AlphaShape(X;α=nothing,search=(0.0, 10.0),MaxSteps=100))[2]
   @info "Runtime for 1000 2D points and 100 optimisation steps", t
end


text:
module CycleGAN

using Flux.Optimise: update!
using Flux.Losses: mse, mae
using Statistics
using Zygote

include("cyclegan.jl")

export Discriminator
export Generator

export Cyclegan

export train_gan

end # end module


text:
using UpdatableCholeskyFactorizations
using LinearAlgebra
using BenchmarkTools

elty = Float64
n = 128
m = 16
A_full = randn(elty, 2n, 2n)
A_full = A_full'A_full
A = A_full[1:n, 1:n]
a = A_full[1:n+1, n+1]
B = A_full[1:2n, n+1:2n]
Bm = A_full[1:n+m, n+1:n+m] # adding m columns
C = updatable_cholesky(A)

set = @benchmarkset "updatable cholesky" begin # for n in 1:16 IDEA: could have loop over different sizes
    @case "$n by $n cholesky" cholesky($A)
    @case "$(2n) by $(2n) cholesky" cholesky($A_full)
    @case "$n by $n updatable_cholesky" updatable_cholesky($A)
    @case "adding vector to $n by $n UpdatableCholesky" add_column!(D, $a) setup=(D=copy($C))
    @case "adding $m vectors to $n by $n UpdatableCholesky" add_column!(D, $Bm) setup=(D=copy($C))
    @case "adding $n x $n matrix to $n by $n UpdatableCholesky" add_column!(D, $B) setup=(D=copy($C))
end

result = run(set)

r = minimum(result["updatable cholesky"])
# times = Vector{Float64}(undef, length(r))
# titles = Vector{String}(undef, length(r))
# i = 1
for (title, trial) in r.data
    println()
    println(title)
    display(trial)
    # times[i] = trial.time
    # titles[i] = title
    # i += 1
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ContentSummaryResult <: SwaggerModel
    ContentSummary::Any # spec type: Union{ Nothing, ContentSummary } # spec name: ContentSummary

    function ContentSummaryResult(;ContentSummary=nothing)
        o = new()
        validate_property(ContentSummaryResult, Symbol("ContentSummary"), ContentSummary)
        setfield!(o, Symbol("ContentSummary"), ContentSummary)
        o
    end
end # type ContentSummaryResult

const _property_map_ContentSummaryResult = Dict{Symbol,Symbol}(Symbol("ContentSummary")=>Symbol("ContentSummary"))
const _property_types_ContentSummaryResult = Dict{Symbol,String}(Symbol("ContentSummary")=>"ContentSummary")
Base.propertynames(::Type{ ContentSummaryResult }) = collect(keys(_property_map_ContentSummaryResult))
Swagger.property_type(::Type{ ContentSummaryResult }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ContentSummaryResult[name]))}
Swagger.field_name(::Type{ ContentSummaryResult }, property_name::Symbol) =  _property_map_ContentSummaryResult[property_name]

function check_required(o::ContentSummaryResult)
    true
end

function validate_property(::Type{ ContentSummaryResult }, name::Symbol, val)
end


text:
# unfished
const _solarized_palette = (:yellow, :orange, :red, :magenta, :violet, :blue,
                               :cyan, :green)

const _solarized_colors = Dict{Symbol, Tuple}(
                          :base03 => (0, 43, 54),
                          :base02 => (7, 54, 66),
                          :base01 => (88, 110, 117),
                          :base00 => (101, 123, 131),
                          :base0  => (131, 148, 150),
                          :base1  => (147, 161, 161),
                          :base2  => (238, 232, 213),
                          :base3  => (253, 246, 227),
                          :yellow => (181, 137, 0),
                          :orange => (203, 75, 22),
                          :red    => (220, 50, 47),
                          :magenta=> (211, 54, 130),
                          :violet => (108, 113, 196),
                          :blue   => (38, 139, 210),
                          :cyan   => (42, 161, 152),
                          :green  => (133, 153, 0))

_get_solarized_color(c) = _255_to_1(c, _solarized_colors)

add_theme(:solarized_base,
    palette  = [map(c->_get_solarized_color(c),_solarized_palette)...],
    bglegend = _invisible,
    fglegend = _invisible)

add_theme(:solarized,
    base = :solarized_base,
    bginside = _get_solarized_color(:base03),
    bgoutside= _get_solarized_color(:base02),
    fg       = _get_solarized_color(:base0),
    fgtext   = _get_solarized_color(:base0),
    fglegend = _get_solarized_color(:base0),
    fgguide  = _get_solarized_color(:base1))

add_theme(:solarized_dark, base = :solarized)
add_theme(:solarized_alldark,
          base = :solarized_dark,
          bgoutside = _get_solarized_color(:base03))

add_theme(:solarized_bright,
          base = :solarized_base,
          bginside = _get_solarized_color(:base3),
          bgoutside= _get_solarized_color(:base2),
          fg       = _get_solarized_color(:base00),
          fgtext   = _get_solarized_color(:base00),
          fgguide  = _get_solarized_color(:base01))

add_theme(:solarized_allbright,
          base = :solarized_bright,
          bgoutside = _get_solarized_color(:base2))


text:
import Base: haskey, getindex, setindex!, start, done, next
# A number of sources are put together
type SourceList <: AbstractSource
  sources::Dict{Any, AbstractSource}
  # helps to build a list of sources
  SourceList(args...) = new(Dict{Any, AbstractSource}(args...))
end

function SourceList(main::AbstractSource, args...)
  result = SourceList(args...)
  result[""] = main
  result
end

getindex(d::SourceList, key) = getindex(d.sources, key)
setindex!(d::SourceList, value::AbstractSource, key) = setindex!(d.sources, value, key)
haskey(d::SourceList, key) = haskey(d.sources, key)
start(d::SourceList) = start(d.sources)
done(d::SourceList, state) = done(d.sources, state)
next(d::SourceList, state) = next(d.sources, state)

type SubPackage <: AbstractPackage
  parent::AbstractPackage
  subdir::String
end

sourcedir(p::SubPackage, workdir::String) = joinpath(sourcedir(p.parent, workdir), p.subdir)

function download(package::AbstractPackage, sources::SourceList, workdir::String; kwargs...)
  if haskey(sources, "")
    download(package, sources[""], workdir; kwargs...)
  end
  for (dir, source) in sources
    path = string(dir)
    if length(path) == 0 continue end
    if path[1] == '/'
      path = path[2:end]
    end
    println("downaload $source")
    download(SubPackage(package, path), source, workdir; kwargs...)
  end
end

max_stars_repo_name: UCL/DebbyPacker

max_issues_repo_name: UCL/DebbyPacker

max_forks_repo_name: UCL/DebbyPacker