text:
const duckdb_database = Ptr{Cvoid}
const duckdb_config = Ptr{Cvoid}
const duckdb_connection = Ptr{Cvoid}
const duckdb_prepared_statement = Ptr{Cvoid}
const duckdb_logical_type = Ptr{Cvoid}
const duckdb_data_chunk = Ptr{Cvoid}
const duckdb_vector = Ptr{Cvoid}
const duckdb_appender = Ptr{Cvoid}
const duckdb_logical_type = Ptr{Cvoid}
const duckdb_value = Ptr{Cvoid}
const duckdb_table_function = Ptr{Cvoid}
const duckdb_bind_info = Ptr{Cvoid}
const duckdb_init_info = Ptr{Cvoid}
const duckdb_function_info = Ptr{Cvoid}
const duckdb_replacement_scan_info = Ptr{Cvoid}
const DuckDBSuccess = 0;
const DuckDBError = 1;
const duckdb_state = Int32;

@enum DUCKDB_TYPE_::UInt32 begin
    DUCKDB_TYPE_INVALID = 0
    DUCKDB_TYPE_BOOLEAN
    DUCKDB_TYPE_TINYINT
    DUCKDB_TYPE_SMALLINT
    DUCKDB_TYPE_INTEGER
    DUCKDB_TYPE_BIGINT
    DUCKDB_TYPE_UTINYINT
    DUCKDB_TYPE_USMALLINT
    DUCKDB_TYPE_UINTEGER
    DUCKDB_TYPE_UBIGINT
    DUCKDB_TYPE_FLOAT
    DUCKDB_TYPE_DOUBLE
    DUCKDB_TYPE_TIMESTAMP
    DUCKDB_TYPE_DATE
    DUCKDB_TYPE_TIME
    DUCKDB_TYPE_INTERVAL
    DUCKDB_TYPE_HUGEINT
    DUCKDB_TYPE_VARCHAR
    DUCKDB_TYPE_BLOB
    DUCKDB_TYPE_DECIMAL
    DUCKDB_TYPE_TIMESTAMP_S
    DUCKDB_TYPE_TIMESTAMP_MS
    DUCKDB_TYPE_TIMESTAMP_NS
    DUCKDB_TYPE_ENUM
    DUCKDB_TYPE_LIST
    DUCKDB_TYPE_STRUCT
    DUCKDB_TYPE_MAP
    DUCKDB_TYPE_UUID
    DUCKDB_TYPE_JSON
end

const DUCKDB_TYPE = DUCKDB_TYPE_

"""
Days are stored as days since 1970-01-01\n
Use the duckdb_from_date/duckdb_to_date function to extract individual information

"""
struct duckdb_date
    days::Int32
end

struct duckdb_date_struct
    year::Int32
    month::Int8
    day::Int8
end

"""
Time is stored as microseconds since 00:00:00\n
Use the duckdb_from_time/duckdb_to_time function to extract individual information

"""
struct duckdb_time
    micros::Int64
end

struct duckdb_time_struct
    hour::Int8
    min::Int8
    sec::Int8
    micros::Int32
end

"""
Timestamps are stored as microseconds since 1970-01-01\n
Use the duckdb_from_timestamp/duckdb_to_timestamp function to extract individual information

"""
struct duckdb_timestamp
    micros::Int64
end

struct duckdb_timestamp_struct
    date::Ref{duckdb_date_struct}
    time::Ref{duckdb_time_struct}
end

struct duckdb_interval
    months::Int32
    days::Int32
    micros::Int64
end

"""
Hugeints are composed in a (lower, upper) component\n
The value of the hugeint is upper * 2^64 + lower\n
For easy usage, the functions duckdb_hugeint_to_double/duckdb_double_to_hugeint are recommended

"""
struct duckdb_hugeint
    lower::UInt64
    upper::Int64
end

struct duckdb_string_t
    length::UInt32
    data::NTuple{12, UInt8}
end

struct duckdb_list_entry_t
    offset::UInt64
    length::UInt64
end

STRING_INLINE_LENGTH = 12

struct duckdb_column
    __deprecated_data::Ptr{Cvoid}
    __deprecated_nullmask::Ptr{UInt8}
    __deprecated_type::Ptr{DUCKDB_TYPE}
    __deprecated_name::Ptr{UInt8}
    internal_data::Ptr{Cvoid}
end

struct duckdb_result
    __deprecated_column_count::Ptr{UInt64}
    __deprecated_row_count::Ptr{UInt64}
    __deprecated_rows_changed::Ptr{UInt64}
    __deprecated_columns::Ptr{duckdb_column}
    __deprecated_error_message::Ptr{UInt8}
    internal_data::Ptr{Cvoid}
end

INTERNAL_TYPE_MAP = Dict(
    DUCKDB_TYPE_BOOLEAN => Bool,
    DUCKDB_TYPE_TINYINT => Int8,
    DUCKDB_TYPE_SMALLINT => Int16,
    DUCKDB_TYPE_INTEGER => Int32,
    DUCKDB_TYPE_BIGINT => Int64,
    DUCKDB_TYPE_UTINYINT => UInt8,
    DUCKDB_TYPE_USMALLINT => UInt16,
    DUCKDB_TYPE_UINTEGER => UInt32,
    DUCKDB_TYPE_UBIGINT => UInt64,
    DUCKDB_TYPE_FLOAT => Float32,
    DUCKDB_TYPE_DOUBLE => Float64,
    DUCKDB_TYPE_TIMESTAMP => Int64,
    DUCKDB_TYPE_TIMESTAMP_S => Int64,
    DUCKDB_TYPE_TIMESTAMP_MS => Int64,
    DUCKDB_TYPE_TIMESTAMP_NS => Int64,
    DUCKDB_TYPE_DATE => Int32,
    DUCKDB_TYPE_TIME => Int64,
    DUCKDB_TYPE_INTERVAL => duckdb_interval,
    DUCKDB_TYPE_HUGEINT => duckdb_hugeint,
    DUCKDB_TYPE_UUID => duckdb_hugeint,
    DUCKDB_TYPE_VARCHAR => duckdb_string_t,
    DUCKDB_TYPE_JSON => duckdb_string_t,
    DUCKDB_TYPE_BLOB => duckdb_string_t,
    DUCKDB_TYPE_UUID => duckdb_hugeint,
    DUCKDB_TYPE_LIST => duckdb_list_entry_t,
    DUCKDB_TYPE_STRUCT => Cvoid,
    DUCKDB_TYPE_MAP => Cvoid
)

JULIA_TYPE_MAP = Dict(
    DUCKDB_TYPE_INVALID => Missing,
    DUCKDB_TYPE_BOOLEAN => Bool,
    DUCKDB_TYPE_TINYINT => Int8,
    DUCKDB_TYPE_SMALLINT => Int16,
    DUCKDB_TYPE_INTEGER => Int32,
    DUCKDB_TYPE_BIGINT => Int64,
    DUCKDB_TYPE_HUGEINT => Int128,
    DUCKDB_TYPE_UTINYINT => UInt8,
    DUCKDB_TYPE_USMALLINT => UInt16,
    DUCKDB_TYPE_UINTEGER => UInt32,
    DUCKDB_TYPE_UBIGINT => UInt64,
    DUCKDB_TYPE_FLOAT => Float32,
    DUCKDB_TYPE_DOUBLE => Float64,
    DUCKDB_TYPE_DATE => Date,
    DUCKDB_TYPE_TIME => Time,
    DUCKDB_TYPE_TIMESTAMP => DateTime,
    DUCKDB_TYPE_TIMESTAMP_S => DateTime,
    DUCKDB_TYPE_TIMESTAMP_MS => DateTime,
    DUCKDB_TYPE_TIMESTAMP_NS => DateTime,
    DUCKDB_TYPE_INTERVAL => Dates.CompoundPeriod,
    DUCKDB_TYPE_UUID => UUID,
    DUCKDB_TYPE_VARCHAR => AbstractString,
    DUCKDB_TYPE_JSON => AbstractString,
    DUCKDB_TYPE_ENUM => AbstractString,
    DUCKDB_TYPE_BLOB => Base.CodeUnits{UInt8, String},
    DUCKDB_TYPE_MAP => Dict
)

# convert a DuckDB type into Julia equivalent
function duckdb_type_to_internal_type(x::DUCKDB_TYPE)
    if !haskey(INTERNAL_TYPE_MAP, x)
        throw(NotImplementedException(string("Unsupported type for duckdb_type_to_internal_type: ", x)))
    end
    return INTERNAL_TYPE_MAP[x]
end

function duckdb_type_to_julia_type(x)
    type_id = get_type_id(x)
    if type_id == DUCKDB_TYPE_DECIMAL
        internal_type_id = get_internal_type_id(x)
        scale = get_decimal_scale(x)
        if internal_type_id == DUCKDB_TYPE_SMALLINT
            return FixedDecimal{Int16, scale}
        elseif internal_type_id == DUCKDB_TYPE_INTEGER
            return FixedDecimal{Int32, scale}
        elseif internal_type_id == DUCKDB_TYPE_BIGINT
            return FixedDecimal{Int64, scale}
        elseif internal_type_id == DUCKDB_TYPE_HUGEINT
            return FixedDecimal{Int128, scale}
        else
            throw(NotImplementedException("Unimplemented internal type for decimal"))
        end
    elseif type_id == DUCKDB_TYPE_LIST
        return Vector{Union{Missing, duckdb_type_to_julia_type(get_list_child_type(x))}}
    elseif type_id == DUCKDB_TYPE_STRUCT
        child_count = get_struct_child_count(x)
        names::Vector{Symbol} = Vector()
        for i in 1:child_count
            child_name::Symbol = Symbol(get_struct_child_name(x, i))
            push!(names, child_name)
        end
        names_tuple = Tuple(x for x in names)
        return Union{Missing, NamedTuple{names_tuple}}
    end
    if !haskey(JULIA_TYPE_MAP, type_id)
        throw(NotImplementedException(string("Unsupported type for duckdb_type_to_julia_type: ", type_id)))
    end
    return JULIA_TYPE_MAP[type_id]
end

const ROUNDING_EPOCH_TO_UNIX_EPOCH_DAYS = 719528
const ROUNDING_EPOCH_TO_UNIX_EPOCH_MS = 62167219200000

sym(ptr) = ccall(:jl_symbol, Ref{Symbol}, (Ptr{UInt8},), ptr)

max_stars_repo_name: lokax/duckdb

max_issues_repo_name: lokax/duckdb

max_forks_repo_name: lokax/duckdb


text:
__precompile__()
module TypographicUnits
export pt, bp, pc, dd, cc, nd, nc, sp, em, ex, px, en

using Unitful
using Unitful: inch, Length
import Base: +, *, -, /, show

## Absolute typographic units ################################################

# (Based on LaTeX values, cf. https://en.wikibooks.org/wiki/LaTeX/Lengths)

@unit pt "pt" Point       100inch//7227  false
@unit bp "bp" BigPoint       inch//72    false
@unit pc "pc" Pica           12pt        false
@unit dd "dd" Didôt        1238pt//1157  false
@unit cc "cc" Cîcero           dd//12    false
@unit nd "nd" NewDidôt      685pt//642   false
@unit nc "nc" NewCicero    1370pt//107   false
@unit sp "sp" ScaledPoint      pt//65536 false

## Relative typographic units ################################################

# The magnitudes of the units em, ex and px are unknown, and so lenghts
# measured in these units are incommensurable with other lengths, and with
# each other. Therefore they are each given a separate length-like dimension.

@dimension 𝐄𝐌 "𝐄𝐌" EmLength
@dimension 𝐄𝐗 "𝐄𝐗" ExLength
@dimension 𝐏𝐗 "𝐏𝐗" PxLength

@refunit   em "em" Em 𝐄𝐌    false
@refunit   ex "ex" Ex 𝐄𝐗    false
@refunit   px "px" Px 𝐏𝐗    false
@unit      en "en" En em//2 false

## A mixture of incommensurable length-like values ###########################

struct MixedLength{T <: Real,
    L<:Length{T}, M<:EmLength{T}, X<:ExLength{T}, P<:PxLength{T}} <: Number
    len::L
    ems::M
    exs::X
    pxs::P
end

MixedLength(len::L, ems::M, exs::X, pxs::P) where {T <: Real,
    L <: Length{T}, M <: EmLength{T}, X <: ExLength{T}, P <: PxLength{T}} =
    MixedLength{T, L, M, X, P}(len, ems, exs, pxs)

MixedLength(len::L, ems::M, exs::X, pxs::P) where
    {L <: Length, M <: EmLength, X <: ExLength, P <: PxLength} =
    MixedLength(_promote(len, ems, exs, pxs)...)

const M = MixedLength
getfields(x::M) = x.len, x.ems, x.exs, x.pxs

M(x::Length)   = M(  x, 0em, 0ex, 0px)
M(x::EmLength) = M(0pt,   x, 0ex, 0px)
M(x::ExLength) = M(0pt, 0em,   x, 0px)
M(x::PxLength) = M(0pt, 0em, 0ex,   x)

## Some limited mixed arithmetic #############################################

+(x::M, y::M)    = M(getfields(x) .+ getfields(y)...)
-(x::M, y::M)    = M(getfields(x) .- getfields(y)...)
-(x::M)          = M(.- getfields(x)...)
*(x::Real, y::M) = M(x .* getfields(y)...)

*(x::M, y::Real) = y * x
/(x::M, y::Real) = x * inv(y)

# Unitful.jl already uses the promotion mechanism on the generated
# unit-related types, so promoting two incommensurable length-like quantities
# to a MixedLength may not be altogether straightforward. Instead, we'll just
# hard-code the operators we need:

const LengthNames = (:Length, :EmLength, :ExLength, :PxLength, :M)
for S in LengthNames, T in LengthNames
    S == T && continue
    @eval +(x::$S, y::$T) = M(x) + M(y)
    @eval -(x::$S, y::$T) = M(x) - M(y)
end

## Mathematical functions ####################################################

