text:
# Goal: Make sure some combinations that raise errors do so
using Test
using OrdinaryDiffEq
using LinearAlgebra
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear, prob_ode_lotkavoltera, prob_ode_fitzhughnagumo

import ODEFilters: remake_prob_with_jac

@testset "EK1 requires Jac" begin
    prob = prob_ode_lotkavoltera
    @test_throws ErrorException solve(prob, EK1())
    @test solve(remake_prob_with_jac(prob), EK1()) isa ODEFilters.ProbODESolution
end

@testset "One-dim problems don't work so far!" begin
    prob = prob_ode_linear
    @test_throws ErrorException solve(prob, EK0())
end

@testset "Fixed-timestep requires dt" begin
    prob = prob_ode_lotkavoltera
    @test_throws ErrorException solve(prob, EK0(), adaptive=false)
    @test solve(prob, EK0(), adaptive=false, dt=0.05) isa ODEFilters.ProbODESolution
end


text:
using Revise
using BenchmarkTools
using UnicodePlots

using jLRS

# How do we want to use this tool?
##################################

# Instantiate a sampling strategy

# Example 1:
# Look at some location and how OCO-x samples in some radius of FoCo
# ---------

s1 = jLRS.OCOSampling(-105.0844, # target lon
                      40.5853, # target lat
                      50, # radius in km
                      "/home/psomkuti/oco3.db")

s2 = jLRS.OCOSampling(-104.0844, # target lon
                      40.5853, # target lat
                      50, # radius in km
                      "/home/psomkuti/oco2.db")

s3 = s1 + s2

# Using the locations and samplings, produce a light response curve given some sampling
weekly = jLRS.WeeklySampling(1)
halfweekly = jLRS.WeeklySampling(0.5)

noagg = jLRS.FullROI()
bigagg = jLRS.RegularGridCells(2, 2)

# Perform the aggregation step
agg1 = jLRS.aggregate_scenes(s3, noagg, weekly)
#agg2 = jLRS.aggregate_scenes(s1, noagg, halfweekly)
#agg3 = jLRS.aggregate_scenes(s1, bigagg, weekly)
#agg4 = jLRS.aggregate_scenes(s1, bigagg, halfweekly)

# Example 2:
# Look at some bounding box (lower left, upper right) with OCO-x samples
# ---------
# s2 = OCOSampling(bounding_box=[-55.0, -14.0, -47.0, -5.0],
#                  "/path/to/OCO3")

ext: jl


text:
module TestSetup
    export @title, @subtitle

    macro title(s)
        return :("■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ " * $s)
    end
    macro subtitle(s)
        return :("████████████████████████████████████████████████ " * $s)
    end
end

module OnlineStatsTest
    include("testfiles/messy_output_test.jl")
    include("testfiles/fit_methods_test.jl")
    include("testfiles/multivariate_test.jl")
    include("testfiles/summary_test.jl")
    include("testfiles/weight_test.jl")
    include("testfiles/distributions_test.jl")
    include("testfiles/modeling_test.jl")
    include("testfiles/statlearn_test.jl")
    include("testfiles/streamstats_test.jl")
end


text:
using TensorFlow
using Distributions
include("mnist_loader.jl")

loader = DataLoader()

session = Session(Graph())

function weight_variable(shape)
    initial = map(Float32, rand(Normal(0, .001), shape...))
    return Variable(initial)
end

function bias_variable(shape)
    initial = fill(Float32(.1), shape...)
    return Variable(initial)
end

function conv2d(x, W)
    nn.conv2d(x, W, [1, 1, 1, 1], "SAME")
end

function max_pool_2x2(x)
    nn.max_pool(x, [1, 2, 2, 1], [1, 2, 2, 1], "SAME")
end

@tf begin

    x = placeholder(Float32)
    y_ = placeholder(Float32)

    W_conv1 = weight_variable([5, 5, 1, 32])
    b_conv1 = bias_variable([32])

    x_image = reshape(x, [-1, 28, 28, 1])

    h_conv1 = nn.relu(conv2d(x_image, W_conv1) + b_conv1)
    h_pool1 = max_pool_2x2(h_conv1)

    W_conv2 = weight_variable([5, 5, 32, 64])
    b_conv2 = bias_variable([64])

    h_conv2 = nn.relu(conv2d(h_pool1, W_conv2) + b_conv2)
    h_pool2 = max_pool_2x2(h_conv2)

    W_fc1 = weight_variable([7*7*64, 1024])
    b_fc1 = bias_variable([1024])

    h_pool2_flat = reshape(h_pool2, [-1, 7*7*64])
    h_fc1 = nn.relu(h_pool2_flat * W_fc1 + b_fc1)

    keep_prob = placeholder(Float32)
    h_fc1_drop = nn.dropout(h_fc1, keep_prob)

    W_fc2 = weight_variable([1024, 10])
    b_fc2 = bias_variable([10])

    global y_conv = nn.softmax(h_fc1_drop * W_fc2 + b_fc2)

    global cross_entropy = reduce_mean(-reduce_sum(y_.*log(y_conv), axis=[2]))
end

train_step = train.minimize(train.AdamOptimizer(1e-4), cross_entropy)

correct_prediction = indmax(y_conv, 2) .== indmax(y_, 2)

accuracy = reduce_mean(cast(correct_prediction, Float32))

run(session, global_variables_initializer())

for i in 1:200
    batch = next_batch(loader, 50)
    if i%100 == 1
        train_accuracy = run(session, accuracy, Dict(x=>batch[1], y_=>batch[2], keep_prob=>1.0))
        info("step $i, training accuracy $train_accuracy")
    end
    run(session, train_step, Dict(x=>batch[1], y_=>batch[2], keep_prob=>.5))
end

testx, testy = load_test_set()
test_accuracy = run(session, accuracy, Dict(x=>testx, y_=>testy, keep_prob=>1.0))
info("test accuracy $test_accuracy")

visualize()

max_stars_repo_name: octonion/tensorflow

max_issues_repo_name: octonion/tensorflow

max_forks_repo_name: octonion/tensorflow


text:
@testset "Decomposing derivative" begin
    cases = [
        ["yy_11", ["y", "yy", "yy_"], ("yy", 11)],
        ["xx_xx_xx_0", ["xx", "x", "xx_xx_xx"], ("xx_xx_xx", 0)],
        ["abc154f", ["ab", "abc"], nothing],
        ["c_1542673", ["a", "b", "c"], ("c", 1542673)],
        ["a", ["a"], nothing]
    ]
    for c in cases
        @test decompose_derivative(c[1], c[2]) == c[3]
    end
end


text:
module TestNetworkRoles
   using Base.Test
   using Brim
   #=testdir = dirname(@__FILE__)=#
   A = [1 0 0; 1 1 0; 0 1 1; 0 1 1]
   S = [1 0; 0 1; 0 1; 1 0; 1 0; 0 1; 0 1]
   M = Modular(A, S)
   roles = network_roles(M)
   ## Approximate values, calculated by hand
   @test abs(roles[:z][1] - 1.1547) <= 0.01
   @test abs(roles[:c][2] - 0.4444) <= 0.01
   @test roles[:c][5] == 0.5
end


text:
#=
Case folding for UTF-16 encoded strings

Copyright 2018 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md
=#

# These are more complex, and maybe belong in a separate UTF16Str.jl package

function _lower(::Type{<:Str{UTF16CSE}}, beg, off, len)
    buf, out = _allocate(UInt16, len)
    unsafe_copyto!(out, beg, len)
    fin = out + (len<<1)
    out += off
    while out < fin
        ch = get_codeunit(out)
        if ch <= 0xff
            _is_upper_al(ch) && set_codeunit!(out, ch += 0x20)
        elseif is_surrogate_trail(ch)
            # pick up previous code unit (lead surrogate)
            cp = get_supplementary(get_codeunit(out - 2), ch)
            if cp < 0x1ffff && _can_lower_slp(cp)
                w1, w2 = get_utf16(_lower_slp(cp))
                set_codeunit!(out - 2, w1)
                set_codeunit!(out,     w2)
            end
        elseif _can_lower_bmp(ch)
            set_codeunit!(out, _lower_bmp(ch))
        end
        out += 2
    end
    Str(UTF16CSE, buf)
end

function lowercase(str::MaybeSub{S}) where {C<:UTF16CSE,S<:Str{C}}
    @preserve str begin
        pnt = beg = pointer(str)
        fin = beg + sizeof(str)
        while pnt < fin
            ch = get_codeunit(pnt)
            prv = pnt
            (ch <= 0xff ? _is_upper_al(ch) :
             (is_surrogate_lead(ch)
              ? (cp = get_supplementary(ch, get_codeunit(pnt += 2));
                 cp <= 0x1ffff && _can_lower_slp(cp))
              : _can_lower_bmp(ch))) &&
                  return _lower(UTF16Str, beg, prv-beg, ncodeunits(str))
            pnt += 2
        end
    end
    str
end

