text:
using IterableTables
using Distributions
using Statistics

"""
    RandomVariable{T}

A RandomVariable can be "assigned" to different model parameters. Its
value can be used directly or applied by adding to, or multiplying by,
reference values for the indicated parameter. Note that the distribution
must be instantiated, e.g., call as RandomVariable(:foo, Normal()), not
RandomVariable(:foo, Normal).
"""
struct RandomVariable{T}
    name::Symbol
    dist::T

    function RandomVariable(name::Symbol, dist::T) where T
        self = new{T}(name, dist)
        return self
    end
end

Base.eltype(rv::RandomVariable) = eltype(rv.dist)

distribution(rv::RandomVariable) = rv.dist

@enum ScenarioLoopPlacement OUTER INNER
@enum SamplingOptions LHS RANDOM

# SampleStore is a faux Distribution that implements base.rand()
# to yield stored values.
mutable struct SampleStore{T}
    values::Vector{T}   # generally Int or Float64
    idx::Int            # index of next value to return

    function SampleStore(values::Vector{T}) where T
        return new{T}(values, 1)
    end
end

Base.eltype(ss::SampleStore) = eltype(ss.values)

distribution(ss::SampleStore) = ss.dist

#
# TBD: maybe have a different SampleStore subtype for values drawn from a dist
# versus those loaded from a file, which would be treated as immutable?
#

function Statistics.quantile(ss::SampleStore{T}, q::Float64) where T
    return quantile(sort(ss.values), q)
end

# TBD: This interpolates values between those in the vector. Is this reasonable?
# Probably shouldn't use correlation on values loaded from a file rather than
# from a proper distribution.
function Statistics.quantile(ss::SampleStore{T}, probs::AbstractArray) where T
    return quantile.(ss, probs)
end

struct TransformSpec
    paramname::Symbol
    op::Symbol
    rvname::Symbol
    dims::Vector{Any}

    function TransformSpec(paramname::Symbol, op::Symbol, rvname::Symbol, dims::Vector{T}=[]) where T
        if ! (op in (:(=), :(+=), :(*=)))
            error("Valid operators are =, +=, and *= (got $op)")
        end

        return new(paramname, op, rvname, dims)
    end
end

"""
    CorrelationSpec

Defines a target rank correlation to establish between the two named random vars.
"""
struct CorrelationSpec
    name1::Symbol
    name2::Symbol
    value::Float64
end

"""
    Base.rand(s::SampleStore{T}, n::Int=1) where T

Pseudo `rand`` function that just returns the next `n` values from the vector of
stored values. Currently does not support cycling through results multiple times,
but you can `reset` the `SampleStore` to reuse it.
"""
function Base.rand(s::SampleStore{T}, n::Int=1) where T
    idx = s.idx
    s.idx = next_idx = idx + n
    return n > 1 ? s.values[idx:(next_idx - 1)] : s.values[idx]
end

function Base.reset(s::SampleStore{T}) where T
    s.idx = 1
    return nothing
end

# debugging tool
# global save_dict = Dict()

"""
    MonteCarloSimulation

Holds all the data that defines a Monte Carlo simulation.
"""
mutable struct MonteCarloSimulation
    trials::Int
    current_trial::Int
    current_data::Any               # holds data for current_trial when current_trial > 0
    rvdict::OrderedDict{Symbol, RandomVariable}
    translist::Vector{TransformSpec}
    corrlist::Vector{CorrelationSpec}
    savelist::Vector{Tuple{Symbol, Symbol}}
    dist_rvs::Vector{RandomVariable}
    nt_type::Any                    # a generated NamedTuple type to hold data for a single trial
    models::Vector{Model}
    results::Vector{Dict{Tuple, DataFrame}}

    function MonteCarloSimulation(rvlist::Vector,
                                  translist::Vector{TransformSpec},
                                  corrlist::Vector{CorrelationSpec},
                                  savelist::Vector{Tuple{Symbol, Symbol}})
        self = new()
        self.trials = 0
        self.current_trial = 0
        self.current_data = nothing
        self.rvdict = OrderedDict([rv.name => rv for rv in rvlist])
        self.translist = translist
        self.corrlist = corrlist
        self.savelist = savelist
        self.dist_rvs = [rv for rv in rvlist]

        names = (keys(self.rvdict)...,)
        types = [eltype(fld) for fld in values(self.rvdict)]
        self.nt_type = NamedTuple{names, Tuple{types...}}

        # These are parallel arrays; each model has a corresponding results dict
        self.models = Vector{Model}(undef, 0)
        self.results = [Dict{Tuple, DataFrame}()]
        # save_dict[:mcs] = self
        return self
    end
end

struct MCSIterator{NT}
    mcs::MonteCarloSimulation

    function MCSIterator{NT}(mcs::MonteCarloSimulation) where NT
        return new{NT}(mcs)
    end
end


text:
"Abstract type for Thunks."
abstract type Think end
abstract type WrappedThink <: Think end
"""
    Thunk(function, args, kwargs)

Type that represents a thunk. Only `evaluated` is a public field.
Use `reify` to evaluate. Evaluates to Weak head normal form, meaning that
non-free symbols must be assigned / defined. Contrast with Unevaluated.
"""
mutable struct Thunk <: Think
    f::Any # usually a Function, but could be any callable
    # args will be passed to f. Needs to support ... (splat), eg Array or Tuple
    args::Any
    kwargs::Dict # kwargs will be passed to f, cleared after evaluation
    evaluated::Bool # false until computed, then true
    result::Any # cache result once computed
    Thunk(f, args) = new(f, args, Dict(), false, nothing)
    Thunk(f, args, kwargs) = new(f, args, kwargs, false, nothing)
end

"""
Keep expression, including args/kwargs, as unevaluated.
"""
mutable struct Unevaluated <: Think
    f::Any # usually a Function, but could be any callable
    # args will be passed to f. Needs to support ... (splat), eg Array or Tuple
    args::Any # a primitive thunk of the form: () -> tuple
    kwargs::Any # a primitive thunk of the form: () -> dict
    evaluated::Bool # false until computed, then true
    result::Any # cache result once computed
    Unevaluated(f, args) = new(f, args, ()->Dict(), false, nothing)
    Unevaluated(f, args, kwargs) = new(f, args, kwargs, false, nothing)
end

"""
A Thunk that can "undo" code evaluation. Does not clear f or args.

Useful for interactive coding. ie if you revise a function definition,
can undo & reify the thunk again.
"""
mutable struct Reversible <: Think
    f::Any # usually a Function, but could be any callable
    # args will be passed to f. Needs to support ... (splat), eg Array or Tuple
    args::Any
    kwargs::Dict # kwargs will be passed to f, cleared after evaluation
    evaluated::Bool # false until computed, then true
    result::Any # cache result once computed
    Reversible(f, args) = new(f, args, Dict(), false, nothing)
    Reversible(f, args, kwargs) = new(f, args, kwargs, false, nothing)
end

"No operation."
function noop(args...; kwargs...)
    nothing
end

"""
A Thunk that can be checkpointed.
"""
mutable struct Checkpointable <: WrappedThink
    wrapped_thunk::Think
    # function that takes value of wrapped_thunk, and stores the result, ie on disk
    checkpoint::Any
    # instead of evaluating wrapped_thunk, try restoring value from disk
    restore::Any # load the result. ideally, no dependencies
    Checkpointable(t,c,r) = new(t,c,r)
    Checkpointable(t,r) = new(t,noop,r)
end

function Base.getindex(self::Think, index)
    thunk(getindex)(self, index)
end

Base.iterate(t::Think, state=1) = (t[state], state+1)

function thunk(f)
    (args...; kwargs...) -> Thunk(f, args, kwargs)
end

"""
    reify(thunk::Think)
    reify(value::Any)

Reify a thunk into a value.

In other words, compute the value of the expression.

We walk through the thunk's arguments and keywords, recursively evaluating each one,
and then evaluating the thunk's function with the evaluated arguments.
"""

function set_result(thunk::WrappedThink)
    set_result(thunk.wrapped_thunk)
end

function set_result(thunk::Think, result)
    thunk.result = result
    thunk.evaluated = true
    # clear to allow garbage collection
    thunk.args = []
    thunk.kwargs = Dict()
end

function set_result(thunk::Reversible, result)
    thunk.result = result
    thunk.evaluated = true
end

function reify(thunk::Think)
    if thunk.evaluated
        return thunk.result
    else
        args = [reify(x) for x in thunk.args]
        kwargs = Dict(k => reify(v) for (k,v) in thunk.kwargs)
        result = thunk.f(args...; kwargs...)
        set_result(thunk, result)
        return result
    end
end

function reify(thunk::Unevaluated)
    if thunk.evaluated
        return thunk.result
    else
        args = [reify(x) for x in thunk.args()]
        kwargs = Dict(k => reify(v) for (k,v) in thunk.kwargs())
        result = thunk.f(args...; kwargs...)
        set_result(thunk, result)
        return result
    end
end

function reify(thunk::Reversible)
    if thunk.evaluated
        return thunk.result
    else
        args = [reify(x) for x in thunk.args]
        kwargs = Dict(k => reify(v) for (k,v) in thunk.kwargs)
        result = thunk.f(args...; kwargs...)
        set_result(thunk, result)
        return result
    end
end

function reify(thunk::Checkpointable)
    if thunk.wrapped_thunk.evaluated
        return thunk.wrapped_thunk.result
    end
    result = nothing
    try
        result = reify(thunk.restore)
    catch
        for (exc, bt) in Base.catch_stack()
            showerror(stdout, exc, bt)
            println(stdout)
        end
    end
    if ~isnothing(result)
        set_result(thunk.wrapped_thunk, result)
        return result
    else
        result = reify(thunk.wrapped_thunk)
        reify(thunk.checkpoint)(result)
        return result
    end
end

function reify(value)
    value
end

"Undo function call (non-recursively)"
function undo(thunk::Reversible)
    thunk.result = nothing
    thunk.evaluated = false
    thunk