# For now, just float. Any other functions that make sense here would be
# easy enough to add.

import Base: float

for func in [:float]
    @eval $func(x::MixedLength) = MixedLength($func.(getfields(x))...)
end

## Pretty-printing ###########################################################

function show(io::IO, x::M)

    compact = get(io, :compact, false)
    plus = compact ? "+" : " + "
    minus = compact ? "-" : " - "

    first = true
    for val in getfields(x)
        if !iszero(val)
            if !first
                print(io, sign(val) < 0 ? minus : plus)
                val = abs(val)
            end
            show(io, val)
            first = false
        end
    end
    first && show(io, x.len)

end

##############################################################################

# Currently (v. 0.7.0), the Unitful.jl promotion implementation doesn't work
# for more than three values [1], so here's a local alternative, sidestepping
# much of the promotion machinery.
# [1]: https://github.com/ajkeller34/Unitful.jl/issues/119
_promote(x...) = (T=promote_type(map(typeof, x)...); map(y->convert(T, y), x))

##############################################################################

const localunits = Unitful.basefactors
function __init__()

    merge!(Unitful.basefactors, localunits)
    Unitful.register(TypographicUnits)

    # Seems to be required (Unitful 0.7.0), even though this is called by
    # @refunit:
    Unitful.preferunits(em, ex, px)

    # We probably don't want to convert to meters if we mix in mm or cm.
    # (Requires that this is imported before, say, Unitful.DefaultSymbols.)
    Unitful.preferunits(pt)

end

end

ext: jl


text:
using Plots

sqrtN = 100
N = sqrtN^2
population = fill(1, (sqrtN, sqrtN))
# 1 == susceptible, 2 == infectious, 3 == recovered

# One random infectious person
population[rand(1:N)] = 2

# Probability of infected individual infecting susceptible neighbour
p = 0.4
# Probability of infected individual recovering
q = 0.01

anim = Animation()
function populationplot(population::Matrix{Int64})
    return heatmap(
        hcat(population, fill(missing, sqrtN), [1,2,3, fill(missing, sqrtN - 3)...]),
        legend = false,
        color = [:gold, :red, :blue],
        size = (4*sqrtN,4*sqrtN),
        showaxis = false,
        ticks = false
    )
end
frame(anim, populationplot(population))

maxsteps = 500
S = zeros(maxsteps+1)
S[1] = N - 1
I = zeros(maxsteps+1)
I[1] = 1

for n ∈ 2:(maxsteps+1)
    reds = CartesianIndex{2}[]
    blues = CartesianIndex{2}[]
    for j ∈ 1:sqrtN, i ∈ 1:sqrtN
        if population[i,j] == 2
            i > 1     && population[i-1,j] == 1 && rand() < p && push!(reds,CartesianIndex(i-1,j))
            j > 1     && population[i,j-1] == 1 && rand() < p && push!(reds,CartesianIndex(i,j-1))
            i < sqrtN && population[i+1,j] == 1 && rand() < p && push!(reds,CartesianIndex(i+1,j))
            j < sqrtN && population[i,j+1] == 1 && rand() < p && push!(reds,CartesianIndex(i,j+1))
            rand() < q && push!(blues,CartesianIndex(i,j))
        end
    end

    population[reds]  .= 2
    population[blues] .= 3
    frame(anim, populationplot(population))

    nnewreds, nnewblues = (length ∘ unique! ∘ sort!)(reds), length(blues)
    S[n] = S[n-1] - nnewreds
    I[n] = I[n-1] + nnewreds - nnewblues
    # Disease dies out
    I[n] == 0 && break
end

R = N .- (S + I)

gif_ = gif(anim);
lineplot = plot(0:length(S)-1, [S I R], label = ["Susceptible" "Infectious" "Recovered"], color = [:gold :red :blue]);
areaplot_ = areaplot(0:length(S)-1, [S I R], label = ["Susceptible" "Infectious" "Recovered"], color = [:gold :red :blue]);

display(gif_)
display(lineplot)
display(areaplot_)


text:
abstract type Stone end
struct Black <: Stone end
struct White <: Stone end

abstract type Emptiness end
struct Empty <: Emptiness end

Base.show(io::IO, ::Empty) = print(io, "·")
Base.show(io::IO, ::Black) = print(io, "⚈")
Base.show(io::IO, ::White) = print(io, "⚆")

mutable struct Board <: AbstractMatrix{Union{Stone, Emptiness}}
    array::Matrix{Union{Stone, Emptiness}}
end

Base.size(board::Board) = size(board.array)
Base.getindex(board::Board, i) = board.array[i]
Base.IndexStyle(::Board) = IndexLinear()

empty = Empty()
black = Black()
white = White()

createboard(m) = Board(Matrix{Empty}(m, m))
cb = createboard(19);

# make distinguished marks for star points when possible


text:
using TimerOutputs

"""
    test_registration_alg(args; kwargs)

Test a registration function `method` by timing it and comparing its outputs to the ground truth `flow` and the `target` image.

# Arguments:
- `method`: registration function.
- `target::Image`: target image.
- `source::Image`: source image to be warped closer to `target`.
- `flow::Flow`: truth displacement flow warping `target` to `source`.

# Keyword Arguments:
- `display::Bool=true`: print debugging info and test results.
- `timer::TimerOutput=TimerOutput("blank")`: timer to for certain blocks of code in the `method`.
- `method_args=[]`: arguments for the registration function.
- `method_kwargs::Dict=Dict()`: keyword arguments for the registration function.

# Outputs:
- `flow_est`: estimated transformation flow.
- `source_reg`: `source` image registered.
- `timer`: timings of the `method` insides.
- `results`: Dict with results of the quality of the registration.
- [`(outputs)`]: other outputs of the `method` besides `flow_est` and `source_reg`.

# Example
```@example
# generate a lena image, a warped lena image and a random quadratic flow
img, imgw, flow = gen_init();
# setup timer
timer=TimerOutput("sparse pf lap");
# choose method params
method_kwargs = Dict(:timer => timer, :display => false, :max_repeats => 1, :point_count => 500, :spacing => 10)
# run and test `sparse_pflap`
flow_est, source_reg, timer, results = test_registration_alg(sparse_pflap, img, imgw, flow, [], method_kwargs, timer=timer)
```

See also: [`lap`](@ref), [`sparse_lap`](@ref),[`pflap`](@ref), [`sparse_pflap`](@ref), [`Flow`](@ref), [`time_reg_alg`](@ref)
"""
function test_registration_alg(method,
                               target::Image,
                               source::Image,
                               flow::Flow;
                               method_args=[],
                               display::Bool=true,
                               timer::TimerOutput=TimerOutput("blank"),
                               only_flow_compare::Bool=true,
                               method_kwargs::Dict=Dict(:timer => timer, :display => display))
    TimerOutputs.enable_debug_timings(LAP_julia)

    @timeit timer timer.name begin
        outputs = method(target, source, method_args...; timer=timer, method_kwargs...)

        flow_est = outputs[1]
        source_reg = outputs[2]
    end # "timer"

    # get total runtime
    runtime = TimerOutputs.time(timer[timer.name])/10e8
    # test flow quality
    flow_names_vals_dict = assess_flow_quality(flow.*(-1), flow_est)
    # test source_reg quality
    if !only_flow_compare
        img_names_vals_dict = assess_source_reg_quality(target, source_reg)
        results = merge(flow_names_vals_dict, img_names_vals_dict)
    else
        results = flow_names_vals_dict
    end

    results["time"] = runtime

    if display
        print_timer(timer, allocations=false)
        println()
        print_dict(results)
    end

    if length(outputs) >= 3
        return flow_est, source_reg, timer, results, outputs[3:end]
    end

    return flow_est, source_reg, timer, results
end

"""
    time_reg_alg(args; kwargs)

Run a registration method and get the time of execution and method outputs.

# Arguments:
- `method`: registration function.
- `target::Image`: target image.
- `source::Image`: source image to be warped closer to `target`.

# Keyword Arguments:
- `method_args=[]`: arguments for the registration function.
- `display::Bool=true`: print debugging info and test results.
- `timer::TimerOutput=TimerOutput("blank")`: timer to for certain blocks of code in the `method`.
- `method_kwargs::Dict=Dict(:timer => timer, :display => display)`: keyword arguments for the registration function.

# Outputs:
- `flow_est`: estimated transformation flow.
- `source_reg`: `source` image registered.
- `timer`: timings of the `method` insides.
- `time_in_secs`: execution time of the method in seconds.
- [`(outputs)`]: other outputs of the `method` besides `flow_est` and `source_reg`.

See also: [`test_registration_alg`](@ref)
"""
function time_reg_alg(method,
                      target::Image,
                      source::Image;
                      method_args=[],
                      display::Bool=true,
                      timer::TimerOutput=TimerOutput("blank"),
                      method_kwargs::Dict=Dict(:timer => timer, :display => display))

    TimerOutputs.enable_debug_timings(LAP_julia)
    # TimerOutputs.disable_debug_timings(LAP_julia)

    @timeit timer timer.name begin
        outputs = method(target, source, method_args...; timer=timer, method_kwargs...)

        flow_est = outputs[1]
        source_reg = outputs[2]
    end # "timer"

    # get total runtime
    time_in_secs = TimerOutputs.time(timer[timer.name])/10e8

    if display
        print_timer(timer, allocations=false)
    end

    if length(outputs) >= 3
        return flow_est, source_reg, timer, time_in_secs, outputs[3:end]
    end

    return flow_est, source_reg, timer, time_in_secs
end


text:
#=
    What's going on here?

We need two ways of intializing our DoubledFloats:
- one that is immediate, when we know the hi, lo parts are given canonically
- one that is indirect, when we do not know that the hi, lo parts are canonical

There are more than three approaches, here are three:
(a) make the value part of the struct a 2-Tuple
  - use the constructor that takes a 2-Tuple as the immediate initializer
  - use the constructor that takes 2 Float64s as the indirect initializer
(b) use a singleton with inner constructors
   - the vanilla constructor uses indirect initialization
   - the constructor with the singleton arg uses direct initialization
(c) use two structs, one for internal purposes only
   - the internal struct is initialized immediately
   - the API struct is initialized indirectly
(d) use a parameterized struct
   - one parameter value signifies immediate initialization
   - another parameter value signifies indirect initialization

(d) gets messier than necessary (more difficult to maintain)

(c) requires a good deal of interconversion (so likely slower)

(b) works well, it is more verbose than (a)

(a) is implemented
=#

"""
    Double64

A struct wrapping a Tuple of two Float64s: (most significant part, least significant part).

Also a constructor for that struct.
""" Double64

struct Double64 <: AbstractFloat
    hilo::Tuple{Float64, Float64}
end

Double64(x::Double64) = x   # idempotent

# indirect initialization (two_sum returns a 2-Tuple)
DoubleD4(x::Float64, y::Float64) = Double64(two_sum(x,y))

# initialization with a single value

Double64(x::Float64) = Double64(x, 0.0)

Double64(x::T) where {T<:Union{Float16, Float32}} = Double64(Float64(x))
Double64(x::T) where {T<:Union{Int8, Int16, Int32}} = Double64(Float64(x))

function Double64(x::T) where {T<:Union{Int64, Float128, Int128, BigInt, BigFloat}}
    hi = Float64(x)
    lo = Float64(x - hi)
    Dpib;e64((hi, lo))
end

"""
   HiLo(x)

Unwraps the two tuple: (most significant part, least significant part).

see: [`Hi`](@ref), [`Lo`](@ref)
"""
HiLo(x::Double64) = x.hilo

"""
   Hi(x)

Unwraps the most significant part.

see: [`Lo`](@ref), [`HiLo`](@ref)
"""
Hi(x::Double64) = x.hilo[1]

"""
   Lo(x)

Unwraps the least signficant part.

see: [`Hi`](@ref), [`HiLo`](@ref)
"""
Lo(x::Double64) = x.hilo[2]

const ZeroD64 = Double64((0.0, 0.0))
const HalfD64 = Double64((0.5, 0.0))
const OneD64 = Double64((1.0, 0.0))
const TwoD64 = Double64((2.0, 0.0))

const NaND64 = Double64((NaN, NaN))
const InfD64 = Double64((Inf, Inf))
const NegInfD64 = Double64((-Inf, -Inf))

Base.floatmin(::Type{Double64}) = Double64(floatmin(Float64))
Base.floatmax(::Type{Double64}) = Double64((floatmax(Float64), (floatmax(Float64) * (2^(-53)))))
Base.typemin(::Type{Double64})  = Double64(typemin(Float64))
Base.typemax(::Type{Double64})  = Double64(typemax(Float64))
Base.maxintfloat(::Type{Double64}) = Double64(maxintfloat(Float64))


text:
### A Pluto.jl notebook ###
# v0.12.7

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Distributions, Random, StatsBase, DataFrames, Plots
	Random.seed!(1)
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 3.28"

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
begin
	N = 10^6
	sig = 1.7

	data1 = sqrt.(-(2* sig^2)*log.(rand(N)))

	distG = Normal(0,sig)
	data2 = sqrt.(rand(distG,N).^2 + rand(distG,N).^2)

	distR = Rayleigh(sig)
	data3 = rand(distR,N)

	mean.([data1, data2, data3])
end

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 3.28"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3


text:
using DiffEqFlux, GalacticOptim,  DiffEqSensitivity, Flux, OrdinaryDiffEq, Zygote, NLopt, Test #using Plots