function lowercase_first(str::MaybeSub{S}) where {C<:UTF16CSE,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = pointer(str)
        ch = get_codeunit(pnt)
        if ch <= 0x7f
            _is_upper_al(ch) || return str
            ch += 0x20
            buf, out = _allocate(codeunit(C), len)
            len > 1 && unsafe_copyto!(out, pnt, len)
        elseif is_surrogate_lead(ch)
            cp = get_supplementary(ch, get_codeunit(pnt + 2))
            (cp <= 0x1ffff && _can_lower_slp(cp)) || return str
            ch, c2 = get_utf16(_lower_slp(cp))
            buf, out = _allocate(codeunit(C), len)
            len > 1 && unsafe_copyto!(out, pnt, len)
            set_codeunit!(out + 2, c2)
        else
            _can_lower_ch(ch) || return str
            ch = _lowercase(ch)
            buf, out = _allocate(codeunit(C), len)
            len > 1 && unsafe_copyto!(out, pnt, len)
        end
        set_codeunit!(out, ch)
        Str(C, buf)
    end
end

function _upper(::Type{<:Str{UTF16CSE}}, beg, off, len)
    buf, out = _allocate(UInt16, len)
    unsafe_copyto!(out, beg, len)
    fin = out + (len<<1)
    out += off
    while out < fin
        ch = get_codeunit(out)
        if ch <= 0x7f
            _is_lower_a(ch) && set_codeunit!(out, ch -= 0x20)
        elseif ch <= 0xff
            set_codeunit!(out, _uppercase_l(ch))
        elseif is_surrogate_trail(ch)
            # pick up previous code unit (lead surrogate)
            cp = get_supplementary(get_codeunit(out - 2), ch)
            if cp <= 0x1ffff && _can_upper_slp(cp)
                w1, w2 = get_utf16(_upper_slp(cp))
                set_codeunit!(out - 2, w1)
                set_codeunit!(out,     w2)
            end
        elseif _can_upper_bmp(ch)
            set_codeunit!(out, _upper_bmp(ch))
        end
        out += 2
    end
    Str(UTF16CSE, buf)
end

function uppercase(str::MaybeSub{S}) where {C<:UTF16CSE,S<:Str{C}}
    @preserve str begin
        pnt = beg = pointer(str)
        fin = beg + sizeof(str)
        while pnt < fin
            ch = get_codeunit(pnt)
            prv = pnt
            if (ch <= 0x7f ? _is_lower_a(ch) : ch <= 0xff ? _can_upper_lat(ch) :
                (is_surrogate_lead(ch)
                 ? (cp = get_supplementary(ch, get_codeunit(pnt += 2));
                    cp <= 0x1ffff && _can_upper_slp(cp))
                 : _can_upper_bmp(ch)))
                return _upper(UTF16Str, beg, prv-beg, ncodeunits(str))
            end
            pnt += 2
        end
    end
    str
end

function uppercase_first(str::MaybeSub{S}) where {C<:UTF16CSE,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = pointer(str)
        ch = get_codeunit(pnt)
        if ch <= 0x7f
            _is_lower_a(ch) || return str
            ch -= 0x20
            buf, out = _allocate(codeunit(C), len)
            len > 1 && unsafe_copyto!(out, pnt, len)
        elseif is_surrogate_lead(ch)
            cp = get_supplementary(ch, get_codeunit(pnt + 2))
            (cp <= 0x1ffff && _can_upper_slp(cp)) || return str
            buf, out = _allocate(codeunit(C), len)
            len > 1 && unsafe_copyto!(out, pnt, len)
            ch, c2 = get_utf16(_upper_slp(cp))
            set_codeunit!(out + 2, c2)
        else
            cp = ch
            (ch = _titlecase(ch)) == cp && return str
            buf, out = _allocate(codeunit(C), len)
            len > 1 && unsafe_copyto!(out, pnt, len)
        end
        set_codeunit!(out, ch)
        Str(C, buf)
    end
end


text:
@recipe function plot(lasso::LassoResult)
    yguide --> "parameter"
    xguide --> "lambda"
    xflip --> true
    lab := hcat(name_regressors(lasso.mdl)...)
    (lasso.λ, lasso.β')
end


text:
module Utilities

using ..Turing: AbstractSampler, Sampler
using Distributions, Bijectors
using StatsFuns, SpecialFunctions
using MCMCChains: AbstractChains, Chains, setinfo
import Distributions: sample

export  vectorize,
        reconstruct,
        reconstruct!,
        Sample,
        Chain,
        init,
        vectorize

include("helper.jl")
include("robustinit.jl")
include("io.jl")           # I/O

end # module


text:
#=
Utilities used by UMAP.jl
=#

@inline fit_ab(_, __, a, b) = a, b

"""
    fit_ab(min_dist, spread, _a, _b) -> a, b

Find a smooth approximation to the membership function of points embedded in ℜᵈ.
This fits a smooth curve that approximates an exponential decay offset by `min_dist`,
returning the parameters `(a, b)`.
"""
function fit_ab(min_dist, spread, ::Nothing, ::Nothing)
    ψ(d) = d >= min_dist ? exp(-(d - min_dist)/spread) : 1.
    xs = LinRange(0., spread*3, 300)
    ys = map(ψ, xs)
    @. curve(x, p) = (1. + p[1]*x^(2*p[2]))^(-1)
    result = curve_fit(curve, xs, ys, [1., 1.], lower=[0., -Inf])
    a, b = result.param
    return a, b
end

knn_search(X::AbstractMatrix, k, metric::Symbol) = knn_search(X, k, Val(metric))

# treat given matrix `X` as distance matrix
knn_search(X::AbstractMatrix, k, ::Val{:precomputed}) = _knn_from_dists(X, k)

"""
    knn_search(X, k, metric) -> knns, dists

Find the `k` nearest neighbors of each point.

`metric` may be of type:
- ::Symbol - `knn_search` is dispatched to one of the following based on the evaluation of `metric`:
- ::Val(:precomputed) - computes neighbors from `X` treated as a precomputed distance matrix.
- ::SemiMetric - computes neighbors from `X` treated as samples, using the given metric.

# Returns
- `knns`: `knns[j, i]` is the index of node i's jth nearest neighbor.
- `dists`: `dists[j, i]` is the distance of node i's jth nearest neighbor.
"""
function knn_search(X::AbstractMatrix,
                    k,
                    metric::SemiMetric)
    if size(X, 2) < 4096
        return knn_search(X, k, metric, Val(:pairwise))
    else
        return knn_search(X, k, metric, Val(:approximate))
    end
end

# compute all pairwise distances
# return the nearest k to each point v, other than v itself
function knn_search(X::AbstractMatrix{S},
                    k,
                    metric,
                    ::Val{:pairwise}) where {S <: Real}
    num_points = size(X, 2)
    dist_mat = Array{S}(undef, num_points, num_points)
    pairwise!(dist_mat, metric, X, dims=2)
    # all_dists is symmetric distance matrix
    return _knn_from_dists(dist_mat, k)
end

# find the approximate k nearest neighbors using NNDescent
function knn_search(X::AbstractMatrix{S},
                    k,
                    metric,
                    ::Val{:approximate}) where {S <: Real}
    knngraph = nndescent(X, k, metric)
    return knn_matrices(knngraph)
end

"""
    knn_search(X, Q, k, metric, knns, dists) -> knns, dists

Given a matrix `X` and a matrix `Q`, use the given metric to compute the `k` nearest neighbors out of the
columns of `X` from the queries (columns in `Q`).
If the matrices are large, reconstruct the approximate nearest neighbors graph of `X` using the given `knns` and `dists`,
representing indices and distances of pairwise neighbors of `X`, and use this to search for approximate nearest
neighbors of `Q`.
If the matrices are small, search for exact nearest neighbors of `Q` by computing all pairwise distances with `X`.

`metric` may be of type:
- ::Symbol - `knn_search` is dispatched to one of the following based on the evaluation of `metric`:
- ::Val(:precomputed) - computes neighbors from `X` treated as a precomputed distance matrix.
- ::SemiMetric - computes neighbors from `X` treated as samples, using the given metric.

# Returns
- `knns`: `knns[j, i]` is the index of node i's jth nearest neighbor.
- `dists`: `dists[j, i]` is the distance of node i's jth nearest neighbor.
"""
function knn_search(X::AbstractMatrix,
                    Q::AbstractMatrix,
                    k::Integer,
                    metric::SemiMetric,
                    knns::AbstractMatrix{<:Integer},
                    dists::AbstractMatrix{<:Real})
    if size(X, 2) < 4096
        return _knn_from_dists(pairwise(metric, X, Q, dims=2), k, ignore_diagonal=false)
    else
        knngraph = HeapKNNGraph(collect(eachcol(X)), metric, knns, dists)
        return search(knngraph, collect(eachcol(Q)), k; max_candidates=8*k)
    end
end

function _knn_from_dists(dist_mat::AbstractMatrix{S}, k::Integer; ignore_diagonal=true) where {S <: Real}
    # Ignore diagonal 0 elements (which will be smallest) when distance matrix represents pairwise distances of the same set
    # If dist_mat represents distances between two different sets, diagonal elements be nontrivial
    range = (1:k) .+ ignore_diagonal
    knns_ = [partialsortperm(view(dist_mat, :, i), range) for i in 1:size(dist_mat, 2)]
    dists_ = [dist_mat[:, i][knns_[i]] for i in eachindex(knns_)]
    knns = hcat(knns_...)::Matrix{Int}
    dists = hcat(dists_...)::Matrix{S}
    return knns, dists
end

# combine local fuzzy simplicial sets
@inline function combine_fuzzy_sets(fs_set::AbstractMatrix{T},
                                    set_op_ratio) where {T}
    return set_op_ratio .* fuzzy_set_union(fs_set) .+
           (one(T) - set_op_ratio) .* fuzzy_set_intersection(fs_set)
end

@inline function fuzzy_set_union(fs_set::AbstractMatrix)
    return fs_set .+ fs_set' .- (fs_set .* fs_set')
end

@inline function fuzzy_set_intersection(fs_set::AbstractMatrix)
    return fs_set .* fs_set'
end


text:
"""
Runs the pagerank algorithm on the given graph. It returns the nodes in the order in which they are ranked, and their respective rank scores.

Input:
-----
* G::SparseMatrixCSC

Output:
------
* top_nodes::Vector{Int32}
* top_ranks:;Vector{Float32}

"""
function pagerank(G::SparseMatrixCSC)
	rows = G.colptr - 1
    cols = G.rowval - 1
    nodes = length(rows) - 1
    edges = length(cols)
    rows = map(Int32, rows)
    cols = map(Int32, cols)
    top_nodes = Array(Int32, nodes)
    top_ranks = Array(Cfloat, nodes)
	ccall((:pagerank, gunrock_lib), Void, (Ptr{Cint}, Ptr{Cfloat}, Int, Int, Ptr{Cint}, Ptr{Cint}), top_nodes, top_ranks, nodes, edges, rows, cols)
    top_nodes += 1
	top_nodes, top_ranks
end


text:
# Types -------

type Organization <: Owner
    id
    email
    name
    login
    blog
    bio
    location
    gravatar_id
    avatar_url
    company

    public_repos
    owned_private_repos
    total_private_repos

    public_gists
    private_gists

    followers
    following
    collaborators

    updated_at
    created_at

    plan
    user_type
    site_admin
    disk_usage

    function Organization(data::Dict)
        new(get(data, "id", nothing),
            get(data, "email", nothing),
            get(data, "name", nothing),
            get(data, "login", nothing),
            get(data, "blog", nothing),
            get(data, "bio", nothing),
            get(data, "location", nothing),
            get(data, "gravatar_id", nothing),
            get(data, "avatar_url", nothing),
            get(data, "company", nothing),
            get(data, "public_repos", nothing),
            get(data, "owned_private_repos", nothing),
            get(data, "total_private_repos", nothing),
            get(data, "public_gists", nothing),
            get(data, "private_gists", nothing),
            get(data, "followers", nothing),
            get(data, "following", nothing),
            get(data, "collaborators", nothing),
            get(data, "updated_at", nothing),
            get(data, "created_at", nothing),
            get(data, "plan", nothing),
            get(data, "type", nothing),
            get(data, "site_admin", nothing),
            get(data, "disk_usage", nothing))
    end
end

function show(io::IO, org::Organization)
    print(io, "$User - $(org.login)")

    if org.name != nothing && !isempty(org.name) && org.blog != nothing && !isempty(org.blog)
        print(io, " ($(org.name), $(org.blog))")
    elseif org.name != nothing
        print(io, " ($(org.name))")
    elseif org.blog != nothing
        print(io, " ($(org.blog))")
    end

    if org.bio != nothing && !isempty(org.bio)
        print(io, "\n\"$(org.bio)\"")
    end
end

# Interface -------

function org(name; auth = AnonymousAuth(), options...)
    org(auth, name; options...)
end

function org(auth::Authorization, name; headers = Dict(), options...)
    authenticate_headers(headers, auth)
    r = get(URI(API_ENDPOINT; path = "/orgs/$name");
            headers = headers,
            options...)

    handle_error(r)

    Organization(JSON.parse(r.data))
end

function orgs(user::String; auth = AnonymousAuth(), options...)
    orgs(auth, user; options...)
end

function orgs(user::User; auth = AnonymousAuth(), options...)
    orgs(auth, user.login; options...)
end

function orgs(auth::Authorization, user; headers = Dict(), options...)
    authenticate_headers(headers, auth)
    r = get(URI(API_ENDPOINT; path = "/users/$user/orgs");
            headers = headers,
            options...)

    handle_error(r)

    map!(o -> Organization(o), JSON.parse(r.data))
end


text:
qmax_default(alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}) = isadaptive(alg) ? 9//8 : 0
qmin_default(alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}) = isadaptive(alg) ? 1//5 : 0

delta_default(alg) = 1//1
delta_default(alg::SRIW1) = 1//6

isadaptive(alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}) = false
isadaptive(alg::Union{StochasticDiffEqAdaptiveAlgorithm,StochasticDiffEqRODEAdaptiveAlgorithm,StochasticDiffEqJumpAdaptiveAlgorithm,StochasticDiffEqJumpDiffusionAdaptiveAlgorithm}) = true
isadaptive(alg::Union{StochasticDiffEqCompositeAlgorithm,StochasticDiffEqRODECompositeAlgorithm}) = all(isadaptive.(alg.algs))
isadaptive(prob,alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}) = isadaptive(alg)
isadaptive(prob::JumpProblem,alg::ImplicitEM) = false

# For whether an algorithm uses a priori dt estimates or utilizes an error estimate
isaposteriori(alg) = false
isaposteriori(alg::CaoTauLeaping) = true

alg_order(alg::EM) = 1//2
alg_order(alg::LambaEM) = 1//2
alg_order(alg::ImplicitEM) = 1//2
alg_order(alg::ImplicitEulerHeun) = 1//2
alg_order(alg::ImplicitRKMil) = 1//1
alg_order(alg::WangLi3SMil_A) = 1//1
alg_order(alg::WangLi3SMil_B) = 1//1
alg_order(alg::WangLi3SMil_C) = 1//1
alg_order(alg::WangLi3SMil_D) = 1//1
alg_order(alg::WangLi3SMil_E) = 1//1
alg_order(alg::WangLi3SMil_F) = 1//1
alg_order(alg::ISSEM) = 1//2
alg_order(alg::ISSEulerHeun) = 1//2
alg_order(alg::SplitEM) = 1//2
alg_order(alg::PCEuler) = 1//2
alg_order(alg::IIF1M) = 1//2
alg_order(alg::IIF2M) = 1//2
alg_order(alg::IIF1Mil) = 1//1
alg_order(alg::EulerHeun) = 1//2
alg_order(alg::LambaEulerHeun) = 1//2
alg_order(alg::RandomEM) = 1//2
alg_order(alg::SimplifiedEM) = 1//2
alg_order(alg::RKMil) = 1//1
alg_order(alg::RKMilCommute) = 1//1
alg_order(alg::RKMil_General) = 1//1

# Generalised version of SROCK1, both Ito ans Stratonovich, will have strong order of 1//2
# and weak order of 1 for Multidimensional Weiner process
# Stratonovich version strong order 1 for 1 dimensional Weiner Process or if noise is commutative
# Ito version can have strong order version for 1 dimensional Weiner Process,
# diagonal noise or commutative noise
alg_order(alg::SROCK1) = 1//2
alg_order(alg::SROCK2) = 1//1
alg_order(alg::KomBurSROCK2) = 1//1
alg_order(alg::SROCKC2) = 1//1
alg_order(alg::SROCKEM) = alg.strong_order_1 ? 1//1 : 1//2
alg_order(alg::SKSROCK) = 1//2
alg_order(alg::TangXiaoSROCK2) = 1//1

alg_order(alg::SRI) = alg.tableau.order
alg_order(alg::SRIW1) = 3//2
alg_order(alg::SRIW2) = 3//2
alg_order(alg::SOSRI) = 3//2
alg_order(alg::SOSRI2) = 3//2
alg_order(alg::SRA) = alg.tableau.order
alg_order(alg::SRA1) = 2//1
alg_order(alg::SRA2) = 2//1
alg_order(alg::SRA3) = 2//1
alg_order(alg::SOSRA) = 2//1
alg_order(alg::SOSRA2) = 2//1