end

"Undo function call (non-recursively)"
function undo(thunk::WrappedThink)
    undo(thunk.wrapped_thunk)
    thunk
end

ext: jl


text:
using FractionalDiffEq
function EnzymeKinetics!(dy, y, ϕ, t)
    dy[1] = 10.5-y[1]/(1+0.0005*ϕ[4]^3)
    dy[2] = y[1]/(1+0.0005*ϕ[4]^3)-y[2]
    dy[3] = y[2]-y[3]
    dy[4] = y[3]-0.5*y[4]
end
q = [60, 10, 10, 20]; α = [0.95, 0.95, 0.95, 0.95]
prob = FDDESystem(EnzymeKinetics!, q, α, 4, 6)
sold, sol = solve(prob, 0.1, DelayABM())
tspan = collect(0:0.1:2)
using Plots
plot(tspan, sol[:, 1])
plot!(tspan, sol[:, 2])
plot!(tspan, sol[:, 3])
plot!(tspan, sol[:, 4])


text:
@testset "Misc" begin
  include("Misc/psibound.jl")
  include("Misc/FiniteField.jl")
  include("Misc/Primes.jl")
  include("Misc/Places.jl")
  include("Misc/meataxe.jl")
  include("Misc/stable_subgroups.jl")
  include("Misc/NumberField.jl")
  include("Misc/Poly.jl")
  include("Misc/Localization.jl")
  include("Misc/OrdLocalization.jl")
  include("Misc/jordan_test.jl")
  include("Misc/RelFinField.jl")
end


text:
# how struct padding affects memory usage
using DataFrames

# define a struct that contains two fields
struct TwoFields{S,T}
    x::S
    y::T
end

# create arrays of TwoFields having different combinations of Int types
N = 100
overhead = 40  # number of bytes for managing the array
results = DataFrame(xtype=Type[], ytype=Type[], n=Int[], size=Int[], unpadded_size=Int[])
for S ∈ (Int8, Int16, Int32, Int64), T ∈ (Int8, Int16, Int32, Int64)
    var = Symbol("two_$(S)_$(T)")
    obj = [TwoFields{S,T}(S(i), T(i)) for i in 1:N]
    push!(results, (xtype=S, ytype=T, n=N,
                    size=Base.summarysize(obj),
                    unpadded_size=overhead + N * (sizeof(S) + sizeof(T))))
end
results.diff = results.size .- results.unpadded_size
results[Symbol("diff%")] = round.(Int, results.diff ./ results.size .* 100)
results

# results
#=
julia> results
16×7 DataFrame
│ Row │ xtype │ ytype │ n     │ size  │ unpadded_size │ diff  │ diff% │
│     │ Type  │ Type  │ Int64 │ Int64 │ Int64         │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┼───────┼───────────────┼───────┼───────┤
│ 1   │ Int8  │ Int8  │ 100   │ 240   │ 240           │ 0     │ 0     │
│ 2   │ Int8  │ Int16 │ 100   │ 440   │ 340           │ 100   │ 23    │
│ 3   │ Int8  │ Int32 │ 100   │ 840   │ 540           │ 300   │ 36    │
│ 4   │ Int8  │ Int64 │ 100   │ 1640  │ 940           │ 700   │ 43    │
│ 5   │ Int16 │ Int8  │ 100   │ 440   │ 340           │ 100   │ 23    │
│ 6   │ Int16 │ Int16 │ 100   │ 440   │ 440           │ 0     │ 0     │
│ 7   │ Int16 │ Int32 │ 100   │ 840   │ 640           │ 200   │ 24    │
│ 8   │ Int16 │ Int64 │ 100   │ 1640  │ 1040          │ 600   │ 37    │
│ 9   │ Int32 │ Int8  │ 100   │ 840   │ 540           │ 300   │ 36    │
│ 10  │ Int32 │ Int16 │ 100   │ 840   │ 640           │ 200   │ 24    │
│ 11  │ Int32 │ Int32 │ 100   │ 840   │ 840           │ 0     │ 0     │
│ 12  │ Int32 │ Int64 │ 100   │ 1640  │ 1240          │ 400   │ 24    │
│ 13  │ Int64 │ Int8  │ 100   │ 1640  │ 940           │ 700   │ 43    │
│ 14  │ Int64 │ Int16 │ 100   │ 1640  │ 1040          │ 600   │ 37    │
│ 15  │ Int64 │ Int32 │ 100   │ 1640  │ 1240          │ 400   │ 24    │
│ 16  │ Int64 │ Int64 │ 100   │ 1640  │ 1640          │ 0     │ 0     │
=#


text:
module Auditor

using BinaryBuilderBase
using Base.BinaryPlatforms
using Pkg
using ObjectFile

using BinaryBuilderBase: march

export audit, collect_files, collapse_symlinks

include("auditor/instruction_set.jl")
include("auditor/dynamic_linkage.jl")
include("auditor/symlink_translator.jl")
include("auditor/compiler_abi.jl")
include("auditor/soname_matching.jl")
include("auditor/filesystems.jl")
include("auditor/extra_checks.jl")
include("auditor/codesigning.jl")

# AUDITOR TODO LIST:
#
# * Build dlopen() clone that inspects and tries to figure out why
#   something can't be opened.  Possibly use that within BinaryProvider too?

"""
    audit(prefix::Prefix, src_name::AbstractString = "";
                          io=stderr,
                          platform::AbstractPlatform = HostPlatform(),
                          verbose::Bool = false,
                          silent::Bool = false,
                          autofix::Bool = false,
                          has_csl::Bool = true,
                          require_license::Bool = true,
          )

Audits a prefix to attempt to find deployability issues with the binary objects
that have been installed within.  This auditing will check for relocatability
issues such as dependencies on libraries outside of the current `prefix`,
usage of advanced instruction sets such as AVX2 that may not be usable on many
platforms, linkage against newer glibc symbols, etc...

This method is still a work in progress, only some of the above list is
actually implemented, be sure to actually inspect `Auditor.jl` to see what is
and is not currently in the realm of fantasy.
"""
function audit(prefix::Prefix, src_name::AbstractString = "";
                               io=stderr,
                               platform::AbstractPlatform = HostPlatform(),
                               verbose::Bool = false,
                               silent::Bool = false,
                               autofix::Bool = false,
                               has_csl::Bool = true,
                               require_license::Bool = true,
               )
    # This would be really weird, but don't let someone set `silent` and `verbose` to true
    if silent
        verbose = false
    end

    # Canonicalize immediately
    prefix = Prefix(realpath(prefix.path))
    if verbose
        @info("Beginning audit of $(prefix.path)")
    end

    # If this is false then it's bedtime for bonzo boy
    all_ok = true

    # Translate absolute symlinks to relative symlinks, if possible
    translate_symlinks(prefix.path; verbose=verbose)

    # Inspect binary files, looking for improper linkage
    predicate = f -> (filemode(f) & 0o111) != 0 || valid_library_path(f, platform)
    bin_files = collect_files(prefix, predicate; exclude_externalities=false)
    for f in collapse_symlinks(bin_files)
        # If `f` is outside of our prefix, ignore it.  This happens with files from our dependencies
        if !startswith(f, prefix.path)
            continue
        end

        # Peel this binary file open like a delicious tangerine
        try
            readmeta(f) do oh
                if !is_for_platform(oh, platform)
                    if verbose
                        @warn("Skipping binary analysis of $(relpath(f, prefix.path)) (incorrect platform)")
                    end
                else
                    # Check that the ISA isn't too high
                    all_ok &= check_isa(oh, platform, prefix; verbose, silent)
                    # Check that the OS ABI is set correctly (often indicates the wrong linker was used)
                    all_ok &= check_os_abi(oh, platform; verbose)
                    # Make sure all binary files are executables, if libraries aren't
                    # executables Julia may not be able to dlopen them:
                    # https://github.com/JuliaLang/julia/issues/38993.  In principle this
                    # should be done when autofix=true, but we have to run this fix on MKL
                    # for Windows, for which however we have to set autofix=false:
                    # https://github.com/JuliaPackaging/Yggdrasil/pull/922.
                    all_ok &= ensure_executability(oh; verbose, silent)

                    # If this is a dynamic object, do the dynamic checks
                    if isdynamic(oh)
                        # Check that the libgfortran version matches
                        all_ok &= check_libgfortran_version(oh, platform; verbose, has_csl)
                        # Check whether the library depends on any of the most common
                        # libraries provided by `CompilerSupportLibraries_jll`.
                        all_ok &= check_csl_libs(oh, platform; verbose, has_csl)
                        # Check that the libstdcxx string ABI matches
                        all_ok &= check_cxxstring_abi(oh, platform; verbose)
                        # Check that this binary file's dynamic linkage works properly.  Note to always
                        # DO THIS ONE LAST as it can actually mutate the file, which causes the previous
                        # checks to freak out a little bit.
                        all_ok &= check_dynamic_linkage(oh, prefix, bin_files;
                                                        platform, silent, verbose, autofix, src_name)
                    end
                end
            end

            # Ensure this file is codesigned (currently only does something on Apple platforms)
            all_ok &= ensure_codesigned(f, prefix, platform; verbose, subdir=src_name)
        catch e
            if !isa(e, ObjectFile.MagicMismatch)
                rethrow(e)
            end

            # If this isn't an actual binary file, skip it
            if verbose
                @info("Skipping binary analysis of $(relpath(f, prefix.path))")
            end
        end
    end

    # Find all dynamic libraries
    shlib_files = filter(f -> startswith(f, prefix.path) && valid_library_path(f, platform), collapse_symlinks(bin_files))

    for f in shlib_files
        # Inspect all shared library files for our platform (but only if we're
        # running native, don't try to load library files from other platforms)
        if platforms_match(platform, HostPlatform())
            if verbose
                @info("Checking shared library $(relpath(f, prefix.path))")
            end

            # dlopen() this library in a separate Julia process so that if we
            # try to do something silly like `dlopen()` a .so file that uses
            # LLVM in interesting ways on startup, it doesn't kill our main
            # Julia process.
            dlopen_cmd = """
                using Libdl
                try
                    dlopen($(repr(f)))
                    exit(0)
                catch e
                    if $(repr(verbose))
                        Base.display_error(e)
                    end
                    exit(1)
                end
            """
            try
                p = open(`$(Base.julia_cmd()) -e $dlopen_cmd`)
                wait(p)
                if p.exitcode != 0
                    throw("Invalid exit code!")
                end
            catch
                # TODO: Use the relevant ObjFileBase packages to inspect why
                # this file is being nasty to us.
                if !silent
                    @warn("$(relpath(f, prefix.path)) cannot be dlopen()'ed")
                end
                all_ok = false
            end
        end

        # Ensure that all libraries have at least some kind of SONAME, if we're
        # on that kind of platform
        if !Sys.iswindows(platform)
            all_ok &= ensure_soname(prefix, f, platform; verbose, autofix, subdir=src_name)
        end

        # Ensure that this library is available at its own SONAME
        all_ok &= symlink_soname_lib(f; verbose=verbose, autofix=autofix)
    end

    if Sys.iswindows(platform)
        # We also cannot allow any symlinks in Windows because it requires
        # Admin privileges to create them.  Orz
        symlinks = collect_files(prefix, islink, exclude_dirs = false)
        for f in symlinks
            try
                src_path = realpath(f)
                if isfile(src_path) || isdir(src_path)
                    rm(f; force=true)
                    cp(src_path, f, follow_symlinks = true)
                end
            catch
            end
        end

        # If we're targeting a windows platform, check to make sure no .dll
        # files are sitting in `$prefix/lib`, as that's a no-no.  This is
        # not a fatal offense, but we'll yell about it.
        lib_dll_files =  filter(f -> valid_library_path(f, platform), collect_files(joinpath(prefix, "lib"), predicate))
        for f in lib_dll_files
            if !silent
                @warn("$(relpath(f, prefix.path)) should be in `bin`!")
            end
        end

        # Even more than yell about it, we're going to automatically move
        # them if there are no `.dll` files outside of `lib`.  This is
        # indicative of a simplistic build system that just don't know any
        # better with regards to windows, rather than a complicated beast.
        outside_dll_files = [f for f in shlib_files if !(f in lib_dll_files)]
        if autofix && !isempty(lib_dll_files) && isempty(outside_dll_files)
            if !silent
                @warn("Simple buildsystem detected; Moving all `.dll` files to `bin`!")
            end

            mkpath(joinpath(prefix, "bin"))
            for f in lib_dll_files
                mv(f, joinpath(prefix, "bin", basename(f)))
            end
        end
    end

    # Check that we're providing a license file
    if require_license
        all_ok &= check_license(prefix, src_name; verbose=verbose, silent=silent)
    end

    # Perform filesystem-related audit passes
    predicate = f -> !startswith(f, joinpath(prefix, "logs"))
    all_files = collect_files(prefix, predicate)

    # Search for absolute paths in this prefix
    all_ok &= check_absolute_paths(prefix, all_files; silent=silent)

    # Search for case-sensitive ambiguities
    all_ok &= check_case_sensitivity(prefix)
    return all_ok