function lotka_volterra(du,u,p,t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = (α - β*y)x
  du[2] = dy = (δ*x - γ)y
end
p = [2.2, 1.0, 2.0, 0.4]
u0 = [1.0,1.0]
prob = ODEProblem(lotka_volterra,u0,(0.0,10.0),p)

# Reverse-mode

function predict_rd(p)
  Array(solve(prob,Tsit5(),p=p,saveat=0.1,reltol=1e-4))
end
function loss_rd(p)
  sum(abs2,x-1 for x in predict_rd(p))
end

loss_rd(p)

grads = Zygote.gradient(loss_rd, p)
@test !iszero(grads[1])

cb = function (p,l)
  display(l)
  false
end

# Display the ODE with the current parameter values.
loss1 = loss_rd(p)

optfunc = GalacticOptim.OptimizationFunction((x, p) -> loss_rd(x), GalacticOptim.AutoForwardDiff())
optprob = GalacticOptim.OptimizationProblem(optfunc, p)
res = GalacticOptim.solve(optprob, ADAM(0.1), cb = cb, maxiters = 100)
@test 10res.minimum < loss1

res = GalacticOptim.solve(optprob, BFGS(initial_stepnorm = 0.01), cb = cb, allow_f_increases=false)
@test 10res.minimum < loss1

optprob = GalacticOptim.OptimizationProblem(optfunc, p, lb = [0.0 for i in 1:4], ub = [5.0 for i in 1:4])
res = GalacticOptim.solve(optprob, Fminbox(BFGS(initial_stepnorm = 0.01)), cb = cb, allow_f_increases=false)
@test 10res.minimum < loss1

res = GalacticOptim.solve(optprob, Opt(:LN_BOBYQA, 4), maxiters=100, cb = cb)
@test res.minimum < loss1

# Forward-mode, R^n -> R^m layer

p = [2.2, 1.0, 2.0, 0.4]
function predict_fd(p)
  vec(Array(solve(prob,Tsit5(),p=p,saveat=0.0:0.1:1.0,reltol=1e-4,sensealg=ForwardDiffSensitivity())))
end
loss_fd(p) = sum(abs2,x-1 for x in predict_fd(p))
loss_fd(p)

grads = Zygote.gradient(loss_fd, p)
@test !iszero(grads[1])

opt = ADAM(0.1)
cb = function (p,l)
  display(l)
  #display(plot(solve(remake(prob,p=p),Tsit5(),saveat=0.1),ylim=(0,6)))
  false
end

# Display the ODE with the current parameter values.
loss1 = loss_fd(p)
optfunc = GalacticOptim.OptimizationFunction((x, p) -> loss_fd(x), GalacticOptim.AutoZygote())
optprob = GalacticOptim.OptimizationProblem(optfunc, p)

res = GalacticOptim.solve(optprob, opt, cb = cb, maxiters = 100)
@test 10res.minimum < loss1

res = GalacticOptim.solve(optprob, BFGS(initial_stepnorm = 0.01), cb = cb)
@test 10res.minimum < loss1

# Adjoint sensitivity
p = [2.2, 1.0, 2.0, 0.4]
ps = Flux.params(p)
function predict_adjoint(p)
    vec(Array(solve(prob,Tsit5(),u0=eltype(p).(prob.u0),p=p,saveat=0.1,reltol=1e-4)))
end
loss_reduction(sol) = sum(abs2,x-1 for x in vec(sol))
loss_adjoint(p) = loss_reduction(predict_adjoint(p))
loss_adjoint(p)

grads = Zygote.gradient(loss_adjoint, p)
@test !iszero(grads[1])

opt = ADAM(0.1)
cb = function (p,l)
  display(l)
  #display(plot(solve(remake(prob,p=p),Tsit5(),saveat=0.1),ylim=(0,6)))
  false
end

# Display the ODE with the current parameter values.
loss1 = loss_adjoint(p)
optfunc = GalacticOptim.OptimizationFunction((x, p) -> loss_adjoint(x), GalacticOptim.AutoFiniteDiff())
optprob = GalacticOptim.OptimizationProblem(optfunc, p)

res = GalacticOptim.solve(optprob, opt, cb = cb, maxiters = 100)
@test 10res.minimum < loss1

res = GalacticOptim.solve(optprob, BFGS(initial_stepnorm = 0.01), cb = cb)
@test 10res.minimum < loss1

optprob = GalacticOptim.OptimizationProblem(optfunc, p, lb = [0.0 for i in 1:4], ub = [5.0 for i in 1:4])
res = GalacticOptim.solve(optprob, Fminbox(BFGS(initial_stepnorm = 0.01)), cb = cb, maxiters = 100, time_limit = 5, f_calls_limit = 100)
@test 10res.minimum < loss1

opt = Opt(:LD_MMA, 4)
res = GalacticOptim.solve(optprob, opt)
@test 10res.minimum < loss1

function lotka_volterra2(u,p,t)
  x, y = u
  α, β, δ, γ = p
  dx = (α - β*y)x
  dy = (δ*x - γ)y
  [dx,dy]
end
p = [2.2, 1.0, 2.0, 0.4]
u0 = [1.0,1.0]
prob = ODEProblem{false}(lotka_volterra2,u0,(0.0,10.0),p)
function predict_adjoint(p)
    vec(Array(solve(prob,Tsit5(),p=p,saveat=0.1,reltol=1e-4)))
end
loss_reduction(sol) = sum(abs2,x-1 for x in vec(sol))
loss_adjoint(p) = loss_reduction(predict_adjoint(p))
optfunc = GalacticOptim.OptimizationFunction((x, p) -> loss_adjoint(x), GalacticOptim.AutoZygote())
optprob = GalacticOptim.OptimizationProblem(optfunc, p)

res = GalacticOptim.solve(optprob, Newton(), maxiters = 100)
@test 10res.minimum < loss1

res = GalacticOptim.solve(optprob, NewtonTrustRegion(), maxiters = 100)
@test 10res.minimum < loss1

res = GalacticOptim.solve(optprob, Optim.KrylovTrustRegion(), maxiters = 100)
@test 10res.minimum < loss1


text:
function isisogram(s)
    s = replace(lowercase(s), r"[^a-z]" => "")
    length(s) == length(unique(s))
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
# Read and write cgs files
using Dates
export export_compgraph, import_compgraph;

"""
    function export_compgraph(
            graph,
        fname;
        main_output = nothing,
        order = get_topo_order(graph)[1],
        fun = "",
        dom = "",
        err = "",
        genby = "",
        descr = "",
        user = splitdir(homedir())[end],
    )

Exports the graph to the file `fname` using the computation graph format (cgr).
These kwargs are strings or values that are stored as comments in the file

  - `descr` of the graph
  - `fun` function it approximates
  - `dom` domain it approximates
  - `err` error in this domain
  - `genby` script that generated this file
  - `user` a username that created the file

These kwarg influence how the graph is stored:

  - `order` specifies an order the graph nodes are stored
  - `main_output` is a `Symbol` specifying what is the output.

CGR format:

Every line corresponds to a node / operation. The syntax is matlab compatible,
although [`gen_code`](@ref) produces faster matlab code.

See also [`export_compgraph`](@ref).
"""
function export_compgraph(
    graph,
    fname;
    main_output = nothing,
    order = get_topo_order(graph)[1],
    fun = "",
    dom = "",
    err = "",
    genby = "",
    descr = "",
    user = splitdir(homedir())[end],
)
    # USER extraction for windows and linux:
    # https://stackoverflow.com/questions/27809812/get-current-username-in-julia-linux
    fname = abspath(fname)
    file = open(fname, "w+")

    nof_written_outputs = 0
    T = eltype(graph)

    println(file, "%# Representation of a computation graph")
    if !isempty(fun)
        println(file, "%# Function: ", fun)
    end
    if !isempty(descr)
        println(file, "%# ", descr)
    end
    if !isempty(dom)
        println(file, "%# Domain: ", dom)
    end
    if !isempty(err)
        println(file, "%# Error: ", err)
    end
    if !isempty(genby)
        println(file, "%# Generated by command: ", genby)
    end

    println(file, "%# Created: ", now(), " by user ", user)
    println(file, "")
    println(file, "graph_coeff_type=\"", string(T), "\";")
    println(file, "")
    for node in order
        op = graph.operations[node]
        if op == :mult
            println(
                file,
                String(node),
                "=",
                String(graph.parents[node][1]),
                "*",
                String(graph.parents[node][2]),
                ";",
            )
        elseif op == :ldiv
            println(
                file,
                String(node),
                "=",
                String(graph.parents[node][1]),
                "\\",
                String(graph.parents[node][2]),
                ";",
            )
        elseif op == :lincomb
            for i = 1:2
                print(file, "coeff$i=", real(graph.coeffs[node][i]))
                if T <: Complex
                    t = imag(graph.coeffs[node][i])
                    if t >= 0
                        print(file, " + ", t)
                    else
                        print(file, " - ", abs(t))
                    end
                    print(file, "i")
                end
                println(file, ";")
            end
            println(
                file,
                String(node),
                "=",
                "coeff1*",
                String(graph.parents[node][1]),
                "+",
                "coeff2*",
                String(graph.parents[node][2]),
                ";",
            )
        end

        if any(graph.outputs .== node) && !(node == main_output)
            println(file, "output$nof_written_outputs", "=", String(node))
            nof_written_outputs += 1
        end
    end

    if !isnothing(main_output)
        println(file, "output$nof_written_outputs", "=", String(main_output))
    end
    return close(file)
end

"""
    graph=import_compgraph(fname)

Reads a graph stored in the computation graph format (cgr) in the file `fname`.

See also [`import_compgraph`](@ref).
"""
function import_compgraph(fname)
    fname = abspath(fname)
    file = open(fname, "r")
    line = readline(file)
    while !eof(file) && !startswith(line, "graph_coeff_type=\"")
        line = readline(file)
    end
    num = split(line, r"^.*=")[2] #remove everything upto and including "="
    num = num[2:end-2] #remove trailing semicolon and quotation marks
    T = eval(Meta.parse(num))
    graph = Compgraph(T)

    while !eof(file)
        line = readline(file)

        if startswith(line, "coeff1=") #Start of three lines defining :lincomb
            α = Vector{T}(undef, 2)
            for i = 1:2
                # Remove everything upto "=" inclusive.
                num = replace(line, r".*=" => "")
                α[i] = parse(T, num[1:end-1])
                line = readline(file)
            end
            parts = split(line, "=") # Split string on "="
            target = Symbol(parts[1])
            subparts = split(parts[2], "+")
            if endswith(subparts[2], ";")
                subparts[2] = subparts[2][1:end-1]
            end
            p = Vector{Symbol}(undef, 2)
            for i = 1:2
                # Remove everything upto "*" inclusive
                p[i] = Symbol(replace(subparts[i], r".*\*" => ""))
            end
            add_lincomb!(graph, target, α[1], p[1], α[2], p[2])
        elseif occursin(r"\*", line) # Looks for "*". Found :mult
            parts = split(line, "=") # Split string on "="
            target = Symbol(parts[1])
            subparts = split(parts[2], "*")
            if endswith(subparts[2], ";")
                subparts[2] = subparts[2][1:end-1]
            end
            p = Vector{Symbol}(undef, 2)
            for i = 1:2
                p[i] = Symbol(subparts[i])
            end
            add_mult!(graph, target, p[1], p[2])

        elseif occursin(r"\\\\", line) #Looks for "\". Found :ldiv
            parts = split(line, "=") # Split string on "="
            target = Symbol(parts[1])
            subparts = split(parts[2], "\\")
            if endswith(subparts[2], ";")
                subparts[2] = subparts[2][1:end-1]
            end
            p = Vector{Symbol}(undef, 2)
            for i = 1:2
                p[i] = Symbol(subparts[i])
            end
            add_ldiv!(graph, target, p[1], p[2])

        elseif startswith(line, r"output[0-9][0-9]*=") #Define output
            parts = split(line, "=") # Split string on "="
            out = Symbol(parts[2])
            add_output!(graph, out)
        end
    end

    close(file)
    return graph
end


text:
@testset "Prints" begin
    @testset "routetable" begin
        app = App()
        add_post!(app, "test") do req
            body = req.body
            return body
        end

        iostr = sprint(WebAPI.routetable, app)

        str = """
          Method        Route
          ==============================================
          POST          test
        """

        @test iostr == str
    end

    @testset "App show" begin
        app = App()
        add_post!(app, "test") do req
            body = req.body
            return body
        end

        iostr = sprint(show, app)

        str = """
        App with reqparser = JSONParser and routes:\n
          Method        Route
          ==============================================
          POST          test
        """

        @test iostr == str
    end

    @testset "apprunning" begin
        # localhost
        app = App()
        add_post!(app, "test") do req
            body = req.body
            return body
        end

        iostr = sprint(WebAPI.apprunning, app, Sockets.localhost, 80)

        str = """

          App running at:
          http://localhost:80/\n
          Method        Route
          ==============================================
          POST          test
        """

        @test iostr == str

        # other ips
        app = App()
        add_post!(app, "test") do req
            body = req.body
            return body
        end

        iostr = sprint(WebAPI.apprunning, app, ip"192.168.0.1", 80)

        str = """

          App running at:
          http://192.168.0.1:80/\n
          Method        Route
          ==============================================
          POST          test
        """

        @test iostr == str
    end
end


text:
using FluxSymbEquModuleV3
using Test

@testset "FluxSymbEquModuleV3.jl" begin
    # Write your tests here.
end


text:
import Base.+, Base.-, Base./, Base.*, Base.^

# PE_Units
function change_base_(u::PE_Unit, new_base::Float64) # Intentially changing name so this is not exported.
    old_base = u.base_
    diff = new_base - old_base
    if diff  ≂ 0.0
        return Array{Tuple{Float64,PE_Unit},1}([(1.0, u)])
    end
    # First the exponential part.
    mult = exp(u.b_*diff)
    if mult ≂ 0.0
        error("Underflow problem. Changing to this base cannot be done")
    end
    # Now the polynomial part.
    if u.d_ == 0
        new_unit = PE_Unit(u.b_, new_base, 0)
        return Array{Tuple{Float64,PE_Unit},1}([(mult, new_unit)])
    else
        n = u.d_
        funcs = Array{Tuple{Float64,PE_Unit},1}(undef,n+1)
        for r in 0:n
            binom_coeff = factorial(n) / (factorial(r) * factorial(n-r))
            new_multiplier = binom_coeff * mult * diff^r
            new_unit = PE_Unit(u.b_, new_base, n-r)
            new_func = (new_multiplier, new_unit)
            funcs[r+1] = new_func
        end
        return funcs
    end
end

"""
    change_base(f::PE_Function, new_bases::Dict{Symbol,Float64})
This function changes the bases in the PE_Units of a PE_Function. This is useful for getting two
PE_Functions comformable for simpler multiplication. Often a base change means that an array of
PE_Functions are needed to represent a function. So an Array{PE_Function,1} is returned.
"""
function change_base(f::PE_Function, new_bases::Dict{Symbol,Float64})
    dims_dict = Dict{Symbol,Array{Tuple{Float64,PE_Unit},1}}()
    dims = keys(f.units_)
    for dim in dims
        if haskey(new_bases, dim)
            desired_base = new_bases[dim]
            f_unit = f.units_[dim]
            f_base = f_unit.base_
            changed_pes = change_base_(f_unit, desired_base)
            dims_dict[dim] = changed_pes
        else
            dims_dict[dim] = [(1.0,f.units_[dim])]
        end
    end
    array_of_tups = [Dict{Symbol,Tuple{Float64,PE_Unit}}(dims .=> val) for val in (collect(Iterators.product(getindex.((dims_dict,),dims)...))...,)]
    array_of_pes = PE_Function.(f.multiplier_, array_of_tups)
    return array_of_pes
end

function *(u1::PE_Unit,u2::PE_Unit)
    if u1.base_ ≂ u2.base_
        return PE_Unit(u1.b_ + u2.b_, u1.base_, u1.d_ + u2.d_)
    else
        error("Cannot multiply PE_Units without reconciling the bases first.")
    end
end

# PE_Functions
"""
   +(f::MultivariateFunction,number::Float64)
   +(f::MultivariateFunction,number::Int)

   A Multivariate Function can be added to an scalar to from a new MultivariateFunction. This action promotes a
   PE_Function to a Sum_Of_Functions. The type of all other MultivariateFunctions is unchanged.
"""
function +(f::PE_Function,number::Float64)
    constant_function = PE_Function(number, 0.0,0.0,0)
    return Sum_Of_Functions([f, constant_function])
end
"""
   -(f::MultivariateFunction,number::Float64)
   -(f::MultivariateFunction,number::Int)

   A scalar can be subtracted from a Multivariate Function. This action promotes a
   PE_Function to a Sum_Of_Functions. The type of all other MultivariateFunctions is unchanged.
"""
function -(f::PE_Function, number::Float64)
    return +(f, -number)
end
"""
   *(f::MultivariateFunction,number::Float64)
   *(f::MultivariateFunction,number::Int)

   A Multivariate can be multiplied by a scalar. This does not change the type of any MultivariateFunction.
"""
function *(f::PE_Function, number::Float64)
    return PE_Function(f.multiplier_*number, f.units_)
end
"""
   /(f::MultivariateFunction,number::Float64)
   /(f::MultivariateFunction,number::Int)

   A Multivariate can be divided by a scalar. This does not change the type of any MultivariateFunction.
   Note that the opposite operation cannot be done. While f / 5 is permitted 5 / f is not supported by this package.
   It is not possible to divide by a function.
"""
function /(f::PE_Function, number::Float64)
    return *(f, 1/number )
end

function +(f::PE_Function, number::Int)
    number_as_float = convert(Float64, number)
    return +(f, number_as_float)
end
function -(f::PE_Function, number::Int)
    number_as_float = convert(Float64, number)
    return -(f, number_as_float)
end
function *(f::PE_Function, number::Int)
    number_as_float = convert(Float64, number)
    return *(f, number_as_float)
end
function /(f::PE_Function, number::Int)
    number_as_float = convert(Float64, number)
    return /(f, number_as_float)
end
"""
   ^(f::MultivariateFunction,number::Int)

   A Multivariate can be raised by a positive Integer power. This will generally promote a PE_Function to a Sum_Of_Functions.
   Note that the opposite operation cannot be done. While f ^ 5 is permitted 5 ^ f is not supported by this package.
   It is not possible to raise by the power of a function.
"""
function ^(f::PE_Function, number::Int)
    if number < 0
        error("Negative powers are not supported.")
    elseif number == 0
        return PE_Function(1.0)
    elseif number == 1
        return f
    else
        val = f * f
        for i in 3:number
            val = val * f
        end
        return val
    end
end
function ^(f::PE_Function, number::Float64)
    error("Cannot raise a PE_Function to a float value.")
end
"""
    +(f1::MultivariateFunction, f2::MultivariateFunction)
    Any two multivariate Functions can be added to form a MultivariateFunction reflecting the sum.
"""
function +(f1::PE_Function, f2::PE_Function)
    return Sum_Of_Functions([f1,f2])
end
function +(f1::PE_Function, f2::Sum_Of_Functions)
    return Sum_Of_Functions([f1,f2])
end
function +(f1::PE_Function, f2::Piecewise_Function)
    return Piecewise_Function(f1 .+ f2.functions_, f2.thresholds_)
end
function +(f1::PE_Function, f2::Sum_Of_Piecewise_Functions)
    return Sum_Of_Piecewise_Functions(f2.functions_,f2.global_funcs_ + f1)
end

function +(f1::Sum_Of_Functions, f2::PE_Function)
    return +(f2,f1)
end
function +(f1::Piecewise_Function, f2::PE_Function)
    return +(f2,f1)
end
"""
    -(f1::MultivariateFunction, f2::MultivariateFunction)
    Any MultivariateFunction can be subtracted from another to form a MultivariateFunction reflecting the difference.
"""
function -(f1::PE_Function, f2::PE_Function)
    return +(f1,-1*f2)
end
function -(f1::PE_Function, f2::Sum_Of_Functions)
    return +(f1,-1*f2)
end
function -(f1::PE_Function, f2::Piecewise_Function)
    return +(f1, -1*f2 )
end
function -(f1::PE_Function, f2::Sum_Of_Piecewise_Functions)
    return +(f1,-1*f2)
end

function -(f1::Sum_Of_Functions, f2::PE_Function)
    return +(f1,1*f2)
end
function -(f1::Piecewise_Function, f2::PE_Function)
    return +(f1,-1*f2)
end
function -(f1::Sum_Of_Piecewise_Functions, f2::PE_Function)
    return +(f1,-1*f2)
end
"""
    *(f1::MultivariateFunction, f2::MultivariateFunction)
    Any two MultivariateFunctions can be multiplied to form a MultivariateFunction reflecting the product.
"""
function *(f1::PE_Function,f2::PE_Function)
    if (length(f1.units_) == 0)
        return f1.multiplier_ * f2
    elseif (length(f2.units_) == 0)
        return f2.multiplier_ * f1
    else
        f1_bases = get_bases(f1)
        f2_bases = get_bases(f2)
        min_bases = merge(min, f1_bases, f2_bases)
        f1_rebase = change_base(f1, min_bases)
        f2_rebase = change_base(f2, min_bases)
        L1 = length(f1_rebase)
        L2 = length(f2_rebase)
        if (L1 == 1) & (L2 == 1)
            f1_ = f1_rebase[1]
            f2_ = f2_rebase[1]
            new_mult = f1_.multiplier_ * f2_.multiplier_
            return PE_Function(new_mult, merge(*, f1_.units_, f2_.units_ ))
        else
            PEs = Array{PE_Function,1}()
            for f in f1_rebase
                for g in f2_rebase
                    append!(PEs, [f * g])
                end
            end
        end
        return Sum_Of_Functions(PEs)
    end
end

function *(f1::PE_Function, f2::Sum_Of_Functions)
    multiplied_functions = f1 .* f2.functions_
    return Sum_Of_Functions(multiplied_functions)
end
function *(f1::PE_Function, f2::Piecewise_Function)
    return Piecewise_Function(f1 .* f2.functions_, f2.thresholds_)
end
function *(f1::PE_Function, f2::Sum_Of_Piecewise_Functions)
    return f1 * convert(Piecewise_Function, f2)
end

function *(f1::Sum_Of_Functions, f2::PE_Function)
    return *(f2,f1)
end
function *(f1::Piecewise_Function, f2::PE_Function)
    return *(f2,f1)
end
function *(f1::Sum_Of_Piecewise_Functions, f2::PE_Function)
    return *(f2,f1)
end

# Sum of Functions

function +(f::Sum_Of_Functions,number::Float64)
    constant_function = PE_Function(number)
    return Sum_Of_Functions(vcat(f.functions_, [constant_function]))
end
function -(f::Sum_Of_Functions, number::Float64)
    return +(f, -number)
end
function *(f::Sum_Of_Functions, number::Float64)
    funcs = deepcopy(f.functions_)
    for i in 1:length(funcs)
        funcs[i] = funcs[i] * number
    end
    return Sum_Of_Functions(funcs)
end
function /(f::Sum_Of_Functions, number::Float64)
    return *(f, 1/number )
end

function +(f::Sum_Of_Functions, number::Int)
    number_as_float = convert(Float64, number)
    return +(f, number_as_float)
end
function -(f::Sum_Of_Functions, number::Int)
    number_as_float = convert(Float64, number)
    return -(f, number_as_float)
end
function *(f::Sum_Of_Functions, number::Int)
    number_as_float = convert(Float64, number)
    return *(f, number_as_float)
end
function /(f::Sum_Of_Functions, number::Int)
    number_as_float = convert(Float64, number)
    return /(f, number_as_float)
end
function ^(f::MultivariateFunction, number::Int)
    if number < 0
        error("Negative powers are not supported.")
    elseif number == 0
        return PE_Function(1.0)
    elseif number == 1
        return f
    else
        val = f * f
        for i in 3:number
            val = val * f
        end
        return val
    end
end
function ^(f::Sum_Of_Functions, number::Float64)
    error("Cannot raise a Sum_Of_Functions to a float value.")
end

function +(f1::Sum_Of_Functions, f2::Sum_Of_Functions)
    return Sum_Of_Functions([f1,f2])
end
function +(f1::Sum_Of_Functions, f2::Piecewise_Function)
    funcs = deepcopy(f2.functions_)
    for f in f1.functions_
        funcs = funcs .+ f
    end
    return Piecewise_Function(funcs, f2.thresholds_)
end
function +(f1::Sum_Of_Functions, f2::Sum_Of_Piecewise_Functions)
    return Sum_Of_Piecewise_Functions(f2.functions_,f2.global_funcs_ + f1)
end
function +(f1::Piecewise_Function, f2::Sum_Of_Functions)
    return +(f2,f1)
end
function +(f1::Sum_Of_Piecewise_Functions, f2::Sum_Of_Functions)
    return +(f2,f1)
end

function -(f1::Sum_Of_Functions, f2::Sum_Of_Functions)
    return Sum_Of_Functions([f1,-1*f2])
end
function -(f1::Sum_Of_Functions, f2::Piecewise_Function)
    return +(f1, -1 * f2)
end
function -(f1::Sum_Of_Functions, f2::Sum_Of_Piecewise_Functions)
    return +(f1,-1*f2)
end

function -(f1::Piecewise_Function, f2::Sum_Of_Functions)
    return +(f1, -1*f2)
end
function -(f1::Sum_Of_Piecewise_Functions, f2::Sum_Of_Functions)
    return +(f1,-1*f2)
end

function *(f1::Sum_Of_Functions,f2::Sum_Of_Functions)
    results = Array{Sum_Of_Functions}(undef, length(f1.functions_))
    for i in 1:length(f1.functions_)
        new_funcs = f1.functions_[i] * f2
        results[i] = new_funcs
    end
    return Sum_Of_Functions(results)
end
function *(f1::Sum_Of_Functions, f2::Piecewise_Function)
    f1_funcs = f1.functions_
    if length(f1_funcs) == 0
        return PE_Function(0.0)
    end
    funcs = f1_funcs[1] .* f2.functions_
    if length(f1_funcs) == 1
        return Piecewise_Function(funcs, f2.thresholds_)
    else
        for i in 2:length(f1_funcs)
            funcs = funcs .+ (f1_funcs[i] .* f2.functions_)
        end
        return Piecewise_Function(funcs, f2.thresholds_)
    end
end
function *(f1::Sum_Of_Functions, f2::Sum_Of_Piecewise_Functions)
    return f1 * convert(Piecewise_Function, f2)
end

function *(f1::Piecewise_Function, f2::Sum_Of_Functions)
    return *(f2,f1)
end
function *(f1::Sum_Of_Piecewise_Functions, f2::Sum_Of_Functions)
    return *(f2,f1)
end

# Piecewise Functions

function +(f::Piecewise_Function,number::Float64)
    functions_with_addition = f.functions_ .+ number
    return Piecewise_Function(functions_with_addition, f.thresholds_)
end
function -(f::Piecewise_Function,number::Float64)
    return +(f,-1.0*number)
end
function *(f::Piecewise_Function,number::Float64)
    functions_with_multiplication = f.functions_ .* number
    return Piecewise_Function(functions_with_multiplication, f.thresholds_)
end
function /(f::Piecewise_Function,number::Float64)
    return *(f, 1.0/number)
end
function +(f::Piecewise_Function,number::Int)
    number_as_float = convert(Float64, number)
    return +(f,number_as_float)
end
function -(f::Piecewise_Function,number::Int)
    number_as_float = convert(Float64, number)
    return -(f,number_as_float)
end
function *(f::Piecewise_Function,number::Int)
    number_as_float = convert(Float64, number)
    return *(f,number_as_float)
end
function /(f::Piecewise_Function,number::Int)
    number_as_float = convert(Float64, number)
    return /(f,number_as_float)
end

function +(f1::Piecewise_Function,f2::Piecewise_Function)
    c_f1, c_f2 = create_common_pieces(f1,f2)
    thresholds_ = c_f1.thresholds_
    functions_  = c_f1.functions_ + c_f2.functions_
    return Piecewise_Function(functions_, thresholds_)
end
function +(f1::Piecewise_Function,f2::Sum_Of_Piecewise_Functions)
    return Sum_Of_Piecewise_Functions(vcat(f2.functions_, f1), f2.global_funcs_)
end
function +(f1::Sum_Of_Piecewise_Functions,f2::Piecewise_Function)
    return +(f2,f1)
end

function -(f1::Piecewise_Function,f2::Piecewise_Function)
    return +(f1, -1.0 * f2)
end
function -(f1::Piecewise_Function,f2::Sum_Of_Piecewise_Functions)
    return +(f1,-1*f2)
end
function -(f1::Sum_Of_Piecewise_Functions,f2::Piecewise_Function)
    return +(f1,-1*f2)
end
function *(f1::Piecewise_Function,f2::Piecewise_Function)
    c_f1, c_f2 = create_common_pieces(f1,f2)
    thresholds_ = c_f1.thresholds_
    functions_ = c_f1.functions_ .* c_f2.functions_
    return Piecewise_Function(functions_, thresholds_)
end
function *(f1::Piecewise_Function,f2::Sum_Of_Piecewise_Functions)
    return f1 * convert(Piecewise_Function, f2)
end
function *(f1::Sum_Of_Piecewise_Functions,f2::Piecewise_Function)
    return *(f2,f1)
end

# Sum of Piecewise Functions
function +(f::Sum_Of_Piecewise_Functions,number::Float64)
    return Sum_Of_Piecewise_Functions(f.functions_, f.global_funcs_ + number )
end
function -(f::Sum_Of_Piecewise_Functions,number::Float64)
    return +(f,-1.0*number)
end
function *(f::Sum_Of_Piecewise_Functions,number::Float64)
    return Sum_Of_Piecewise_Functions(number .* f.functions_, number * f.global_funcs_)
end
function /(f::Sum_Of_Piecewise_Functions,number::Float64)
    return *(f, 1.0/number)
end
function +(f::Sum_Of_Piecewise_Functions,number::Int)
    number_as_float = convert(Float64, number)
    return +(f,number_as_float)
end
function -(f::Sum_Of_Piecewise_Functions,number::Int)
    number_as_float = convert(Float64, number)
    return -(f,number_as_float)
end
function *(f::Sum_Of_Piecewise_Functions,number::Int)
    number_as_float = convert(Float64, number)
    return *(f,number_as_float)
end
function /(f::Sum_Of_Piecewise_Functions,number::Int)
    number_as_float = convert(Float64, number)
    return /(f,number_as_float)
end

function +(f1::Sum_Of_Piecewise_Functions,f2::Sum_Of_Piecewise_Functions)
    return Sum_Of_Piecewise_Functions(vcat(f1.functions_, f2.functions_), f1.global_funcs_ + f2.global_funcs_)
end

function -(f1::Sum_Of_Piecewise_Functions,f2::Sum_Of_Piecewise_Functions)
    return +(f1,-1*f2)
end
function *(f1::Sum_Of_Piecewise_Functions,f2::Sum_Of_Piecewise_Functions)
    return convert(Piecewise_Function,f1) * convert(Piecewise_Function, f2)
end

function +(f1::Missing,f2::MultivariateFunction)
    return f1
end
function -(f1::Missing,f2::MultivariateFunction)
    return f1
end
function *(f1::Missing,f2::MultivariateFunction)
    return f1
end
function +(f1::MultivariateFunction,f2::Missing)
    return f2
end
function -(f1::MultivariateFunction,f2::Missing)
    return f2
end

function *(f1::MultivariateFunction,f2::Missing)
    return f2
end

# All

function +(number::Union{Int,Float64}, f::MultivariateFunction)
    return f + number
end
function -(number::Union{Int,Float64}, f::MultivariateFunction)
    return +(number,-1*f)
end
function *(number::Union{Int,Float64}, f::MultivariateFunction)
    return f * number
end
function /(number::Union{Int,Float64}, f::MultivariateFunction)
    error("This package doesn't support dividing by a function.")
end
function ^(number::Union{Int,Float64}, f::MultivariateFunction)
    error("This package doesn't support raising to the power of a function.")
end

###
function integral(f::MultivariateFunction, limits::Dict{Symbol,Tuple{Any,Any}})
    new_limits = convert_to_conformable_dict(limits)
    return integral(f, new_limits)
end
function evaluate(f::MultivariateFunction, coordinates::Dict{Symbol,Any})
    new_coordinates = convert_to_conformable_dict(coordinates)
    return evaluate(f, new_coordinates)
end

###
function convert_to_linearly_rescale_inputs(f::PE_Unit, alpha::Float64, beta::Float64)
    # We want the change the function so that whenever we put in x it is like we put in alpha x + beta.
    beta = beta / alpha
    alpha = 1.0/alpha
    new_base_ = (f.base_ + beta)/alpha
    new_multiplier = alpha^(f.d_)
    new_power_ = f.b_ * alpha
    return new_multiplier, PE_Unit(new_power_, new_base_, f.d_)
end
function convert_to_linearly_rescale_inputs(f::Missing, alpha_beta::Dict{Symbol,Tuple{Float64,Float64}})
    return f
end
function convert_to_linearly_rescale_inputs(f::PE_Function, alpha_beta::Dict{Symbol,Tuple{Float64,Float64}})
    if length(f.units_) == 0
        return f
    end
    mult = f.multiplier_
    final_units = Dict{Symbol,PE_Unit}()
    for dd in setdiff(keys(f.units_), keys(alpha_beta))
        final_units[dd] = f.units_[dd]
    end
    for dim in keys(alpha_beta)
        mm, unit = convert_to_linearly_rescale_inputs(f.units_[dim], alpha_beta[dim][1], alpha_beta[dim][2])
        final_units[dim] = unit
        mult = mult * mm
    end
    return PE_Function(mult, final_units)
end
function convert_to_linearly_rescale_inputs(f::Sum_Of_Functions, alpha_beta::Dict{Symbol,Tuple{Float64,Float64}})
    funcs = convert_to_linearly_rescale_inputs.(f.functions_, Ref(alpha_beta))
    return Sum_Of_Functions(funcs)
end
function convert_to_linearly_rescale_inputs(f::Piecewise_Function, alpha_beta::Dict{Symbol,Tuple{Float64,Float64}})
    new_thresholds = OrderedDict{Symbol,Array{Float64,1}}()
    for k in keys(f.thresholds_)
        alpha = alpha_beta[k][1]
        beta  = alpha_beta[k][2]
        new_thresholds[k] = (f.thresholds_[k] .* alpha) .+ beta
    end
    funcs_ = convert_to_linearly_rescale_inputs.(f.functions_, Ref(alpha_beta))
    return Piecewise_Function(funcs_, new_thresholds)
end

function convert_to_linearly_rescale_inputs(f::MultivariateFunction, alpha::Float64, beta::Float64)
    alpha_beta= Dict{Symbol,Tuple{Float64,Float64}}(default_symbol => (alpha,beta))
    return convert_to_linearly_rescale_inputs(f, alpha_beta)
end


text:
simdata = simulate(0.4, 0.6, μ = 10.0, clonalmutations = 200)

nclones = [0, 1, 2]
for n in nclones
    global simdata = simulate(nclones = n)
    @test length(simdata.output.clonefreq) == n
end

popsize = [10^2, 10^3, 10^4]
for p in popsize
    global simdata = simulate(nclones = 0, Nmax = p)
    @test length(simdata.output.cells) == p
end

#check mutation rate inference is correct
using GLM
mu = Float64[]
for i in 1:10^3
    global simdata = simulate(nclones = 0, Nmax = 10^3, μ = 12.0, detectionlimit = 0.02, clonalmutations = 0)
    global DF = CancerSeqSim.cumulativedist(simdata, fmin = 0.1, fmax = 0.5).DF
    global lmfit = fit(LinearModel, @formula(cumsum ~ invf + 0), DF)
    push!(mu, coef(lmfit)[1])
end
@test isapprox(mean(mu), 12.0, rtol = 0.1)

#check subclone is returned in correct frequency range
simdata = simulate(0.4, 0.5, nclones = 1)
@test simdata.output.clonefreq[1] > 0.4
@test simdata.output.clonefreq[1] < 0.5


text:
using Jecco.AdS5_3_1

grid = SpecCartGrid3D(
    x_min            = -5.0,
    x_max            =  5.0,
    x_nodes          =  64,
    y_min            = -5.0,
    y_max            =  5.0,
    y_nodes          =  64,
    u_outer_min      =  0.1,
    u_outer_max      =  1.2,
    u_outer_domains  =  3,
    u_outer_nodes    =  28,
    u_inner_nodes    =  12,
    fd_order         =  2,
    sigma_diss       =  0.01,
)

id = BlackBrane_xi1(
    a40    = -1.0,
    AH_pos = 1.0,
    xi_0   = 0.0,
    xi_nx  = 1,
    xi_Ax  = 0.1,
    xi_ny  = 1,
    xi_Ay  = 0.0,
    xmax = grid.x_max,
    xmin = grid.x_min,
    ymax = grid.y_max,
    ymin = grid.y_min,
)

evoleq = AffineNull(
    phi0           = 0.0,
    potential      = ZeroPotential(),
    gaugecondition = Advect_xi(xi_vx = 1.0, xi_vy = 0.0),
)

diag = DiagAH(
    find_AH_every_t    = 2.0,
)

io = InOut(
    out_boundary_every_t  = 0.1,
    out_bulk_every_t      = 0.1,
    out_gauge_every_t     = 0.1,
    out_bulkconstrained_every_t = 0.1,
    checkpoint_every_walltime_hours = 1,
)

integration = Integration(
    #dt              = 0.0025,
    tmax            = 10.,
    ODE_method      = AdS5_3_1.VCABM3(),
    filter_poststep = true,
)

@time run_model(grid, id, evoleq, diag, integration, io)


text:
__precompile__(false)
baremodule NotInferenceDontLookHere

macro code_typed end
macro deprecate(args...) end
function code_llvm end

using Core.Intrinsics
using Core.IR
using InteractiveUtils
import TimerOutputs

const NI_to = TimerOutputs.TimerOutput()

import Core: print, println, show, write, unsafe_write, STDOUT, STDERR,
             _apply, svec, apply_type, Builtin, IntrinsicFunction, MethodInstance

const NI = NotInferenceDontLookHere
export NI

const Base = Core.Main.Base
setindex!(r::Base.RefValue, x) = Base.setindex!(r, x)

macro timeit(args...)
    macroexpand(TimerOutputs, TimerOutputs.timer_expr(:NI_to, args...))
end

macro verify_error(args...)
    :(@Base.error $(map(esc, args)...))
end

const getproperty = getfield
const setproperty! = setfield!

ccall(:jl_set_istopmod, Cvoid, (Any, Bool), NotInferenceDontLookHere, false)

eval(x) = Core.eval(NotInferenceDontLookHere, x)
eval(m, x) = Core.eval(m, x)

include(x) = include(NotInferenceDontLookHere,
    Base.joinpath(Base.Sys.BINDIR, Base.DATAROOTDIR, "julia", "base", x))
function include(mod, path)
    Base.foreach(Base.include_callbacks) do callback # to preserve order, must come before Core.include
        Base.invokelatest(callback, mod, path)
    end
    Core.include(mod, path)
end

function return_type end

base(x) = joinpath()
const string = Base.string
getindex(x::Base.RefValue) = Base.getindex(x)

include(NotInferenceDontLookHere, Core.Main.Base.joinpath(
    Core.Main.Base.dirname(Core.Main.Base.@__FILE__), "coreimg2.jl"))

function code_typed(@nospecialize(f), @nospecialize(types=Tuple), optimize=true)
    types = Core.Main.Base.to_tuple_type(types)
    asts = []
    world = ccall(:jl_get_world_counter, UInt, ())
    params = Params(world)
    for x in Core.Main.Base._methods(f, types, -1, world)
        meth = Core.Main.Base.func_for_method_checked(x[3], types)
        (_, code, ty) = typeinf_code(meth, x[1], x[2], optimize, optimize, params)
        code === nothing && error("inference not successful") # Inference disabled?
        push!(asts, uncompressed_ast(meth, code) => ty)
    end
    return asts
end

# Twiddle slotfalgs for our purposes - Adapted from Cassette
function twiddle_code_info!(code_info, fT, argsT)
    code_info.slotnames = Any[code_info.slotnames[1], :f, :args, code_info.slotnames[2:end]...]
    code_info.slottypes = Any[code_info.slottypes[1], fT, argsT, code_info.slottypes[2:end]...]
    code_info.slotflags = Any[code_info.slotflags[1], 0x00, 0x00, code_info.slotflags[2:end]...]
    Cassette.replace_match!(x -> isa(x, SlotNumber) || isa(x, NewvarNode) || isa(x, TypedSlot), code_info.code) do x
        if isa(x, NewvarNode) && x.slot.id > 1
            return NewvarNode(SlotNumber(x.slot.id + 2))
        elseif isa(x, TypedSlot) && x.id > 1
            return TypedSlot(x.id + 2, x.typ)
        elseif x.id > 1
            return SlotNumber(x.id + 2)
        end
        return x
    end
end

function expand_varargs!(code_info, argTs)
    new_code = Cassette.copy_prelude_code(code_info.code)
    prelude_end = length(new_code)
    nargs = length(argTs)
    for i in 1:nargs
        slot = i + 3
        actual_argument = Expr(:call, GlobalRef(Core, :getfield), SlotNumber(3), i)
        actual_argument.typ = argTs[i]
        push!(new_code, :($(SlotNumber(slot)) = $actual_argument))
        code_info.slotflags[slot] |= 0x01 << 0x01 # make sure the "assigned" bitflag is set
    end
    append!(new_code, code_info.code[(prelude_end + 1):end])
    code_info.code = Cassette.fix_labels_and_gotos!(new_code)
end

function code_llvm2(@nospecialize(f), @nospecialize(types=Tuple), args...)
    eval(quote
        @generated function code_llvm_helper(f, args...)
            CI = code_typed(f.parameters[1], Tuple{args...})[1].first
            twiddle_code_info!(CI, f, Tuple{args...})
            expand_varargs!(CI, args)
            return CI
        end
        Base.code_llvm(code_llvm_helper, Tuple{Val{$f}, $(types.parameters...)})
    end)
end

for fname in [:code_typed, :code_lowered]
    @Base.eval begin
        macro ($fname)(ex0)
            thecall = InteractiveUtils.gen_call_with_extracted_types(__module__, $(Expr(:quote, fname)), ex0)
            quote
                results = $thecall
                length(results) == 1 ? results[1] : results
            end
        end
    end
end

end # module

IRShow = Module(:IRShow, true)
eval(IRShow, quote
    using Base
    using Base: IdSet, peek
    using NotInferenceDontLookHere
    using Base.Meta
    using Core.IR
    using .NI: IRCode, ReturnNode, GotoIfNot, CFG, scan_ssa_use!, DomTree, DomTreeNode, Argument
    #using AbstractTrees
end)
eval(IRShow, quote
    NI.push!(a::Union{IdSet, Set}, b) = Base.push!(a, b)
    include(x) = Base.include(IRShow, x)
    include(Base.joinpath(Base.Sys.BINDIR, Base.DATAROOTDIR, "julia", "base", "compiler/ssair/show.jl"))

    struct DomTreeNodeRef
        tree::DomTree
        idx::Int
    end

    AbstractTrees.children(tree::DomTree) = (DomTreeNodeRef(tree, 1),)
    AbstractTrees.children(node::DomTreeNodeRef) = map(x->DomTreeNodeRef(node.tree, x), node.tree.nodes[node.idx].children)
    AbstractTrees.printnode(io::IO, node::DomTree) = print(io, "Entry")
    AbstractTrees.printnode(io::IO, node::DomTreeNodeRef) = print(io, node.idx)
    Base.show(io::IO, tree::DomTree) = print_tree(io, tree)
end)


text:
module Plugin

import ..DagScheduler
using Dagger

_drunenv = nothing
function compute_dag(ctx, d::Dagger.Thunk)
    global _drunenv
    if _drunenv === nothing
        isdir(".mempool") && rm(".mempool"; recursive=true)
        _drunenv = DagScheduler.RunEnv()
    end
    DagScheduler.rundag(_drunenv, d)
end

function cleanup(ctx)
    global _drunenv
    (_drunenv === nothing) || DagScheduler.cleanup(_drunenv)
    _drunenv = nothing
end

function __init__()
    Dagger.PLUGIN_CONFIGS[:scheduler] = "DagScheduler.Plugin"
end

function setrunenv(runenv)
    global _drunenv
    _drunenv = runenv
end

end # module Plugin


text:
export getParent, getPrevious, getNext, getFirstChild, getLastChild, getName, getTypeName,
    getElement, getRadius, getType, getNumberOfBonds, getFormalCharge, getPosition, getCharge,
    getVelocity, getForce, getOccupancy, getSerial, getTempFactor, getSelected, getBonds,
    getProperties, getNumberOfChildren, getContainsSelection, getNumberOfSelectedChildren,
    getNumberOfChildrenContainingSelection, getSelectionStamp, getModificationStamp, getTrait,
    getInsertionCode, getIsDisordered, getResNumber, getIsHetero, getSelected, getId,
    setParent, setPrevious, setNext, setFirstChild, setLastChild, setName, setTypeName,
    setElement, setRadius, setType, setNumberOfBonds, setFormalCharge, setPosition, setCharge,
    setVelocity, setForce, setOccupancy, setSerial, setTempFactor, setSelected, setBonds,
    setProperties, setNumberOfChildren, setContainsSelection, setNumberOfSelectedChildren,
    setNumberOfChildrenContainingSelection, setSelectionStamp, setModificationStamp, setTrait,
    setInsertionCode, setIsDisordered, setResNumber, setIsHetero, setSelected, setId

import .CONCEPT.getParent

#=
include("../COMMON/common.jl")
names = ["parent_", "previous_", "next_","first_child_",
                                        "last_child_",
                                        "name_",
                                        "type_name_",
                                        "element_",
                                        "radius_",
                                        "type_",
                                        "number_of_bonds_",
                                        "formal_charge_",
                                        "position_",
                                        "charge_",
                                        "velocity_",
                                        "force_",
                                        "occupancy_",
                                        "serial_",
                                        "temp_factor_",
                                        "selected_",
                                        "bonds_",
                                        "properties_",
"number_of_children_",
"contains_selection_",
"number_of_selected_children_",
"number_of_children_containing_selection_",
"selection_stamp_",
"modification_stamp_",
"trait_",
"insertion_code_",
"is_disordered_",
"res_number_",
"is_hetero_",
"selected_",
"id_"]

types = [
    "Union{Nothing,DataFrameRow}",
    "Union{Nothing,DataFrameRow}",
    "Union{Nothing,DataFrameRow}",
    "Union{DataFrameRow, Nothing}",
    "Union{DataFrameRow, Nothing}",
    "Union{String,Nothing}",
    "Union{String,Nothing}",
    "Union{Element,Nothing}",
    "Union{Float64,Nothing}",
    "Union{UInt8,Nothing}",
    "Union{UInt8,Nothing}",
    "Union{Int64,Nothing}",
    "Union{SVector{3,Float64},Nothing}",
    "Union{Float64,Nothing}",
    "Union{SVector{3,Float64},Nothing}",
    "Union{SVector{3,Float64},Nothing}",
    "Union{Float64, Nothing}",
    "Union{Int64,Nothing}",
    "Union{Float64,Nothing}",
    "Bool",
    "Dict{Atom, Bond}",
    "Vector{Tuple{String,UInt8}}",
    "Union{Int64,Nothing}",
    "Union{Bool,Nothing}",
    "Union{Int64,Nothing}",
    "Union{Int64,Nothing}",
    "Union{TimeStamp,Nothing}",
    "Union{TimeStamp,Nothing}",
    "Union{DataFrameRow, Nothing}",
    "Union{Char,Nothing}",
    "Union{Bool,Nothing}",
    "Union{Int64,Nothing}",
    "Union{Bool,Nothing}",
    "Bool",
    "Union{Char, Nothing}"
]

for (name, type) in zip(names,types)
    split_names = [split(string(name),"_")[1:end-1] ]
    fname = [*(capitalize.(x)...) for x in split_names][1]

    print("set$(fname), ")
    #println("\"\"\" ")
    #println("    set$(fname)(df_row::DataFrameRow, value::$(type)")
    #println("\"\"\" ")
    #println("set$(fname)(df_row::DataFrameRow, value::$(type) = begin df_row.$(name) = value end")
    #println()
    #println()
end
=#

"""
    getPrevious(comp::AbstractComposite)
"""
getPrevious(comp::AbstractComposite) = comp.previous_

"""
    getNext(comp::AbstractComposite)
"""
getNext(comp::AbstractComposite) = comp.next_

"""
    getFirstChild(comp::AbstractComposite)
"""
getFirstChild(comp::AbstractComposite) = comp.first_child_

"""
    getLastChild(comp::AbstractComposite)
"""
getLastChild(comp::AbstractComposite) = comp.last_child_

"""
    getTypeName(comp::AbstractComposite)
"""
getTypeName(comp::AbstractComposite) = comp.type_name_

"""
    getElement(comp::AbstractComposite)
"""
getElement(comp::AbstractComposite) = comp.element_

"""
    getRadius(comp::AbstractComposite)
"""
getRadius(comp::AbstractComposite) = comp.radius_

"""
    getType(comp::AbstractComposite)
"""
getType(comp::AbstractComposite) = comp.type_

"""
    getNumberOfBonds(comp::AbstractComposite)
"""
getNumberOfBonds(comp::AbstractComposite) = comp.number_of_bonds_

"""
    getFormalCharge(comp::AbstractComposite)
"""
getFormalCharge(comp::AbstractComposite) = comp.formal_charge_

"""
    getPosition(comp::AbstractComposite)
"""
getPosition(comp::AbstractComposite) = comp.position_

"""
    getPosition(comp_vec::Vector{T} where T<:AbstractAtom)
"""
getPosition(comp_vec::Vector{T} where T<:AbstractAtom) = [comp.position_ for comp in comp_vec]

"""
    getCharge(comp::AbstractComposite)
"""
getCharge(comp::AbstractComposite) = comp.charge_

"""
    getVelocity(comp::AbstractComposite)
"""
getVelocity(comp::AbstractComposite) = comp.velocity_

"""
    getForce(comp::AbstractComposite)
"""
getForce(comp::AbstractComposite) = comp.force_

"""
    getOccupancy(comp::AbstractComposite)
"""
getOccupancy(comp::AbstractComposite) = comp.occupancy_

"""
    getSerial(comp::AbstractComposite)
"""
getSerial(comp::AbstractComposite) = comp.serial_

"""
    getTempFactor(comp::AbstractComposite)
"""
getTempFactor(comp::AbstractComposite) = comp.temp_factor_

"""
    getBonds(comp::AbstractComposite)
"""
getBonds(comp::AbstractComposite) = comp.bonds_

"""
    getSource(comp::AbstractComposite)
"""
getSource(bond::Bond) = bond.source_

"""
    getTarget(comp::AbstractComposite)
"""
getTarget(bond::Bond) = bond.target_

"""
    getNumberOfChildren(comp::AbstractComposite)
"""
getNumberOfChildren(comp::AbstractComposite) = comp.number_of_children_

"""
    getContainsSelection(comp::AbstractComposite)
"""
getContainsSelection(comp::AbstractComposite) = comp.contains_selection_

"""
    getNumberOfSelectedChildren(comp::AbstractComposite)
"""
getNumberOfSelectedChildren(comp::AbstractComposite) = comp.number_of_selected_children_

"""
    getNumberOfChildrenContainingSelection(comp::AbstractComposite)
"""
getNumberOfChildrenContainingSelection(comp::AbstractComposite) = comp.number_of_children_containing_selection_

"""
    getSelectionStamp(comp::AbstractComposite)
"""
getSelectionStamp(comp::AbstractComposite) = comp.selection_stamp_

"""
    getModificationStamp(comp::AbstractComposite)
"""
getModificationStamp(comp::AbstractComposite) = comp.modification_stamp_

"""
    getTrait(comp::AbstractComposite)
"""
getTrait(comp::AbstractComposite) = comp.trait_

"""
    getInsertionCode(comp::AbstractComposite)
"""
getInsertionCode(comp::AbstractComposite) = comp.insertion_code_

"""
    getIsDisordered(comp::AbstractComposite)
"""
getIsDisordered(comp::AbstractComposite) = comp.is_disordered_

"""
    getResNumber(comp::AbstractComposite)
"""
getResNumber(comp::AbstractComposite) = comp.res_number_

"""
    getId(comp::AbstractComposite)
"""
getId(comp::AbstractComposite) = comp.id_

"""
    setParent(comp::AbstractComposite, value::Union{Nothing,AbstractComposite})
"""
setParent(comp::AbstractComposite, value::Union{Nothing,AbstractComposite}) = begin comp.parent_ = value end

"""
    setPrevious(comp::AbstractComposite, value::Union{Nothing,AbstractComposite})
"""
setPrevious(comp::AbstractComposite, value::Union{Nothing,AbstractComposite}) = begin comp.previous_ = value end

"""
    setNext(comp::AbstractComposite, value::Union{Nothing,AbstractComposite})
"""
setNext(comp::AbstractComposite, value::Union{Nothing,AbstractComposite}) = begin comp.next_ = value end

"""
    setFirstChild(comp::AbstractComposite, value::Union{AbstractComposite, Nothing})
"""
setFirstChild(comp::AbstractComposite, value::Union{AbstractComposite, Nothing}) = begin comp.first_child_ = value end

"""
    setLastChild(comp::AbstractComposite, value::Union{AbstractComposite, Nothing})
"""
setLastChild(comp::AbstractComposite, value::Union{AbstractComposite, Nothing}) = begin comp.last_child_ = value end

"""
    setName(comp::AbstractComposite, value::Union{String,Nothing})
"""
setName(comp::AbstractComposite, value::Union{String,Nothing}) = begin comp.name_ = value end

"""
    setTypeName(comp::AbstractComposite, value::Union{String,Nothing})
"""
setTypeName(comp::AbstractComposite, value::Union{String,Nothing}) = begin comp.type_name_ = value end

"""
    setElement(comp::AbstractComposite, value::Union{Element,Nothing})
"""
setElement(comp::AbstractComposite, value::Union{Element,Nothing}) = begin comp.element_ = value end

"""
    setRadius(comp::AbstractComposite, value::Union{Float64,Nothing})
"""
setRadius(comp::AbstractComposite, value::Union{Float64,Nothing}) = begin comp.radius_ = value end

"""
    setType(comp::AbstractComposite, value::Union{UInt8,Nothing})
"""
setType(comp::AbstractComposite, value::Union{UInt8,Nothing}) = begin comp.type_ = value end

"""
    setNumberOfBonds(comp::AbstractComposite, value::Union{UInt8,Nothing})
"""
setNumberOfBonds(comp::AbstractComposite, value::Union{UInt8,Nothing}) = begin comp.number_of_bonds_ = value end

"""
    setFormalCharge(comp::AbstractComposite, value::Union{Int64,Nothing})
"""
setFormalCharge(comp::AbstractComposite, value::Union{Int64,Nothing}) = begin comp.formal_charge_ = value end

"""
    setPosition(comp::AbstractComposite, value::Union{SVector{3,Float64},Nothing})
"""
setPosition(comp::AbstractComposite, value::Union{SVector{3,Float64},Nothing}) = begin comp.position_ = value end

"""
    setCharge(comp::AbstractComposite, value::Union{Float64,Nothing})
"""
setCharge(comp::AbstractComposite, value::Union{Float64,Nothing}) = begin comp.charge_ = value end

"""
    setVelocity(comp::AbstractComposite, value::Union{SVector{3,Float64},Nothing})
"""
setVelocity(comp::AbstractComposite, value::Union{SVector{3,Float64},Nothing}) = begin comp.velocity_ = value end

"""
    setForce(comp::AbstractComposite, value::Union{SVector{3,Float64},Nothing})
"""
setForce(comp::AbstractComposite, value::Union{SVector{3,Float64},Nothing}) = begin comp.force_ = value end

"""
    setOccupancy(comp::AbstractComposite, value::Union{Float64, Nothing})
"""
setOccupancy(comp::AbstractComposite, value::Union{Float64, Nothing}) = begin comp.occupancy_ = value end

"""
    setSerial(comp::AbstractComposite, value::Union{Int64,Nothing})
"""
setSerial(comp::AbstractComposite, value::Union{Int64,Nothing}) = begin comp.serial_ = value end

"""
    setTempFactor(comp::AbstractComposite, value::Union{Float64,Nothing})
"""
setTempFactor(comp::AbstractComposite, value::Union{Float64,Nothing}) = begin comp.temp_factor_ = value end

"""
    setBonds(comp::AbstractComposite, value::Dict{Atom, Bond})
"""
setBonds(comp::AbstractComposite, value::Dict{Atom, Bond}) = begin comp.bonds_ = value end

"""
    setProperties(comp::AbstractComposite, value::Vector{Tuple{String,UInt8}})
"""
setProperties(comp::AbstractComposite, value::Vector{Tuple{String,UInt8}}) = begin comp.properties_ = value end

"""
    setNumberOfChildren(comp::AbstractComposite, value::Union{Int64,Nothing})
"""
setNumberOfChildren(comp::AbstractComposite, value::Union{Int64,Nothing}) = begin comp.number_of_children_ = value end

"""
    setContainsSelection(comp::AbstractComposite, value::Union{Bool,Nothing})
"""
setContainsSelection(comp::AbstractComposite, value::Union{Bool,Nothing}) = begin comp.contains_selection_ = value end

"""
    setNumberOfSelectedChildren(comp::AbstractComposite, value::Union{Int64,Nothing})
"""
setNumberOfSelectedChildren(comp::AbstractComposite, value::Union{Int64,Nothing}) = begin comp.number_of_selected_children_ = value end

"""
    setNumberOfChildrenContainingSelection(comp::AbstractComposite, value::Union{Int64,Nothing})
"""
setNumberOfChildrenContainingSelection(comp::AbstractComposite, value::Union{Int64,Nothing}) = begin comp.number_of_children_containing_selection_ = value end

"""
    setSelectionStamp(comp::AbstractComposite, value::Union{TimeStamp,Nothing})
"""
setSelectionStamp(comp::AbstractComposite, value::Union{TimeStamp,Nothing}) = begin comp.selection_stamp_ = value end

"""
    setModificationStamp(comp::AbstractComposite, value::Union{TimeStamp,Nothing})
"""
setModificationStamp(comp::AbstractComposite, value::Union{TimeStamp,Nothing}) = begin comp.modification_stamp_ = value end

"""
    setTrait(comp::AbstractComposite, value::Union{AbstractComposite, Nothing})
"""
setTrait(comp::AbstractComposite, value::Union{AbstractComposite, Nothing}) = begin comp.trait_ = value end

"""
    setInsertionCode(comp::AbstractComposite, value::Union{Char,Nothing})
"""
setInsertionCode(comp::AbstractComposite, value::Union{Char,Nothing}) = begin comp.insertion_code_ = value end

"""
    setIsDisordered(comp::AbstractComposite, value::Union{Bool,Nothing})
"""
setIsDisordered(comp::AbstractComposite, value::Union{Bool,Nothing}) = begin comp.is_disordered_ = value end

"""
    setResNumber(comp::AbstractComposite, value::Union{Int64,Nothing})
"""
setResNumber(comp::AbstractComposite, value::Union{Int64,Nothing}) = begin comp.res_number_ = value end

"""
    setIsHetero(comp::AbstractComposite, value::Union{Bool,Nothing})
"""
setIsHetero(comp::AbstractComposite, value::Union{Bool,Nothing}) = begin comp.is_hetero_ = value end

"""
    setId(comp::AbstractComposite, value::Union{Char, Nothing})
"""
setId(comp::AbstractComposite, value::Union{Char, Nothing}) = begin comp.id_ = value end

#----------------DATAFRAMEROW getter setter-----------------------

"""
    getPrevious(df_row::DataFrameRow)
"""
getPrevious(df_row::DataFrameRow) = df_row.previous_

"""
    getName(df_row::DataFrameRow)
"""
getName(df_row::DataFrameRow) = df_row[1]

"""
    getTypeName(df_row::DataFrameRow)
"""
getTypeName(df_row::DataFrameRow) = df_row.type_name_

"""
    getElement(df_row::DataFrameRow)
"""
getElement(df_row::DataFrameRow) = df_row.element_

"""
    getRadius(df_row::DataFrameRow)
"""
getRadius(df_row::DataFrameRow) = df_row.radius_

"""
    getType(df_row::DataFrameRow)
"""
getType(df_row::DataFrameRow) = df_row.type_

"""
    getNumberOfBonds(df_row::DataFrameRow)
"""
getNumberOfBonds(df_row::DataFrameRow) = df_row.number_of_bonds_

"""
    getFormalCharge(df_row::DataFrameRow)
"""
getFormalCharge(df_row::DataFrameRow) = df_row.formal_charge_

"""
    getPosition(df_row::DataFrameRow)
"""
getPosition(df_row::DataFrameRow) = df_row.position_

"""
    getPosition(df::DataFrame)
"""
getPosition(df::DataFrame) = df.position_

"""
    getPosition(df::DataFrame)
"""
getPosition(df::SubDataFrame) = df.position_

"""
    getCharge(df_row::DataFrameRow)
"""
getCharge(df_row::DataFrameRow) = df_row.charge_

"""
    getVelocity(df_row::DataFrameRow)
"""
getVelocity(df_row::DataFrameRow) = df_row.velocity_

"""
    getForce(df_row::DataFrameRow)
"""
getForce(df_row::DataFrameRow) = df_row.force_

"""
    getOccupancy(df_row::DataFrameRow)
"""
getOccupancy(df_row::DataFrameRow) = df_row.occupancy_

"""
    getSerial(df_row::DataFrameRow)
"""
getSerial(df_row::DataFrameRow) = df_row.serial_

"""
    getTempFactor(df_row::DataFrameRow)
"""
getTempFactor(df_row::DataFrameRow) = df_row.temp_factor_

"""
    getProperties(df_row::DataFrameRow)
"""
getProperties(df_row::DataFrameRow) = df_row.properties_

"""
    getNumberOfChildren(df_row::DataFrameRow)
"""
getNumberOfChildren(df_row::DataFrameRow) = df_row.number_of_children_

"""
    getContainsSelection(df_row::DataFrameRow)
"""
getContainsSelection(df_row::DataFrameRow) = df_row.contains_selection_

"""
    getNumberOfSelectedChildren(df_row::DataFrameRow)
"""
getNumberOfSelectedChildren(df_row::DataFrameRow) = df_row.number_of_selected_children_

"""
    getNumberOfChildrenContainingSelection(df_row::DataFrameRow)
"""
getNumberOfChildrenContainingSelection(df_row::DataFrameRow) = df_row.number_of_children_containing_selection_

"""
    getSelectionStamp(df_row::DataFrameRow)
"""
getSelectionStamp(df_row::DataFrameRow) = df_row.selection_stamp_

"""
    getModificationStamp(df_row::DataFrameRow)
"""
getModificationStamp(df_row::DataFrameRow) = df_row.modification_stamp_

"""
    getInsertionCode(df_row::DataFrameRow)
"""
getInsertionCode(df_row::DataFrameRow) = df_row.insertion_code_

"""
    getIsDisordered(df_row::DataFrameRow)
"""
getIsDisordered(df_row::DataFrameRow) = df_row.is_disordered_

"""
    getResNumber(df_row::DataFrameRow)
"""
getResNumber(df_row::DataFrameRow) = df_row.res_number_

"""
    getIsHetero(df_row::DataFrameRow)
"""
getIsHetero(df_row::DataFrameRow) = df_row.is_hetero_

"""
    getSelected(df_row::DataFrameRow)
"""
getSelected(df_row::DataFrameRow) = df_row.selected_

"""
    getId(df_row::DataFrameRow)
"""
getId(df_row::DataFrameRow) = df_row.id_

"""
    setParent(df_row::DataFrameRow, value::Union{Nothing,DataFrameRow}
"""
setParent(df_row::DataFrameRow, value::Union{Nothing,DataFrameRow}) = begin df_row.parent_ = value end

"""
    setPrevious(df_row::DataFrameRow, value::Union{Nothing,DataFrameRow}
"""
setPrevious(df_row::DataFrameRow, value::Union{Nothing,DataFrameRow}) = begin df_row.previous_ = value end

"""
    setNext(df_row::DataFrameRow, value::Union{Nothing,DataFrameRow}
"""
setNext(df_row::DataFrameRow, value::Union{Nothing,DataFrameRow}) = begin df_row.next_ = value end

"""
    setFirstChild(df_row::DataFrameRow, value::Union{DataFrameRow, Nothing}
"""
setFirstChild(df_row::DataFrameRow, value::Union{DataFrameRow, Nothing}) = begin df_row.first_child_ = value end

"""
    setLastChild(df_row::DataFrameRow, value::Union{DataFrameRow, Nothing}
"""
setLastChild(df_row::DataFrameRow, value::Union{DataFrameRow, Nothing}) = begin df_row.last_child_ = value end

"""
    setName(df_row::DataFrameRow, value::Union{String,Nothing}
"""
setName(df_row::DataFrameRow, value::Union{String,Nothing}) = begin df_row.name_ = value end

"""
    setTypeName(df_row::DataFrameRow, value::Union{String,Nothing}
"""
setTypeName(df_row::DataFrameRow, value::Union{String,Nothing}) = begin df_row.type_name_ = value end

"""
    setElement(df_row::DataFrameRow, value::Union{Element,Nothing}
"""
setElement(df_row::DataFrameRow, value::Union{Element,Nothing}) = begin df_row.element_ = value end

"""
    setRadius(df_row::DataFrameRow, value::Union{Float64,Nothing}
"""
setRadius(df_row::DataFrameRow, value::Union{Float64,Nothing}) = begin df_row.radius_ = value end

"""
    setType(df_row::DataFrameRow, value::Union{UInt8,Nothing}
"""
setType(df_row::DataFrameRow, value::Union{UInt8,Nothing}) = begin df_row.type_ = value end

"""
    setNumberOfBonds(df_row::DataFrameRow, value::Union{UInt8,Nothing}
"""
setNumberOfBonds(df_row::DataFrameRow, value::Union{UInt8,Nothing}) = begin df_row.number_of_bonds_ = value end

"""
    setFormalCharge(df_row::DataFrameRow, value::Union{Int64,Nothing}
"""
setFormalCharge(df_row::DataFrameRow, value::Union{Int64,Nothing}) = begin df_row.formal_charge_ = value end

"""
    setPosition(df_row::DataFrameRow, value::Union{SVector{3,Float64},Nothing}
"""
setPosition(df_row::DataFrameRow, value::Union{SVector{3,Float64},Nothing}) = begin df_row.position_ = value end

"""
    setCharge(df_row::DataFrameRow, value::Union{Float64,Nothing}
"""
setCharge(df_row::DataFrameRow, value::Union{Float64,Nothing}) = begin df_row.charge_ = value end

"""
    setVelocity(df_row::DataFrameRow, value::Union{SVector{3,Float64},Nothing}
"""
setVelocity(df_row::DataFrameRow, value::Union{SVector{3,Float64},Nothing}) = begin df_row.velocity_ = value end

"""
    setForce(df_row::DataFrameRow, value::Union{SVector{3,Float64},Nothing}
"""
setForce(df_row::DataFrameRow, value::Union{SVector{3,Float64},Nothing}) = begin df_row.force_ = value end

"""
    setOccupancy(df_row::DataFrameRow, value::Union{Float64, Nothing}
"""
setOccupancy(df_row::DataFrameRow, value::Union{Float64, Nothing}) = begin df_row.occupancy_ = value end

"""
    setSerial(df_row::DataFrameRow, value::Union{Int64,Nothing}
"""
setSerial(df_row::DataFrameRow, value::Union{Int64,Nothing}) = begin df_row.serial_ = value end

"""
    setTempFactor(df_row::DataFrameRow, value::Union{Float64,Nothing}
"""
setTempFactor(df_row::DataFrameRow, value::Union{Float64,Nothing}) = begin df_row.temp_factor_ = value end

"""
    setBonds(df_row::DataFrameRow, value::Dict{Atom, Bond}
"""
setBonds(df_row::DataFrameRow, value::Dict{Atom, Bond}) = begin df_row.bonds_ = value end

"""
    setProperties(df_row::DataFrameRow, value::Vector{Tuple{String,UInt8}}
"""
setProperties(df_row::DataFrameRow, value::Vector{Tuple{String,UInt8}}) = begin df_row.properties_ = value end

"""
    setNumberOfChildren(df_row::DataFrameRow, value::Union{Int64,Nothing}
"""
setNumberOfChildren(df_row::DataFrameRow, value::Union{Int64,Nothing}) = begin df_row.number_of_children_ = value end

"""
    setContainsSelection(df_row::DataFrameRow, value::Union{Bool,Nothing}
"""
setContainsSelection(df_row::DataFrameRow, value::Union{Bool,Nothing}) = begin df_row.contains_selection_ = value end

"""
    setNumberOfSelectedChildren(df_row::DataFrameRow, value::Union{Int64,Nothing}
"""
setNumberOfSelectedChildren(df_row::DataFrameRow, value::Union{Int64,Nothing}) = begin df_row.number_of_selected_children_ = value end

"""
    setNumberOfChildrenContainingSelection(df_row::DataFrameRow, value::Union{Int64,Nothing}
"""
setNumberOfChildrenContainingSelection(df_row::DataFrameRow, value::Union{Int64,Nothing}) = begin df_row.number_of_children_containing_selection_ = value end

"""
    setSelectionStamp(df_row::DataFrameRow, value::Union{TimeStamp,Nothing}
"""
setSelectionStamp(df_row::DataFrameRow, value::Union{TimeStamp,Nothing}) = begin df_row.selection_stamp_ = value end

"""
    setModificationStamp(df_row::DataFrameRow, value::Union{TimeStamp,Nothing}
"""
setModificationStamp(df_row::DataFrameRow, value::Union{TimeStamp,Nothing}) = begin df_row.modification_stamp_ = value end

"""
    setTrait(df_row::DataFrameRow, value::Union{DataFrameRow, Nothing}
"""
setTrait(df_row::DataFrameRow, value::Union{DataFrameRow, Nothing}) = begin df_row.trait_ = value end

"""
    setInsertionCode(df_row::DataFrameRow, value::Union{Char,Nothing}
"""
setInsertionCode(df_row::DataFrameRow, value::Union{Char,Nothing}) = begin df_row.insertion_code_ = value end

"""
    setIsDisordered(df_row::DataFrameRow, value::Union{Bool,Nothing}
"""
setIsDisordered(df_row::DataFrameRow, value::Union{Bool,Nothing}) = begin df_row.is_disordered_ = value end

"""
    setResNumber(df_row::DataFrameRow, value::Union{Int64,Nothing}
"""
setResNumber(df_row::DataFrameRow, value::Union{Int64,Nothing}) = begin df_row.res_number_ = value end

"""
    setIsHetero(df_row::DataFrameRow, value::Union{Bool,Nothing}
"""
setIsHetero(df_row::DataFrameRow, value::Union{Bool,Nothing}) = begin df_row.is_hetero_ = value end

"""
    setSelected(df_row::DataFrameRow, value::Bool
"""
setSelected(df_row::DataFrameRow, value::Bool) = begin df_row.selected_ = value end

"""
    setId(df_row::DataFrameRow, value::Union{Char, Nothing}
"""
setId(df_row::DataFrameRow, value::Union{Char, Nothing}) = begin df_row.id_ = value end


text:
using Distributions
using GaussianProcesses
using Plots, StatPlots
using LaTeXStrings
#  include("X3_empirical.jl")
include("utilities.jl")
include("self_confidence.jl")

function annotate_pct(coords::Tuple{Float64,Float64},txt::Plots.PlotText,plt::Plots.Plot)
    xmin = plt.subplots[1].attr[:xaxis][:extrema].emin
    xmax = plt.subplots[1].attr[:xaxis][:extrema].emax
    ymin = plt.subplots[1].attr[:yaxis][:extrema].emin
    ymax = plt.subplots[1].attr[:yaxis][:extrema].emax
    #  println("$xmin, $xmax, $ymin, $ymax")
    Plots.annotate!(coords[1]*(xmax-xmin)+xmin,coords[2]*(ymax-ymin)+ymin,txt)
end

srand(12345)
trusted_n=5;                          #number of training points
trusted_x = squeeze([0. 1. rand(1,trusted_n-2)],1);              #predictors
#  trusted_y = ones(trusted_n);    #regressors
trusted_y = zeros(trusted_n);    #regressors

candidate_n=5;
candidate_x = squeeze([0. 1. 0.1464466 rand(1,candidate_n-3)],1);              #predictors
#  candidate_y = (candidate_x-0.5).^2 - 0.125 + 1.;    #regressors
candidate_y = (candidate_x-0.5).^2 - 0.125;    #regressors

#Select mean and covariance function
#  mZero = MeanConst(1.0)
mZero = MeanZero()
kern = SE(-0.5,0.0)                   #Sqaured exponential kernel (note that hyperparameters are on the log scale)

logObsNoise_trust = -4.0                        # log standard deviation of observation noise (this is optional)
logObsNoise_cand = -4.0                        # log standard deviation of observation noise (this is optional)
trusted_gp = GP(trusted_x,trusted_y,mZero,kern,logObsNoise_trust)       #Fit the GP
candidate_gp = GP(candidate_x,candidate_y,mZero,kern,logObsNoise_cand)       #Fit the GP

#  points_of_interest_x = [0.150,0.5,0.84,0.985]
points_of_interest_x = [0.15,0.5,0.84,0.985]
#  points_of_interest_y = [1.,0.89,1.04,1.116]
points_of_interest_y = [0.,-0.125,0.03,0.116]

μ_trust, σ²_trust = predict_y(trusted_gp,points_of_interest_x);
μ_cand, σ²_cand = predict_y(candidate_gp,points_of_interest_x);

dist_trust = []
for (m,s) in zip(μ_trust,sqrt.(σ²_trust))
    push!(dist_trust,Normal(m,s))
end
dist_cand = []
for (m,s) in zip(μ_cand,sqrt.(σ²_cand))
    push!(dist_cand,Normal(m,s))
end

SQ = []
SQ_raw = []
solver_reward_range = [[-0.,5.], [-0.,0.5], [-0.,0.05], [-0.,0.005]]
for (c,t) in zip(dist_cand,dist_trust)

    sqs = []
    sq_raw = []
    for i = 1:length(solver_reward_range)
        println("#################################")
        println("i: $i")
        (sq,sqr) = X3(c,t,global_rwd_range=solver_reward_range[i],return_raw_sq=true)
        println("using $(solver_reward_range[i]) reward range, got $sq")
        push!(sqs,sq)
        push!(sq_raw,sqr)
    end
    push!(SQ,sqs)
    push!(SQ_raw,sq_raw)
end

show_obs = false
beta_val = 0.999

Plots.scalefontsizes(1.0)

#  p1 = Plots.plot(trusted_gp,β=beta_val,linecolor=:blue,fillcolor=:blue,lw=3,obsv=show_obs,fillalpha=0.2,label=L"test",legend=:topleft)
#  Plots.plot!(candidate_gp,β=beta_val,linecolor=:red,fillcolor=:red,lw=2,obsv=show_obs,fillalpha=0.2,label=L"test2")
p1 = Plots.plot(trusted_gp,β=beta_val,linecolor=:blue,fillcolor=:blue,lw=3,obsv=show_obs,fillalpha=0.2,label=L"trusted $\mathcal{S}^T$",legend=:topleft)
Plots.plot!(candidate_gp,β=beta_val,linecolor=:red,fillcolor=:red,lw=2,obsv=show_obs,fillalpha=0.2,label=L"candidate $\mathcal{S}^C$")
Plots.scatter!(points_of_interest_x,points_of_interest_y,marker=:star,markersize=10,label="points of interest")
xlabel!("task parameter")
ylabel!("reward")
for i = 1:length(points_of_interest_x)
    point_label = '@' + i
    Plots.annotate!(points_of_interest_x[i],points_of_interest_y[i],Plots.text("$point_label",:black,:right,:bottom))
end

p2 = Plots.plot(dist_trust[1],linecolor=:blue,fill=(0,0.2,:blue),legend=false,title="A")
Plots.plot!(dist_cand[1],linecolor=:red,fill=(0,0.2,:red))
annotate_pct((0.65,0.8),Plots.text(LaTeXString(@sprintf("\$x_Q(r=5):%0.3f\$\n\$x_Q(r=0.5):%0.3f\$\n\$x_Q(r=0.05):%0.3f\$\n\$x_Q(r=0.005)%0.3f\$",SQ[1][1],SQ[1][2],SQ[1][3],SQ[1][4])),14,:center,:left),p2)
#  annotate_pct((0.65,0.8),Plots.text(LaTeXString(@sprintf("\$x^2\$")),14,:center,:left),p2)
xlabel!("reward")
ylabel!("p(rwd)")
#  annotate_pct((0.8,0.8),text("SQ: $(@sprintf("%0.3f",SQ[1][1]))",:center,:center),p2)
#  annotate!(0.07,10.,text(L"SQ\rightarrow High?",:black,:center))

p3 = Plots.plot(dist_trust[2],linecolor=:blue,fill=(0,0.2,:blue),legend=false,title="B")
Plots.plot!(dist_cand[2],linecolor=:red,fill=(0,0.2,:red))
annotate_pct((0.15,0.8),Plots.text(LaTeXString(@sprintf("\$x_Q(r=5):%0.3f\$\n\$x_Q(r=0.5):%0.3f\$\n\$x_Q(r=0.05):%0.3f\$\n\$x_Q(r=0.005)%0.3f\$",SQ[2][1],SQ[2][2],SQ[2][3],SQ[2][4])),14,:center,:left),p3)
xlabel!("reward")
ylabel!("p(rwd)")
#  annotate_pct((0.8,0.8),text("SQ: $(@sprintf("%0.3f",SQ[2][1]))",:center,:center),p3)
#  annotate!(-0.15,12.,text(L"SQ\rightarrow Very Low",:black,:center))

p4 = Plots.plot(dist_trust[3],linecolor=:blue,fill=(0,0.2,:blue),legend=false,title="C")
Plots.plot!(dist_cand[3],linecolor=:red,fill=(0,0.2,:red))
annotate_pct((0.65,0.8),Plots.text(LaTeXString(@sprintf("\$x_Q(r=5):%0.3f\$\n\$x_Q(r=0.5):%0.3f\$\n\$x_Q(r=0.05):%0.3f\$\n\$x_Q(r=0.005)%0.3f\$",SQ[3][1],SQ[3][2],SQ[3][3],SQ[3][4])),14,:center,:left),p4)
xlabel!("reward")
ylabel!("p(rwd)")
#  annotate_pct((0.8,0.8),text("SQ: $(@sprintf("%0.3f",SQ[3][1]))",:center,:center),p4)
#  annotate!(0.15,12.,text(L"SQ\rightarrow Low?",:black,:center))

p5 = Plots.plot(dist_trust[4],linecolor=:blue,fill=(0,0.2,:blue),legend=false,title="D")
Plots.plot!(dist_cand[4],linecolor=:red,fill=(0,0.2,:red))
annotate_pct((0.65,0.8),Plots.text(LaTeXString(@sprintf("\$x_Q(r=5):%0.3f\$\n\$x_Q(r=0.5):%0.3f\$\n\$x_Q(r=0.05):%0.3f\$\n\$x_Q(r=0.005)%0.3f\$",SQ[4][1],SQ[4][2],SQ[4][3],SQ[4][4])),14,:center,:left),p5)
#  annotate_pct((0.65,0.8),Plots.text(@sprintf("SQ(r=5):%0.3f\nSQ(r=0.5):%0.3f\nSQ(r=0.05):%0.3f\nSQ(r=0.005)%0.3f",SQ[4][1],SQ[4][2],SQ[4][3],SQ[4][4]),14,:center,:left),p5)
xlabel!("reward")
ylabel!("p(rwd)")
#  annotate_pct((0.8,0.8),text("SQ: $(@sprintf("%0.3f",SQ[4][1]))",:center,:center),p5)
#  annotate!(0.15,12.,text(L"SQ\rightarrow Very High",:black,:center))

# empty filler frames
pe1 = Plots.plot(framestyle = :none)
annotate!(0.1,0.5,text("Solver Quality:\n  0->Poor Quality\n  1-> Equal to Trusted\n  2->Better than Trusted\nThink of as fraction of quality\nthresholded at 2",10,:center,:left))
#  annotate!(0.,0.,text("empty",:black,:center))
pe2 = Plots.plot(framestyle = :none)
Plots.annotate!(0.1,0.65,Plots.text("* Calculate mean difference weighted by std",10,:center,:left))
Plots.annotate!(0.5,0.5,Plots.text(L"wd=\frac{m(c)/s(c)-m(t)/s(t)}{range(rwd)}",10,:center,:center))
Plots.annotate!(0.1,0.3,Plots.text("* SQ is a scaled Sigmoid on wd",10,:center,:left))
Plots.annotate!(0.5,0.15,Plots.text(L"SQ=\frac{2exp(wd)}{(exp(wd)+1)}",10,:center,:center))
#  annotate!(0.,0.,text("empty",:black,:center))
pe3 = Plots.plot(framestyle = :none)
Plots.annotate!(0.1,0.65,Plots.text("* In these plots:",10,:center,:left))
Plots.annotate!(0.5,0.5,Plots.text(LaTeXString(@sprintf("range(rwd) = %s",solver_reward_range)),10,:center,:center))
#  annotate!(0.,0.,text("empty",:black,:center))

p1notes = Plots.plot(framestyle = :none)
Plots.annotate!(0.5,0.5,Plots.text(@sprintf("SQ: %0.3f\n%0.3f\n%0.3f\n%0.3f",SQ[1][1],SQ[1][2],SQ[1][3],SQ[1][4]),10,:center,:center))

p2notes = Plots.plot(framestyle = :none)
Plots.annotate!(0.5,0.5,Plots.text("SQ: $(@sprintf("%0.3f",SQ[2][1]))",10,:center,:center))

p3notes = Plots.plot(framestyle = :none)
Plots.annotate!(0.5,0.5,Plots.text("SQ: $(@sprintf("%0.3f",SQ[3][1]))",10,:center,:center))

p4notes = Plots.plot(framestyle = :none)
Plots.annotate!(0.5,0.5,Plots.text("SQ: $(@sprintf("%0.3f",SQ[4][1]))",10,:center,:center))

savefig(p1,"figs/p1.pdf")
savefig(p2,"figs/p2.pdf")
savefig(p3,"figs/p3.pdf")
savefig(p4,"figs/p4.pdf")
savefig(p5,"figs/p5.pdf")
paper_layout = Plots.plot(p2,p3,p4,p5,size=(1000,800))
savefig(paper_layout,"figs/point_compare.pdf")
Plots.plot(p1,p2,pe1,p3,pe2,p4,pe3,p5,layout=(4,2),size=(800,800))

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
using SimTools
using Test

@testset "SimTools.jl" begin
    include("bisection_search_tests.jl")
    include("grad_desc_tests.jl")
    include("OnDiskMonitor_tests.jl")
end