alg_order(alg::DRI1) = 1//1
alg_order(alg::RI1) = 1//1

alg_order(alg::TauLeaping) = 1//1
alg_order(alg::CaoTauLeaping) = 1//1

alg_order(alg::SKenCarp) = 2//1
alg_order(alg::Union{StochasticDiffEqCompositeAlgorithm,StochasticDiffEqRODECompositeAlgorithm}) = maximum(alg_order.(alg.algs))
get_current_alg_order(alg::StochasticDiffEqAlgorithm,cache) = alg_order(alg)
get_current_alg_order(alg::Union{StochasticDiffEqCompositeAlgorithm,StochasticDiffEqRODECompositeAlgorithm},cache) = alg_order(alg.algs[cache.current])

beta2_default(alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}) = isadaptive(alg) ? 2//(5alg_order(alg)) : 0
beta1_default(alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm},beta2) = isadaptive(alg) ? 7//(10alg_order(alg)) : 0

isdtchangeable(alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}) = true

alg_interpretation(alg::StochasticDiffEqAlgorithm) = :Ito
alg_interpretation(alg::EulerHeun) = :Stratonovich
alg_interpretation(alg::LambaEulerHeun) = :Stratonovich
alg_interpretation(alg::KomBurSROCK2) = :Stratonovich
alg_interpretation(alg::RKMil{interpretation}) where {interpretation} = interpretation
alg_interpretation(alg::SROCK1{interpretation}) where {interpretation} = interpretation
alg_interpretation(alg::RKMilCommute{interpretation}) where {interpretation} = interpretation
alg_interpretation(alg::RKMil_General) = alg.interpretation
alg_interpretation(alg::ImplicitRKMil{CS,AD,F,S,N,T2,Controller,interpretation}) where {CS,AD,F,S,N,T2,Controller,interpretation} = interpretation

alg_compatible(prob,alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}) = true
alg_compatible(prob,alg::StochasticDiffEqAlgorithm) = false

function alg_compatible(prob::JumpProblem,alg::StochasticDiffEqAlgorithm)
    alg_compatible(prob.prob,alg) && prob.regular_jump === nothing &&
    typeof(prob.prob) <: DiffEqBase.AbstractSDEProblem
end
alg_compatible(prob::JumpProblem,alg::EM) = alg_compatible(prob.prob,alg)
alg_compatible(prob::JumpProblem,alg::ImplicitEM) = alg_compatible(prob.prob,alg)

alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::SRI) = is_diagonal_noise(prob)
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::SRIW1) = is_diagonal_noise(prob)
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::SRIW2) = is_diagonal_noise(prob)
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::SOSRI) = is_diagonal_noise(prob)
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::SOSRI2) = is_diagonal_noise(prob)
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::SRA) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::SRA1) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::SRA2) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::SRA3) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::SOSRA) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::SOSRA2) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::DRI1) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::RI1) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::SKenCarp) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::EM) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::LambaEM) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::WangLi3SMil_A) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::WangLi3SMil_B) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::WangLi3SMil_C) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::WangLi3SMil_D) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::WangLi3SMil_E) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::WangLi3SMil_F) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::SROCK1) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::SROCK2) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::KomBurSROCK2) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::SROCKC2) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::SROCKEM) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::SKSROCK) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::TangXiaoSROCK2) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::EulerHeun) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::LambaEulerHeun) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::SplitEM) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::PCEuler) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::ImplicitEM) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::ImplicitEulerHeun) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::ISSEM) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::ISSEulerHeun) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::SimplifiedEM) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::RKMil) = is_diagonal_noise(prob)
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::ImplicitRKMil) = is_diagonal_noise(prob)
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::RKMilCommute) = true # No good check for commutative noise
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::RKMil_General) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::IIF1M) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::IIF2M) = true
alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::Union{StochasticDiffEqCompositeAlgorithm,StochasticDiffEqRODECompositeAlgorithm}) = max((alg_compatible(prob,a) for a in alg.algs)...)

function alg_compatible(prob::JumpProblem,alg::Union{StochasticDiffEqJumpAdaptiveAlgorithm,StochasticDiffEqJumpAlgorithm})
    prob.prob isa DiscreteProblem
end

alg_needs_extra_process(alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}) = false
alg_needs_extra_process(alg::Union{StochasticDiffEqCompositeAlgorithm,StochasticDiffEqRODECompositeAlgorithm}) = max((alg_needs_extra_process(a) for a in alg.algs)...)
alg_needs_extra_process(alg::SRI) = true
alg_needs_extra_process(alg::SRIW1) = true
alg_needs_extra_process(alg::SRIW2) = true
alg_needs_extra_process(alg::SOSRI) = true
alg_needs_extra_process(alg::SOSRI2) = true
alg_needs_extra_process(alg::SRA) = true
alg_needs_extra_process(alg::SRA1) = true
alg_needs_extra_process(alg::SRA2) = true
alg_needs_extra_process(alg::SRA3) = true
alg_needs_extra_process(alg::SOSRA) = true
alg_needs_extra_process(alg::SOSRA2) = true
alg_needs_extra_process(alg::SKenCarp) = true
alg_needs_extra_process(alg::DRI1) = true
alg_needs_extra_process(alg::RI1) = true

OrdinaryDiffEq.alg_autodiff(alg::StochasticDiffEqNewtonAlgorithm{CS,AD,Controller}) where {CS,AD,Controller} = AD
OrdinaryDiffEq.alg_autodiff(alg::StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}) where {CS,AD,Controller} = AD
OrdinaryDiffEq.alg_autodiff(alg::StochasticDiffEqJumpNewtonAdaptiveAlgorithm{CS,AD,Controller}) where {CS,AD,Controller} = AD
OrdinaryDiffEq.alg_autodiff(alg::StochasticDiffEqJumpNewtonDiffusionAdaptiveAlgorithm{CS,AD,Controller}) where {CS,AD,Controller} = AD

OrdinaryDiffEq.get_current_alg_autodiff(alg::StochasticDiffEqCompositeAlgorithm, cache) = OrdinaryDiffEq.alg_autodiff(alg.algs[cache.current])

OrdinaryDiffEq.get_chunksize(alg::StochasticDiffEqNewtonAlgorithm{CS,AD,Controller}) where {CS,AD,Controller} = CS
OrdinaryDiffEq.get_chunksize(alg::StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}) where {CS,AD,Controller} = CS
OrdinaryDiffEq.get_chunksize(alg::StochasticDiffEqJumpNewtonAdaptiveAlgorithm{CS,AD,Controller}) where {CS,AD,Controller} = CS
OrdinaryDiffEq.get_chunksize(alg::StochasticDiffEqJumpNewtonDiffusionAdaptiveAlgorithm{CS,AD,Controller}) where {CS,AD,Controller} = CS

alg_mass_matrix_compatible(alg::StochasticDiffEqAlgorithm) = false

alg_can_repeat_jac(alg::StochasticDiffEqAlgorithm) = true

function alg_mass_matrix_compatible(alg::Union{StochasticDiffEqNewtonAlgorithm,StochasticDiffEqNewtonAdaptiveAlgorithm})
    if alg.symplectic
        return true
    elseif alg.theta == 1
        return true
    else
        error("Algorithm must be set as symplectic or theta=1 for mass matrices")
    end
end

is_split_step(::StochasticDiffEqAlgorithm) = false
is_split_step(::EM{split}) = split
is_split_step(::LambaEM{split}) = split

alg_stability_size(alg::SOSRI2) = 10.6
alg_stability_size(alg::SOSRA2) = 5.3

is_composite(alg) = false
is_composite(alg::StochasticDiffEqCompositeAlgorithm) = true
is_composite(alg::StochasticDiffEqRODECompositeAlgorithm) = true
function unwrap_alg(integrator, is_nlsolve)
  alg = integrator.alg
  if !is_composite(alg)
    return alg
  elseif typeof(alg.choice_function) <: AutoSwitch
    num = is_nlsolve ? 2 : 1
    return alg.algs[num]
  else
    return alg.algs[integrator.cache.current]
  end
end

issplit(::StochasticDiffEqAlgorithm) = false
issplit(::SplitSDEAlgorithms) = true

function OrdinaryDiffEq.unwrap_alg(integrator::SDEIntegrator, is_stiff)
  alg = integrator.alg
  if !is_composite(alg)
    return alg
  elseif typeof(alg.choice_function) <: AutoSwitch
    num = is_stiff ? 2 : 1
    return alg.algs[num]
  else
    return alg.algs[integrator.cache.current]
  end
end

alg_control_rate(::StochasticDiffEqAlgorithm) = false
alg_control_rate(::TauLeaping) = true