end

"""
    compatible_marchs(p::AbstractPlatform)

Return a (sorted) list of compatible microarchitectures, starting from the most
compatible to the most highly specialized.  If no microarchitecture is specified within
`p`, returns the most generic microarchitecture possible for the given architecture.
"""
function compatible_marchs(platform::AbstractPlatform)
    if !haskey(BinaryPlatforms.arch_march_isa_mapping, arch(platform))
        throw(ArgumentError("Architecture \"$(arch(platform))\" does not contain any known microarchitectures!"))
    end

    # This is the list of microarchitecture names for this architecture
    march_list = String[march for (march, features) in BinaryPlatforms.arch_march_isa_mapping[arch(platform)]]

    # Fast-path a `march()` of `nothing`, which defaults to only the most compatible microarchitecture
    if march(platform) === nothing
        return march_list[begin:begin]
    end

    # Search for this platform's march in the march list
    idx = findfirst(m -> m == march(platform), march_list)
    if idx === nothing
        throw(ArgumentError("Microarchitecture \"$(march(platform))\" not valid for architecture \"$(arch(platform))\""))
    end

    # Return all up to that index
    return march_list[begin:idx]
end

function check_isa(oh, platform, prefix;
                   verbose::Bool = false,
                   silent::Bool = false)
    detected_march = analyze_instruction_set(oh, platform; verbose=verbose)
    platform_marchs = compatible_marchs(platform)
    if detected_march ∉ platform_marchs
        # The object file is for a microarchitecture not compatible with
        # the desired one
        if !silent
            msg = replace("""
            Minimum instruction set detected for $(relpath(path(oh), prefix.path)) is
            $(detected_march), not $(last(platform_marchs)) as desired.
            """, '\n' => ' ')
            @warn(strip(msg))
        end
        return false
    elseif detected_march != last(platform_marchs)
        # The object file is compatible with the desired
        # microarchitecture, but using a lower instruction set: inform
        # the user that they may be missing some optimisation, but there
        # is no incompatibility, so return true.
        if !silent
            msg = replace("""
            Minimum instruction set detected for $(relpath(path(oh), prefix.path)) is
            $(detected_march), not $(last(platform_marchs)) as desired.
            You may be missing some optimization flags during compilation.
            """, '\n' => ' ')
            @warn(strip(msg))
        end
    end
    return true
end

function check_dynamic_linkage(oh, prefix, bin_files;
                               platform::AbstractPlatform = HostPlatform(),
                               verbose::Bool = false,
                               silent::Bool = false,
                               autofix::Bool = true,
                               src_name::AbstractString = "",
                               )
    all_ok = true
    # If it's a dynamic binary, check its linkage
    if isdynamic(oh)
        rp = RPath(oh)

        if verbose
            @info("Checking $(relpath(path(oh), prefix.path)) with RPath list $(rpaths(rp))")
        end

        # Look at every dynamic link, and see if we should do anything about that link...
        libs = find_libraries(oh)
        ignored_libraries = String[]
        for libname in keys(libs)
            if should_ignore_lib(libname, oh)
                push!(ignored_libraries, libname)
                continue
            end

            # If this is a default dynamic link, then just rewrite to use rpath and call it good.
            if is_default_lib(libname, oh)
                if autofix
                    if verbose
                        @info("Rpathify'ing default library $(libname)")
                    end
                    relink_to_rpath(prefix, platform, path(oh), libs[libname]; verbose, subdir=src_name)
                end
                continue
            end

            if !isfile(libs[libname])
                # If we couldn't resolve this library, let's try autofixing,
                # if we're allowed to by the user
                if autofix
                    # First, is this a library that we already know about?
                    known_bins = lowercase.(basename.(bin_files))
                    kidx = findfirst(known_bins .== lowercase(basename(libname)))
                    if kidx !== nothing
                        # If it is, point to that file instead!
                        new_link = update_linkage(prefix, platform, path(oh), libs[libname], bin_files[kidx]; verbose, subdir=src_name)

                        if verbose && new_link !== nothing
                            @info("Linked library $(libname) has been auto-mapped to $(new_link)")
                        end
                    else
                        if !silent
                            @warn("Linked library $(libname) could not be resolved and could not be auto-mapped")
                            if is_troublesome_library_link(libname, platform)
                                @warn("Depending on $(libname) is known to cause problems at runtime, make sure to link against the JLL library instead")
                            end
                        end
                        all_ok = false
                    end
                else
                    if !silent
                        @warn("Linked library $(libname) could not be resolved within the given prefix")
                    end
                    all_ok = false
                end
            elseif !startswith(libs[libname], prefix.path)
                if !silent
                    @warn("Linked library $(libname) (resolved path $(libs[libname])) is not within the given prefix")
                end
                all_ok = false
            end
        end

        if verbose && !isempty(ignored_libraries)
            @info("Ignored system libraries $(join(ignored_libraries, ", "))")
        end

        # If there is an identity mismatch (which only happens on macOS) fix it
        if autofix
            fix_identity_mismatch(prefix, platform, path(oh), oh; verbose, subdir=src_name)
        end
    end
    return all_ok
end

if VERSION < v"1.5-DEV"
function walkdir_nosymlinks(path::String)
    function adjuster(out_c::Channel)
        for (root, dirs, files) in walkdir(path)
            for d in dirs
                if islink(joinpath(root, d))
                    push!(files, d)
                end
            end
            filter!(d -> !islink(joinpath(root, d)), dirs)
            put!(out_c, (root, dirs, files))
        end
    end
    return Channel(adjuster)
end
else
# No adjustment necessary on 1.5+
walkdir_nosymlinks(path) = walkdir(path)
end
"""
    collect_files(path::AbstractString, predicate::Function = f -> true)

Find all files that satisfy `predicate()` when the full path to that file is
passed in, returning the list of file paths.
"""
function collect_files(path::AbstractString, predicate::Function = f -> true;
                       exclude_externalities::Bool = true, exclude_dirs::Bool = true)
    # Sometimes `path` doesn't actually live where we think it does, so canonicalize it immediately
    path = Pkg.Types.safe_realpath(path)

    if !isdir(path)
        return String[]
    end
    # If we are set to exclude externalities, then filter out symlinks that point
    # outside of our given `path`.
    if exclude_externalities
        old_predicate = predicate
        predicate = f -> old_predicate(f) && !(islink(f) && !startswith(Pkg.Types.safe_realpath(f), path))
    end
    collected = String[]
    for (root, dirs, files) in walkdir_nosymlinks(path)
        if exclude_dirs
            list = files
        else
            list = append!(files, dirs)
        end
        for f in list
            f_path = joinpath(root, f)

            if predicate(f_path)
                push!(collected, f_path)
            end
        end
    end
    return collected
