text:
## This is an example on the usage of PolyChaos for non-intrusive PCE using
## the projection and the regression approach and a comarison of them in a multi-
## variate setting.
## Model: f = ΣX^2, X = {X_1, ..., X_k}, X_i ~ N(0, 1), Y ~ Χ²(k)

using PolyChaos
using LinearAlgebra

### Parameters ###
μ, σ = 0, 1
maxdeg = 3
Nrec = 20
k = 3 # degrees of freedom

### Model equation ###
function model(X)
    @assert length(X) == k "Degrees of freedom is $k, but $(length(X)) inputs given"
    println(X)
    sum(x^2 for x in X)
end

### Setup and compute PCE coefficients of x ###
op = GaussOrthoPoly(maxdeg, Nrec=Nrec, addQuadrature=true)
mop = MultiOrthoPoly([op for i in 1:k], maxdeg)
ind = mop.ind

# Find PCE for all X_i
L = dim(mop)
# assign2multi(): take coefficients x for univariate elements and assign then to their accoring multi index.
# If x has 1 element, only take degree 1 basis elements, if it has 2 take the first 2 etc.
x = [ assign2multi(convert2affinePCE(μ, σ, op), i, mop.ind) for i in 1:k ]
# we can describe a gaussian RV with ~N(0,1) by PCE coefficients x_0 = 1 and x_1 = 1

# Compute tensors <ϕ_m,ϕ_m> and <ϕ_1,ϕ_2,ϕ_m>
t2 = Tensor(2, mop)
T2 = [t2.get([i,i]) for i = 0:L-1]
t3 = Tensor(3, mop)

# ------- Reference -------
### Intrusive PCE ###
println("\t == Galerking Projection ==")
# x[i][j] - i = which variable, j = coefficient for degree j
y_intr = [ sum( x[i][j1] * x[i][j2] * t3.get([j1-1, j2-1, m-1]) for i = 1:k, j1 = 1:L, j2 = 1:L ) / t2.get([m-1, m-1]) for m in 1:L ]
# y_intr = [ sum( x[i][j1] * x[i][j2] * t3.get([j1-1, j2-1, m-1]) for i = 1:k, j1 = 1:L, j2 = 1:L ) for m in 1:L ]
println(y_intr)
println()

### Analytic ###
# ...

# ------- Non-intrusive PCE -------
### Projection
println("\t == Projection approach ==")

# Perform projection via numerical integration
function computeProjection(i::Int, index::AbstractVector{<:Int})
    g(v::AbstractVector{<:Real}) = model(v) * evaluate(index, v, mop)
    γ = t2.get([i,i]) # Normalization of coefficients (-> ortho normal) TODO: Why?
    # integ = integrate(x1 -> integrate(x2 -> integrate(x3 -> g([x1,x2,x3]), mop.uni[3]), mop.uni[2]), mop.uni[1])
    println("foo")
    v = Vector{Real}()
    integ = multi_integral(v, k, g)
    integ / γ
end

function multi_integral(v::AbstractVector{<:Real}, n::Int, g)
    if n == 1
        integrate(x -> g(push!(v,x)), mop.uni[n])
        v = Vector{Real}()
    else
        integrate(x -> multi_integral(push!(v,x), n-1, g), mop.uni[n])
    end
end

y_proj = [computeProjection(i-1, index) for (i, index) in enumerate(eachrow(ind))]
println.(y_proj)
println()

# # Comparison to intrusive solution
print("Comapre coefficients intrusive PCE <-> projection:")
println(norm(y_intr - y_proj, Inf), "\n")
println()

### Regression
println("\t == Regression approach ==")
# Draw n samples, where N > P has to hold, P = basis size
nSamples = maxdeg * k * 10
X = sampleMeasure(nSamples, mop)

# Evaluate model -> vector Y
Y = model.(X)

# Build matrix Φ with ϕ(x(i))
# Φ = Array{Float64}(undef, nSamples, maxdegree+1)
# Φ = [ evaluate(j, X[i], op) for i = 1:nSamples, j = 0:maxdeg]

# Ordinary least squares regression
# y_reg = leastSquares(Φ, Y)
# println.(y_reg)
# println()

# Validation of PCE model
# println("Comapre coefficients analytic <-> regression:")
# println(norm(y_ana - y_reg, Inf))

# genError = empError(Y, Φ, y_reg)
# println("Determination coefficient R² (normalized empicial error): ", 1 - genError)

# ϵLoo = looError(Y, Φ, y_reg)
# println("Determination coefficient Q² (leave-one-out error): ", 1- ϵLoo)
# println()

### Monte-Carlo
# Evaluate model function on same set of previously drawn samples X
# y_mc = model.(X)
#println()

# ------- Comparison of moments -------
println("\n\t Comparison of moments to analytic solution:")

# Analytic moments for y
mean_ana = k
std_ana = sqrt(2*k)
skew_ana = sqrt(8/k)

# PCE skewness
function skew(y)
    e3 = sum( y[i] * y[j] * y[k] * t3.get([i-1, j-1, k-1]) for i = 1:L, j = 1:L, k = 1:L )
    μ = y[1]
    σ = std(y, mop)
    (e3 - 3 * μ * σ^2 - μ^3) / (σ^3)
end

# Intrusive PCE moments
println("= Error Intrusive PCE vs Analytic=")
mean_intr = mean(y_intr, mop)
std_intr  = std(y_intr, mop)
skew_intr = skew(y_intr)
error_mean_intr = abs(mean_ana - mean_intr)
error_std_intr  = abs(std_ana - std_intr)
error_skew_intr = abs(skew_ana - skew_intr)
println("\t\t\t error mean: \t $(error_mean_intr)")
println("\t\t\t error std: \t $(error_std_intr)")
println("\t\t\t error skew: \t $(error_skew_intr)\n")

# # MC moments
# mean_mc, std_mc, skew_mc = mean(y_mc), std(y_mc), skewness(y_mc)
# error_mean_mc = abs(mean_ana - mean_mc)
# error_std_mc = abs(std_ana - std_mc)
# println("\t\t\t error MC, mean: \t $(error_mean_mc)")
# println("\t\t\t error MC, std: \t $(error_std_mc)")

# # Projection moments
# mean_proj = mean(y_proj, op)
# std_proj = std(y_proj, op)
# error_mean_proj = abs(mean_ana - mean_proj)
# error_std_proj = abs(std_ana - std_proj)
# println("\t\t\t error proj, mean: \t $(error_mean_proj)")
# println("\t\t\t error proj, std: \t $(error_std_proj)")

# # Regression moments
# mean_reg = mean(y_reg, op)
# std_reg = std(y_reg, op)
# error_mean_reg = abs(mean_ana - mean_reg)
# error_std_reg = abs(std_ana - std_reg)
# println("\t\t\t error reg, mean: \t $(error_mean_reg)")
# println("\t\t\t error reg, std: \t $(error_std_reg)")

# ------- Plotting of PDF -------
using Plots
Nsmpl = 10000
ξ = sampleMeasure(Nsmpl, mop)

import SpecialFunctions: gamma

### Intrusive PCE ###
samp_intr = evaluatePCE(y_intr, ξ, mop)
histogram(samp_intr; normalize=true, xlabel="t", ylabel="ρ(t)")
# Analytic comparison
ρ(t) = 1  / (2^(0.5*k) * gamma(.5*k)) * t^(.5*k-1) * exp(-.5*t) # analytic pdf
t = range(.1; stop=maximum(samp_intr), length=100)
plot!(t, ρ.(t), w=4)

ext: jl

max_stars_repo_name: adriangrupp/sparsePolyChaos

max_issues_repo_name: adriangrupp/sparsePolyChaos

max_forks_repo_name: adriangrupp/sparsePolyChaos


text:
using RealSense.LibRealSense
using Test

err = Ref{Ptr{rs2_error}}(0)

ctx = rs2_create_context(RS2_API_VERSION, err)
@test err[] == C_NULL

device_list = rs2_query_devices(ctx, err)
@test err[] == C_NULL

dev_count = rs2_get_device_count(device_list, err)
@test err[] == C_NULL
println("There are $dev_count connected RealSense devices.")

rs2_delete_device_list(device_list)
rs2_delete_context(ctx)


text:
# Autogenerated wrapper script for LibUV_jll for x86_64-w64-mingw32
export libuv