text:
int_rules_4_5_2_1 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.5.2.1*(a+b*sec)^m*(c+d*sec)^n =#
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~c) ^ ~n * Antiderivative(ExpandTrig((1 + (~d / ~c) * csc(~e + ~f * ~x)) ^ ~n, (~a + ~b * csc(~e + ~f * ~x)) ^ ~m, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (IGtQ(~m, 0) && (ILtQ(~n, 0) && LtQ(~m + ~n, 2)))))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (-(~a) * ~c) ^ ~m * Antiderivative(cot(~e + ~f * ~x) ^ (2 * ~m) * (~c + ~d * csc(~e + ~f * ~x)) ^ (~n - ~m), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (IntegerQ(~m) && (RationalQ(~n) && Not(IntegerQ(~n) && GtQ(~m - ~n, 0))))))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (((-(~a) * ~c) ^ (~m + 1 / 2) * cot(~e + ~f * ~x)) / (sqrt(~a + ~b * csc(~e + ~f * ~x)) * sqrt(~c + ~d * csc(~e + ~f * ~x)))) * Antiderivative(cot(~e + ~f * ~x) ^ (2 * ~m), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && IntegerQ(~m + 1 / 2)))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (2 * ~a * ~c * cot(~e + ~f * ~x) * (~c + ~d * csc(~e + ~f * ~x)) ^ (~n - 1)) / (~f * (2 * ~n - 1) * sqrt(~a + ~b * csc(~e + ~f * ~x))) + ~c * Antiderivative(sqrt(~a + ~b * csc(~e + ~f * ~x)) * (~c + ~d * csc(~e + ~f * ~x)) ^ (~n - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && GtQ(~n, 1 / 2)))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (-2 * ~a * cot(~e + ~f * ~x) * (~c + ~d * csc(~e + ~f * ~x)) ^ ~n) / (~f * (2 * ~n + 1) * sqrt(~a + ~b * csc(~e + ~f * ~x))) + (1 / ~c) * Antiderivative(sqrt(~a + ~b * csc(~e + ~f * ~x)) * (~c + ~d * csc(~e + ~f * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && LtQ(~n, -1 / 2)))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ (3 / 2) * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (-4 * (~a) ^ 2 * cot(~e + ~f * ~x) * (~c + ~d * csc(~e + ~f * ~x)) ^ ~n) / (~f * (2 * ~n + 1) * sqrt(~a + ~b * csc(~e + ~f * ~x))) + (~a / ~c) * Antiderivative(sqrt(~a + ~b * csc(~e + ~f * ~x)) * (~c + ~d * csc(~e + ~f * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && LtQ(~n, -1 / 2)))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ (3 / 2) * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (-2 * (~a) ^ 2 * cot(~e + ~f * ~x) * (~c + ~d * csc(~e + ~f * ~x)) ^ ~n) / (~f * (2 * ~n + 1) * sqrt(~a + ~b * csc(~e + ~f * ~x))) + ~a * Antiderivative(sqrt(~a + ~b * csc(~e + ~f * ~x)) * (~c + ~d * csc(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && Not(LeQ(~n, -1 / 2))))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ (5 / 2) * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (-8 * (~a) ^ 3 * cot(~e + ~f * ~x) * (~c + ~d * csc(~e + ~f * ~x)) ^ ~n) / (~f * (2 * ~n + 1) * sqrt(~a + ~b * csc(~e + ~f * ~x))) + ((~a) ^ 2 / (~c) ^ 2) * Antiderivative(sqrt(~a + ~b * csc(~e + ~f * ~x)) * (~c + ~d * csc(~e + ~f * ~x)) ^ (~n + 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && LtQ(~n, -1 / 2)))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((-(~a) * ~c * cot(~e + ~f * ~x)) / (~f * sqrt(~a + ~b * csc(~e + ~f * ~x)) * sqrt(~c + ~d * csc(~e + ~f * ~x)))) * Subst(Antiderivative(((~b + ~a * ~x) ^ (~m - 1 / 2) * (~d + ~c * ~x) ^ (~n - 1 / 2)) / (~x) ^ (~m + ~n), ~x), ~x, sin(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (IntegerQ(~m - 1 / 2) && EqQ(~m + ~n, 0))))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => ((~a * ~c * cot(~e + ~f * ~x)) / (~f * sqrt(~a + ~b * csc(~e + ~f * ~x)) * sqrt(~c + ~d * csc(~e + ~f * ~x)))) * Subst(Antiderivative(((~a + ~b * ~x) ^ (~m - 1 / 2) * (~c + ~d * ~x) ^ (~n - 1 / 2)) / ~x, ~x), ~x, csc(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && EqQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) * (~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => ~a * ~c * ~x + ~b * ~d * Antiderivative(csc(~e + ~f * ~x) ^ 2, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && EqQ(~b * ~c + ~a * ~d, 0)
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) * (~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => ~a * ~c * ~x + (~b * ~c + ~a * ~d) * Antiderivative(csc(~e + ~f * ~x), ~x) + ~b * ~d * Antiderivative(csc(~e + ~f * ~x) ^ 2, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && NeQ(~b * ~c + ~a * ~d, 0))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) * (~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => ~c * Antiderivative(sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) + ~d * Antiderivative(sqrt(~a + ~b * csc(~e + ~f * ~x)) * csc(~e + ~f * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && EqQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) * (~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => ~a * ~c * Antiderivative(1 / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) + Antiderivative((csc(~e + ~f * ~x) * (~b * ~c + ~a * ~d + ~b * ~d * csc(~e + ~f * ~x))) / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && NeQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => (-(~b) * ~d * cot(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ (~m - 1)) / (~f * ~m) + (1 / ~m) * Antiderivative((~a + ~b * csc(~e + ~f * ~x)) ^ (~m - 1) * Simp(~a * ~c * ~m + (~b * ~c * ~m + ~a * ~d * (2 * ~m - 1)) * csc(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (GtQ(~m, 1) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && IntegerQ(2 * ~m))))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => (-(~b) * ~d * cot(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ (~m - 1)) / (~f * ~m) + (1 / ~m) * Antiderivative((~a + ~b * csc(~e + ~f * ~x)) ^ (~m - 2) * Simp((~a) ^ 2 * ~c * ~m + ((~b) ^ 2 * ~d * (~m - 1) + 2 * ~a * ~b * ~c * ~m + (~a) ^ 2 * ~d * ~m) * csc(~e + ~f * ~x) + ~b * (~b * ~c * ~m + ~a * ~d * (2 * ~m - 1)) * csc(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (GtQ(~m, 1) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && IntegerQ(2 * ~m))))
    @apply_utils Antiderivative((~c + ~(d') * csc(~(e') + ~(f') * ~x)) / (~a + ~(b') * csc(~(e') + ~(f') * ~x)), ~x) => (~c * ~x) / ~a - ((~b * ~c - ~a * ~d) / ~a) * Antiderivative(csc(~e + ~f * ~x) / (~a + ~b * csc(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && NeQ(~b * ~c - ~a * ~d, 0)
    @apply_utils Antiderivative((~c + ~(d') * csc(~(e') + ~(f') * ~x)) / sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)), ~x) => (~c / ~a) * Antiderivative(sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) - ((~b * ~c - ~a * ~d) / ~a) * Antiderivative(csc(~e + ~f * ~x) / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && EqQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative((~c + ~(d') * csc(~(e') + ~(f') * ~x)) / sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)), ~x) => ~c * Antiderivative(1 / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) + ~d * Antiderivative(csc(~e + ~f * ~x) / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && NeQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => (-((~b * ~c - ~a * ~d)) * cot(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ ~m) / (~b * ~f * (2 * ~m + 1)) + (1 / ((~a) ^ 2 * (2 * ~m + 1))) * Antiderivative((~a + ~b * csc(~e + ~f * ~x)) ^ (~m + 1) * Simp(~a * ~c * (2 * ~m + 1) - (~b * ~c - ~a * ~d) * (~m + 1) * csc(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (LtQ(~m, -1) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && IntegerQ(2 * ~m))))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => (~b * (~b * ~c - ~a * ~d) * cot(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ (~m + 1)) / (~a * ~f * (~m + 1) * ((~a) ^ 2 - (~b) ^ 2)) + (1 / (~a * (~m + 1) * ((~a) ^ 2 - (~b) ^ 2))) * Antiderivative((~a + ~b * csc(~e + ~f * ~x)) ^ (~m + 1) * Simp((~c * ((~a) ^ 2 - (~b) ^ 2) * (~m + 1) - (~a * (~b * ~c - ~a * ~d) * (~m + 1)) * csc(~e + ~f * ~x)) + ~b * (~b * ~c - ~a * ~d) * (~m + 2) * csc(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (LtQ(~m, -1) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && IntegerQ(2 * ~m))))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => ~c * Antiderivative((~a + ~b * csc(~e + ~f * ~x)) ^ ~m, ~x) + ~d * Antiderivative((~a + ~b * csc(~e + ~f * ~x)) ^ ~m * csc(~e + ~f * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && Not(IntegerQ(2 * ~m)))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) / (~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => (1 / ~c) * Antiderivative(sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) - (~d / ~c) * Antiderivative((csc(~e + ~f * ~x) * sqrt(~a + ~b * csc(~e + ~f * ~x))) / (~c + ~d * csc(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) || EqQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) / (~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => (~a / ~c) * Antiderivative(1 / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) + ((~b * ~c - ~a * ~d) / ~c) * Antiderivative(csc(~e + ~f * ~x) / (sqrt(~a + ~b * csc(~e + ~f * ~x)) * (~c + ~d * csc(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ (3 / 2) / (~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => (~a / ~c) * Antiderivative(sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) + ((~b * ~c - ~a * ~d) / ~c) * Antiderivative((csc(~e + ~f * ~x) * sqrt(~a + ~b * csc(~e + ~f * ~x))) / (~c + ~d * csc(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) || EqQ((~c) ^ 2 - (~d) ^ 2, 0)))
    #= Antiderivative(((~a)+(~b')*csc((~e')+(~f')*(~x)))^(3/2)/((~c)+(~d')*csc((~e')+(~f')*(~x))),(~x)*Symbol) := b/d*Antiderivative(sqrt(a+b*csc(e+f*x)),x) -  (b*c-a*d)/d*Antiderivative(sqrt(a+b*csc(e+f*x))/(c+d*csc(e+f*x)),x) <--  FreeQ([a,b,c,d,e,f],x) && NeQ(b*c-a*d,0) && NeQ(a^2-b^2,0) &&  NeQ(c^2-d^2,0) =#
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ (3 / 2) / (~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => (1 / (~c * ~d)) * Antiderivative(((~a) ^ 2 * ~d + (~b) ^ 2 * ~c * csc(~e + ~f * ~x)) / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) - ((~b * ~c - ~a * ~d) ^ 2 / (~c * ~d)) * Antiderivative(csc(~e + ~f * ~x) / (sqrt(~a + ~b * csc(~e + ~f * ~x)) * (~c + ~d * csc(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(1 / (sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) * (~c + ~(d') * csc(~(e') + ~(f') * ~x))), ~x) => (1 / (~c * (~b * ~c - ~a * ~d))) * Antiderivative(((~b * ~c - ~a * ~d) - ~b * ~d * csc(~e + ~f * ~x)) / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) + ((~d) ^ 2 / (~c * (~b * ~c - ~a * ~d))) * Antiderivative((csc(~e + ~f * ~x) * sqrt(~a + ~b * csc(~e + ~f * ~x))) / (~c + ~d * csc(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) || EqQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(1 / (sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) * (~c + ~(d') * csc(~(e') + ~(f') * ~x))), ~x) => (1 / ~c) * Antiderivative(1 / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) - (~d / ~c) * Antiderivative(csc(~e + ~f * ~x) / (sqrt(~a + ~b * csc(~e + ~f * ~x)) * (~c + ~d * csc(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && NeQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) * sqrt(~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => ((sqrt(~a + ~b * csc(~e + ~f * ~x)) * sqrt(~c + ~d * csc(~e + ~f * ~x))) / cot(~e + ~f * ~x)) * Antiderivative(cot(~e + ~f * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && EqQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) * sqrt(~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => ~c * Antiderivative(sqrt(~a + ~b * csc(~e + ~f * ~x)) / sqrt(~c + ~d * csc(~e + ~f * ~x)), ~x) + ~d * Antiderivative((csc(~e + ~f * ~x) * sqrt(~a + ~b * csc(~e + ~f * ~x))) / sqrt(~c + ~d * csc(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && NeQ(~b * ~c - ~a * ~d, 0)
    @apply_utils Antiderivative(sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) / sqrt(~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => (1 / ~c) * Antiderivative(sqrt(~a + ~b * csc(~e + ~f * ~x)) * sqrt(~c + ~d * csc(~e + ~f * ~x)), ~x) - (~d / ~c) * Antiderivative((csc(~e + ~f * ~x) * sqrt(~a + ~b * csc(~e + ~f * ~x))) / sqrt(~c + ~d * csc(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && EqQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) / sqrt(~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => ((-2 * ~a) / ~f) * Subst(Antiderivative(1 / (1 + ~a * ~c * (~x) ^ 2), ~x), ~x, cot(~e + ~f * ~x) / (sqrt(~a + ~b * csc(~e + ~f * ~x)) * sqrt(~c + ~d * csc(~e + ~f * ~x)))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) / sqrt(~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => (~a / ~c) * Antiderivative(sqrt(~c + ~d * csc(~e + ~f * ~x)) / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) + ((~b * ~c - ~a * ~d) / ~c) * Antiderivative(csc(~e + ~f * ~x) / (sqrt(~a + ~b * csc(~e + ~f * ~x)) * sqrt(~c + ~d * csc(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && EqQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) / sqrt(~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => ((2 * (~a + ~b * csc(~e + ~f * ~x))) / (~c * ~f * Rt((~a + ~b) / (~c + ~d), 2) * cot(~e + ~f * ~x))) * sqrt(((~b * ~c - ~a * ~d) * (1 + csc(~e + ~f * ~x))) / ((~c - ~d) * (~a + ~b * csc(~e + ~f * ~x)))) * sqrt((-((~b * ~c - ~a * ~d)) * (1 - csc(~e + ~f * ~x))) / ((~c + ~d) * (~a + ~b * csc(~e + ~f * ~x)))) * Elliptic.Pi((~a * (~c + ~d)) / (~c * (~a + ~b)), asin((Rt((~a + ~b) / (~c + ~d), 2) * sqrt(~c + ~d * csc(~e + ~f * ~x))) / sqrt(~a + ~b * csc(~e + ~f * ~x))), ((~a - ~b) * (~c + ~d)) / ((~a + ~b) * (~c - ~d))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(1 / (sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) * sqrt(~c + ~(d') * csc(~(e') + ~(f') * ~x))), ~x) => (cot(~e + ~f * ~x) / (sqrt(~a + ~b * csc(~e + ~f * ~x)) * sqrt(~c + ~d * csc(~e + ~f * ~x)))) * Antiderivative(1 / cot(~e + ~f * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && EqQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(1 / (sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) * sqrt(~c + ~(d') * csc(~(e') + ~(f') * ~x))), ~x) => (1 / ~a) * Antiderivative(sqrt(~a + ~b * csc(~e + ~f * ~x)) / sqrt(~c + ~d * csc(~e + ~f * ~x)), ~x) - (~b / ~a) * Antiderivative(csc(~e + ~f * ~x) / (sqrt(~a + ~b * csc(~e + ~f * ~x)) * sqrt(~c + ~d * csc(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && NeQ(~b * ~c - ~a * ~d, 0)
    @apply_utils Antiderivative(sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) / (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ (3 / 2), ~x) => (1 / ~c) * Antiderivative(sqrt(~a + ~b * csc(~e + ~f * ~x)) / sqrt(~c + ~d * csc(~e + ~f * ~x)), ~x) - (~d / ~c) * Antiderivative((csc(~e + ~f * ~x) * sqrt(~a + ~b * csc(~e + ~f * ~x))) / (~c + ~d * csc(~e + ~f * ~x)) ^ (3 / 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (((~a) ^ 2 * cot(~e + ~f * ~x)) / (~f * sqrt(~a + ~b * csc(~e + ~f * ~x)) * sqrt(~a - ~b * csc(~e + ~f * ~x)))) * Subst(Antiderivative(((~a + ~b * ~x) ^ (~m - 1 / 2) * (~c + ~d * ~x) ^ ~n) / (~x * sqrt(~a - ~b * ~x)), ~x), ~x, csc(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && IntegerQ(~m - 1 / 2))))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => Antiderivative(((~b + ~a * sin(~e + ~f * ~x)) ^ ~m * (~d + ~c * sin(~e + ~f * ~x)) ^ ~n) / sin(~e + ~f * ~x) ^ (~m + ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (IntegerQ(~m) && (IntegerQ(~n) && LeQ(-2, ~m + ~n, 0))))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((sqrt(~d + ~c * sin(~e + ~f * ~x)) * sqrt(~a + ~b * csc(~e + ~f * ~x))) / (sqrt(~b + ~a * sin(~e + ~f * ~x)) * sqrt(~c + ~d * csc(~e + ~f * ~x)))) * Antiderivative(((~b + ~a * sin(~e + ~f * ~x)) ^ ~m * (~d + ~c * sin(~e + ~f * ~x)) ^ ~n) / sin(~e + ~f * ~x) ^ (~m + ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (IntegerQ(~m + 1 / 2) && (IntegerQ(~n + 1 / 2) && LeQ(-2, ~m + ~n, 0))))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((sin(~e + ~f * ~x) ^ (~m + ~n) * (~a + ~b * csc(~e + ~f * ~x)) ^ ~m * (~c + ~d * csc(~e + ~f * ~x)) ^ ~n) / ((~b + ~a * sin(~e + ~f * ~x)) ^ ~m * (~d + ~c * sin(~e + ~f * ~x)) ^ ~n)) * Antiderivative(((~b + ~a * sin(~e + ~f * ~x)) ^ ~m * (~d + ~c * sin(~e + ~f * ~x)) ^ ~n) / sin(~e + ~f * ~x) ^ Simplify(~m + ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ(~m + ~n, 0) && Not(IntegerQ(2 * ~m))))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => Antiderivative(ExpandTrig((~a + ~b * csc(~e + ~f * ~x)) ^ ~m, (~c + ~d * csc(~e + ~f * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => Unintegrable((~a + ~b * csc(~e + ~f * ~x)) ^ ~m * (~c + ~d * csc(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * sec(~(e') + ~(f') * ~x)) ^ ~(m') * (~(d') / sec(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~d) ^ ~m * Antiderivative((~b + ~a * cos(~e + ~f * ~x)) ^ ~m * (~d * cos(~e + ~f * ~x)) ^ (~n - ~m), ~x) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~n], ~x) && (Not(IntegerQ(~n)) && IntegerQ(~m))
    @apply_utils Antiderivative((~(a') + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~(d') / csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~d) ^ ~m * Antiderivative((~b + ~a * sin(~e + ~f * ~x)) ^ ~m * (~d * sin(~e + ~f * ~x)) ^ (~n - ~m), ~x) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~n], ~x) && (Not(IntegerQ(~n)) && IntegerQ(~m))
    @apply_utils Antiderivative((~(a') + ~(b') * sec(~(e') + ~(f') * ~x)) ^ ~(m') * (~(c') * (~(d') * sec(~(e') + ~(f') * ~x)) ^ ~p) ^ ~n, ~x) => (((~c) ^ IntPart(~n) * (~c * (~d * sec(~e + ~f * ~x)) ^ ~p) ^ FracPart(~n)) / (~d * sec(~e + ~f * ~x)) ^ (~p * FracPart(~n))) * Antiderivative((~a + ~b * sec(~e + ~f * ~x)) ^ ~m * (~d * sec(~e + ~f * ~x)) ^ (~n * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && Not(IntegerQ(~n))
    @apply_utils Antiderivative((~(a') + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~(c') * (~(d') * csc(~(e') + ~(f') * ~x)) ^ ~p) ^ ~n, ~x) => (((~c) ^ IntPart(~n) * (~c * (~d * csc(~e + ~f * ~x)) ^ ~p) ^ FracPart(~n)) / (~d * csc(~e + ~f * ~x)) ^ (~p * FracPart(~n))) * Antiderivative((~a + ~b * cos(~e + ~f * ~x)) ^ ~m * (~d * cos(~e + ~f * ~x)) ^ (~n * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && Not(IntegerQ(~n))
end

max_stars_repo_path: src/rules/4 Trig functions/4.5 Secant/4.5.2.1 (a+b sec)^m (c+d sec)^n.jl


text:
const βc = log(1+sqrt(2))/2

βofh(h) = βc/asinh(1) * asinh(sqrt(1/h))
hofβ(β) = sinh(β * asinh(1)/βc)^(-2)

function magofβ(β)
    if β > βc
        (1-sinh(2*β)^-4)^(1/8)
    else
        0
    end
end

#Crone Thesis tensors:
function isingtensor(β)
    a = DTensor{Float64}((2,2,2,2))
    initwithzero!(a)
    a[1,1,1,1] = a[2,2,2,2] = 1
    cβ, sβ = sqrt(cosh(β)), sqrt(sinh(β))
    Q = DTensor(1/sqrt(2) * [cβ+sβ cβ-sβ; cβ-sβ cβ+sβ])
    @tensor a[1,2,3,4] := a[-1,-2,-3,-4] * Q[-1,1] * Q[-2,2] * Q[-3,3] * Q[-4,4]
    return a
end

function isingmagtensor(β)
    a = DTensor{Float64}((2,2,2,2))
    initwithzero!(a)
    a[1,1,1,1] = 1
    a[2,2,2,2] = -1
    cβ, sβ = sqrt(cosh(β)), sqrt(sinh(β))
    Q = DTensor(1/sqrt(2) * [cβ+sβ cβ-sβ; cβ-sβ cβ+sβ])
    @tensor a[1,2,3,4] := a[-1,-2,-3,-4] * Q[-1,1] * Q[-2,2] * Q[-3,3] * Q[-4,4]
    return a
end

function isingenvironment(β, χ; fixed::Bool = true)
    C, T = zeros(Float64,χ,χ), zeros(Float64, χ, 2, χ)
    C[1,1] = T[1,1,1] =  1
    !fixed && (C[2,2] = T[2,2,2] = 1)

    cβ, sβ = sqrt(cosh(β)), sqrt(sinh(β))
    q = 1/sqrt(2) * [cβ+sβ cβ-sβ; cβ-sβ cβ+sβ]
    Q = zeros(Float64, χ, χ)
    Q[1:2,1:2] = q
    @tensor T[1,2,3] := T[-1,-2,-3] * Q[1,-1] * q[2,-2] * Q[3,-3]
    @tensor C[1,2] := C[-1,-2] * Q[1,-1] * Q[2,-2]
    return (DTensor(C), DTensor(T))
end

function isingenvironmentz2(β::Float64, χ::Int)
    C = DASTensor{Float64,2}(Z2(),ntuple(i -> Z2Charges(), 2),
            ([χ,χ],[χ,χ]), InOut(-1,-1))
    T = DASTensor{Float64,3}(Z2(),ntuple(i -> Z2Charges(),3),
            ([χ,χ],[1,1],[χ,χ]),InOut(1,-1,1))
    initwithzero!(C)
    initwithzero!(T)
    foreach(k -> C[k][1] = 1, keys(C))
    foreach(k -> T[k][1] = 1, keys(T))

    q =  DASTensor{Float64,2}(Z2(),(Z2Charges(), Z2Charges()),
            ([1,1],[1,1]), InOut(1,-1))
    qχ = DASTensor{Float64,2}(Z2(),(Z2Charges(), Z2Charges()),
            ([χ,χ],[χ,χ]), InOut(1,-1))
    cβ, sβ = sqrt(2) * [sqrt(cosh(β))], sqrt(2) * [sqrt(sinh(β))]
    k00 = DASSector(Z2Charge(0),Z2Charge(0))
    k11 = DASSector(Z2Charge(1),Z2Charge(1))
    q[k00]  = reshape(cβ, 1,1)
    q[k11]  = reshape(sβ, 1,1)
    qχ[k00] = zeros(Float64,χ,χ)
    qχ[k11] = zeros(Float64,χ,χ)
    qχ[k00][1] = cβ[1]
    qχ[k11][1] = sβ[1]
    @tensor T[1,2,3] := T[-1,-2,-3] * qχ[1,-1] * q[2,-2] * qχ[3,-3]
    @tensor C[1,2] := C[-1,-2] * qχ[1,-1] * qχ[2,-2]
    return (C, T)
end

function isingtensorz2(β)
    a = DASTensor{Float64,4}(Z2(),
        ntuple(i -> Z2Charges(), 4),
        ntuple(i -> [1,1], 4),
        InOut(1,1,1,1))
    for k in invariantsectors(charges(a), in_out(a))
        a[k] = reshape([1],1,1,1,1)
    end
    b = DASTensor{Float64,2}(Z2(),
        ntuple(i -> Z2Charges(), 2),
        ntuple(i -> [1,1], 2),
        InOut(1,-1))
    cβ, sβ = sqrt(2) * [sqrt(cosh(β))], sqrt(2) * [sqrt(sinh(β))]
    b[DASSector(Z2Charge(0),Z2Charge(0))] = reshape(cβ,1,1)
    b[DASSector(Z2Charge(1),Z2Charge(1))] = reshape(sβ,1,1)
    @tensor a[o1,o2,o3,o4] := a[c1,c2,c3,c4] * b[o1,c1] * b[o2,c2] * b[o3,c3] * b[o4,c4]
    return a
end

isingmagtensorz2(β) = isingtensorz2(β)

isingtensors(β) = (isingtensor(β), isingmagtensor(β))
isingtensorsz2(β) = (isingtensorz2(β), isingmagtensorz2(β))

magnetisation(state, a, asz) = magnetisation(environment(state), a, asz)
magnetisation((C,T,a,asz)) = magnetisation(C,T,a,asz)

function magnetisation(env::AbstractTensor{<:Any,4}, a, asz)
    @tensor mag  = env[-1,-2,-3,-4] * asz[-1,-2,-3,-4]
    @tensor norm = env[-1,-2,-3,-4] * a[-1,-2,-3,-4]
    return mag/norm
end

function magnetisation(C::AbstractTensor{S,2}, T::AbstractTensor{S,3},
        a::AbstractTensor{S,4}, asz::AbstractTensor{S,4}) where S
    env = environment(ntuple(i -> C,4), ntuple(i -> T,4))
    return magnetisation(env, a, asz)
end

function isingctm(β, χ, fixed::Bool = true; log::Bool = false, kwargs...)
    cinit, tinit = isingenvironment(β, χ, fixed = fixed)
    a, asz = isingtensors(β)
    if log
        st = ctm(a, asz, χ, Cinit = cinit, Tinit = tinit; kwargs...)
        info = (β = β, χ = χ, n_it = st.n_it[], diffs = st.diffs)
        return (C = st.C, T = st.T, A = a, M = asz, ttype = :DTensor, info...)
    end
    st = ctm(a, asz, χ, Cinit = cinit, Tinit = tinit ; kwargs...)
    c, t = st.C, st.T
    return (C = st.C, T = st.T, A = a, M = asz)
end

function isingctmz2(β, χ; log::Bool = false, kwargs...)
    cinit, tinit = isingenvironmentz2(β, χ)
    a, asz = isingtensorsz2(β)
    if log
        st = ctm(a, asz, χ, Cinit = cinit, Tinit = tinit; kwargs...)
        info = (β = β, χ = χ, n_it = st.n_it[], diffs = st.diffs)
        return (C = st.C, T = st.T, A = a, M = asz, ttype = :ZNTensor, info...)
    end
    st = ctm(a, asz, χ, Cinit = cinit, Tinit = tinit; kwargs...)
    c, t = st.C, st.T
    return (C = c, T = t, A = a, M = asz)
end

function σxtensor()
    σx = DASTensor{Float64,2}(
            ZN{2}(),
            (Z2Charges(), Z2Charges()),
            ([1,1], [1,1]),
            InOut(1,-1))
    σx[DASSector(Z2Charge(1), Z2Charge(0))] = reshape([1], 1, 1)
    σx[DASSector(Z2Charge(0), Z2Charge(1))] = reshape([1], 1, 1)
    return σx
end

function σztensor()
    σz = DASTensor{Float64,2}(
            ZN{2}(),
            (Z2Charges(), Z2Charges()),
            ([1,1], [1,1]),
            InOut(1,-1))
    σz[DASSector(Z2Charge(0), Z2Charge(0))] = reshape([1], 1, 1)
    σz[DASSector(Z2Charge(1), Z2Charge(1))] = reshape([-1], 1, 1)
    return σz
end

environment(state::CTMState) = environment(state.Cs, state.Ts)
environment(state::rotsymCTMState) =
    environment(ntuple(i -> state.C,4), ntuple(i->state.T,4))
function environment(state::transconjCTMState)
    C = state.C
    T1, T2 = state.Ts
    environment((C,C',C,C'), (T1, T2', T1', T2))
end
function environment((C1,C2,C3,C4), (T1,T2,T3,T4))
    @tensor env[1,2,3,4] := C1[-8, -1] * T1[-1,3,-2] *  C2[-2,-3] *
                            T4[-7,4,-8] * #=A[1,2,3,4]=#T2[-3,2,-4] *
                            C4[-6,-7] *  T3[-5,1,-6] *  C3[-4,-5]
    return env
end


text:
InputFolder = './Images/CrazyCells/';
OutputFolder = './Results/Images/CrazyCells/';

@iA = '*_C00*.tif';		% Image filter
@iS = '*_C01*.tif';		% Image seeds filter

@fxgs_lSeededPolarGeod [iA, iS] > [L];
params.ThDiv = 180;
params.ObjRad = 35;
params.MinArea = 175;
params.Method = 'simple';
/endf

/show iA > iS;
/show iA > L;
/keep L > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
"""
    struct Plane{T} <: AbstractPlanarSurfacePrimitive{T}

Surface primitive describing a two-dimensional flat plane in three-dimensional space.

## Fields
* `origin::CartesianPoint{T}`: Point in the `Plane`.
* `normal::CartesianVector{T}`: Normal vector of the `Plane`, normalized to length 1.
"""
struct Plane{T} <: AbstractPlanarSurfacePrimitive{T}
    origin::CartesianPoint{T}
    normal::CartesianVector{T}
    Plane{T}(o, n) where {T} = new{T}(o, normalize(n))
end
Plane(origin::CartesianPoint{T}, normal::CartesianVector{T}) where {T} = Plane{T}(origin, normal)

normal(p::Plane) = p.normal
origin(p::Plane) = p.origin

isinfront(pt::AbstractCoordinatePoint, p::Plane) = (pt - origin(p)) ⋅ normal(p) > 0
isbehind(pt::AbstractCoordinatePoint, p::Plane) = (pt - origin(p)) ⋅ normal(p) < 0

_distance(pt::AbstractCoordinatePoint, p::Plane) = (pt - origin(p)) ⋅ normal(p)
distance(pt::AbstractCoordinatePoint, p::Plane) = abs(_distance(pt, p))

"""
    intersection(p::Plane{T}, line::Line{T}) where {T}

Calculates the intersections of a `Line` with a `Plane`.

## Arguments
* `cm::Plane{T}`: The `Plane`.
* `l::Line{T}`: The `Line`.

!!! note
    The function will always return one Point as a Tuple.
    If the line is parallel to the plane, the point will have `NaN`'s/`Inf`'s as values.
"""
function intersection(p::Plane{T}, line::Line{T}) where {T}
    ndir = normalize(line.direction)
    λ = (p.normal ⋅ p.origin - p.normal ⋅ line.origin) / (p.normal ⋅ ndir)
    (line.origin + λ * ndir,)
end

intersection(p::AbstractPlanarSurfacePrimitive{T}, line::Line{T}) where {T} = intersection(Plane(p), line)


text:
mutable struct quadrotor_input_t <: LCMType
    timestamp::Int64
    motors::SVector{4, Float64}
end

@lcmtypesetup(quadrotor_input_t)


text:
__precompile__(true)
module DataStreams

export Data, DataFrame

module Data

if !isdefined(Core, :String)
    typealias String UTF8String
end

# Data.Schema
"""
A `Data.Schema` describes a tabular dataset (i.e. a set of optionally named, typed columns with records as rows)
`Data.Schema` allow `Data.Source` and `Data.Sink` to talk to each other and prepare to provide/receive data through streaming.
`Data.Schema` fields include:

 * `Data.header(schema)` to return the header/column names in a `Data.Schema`
 * `Data.types(schema)` to return the column types in a `Data.Schema`; `Nullable{T}` indicates columns that may contain missing data (null values)
 * `Data.size(schema)` to return the (# of rows, # of columns) in a `Data.Schema`

`Data.Source` and `Data.Sink` interfaces both require that `Data.schema(source_or_sink)` be defined to ensure
that other `Data.Source`/`Data.Sink` can work appropriately.
"""
type Schema
    header::Vector{String}       # column names
    types::Vector{DataType}      # Julia types of columns
    rows::Int                    # number of rows in the dataset
    cols::Int                    # number of columns in a dataset
    metadata::Dict{Any, Any}     # for any other metadata we'd like to keep around (not used for '==' operation)
    function Schema(header::Vector, types::Vector{DataType}, rows::Integer=0, metadata::Dict=Dict())
        cols = length(header)
        cols != length(types) && throw(ArgumentError("length(header): $(length(header)) must == length(types): $(length(types))"))
        header = String[string(x) for x in header]
        return new(header, types, rows, cols, metadata)
    end
end

Schema(header, types::Vector{DataType}, rows::Integer=0, meta::Dict=Dict()) = Schema(String[i for i in header], types, rows, meta)
Schema(types::Vector{DataType}, rows::Integer=0, meta::Dict=Dict()) = Schema(String["Column$i" for i = 1:length(types)], types, rows, meta)
const EMPTYSCHEMA = Schema(String[], DataType[], 0, Dict())
Schema() = EMPTYSCHEMA

header(sch::Schema) = sch.header
types(sch::Schema) = sch.types
Base.size(sch::Schema) = (sch.rows, sch.cols)
Base.size(sch::Schema, i::Int) = ifelse(i == 1, sch.rows, ifelse(i == 2, sch.cols, 0))

function Base.show(io::IO, schema::Schema)
    println(io, "Data.Schema:")
    println(io, "rows: $(schema.rows)\tcols: $(schema.cols)")
    if schema.cols <= 0
        println(io)
    else
        println(io, "Columns:")
        Base.print_matrix(io, hcat(schema.header, schema.types))
    end
end

# Data.Source / Data.Sink
abstract Source

function reset! end
function isdone end
reference(x) = UInt8[]

function getfield end
function getcolumn end

abstract Sink

open!(sink, source) = nothing
function streamfield! end
streamfield!{T}(sink, source, ::Type{T}, row, col, cols, sinkrows) = streamfield!(sink, source, T, row, col, cols)
cleanup!(sink) = nothing
flush!(sink) = nothing
close!(sink) = nothing

function streamcolumn! end

abstract StreamType
immutable Field <: StreamType end
immutable Column <: StreamType end

"""
`Data.streamtype{T<:Data.Source, S<:Data.StreamType}(::Type{T}, ::Type{S})` => Bool

Indicates whether the source `T` supports streaming of type `S`. To be overloaded by individual sources according to supported `Data.StreamType`s
"""
function streamtype end

# generic fallback for all Sources
Data.streamtype{T<:StreamType}(source, ::Type{T}) = false

"""
`Data.streamtypes{T<:Data.Sink}(::Type{T})` => Vector{StreamType}

Returns a list of `Data.StreamType`s that the sink is able to receive; the order of elements indicates the sink's streaming preference
"""
function streamtypes end

schema(source_or_sink) = isdefined(source_or_sink, :schema) ? deepcopy(source_or_sink.schema) : throw(ArgumentError("unable to get schema of $source_or_sink"))
header(source_or_sink) = header(schema(source_or_sink))
types(source_or_sink) = types(schema(source_or_sink))
Base.size(source::Source) = size(schema(source))
Base.size(source::Source, i) = size(schema(source),i)
Base.size(sink::Sink) = size(schema(sink))
Base.size(sink::Sink, i) = size(schema(sink),i)
setrows!(source, rows) = isdefined(source, :schema) ? (source.schema.rows = rows; nothing) : nothing
setcols!(source, cols) = isdefined(source, :schema) ? (source.schema.cols = cols; nothing) : nothing

# Data.stream!
function stream! end

# generic definitions
function Data.stream!{T, TT}(source::T, ::Type{TT}, append::Bool, args...)
    typs = Data.streamtypes(TT)
    for typ in typs
        if Data.streamtype(T, typ)
            sink = TT(source, typ, append, args...)
            return Data.stream!(source, typ, sink)
        end
    end
    throw(ArgumentError("`source` doesn't support the supported streaming types of `sink`: $typs"))
end
# for backwards compatibility
Data.stream!{T, TT}(source::T, ::Type{TT}) = Data.stream!(source, TT, false, ())

function Data.stream!{T, TT}(source::T, sink::TT, append::Bool)
    typs = Data.streamtypes(TT)
    for typ in typs
        if Data.streamtype(T, typ)
            sink = TT(sink, source, typ, append)
            return Data.stream!(source, typ, sink)
        end
    end
    throw(ArgumentError("`source` doesn't support the supported streaming types of `sink`: $typs"))
end
Data.stream!{T, TT <: Data.Sink}(source::T, sink::TT) = Data.stream!(source, sink, false)

# Generic Data.stream! method for Data.Field
function Data.stream!{T1, T2}(source::T1, ::Type{Data.Field}, sink::T2)
    Data.types(source) == Data.types(sink) || throw(ArgumentError("schema mismatch: \n$(Data.schema(source))\nvs.\n$(Data.schema(sink))"))
    Data.isdone(source, 1, 1) && return sink
    rows, cols = size(source)
    sinkrows = max(0, size(sink, 1) - rows)
    types = Data.types(source)
    row = 1
    Data.open!(sink, source)
    try
        if rows == -1
            while true
                for col = 1:cols
                    @inbounds T = types[col]
                    Data.streamfield!(sink, source, T, row, col, cols)
                end
                row += 1
                Data.isdone(source, row, cols) && break
            end
            Data.setrows!(source, row - 1)
        else
            for row = 1:rows, col = 1:cols
                @inbounds T = types[col]
                Data.streamfield!(sink, source, T, row, col, cols, sinkrows)
            end
        end
    catch e
        Data.cleanup!(sink)
        rethrow(e)
    end
    Data.flush!(sink)
    return sink
end

# DataFrames DataStreams definitions
using DataFrames, NullableArrays, CategoricalArrays, WeakRefStrings

# because there's currently not a better place for this to live
import Base.==
=={T}(x::WeakRefString{T}, y::CategoricalArrays.CategoricalValue) = String(x) == String(y)
=={T}(y::CategoricalArrays.CategoricalValue, x::WeakRefString{T}) = String(x) == String(y)

# AbstractColumn definitions
nullcount(A::NullableVector) = sum(A.isnull)
nullcount(A::Vector) = 0
nullcount(A::CategoricalArray) = 0
nullcount(A::NullableCategoricalArray) = sum(A.refs .== 0)
if isdefined(Main, :DataArray)
    nullcount(A::DataArray) = sum(A.na)
end

allocate{T}(::Type{T}, rows, ref) = Array{T}(rows)
function allocate{T}(::Type{Nullable{T}}, rows, ref)
    A = Array{T}(rows)
    return NullableArray{T, 1}(A, fill(true, rows), isempty(ref) ? UInt8[] : ref)
end
allocate{S,R}(::Type{CategoricalArrays.CategoricalValue{S,R}}, rows, ref) = CategoricalArray{S,1,R}(rows)
allocate{S,R}(::Type{Nullable{CategoricalArrays.CategoricalValue{S,R}}}, rows, ref) = NullableCategoricalArray{S,1,R}(rows)

# DataFrames DataStreams implementation
function Data.schema(df::DataFrame)
    return Data.Schema(map(string, names(df)),
            DataType[eltype(A) for A in df.columns], size(df, 1))
end

# DataFrame as a Data.Source
function Data.isdone(source::DataFrame, row, col)
    rows, cols = size(source)
    return row > rows || col > cols
end

Data.streamtype(::Type{DataFrame}, ::Type{Data.Column}) = true
Data.streamtype(::Type{DataFrame}, ::Type{Data.Field}) = true

Data.getcolumn{T}(source::DataFrame, ::Type{T}, col) = (@inbounds A = source.columns[col]; return A)
Data.getcolumn{T}(source::DataFrame, ::Type{Nullable{T}}, col) = (@inbounds A = source.columns[col]::NullableVector{T}; return A)
Data.getcolumn{T,R}(source::DataFrame, ::Type{CategoricalValue{T,R}}, col) = (@inbounds A = source.columns[col]::CategoricalVector{T,R}; return A)
Data.getcolumn{T,R}(source::DataFrame, ::Type{Nullable{CategoricalValue{T,R}}}, col) = (@inbounds A = source.columns[col]::NullableCategoricalVector{T,R}; return A)

Data.getfield{T}(source::DataFrame, ::Type{T}, row, col) = (@inbounds A = Data.getcolumn(source, T, col); return A[row]::T)

# DataFrame as a Data.Sink
DataFrame{T<:Data.StreamType}(so, ::Type{T}, append::Bool, args...) = DataFrame(Data.schema(so), T, Data.reference(so))

function DataFrame{T<:Data.StreamType}(sch::Schema, ::Type{T}=Data.Field, ref::Vector{UInt8}=UInt8[])
    rows, cols = size(sch)
    rows = T === Data.Column || rows < 0 ? 0 : rows # don't pre-allocate for Column streaming
    columns = Vector{Any}(cols)
    types = Data.types(sch)
    for i = 1:cols
        columns[i] = allocate(types[i], rows, ref)
    end
    return DataFrame(columns, map(Symbol, Data.header(sch)))
end

# given an existing DataFrame (`sink`), make any necessary changes for streaming `source`
# to it, given we know if we'll be `appending` or not
function DataFrame{T<:Data.StreamType}(sink, source, ::Type{T}, append)
    sch = Data.schema(source)
    rows, cols = size(sch)
    if T === Data.Column
        if !append
            for col in sink.columns
                empty!(col)
            end
        end
    else
        if rows > 0
            newlen = rows + (append ? size(sink, 1) : 0)
            for col in sink.columns
                resize!(col, newlen)
            end
        end
    end
    return sink
end

Data.streamtypes(::Type{DataFrame}) = [Data.Column, Data.Field]

function Data.streamfield!{T}(sink::DataFrame, source, ::Type{T}, row, col, cols)
    push!(sink.columns[col], Data.getfield(source, T, row, col))
    return nothing
end

function streamfield{T}(dest, ::Type{T}, row, val)
    @inbounds dest[row] = val
    return nothing
end
function Data.streamfield!{T}(sink::DataFrame, source, ::Type{T}, row, col, cols, sinkrows)
    streamfield(sink.columns[col], T, sinkrows + row, Data.getfield(source, T, row, col))
    return nothing
end

function appendcolumn!{T}(source, ::Type{T}, dest, col)
    column = Data.getcolumn(source, T, col)
    append!(dest, column)
    return length(dest)
end

function appendcolumn!{T}(source, ::Type{Nullable{WeakRefString{T}}}, dest, col)
    column = Data.getcolumn(source, Nullable{WeakRefString{T}}, col)
    offset = length(dest.values)
    parentoffset = length(dest.parent)
    append!(dest.isnull, column.isnull)
    append!(dest.parent, column.parent)
    # appending new data to `dest` would invalid all existing WeakRefString pointers
    resize!(dest.values, length(dest) + length(column))
    for i = 1:offset
        old = dest.values[i]
        dest.values[i] = WeakRefString{T}(pointer(dest.parent, old.ind), old.len, old.ind)
    end
    for i = 1:length(column)
        old = column.values[i]
        dest.values[offset + i] = WeakRefString{T}(pointer(dest.parent, parentoffset + old.ind), old.len, parentoffset + old.ind)
    end
    return length(dest)
end

function Data.streamcolumn!{T}(sink::DataFrame, source, ::Type{T}, col, row)
    if row == 0
        sink.columns[col] = Data.getcolumn(source, T, col)
        len = length(sink.columns[col])
    else
        len = appendcolumn!(source, T, sink.columns[col], col)
    end
    return len
end

function Data.stream!{T}(source::T, ::Type{Data.Column}, sink)
    Data.types(source) == Data.types(sink) || throw(ArgumentError("\n\nschema mismatch: \n\nSOURCE: $(Data.schema(source))\n\nvs.\n\nSINK: $(Data.schema(sink))"))
    rows, cols = size(source)
    Data.isdone(source, 1, 1) && return sink
    types = Data.types(source)
    row = cur_row = 0
    for col = 1:cols
        row = Data.streamcolumn!(sink, source, types[col], col, row)
    end
    while !Data.isdone(source, row+1, cols)
        for col = 1:cols
            cur_row = Data.streamcolumn!(sink, source, types[col], col, row)
        end
        row += cur_row
    end
    Data.setrows!(source, row)
    return sink
end

end # module Data

end # module DataStreams


text:
function ffd(f::Function, x=0.0, h=1e-8)
    val = (f(x+h)-f(x))/h
    return val
end

function bfd(f::Function, x=0.0, h=1e-8)
    val = (f(x)-f(x-h))/h
    return val
end

function cfd(f::Function, x=0.0, h=1e-5)
    val = (f(x+h)-f(x-h))/2/h
    return val
end


text:
using Documenter
using Rimu

DocMeta.setdocmeta!(
    Rimu,
    :DocTestSetup,
    :(using Rimu; using Rimu.RMPI; using Rimu.StatsTools; using DataFrames; using Random;
      using LinearAlgebra);
    recursive=true
)
# Run with fix=true to fix docstrings
doctest(Rimu)


text:
###############################################################################
#
# Create IncomingRadiation
#
###############################################################################
"""
    create_incoming_radiation(
                wls::WaveLengths{FT},
                wlfn::String = FILE_SUN
    ) where {FT<:AbstractFloat}

Create an `AbstractIncomingRadiation` struct, given
- `wls` [`WaveLengths`](@ref) type struct
- `wlfn` File that saves incoming wave information
"""
function create_incoming_radiation(
            wls::WaveLengths{FT},
            wlfn::String = FILE_SUN
) where {FT<:AbstractFloat}
    @unpack sWL, nWL = wls;

    # Read data
    _suni  = matread(wlfn)["sun"]
    _wl    = _suni["wl"      ]
    _Edir  = _suni["Edirect" ]
    _Ediff = _suni["Ediffuse"]

    # create arrays
    # WL  = zeros(FT, N)
    Edir  = zeros(FT, nWL)
    Ediff = zeros(FT, nWL)

    # fill in the arrays
    # println("Reading Optical Parameters from ", sWL[1], " to ", sWL[end])
    for i in 1:nWL
        wo = findall( (_wl.>=sWL[i]) .& (_wl.<sWL[i+1]) )
        if length(wo)==0
            @warn "Some wavelengths out of bounds $(string(sWL[i]))";
        end
        # remove WL to avoid duplication
        # WL[i]  = mean(   _wl[wo])
        Edir[i]  = mean( _Edir[wo])
        Ediff[i] = mean(_Ediff[wo])
    end

    # create struct from the arrays
    return IncomingRadiation{FT}(Edir, Ediff)
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
using Base: Indices
using Test
using Plots
using FindPeaks

@testset "Tests" begin

    @testset "functions" begin
        include("test_functions.jl")
    end

    @testset "Peaks" begin
        include("test_PeakResults.jl")
    end

end

max_stars_repo_name: vkumpost/FindPeaks

max_issues_repo_name: vkumpost/FindPeaks

max_forks_repo_name: vkumpost/FindPeaks


text:
@testset "MerCount" begin
    v = DNAMer{4}[
        mer"AAAA", mer"AAAA", mer"AAAA", mer"ATAG",
        mer"GGGG", mer"GGGG", mer"GGGT", mer"AGGT",
        mer"AGGT", mer"AGGT", mer"AGGT"
    ]

    v2 = DNAMer{4}[
        mer"AAAA", mer"AAAA", mer"AAAA", mer"ATAG",
        mer"GGGG", mer"GGGG", mer"GGGT", mer"AGGT",
        mer"ACGT", mer"AAGT", mer"AAGT"
    ]

    function test_merge_into(a, b, a′)
        merge_into!(a, b)
        return isempty(b) && a == a′
    end

    @test collapse_into_counts(v) == [
        MerCount{DNAMer{4}}(mer"AAAA", 3),
        MerCount{DNAMer{4}}(mer"AGGT", 4),
        MerCount{DNAMer{4}}(mer"ATAG", 1),
        MerCount{DNAMer{4}}(mer"GGGG", 2),
        MerCount{DNAMer{4}}(mer"GGGT", 1)
    ]

    @test collapse_into_counts(v2) == [
        MerCount{DNAMer{4}}(mer"AAAA", 3),
        MerCount{DNAMer{4}}(mer"AAGT", 2),
        MerCount{DNAMer{4}}(mer"ACGT", 1),
        MerCount{DNAMer{4}}(mer"AGGT", 1),
        MerCount{DNAMer{4}}(mer"ATAG", 1),
        MerCount{DNAMer{4}}(mer"GGGG", 2),
        MerCount{DNAMer{4}}(mer"GGGT", 1)
    ]

    @test test_merge_into(collapse_into_counts(v), collapse_into_counts(v2), [
        MerCount{DNAMer{4}}(mer"AAAA", 6),
        MerCount{DNAMer{4}}(mer"AAGT", 2),
        MerCount{DNAMer{4}}(mer"ACGT", 1),
        MerCount{DNAMer{4}}(mer"AGGT", 5),
        MerCount{DNAMer{4}}(mer"ATAG", 2),
        MerCount{DNAMer{4}}(mer"GGGG", 4),
        MerCount{DNAMer{4}}(mer"GGGT", 2)
    ])
end


text:
"""
    serialize(val) :: BSON

Encodes a Julia value `val` into BSON format.

    serialize(io::IO, val)

Writes a Julia value `val` into `io` using BSON format.
"""
function serialize end

"""
    deserialize(doc::Union{BSON, Dict})

Decodes a BSON that was previously encoded using `BSONSerializer.serialize` method
into a Julia value.

This method also accepts a `Dict` as input, given that
you may want to decode just a portion of a BSON document.

Type information is encoded as a String in the BSON document,
so the `deserialize` method uses `m.eval` to convert that to a Julia type,
where `m` is the 2nd argument of the `deserialize` method.
"""
function deserialize end

function serialize(val::T) :: BSON where {T}
    return serialize(Serializable(val))
end

function serialize(val::Serializable{T}) where {T}
    error("Call @BSONSerializable($T) to generate serialize code for $T.")
end

# based on BSON.jl
resolve_typepath(fs::Vector) = foldl((m, f) -> getfield(m, Symbol(f)), fs; init = Main)

function deserialize(bson::Union{BSON, Dict})
    @assert haskey(bson, "type") && haskey(bson, "args")
    datatype = resolve_typepath(bson["type"])
    @assert isa(datatype, DataType)
    return deserialize(bson, Serializable{datatype})
end

"""
    roundtrip(val:::T) :: T

Applies `serialize` and `deserialize` to `val`.
Useful for testing.
"""
function roundtrip(val::T) where {T}
    return deserialize(serialize(Serializable(val)), Serializable{T})
end

function serialize(io::IO, val::T) where {T}
	write(io, serialize(val))
end

function deserialize(io::IO)
	return deserialize(read(io, BSON))
end

@BSONSerializable(Missing)
@BSONSerializable(Nothing)


text:
# Testing discrete univariate distributions

using Distributions
import JSON
using Base.Test

function verify_and_test_drive(jsonfile, n_tsamples::Int)
    R = JSON.parsefile(jsonfile)
    for (ex, dct) in R
        println("    testing $(ex)")

        # check type
        dtype = eval(symbol(dct["dtype"]))
        d = eval(parse(ex))
        if dtype == TruncatedNormal
            @test isa(d, Truncated{Normal})
        else
            @assert isa(dtype, Type) && dtype <: UnivariateDistribution
            @test isa(d, dtype)
        end

        # verification and testing
        verify_and_test(d, dct, n_tsamples)
    end
end

_parse_x(d::DiscreteUnivariateDistribution, x) = int(x)
_parse_x(d::ContinuousUnivariateDistribution, x) = float64(x)

function verify_and_test(d::UnivariateDistribution, dct::Dict, n_tsamples::Int)
    # verify parameters
    pdct = dct["params"]
    for (fname, val) in pdct
        f = eval(symbol(fname))
        @assert isa(f, Function)
        Base.Test.test_approx_eq(f(d), val, "$fname(d)", "val")
    end

    # verify stats
    r_min = float64(dct["minimum"])
    r_max = float64(dct["maximum"])
    @assert r_min < r_max

    @test_approx_eq minimum(d) r_min
    @test_approx_eq maximum(d) r_max
    @test_approx_eq_eps mean(d) float64(dct["mean"]) 1.0e-8
    @test_approx_eq_eps var(d) float64(dct["var"]) 1.0e-8
    @test_approx_eq_eps median(d) dct["median"] 1.0

    if applicable(entropy, d)
        @test_approx_eq_eps entropy(d) dct["entropy"] 1.0e-7
    end

    # verify quantiles
    @test_approx_eq_eps quantile(d, 0.10) dct["q10"] 1.0e-8
    @test_approx_eq_eps quantile(d, 0.25) dct["q25"] 1.0e-8
    @test_approx_eq_eps quantile(d, 0.50) dct["q50"] 1.0e-8
    @test_approx_eq_eps quantile(d, 0.75) dct["q75"] 1.0e-8
    @test_approx_eq_eps quantile(d, 0.90) dct["q90"] 1.0e-8

    # verify logpdf and cdf at certain points
    pts = dct["points"]
    for pt in pts
        x = _parse_x(d, pt["x"])
        lp = float64(pt["logpdf"])
        cf = float64(pt["cdf"])
        Base.Test.test_approx_eq(logpdf(d, x), lp, "logpdf(d, $x)", "lp")
        Base.Test.test_approx_eq(cdf(d, x), cf, "cdf(d, $x)", "cf")
    end

    # generic testing
    test_distr(d, n_tsamples)
end

## main

for c in ["discrete",
          "continuous"]

    title = string(uppercase(c[1]), c[2:end])
    println("    [$title]")
    println("    ------------")
    jsonfile = joinpath(dirname(@__FILE__), "$(c)_test.json")
    verify_and_test_drive(jsonfile, 10^6)
    println()
end


text:
module StochasticGroundMotionSimulation

using Interpolations
using Roots
using ForwardDiff
using ForwardDiff: Dual
using QuadGK
using FastGaussQuadrature
using LinearAlgebra

export Oscillator,
	FourierParameters,
	SourceParameters,
	GeometricSpreadingParameters,
	NearSourceSaturationParameters,
	AnelasticAttenuationParameters,
	PathParameters,
	SiteParameters,
	RandomVibrationParameters,
    period,
    transfer,
    transfer!,
    squared_transfer,
    squared_transfer!,
    site_amplification,
	kappa_filter,
    corner_frequency,
	geometric_spreading,
	near_source_saturation,
	equivalent_point_source_distance,
	anelastic_attenuation,
	fourier_constant,
	fourier_source,
	fourier_source_shape,
	fourier_path,
	fourier_attenuation,
	fourier_site,
	combined_kappa_frequency,
	fourier_spectral_ordinate,
	fourier_spectrum,
	fourier_spectrum!,
	spectral_moment,
	spectral_moments,
	excitation_duration,
	rms_duration,
	peak_factor,
	rvt_response_spectral_ordinate,
	rvt_response_spectrum,
	rvt_response_spectrum!

# Write your package code here.
include("oscillator/PJSoscillator.jl")
include("fourier/PJSfourierParameters.jl")
include("rvt/PJSrandomVibrationParameters.jl")
include("fourier/PJSsite.jl")
include("fourier/PJSsource.jl")
include("fourier/PJSpath.jl")
include("fourier/PJSfourierSpectrum.jl")
include("duration/PJSduration.jl")
include("rvt/PJSintegration.jl")
include("rvt/PJSpeakFactor.jl")
include("rvt/PJSrandomVibration.jl")

end


text:
using InverseGaussTransform

# ==========================
# TESTS
# ==========================

function trap_rule_test()
    Narr = [20,40,80,160,320,320*2,320*4,320*8,320*16,320*32]
    exact = 4^3/3
    errorArr = zeros(length(Narr))
    for ti=1:length(Narr)
        xarr = collect(range(0,stop=4.0,length=Narr[ti]))
        farr = xarr.^2
        quad = trap_rule(xarr,farr)
        errorArr[ti] = abs((quad - exact)/exact)
    end

    for tj=2:length(Narr)
        println("ratio: ",errorArr[tj-1]/errorArr[tj])
    end
    nothing
end

function cost_test()
    x1 = 0.7
    s = 1.2
    h = 0.5

    N = 200
    y = collect(range(-2,stop=2,length=N))
    F = zeros(N)
    for ti=1:N
        F[ti] = s*GaussKernel(y[ti] - x1,h)
    end

    # initial conditions [x0,s0,h0]
    #θ = [x1,s,h]
    #θ = [-1.0,1.0,0.5,0.2,0.2]
    θ = [-1.0,0.5,0.1]

    plot_fields(θ,y,F)
    J = cost(θ,y,F)

    nothing
end

function grad_cost_test()
    x1 = 0.7
    s = 1.2
    h = 0.5

    N = 200
    y = collect(range(-2,stop=2,length=N))
    F = zeros(N)
    for ti=1:N
        F[ti] = s*GaussKernel(y[ti] - x1,h)
    end

    # initial conditions [x0,s0,h0]
    #θ = [x1,s,h]
    #θ = [-1.0,1.0,0.5,0.2,0.2]
    θ = [0.0,0.0,0.1]

    plot_fields(θ,y,F)
    ∇J = grad_cost_FD(θ,y,F)
    println(∇J)

    nothing
end

function gradient_descent_test()
    #x1 = [0.7,-0.2]
    #s = [1.2,0.4]
    #x1 = 0.7
    #s  = 1.2
    #h = 0.5
    M = 100
    x1 = 2*rand(M) .- 1.0
    s1  = rand(M)
    h = 0.5

    N = 200
    y = collect(range(-2,stop=2,length=N))
    F = zeros(N)
    for ti=1:N
        for tj = 1:M
            F[ti] += s[tj]*GaussKernel(y[ti] - x1[tj],h)
        end
    end

    # initial conditions [x0,s0,h0]
    #θ = [x1,s,h]
    θ0 = [-0.3,0.0,0.1]
    #θ0 = [0.0,0.0,0.1]
    #θ0 = #[-1.0,1.0,0.5,0.2,0.2]
    #θ0 = [0.0,0.5,0.1]

    θ = gradient_descent(θ0,y,F;NiterMAX=1000,TOL=1e-8,pause=true)

    println(θ)
end

function gradient_descent_test_print_for_dissertation()
    M = 5
    x1 = 2*rand(M) .- 1.0
    s1  = rand(M)
    h = 0.5

    N = 200
    y = collect(range(-3,stop=3,length=N))
    F = zeros(N)
    for ti=1:N
        for tj = 1:M
            F[ti] += s1[tj]*GaussKernel(y[ti] - x1[tj],h)
        end
    end

    #θ0 = [-0.3,0.0,0.1]
    K= M+20
    θ0 = vcat(collect(range(-1,stop=1,length=K)),0.01*ones(K),0.1)

    error = 100.0
    α = 0.04

    p = plot(y,F,label="true",xlabel="y",xaxis=[-2,2])

    iter = 1
    θ = θ0
    Jtemp = 100.0
    NiterMAX = 2000
    TOL = 1e-8
    while iter<NiterMAX && error > TOL
        println("\riteration ",iter,"          ")
        plot_fields(θ,y,F;iter=iter)

        #readline()

        J = cost(θ,y,F)
        error = abs(J - Jtemp)

        ∇J = grad_cost(θ,y,F)

        gradient_descent_step!(θ,∇J,α)

        G = get_gauss_field(θ,y)
        if mod(iter,40) == 0
            plot!(p,y,G,linestyle=:dash,label="iter $(iter)")
        end

        Jtemp = J
        iter += 1
    end

   savefig(p,"test.png")

    println(θ)
end

function gradient_descent_test_print_for_dissertation_COST()
    x1 = 0.7
    s1  = 1.2
    h = 0.5

    N = 200
    y = collect(range(-3,stop=3,length=N))
    F = zeros(N)
    for ti=1:N
        F[ti] += s1*GaussKernel(y[ti] - x1,h)
    end

    θ0 = [-0.3,0.0,0.1]

    error = 100.0
    α = 0.2

    Jarr = [cost(θ0,y,F)]

    iter = 1
    θ = θ0
    Jtemp = 100.0
    NiterMAX = 5000
    TOL = 1e-8
    while iter<NiterMAX && error > TOL
        println("\riteration ",iter,"          ")
        #plot_fields(θ,y,F;iter=iter)

        J = cost(θ,y,F)
        push!(Jarr,J)
        error = abs(J - Jtemp)

        ∇J = grad_cost(θ,y,F)

        gradient_descent_step!(θ,∇J,α)

        G = get_gauss_field(θ,y)

        Jtemp = J
        iter += 1
    end

    p = plot(0:(length(Jarr)-1),Jarr,ylabel="J (cost)",label="",xlabel="iteration #")

    savefig(p,"cost.png")

    println(θ)
end

function dissertation_STEP()
    N_pts = 200
    y = collect(range(-3,stop=3,length=N_pts))
    F = zeros(N_pts)
    for ti=1:N_pts
        F[ti] = step_function(y[ti])
    end

    iterArr = Array{Int}(undef,0)
    indArr  = Array{Int}(undef,0)
    minJArr = Array{Float64}(undef,0)

    for tj=19:20
        println("\r calculating N=$(tj)        ")
        push!(indArr,tj)
        θ0 = vcat(collect(range(-1,stop=1,length=tj)),0.01*ones(tj),0.1)

        error = 100.0
        α = 0.04

        Jarr = [cost(θ0,y,F)]
        J = zero(Float64)

        iter = 1
        θ = θ0
        Jtemp = 100.0
        NiterMAX = 10000
        TOL = 1e-16
        while iter<NiterMAX && error > TOL
            #println("\riteration ",iter,"          ")
            #plot_fields(θ,y,F;iter=iter)

            J = cost(θ,y,F)
            #println(J)
            #push!(Jarr,J)
            error = abs(J - Jtemp)

            ∇J = grad_cost(θ,y,F)

            gradient_descent_step!(θ,∇J,α)

            G = get_gauss_field(θ,y)

            Jtemp = J
            iter += 1
        end
        plot_fields(θ,y,F)
        plot(p,"test.png")
        readline()
        push!(iterArr,iter-1)
        push!(minJArr,J)
    end

    #p1 = plot(indArr,minJArr,ylabel="minimum J",label="",xlabel="M")
    #savefig(p1,"cost_multiN_J.png")
    #p2 = plot(indArr,iterArr,ylabel="iterations to convergence",label="",xlabel="M")
    #savefig(p2,"cost_multiN_iter.png")
    nothing
end

function step_function(x)
    if x > -0.5 && x < 0.5
        return 1.0
    else
        return 0.0
    end
end


text:
mexpr(x::WExpr) = x

mexpr(x::Float32) = x
mexpr(x::Float64) = x
mexpr(x::Int32) = x
mexpr(x::Int64) = x

mexpr(x::Rational) = WExpr(:Rational, mexpr(x.num), mexpr(x.den))
mexpr(x::Complex) = WExpr(:Complex, mexpr(real(x)), mexpr(imag(x)))

mexpr(x::String) = x

mexpr(x::Bool) = x ? WSymbol(:True) : WSymbol(:False)
mexpr(x::Nothing) = WSymbol(:Nothing)

macro mdef(jfun, mfun)
    quote
        MathLink.mexpr(::typeof($(esc(jfun)))) = MSymbol($(QuoteNode(mfun)))
    end
end

@mdef(+, Plus)
@mdef(-, Subtract)
@mdef(*, Times)
@mdef(/, Divide)
@mdef(^, Power)

@mdef log Log
@mdef exp Exp
@mdef sin Sin
@mdef cos Cos
@mdef tan Tan
@mdef mod Mod

@mdef pi Pi
@mdef MathConstants.e E

macro mexpr(ex)
    :(mexpr($ex))
end

macro mexpr(ex::Symbol)
    :($(Expr(:isdefined, esc(ex))) ? mexpr($(esc(ex))) : MSymbol($(QuoteNode(ex))))
end

macro mexpr(ex::Expr)
    if ex.head == :call
        :(MFunc(@mexpr($(ex.args[1])),
                Any[$([:(@mexpr($arg)) for arg in ex.args[2:end]]...)]))
    elseif ex.head == :block
        :(MFunc(MSymbol(:CompoundExpression),
                Any[$([:(@mexpr($arg)) for arg in ex.args if !isa(arg, LineNumberNode)]...)]))
    elseif ex.head == :ref
        :(MFunc(MSymbol(:Part),
                Any[$([:(@mexpr($arg)) for arg in ex.args]...)]))
    else
        ex
    end
end

#=

# Conversion

from_mma(x) = x
const symbols = Dict(:True => true, :False => false, :Null => nothing)
from_mma(s::Symbol) = haskey(symbols, s) ? symbols[s] : s

to_mma(x::T) where {T<:Union{Int64,Int32,Float64,Float32,Symbol,AbstractString}}= x

function to_mma(x::Expr)
  if x.head == :call
    head = haskey(aliases, x.args[1]) ? aliases[x.args[1]] : x.args[1]
    MExpr{head}(map(to_mma, x.args[2:end]))
  elseif x.head == :block
    MExpr{:CompoundExpression}(map(to_mma, x.args))
  elseif x.head == :ref
    MExpr{:Part}(map(to_mma, x.args))
  elseif x.head == :cell1d
    MExpr{:List}(x.args)
  else
    error("Unsupported $(x.head) expression.")
  end
end

function to_mma(x::QuoteNode)
  typeof(x.value) == Symbol && return x.value
  error("Cannot call to_mma on QuoteNode($(x.value))")
end

# Other types

to_mma(x::Bool) = x ? :True : :False

to_mma(x::MExpr{T}) where {T} = MExpr{T}(map(to_mma, x.args))
from_mma(x::MExpr{T}) where {T} = MExpr{T}(map(from_mma, x.args))

to_mma(x::Rational) = MExpr(:Rational, x.num, x.den)
from_mma(f::MExpr{:Rational}) = f.args[1]//f.args[2]

to_mma(x::Complex) = MExpr(:Complex, to_mma(real(x)), to_mma(imag(x)))
from_mma(x::MExpr{:Complex}) = Complex(x.args[1], x.args[2])

to_mma(xs::Vector) = MExpr{:List}(map(to_mma, xs))
from_mma(l::MExpr{:List}) = map(from_mma, l.args)

# Julia Expression Conversion

to_expr(x) = x
to_expr(x::MExpr{T}) where {T} = Expr(:call, T, map(to_expr, x.args)...)

for (j, m) in aliases
  j = Expr(:quote, j)
  m = Expr(:quote, m)
  @eval to_expr(x::MExpr{$m}) = Expr(:call, $j, map(to_expr,x.args)...)
end

import Base.convert
convert(::Type{Expr}, x::MExpr) = to_expr(from_mma(x))
convert(::Type{MExpr}, x::Expr) = to_mma(x)
=#


text:
struct TimeoutException
    duration
end

function Base.showerror(io::IO, te::TimeoutException)
    print(io, "TimeoutException: Operation did not finish in ", te.duration)

    if !isa(te.duration, Period)
        print(io, " seconds")
    end
end

function asynctimedwait(fn, secs; kill=false)
    t = @async fn()
    timedwait(() -> istaskdone(t), secs)

    if istaskdone(t)
        fetch(t)
        return true
    else
        if kill
            Base.throwto(t, TimeoutException(secs))
        end

        return false
    end
end


text:
using WindMixing: animate_LES_3D
using Images, FileIO
using ImageTransformations
using CairoMakie

u_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\u.png"))))
v_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\v.png"))))
T_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\T.png"))))
uw_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\uw.png"))))
vw_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\vw.png"))))
wT_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\wT.png"))))
Ri_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\Ri.png"))))
z_img = load(assetpath(joinpath(pwd(), "images_plot\\z.png")))
T_3D_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\T_3D.png"))))

axis_images = (
    u = u_img,
    v = v_img,
    T = T_img,
    z = z_img,
    uw = uw_img,
    vw = vw_img,
    wT = wT_img,
    Ri = Ri_img,
    T_3D = T_3D_img,
)

FILE_DIR = "C:\\Users\\xinle\\Downloads\\three_layer_constant_fluxes_linear_hr144_Qu5.5e-04_Qb3.0e-08_f1.0e-04_Nh256_Nz128_strong_wind_weak_cooling"
OUTPUT_PATH = "test_LES"
animate_LES_3D(FILE_DIR, OUTPUT_PATH, axis_images, fps=5)