end
# Unwrap Prefix objects automatically
collect_files(prefix::Prefix, args...; kwargs...) = collect_files(prefix.path, args...; kwargs...)

"""
    collapse_symlinks(files::Vector{String})

Given a list of files, prune those that are symlinks pointing to other files
within the list.
"""
function collapse_symlinks(files::Vector{String})
    abs_files = String[]
    # Collapse symlinks down to real files, but don't die if we've got a broken symlink
    for f in files
        try
            push!(abs_files, realpath(f))
        catch
        end
    end

    # Return true if it's not a link and the real path exists in abs_files
    predicate = f -> begin
        try
            return !(islink(f) && realpath(f) in abs_files)
        catch
            return false
        end
    end
    return filter(predicate, files)
end

"""
    check_license(prefix, src_name; verbose::Bool = false,, silent::Bool = false)

Check that there are license files for the project called `src_name` in the `prefix`.
"""
function check_license(prefix::Prefix, src_name::AbstractString = "";
                       verbose::Bool = false, silent::Bool = false)
    if verbose
        @info("Checking license file")
    end
    license_dir = joinpath(prefix.path, "share", "licenses", src_name)
    if isdir(license_dir) && length(readdir(license_dir)) >= 1
        if verbose
            @info("Found license file(s): " * join(readdir(license_dir), ", "))
        end
        return true
    else
        if !silent
            @error("Unable to find valid license file in \"\${prefix}/share/licenses/$(src_name)\"")
        end
        # This is pretty serious; don't let us get through without a license
        return false
    end
end

end # module Auditor


text:
module TestUtils

using FLoops: _any
using Test

function test_any()
    tf = [false, true]
    @testset for arity in 0:12
        combinations = Iterators.product((tf for _ in 1:arity)...)
        for bools in combinations
            @test _any(bools...) == any(bools)
        end
    end
end

end  # module


text:
# This file was generated, do not modify it. # hide
mutable struct CompositeModel2 <: DeterministicNetwork
    std_model::Standardizer
    box_model::UnivariateBoxCoxTransformer
    ridge_model::RidgeRegressor
end

function MLJ.fit(m::CompositeModel2, verbosity::Int, X, y)
    Xs = source(X)
    ys = source(y)
    W = MLJ.transform(machine(m.std_model, Xs), Xs)
    box = machine(m.box_model, ys)
    z = MLJ.transform(box, ys)
    ẑ = predict(machine(m.ridge_model, W, z), W)
    ŷ = inverse_transform(box, ẑ)
    mach = machine(Deterministic(), Xs, ys; predict=ŷ)
    fit!(mach, verbosity=verbosity - 1)
    return mach()
end

mdl = CompositeModel2(Standardizer(), UnivariateBoxCoxTransformer(),
                      RidgeRegressor(lambda=0.1))
cm = machine(mdl, X, y)
res = evaluate!(cm, resampling=Holdout(fraction_train=0.8), measure=rms)
round(res.measurement[1], sigdigits=3)


text:
ENV["GKSwstype"] = "100" #src

# # Getting Started with Data Science in Julia
#
# > [A]lways bite the bullet with regard to understanding the basics of the data first before you do anything else...
# >
# > \- Caitlin Smallwood, Vice President of Science and Algorithms at Netflix

println("Hello, Data Science World!")

# # What does your data look like?
#
# Before you can solve complex problems with data, you should have a firm grasp on what your data looks like.  Are your variables continuous or categorical?  What do the distributions look like?  Are there missing observations?  Are variables correlated?  All of these questions arise in the data science workflow, so it's helpful to know the answers from the start of your analysis.

# Let's start with loading the [**`Statistics`**](https://docs.julialang.org/en/latest/stdlib/Statistics/) package, which provides basic statistical operations like means, variances, etc.

using Statistics

# You can examine the items that a package exports with the `names` function.

names(Statistics)

# To start, let's try some of the functions from **`Statistics`** on simulated data.

x = randn(100, 3)

mean(x, dims=1)  # calculate over first dimension (column means)

#-

cor(x)  # correlation matrix

# For a more realistic example, let's load some data from the **`RDatasets`** package, which has a large collection of datasets that get loaded as a `DataFrame`.

# The [`"iris"` dataset](https://en.wikipedia.org/wiki/Iris_flower_data_set) is a collection of measurements for three different species of the iris flower: Iris Setosa, Iris Virginica, and Iris Versicolor (shown below).  The measurements consist of length and width of the petal and sepal (part underneath the flower).

#

using RDatasets

iris = dataset("datasets", "iris")

# DataFrames are discussed in another course module.  For now, we'll just use the fact that data vectors from the DataFrame can be selected with `mydf.mycol`:

iris.SepalLength

# What is the average sepal width across all three species?

mean(iris.SepalWidth)

# Maximum and minimum?

@show minimum(iris.SepalWidth)
@show maximum(iris.SepalWidth)
extrema(iris.SepalWidth)

# Is petal width correlated with petal length?

cor(iris.PetalWidth, iris.PetalLength)

# While we could examine each column separately, a much quicker way to summarize our variables is with the `describe` function:

describe(iris)

# # Random Sampling
#
#
#
# [Random Sampling](https://en.wikipedia.org/wiki/Sampling_(statistics)) plays an integral part in many data science tasks, such as:
#
# - Splitting data into multiple datasets for [cross validation](https://en.wikipedia.org/wiki/Cross-validation_(statistics)).
# - Subsampling a large dataset to something more manageable.
# - Running statistical simulations.
# - [Statistical Bootstrap](https://en.wikipedia.org/wiki/Bootstrapping_(statistics))
#
# The **`StatsBase`** package has the `sample` function to factilitate sampling [with or without replacement](https://en.wikipedia.org/wiki/Sampling_(statistics)#Replacement_of_selected_units).
# Imagine you had a bag with 20 numbers in it.
# Sampling with replacement means that you put the number you drew back into the bag before the next draw.
# Just like rolling a dice, each draw is _independent_ and a given number might get drawn again.
# On the other hand, if you were to _not_ put the drawn number back in the bag, that would be representative of sampling without replacement.
# Each draw is _dependent_ upon all previous draws as those numbers are no longer available for selection.
# For bootstrapping and statistical simulations, it's typical to use sampling with replacement.
# While the two behave nearly the same if the sample is a very small fraction of the overall population, it becomes important to sample without replacement if the sample is large or if using uneven importance weights.

using StatsBase

y1 = sample(1:20, 20, replace=true)  # replaces units after being selected

#-

y2 = sample(1:20, 20, replace=false)  # DOES NOT replace units after being selected

# Next we'll use the [`countmap`](http://juliastats.github.io/StatsBase.jl/latest/search.html?q=countmap) function to get the count of each unique value in our samples.  Notice that when we do sampling *without* replacement, the values will only appear one time.  When sampling *with* replacement, any given selected item is then placed back into the pool of possible items to select from.

countmap(y1)

#-

countmap(y2)

# This is just a small sample (pun intended) of the features in **`StatsBase`**.  For a more complete look of what is possible, check out the [documentation](http://juliastats.github.io/StatsBase.jl/latest/index.html).

# # Parametric Distributions
#
# The [**`Distributions`**](https://github.com/JuliaStats/Distributions.jl) package provides an interface for working with probability distributions.  The full documentation is [here](https://juliastats.github.io/Distributions.jl/stable/).

# Here we'll also load the **`StatsPlots`** package (more on this in the next course module) to visualize the probability distributions.

using Distributions, StatsPlots

plot(Normal(), label = "Normal(0, 1)")

#-

plot!(Gamma(5, 1), label = "Gamma(5, 1)")

# The package defines a number of functions that together create a consistent "grammar" for discussing distributions:
#
# - Probability density function: `pdf`
# - Cumulative distribution function: `cdf`
# - The inverse cumulative distribution function with the standard `quantile` statistical function
# - etc.

d = Normal()

pdf(d, 0), cdf(d, 0), mean(d), var(d), quantile(d, .5), mode(d)

#-

d = Gamma(5, 1)

pdf(d, 0), cdf(d, 0), mean(d), var(d), quantile(d, .5), mode(d)

# As an example, let's write a function that uses [Newton's Method](https://en.wikipedia.org/wiki/Newton%27s_method) that can find a given quantile for any continuous univariate distribution.  Newton's method attempts to find the root for a function $f$ by performing iterations of the form:
#
# $$\theta^{(t)} = \theta^{(t-1)} - \frac{f(\theta^{(t-1)})}{f'(\theta^{(t-1)})}.$$
#
# For quantiles, we want to find the root of the function $F(\theta) - q$ where $F$ is the cumulative density function of the distribtuion and $q \in (0, 1)$.  Using **`Distributions`**, this looks something like

function myquantile(d::Distribution, q::Number)
    θ = mean(d)
    for i in 1:20
        θ -= (cdf(d, θ) - q) / pdf(d, θ)  # θ = θ - (F(θ) - q) / F'(θ)
    end
    θ
end

# Does our `myquantile` function work as expected?  Let's try out our function on several distributions.

d = Normal()
myquantile(d, .5), quantile(d, .5)

#-

d = Gamma(4,3)
myquantile(d, .7), quantile(d, .7)

# The above example shows off the power of generic functions.  Instead of hard-coding the distribution (as would be necessary in R), we can write functions in terms of an arbitrary distribution (without extra effort).  This gives us a lot of flexibility for tasks such as writing [Gibbs Samplers](https://en.wikipedia.org/wiki/Gibbs_sampling) that can swap out distributions with ease.

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
using Test

@testset "./src/LocalPolytope/utils.jl" begin

using BellScenario

@testset "dimension()" begin

    @testset "vectors" begin
        vertices = [[1,0,0,0],[0,1,0,0],[0,0,1,0],[1,1,0,0]]
        @test LocalPolytope.dimension(vertices) == 3
    end

    @testset "matrices" begin
        vertices = [[1 0;0 0],[0 0;0 1],[1 1;0 0],[1 1;0 1]]
        @test LocalPolytope.dimension(vertices) == 3
    end

    @testset "Deterministic Strategies" begin
        vertices = DeterministicStrategy.([[1 1;0 0],[0 0;1 1],[1 0;0 1]])
        @test LocalPolytope.dimension(vertices) == 2
    end