JLLWrappers.@generate_wrapper_header("LibUV")
JLLWrappers.@declare_library_product(libuv, "libuv-2.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libuv,
        "bin\\libuv-2.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
import Pkg
Pkg.activate(joinpath(pkgdir(EchelleCCFs),"examples"))

using EchelleCCFs
using HDF5, FileIO

include("../deps/download_soap_example_from_gilbertson_etal_2020.jl")

function read_soap_output_gilbertson_ford_dumuseque_2020(filename::String, obs_idx::Integer)
    local time, λ, flux
    c = h5open(filename, "r") do file
        flux_size = size(file["active"])
        @assert 1<= obs_idx <= flux_size[2]
        flux = file["active"][1:flux_size[1],obs_idx]
        λair =  read(file,"lambdas")
        λ = EchelleCCFs.λ_air_to_vac.(λair)
        idx_min = findfirst(x->x>0, flux)
        idx_max = findlast(x->x>0, flux)
        λ = view(λ,idx_min:idx_max)
        flux = view(flux,idx_min:idx_max)
        p_rot_eq = 25.05 # days
        time = file["phases"][obs_idx] * p_rot_eq
    end
    return (time=time, λ=λ, flux = flux)
end

filename = joinpath(pkgdir(EchelleCCFs),"data","spectra","res-1000-1years_full_id1.h5")
(time, λ, flux) = read_soap_output_gilbertson_ford_dumuseque_2020(filename, 1)

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
# Transformer definitions and implementations.
module Transformers

export Transformer,
       Learner,
       OneHotEncoder,
       Imputer,
       Pipeline,
       Wrapper,
       Identity,
       Baseline,
       PrunedTree,
       RandomForest,
       DecisionStumpAdaboost,
       StandardScaler,
       PCA,
       VoteEnsemble,
       StackEnsemble,
       BestLearner,
       SKLLearner,
       CRTLearner,
       fit!,
       transform!

# Obtain system details
import Orchestra.System: LIB_SKL_AVAILABLE, LIB_CRT_AVAILABLE

# Include abstract types as convenience
importall Orchestra.Types

# Include atomic Orchestra transformers
include(joinpath("orchestra", "baseline.jl"))
importall .BaselineMethods
include(joinpath("orchestra", "transformers.jl"))
importall .OrchestraTransformers

# Include Julia transformers
include(joinpath("julia", "decisiontree.jl"))
importall .DecisionTreeWrapper
include(joinpath("julia", "mlbase.jl"))
importall .MLBaseWrapper
include(joinpath("julia", "dimensionalityreduction.jl"))
importall .DimensionalityReductionWrapper

# Include Python transformers
if LIB_SKL_AVAILABLE
  include(joinpath("python", "scikit_learn.jl"))
  importall .ScikitLearnWrapper
end

# Include R transformers
if LIB_CRT_AVAILABLE
  include(joinpath("r", "caret.jl"))
  importall .CaretWrapper
end

# Include aggregate transformers last, dependent on atomic transformers
include(joinpath("orchestra", "ensemble.jl"))
importall .EnsembleMethods

end # module


text:
"""
This function reads the test system time series and returns representative hour weight.
"""
function read_test_system(data_dir::String,
                          test_system_dir::String,
                          base_dir::String,
                          annual_growth::AxisArrays.AxisArray{Float64, 2},
                          start_year::Int64,
                          n_rep_days::Int64)

    test_sys_hour_weight = nothing
    zones = nothing
    representative_days = nothing
    rep_hour_weight = nothing
    system_peak_load = nothing
    test_sys_hour_weight = nothing
    zonal_lines = nothing

    if occursin("RTS-GMLC", test_system_dir)
        zones,
        representative_days,
        rep_hour_weight,
        system_peak_load,
        test_sys_hour_weight,
        zonal_lines = read_rts(data_dir,
                               test_system_dir,
                               base_dir,
                               annual_growth,
                               start_year,
                               n_rep_days)
    end

    return zones, representative_days, rep_hour_weight, system_peak_load, test_sys_hour_weight, zonal_lines
end


text:
using CUDA, Base.Test

#Initialize device
dev = CuDevice(0)

#Create context
ctx = create_context(dev)

#Modules include all functions
md = CuModule("thing.ptx")

#Load function
vadd = CuFunction(md, "vadd")
vmul = CuFunction(md, "vmul")
vsub = CuFunction(md, "vsub")
vdiv = CuFunction(md, "vdiv")

#Init
a = rand(Float32, 10)
b = rand(Float32, 10)
ad = CuArray(a)
bd = CuArray(b)

#Addition
c = zeros(Float32, 10)
cd = CuArray(c)
CUDA.launch(vadd, 10, 1, (ad,bd,cd))
c = to_host(cd)
@test_approx_eq c a+b

#Subtraction
c = zeros(Float32, 10)
cd = CuArray(c)
CUDA.launch(vsub, 10, 1, (ad,bd,cd))
c = to_host(cd)
@test_approx_eq c a-b

#Multiplication
c = zeros(Float32, 10)
cd = CuArray(c)
CUDA.launch(vmul, 10, 1, (ad,bd,cd))
c = to_host(cd)
@test_approx_eq c a.*b

#Division
c = zeros(Float32, 10)
cd = CuArray(c)
CUDA.launch(vdiv, 10, 1, (ad,bd,cd))
c = to_host(cd)
@test_approx_eq c a./b

#Negative test cases
a = rand(Float32, 10)
ad = CuArray(Float32, 5)
@test_throws ArgumentError copy!(ad, a)
@test_throws ArgumentError copy!(a, ad)

#Utility
@test CUDA.ndims(ad) == 1
@test CUDA.eltype(ad) == Float32

#Not sure what this does
push(ctx)
pop(ctx)

#Count number of devices
CUDA.devcount()

#Device info
println("Name of device = $(name(dev))")
println("Device capability = $(capability(dev))")
println("Device memory = $(totalmem(dev))")

#Free memory and finish up
free(ad)
free(bd)
free(cd)
unload(md)
destroy(ctx)


text:
#
# Utilities.
#

#
# Method grouping.
#

"""
$(:SIGNATURES)

Group all methods of function `func` with type signatures `typesig` in module `modname`.
Keyword argument `exact = true` matches signatures "exactly" with `==` rather than `<:`.

# Examples

```julia
groups = methodgroups(f, Union{Tuple{Any}, Tuple{Any, Integer}}, Main; exact = false)
```
"""
function methodgroups(func, typesig, modname; exact = true)
    # Group methods by file and line number.
    local methods = Set{Method}(getmethods(func, typesig))
    local groups = groupby(Tuple{Symbol, Int}, Vector{Method}, methods) do m
        (m.file, m.line), m
    end

    # Filter out methods from other modules and with non-matching signatures.
    local typesigs = alltypesigs(typesig)
    local results = Vector{Method}[]
    for (key, group) in groups
        filter!(group) do m
            local ismod = m.module == modname
            exact ? (ismod && Base.tuple_type_tail(m.sig) in typesigs) : ismod
        end
        isempty(group) || push!(results, group)
    end

    # Sort the groups by file and line.
    sort!(results, lt = comparemethods, by = first)

    return results
end

"""
$(:SIGNATURES)

Compare methods `a` and `b` by file and line number.
"""
function comparemethods(a::Method, b::Method)
    comp = a.file < b.file ? -1 : a.file > b.file ? 1 : 0
    comp == 0 ? a.line < b.line : comp < 0
end

"""
$(:SIGNATURES)

A helper method for [`getmethods`](@ref) that collects methods in `results`.
"""
function getmethods!(results, f, sig)
    if sig == Union{}
        append!(results, methods(f))
    elseif isa(sig, Union)
        for each in sig.types
            append!(results, getmethods(f, each))
        end
    else
        append!(results, methods(f, sig))
    end
    return results
end
"""
$(:SIGNATURES)

Collect and return all methods of function `f` matching signature `sig`.

This is similar to `methods(f, sig)`, but handles type signatures found in `DocStr` objects
more consistently that `methods`.
"""
getmethods(f, sig) = getmethods!(Method[], f, sig)

"""
$(:SIGNATURES)

Is the type `t` a `bitstype`?
"""
isbitstype(t::ANY) = isleaftype(t) && sizeof(t) > 0 && isbits(t)

"""
$(:SIGNATURES)

Is the type `t` an `abstract` type?
"""
isabstracttype(t::ANY) = isa(t, DataType) && getfield(t, :abstract)

"""
$(:SIGNATURES)

Returns a `SimpleVector` of the `Tuple` types contained in `sig`.
"""
alltypesigs(sig) = isa(sig, Union) ? sig.types : Core.svec(sig)

"""
$(:SIGNATURES)

A helper method for [`groupby`](@ref) that uses a pre-allocated `groups` `Dict`.
"""
function groupby!(f, groups, data)
    for each in data
        key, value = f(each)
        push!(get!(groups, key, []), value)
    end
    return sort!(collect(groups), by = first)
end

"""
$(:SIGNATURES)

Group `data` using function `f` where key type is specified by `K` and group type by `V`.

The function `f` takes a single argument, an element of `data`, and should return a 2-tuple
of `(computed_key, element)`. See the example below for details.

# Examples

```julia
groupby(Int, Vector{Int}, collect(1:10)) do num
    mod(num, 3), num
end
```
"""
groupby(f, K, V, data) = groupby!(f, Dict{K, V}(), data)

"""
$(:SIGNATURES)

Remove the `Pkg.dir` part of a file `path` if it exists.
"""
function cleanpath(path::AbstractString)
    local pkgdir = joinpath(Pkg.dir(), "")
    return startswith(path, pkgdir) ? first(split(path, pkgdir; keep = false)) : path
end

"""
$(:SIGNATURES)

Parse all docstrings defined within a module `mod`.
"""
function parsedocs(mod::Module)
    for (binding, multidoc) in Docs.meta(mod)
        for (typesig, docstr) in multidoc.docs
            Docs.parsedoc(docstr)
        end
    end
end

"""
$(:SIGNATURES)

Print a simplified representation of a method signature to `buffer`. Some of these
simplifications include:

  * no `TypeVar`s;
  * no types;
  * no keyword default values;
  * `?` printed where `#unused#` arguments are found.

# Examples

```julia
f(x; a = 1, b...) = x
sig = printmethod(Docs.Binding(Main, :f), f, first(methods(f)))
```
"""
function printmethod(buffer::IOBuffer, binding::Docs.Binding, func, method::Method)
    # TODO: print qualified?
    print(buffer, binding.var)
    print(buffer, "(")
    join(buffer, arguments(method), ", ")
    local kws = keywords(func, method)
    if !isempty(kws)
        print(buffer, "; ")
        join(buffer, kws, ", ")
    end
    print(buffer, ")")
    return buffer
end

printmethod(b, f, m) = takebuf_string(printmethod(IOBuffer(), b, f, m))

# Handle differences between Julia 0.5/0.6 where `LambdaInfo` is replaced with `CodeInfo`.
if isdefined(Base, :LambdaInfo)
    has_method_source(m::Method) = isdefined(m, :lambda_template)
    get_method_source(m::Method) = m.lambda_template
    nargs(m::Method) = m.lambda_template.nargs
else
    has_method_source(m::Method) = isdefined(m, :source)
    get_method_source(m::Method) = m.source
    nargs(m::Method) = m.nargs
end

"""
$(:SIGNATURES)

Returns the list of keywords for a particular method `m` of a function `func`.

# Examples

```julia
f(x; a = 1, b...) = x
kws = keywords(f, first(methods(f)))
```
"""
function keywords(func, m::Method)
    local table = methods(func).mt
    if isdefined(table, :kwsorter)
        local kwsorter = table.kwsorter
        local signature = Base.tuple_type_cons(Vector{Any}, m.sig)
        if method_exists(kwsorter, signature)
            local method = which(kwsorter, signature)
            if has_method_source(method)
                local template = get_method_source(method)
                # `.slotnames` is a `Vector{Any}`. Convert it to the right type.
                local args = map(Symbol, template.slotnames[(nargs(method) + 1):end])
                # Only return the usable symbols, not ones that aren't identifiers.
                filter!(arg -> !contains(string(arg), "#"), args)
                # Keywords *may* not be sorted correctly. We move the vararg one to the end.
                local index = findfirst(arg -> endswith(string(arg), "..."), args)
                if index > 0
                    args[index], args[end] = args[end], args[index]
                end
                return args
            end
        end
    end
    return Symbol[]
end

"""
$(:SIGNATURES)

Returns the list of arguments for a particular method `m`.

# Examples

```julia
f(x; a = 1, b...) = x
args = arguments(first(methods(f)))
```
"""
function arguments(m::Method)
    if has_method_source(m)
        local template = get_method_source(m)
        if isdefined(template, :slotnames)
            local args = map(template.slotnames[1:nargs(m)]) do arg
                arg === Symbol("#unused#") ? "?" : arg
            end
            return filter(arg -> arg !== Symbol("#self#"), args)
        end
    end
    return Symbol[]
end

#
# Source URLs.
#
# Based on code from https://github.com/JuliaLang/julia/blob/master/base/methodshow.jl.
#
# Customised to handle URLs on travis since the directory is not a Git repo and we must
# instead rely on `TRAVIS_REPO_SLUG` to get the remote repo.
#

"""
$(:SIGNATURES)

Get the URL (file and line number) where a method `m` is defined.

Note that this is based on the implementation of `Base.url`, but handles URLs correctly
on TravisCI as well.
"""
url(m::Method) = url(m.module, string(m.file), m.line)

if VERSION < v"0.5.0-dev"
    url(mod, file, line) = ""
else
    function url(mod::Module, file::AbstractString, line::Integer)
        file = is_windows() ? replace(file, '\\', '/') : file
        if Base.inbase(mod) && !isabspath(file)
            local base = "https://github.com/JuliaLang/julia/tree"
            if isempty(Base.GIT_VERSION_INFO.commit)
                return "$base/v$VERSION/base/$file#L$line"
            else
                local commit = Base.GIT_VERSION_INFO.commit
                return "$base/$commit/base/$file#L$line"
            end
        else
            if isfile(file)
                local d = dirname(file)
                return LibGit2.with(LibGit2.GitRepoExt(d)) do repo
                    LibGit2.with(LibGit2.GitConfig(repo)) do cfg
                        local u = LibGit2.get(cfg, "remote.origin.url", "")
                        local m = match(LibGit2.GITHUB_REGEX, u)
                        u = m === nothing ? get(ENV, "TRAVIS_REPO_SLUG", "") : m.captures[1]
                        local commit = string(LibGit2.head_oid(repo))
                        local root = LibGit2.path(repo)
                        if startswith(file, root) || startswith(realpath(file), root)
                            local base = "https://github.com/$u/tree"
                            local filename = file[(length(root) + 1):end]
                            return "$base/$commit/$filename#L$line"
                        else
                            return ""
                        end
                    end
                end
            else
                return ""
            end
        end
    end
end


text:
using LinearAlgebra

function generic_matmatmul!(C::AbstractVecOrMat{R}, A::AbstractVecOrMat{T}, B::AbstractVecOrMat{S}) where {T,S,R}
    if size(A,2) != size(B,1)
        throw(DimensionMismatch("matrix A has dimensions $(size(A)), matrix B has dimensions $(size(B))"))
    end
    if size(C,1) != size(A,1) || size(C,2) != size(B,2)
        throw(DimensionMismatch("result C has dimensions $(size(C)), needs $((size(A,1),size(B,2)))"))
    end
    if isempty(A) || isempty(B)
        return fill!(C, zero(R))
    end

    function kernel(C, A, B)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        j = (blockIdx().y-1) * blockDim().y + threadIdx().y

        if i <= size(A,1) && j <= size(B,2)
            z2 = zero(A[i, 1]*B[1, j] + A[i, 1]*B[1, j])
            Ctmp = convert(promote_type(R, typeof(z2)), z2)
            for k in 1:size(A,2)
                Ctmp += A[i, k]*B[k, j]
            end
            C[i,j] = Ctmp
        end

        return
    end

    max_threads = 256
    threads_x = min(max_threads, size(C,1))
    threads_y = min(max_threads ÷ threads_x, size(C,2))
    threads = (threads_x, threads_y)
    blocks = ceil.(Int, (size(C,1), size(C,2)) ./ threads)

    @cuda threads=threads blocks=blocks kernel(C, A, B)

    C
end

LinearAlgebra.mul!(C::CuVecOrMat, A::CuVecOrMat, B::CuVecOrMat) = generic_matmatmul!(C, A, B)
LinearAlgebra.mul!(C::CuVecOrMat, A::CuVecOrMat, B::LinearAlgebra.Adjoint{<:Any, <:CuVecOrMat}) = generic_matmatmul!(C, A, B)
LinearAlgebra.mul!(C::CuVecOrMat, A::CuVecOrMat, B::LinearAlgebra.Transpose{<:Any, <:CuVecOrMat}) = generic_matmatmul!(C, A, B)
LinearAlgebra.mul!(C::CuVecOrMat, A::LinearAlgebra.Adjoint{<:Any, <:CuVecOrMat}, B::CuVecOrMat) = generic_matmatmul!(C, A, B)
LinearAlgebra.mul!(C::CuVecOrMat, A::LinearAlgebra.Transpose{<:Any, <:CuVecOrMat}, B::CuVecOrMat) = generic_matmatmul!(C, A, B)
LinearAlgebra.mul!(C::CuVecOrMat, A::LinearAlgebra.Transpose{<:Any, <:CuVecOrMat}, B::LinearAlgebra.Adjoint{<:Any, <:CuVecOrMat}) = generic_matmatmul!(C, A, B)
LinearAlgebra.mul!(C::CuVecOrMat, A::LinearAlgebra.Adjoint{<:Any, <:CuVecOrMat}, B::LinearAlgebra.Transpose{<:Any, <:CuVecOrMat}) = generic_matmatmul!(C, A, B)
LinearAlgebra.mul!(C::CuVecOrMat, A::LinearAlgebra.Adjoint{<:Any, <:CuVecOrMat}, B::LinearAlgebra.Adjoint{<:Any, <:CuVecOrMat}) = generic_matmatmul!(C, A, B)
LinearAlgebra.mul!(C::CuVecOrMat, A::LinearAlgebra.Transpose{<:Any, <:CuVecOrMat}, B::LinearAlgebra.Transpose{<:Any, <:CuVecOrMat}) = generic_matmatmul!(C, A, B)

function generic_rmul!(X::CuArray, s::Number)
    function kernel(X, s)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        @inbounds X[i] *= s
        return
    end
    @cuda blocks=length(X) kernel(X, s)
    X
end

LinearAlgebra.rmul!(A::CuArray, b::Number) = generic_rmul!(A, b)

function generic_lmul!(s::Number, X::CuArray)
    function kernel(s, X)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        @inbounds X[i] = s*X[i]
        return
    end
    @cuda blocks=length(X) kernel(s, X)
    X
end

LinearAlgebra.lmul!(a::Number, B::CuArray) = generic_lmul!(a, B)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Center()

Applies the center transform to all columns of the table.
The center transform of the column `x`, with mean `μ`,
is defined by `x .- μ`.
"""
struct Center <: Colwise end

assertions(::Type{Center}) = [assert_continuous]

isrevertible(::Type{Center}) = true

colcache(::Center, x) = mean(x)

colapply(::Center, x, μ)  = @. x - μ

colrevert(::Center, y, μ) = @. y + μ


text:
# Randomly convert an operator into another one (binary->binary;
# unary->unary)
function mutateOperator(tree::Node, options::Options)::Node
    if countOperators(tree) == 0
        return tree
    end
    node = randomNode(tree)
    while node.degree == 0
        node = randomNode(tree)
    end
    if node.degree == 1
        node.op = rand(1:options.nuna)
    else
        node.op = rand(1:options.nbin)
    end
    return tree
end

# Randomly perturb a constant
function mutateConstant(
        tree::Node, temperature::T,
        options::Options)::Node where {T<:Real}
    # T is between 0 and 1.

    if countConstants(tree) == 0
        return tree
    end
    node = randomNode(tree)
    while node.degree != 0 || node.constant == false
        node = randomNode(tree)
    end

    bottom = convert(T, 1//10)
    maxChange = options.perturbationFactor * temperature + convert(T, 1) + bottom
    factor = maxChange^Float32(rand())
    makeConstBigger = rand() > 0.5

    if makeConstBigger
        node.val *= convert(CONST_TYPE, factor)
    else
        node.val /= convert(CONST_TYPE, factor)
    end

    if rand() > options.probNegate
        node.val *= -1
    end

    return tree
end

# Add a random unary/binary operation to the end of a tree
function appendRandomOp(tree::Node, options::Options, nfeatures::Int)::Node
    node = randomNode(tree)
    while node.degree != 0
        node = randomNode(tree)
    end

    choice = rand()
    makeNewBinOp = choice < options.nbin/(options.nuna + options.nbin)

    if makeNewBinOp
        newnode = Node(
            rand(1:options.nbin),
            makeRandomLeaf(nfeatures),
            makeRandomLeaf(nfeatures)
        )
    else
        newnode = Node(
            rand(1:options.nuna),
            makeRandomLeaf(nfeatures)
        )
    end

    if newnode.degree == 2
        node.r = newnode.r
    end
    node.l = newnode.l
    node.op = newnode.op
    node.degree = newnode.degree
    node.val = newnode.val
    node.feature = newnode.feature
    node.constant = newnode.constant

    return tree
end

# Insert random node
function insertRandomOp(tree::Node, options::Options, nfeatures::Int)::Node
    node = randomNode(tree)
    choice = rand()
    makeNewBinOp = choice < options.nbin/(options.nuna + options.nbin)
    left = copyNode(node)

    if makeNewBinOp
        right = makeRandomLeaf(nfeatures)
        newnode = Node(
            rand(1:options.nbin),
            left,
            right
        )
    else
        newnode = Node(
            rand(1:options.nuna),
            left
        )
    end
    if newnode.degree == 2
        node.r = newnode.r
    end
    node.l = newnode.l
    node.op = newnode.op
    node.degree = newnode.degree
    node.val = newnode.val
    node.feature = newnode.feature
    node.constant = newnode.constant
    return tree
end

# Add random node to the top of a tree
function prependRandomOp(tree::Node, options::Options, nfeatures::Int)::Node
    node = tree
    choice = rand()
    makeNewBinOp = choice < options.nbin/(options.nuna + options.nbin)
    left = copyNode(tree)

    if makeNewBinOp
        right = makeRandomLeaf(nfeatures)
        newnode = Node(
            rand(1:options.nbin),
            left,
            right
        )
    else
        newnode = Node(
            rand(1:options.nuna),
            left
        )
    end
    if newnode.degree == 2
        node.r = newnode.r
    end
    node.l = newnode.l
    node.op = newnode.op
    node.degree = newnode.degree
    node.val = newnode.val
    node.feature = newnode.feature
    node.constant = newnode.constant
    return node
end

function makeRandomLeaf(nfeatures::Int)::Node
    if rand() > 0.5
        return Node(randn(CONST_TYPE))
    else
        return Node(rand(1:nfeatures))
    end
end

# Return a random node from the tree with parent
function randomNodeAndParent(tree::Node, parent::Union{Node, Nothing})::Tuple{Node, Union{Node, Nothing}}
    if tree.degree == 0
        return tree, parent
    end
    a = countNodes(tree)
    b = 0
    c = 0
    if tree.degree >= 1
        b = countNodes(tree.l)
    end
    if tree.degree == 2
        c = countNodes(tree.r)
    end

    i = rand(1:1+b+c)
    if i <= b
        return randomNodeAndParent(tree.l, tree)
    elseif i == b + 1
        return tree, parent
    end

    return randomNodeAndParent(tree.r, tree)
end

# Select a random node, and replace it an the subtree
# with a variable or constant
function deleteRandomOp(tree::Node, options::Options, nfeatures::Int)::Node
    node, parent = randomNodeAndParent(tree, nothing)
    isroot = (parent === nothing)

    if node.degree == 0
        # Replace with new constant
        newnode = makeRandomLeaf(nfeatures)
        node.degree = newnode.degree
        node.val = newnode.val
        node.constant = newnode.constant
        if !newnode.constant
            node.feature = newnode.feature
        end
    elseif node.degree == 1
        # Join one of the children with the parent
        if isroot
            return node.l
        elseif parent.l == node
            parent.l = node.l
        else
            parent.r = node.l
        end
    else
        # Join one of the children with the parent
        if rand() < 0.5
            if isroot
                return node.l
            elseif parent.l == node
                parent.l = node.l
            else
                parent.r = node.l
            end
        else
            if isroot
                return node.r
            elseif parent.l == node
                parent.l = node.r
            else
                parent.r = node.r
            end
        end
    end
    return tree
end

# Create a random equation by appending random operators
function genRandomTree(length::Int, options::Options, nfeatures::Int)::Node
    tree = Node(convert(CONST_TYPE, 1))
    for i=1:length
        tree = appendRandomOp(tree, options, nfeatures)
    end
    return tree
end


text:
int_rules_6_1_11 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 6.1.11*(e*x)^m*(a+b*x^n)^p*sinh =#
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ ~n) ^ ~(p') * sinh(~(c') + ~(d') * ~x), ~x) => Antiderivative(ExpandIntegrand(sinh(~c + ~d * ~x), (~a + ~b * (~x) ^ ~n) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ ~n) ^ ~(p') * cosh(~(c') + ~(d') * ~x), ~x) => Antiderivative(ExpandIntegrand(cosh(~c + ~d * ~x), (~a + ~b * (~x) ^ ~n) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ ~n) ^ ~p * sinh(~(c') + ~(d') * ~x), ~x) => (((~x) ^ (-(~n) + 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * sinh(~c + ~d * ~x)) / (~b * ~n * (~p + 1)) - ((-(~n) + 1) / (~b * ~n * (~p + 1))) * Antiderivative((~x) ^ -(~n) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * sinh(~c + ~d * ~x), ~x)) - (~d / (~b * ~n * (~p + 1))) * Antiderivative((~x) ^ (-(~n) + 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * cosh(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (IntegerQ(~p) && (IGtQ(~n, 0) && (LtQ(~p, -1) && GtQ(~n, 2))))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ ~n) ^ ~p * cosh(~(c') + ~(d') * ~x), ~x) => (((~x) ^ (-(~n) + 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * cosh(~c + ~d * ~x)) / (~b * ~n * (~p + 1)) - ((-(~n) + 1) / (~b * ~n * (~p + 1))) * Antiderivative((~x) ^ -(~n) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * cosh(~c + ~d * ~x), ~x)) - (~d / (~b * ~n * (~p + 1))) * Antiderivative((~x) ^ (-(~n) + 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * sinh(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (IntegerQ(~p) && (IGtQ(~n, 0) && (LtQ(~p, -1) && GtQ(~n, 2))))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ ~n) ^ ~p * sinh(~(c') + ~(d') * ~x), ~x) => Antiderivative(ExpandIntegrand(sinh(~c + ~d * ~x), (~a + ~b * (~x) ^ ~n) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (ILtQ(~p, 0) && (IGtQ(~n, 0) && (EqQ(~n, 2) || EqQ(~p, -1))))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ ~n) ^ ~p * cosh(~(c') + ~(d') * ~x), ~x) => Antiderivative(ExpandIntegrand(cosh(~c + ~d * ~x), (~a + ~b * (~x) ^ ~n) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (ILtQ(~p, 0) && (IGtQ(~n, 0) && (EqQ(~n, 2) || EqQ(~p, -1))))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ ~n) ^ ~p * sinh(~(c') + ~(d') * ~x), ~x) => Antiderivative((~x) ^ (~n * ~p) * (~b + ~a * (~x) ^ -(~n)) ^ ~p * sinh(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (ILtQ(~p, 0) && ILtQ(~n, 0))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ ~n) ^ ~p * cosh(~(c') + ~(d') * ~x), ~x) => Antiderivative((~x) ^ (~n * ~p) * (~b + ~a * (~x) ^ -(~n)) ^ ~p * cosh(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (ILtQ(~p, 0) && ILtQ(~n, 0))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ ~n) ^ ~p * sinh(~(c') + ~(d') * ~x), ~x) => Unintegrable((~a + ~b * (~x) ^ ~n) ^ ~p * sinh(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x)
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ ~n) ^ ~p * cosh(~(c') + ~(d') * ~x), ~x) => Unintegrable((~a + ~b * (~x) ^ ~n) ^ ~p * cosh(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * sinh(~(c') + ~(d') * ~x), ~x) => Antiderivative(ExpandIntegrand(sinh(~c + ~d * ~x), (~e * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * cosh(~(c') + ~(d') * ~x), ~x) => Antiderivative(ExpandIntegrand(cosh(~c + ~d * ~x), (~e * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~p * sinh(~(c') + ~(d') * ~x), ~x) => ((~e) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * sinh(~c + ~d * ~x)) / (~b * ~n * (~p + 1)) - ((~d * (~e) ^ ~m) / (~b * ~n * (~p + 1))) * Antiderivative((~a + ~b * (~x) ^ ~n) ^ (~p + 1) * cosh(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && (IntegerQ(~p) && (EqQ((~m - ~n) + 1, 0) && (LtQ(~p, -1) && (IntegerQ(~n) || GtQ(~e, 0)))))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~p * cosh(~(c') + ~(d') * ~x), ~x) => ((~e) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * cosh(~c + ~d * ~x)) / (~b * ~n * (~p + 1)) - ((~d * (~e) ^ ~m) / (~b * ~n * (~p + 1))) * Antiderivative((~a + ~b * (~x) ^ ~n) ^ (~p + 1) * sinh(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && (IntegerQ(~p) && (EqQ((~m - ~n) + 1, 0) && (LtQ(~p, -1) && (IntegerQ(~n) || GtQ(~e, 0)))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~p * sinh(~(c') + ~(d') * ~x), ~x) => (((~x) ^ ((~m - ~n) + 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * sinh(~c + ~d * ~x)) / (~b * ~n * (~p + 1)) - (((~m - ~n) + 1) / (~b * ~n * (~p + 1))) * Antiderivative((~x) ^ (~m - ~n) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * sinh(~c + ~d * ~x), ~x)) - (~d / (~b * ~n * (~p + 1))) * Antiderivative((~x) ^ ((~m - ~n) + 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * cosh(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (ILtQ(~p, -1) && (IGtQ(~n, 0) && (RationalQ(~m) && (GtQ((~m - ~n) + 1, 0) || GtQ(~n, 2)))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~p * cosh(~(c') + ~(d') * ~x), ~x) => (((~x) ^ ((~m - ~n) + 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * cosh(~c + ~d * ~x)) / (~b * ~n * (~p + 1)) - (((~m - ~n) + 1) / (~b * ~n * (~p + 1))) * Antiderivative((~x) ^ (~m - ~n) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * cosh(~c + ~d * ~x), ~x)) - (~d / (~b * ~n * (~p + 1))) * Antiderivative((~x) ^ ((~m - ~n) + 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * sinh(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (ILtQ(~p, -1) && (IGtQ(~n, 0) && (RationalQ(~m) && (GtQ((~m - ~n) + 1, 0) || GtQ(~n, 2)))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~p * sinh(~(c') + ~(d') * ~x), ~x) => Antiderivative(ExpandIntegrand(sinh(~c + ~d * ~x), (~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (ILtQ(~p, 0) && (IntegerQ(~m) && (IGtQ(~n, 0) && (EqQ(~n, 2) || EqQ(~p, -1)))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~p * cosh(~(c') + ~(d') * ~x), ~x) => Antiderivative(ExpandIntegrand(cosh(~c + ~d * ~x), (~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (ILtQ(~p, 0) && (IntegerQ(~m) && (IGtQ(~n, 0) && (EqQ(~n, 2) || EqQ(~p, -1)))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~p * sinh(~(c') + ~(d') * ~x), ~x) => Antiderivative((~x) ^ (~m + ~n * ~p) * (~b + ~a * (~x) ^ -(~n)) ^ ~p * sinh(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (ILtQ(~p, 0) && ILtQ(~n, 0))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~p * cosh(~(c') + ~(d') * ~x), ~x) => Antiderivative((~x) ^ (~m + ~n * ~p) * (~b + ~a * (~x) ^ -(~n)) ^ ~p * cosh(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (ILtQ(~p, 0) && ILtQ(~n, 0))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * sinh(~(c') + ~(d') * ~x), ~x) => Unintegrable((~e * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p * sinh(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n) ^ ~(p') * cosh(~(c') + ~(d') * ~x), ~x) => Unintegrable((~e * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p * cosh(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x)
end

max_stars_repo_path: src/rules/6 Hyperbolic functions/6.1 Hyperbolic sine/6.1.11 (e x)^m (a+b x^n)^p sinh.jl


text:
using FastGaussQuadrature

struct IntegratedLimbDark{LD<:AbstractLimbDark,WT,NT} <: AbstractLimbDark
    driver::LD
    order::Int
    nodes::NT
    weights::WT
end

@doc raw"""
    IntegratedLimbDark(limbdark; N=21, basis=:legendre)
    IntegratedLimbDark(u; kwargs...)

Computes the time-averaged flux in the middle of an exposure by wrapping a limb darkening law `limbdark` with a quadrature scheme. For each time step `t`, `N` extra points are *super-sampled* from `t-texp/2` to `t+texp/2`and the time-averaged flux is calculated via quadrature.

If a set of limb darkening coefficients, `u`, is provided, a [`PolynomialLimbDark`](@ref) law will be used by default.

# Mathematical form
```math
\bar{F}(t) = \frac{1}{\Delta t}\int_{t-\Delta t / 2}^{t+\Delta t / 2}{F(t')dt'}
```
where $F$ is the wrapped limb darkening law and $\Delta t$ is the exposure time.

## Quadrature
The integration is approximated via [Guassian quadrature](https://en.wikipedia.org/wiki/Gaussian_quadrature)
```math
\frac{1}{\Delta t} \int{F(t')dt'} \approx \frac12\sum_i^N{w_i * F(\frac{\Delta t}{2}\xi_i + t)}
```
where the weights `w_i` and nodes `ξ_i` are defined by the given quadrature rule. The nodes are defined by evaluating orthogonal polynomials `N` times between -1 and 1. Notice the change of interval required to go from the natural bounds of the orthogonal polynomial basis, `-1, 1`, to the range defined by the exposure time.

The following bases are available from [FastGaussQuadrature.jl](https://github.com/JuliaApproximation/FastGaussQuadrature.jl). In addition, a function can be passed which calculates `nodes, weights = f(N)`.
* `:legendre` - Legendre polynomial base on the open `(-1, 1)`
* `:radau` - Legendre polynomial base on the semi-open `[-1, 1)` interval
* `:lobatto` - Legendre polynomial base on the closed `[-1, 1]` interval
"""
function IntegratedLimbDark(driver::AbstractLimbDark; N=21, basis=:legendre)
    nodes, weights = get_nodes_weights(basis, N)
    return IntegratedLimbDark(driver, N, nodes, weights)
end

get_nodes_weights(basis, N) = basis(N)

function get_nodes_weights(basis::Symbol, N)
    if basis === :legendre
        return gausslegendre(N)
    elseif basis === :radau
        return gaussradau(N)
    elseif basis === :lobatto
        return gausslobatto(N)
    else
        error("$basis not a support quadrature basis")
    end
end

function IntegratedLimbDark(u::AbstractVector; kwargs...)
    driver = PolynomialLimbDark(u)
    return IntegratedLimbDark(driver; kwargs...)
end

compute(ld::IntegratedLimbDark, orbit::AbstractOrbit, t, r; texp=nothing) =
    compute(ld, orbit, t, r, texp)

function compute(ld::IntegratedLimbDark, orbit::AbstractOrbit, t, r, texp)
    # perform change of interval
    half_texp = 0.5 * texp
    # perform quadrature
    flux = sum(zip(ld.weights, ld.nodes)) do (w, ξ)
        f = compute(ld.driver, orbit, half_texp * ξ + t, r)
        return w * f
    end
    return 0.5 * flux
end

compute(ld::IntegratedLimbDark, orbit::AbstractOrbit, t, r, ::Nothing) =
    compute(ld.driver, orbit, t, r)


text:
function main(dataPath)
	lines = parse(readall(pipeline(`cat $[dataPath]ratings.csv`, `wc -l`)))
	@time ratings = readdlm(dataPath * "ratings.csv", ',', header=true, dims=(lines, 4))
	ratingsHeader = ratings[2]
	ratings = ratings[1]

	order = shuffle(collect(1:size(ratings)[1]))
	ratings = ratings[order, :]

	trainEnd = Int(floor(size(ratings)[1] * 0.8))
	train = ratings[1 : trainEnd, :]
	crossEnd = Int(floor(size(ratings)[1] * 0.9))
	cross = ratings[trainEnd+1 : crossEnd, :]
	test = ratings[crossEnd+1 : end, :]

	train = [ratingsHeader; train]
	cross = [ratingsHeader; cross]
	test = [ratingsHeader; test]

	@printf "Finished moving, begin saving\n"

	writecsv(dataPath * "train.csv", train)
	writecsv(dataPath * "cross.csv", cross)
	writecsv(dataPath * "test.csv", test)

	@printf "Done :)\n"
end

main("ml-latest-small/fixed/")


text:
using Test, Suppressor

# note, capturing test results does not work in a testset, this functionality is not tested
@testset "./script/utilities.jl" begin

include("../../script/utilities.jl")
using .ScriptUtilities: print_test_results, capture_test

test_pass() = @testset "test_method_pass()" begin @test true end
test_fail() = @testset "test_method_fail()" begin @test false end
test_args(x) = @testset "test_method_args()" begin @test x == 1 end
not_test_pass() = println("Hello World!")
not_test_fail() = throw(ErrorException("Hello World Failure!"))
not_test_args(x) = println(x)

@testset "print_test_results()" begin
    @testset "Boolean pass return values" begin
        pass = @suppress print_test_results(not_test_pass)
        @test pass

        pass = @suppress print_test_results(not_test_fail)
        @test !pass
    end

    @testset "test_pass() results print to STDOUT" begin
        output = @capture_out print_test_results(test_pass)

        @test occursin(r"Version :", output)
        @test occursin(r"(\[96aab1c2\])? SignalingDimension", output)
        @test occursin(r"Test Completed : \d\d\d\d-\d\d-\d\dT\d\d:\d\d:\d\d\n", output)
        @test occursin(r"Elapsed Time : \d+\.\d+\n", output)
        @test occursin(r"Test Method : test_pass\n", output)
        @test occursin(r"Test Args : \(\)\n", output)
        @test occursin(r"Test Pass : true\n", output)
    end

    @testset "test_args(1) results print to STDOUT" begin
        output = @capture_out print_test_results(test_args, params=[1])

        @test occursin(r"Version :", output)
        @test occursin(r"(\[96aab1c2\])? SignalingDimension", output)
        @test occursin(r"Test Completed : \d\d\d\d-\d\d-\d\dT\d\d:\d\d:\d\d\n", output)
        @test occursin(r"Elapsed Time : \d+\.\d+\n", output)
        @test occursin(r"Test Method : test_args\n", output)
        @test occursin(r"Test Args : \(1\)\n", output)
        @test occursin(r"Test Pass : true\n", output)
    end

    @testset "not_test_fail() results print to STDOUT" begin
        output = @capture_out print_test_results(not_test_fail)

        @test occursin(r"Version :", output)
        @test occursin(r"(\[96aab1c2\])? SignalingDimension", output)
        @test occursin(r"Test Completed : \d\d\d\d-\d\d-\d\dT\d\d:\d\d:\d\d\n", output)
        @test occursin(r"Elapsed Time : \d+\.\d+\n", output)
        @test occursin(r"Test Method : not_test_fail\n", output)
        @test occursin(r"Test Args : \(\)\n", output)
        @test occursin(r"Test Pass : false\n", output)
        @test occursin(r"Hello World Failure!", output)
    end

    @testset "prints to file" begin
        test_dir = mkdir("./test/unit/script_utilities_temp")

        print_test_results(test_pass, stdout=false, dir=test_dir)

        files = readdir(test_dir)

        @test length(files) == 1
        @test occursin(r"^test_pass_\d\d\d\d-\d\d-\d\dT\d\d-\d\d-\d\d\.txt$", files[1])

        output = read(test_dir * "/" * files[1], String)

        @test occursin(r"Version :", output)
        @test occursin(r"(\[96aab1c2\])? SignalingDimension", output)
        @test occursin(r"Test Completed : \d\d\d\d-\d\d-\d\dT\d\d:\d\d:\d\d\n", output)
        @test occursin(r"Elapsed Time : \d+\.\d+\n", output)
        @test occursin(r"Test Method : test_pass\n", output)
        @test occursin(r"Test Args : \(\)\n", output)
        @test occursin(r"Test Pass : true\n", output)

        rm(test_dir, recursive = true)
    end
end

@testset "capture_test()" begin
    @testset "capture stdout from pass" begin
        output = capture_test(not_test_pass)
        @test output == (true, "Hello World!\n")
    end

    @testset "args passed to test method" begin
        output = capture_test(not_test_args, params=[57])
        @test output == (true, "57\n")
    end

    @testset "capture exception from fail" begin
        output = capture_test(not_test_fail)
        @test output == (false, "ErrorException(\"Hello World Failure!\")\n")
    end
end

end


text:
using QSM: padfastfft, padarray!, unpadarray!

@testset "LSMR" begin
    include("lsmr.jl")
end

@testset "Poisson solver" begin
    include("poisson_solver/poisson_solver.jl")
end

@testset "padfastfft and psf2otf" for T in (Float32, Float64)
    Random.seed!(501)
    for sz in [(8, 9, 10), (11, 12, 13)]
        for ksz in [(3, 4, 3), (4, 3, 4)]
            x = randn(T, sz)
            k = randn(T, ksz)
            y = zeros(T, sz)
            z = zeros(T, sz)

            y = unpadarray!(y, conv(k, x))

            for _rfft in (false, true)
                xp = padfastfft(x, ksz, rfft=_rfft)
                K = psf2otf(k, size(xp), rfft=_rfft)

                if _rfft
                    @test iseven(size(xp, 1))
                    P = plan_rfft(xp)
                    zp = P\(K.*(P*xp))
                else
                    P = plan_fft(xp)
                    zp = real(P\(K.*(P*xp)))
                end

                @test unpadarray!(z, zp) ≈ y
            end
        end
    end
end

@testset "padarray!" begin
    A = reshape(1:25, 5, 5)

    @test padarray!(rand(Int, (8, 8)), A, :fill, 0) ==
        [
            0  0  0   0   0   0   0  0
            0  0  0   0   0   0   0  0
            0  0  1   6  11  16  21  0
            0  0  2   7  12  17  22  0
            0  0  3   8  13  18  23  0
            0  0  4   9  14  19  24  0
            0  0  5  10  15  20  25  0
            0  0  0   0   0   0   0  0
        ]

    @test padarray!(rand(Int, (8, 9)), A, :circular) ==
        [
            19  24  4   9  14  19  24  4   9
            20  25  5  10  15  20  25  5  10
            16  21  1   6  11  16  21  1   6
            17  22  2   7  12  17  22  2   7
            18  23  3   8  13  18  23  3   8
            19  24  4   9  14  19  24  4   9
            20  25  5  10  15  20  25  5  10
            16  21  1   6  11  16  21  1   6
        ]

    @test padarray!(rand(Int, (9, 8)), A, :replicate) ==
        [
            1  1  1   6  11  16  21  21
            1  1  1   6  11  16  21  21
            1  1  1   6  11  16  21  21
            2  2  2   7  12  17  22  22
            3  3  3   8  13  18  23  23
            4  4  4   9  14  19  24  24
            5  5  5  10  15  20  25  25
            5  5  5  10  15  20  25  25
            5  5  5  10  15  20  25  25
        ]

    @test padarray!(rand(Int, (9, 9)), A, :symmetric) ==
        [
             7  2  2   7  12  17  22  22  17
             6  1  1   6  11  16  21  21  16
             6  1  1   6  11  16  21  21  16
             7  2  2   7  12  17  22  22  17
             8  3  3   8  13  18  23  23  18
             9  4  4   9  14  19  24  24  19
            10  5  5  10  15  20  25  25  20
            10  5  5  10  15  20  25  25  20
             9  4  4   9  14  19  24  24  19
        ]

    @test padarray!(rand(Int, (9, 9)), A, :reflect) ==
        [
            13   8  3   8  13  18  23  18  13
            12   7  2   7  12  17  22  17  12
            11   6  1   6  11  16  21  16  11
            12   7  2   7  12  17  22  17  12
            13   8  3   8  13  18  23  18  13
            14   9  4   9  14  19  24  19  14
            15  10  5  10  15  20  25  20  15
            14   9  4   9  14  19  24  19  14
            13   8  3   8  13  18  23  18  13
        ]
end


text:
using CMDC

c = Client()

function multiple_dataset_all(c=Client())
    request!(c, demographics(), covid_us())
    fetch(c)
end

function multiple_dataset_states_only(c=Client())
    request!(c, demographics(), covid_us(location="<100"))
    fetch(c)
end

function multiple_dataset_counties_only(c=Client())
    request!(c, demographics(), covid_us(location=">1000"))
    fetch(c)
end

function single_dataset_all(c=Client())
    request!(c, mobility_devices())
    fetch(c)
end

function single_dataset_deaths_filter(c=Client())
    request!(c, covid_us(location="<100", variable="deaths_total", value=">100"))

    fetch(c)
end

function single_dataset_multiplestatesallcounties(c=Client())
    request!(c, mobility_devices(state=[6, 48]))

    fetch(c)
end

function single_dataset_onestateallcounties(c=Client())
    request!(c, mobility_devices(state=6))
    fetch(c)
end

function single_dataset_variableselect(c=Client())
    request!(
        c, demographics(
            variable=[
                "Total population",
                "Fraction of population over 65",
                "Median age"
            ]
        )
    )

    fetch(c)
end


text:
using DrWatson
@quickactivate "PECUZAL_Julia"

using DelayEmbeddings, DynamicalSystemsBase
using DifferentialEquations
using Random
using Test
using DelimitedFiles

include("../../src/pecuzal_method.jl")

## In this script we evaluate the dependence of the returned reconstruction
# parameters on the parameter determining the δ-neighborhood, needed for the
# computation of the continuity-statistic

## Dependence on δ-neighborhood

# The time series have been computed in the script `/test/timeseries/produce_timeseries.jl`
s = readdlm("./test/timeseries/lorenz_pecora_uni_x.csv")
s = vec(s[1:5000]) # input timeseries = x component of lorenz
w = estimate_delay(s, "mi_min")
Tmax = 100
KNN = 3
samplesize = 1
Tw = 20

Ks = 8:20
τ_vals = []
Ls = []
sizes = []
for K in Ks
    display(K)
    YY, τ_valss, _, Lss , _ = pecuzal_embedding_update(s;
                                τs = 0:Tmax , w = w, samplesize = samplesize,
                                K = K, KNN = KNN, Tw = Tw)
    push!(sizes,size(YY,2))
    push!(τ_vals, τ_valss)
    push!(Ls, Lss)
end

writedlm("./scripts/computed data/dependence_on_delta_neighborhood_Ks.csv",Ks)
writedlm("./scripts/computed data/dependence_on_delta_neighborhood_sizes.csv",sizes)
writedlm("./scripts/computed data/dependence_on_delta_neighborhood_Ls.csv",Ls)
writedlm("./scripts/computed data/dependence_on_delta_neighborhood_tau_vals.csv",τ_vals)

## Plot results (see script `plot_results_params_dependencies.jl`)

# deltas = readdlm("./scripts/computed data/dependence_on_delta_neighborhood_Ks.csv")
# sizes = readdlm("./scripts/computed data/dependence_on_delta_neighborhood_sizes.csv")
# Ls = readdlm("./scripts/computed data/dependence_on_delta_neighborhood_Ls.csv")
# τ_vals = readdlm("./scripts/computed data/dependence_on_delta_neighborhood_tau_vals.csv")
#
#
# using PyPlot
# pygui(true)
#
# fig = figure(figsize=[15,8])
#
# panelnames = ["A" "B"]
#
# lwg = 3         # linewidth of the graph
# lwa = 2         # linewidth of the axis
# fsp = 20        # Fontsize of panelnames
# fsa = 16        # Fontsize of the axis
# fsl = 12        # Fontsize of the legendentries
# fst = 16        # Fontsize of title
# axislabelsize = 16 # axislabelsize
# ticklabelsize = 12  # labelsize of ticks
# ms = 10         # markersize of maxima
# ms_style = "*"  # markerstyle
#
# subplots_adjust(left = 0.07,
# right = 0.96,
# bottom = 0.08,
# top = 0.9,
# wspace = 0.25,
# hspace = 0.25)
#
# rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
# font0 = Dict(
#         "font.size" => fsa,
#         "axes.labelweight" => "normal",
#         "axes.labelsize" => axislabelsize,
#         "axes.linewidth" => lwa,
#         "xtick.labelsize" => ticklabelsize,
#         "ytick.labelsize" => ticklabelsize,
#         "legend.fontsize" => fsl)
# merge!(rcParams, font0)
#
# ax = fig.add_subplot(1, 2, 1)
# p1 = plot(deltas, Ls[:,1], label="embedding cylce 1", linewidth=lwg)
# color1 = p1[1].get_color()
# scatter(deltas, Ls[:,1], label="")
# p2 = plot(deltas, Ls[:,2], label="embedding cylce 2", linewidth=lwg)
# color2 = p2[1].get_color()
# scatter(deltas, Ls[:,2], label="")
# legend()
# ax.text(-0.15, 1.1, panelnames[1], transform=ax.transAxes,
#      fontsize=fsp, fontweight="bold", va="top")
# title("minimum L")
# xlabel(L"\delta-Neighborhood size")
# ylabel(L"L_{min}")
# grid()
#
# ax = fig.add_subplot(1, 2, 2)
# plot(deltas, τ_vals[:,1], label="τ₁, embedding cylce 0", linewidth=lwg, color="r")
# scatter(deltas, τ_vals[:,1], label="", color="r")
# plot(deltas, τ_vals[:,2], label="τ₂, embedding cylce 1", linewidth=lwg, color=color1)
# scatter(deltas, τ_vals[:,2], label="", color=color1)
# plot(deltas, τ_vals[:,3], label="τ₃, embedding cylce 2", linewidth=lwg, color=color2)
# scatter(deltas, τ_vals[:,3], label="", color=color2)
# legend()
# ax.set_ylim([-1,24])
# ax.text(-0.15, 1.1, panelnames[2], transform=ax.transAxes,
#      fontsize=fsp, fontweight="bold", va="top")
# title("chosen delays")
# xlabel(L"\delta-Neighborhood size")
# ylabel("τ")
# grid()

max_stars_repo_path: scripts/dependence on input parameters/dependence_on_delta_neighborhood.jl


text:
using PyPlot

"""
Plot the mole fractions of the simulation bsol from t0 to tf
using N logarithmically spaced time points
only plots species who have mole fractions > tol at some point
in the simulation
"""
function plotmolefractions(bsol::Q, tf::V; t0::Z=1e-15,N::Z2=1000,tol::Z3=0.01,exclude::M=Array{String,1}()) where {Q<:Simulation, V<:Real, Z<:Real, Z2<:Real, Z3<:Real, M<:AbstractArray{String,1}}
    ts = exp.(range(log(t0),length=N,stop=log(tf)))
    xs = hcat(molefractions.(bsol,ts)...)
    maxes = maximum(xs,dims=2)
    spnames = []
    for i = 1:length(maxes)
        if maxes[i] > tol && !(bsol.domain.phase.species[i].name in exclude)
            plot(ts,xs[i,:])
            push!(spnames,bsol.domain.phase.species[i].name)
        end
    end
    legend(spnames)
    xlabel("Time in sec")
    ylabel("Mole Fraction")
end

"""
Plot the mole fractions of the simulation bsol at the time points solved for
only plots species who have mole fractions > tol at some point
in the simulation
"""
function plotmolefractions(bsol::Q; tol::V=0.01, exclude::M=Array{String,1}()) where {Q<:Simulation, V<:Real, M<:AbstractArray{String,1}}
    xs = molefractions(bsol)
    maxes = maximum(xs,dims=2)
    spnames = []
    for i = 1:length(maxes)
        if maxes[i] > tol && !(bsol.domain.phase.species[i].name in exclude)
            plot(bsol.sol.t,xs[i,:])
            push!(spnames,bsol.domain.phase.species[i].name)
        end
    end
    legend(spnames)
    xlabel("Time in sec")
    ylabel("Mole Fraction")
end

"""
Plot the molefractions of the species with names in spcnames over
the bsol time interval
"""
function plotmolefractions(bsol::Q,spcnames::V) where {Q<:Simulation,V<:AbstractArray}
    for name in spcnames
        plot(bsol.sol.t,[molefractions(bsol,name,t) for t in bsol.sol.t])
    end
    legend(spcnames)
end

export plotmolefractions

function plotradicalmolefraction(bsol;t=0.0,ts=[])
    if t == 0.0;
        t = bsol.sol.t[end]
    end
    if ts == []
        ts = range(0.0,stop=t, length=100);
    end
    radmf = zeros(length(ts))
    rads = [x.radicalelectrons for x in bsol.domain.phase.species]
    for (i,t) in enumerate(ts)
        radmf[i] = dot(rads,molefractions(bsol,t))
    end
    plot(ts,radmf)
    ylabel("Radical Mole Fraction")
    xlabel("Time in sec")
end

export plotradicalmolefraction

function plotmaxthermoforwardsensitivity(bsol, spcname; N=0, tol= 1e-2)
    spnames = getfield.(bsol.domain.phase.species,:name)
    values = Array{Float64,1}()
    outnames = Array{String,1}()
    for spn in spnames
        ind = argmax(abs.(getconcentrationsensitivity.(bsol,spcname,spn,bsol.sol.t)))
        val = getconcentrationsensitivity(bsol,spcname,spn,bsol.sol.t[ind])
        if abs(val)*4184.0 > tol
            push!(values,val)
            push!(outnames,spn)
        end
    end
    inds = reverse(sortperm(abs.(values)))
    if N == 0
        N = length(inds)
    elseif N > length(inds)
        N = length(inds)
    end
    inds = inds[1:N]
    xs = Array{Float64,1}(1:length(inds))
    barh(xs,reverse(values[inds].*4184.0))
    yticks(xs,reverse(outnames[inds]))
    xlabel("dLn([$spcname])/d(G_i) mol/kcal")
end
export plotmaxthermoforwardsensitivity

function plotmaxrateforwardsensitivity(bsol, spcname; N=0, tol= 1e-2)
    Nrxns = length(bsol.domain.phase.reactions)
    values = Array{Float64,1}()
    outinds = Array{Int64,1}()
    for i in 1:Nrxns
        ind = argmax(abs.(getconcentrationsensitivity.(bsol,spcname,i,bsol.sol.t)))
        val = getconcentrationsensitivity(bsol,spcname,i,bsol.sol.t[ind])
        if abs(val) > tol
            push!(values,val)
            push!(outinds,i)
        end
    end
    inds = reverse(sortperm(abs.(values)))
    if N == 0
        N = length(inds)
    elseif N > length(inds)
        N = length(inds)
    end
    inds = inds[1:N]
    xs = Array{Float64,1}(1:length(inds))
    barh(xs,reverse(values[inds]))
    yticks(xs,reverse(getrxnstr.(bsol.domain.phase.reactions[outinds[inds]])))
    xlabel("dLn([$spcname])/d(Ln(k_i))")
end
export plotmaxrateforwardsensitivity

"""
make a bar graph of the production/loss for the given species
associated with each reaction
N reactions are included all of which must have absolute value greater than abs(maximum prod or loss rate)*tol
"""
function plotrops(bsol::Y,name::X,t::Z;N=0,tol=0.01) where {Y<:Simulation, X<:AbstractString, Z<:Real}
    if !(name in getfield.(bsol.domain.phase.species,:name))
        error("Species $name not in domain")
    end
    rop = rops(bsol,name,t)
    inds = rop.nzind[reverse(sortperm(abs.(rop.nzval)))]
    if N == 0
        N = length(inds)
    elseif N > length(inds)
        N = length(inds)
    end
    inds = inds[1:N]
    mval = abs(rop[inds[1]])
    minval = mval*tol
    k = 1
    while k < length(inds) && abs(rop[inds[k]]) >= minval
        k += 1
    end
    inds = inds[1:k]
    xs = Array{Float64,1}(1:length(inds))
    barh(xs,reverse(rop[inds]))
    yticks(xs,reverse(getrxnstr.(bsol.domain.phase.reactions[inds])))
    xlabel("Production/Loss Rate mol/s")
    return
end

"""
make a bar graph of the production/loss for the given species
associated with each reaction
N reactions are included all of which must have absolute value greater than abs(maximum prod or loss rate)*tol
"""
function plotrops(ssys::Y,name::X,t::Z;N=0,tol=0.01) where {Y<:SystemSimulation, X<:AbstractString, Z<:Real}
    if !(name in ssys.names)
        error("Species $name not in domain")
    end
    rop = rops(ssys,name,t)
    inds = rop.nzind[reverse(sortperm(abs.(rop.nzval)))]
    if N == 0
        N = length(inds)
    elseif N > length(inds)
        N = length(inds)
    end
    inds = inds[1:N]
    mval = abs(rop[inds[1]])
    minval = mval*tol
    k = 1
    while k < length(inds) && abs(rop[inds[k]]) >= minval
        k += 1
    end
    inds = inds[1:k]
    xs = Array{Float64,1}(1:length(inds))
    barh(xs,reverse(rop[inds]))
    yticks(xs,reverse(getrxnstr.(ssys.reactions[inds])))
    xlabel("Production/Loss Rate mol/s")
    return
end

"""
make a line graph of the production/loss for the given species
associated with each reaction across a time domain
reactions with maximum (over time) production value greater than max production*tol or
maximum (over time) loss value greater than maximum loss*tol are included
"""
function plotrops(bsol::Y,name::X;rxnrates=Array{Float64,1}(),ts=Array{Float64,1}(),tol=0.05) where {Y<:Simulation, X<:AbstractString}
    if !(name in getfield.(bsol.domain.phase.species,:name))
        error("Species $name not in domain")
    end
    if length(rxnrates) == 0 || length(ts) == 0
        rxnrates = rates(bsol)
        ts = bsol.sol.t
    end
    ind = spcindex(bsol,name)

    cs = [count(isequal(ind),bsol.domain.phase.reactions[i].productinds)-count(isequal(ind),bsol.domain.phase.reactions[i].reactantinds) for i in 1:length(bsol.domain.phase.reactions)]
    @views maxrates = maximum(cs.*rxnrates,dims=2)[:,1]
    @views minrates = minimum(cs.*rxnrates,dims=2)[:,1]
    maxthresh = maximum(maxrates)*tol
    minthresh = minimum(minrates)*tol
    maxperm = reverse(sortperm(maxrates))
    minperm = sortperm(minrates)

    leg = Array{String,1}()
    for i in 1:length(maxperm)
        ind = maxperm[i]
        if maxrates[ind] >= maxthresh
            push!(leg,getrxnstr(bsol.domain.phase.reactions[ind]))
            plot(ts,rxnrates[ind,:]*cs[ind])
        else
            break
        end
    end
    for i in 1:length(minperm)
        ind = minperm[i]
        if minrates[ind] <= minthresh
            push!(leg,getrxnstr(bsol.domain.phase.reactions[ind]))
            plot(ts,rxnrates[ind,:]*cs[ind])
        else
            break
        end
    end
    legend(leg,loc="upper left", bbox_to_anchor=(1,1))
    ylabel("Flux in mol/s")
    xlabel("Time in sec")
end

"""
make a line graph of the production/loss for the given species
associated with each reaction across a time domain
reactions with maximum (over time) production value greater than max production*tol or
maximum (over time) loss value greater than maximum loss*tol are included
"""
function plotrops(bsol::Y,name::X;rxnrates=Array{Float64,1}(),ts=Array{Float64,1}(),tol=0.05) where {Y<:Simulation, X<:AbstractString}
    if !(name in getfield.(bsol.domain.phase.species,:name))
        error("Species $name not in domain")
    end
    if length(rxnrates) == 0 || length(ts) == 0
        rxnrates = rates(bsol)
        ts = bsol.sol.t
    end
    ind = spcindex(bsol,name)

    cs = [count(isequal(ind),bsol.domain.phase.reactions[i].productinds)-count(isequal(ind),bsol.domain.phase.reactions[i].reactantinds) for i in 1:length(bsol.domain.phase.reactions)]
    @views maxrates = maximum(cs.*rxnrates,dims=2)[:,1]
    @views minrates = minimum(cs.*rxnrates,dims=2)[:,1]
    maxthresh = maximum(maxrates)*tol
    minthresh = minimum(minrates)*tol
    maxperm = reverse(sortperm(maxrates))
    minperm = sortperm(minrates)

    leg = Array{String,1}()
    for i in 1:length(maxperm)
        ind = maxperm[i]
        if maxrates[ind] >= maxthresh
            push!(leg,getrxnstr(bsol.domain.phase.reactions[ind]))
            plot(ts,rxnrates[ind,:]*cs[ind])
        else
            break
        end
    end
    for i in 1:length(minperm)
        ind = minperm[i]
        if minrates[ind] <= minthresh
            push!(leg,getrxnstr(bsol.domain.phase.reactions[ind]))
            plot(ts,rxnrates[ind,:]*cs[ind])
        else
            break
        end
    end
    legend(leg,loc="upper left", bbox_to_anchor=(1,1))
    ylabel("Flux in mol/s")
    xlabel("Time in sec")
end

"""
make a line graph of the production/loss for the given species
associated with each reaction across a time domain
reactions with maximum (over time) production value greater than max production*tol or
maximum (over time) loss value greater than maximum loss*tol are included
"""
function plotrops(bsol::Y,name::X;rxnrates=Array{Float64,1}(),ts=Array{Float64,1}(),tol=0.05) where {Y<:Simulation, X<:AbstractString}
    if !(name in getfield.(bsol.domain.phase.species,:name))
        error("Species $name not in domain")
    end
    if length(rxnrates) == 0 || length(ts) == 0
        rxnrates = rates(bsol)
        ts = bsol.sol.t
    end
    ind = spcindex(bsol,name)

    cs = [count(isequal(ind),bsol.domain.phase.reactions[i].productinds)-count(isequal(ind),bsol.domain.phase.reactions[i].reactantinds) for i in 1:length(bsol.domain.phase.reactions)]
    @views maxrates = maximum(cs.*rxnrates,dims=2)[:,1]
    @views minrates = minimum(cs.*rxnrates,dims=2)[:,1]
    maxthresh = maximum(maxrates)*tol
    minthresh = minimum(minrates)*tol
    maxperm = reverse(sortperm(maxrates))
    minperm = sortperm(minrates)

    leg = Array{String,1}()
    for i in 1:length(maxperm)
        ind = maxperm[i]
        if maxrates[ind] >= maxthresh
            push!(leg,getrxnstr(bsol.domain.phase.reactions[ind]))
            plot(ts,rxnrates[ind,:]*cs[ind])
        else
            break
        end
    end
    for i in 1:length(minperm)
        ind = minperm[i]
        if minrates[ind] <= minthresh
            push!(leg,getrxnstr(bsol.domain.phase.reactions[ind]))
            plot(ts,rxnrates[ind,:]*cs[ind])
        else
            break
        end
    end
    legend(leg,loc="upper left", bbox_to_anchor=(1,1))
    ylabel("Flux in mol/s")
    xlabel("Time in sec")
end

export plotrops

"""
make a bar graph of the production/loss for all radicals
associated with each reaction
N reactions are included all of which must have absolute value greater than abs(maximum prod or loss rate)*tol
"""
function plotradicalrops(bsol::Y,t::Z;N=0,tol=0.01) where {Y<:Simulation, Z<:Real}
    rop = rates(bsol,t).*getfield.(bsol.domain.phase.reactions,:radicalchange)
    inds = reverse(sortperm(abs.(rop)))
    if N == 0
        N = length(inds)
    elseif N > length(inds)
        N = length(inds)
    end
    inds = inds[1:N]
    mval = abs(rop[inds[1]])
    minval = mval*tol
    k = 1
    while k < length(inds) && abs(rop[inds[k]]) >= minval
        k += 1
    end
    inds = inds[1:k]
    xs = Array{Float64,1}(1:length(inds))
    barh(xs,reverse(rop[inds]))
    yticks(xs,reverse(getrxnstr.(bsol.domain.phase.reactions[inds])))
    xlabel("Production/Loss Rate mol/(m^3*s)")
    return
end

"""
make a line graph of the production/loss for all radicals
associated with each reaction across a time domain
reactions with maximum (over time) production value greater than max production*tol or
maximum (over time) loss value greater than maximum loss*tol are included
"""
function plotradicalrops(bsol::Y;rxnrates=Array{Float64,1}(),ts=Array{Float64,1}(),tol=0.05) where {Y<:Simulation, X<:AbstractString}
    if length(rxnrates) == 0 || length(ts) == 0
        rxnrates = rates(bsol)
        ts = bsol.sol.t
    end

    cs = getfield.(bsol.domain.phase.reactions,:radicalchange)
    @views maxrates = maximum(cs.*rxnrates,dims=2)[:,1]
    @views minrates = minimum(cs.*rxnrates,dims=2)[:,1]
    maxthresh = maximum(maxrates)*tol
    minthresh = minimum(minrates)*tol
    maxperm = reverse(sortperm(maxrates))
    minperm = sortperm(minrates)

    leg = Array{String,1}()
    for i in 1:length(maxperm)
        ind = maxperm[i]
        if maxrates[ind] >= maxthresh
            push!(leg,getrxnstr(bsol.domain.phase.reactions[ind]))
            plot(ts,rxnrates[ind,:]*cs[ind])
        else
            break
        end
    end
    for i in 1:length(minperm)
        ind = minperm[i]
        if minrates[ind] <= minthresh
            push!(leg,getrxnstr(bsol.domain.phase.reactions[ind]))
            plot(ts,rxnrates[ind,:]*cs[ind])
        else
            break
        end
    end
    legend(leg,loc="upper left", bbox_to_anchor=(1,1))
    ylabel("Radical Flux in mol/s")
    xlabel("Time in sec")
end

export plotradicalrops

function plotthermoadjointsensitivities(bsol::Y,name::X,dps::Z;N=0,tol=0.01) where {Y<:Simulation, X<:AbstractString, Z}
    t = bsol.sol.t[end]
    if name in ["T","V"] || name in bsol.names
        dpvals = dps[1:length(bsol.domain.phase.species)].*4184.0
    else
        error("Name $name not in domain")
    end
    inds = reverse(sortperm(abs.(dpvals)))
    if N == 0
        N = length(inds)
    elseif N > length(inds)
        N = length(inds)
    end
    inds = inds[1:N]
    mval = abs(dpvals[inds[1]])
    minval = mval*tol
    k = 1
    while k < length(inds) && abs(dpvals[inds[k]]) >= minval
        k += 1
    end
    inds = inds[1:k]
    xs = Array{Float64,1}(1:length(inds))
    barh(xs,reverse(dpvals[inds]))
    yticks(xs,reverse(getfield.(bsol.domain.phase.species[inds],:name)))
    if name in ["T","V"]
        xlabel("d$name/dG mol/kcal")
    else
        xlabel("dLn(N$name)/dG mol/kcal")
    end
    return
end
export plotthermoadjointsensitivities

function plotrateadjointsensitivities(bsol::Y,name::X,dps::Z;N=0,tol=0.01) where {Y<:Simulation, X<:AbstractString, Z}
    if !(name in ["T","V"] || name in bsol.names)
        error("Name $name not in domain")
    end
    dpvals = dps[length(bsol.domain.phase.species)+1:end]
    inds = reverse(sortperm(abs.(dpvals)))
    if N == 0
        N = length(inds)
    elseif N > length(inds)
        N = length(inds)
    end
    inds = inds[1:N]
    mval = abs(dpvals[inds[1]])
    minval = mval*tol
    k = 1
    while k < length(inds) && abs(dpvals[inds[k]]) >= minval
        k += 1
    end
    inds = inds[1:k]
    xs = Array{Float64,1}(1:length(inds))
    barh(xs,reverse(dpvals[inds]))
    yticks(xs,reverse(getrxnstr.(bsol.domain.phase.reactions[inds])))

    if name in ["T","V"]
        xlabel("d$name/dLn(kf)")
    else
        xlabel("dLn(N$name)/dLn(kf)")
    end
    return
end
export plotrateadjointsensitivities

function plottimescales(sim,t;taumax=1e6,taumin=1e-18,taures=10.0^0.5,usediag=true)
    Jy = jacobiany(sim.sol(t),sim.domain.p,t,sim.domain,[]);
    return plottimescales(Jy;taumax=taumax,taumin=taumin,taures=taures,usediag=usediag)
end

function plottimescales(Jy;taumax=1e6,taumin=1e-18,taures=10.0^0.5,usediag=true)
    if usediag
        taus = 1.0./abs.(diag(Jy))
    else
        taus = 1.0./abs.(eigvals(Jy))
    end
    PyPlot.hist([x==Inf ? 0.0 : x for x in taus],bins=10.0.^(log10(taumin):log10(taures):log10(taumax)))
    xscale("log")
    xlabel("Species Time Scale [sec]")
    ylabel("Counts")
end
export plottimescales

function plotrxntransitorysensitivities(bsol,name,t;dSdt=nothing,tau=nothing,tol=1e-3,N=0,rxntol=1e-6)
    if !(name in getfield.(bsol.domain.phase.species,:name))
        error("Species $name not in domain")
    elseif !isnothing(dSdt) && (sum(dim > 1 for dim in size(dSdt)) > 1 || maximum(size(dSdt)) != length(bsol.p))
        error("dSdt must be a vector of length number of parameters")
    end

    ind = findfirst(isequal(name),bsol.names)

    rts = rates(bsol,t)
    Rchar = norm(rts)
    Rthresh = rxntol*Rchar

    if dSdt === nothing
        if tau === nothing
            dSdt = transitorysensitivitiesfulltrapezoidal(bsol,t)[ind,length(bsol.names)+1:end]
        else
            dSdt = transitorysensitivitiesfulltrapezoidal(bsol,t,tau)[ind,length(bsol.names)+1:end]
        end
    else
        dSdt = dSdt[length(bsol.names)+1:end]
    end

    inds = reverse(sortperm(abs.(dSdt)))
    minval = 0.0
    dSdtmax = maximum(abs.(dSdt))
    maxthresh = dSdtmax*tol

    inds = [i for i in inds if abs(rts[i]) > Rthresh] #weak filter based on reaction flux

    if N == 0
        N = length(inds)
    elseif N > length(inds)
        N = length(inds)
    end
    inds = inds[1:N]
    mval = abs(dSdt[inds[1]])
    minval = mval*tol
    k = 1
    while k < length(inds) && abs(dSdt[inds[k]]) >= minval
        k += 1
    end
    inds = inds[1:k]
    xs = Array{Float64,1}(1:length(inds))
    barh(xs,reverse(dSdt[inds]))
    yticks(xs,reverse(getrxnstr.(bsol.domain.phase.reactions[inds])))
    xlabel("d/dt dLn([$name])/d(Ln(k_i))")
    xscale("symlog")
    return
end
export plotrxntransitorysensitivities

function plotthermotransitorysensitivities(bsol,name,t;dSdt=nothing,tau=nothing,tol=1e-3,N=0)
    if !(name in getfield.(bsol.domain.phase.species,:name))
        error("Species $name not in domain")
    end
    ind = findfirst(isequal(name),bsol.names)

    if dSdt === nothing
        if tau === nothing
            dSdt = transitorysensitivitiesfulltrapezoidal(bsol,t)[ind,1:length(bsol.names)]
        else
            dSdt = transitorysensitivitiesfulltrapezoidal(bsol,t,tau)[ind,1:length(bsol.names)]
        end
    else
        dSdt = dSdt[1:length(bsol.names)]
    end

    inds = reverse(sortperm(abs.(dSdt)))
    dSdtmax = maximum(abs.(dSdt))
    maxthresh = dSdtmax*tol

    if N == 0
        N = length(inds)
    elseif N > length(inds)
        N = length(inds)
    end
    inds = inds[1:N]
    mval = abs(dSdt[inds[1]])
    minval = mval*tol
    k = 1
    while k < length(inds) && abs(dSdt[inds[k]]) >= minval
        k += 1
    end
    inds = inds[1:k]
    xs = Array{Float64,1}(1:length(inds))
    barh(xs,reverse(dSdt[inds]))
    yticks(xs,reverse(getfield.(bsol.domain.phase.species[inds],:name)))
    xlabel("Normalized Transitory Sensitivities for $name")
    return
end
export plotthermotransitorysensitivities


text:
## runtests.jl  Unit tests for MFCC
## (c) 2015 David A. van Leeuwen
##
## Licensed under the MIT software license, see LICENSE.md

using WAV
include("../src/MFCC.jl")
using .MFCC
using SpecialFunctions
using Statistics

# x, meta, params = MFCC.feacalc("../test/bl2.wav", normtype=:none, method=:wav, augtype=:none, sadtype=:none)
# y = MFCC.feaload("bl2.mfcc")

# @assert x == y

x = randn(1_000)

@show MFCC.mfcc(x)

println("Tests passed")


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

"""
    glasscatalogs()

Returns the complete list of glass catalogs available from GlassCat.

## Example
```julia-repl
julia> glasscatalogs()
41-element Array{Any,1}:
 OpticSim.GlassCat.AMTIR
 OpticSim.GlassCat.ANGSTROMLINK
 OpticSim.GlassCat.APEL
 OpticSim.GlassCat.ARCHER
 OpticSim.GlassCat.ARTON
 OpticSim.GlassCat.AUER_LIGHTING
 OpticSim.GlassCat.BIREFRINGENT
 ⋮
```
"""
glasscatalogs() = _child_modules(GlassCat)

"""
    glassnames(catalog::Module)

Returns the glass names available from a given catalog.

# Example
```julia-repl
julia> glassnames(GlassCat.CARGILLE)
3-element Array{Any,1}:
 "OG0607"
 "OG0608"
 "OG081160"
```
"""
function glassnames(catalog::Module)
    glass_names = names(catalog, all = true, imported = false)
    glasses = []
    for glass_name in glass_names
        glass_name_str = string(glass_name)
        if !occursin("#", glass_name_str) && glass_name_str != "eval" && glass_name_str != "include" && glass_name != nameof(catalog)
            push!(glasses, glass_name)
        end
    end
    return glasses
end

"""
    glassnames()

Returns the glass names available from all catalogs.

# Example
```julia-repl
julia> glassnames()
6-element Array{Pair{Module,Array{Any,1}},1}:
 OpticSim.GlassCat.CARGILLE => ["OG0607", "OG0608", "OG081160"]
     OpticSim.GlassCat.HOYA => ["BAC4", "BACD11"  …  "TAFD65"]
    OpticSim.GlassCat.NIKON => ["BAF10", "BAF11"  …  "_7054"]
    OpticSim.GlassCat.OHARA => ["L_BAL35", "L_BAL35P"  …  "S_TIM8"]
   OpticSim.GlassCat.SCHOTT => ["AF32ECO", "BAFN6"  …  "SFL6"]
   OpticSim.GlassCat.SUMITA => ["BAF1", "BAF10"  …  "ZNSF8"]
```
"""
glassnames() = [m => glassnames(m) for m in _child_modules(GlassCat)]

"""
    findglass(condition::Function) -> Vector{Glass}

Returns the list of glasses which satisfy `condition` where `condition::(Glass -> Bool)`.

# Example
```julia-repl
julia> findglass(x -> (x.Nd > 2.3 && x.λmin < 0.5 && x.λmax > 0.9))
8-element Array{GlassCat.Glass,1}:
 BIREFRINGENT.TEO2_E
 BIREFRINGENT.PBMOO4
 BIREFRINGENT.LINBO3
 INFRARED.CLEARTRAN_OLD
 INFRARED.CLEARTRAN
 INFRARED.SRTIO3
 INFRARED.ZNS_BROAD
 INFRARED.ZNS_VIS
```
"""
function findglass(condition::Function)
    # TODO - make the condition easier to specify (accessor functions for fields?)
    out = Vector{Glass}(undef, 0)
    for g in AGF_GLASSES
        if condition(g)
            push!(out, g)
        end
    end
    for g in OTHER_GLASSES
        if condition(g)
            push!(out, g)
        end
    end
    return out
end

function _child_modules(m::Module)
    ns = names(m, imported = false, all = true)
    ms = []
    for n in ns
        if n != nameof(m)
            try
                x = Core.eval(m, n)
                if x isa Module
                    push!(ms, x)
                end
            catch
            end
        end
    end
    return ms
end


text:
include("Include.jl")

# what training index?
index = 10

# go -
(p, Yₘ, Y) = learn_optim(index; pₒ = nothing)

max_stars_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation

max_issues_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation

max_forks_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation


text:
using Revise
using Soss

m = @model begin
    n = 10 # number of samples
    m = 10 # reviewers for each example

    p_bad ~ Beta(1,3) |> iid(n)

    fpr ~ Beta(2,5) |> iid(m)
    tpr ~ Beta(5,2) |> iid(m)

    y ~ For(n, m) do nj,mj
            Mix([Bernoulli(fpr[mj]), Bernoulli(tpr[mj])]
                , [1 - p_bad[nj], p_bad[nj]])
        end
end;

using Random
truth = rand(m());
symlogpdf(m())

@time logpdf(m(),truth, codegen)
@time logpdf(m(),truth, codegen)

codegen(m(),truth)

logpdf(m(), merge(truth, (p_bad=shuffle(truth.p_bad),)))

@time result = dynamicHMC(m(), (y=truth.y,), codegen) ;

# result = @time advancedHMC(m(), (y=truth.y,))

pairs(truth)
result |> particles |> pairs


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

module Terminals

export
    TextTerminal,
    UnixTerminal,
    TerminalBuffer,
    TTYTerminal,
    cmove,
    cmove_col,
    cmove_down,
    cmove_left,
    cmove_line_down,
    cmove_line_up,
    cmove_right,
    cmove_up,
    disable_bracketed_paste,
    enable_bracketed_paste,
    end_keypad_transmit_mode,
    getX,
    getY,
    hascolor,
    pos,
    raw!

import Base:
    check_open, # stream.jl
    displaysize,
    flush,
    pipe_reader,
    pipe_writer,
    read,
    readuntil

## TextTerminal ##

abstract TextTerminal <: Base.AbstractPipe

# INTERFACE
pipe_reader(::TextTerminal) = error("Unimplemented")
pipe_writer(::TextTerminal) = error("Unimplemented")
displaysize(::TextTerminal) = error("Unimplemented")
cmove(t::TextTerminal, x, y) = error("Unimplemented")
getX(t::TextTerminal) = error("Unimplemented")
getY(t::TextTerminal) = error("Unimplemented")
pos(t::TextTerminal) = (getX(t), getY(t))

# Relative moves (Absolute position fallbacks)
cmove_up(t::TextTerminal, n) = cmove(getX(t), max(1, getY(t)-n))
cmove_up(t) = cmove_up(t, 1)

cmove_down(t::TextTerminal, n) = cmove(getX(t), max(height(t), getY(t)+n))
cmove_down(t) = cmove_down(t, 1)

cmove_left(t::TextTerminal, n) = cmove(max(1, getX(t)-n), getY(t))
cmove_left(t) = cmove_left(t, 1)

cmove_right(t::TextTerminal, n) = cmove(max(width(t), getX(t)+n), getY(t))
cmove_right(t) = cmove_right(t, 1)

cmove_line_up(t::TextTerminal, n) = cmove(1, max(1, getY(t)-n))
cmove_line_up(t) = cmove_line_up(t, 1)

cmove_line_down(t::TextTerminal, n) = cmove(1, max(height(t), getY(t)+n))
cmove_line_down(t) = cmove_line_down(t, 1)

cmove_col(t::TextTerminal, c) = cmove(c, getY(t))

# Defaults
hascolor(::TextTerminal) = false

# Utility Functions
width(t::TextTerminal) = displaysize(t)[2]
height(t::TextTerminal) = displaysize(t)[1]

# For terminals with buffers
flush(t::TextTerminal) = nothing

clear(t::TextTerminal) = error("Unimplemented")
clear_line(t::TextTerminal, row) = error("Unimplemented")
clear_line(t::TextTerminal) = error("Unimplemented")

raw!(t::TextTerminal, raw::Bool) = error("Unimplemented")

beep(t::TextTerminal) = nothing
enable_bracketed_paste(t::TextTerminal) = nothing
disable_bracketed_paste(t::TextTerminal) = nothing

## UnixTerminal ##

abstract UnixTerminal <: TextTerminal

pipe_reader(t::UnixTerminal) = t.in_stream
pipe_writer(t::UnixTerminal) = t.out_stream

type TerminalBuffer <: UnixTerminal
    out_stream::Base.IO
end

type TTYTerminal <: UnixTerminal
    term_type::String
    in_stream::Base.TTY
    out_stream::Base.TTY
    err_stream::Base.TTY
end

const CSI = "\x1b["

cmove_up(t::UnixTerminal, n) = write(t.out_stream, "$(CSI)$(n)A")
cmove_down(t::UnixTerminal, n) = write(t.out_stream, "$(CSI)$(n)B")
cmove_right(t::UnixTerminal, n) = write(t.out_stream, "$(CSI)$(n)C")
cmove_left(t::UnixTerminal, n) = write(t.out_stream, "$(CSI)$(n)D")
cmove_line_up(t::UnixTerminal, n) = (cmove_up(t, n); cmove_col(t, 1))
cmove_line_down(t::UnixTerminal, n) = (cmove_down(t, n); cmove_col(t, 1))
cmove_col(t::UnixTerminal, n) = write(t.out_stream, "$(CSI)$(n)G")

if is_windows()
    function raw!(t::TTYTerminal,raw::Bool)
        check_open(t.in_stream)
        if Base.ispty(t.in_stream)
            run(if raw
                    `stty raw -echo onlcr -ocrnl opost`
                else
                    `stty sane`
                end,t.in_stream,t.out_stream,t.err_stream)
            true
        else
            ccall(:jl_tty_set_mode,
                 Int32, (Ptr{Void},Int32),
                 t.in_stream.handle, raw) != -1
        end
    end
else
    function raw!(t::TTYTerminal, raw::Bool)
        check_open(t.in_stream)
        ccall(:jl_tty_set_mode, Int32, (Ptr{Void},Int32), t.in_stream.handle, raw) != -1
    end
end
enable_bracketed_paste(t::UnixTerminal) = write(t.out_stream, "$(CSI)?2004h")
disable_bracketed_paste(t::UnixTerminal) = write(t.out_stream, "$(CSI)?2004l")
end_keypad_transmit_mode(t::UnixTerminal) = # tput rmkx
    write(t.out_stream, "$(CSI)?1l\x1b>")

function Base.displaysize(t::UnixTerminal)
    return displaysize(t.out_stream)
end

clear(t::UnixTerminal) = write(t.out_stream, "\x1b[H\x1b[2J")
clear_line(t::UnixTerminal) = write(t.out_stream, "\x1b[0G\x1b[0K")
#beep(t::UnixTerminal) = write(t.err_stream,"\x7")

if is_windows()
    hascolor(t::TTYTerminal) = true
else
    function hascolor(t::TTYTerminal)
        startswith(t.term_type, "xterm") && return true
        try
            return success(`tput setaf 0`)
        catch
            return false
        end
    end
end

end # module


text:
# Inspired from MathOptInterface/src/Test/config.jl
"""
    @polytestset setname subsets

Defines a function `setnametest(lib, exclude)` that runs the tests defined in the dictionary `setnametests`
with the library `lib` except the tests whose dictionary key is in `exclude`.
If `subsets` is `true` then each test runs in fact multiple tests hence the `exclude` argument is passed
as it can also contains test to be excluded from these subsets of tests.
"""
macro polytestset(setname, subsets=false)
    testname = Symbol(string(setname) * "test")
    testdict = Symbol(string(testname) * "s")
    if subsets
        runtest = :( f(lib, exclude) )
    else
        runtest = :( f(lib) )
    end
    esc(:(
        function $testname(lib::Polyhedra.Library, exclude::Vector{String} = String[])
            for (name,f) in $testdict
                if name in exclude
                    continue
                end
                @testset "$name" begin
                    $runtest
                end
            end
        end
    ))
end


text:
module DirectTrajectoryOptimization

using LinearAlgebra
using SparseArrays
using Symbolics, IfElse
using Parameters
using Ipopt
using MathOptInterface
const MOI = MathOptInterface

include(joinpath(pwd(), "src/objective.jl"))
include(joinpath(pwd(), "src/constraints.jl"))
include(joinpath(pwd(), "src/dynamics.jl"))
include(joinpath(pwd(), "src/problem.jl"))
include(joinpath(pwd(), "src/moi.jl"))
include(joinpath(pwd(), "src/solver.jl"))
include(joinpath(pwd(), "src/utils.jl"))

# objective
export Cost

# constraints
export Bound, StageConstraint, ConstraintSet

# dynamics
export Dynamics, DynamicsModel

# problem
export TrajectoryOptimizationProblem

# solver
export Solver, initialize!, solve!

# utils
export linear_interpolation

end # module


text:
export Cell

mutable struct Cell
    color::Float64
    obj::Float64
end

function Cell(color::String, obj::String)
    Cell(color_to_float(color), type_to_float(obj))
end

WallCell = Cell("black", "wall")

function Cell()
    Cell("white", "empty")
end

function object(c::Cell)
    float_to_type(c.obj)
end

function color(c::Cell)
    float_to_color(c.color)
end


text:
using RegisterQD
using LinearAlgebra
using RegisterQD.StaticArrays
using RegisterQD.Interpolations
using RegisterQD.Images
using RegisterQD.CoordinateTransformations
using RegisterQD.Rotations

g = 0.2:0.2:1.2
gradcube = g .* reshape(g, 1, 6) .* reshape(g, 1, 1, 6)
gradcube[1,4,2] = 0
gradcube[4,5,2] = 0 #break up any potential rotational symmetry

mxshift = (1,1,1);
mxrot = (0.01,0.01,0.01);

EYE =Matrix(1.0*I, 3,3)

@testset "initial_tfm improves rigid rotational alignment" begin
    testimage1 = zeros(Float64,10,10,10)
    testimage1 .= NaN
    testimage1[3:8,3:8,3:8] = gradcube
    testimage1 = centered(testimage1)

    tform = RotXYZ(0.1,0.1,0.1)
    mytform = AffineMap(tform, [0,0,0])
    testimage2 = warp(testimage1, mytform, axes(testimage1))

    # insufficient parameters + no initiat_tfm causes large misalignment

    tformtest1, mm1 = qd_rigid(testimage2, testimage1, mxshift, mxrot; print_interval=typemax(Int)) #
    @test !(mm1 < 1e-8)
    @test !isapprox(tformtest1, mytform, atol = 0.1)

    #initial_tfm improves alignment

    tformtest2, mm2 = qd_rigid(testimage2, testimage1, mxshift, mxrot; print_interval=typemax(Int), initial_tfm = mytform) #
 #with the initial_tfm being the true tfm, this should give back the true rotation

    @test mm2 < 1e-8
    @test isapprox(tformtest2, mytform, atol = 0.0001)
    @test isrotation(tformtest2.linear)
end #Test initial_tfm improves rotational alignment for rigid

@testset "initial_tfm improves rigid translational alignment" begin
    testimage3 = zeros(10,20,10)
    testimage3 .= NaN
    testimage3[3:8,3:8,3:8] = gradcube
    testimage3 = centered(testimage3)

    testimage4 = zeros(10,20,10)
    testimage4 .= NaN
    testimage4[3:8,13:18,3:8] = gradcube
    testimage4 = centered(testimage4)

    init_tfm = AffineMap(diagm(0 => fill(1.0,3)), [0.0,10.0,0.0])

    tformtest3 = nothing
    mm3 = nothing

    #no initial_tfm should fail due to total lack of overlap
    @test try
        tformtest3, mm3 = qd_rigid(testimage3, testimage4, mxshift, mxrot; print_interval=typemax(Int)) #
        false
    catch err
        true  # this should break
    end

    tformtest4, mm4 = qd_rigid(testimage3, testimage4, mxshift, mxrot; print_interval=typemax(Int), initial_tfm = init_tfm) #

    @test mm4 < 1e-8 #this mismatch should be lower!
    @test isapprox(tformtest4.translation, [0, 10, 0])
    @test isapprox(tformtest4.linear, EYE)
end

@testset "Test initial_tfm creates real rotations." begin #some of this seems redundant
    testimage1 = zeros(Float64,10,10,10)
    testimage1 .= NaN
    testimage1[3:8,3:8,3:8] = gradcube
    testimage1 = centered(testimage1)

    tform = RotXYZ(0.1,0.1,0.1)
    mytform = AffineMap(tform, [0,0,0])
    testimage2 = warp(testimage1, mytform, axes(testimage1))

    mxrot2 = (0.2,0.2,0.2);
    minwidth_rot = RegisterQD.default_minwidth_rot(CartesianIndices(testimage2), EYE)

    # tests with equal spaces produces real rotations
    # TODO: the next test is machine-dependent, figure out why
    tformtest0, mm0 = qd_rigid(testimage2, testimage1, mxshift, mxrot2; fvalue=1e-5, rtol=0, print_interval=typemax(Int))

    @test mm0 < 1e-5
    @test isapprox(tformtest0, mytform, atol = 0.1)
    @test isrotation(tformtest0.linear)

    tformtest01, mm01 = RegisterQD.qd_rigid_coarse(testimage2, testimage1, mxshift, [mxrot2...], minwidth_rot; SD=EYE, print_interval=typemax(Int))
    @test isrotation(tformtest01.linear)

    tformtest02 = nothing
    mm02 = nothing
    tformtest02, mm02 = RegisterQD.qd_rigid_fine(testimage2, testimage1, [mxrot2...]./2, minwidth_rot; SD=EYE, print_interval=typemax(Int))
    @test isrotation(tformtest02.linear)

    #with skewed spacing produces real rotations.
    testimage5 = testimage1[:,:,-4:2:5]
    testimage6 = testimage2[:,:,-4:2:5]
    ps = (1, 1, 2)
    SD =SArray{Tuple{3,3}}(Diagonal(SVector(ps) ./ minimum(ps)))

    tformtest5, mm5 = qd_rigid(testimage6, testimage5, mxshift, mxrot2; SD=SD, print_interval=typemax(Int))
    @test mm5 <1e-4
    @test isapprox(tformtest5, mytform, atol = 1)
    @test isrotation(tformtest5.linear)
    @test !isrotation(SD*tformtest5.linear*inv(SD))

    tformtest55, mm55 = qd_rigid(testimage6, testimage5, mxshift, mxrot2; SD=SD, print_interval=typemax(Int), initial_tfm = tformtest5)
    @test mm55 <1e-4
    @test isapprox(tformtest55, mytform, atol = 1)
    @test isrotation(tformtest55.linear)
    @test !isrotation(SD*tformtest55.linear*inv(SD))

    #coarse and fine produce real rotations
    tformtest6, mm6 = RegisterQD.qd_rigid_coarse(testimage6, testimage5, mxshift, mxrot2, minwidth_rot; SD=SD, print_interval=typemax(Int))
    @test mm6 <1e-4
    @test isrotation(tformtest6.linear)
    @test !isrotation(SD*tformtest6.linear*inv(SD))

    tformtest66, mm66 = RegisterQD.qd_rigid_coarse(testimage6, testimage5, mxshift, mxrot2, minwidth_rot; SD=SD, print_interval=typemax(Int), initial_tfm = tformtest6)
    @test mm66 <1e-4
    @test isrotation(tformtest66.linear)
    @test !isrotation(SD*tformtest66.linear*inv(SD))

    tformtest7, mm7 = RegisterQD.qd_rigid_fine(testimage6, testimage5, [mxrot2...]./2, minwidth_rot; SD=SD, print_interval=typemax(Int))
    @test mm7 <1e-4
    @test isrotation(tformtest7.linear)
    @test !isrotation(SD*tformtest7.linear*inv(SD))

    tformtest77, mm77 = RegisterQD.qd_rigid_fine(testimage6, testimage5, [mxrot2...]./2, minwidth_rot; SD=SD, print_interval=typemax(Int), initial_tfm = tformtest7)
    @test mm77 <1e-4
    @test isrotation(tformtest77.linear)
    @test !isrotation(SD*tformtest77.linear*inv(SD))
    #fails due to specific error in qd_rigid_fine

    # a non rigid initial_tfm does not return a rigid transformation and prints an error message (?)
    tformtest8, mm8 = qd_rigid(testimage6, testimage5, mxshift, mxrot2; SD=SD, print_interval=typemax(Int), initial_tfm = AffineMap(SD\tformtest5.linear*SD, SD\tformtest5.translation))
    @test mm8 <1e-4
    @test isapprox(tformtest8, mytform, atol = 1)
    @test !isrotation(tformtest8.linear)

    mktemp() do path, io
        redirect_stdout(io) do
            tformtest8, mm8 = qd_rigid(testimage6, testimage5, mxshift, mxrot2; SD=SD, print_interval=typemax(Int), initial_tfm = AffineMap(SD\tformtest5.linear*SD, SD\tformtest5.translation))
        end
        flush(io)
        str = read(path, String)
        @test !isempty(str)
    end

end

@testset "Test initial_tfm improves translational alignment for affine" begin
    testimage3 = zeros(10,20,10)
    testimage3 .= NaN
    testimage3[3:8,3:8,3:8] = gradcube
    testimage3 = centered(testimage3)

    testimage4 = zeros(10,20,10)
    testimage4 .= NaN
    testimage4[3:8,13:18,3:8] = gradcube
    testimage4 = centered(testimage4)

    init_tfm = AffineMap(diagm(0 => fill(1.0,3)), [0.0,10.0,0.0])

    tformtest3 = nothing
    mm3 = nothing

    @test try
        tformtest3, mm3 = qd_affine(testimage3, testimage4, mxshift; print_interval=typemax(Int)) #
        false
    catch err
        true  # this should break
    end

    tformtest4, mm4 = qd_affine(testimage3, testimage4, mxshift;  print_interval = typemax(Int), initial_tfm = init_tfm, fvalue = 1e-5) #
    @test mm4 < 1e-4#this mismatch should be lower!
    @test isapprox(tformtest4.translation, [0.0, 10.0, 0.0])

end

@testset "Test initial_tfm improves general alignment for affine" begin
    testimage1 = zeros(Float64,10,10,10)
    testimage1 .= NaN
    testimage1[3:8,3:8,3:8] = gradcube
    testimage1 = centered(testimage1)

    tform = RotXYZ(0.1,0.1,0.1)
    mytform = AffineMap(tform, [0.0,0.0,0.0])
    testimage2 = warp(testimage1, mytform, axes(testimage1))

    tformtest1, mm1 = qd_affine(testimage2, testimage1, mxshift; print_interval=typemax(Int)) #
    #as the max-rotation is set too low, this should give a bad mismatch, but should still work
    @test !(mm1 < 1e-8)
    @test !isapprox(tformtest1, mytform, atol = 0.1)

    tformtest2, mm2 = qd_affine(testimage2, testimage1, mxshift; print_interval=typemax(Int), initial_tfm = mytform, fvalue = 1e-5) #
    #with the initial_tfm being the true tfm, this should give back the true rotation

    @test mm2 < 1e-8
    @test isapprox(tformtest2, mytform, atol = 0.0001) #this may not work so well. see qd_standard.
end #Test initial_tfm improves rotational alignment for rigid


text:
immutable LinearCombination{S,P,T}<:AbstractTensorMap{S,P,T}
    maps::Vector{AbstractTensorMap{S,P}}
    coeffs::Vector{T}
    function LinearCombination(maps::Vector{AbstractTensorMap{S,P}},coeffs::Vector{T})
        N=length(maps)
        N==length(coeffs) || throw(ArgumentError("Number of coefficients doesn't match number of terms"))
        dom=domain(maps[1])
        codom=codomain(maps[1])
        for n=1:N
            domain(maps[n])==dom || throw(SpaceError("Domain mismatch"))
            codomain(maps[n])==codom || throw(SpaceError("Codomain mismatch"))
            promote_type(T,eltype(maps[n]))==T || throw(InexactError())
        end
        new(maps,coeffs)
    end
end

# basic methods
codomain(A::LinearCombination)=codomain(A.maps[1])
domain(A::LinearCombination)=domain(A.maps[1])
Base.isreal(A::LinearCombination)=all(isreal,A.maps) && all(isreal,A.coeffs) # sufficient but not necessary
Base.issym(A::LinearCombination)=all(issym,A.maps) # sufficient but not necessary
Base.ishermitian(A::LinearCombination)=all(ishermitian,A.maps) && all(isreal,A.coeffs) # sufficient but not necessary
Base.isposdef(A::LinearCombination)=all(isposdef,A.maps) && all(isposdef,A.coeffs) # sufficient but not necessary

# adding linear maps
function +{S,P}(A1::LinearCombination{S,P},A2::LinearCombination{S,P})
    domain(A1)==domain(A2) || throw(SpaceError("Domain mismatch"))
    codomain(A1)==codomain(A2) || throw(SpaceError("Coomain mismatch"))
    T=promote_type(eltype(A1),eltype(A2))
    return LinearCombination{S,P,T}(AbstractTensorMap{S,P}[A1.maps...,A2.maps...],T[A1.coeffs...,A2.coeffs...])
end
function +{S,P}(A1::AbstractTensorMap{S,P},A2::LinearCombination{S,P})
    domain(A1)==domain(A2) || throw(SpaceError("Domain mismatch"))
    codomain(A1)==codomain(A2) || throw(SpaceError("Coomain mismatch"))
    T=promote_type(eltype(A1),eltype(A2))
    return LinearCombination{S,P,T}(AbstractTensorMap{S,P}[A1,A2.maps...],T[one(T),A2.coeffs...])
end
+{S,P}(A1::LinearCombination{S,P},A2::AbstractTensorMap{S,P})=+{S,P}(A2,A1)
function +{S,P}(A1::AbstractTensorMap{S,P},A2::AbstractTensorMap{S,P})
    domain(A1)==domain(A2) || throw(SpaceError("Domain mismatch"))
    codomain(A1)==codomain(A2) || throw(SpaceError("Coomain mismatch"))
    T=promote_type(eltype(A1),eltype(A2))
    return LinearCombination{S,P,T}(AbstractTensorMap{S,P}[A1,A2],T[one(T),one(T)])
end
function -{S,P}(A1::LinearCombination{S,P},A2::LinearCombination{S,P})
    domain(A1)==domain(A2) || throw(SpaceError("Domain mismatch"))
    codomain(A1)==codomain(A2) || throw(SpaceError("Coomain mismatch"))
    T=promote_type(eltype(A1),eltype(A2))
    return LinearCombination{S,P,T}(AbstractTensorMap{S,P}[A1.maps...,A2.maps...],T[A1.coeffs...,map(-,A2.coeffs)...])
end
function -{S,P}(A1::AbstractTensorMap{S,P},A2::LinearCombination{S,P})
    domain(A1)==domain(A2) || throw(SpaceError("Domain mismatch"))
    codomain(A1)==codomain(A2) || throw(SpaceError("Coomain mismatch"))
    T=promote_type(eltype(A1),eltype(A2))
    return LinearCombination{S,P,T}(AbstractTensorMap{S,P}[A1,A2.maps...],T[one(T),map(-,A2.coeffs)...])
end
function -{S,P}(A1::LinearCombination{S,P},A2::AbstractTensorMap{S,P})
    domain(A1)==domain(A2) || throw(SpaceError("Domain mismatch"))
    codomain(A1)==codomain(A2) || throw(SpaceError("Coomain mismatch"))
    T=promote_type(eltype(A1),eltype(A2))
    return LinearCombination{S,P,T}(AbstractTensorMap{S,P}[A1.maps...,A2],T[A1.coeffs...,-one(T)])
end
function -{S,P}(A1::AbstractTensorMap{S,P},A2::AbstractTensorMap{S,P})
    domain(A1)==domain(A2) || throw(SpaceError("Domain mismatch"))
    codomain(A1)==codomain(A2) || throw(SpaceError("Coomain mismatch"))
    T=promote_type(eltype(A1),eltype(A2))
    return LinearCombination{S,P,T}(AbstractTensorMap{S,P}[A1,A2],T[one(T),-one(T)])
end

# scalar multiplication
-{S,P,T}(A::AbstractTensorMap{S,P,T})=LinearCombination{S,P,T}(AbstractTensorMap{S,P}[A],[-one(T)])
-(A::LinearCombination)=LinearCombination(A.maps,-A.coeffs)

function *{S,P}(alpha::Number,A::AbstractTensorMap{S,P})
    T=promote_type(eltype(alpha),eltype(A))
    return LinearCombination{S,P,T}(AbstractTensorMap{S,P}[A],T[alpha])
end
*(A::AbstractTensorMap,alpha::Number)=*(alpha,A)

*(alpha::Number,A::LinearCombination)=LinearCombination(A.maps,alpha*A.coeffs)
*(A::LinearCombination,alpha::Number)=*(alpha,A)

\(alpha::Number,A::AbstractTensorMap)=*(1/alpha,A)
/(A::AbstractTensorMap,alpha::Number)=*(1/alpha,A)

# comparison of LinearCombination objects
==(A::LinearCombination,B::LinearCombination)=(eltype(A)==eltype(B) && A.maps==B.maps && A.coeffs==B.coeffs)

# special transposition behavior
transpose(A::LinearCombination)=LinearCombination{eltype(A)}(AbstractLinearMap[transpose(l) for l in A.maps],A.coeffs)
ctranspose(A::LinearCombination)=LinearCombination{eltype(A)}(AbstractLinearMap[ctranspose(l) for l in A.maps],conj(A.coeffs))

# multiplication with vectors
function Base.A_mul_B!{S,P}(y::AbstractTensor{S,P},A::LinearCombination,x::AbstractTensor{S,P})
    # no size checking, will be done by individual maps
    Base.A_mul_B!(y,A.maps[1],x)
    scale!(A.coeffs[1],y)
    if length(A.maps)>1
        z=similar(y)
    end
    for n=2:length(A.maps)
        Base.A_mul_B!(z,A.maps[n],x)
        Base.LinAlg.axpy!(A.coeffs[n],z,y)
    end
    return y
end
function Base.At_mul_B!{S,P}(y::AbstractTensor{S,P},A::LinearCombination,x::AbstractTensor{S,P})
    # no size checking, will be done by individual maps
    Base.At_mul_B!(y,A.maps[1],x)
    scale!(A.coeffs[1],y)
    if length(A.maps)>1
        z=similar(y)
    end
    for n=2:length(A.maps)
        Base.At_mul_B!(z,A.maps[n],x)
        Base.LinAlg.axpy!(A.coeffs[n],z,y)
    end
    return y
end
function Base.Ac_mul_B!{S,P}(y::AbstractTensor{S,P},A::LinearCombination,x::AbstractTensor{S,P})
    # no size checking, will be done by individual maps
    Base.Ac_mul_B!(y,A.maps[1],x)
    scale!(conj(A.coeffs[1]),y)
    if length(A.maps)>1
        z=similar(y)
    end
    for n=2:length(A.maps)
        Base.Ac_mul_B!(z,A.maps[n],x)
        Base.LinAlg.axpy!(conj(A.coeffs[n]),z,y)
    end
    return y
end


text: using FuzzySets, Test