end

end


text:
#Iterate through samps and circ list, recording the indices at which the samp vector crosses the thresholds in circ_list
#circ_list is expected to be sorted in either forward or reverse order, and the first index of samps is expected to be outside the range spanned by first(circ_list) and last(circ_list)
#After fully iterating through circ_list the search will continue on the reversed circ_list.
#Since the endpoints of circ_list get repeated they can be counted too often due to noise.  For this reason the pad_samps argument determines how long to wait before counting another threshold crossing each time the list is reversed
#Returns a vector-of-vectors where each vector contains the crossing indices of a particular circ_list value
function find_circular(samps::AbstractVector, circ_list::AbstractVector, pad_nsamps::Int)
    is_incr = true
    reversed = false
    rlist = reverse(circ_list)
    if !issorted(circ_list)
        if !issorted(rlist)
            error("circ_list must be in increasing or decreasing order")
        end
        circ_list = rlist
        is_incr = false
        reversed = true
    end
    samp1 = first(samps)
    if (is_incr && (samp1 > circ_list[1])) || (!is_incr && (samp1 < circ_list[end]))
        error("The sample vector begins in the middle of a cycle")
    end
    outputs = Vector{Vector{Int}}(undef, length(circ_list))
    for i = 1:length(circ_list)
        outputs[i] = Int[]
    end
    wait_stop = 0
    circ_i = reversed ? length(circ_list) : 1
    for i = 1:length(samps)
        if i == wait_stop
            if (is_incr && (samps[i] >= circ_list[circ_i])) || (!is_incr && (samps[i] <= circ_list[1]))
                error("The signal had already passed threshold after waiting for pad_nsamps samples")
            end
        elseif i > wait_stop
            if is_incr && (samps[i] >= circ_list[circ_i])
                push!(outputs[circ_i], i)
                if circ_i == length(circ_list)
                    is_incr = false
                    wait_stop = i + pad_nsamps
                else
                    circ_i += 1
                end
            elseif !is_incr && (samps[i] <= circ_list[circ_i])
                push!(outputs[circ_i], i)
                if circ_i == 1
                    is_incr = true
                    wait_stop = i + pad_nsamps
                else
                    circ_i -= 1
                end
            end
        end
    end
    if reversed
        outputs = reverse(outputs)
    end
    return outputs
end

#Returns a vector that is the mean of each cycle in "samps" beginning at start_idx
#Setting start_idx greater than one allows skipping of early samples while the piezo warms up
function mean_cycle(samps::AbstractVector, nsamps_cycle::Int; start_idx = 1)
    lsamps = length(samps)
    @assert start_idx >=1 && start_idx <= lsamps
    @assert (lsamps-start_idx+1) >= nsamps_cycle
    rng_analyzed = start_idx:lsamps
    samps_used = view(samps, rng_analyzed)
    if mod(length(samps_used), nsamps_cycle) != 0
        @warn "Found a fractional number of cycles.  Truncating to analyze a whole number of cycles."
        rng_analyzed = start_idx:start_idx+(nsamps_cycle*div(length(rng_analyzed),nsamps_cycle)-1)
        samps_used = view(samps, rng_analyzed)
    end
    cyc_mat = get_cycles(samps_used, nsamps_cycle)
    print("Averaging $(size(cyc_mat,2)) cycles...\n")
    return reshape(mean(ustrip.(cyc_mat), dims=2), nsamps_cycle)  * unit(samps[1]) #work around limitation of mean function with Unitful Quantities
end

#find integer sample shift to align mod and mon.  (Even with calibration there is often a phase delay between the MOD and MON signal)
#We do this so that we know when the MON signal reaches its minimum
#indmax(xcorr(a,b)) is displaced from the center of the xcorr vector by the amount that b would need to be shifted to align with a
function mon_delay(mod_cyc, mon_cyc)
    xc = xcorr(repeat(ustrip.(mon_cyc), 3), ustrip.(mod_cyc)) #replicate it so that edges don't corrupt results
    ctr_i = div(length(xc),2)+1
    half_window_sz = div(length(mon_cyc),4)
    inner_half = xc[(ctr_i-half_window_sz):(ctr_i+half_window_sz)]
    #amount that mod_cyc needs to be shifted to align with mon_cyc (should be positive)
    #we examine only the inner half of shifts because we don't want to shift more than a half cycle
    return argmax(inner_half) - (div(length(inner_half),2) + 1)
end

#Returns a set of indices marking when the mon signal crosses the values in the crossings vector.
#Useful for deciding when to time camera exposures / laser pulses
#note: if calibration is needed it must be done before calling this function. i.e. sig should be calibrated already
function pulse_timings(mod_cyc::AbstractVector, mon_cyc::AbstractVector, targets::AbstractVector, pad_nsamps::Int)
    @assert length(mod_cyc) == length(mon_cyc)
    nsamps_cycle = length(mod_cyc)
    mon_idxs = find_circular(mon_cyc, targets, pad_nsamps)
    @assert all(map(length, mon_idxs).==2) #should have exactly two crossings for each target with a cyclical waveform
    return mon_idxs
end

function flash_cam_cycs(mod_cyc::AbstractVector, mon_cyc::AbstractVector, slice_zs, pad_nsamps::Int, flash_time::HasTimeUnits, exp_time::HasTimeUnits, sample_rate::HasInverseTimeUnits; is_bidi = true)
    offset = mon_delay(ustrip.(mod_cyc), ustrip.(mon_cyc))
    mon_cyc = circshift(mon_cyc, -offset)
    @assert flash_time <= exp_time
    flash_ctr_is = pulse_timings(mod_cyc, mon_cyc, slice_zs, pad_nsamps)
    flash_ctr_fwd = [x[1] for x in flash_ctr_is]
    flash_ctr_back = reverse([x[2] for x in flash_ctr_is]) #reversed so that corresponding z locations line up
    if is_bidi
        flash_ctr_is = vcat(flash_ctr_fwd, flash_ctr_back)
    else
        flash_ctr_is = flash_ctr_fwd
    end
    flash_cam_cycs(flash_ctr_is, flash_time, exp_time, sample_rate, length(mon_cyc), offset)
end

function flash_cam_cycs(flash_ctr_is, flash_time::HasTimeUnits, exp_time::HasTimeUnits, sample_rate::HasInverseTimeUnits, nsamps_cyc::Int, offset::Int)
    #@assert all(map((x,y) -> isapprox(x, y, rtol = (1/typemax(Int16) * 800μm)), (mean_cyc[flash_ctr_fwd], mean_cyc[flash_ctr_back])))
    min_flash_sep = minimum(diff(flash_ctr_is))
    max_exp_time = ((min_flash_sep-1) / sample_rate) * 0.97
    #@show max_exp_time = (1/stack_rate - (1/samprate(pos) *length(slice_zs))) / length(slice_zs) #leaves one low sample between exposure pulses
    if exp_time >= max_exp_time
        @warn "Could not use the requested exposure time while keeping sufficient separation between exposures.  Using an exposure time of $(max_exp_time) instead.  This means that the maximum framerate required from the camera will be $(100.0 * inv(max_exp_time)/inv(exp_time))% of the requested rate. (choose the ROI size to meet this requirement)"
        exp_time = max_exp_time
    end
    #now place flashes, centered on the indices above
    nsamps_flash = ImagineInterface.calc_num_samps(flash_time, sample_rate)
    if iseven(nsamps_flash)
	    nsamps_flash -= 1 #force odd so that fwd and reverse flashes line up in bidi recordings
    end
    nsamps_from_ctr = div(nsamps_flash,2)
    flash_ivs = map(x-> ClosedInterval(x-nsamps_from_ctr, x+nsamps_from_ctr), flash_ctr_is)
    #now place exposures, aligning exposure ends with flash ends
    nsamps_exp = ImagineInterface.calc_num_samps(exp_time, sample_rate)
    exp_ivs= map(x-> ClosedInterval(x+nsamps_from_ctr-nsamps_exp+1, x+nsamps_from_ctr), flash_ctr_is)
    flash_cyc = ImagineInterface.gen_pulses(nsamps_cyc, flash_ivs)
    cam_cyc = ImagineInterface.gen_pulses(nsamps_cyc, exp_ivs)

    return flash_cyc, cam_cyc, offset
end

#empirical max acceleration
max_acc(v, sr) = maximum(abs.(upreferred.(diff(diff(v)./(1/sr))./(1/sr))))
max_vel(v, sr) = maximum(abs.(upreferred.(diff(v)./(1/sr))))

#Returns a triangle wave equivalent to the input (same max, min, freq), useful for comparisons
function eq_triangle(v, sr::HasInverseTimeUnits)
    mx = maximum(v)
    mn = minimum(v)
    freq = upreferred(sr/length(v))
    fwd, bck = ImagineInterface.gen_bidi_pos(mn, mx, upreferred(1/(2.0*freq)), sr)
    return vcat(fwd,bck)
end

#calculate the fraction of the camera's max framerate that can be utilized with the given sample vector
#by comparing the maximum velocity reached with the ideal constant velocity of a triangle wave
function framerate_efficiency(v, sr::HasInverseTimeUnits)
    tri = eq_triangle(v, sr)
    vin = max_vel(v,sr)
    vtri = max_vel(tri,sr)
    return upreferred(vtri/vin)
end

#calculate the fraction of the camera's max framerate that can be utilized with the given sample vector and slice locations
#This version is more practical than above (above is more of a worst-case efficiency)
function framerate_efficiency(v, sr::HasInverseTimeUnits, slice_zs; pad_nsamps=ImagineInterface.calc_num_samps(0.0005s, sr), toffsets_fwd=fill(0.0s, length(slice_zs)), toffsets_bck=fill(0.0s, length(slice_zs)))
    mexp = max_exp(v,sr,slice_zs; pad_nsamps=pad_nsamps, toffsets_fwd=toffsets_fwd, toffsets_bck=toffsets_bck)
    tri = eq_triangle(v, sr)
    mexp_tri = max_exp(tri, sr, slice_zs; pad_nsamps=2)
    @assert eltype(v) == eltype(slice_zs)
    return upreferred(mexp/mexp_tri)
end

#calculate longest exposure time possible given sample vector v
#Also considers temporal offsets that have been found empirically (toffsets_fwd and toffsets_bck keyword args)
function max_exp(v, sr::HasInverseTimeUnits, slice_zs; pad_nsamps = ImagineInterface.calc_num_samps(0.0005s, sr), toffsets_fwd=fill(0.0s, length(slice_zs)), toffsets_bck = fill(0.0s,length(slice_zs)))
    first_idx = argmin(v)
    v = circshift(v, -first_idx) #assumes slice_zs are sorted and increasing
    timings = ImagineAnalyses.find_circular(v, slice_zs, pad_nsamps)
    fwdt = [x[1] for x in timings]
    bckt = reverse([x[2] for x in timings])
    toffsets_fwd_nsamps = map(x->ImagineInterface.calc_num_samps(x, sr), toffsets_fwd)
    toffsets_bck_nsamps = map(x->ImagineInterface.calc_num_samps(x, sr), toffsets_bck)
    fwdt = fwdt.-toffsets_fwd_nsamps
    bckt = bckt.-reverse(toffsets_bck_nsamps)
    diffi = minimum(diff(vcat(fwdt, bckt))) - 1 #subtract one because we need a sample between exposures (to send TTL down)
    diffi = min(diffi, first(fwdt)+(length(v)-last(bckt)) - 1) #handle first/last timing
    return upreferred(diffi / sr)
end


text:
mutable struct Lattice
    n::Int64 # number of rows in basis matrix
    m::Int64 # number of columns in basis matrix
    basis::Array{Float64} # matrix whose columns form basis of lattice
    μ::Array{Float64} # lower triangular matrix of gram-schmidt coefficients
    r::Array{Float64} # squared norm of gram-schmidt basis vectors
    gs::Array{Float64} # columns are gram-schmidt basis vectors

    """
        Lattice(basis)

    Construct a Lattice generated by the input basis matrix.
    """
    Lattice(basis) = gram_schmidt(basis);
    Lattice(n, m, basis, μ, r, gs) = new(n, m, basis, μ, r, gs);
end

"""
    gram_schmidt(M)

Compute Gram-Schmidt orthogonalization on columns of input matrix.
"""
function gram_schmidt(M::Array{<:Real})
    n = size(M, 1); # number of rows of matrix
    m = size(M, 2); # number of columns of matrix
    μ = zeros(Float64,m,m); # lower triangular matrix of gram-schmidt coefficients
    r = zeros(Float64,m); # squared norm of gram-schmidt basis vectors
    gs = zeros(Float64,n,m); # columns are gram-schmidt basis vectors
    for i = 1:m
        b = M[:,i];
        for j = 1:i-1 # compute gram-schmidt coefficients for j < i
            μ[i,j] = (transpose(b)*gs[:,j])/r[j];
        end
        for j = 1:i-1 # compute i-th gram-schmidt basis vector
            b -= μ[i,j]*gs[:,j];
        end
        gs[:,i] = b;
        r[i] = norm(b)^2; # squared norm of gram-schmidt vector
        μ[i,i] = (transpose(b)*M[:,i])/r[i];
    end
    return Lattice(n, m, M, μ, r, gs);
end

"""
    update_gram_schmidt!(L)

Update Gram-Schmidt coefficients and vectors.
"""
function update_gram_schmidt!(L::Lattice)
    for i = 1:L.m
        b = L.basis[:,i];
        for j = 1:i-1 # compute gram-schmidt coefficients for j < i
            L.μ[i,j] = (transpose(b)*L.gs[:,j])/L.r[j];
        end
        for j = 1:i-1 # compute i-th gram-schmidt basis vector
            b -= L.μ[i,j]*L.gs[:,j];
        end
        L.gs[:,i] = b;
        L.r[i] = norm(b)^2; # squared norm of gram-schmidt vector
        L.μ[i,i] = (transpose(b)*L.basis[:,i])/L.r[i];
    end
    return L;
end

"""
    size_reduce(L)

Compute size reduction of basis vectors of lattice, such that all Gram-Schmidt coefficients have absolute value at most 1/2. Modify the Gram-Schmidt coefficients of input lattice accordingly.
"""
function size_reduce!(L::Lattice)
    for i =2:L.m
        for j = i-1:-1:1
            L.basis[:,i] -= round(L.μ[i,j])*L.basis[:,j];
            for k = 1:j
                L.μ[i,k] -= round(L.μ[i,j])*L.μ[j,k];
            end
        end
    end
    return L
end

"""
    check_size_reduce(L, tol=1e-5)

Check whether basis vectors of input lattice are size-reduced.
"""
function check_size_reduce(L::Lattice, tol::Real=1e-5)
    for i = 1:L.m, j = 1:i-1
        if abs(L.μ[i,j]) > 0.5 + tol
            return false
        end
    end
    return true
end

"""
    lll!(L, δ=0.75)

Compute LLL reduction of input lattice, with parameter δ such that ``0.25 < δ < 1``.
"""
function lll!(L::Lattice, δ::Real=0.75)
    k = 2;
    while k <= L.m
        size_reduce!(L);
        if norm(L.gs[:,k])^2 >= (δ - L.μ[k,k-1]^2)*norm(L.gs[:,k-1])^2
            k += 1;
        else
            b = copy(L.basis[:,k]);
            L.basis[:,k] = copy(L.basis[:,k-1]);
            L.basis[:,k-1] = b;
            update_gram_schmidt!(L);
            k = max(2, k-1);
        end
    end
    return L;
end

"""
    lll(M, δ=0.75)

Compute LLL reduction of input matrix, with parameter δ such that ``0.25 < δ < 1``.
"""
function lll(M::Array{<:Real}, δ::Real=0.75)
    return lll(Lattice(M)).basis;
end

"""
    check_lll(L, δ=0.75)

Check whether input lattice has an LLL-reduced basis, with parameter δ such that ``0.25 < δ < 1``.
"""
function check_lll(L::Lattice, δ::Real=0.75)
    if !check_size_reduce(L)
        return false;
    end
    for i=1:L.m-1
        if L.r[i+1] < (δ-L.μ[i+1,i]^2)*L.r[i]
            return false;
        end
    end
    return true;
end

"""
    svp(M)

Compute shortest vector in lattice generated by M, using Schnorr-Euchner enumeration.
"""
function svp(M::Array{<:Real})
    L = Lattice(M);
    x = zeros(Int64, L.m);
    x[L.m] = 1;
    A = norm(M[:,L.m])^2
    sol, length =  enum_svp(L, x, L.m, A, 0.0);
    vect = sum([sol[j]*L.basis[:,j] for j=1:L.m]);
    return [vect, length, sol];
end

"""
    enum_svp(L, x, i, radius, l)

Enumerate all lattice points in subtree at height i, where parent nodes are defined by x. Helper method for [`svp`](@ref).
"""
function enum_svp(L::Lattice, x::Array{<:Integer}, i::Integer, radius::Real, l::Real)
    sol = copy(x);
    c = -sum(Float64[x[j]*L.μ[j,i] for j=i+1:L.m]);
    interval = sqrt((radius-l)/L.r[i]); # interval of possibilities for x[i]
    k = floor(Int64, c); # stores value of x[i]
    Δk = 0; # used to do search in zig-zag path
    Δ²k = -1; # used to do search in zig-zag path
    lbound = floor(Int64, c - interval);
    ubound = ceil(Int64, c + interval);
    while k>= lbound & k <= ubound
        l_new = l + (k-c)^2*L.r[i];
        if l_new <= radius # possible new solution found
            if i == 1 & ~(l_new == 0) # leaf of enumeration tree
                # l_new is approximation of length of new solution. Do exact computation of length to avoid errors
                old_sol = sol[i];
                sol[i] = k;
                new_length = norm(sum([sol[j]*L.basis[:,j] for j=1:L.m]))^2;
                if radius > new_length # prune enumeration tree
                    radius = new_length;
                else
                    sol[i] = old_sol;
                end
            elseif i > 1 # solve subtree enumeration recursively
                old_sol = sol[i];
                sol[i] = k;
                subtree_sol, subtree_length = enum_svp(L, sol, i-1, radius, l_new);
                if subtree_length < radius
                    sol[1:i-1] = subtree_sol[1:i-1];
                    radius = subtree_length;
                else
                    sol[i] = old_sol;
                end
            end
        end
        Δ²k = - Δ²k;
        Δk = Δ²k - Δk;
        k += Δk;
    end
    return [sol, radius];
end


text:
"""
	_start_timer()

Same as `using Dates; now()`.
"""
_start_timer() = Dates.now()

"""
	_stop_timer(d)

Same as `using Dates; now() - d` where `d` is a DateTime object.
"""
_stop_timer(d::DateTime) = Dates.now() - d

"""
	get_hash_sha256(var)

Get the hash of the content of `var` using SHA256 algorithm.
"""
function get_hash_sha256(var)::String
	io = IOBuffer();
	serialize(io, var)
	result = bytes2hex(sha256(take!(io)))
	close(io)

	result
end

"""
	printcheckpoint([io::IO], any...; format)

Print `any...` to `io` and flush preceeded by the current DateTime in `format`.

A line starter can be specified by the keyword argument `start`. Default is "● ".
"""
function printcheckpoint(
	io::IO, string::AbstractString...;
	format = "dd/mm/yyyy HH:MM:SS",
	start = "● "
)
	println(io, start, Dates.format(now(), "[ $(format) ] "), string...)
	flush(io)
end
printcheckpoint(string::AbstractString...; kwargs...) = printcheckpoint(stdout, string...; kwargs...)
printcheckpoint(io::IO, any::Any...; kwargs...) = printcheckpoint(io, string(any...); kwargs...)
printcheckpoint(any::Any...; kwargs...) = printcheckpoint(stdout, any...; kwargs...)

"""
	human_readable_time(milliseconds)

Print elapsed time in format "HH:MM:SS".
"""
function human_readable_time(ms::Dates.Millisecond)::String
	result = ms.value / 1000
	seconds = round(Int64, result % 60)
	result /= 60
	minutes = round(Int64, result % 60)
	result /= 60
	hours = round(Int64, result % 24)
	string(string(hours; pad=2), ":", string(minutes; pad=2), ":", string(seconds; pad=2))
end

"""
	human_readable_time_s(milliseconds)

Print elapsed time in seconds.
"""
function human_readable_time_s(ms::Dates.Millisecond)::String
	string(ms.value / 1000)
end

"""
	throw_n_log(message; error_type = ErrorException)

Call `@error msg` before throwning.
"""
function throw_n_log(msg::AbstractString; error_type::Type{<:Exception} = ErrorException)
	@error msg
	throw(error_type(msg))
end

## conversion utilities ##

"""
	_strarg(arg)

A utility function used to generate a String representing the final Expr that will cause
the cached file to be loaded.
This means that the following calls will produce the same strings:

```
x = 10
@scache "my-cute-vector" "./" fill(0, 10)
@scache "my-cute-vector" "./" fill(0, x)
@scache "my-cute-vector" "./" fill(0, 5 + 5)
```

This function is called inside [`@_scache`](@ref) to generate the string that will fill the
column `COMMAND` in the cache record (if generated).
"""
_strarg(arg::Any, top_level::Bool = false) = Expr(:call, :string, arg)
function _strarg(arg::Expr, top_level::Bool = false)
	if top_level && arg.head == :escape
		return _strarg(arg.args[1], true)
	elseif top_level && arg.head == :call
		_args = filter(x -> !(typeof(x) == Expr && (x.head == :parameters || x.head == :kw)), arg.args[2:end])

		_params = filter(x -> typeof(x) == Expr && x.head == :parameters, arg.args[2:end])
		_kw =
			if length(_params) > 0
				vcat([p.args for p in _params]...)
			else
				[]
			end
		append!(_kw, filter(x -> typeof(x) == Expr && x.head == :kw, arg.args[2:end]))

		return Expr(:call, :string,
			"$(arg.args[1])(",
			Expr(:call, :join, _toexpr(_strarg.(_args)), ", "),
			length(_kw) > 0 ? "; " : "",
			length(_kw) > 0 ? Expr(:call, :join, _toexpr(_strarg.(_kw, true)), ", ") : "",
			")"
		)
	elseif top_level && arg.head == :kw
		return Expr(:call, :string, "$(arg.args[1]) = ", _strarg(arg.args[2]))
	elseif top_level && arg.head == :parameters
		return join(_strarg.(arg.args), ", ")
	elseif top_level && arg.head == :...
		return _rem(_rem(_rem(_strarg(arg.args[1]), "^\\("), "\\)\$"), ",\$")
	elseif arg.head == :...
		return _strarg(Expr(:call, :join, Expr(:..., _toexpr(arg.args)), ", "))
	else
		Expr(:call, :string, arg)
	end
end

function _rem(s, reg::AbstractString)
	return Expr(:call, :replace, _strarg(s), Expr(:call, :(=>), Expr(:call, :Regex, reg), ""))
end

"""
	_convert_input(arg)

A utility function used to generate a stable Tuple containing the information that will be
used to generate the hash of the cached file.

The resulting object is a NamedTuple{(:args,:kwargs),Tuple{T,D}} where T is a ordered
Tuple containing the arguments passed to the function called and D is a Dict{Symbol,Any}
containing the keyword arguments passed.

Note: a dictionary is used for the keyword arguments because otherwise the hash would change
based on the their order.
"""
_convert_input(arg::Any, top_level::Bool = false) = arg
function _convert_input(arg::Expr, top_level::Bool = false)
	_splat2pairs(v::AbstractVector) = length(v) == 0 ? [] : _splat2pairs(v[1])
	_splat2pairs(ex::Expr) = Expr(:call, :pairs, ex.args[1])

	if top_level && arg.head == :escape
		return _convert_input(arg.args[1], true)
	elseif top_level && arg.head == :call
		_args = filter(x -> !(typeof(x) == Expr && (x.head == :parameters || x.head == :kw)), arg.args[2:end])

		_params = filter(x -> typeof(x) == Expr && x.head == :parameters, arg.args[2:end])
		_kw =
			if length(_params) > 0
				vcat([p.args for p in _params]...)
			else
				[]
			end
		append!(_kw, filter(x -> typeof(x) == Expr && x.head == :kw, arg.args[2:end]))

		_res = filter(x -> typeof(x) == Expr && x.head == :..., _kw)
		_kw = filter(x -> !(typeof(x) == Expr && x.head == :...), _kw)

		return (args = [_convert_input(a, true) for a in _args], kwargs = Dict{Symbol,Any}(_convert_input.(_kw, true)...), res = _splat2pairs(_res) )
	elseif top_level && arg.head == :kw
		return arg.args[1] => _convert_input(arg.args[2])
	elseif top_level && arg.head == :parameters
		return _convert_input.(arg.args)
	else
		return arg
	end
end

## fast "toexpr" conversions

_toexpr(vec::AbstractVector) = Expr(:vect, vec...)
function _toexpr(d::AbstractDict{Symbol,Any})
	n = length(d)

	ks = Vector{Symbol}(undef, n)
	vs = Vector{Any}(undef, n)

	Threads.@threads for (i, (k, v)) in collect(enumerate(d))
		ks[i] = k
		vs[i] = v
	end

	_toexpr(Expr.(:quote, ks)), _toexpr(vs)
end

## escaping utilities ##

isescaped(ex::Any) = false
isescaped(ex::Expr) = typeof(ex) == Expr && ex.head == :escape

unesc(ex::Any) = ex
unesc(ex::Expr) = ex.args[1]

safeunesc(ex::Any) = ex
safeunesc(ex::Expr) = isescaped(ex) ? unesc(ex) : ex

escdepth(ex::Any) = 0
function escdepth(ex::Expr)
	res::Int = 0
	curr_ex = ex
	while (isescaped(curr_ex))
		curr_ex = unesc(curr_ex)
		res += 1
	end
	return res
end

function esc(ex::Any, n::Integer)
	res = ex
	for i in 1:n
		res = esc(res)
	end
	return res
end

unesc_comp(ex::Any) = ex
function unesc_comp(ex::Expr)
	res = ex
	while (isescaped(res))
		res = unesc(res)
	end
	return res
end


text:
module StochasticIntegrals

using DataFrames
using Dates
using Distributions: Normal, quantile
using LinearAlgebra: Symmetric, cholesky, inv, det, LowerTriangular, diag
using UnivariateFunctions
using Random
using Sobol: SobolSeq, next!
using FixedPointAcceleration

include("4_number_generators.jl")
export NumberGenerator, Mersenne, SobolGen, next!

include("1_main_functions.jl")
export ItoIntegral, get_volatility, get_variance, get_covariance, get_correlation, ItoSet, CovarianceAtDate
export get_draws, get_zero_draws, pdf, make_covariance, log_likelihood, brownians_in_use
include("2_data_conversions.jl")
export to_draws, to_dataframe, to_array
include("3_getConfidenceHypercube.jl")
export get_confidence_hypercube
include("5_conditional_distributions.jl")
export generate_conditioned_distribution

end


text:
using Test
using DifferentialEquations

#Constant T and V Ideal Dilute Liquid
phaseDict = readinput("../src/testing/liquid_phase.yml")
spcs = phaseDict["phase"]["Species"]; #mechanism dictionaries index:  phaseDict[phasename]["Species" or "Reactions"]
rxns = phaseDict["phase"]["Reactions"];
solv = phaseDict["Solvents"][0]
liq = IdealDiluteSolution(spcs,rxns,solv;name="phase",diffusionlimited=true) #Define the phase (how species thermodynamic and kinetic properties calculated)
initialconds = Dict(["T"=>450.0,"P"=>1e5,"V"=>1.0e-6*1e6,"octane"=>6.154e-3*1e6,"oxygen"=>4.953e-6*1e6]) #Set simulation Initial Temp and Pressure
domain,y0 = ConstantTVDomain(phase=liq,initialconds=initialconds,constantspecies=["oxygen"]) #Define the domain (encodes how system thermodynamic properties calculated)
react = Reactor(domain,y0,(0.0,140000.01)) #Create the reactor object

sol = solve(react.ode,CVODE_BDF(),abstol=1e-20,reltol=1e-8); #solve the ode associated with the reactor

spcnames = getfield.(liq.species,:name)
octaneind = findfirst(isequal("octane"),spcnames)
y = sol(32977.61568)
@test y[octaneind]/sum(y) ≈ 0.461599061 rtol=3e-2 #from RMG simulator I believe the slight difference is due to better calculation of diffusion limits in RMS

#Constant T and P Ideal Gas
phaseDict = readinput("../src/testing/superminimal.yml") #load mechanism dictionary
spcs = phaseDict["gas"]["Species"]; #mechanism dictionaries index:  phaseDict[phasename]["Species" or "Reactions"]
rxns = phaseDict["gas"]["Reactions"];

ig = IdealGas(spcs,rxns,name="gas") #Define the phase (how species thermodynamic and kinetic properties calculated)
initialconds = Dict(["T"=>1000.0,"P"=>1e5,"H2"=>0.67,"O2"=>0.33]) #Set simulation Initial Temp and Pressure
domain,y0 = ConstantTPDomain(phase=ig,initialconds=initialconds) #Define the domain (encodes how system thermodynamic properties calculated)

react = Reactor(domain,y0,(0.0,150.11094)) #Create the reactor object
sol = solve(react.ode,CVODE_BDF(),abstol=1e-20,reltol=1e-12); #solve the ode associated with the reactor

spcnames = getfield.(ig.species,:name)
h2ind = findfirst(isequal("H2"),spcnames)
o2ind = findfirst(isequal("O2"),spcnames)
h2oind = findfirst(isequal("H2O"),spcnames)
y = sol(20.44002454)
N = sum(y)
@test y[h2ind]/N ≈ 0.412883111 rtol=1e-5 #from RMG simulator
@test y[o2ind]/N ≈ 0.200419093 rtol=1e-5
@test y[h2oind]/N ≈ 0.386618602 rtol=1e-5

#Constant V adiabatic Ideal Gas
#uses superminimal.yml mechanism
initialconds = Dict(["T"=>1000.0,"P"=>10.0e5,"H2"=>0.67,"O2"=>0.33]) #Set simulation Initial Temp and Pressure
domain,y0 = ConstantVDomain(phase=ig,initialconds=initialconds) #Define the domain (encodes how system thermodynamic properties calculated)

react = Reactor(domain,y0,(0.0,0.101)) #Create the reactor object
sol = solve(react.ode,CVODE_BDF(),abstol=1e-20,reltol=1e-12); #solve the ode associated with the reactor

ts = exp.(range(log(1e-15),length=10000,stop=log(0.1)))
IDT = ts[argmax(diff([sol(t)[end] for t in ts]))] #Ignition Delay Time based on argmax(dTdt(t))

@test IDT ≈ 0.038384723436228063 rtol=1e-5 #from Cantera simulation

#test sensitivity analysis
initialconds = Dict(["T"=>1000.0,"P"=>1e5,"H2"=>0.67,"O2"=>0.33]) #Set simulation Initial Temp and Pressure
domain,y0 = ConstantTPDomain(phase=ig,initialconds=initialconds,sensitivity=true) #Define the domain (encodes how system thermodynamic properties calculated)

rmgdgdk = [[  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00],
 [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,  0.00000000e+00,  0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,  0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,  0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,  0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,  0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,  0.00000000e+00,   0.00000000e+00,
    0.00000000e+00],
 [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,  0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00],
 [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00],
 [ -4.17206256e-17,  -2.65921877e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,  0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,  0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,  -1.28556728e-16,  0.00000000e+00,   2.57113455e-16,
    0.00000000e+00,   0.00000000e+00,  0.00000000e+00,   0.00000000e+00,
    0.00000000e+00],
 [  0.00000000e+00,  -2.65921877e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00],
 [  8.34412512e-17,   2.65921877e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   2.57113455e-16,   0.00000000e+00,  -5.14226910e-16,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,],
 [  0.00000000e+00,   2.65921877e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,],
 [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,],
 [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,],
 [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,],
 [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,   0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
    0.00000000e+00,]] #from RMG

rmgdgdk = hcat(rmgdgdk...)'
N = sum(y0)
V= R*N*domain.T/domain.P
@views Gs = calcenthalpyinternalgibbs(domain.phase,domain.T,domain.P,V)[3]
cs = y0[1:length(domain.phase.species)]/V
C = sum(cs)
kfs,krevs = getkfkrevs(phase=domain.phase,T=domain.T,P=domain.P,C=C,N=N,ns=y0[1:length(domain.phase.species)],Gs=Gs,diffs=[])

dgdk = ratederivative(domain;cs=cs,V=V,T=domain.T,kfs=kfs,krevs=krevs,Us=Array{Float64,1}(),N=1.0,wV=Array{Float64,1}(),Cvave=0.0,sparse=false)

dgdkdif = (dgdk-rmgdgdk)./rmgdgdk
@test all((dgdkdif .< 1e-4) .| isnan.(dgdkdif))


text:
    catch_phrase()::String  = join(map( x -> select_random(x)*" " , data["faker"]["company"]["buzzwords"] ))
    bs()::String = join(map( x -> select_random(x)*" " , data["faker"]["company"]["bs"] ))
    company_suffix()::String = select_random(data["faker"]["company"]["suffix"])
    company_prefix()::String  = select_random(data["faker"]["company"]["prefix"])
    company()::String =executor(data["faker"]["company"]["name"])


text:
using SINDY
ds = ODE1()


text:
# include("../src/LarSurf.jl")

using LarSurf
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm, SparseArrays

threshold = 4000
pth = Pio3d.datasets_join_path("medical/orig/sample-data/nrn4.pklz")
datap = Pio3d.read3d(pth)

data3d = datap["data3d"]
# segmentation = convert(Array{Int8, 2}, data3d .> threshold)
segmentation = data3d .> threshold

(V, filtered_topology), Flin, (Vfull, tmodel) = LarSurf.get_surface_grid(segmentation)
(VV, EV, FV, CV) = tmodel
filteredFV = filtered_topology[1]
Plasm.View((V,[VV, EV, filteredFV]))


text:
println("Fitting Boon model")
M1 = BoonModel(N1)
boon = estfun(M1,p1,Q1)


text: # # No Makie recipes yet, text is better I guess


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

using Markdown
using InteractiveUtils

# ╔═╡ 15b8afb9-ace9-4b5b-894e-dce9dda5bf5f
begin
	using Pkg

	try
		using UnicodePlots

	catch ArgumentError
		Pkg.add("UnicodePlots")
		using UnicodePlots
	end

	using Test
	using PlutoUI
end

# ╔═╡ 45d72214-919c-11eb-2e69-23c0786baacb
md"""
### Functional Programming

Src: [Julia 1.0 Cookbook -  Bogumił Kamiński , Przemysław Szufel - 2018](https://www.packtpub.com/product/julia-1-0-programming-cookbook/9781788998369)
"""

# ╔═╡ 6d80a12a-9649-4438-9a12-36003a57bf86
function deriv(f::Function)::Function
	h = √eps()
	df(x) = (f(x + h) - f(x)) / h
	df
end

# ╔═╡ 5e5c6b37-f3ee-42b1-94cc-686d54171833
md"""
Let us test the `deriv` function. We will use the `UnicodePlots.jl` package to plot a $2×x×x + 5×x - 4$ function and its derivative.
"""

# ╔═╡ b3faae61-dcdc-473d-9a06-cc890ae74bd7
begin
	f(x) = 2x*x + 5x - 4
	x = -5.:3.;
end

# ╔═╡ 75442ce9-cf92-4800-b1aa-bfdd4229224c
begin
	plot = lineplot(x, f.(x), width=45, height=15, canvas=DotCanvas, name="f(x)");
	plot = lineplot!(plot, x, deriv(f).(x), name="f'(x)")
end

# ╔═╡ a2db212a-02f8-493c-a1f2-43c0e63c31fa
md"""
Next, we can construct a function that solves any quadratic equation in the form

$$ax^2 + bx +c = 0$$
"""

# ╔═╡ 415677ce-d5d6-4808-8789-f5fa66d7f7df
function q_solve(f)
	c = f(0.)
	f1 = deriv(f)
	b = f1(0.)
	a = f(1,) - b - c
	d = √(b * b - 4 * a * c)
	((-b - d) / 2a, -(-b + d) / 2a)
end

# ╔═╡ bd8d4613-c666-41e3-b877-e9ed6882fe04
@test q_solve(x -> (x - 1) * (x + 7)) == (-7., -1.)

# ╔═╡ b3ee6944-32c7-4dca-ab51-9f1dc0accd91
@test_throws DomainError q_solve(x -> x * x + 1)  ## no solution in the real domain

# ╔═╡ 6486f3b7-9a29-41d1-8d80-514a9212a174
with_terminal() do
	println(q_solve(x -> x * x + 1 + 0im))
	## (-7.450580707946e-9 - 1.0000000074505im, 7.450580707946e-9 - 1.0000000074505im)
end

# ╔═╡ 75a80d9a-953e-41c0-9569-2f1c08ca16cd
md"""
Our quadratic solver function takes any quadratic function as its argument. Such functions can be presented as $f(x) = a×x^2 + b×x + c$, and we calculate $f(0)$ to get the $c$ value.
Next, we calculate `deriv(f)(0)`, which is $(2×a×x + b)(0)$, and hence obtain a numerical approximation of the value of b.
Finally, we calculate a, and now we can calculate root values of the `f` function.

Since we did not specify types, our function works not only for real numbers but also for complex numbers.
This can be observed when comparing the calls `q_solve(x -> x*x + 1)` and
`q_solve(x -> x*x + 1 + 0im)`.
The first one fails because `x × x + 1` does not have a solution in the real domain. On the other hand - the second one - using complex numbers, works.
"""

# ╔═╡ Cell order:
# ╟─45d72214-919c-11eb-2e69-23c0786baacb
# ╠═15b8afb9-ace9-4b5b-894e-dce9dda5bf5f
# ╠═6d80a12a-9649-4438-9a12-36003a57bf86
# ╟─5e5c6b37-f3ee-42b1-94cc-686d54171833
# ╠═b3faae61-dcdc-473d-9a06-cc890ae74bd7
# ╠═75442ce9-cf92-4800-b1aa-bfdd4229224c
# ╟─a2db212a-02f8-493c-a1f2-43c0e63c31fa
# ╠═415677ce-d5d6-4808-8789-f5fa66d7f7df
# ╠═bd8d4613-c666-41e3-b877-e9ed6882fe04
# ╠═b3ee6944-32c7-4dca-ab51-9f1dc0accd91
# ╠═6486f3b7-9a29-41d1-8d80-514a9212a174
# ╟─75a80d9a-953e-41c0-9569-2f1c08ca16cd

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
module SpikesCorrelations

"CorrelationIndex.jl" |> include
# "SpikesCorrelations.jl" |> include
"TilingCoefficient.jl" |> include
"VPDistance.jl" |> include

# Write your package code here.
export correlation_index, tailing_coefficient, spkd

end

# TODO:
# - set up documentation
# - add tests
# - refactor vp distance, so that it is similar to other functions
# - check out how to write a function that takes other functions as argument, just like in distances library