text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Generic function to inverse rotations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export inv_rotation

################################################################################
#                              Inverse Rotations
################################################################################

"""
    inv_rotation(R)

Compute the inverse rotation of `R`, which can be:

- A direction cosine matrix (`DCM`);
- An Euler angle and axis (`EulerAngleAxis`);
- A set of Euler anlges (`EulerAngles`); or
- A quaternion (`Quaternion`).

The output will have the same type as `R`.

!!! note
    If `R` is a DCM, than its transpose is computed instead of its inverse to
    reduce the computational burden. The both are equal if the DCM has unit
    norm. This must be verified by the user.

!!! note
    If `R` is a quaternion, than its conjugate is computed instead of its
    inverse to reduce the computational burden. The both are equal if the
    quaternion has unit norm. This must be verified by the used.

# Example

```jldoctest
julia> D = angle_to_dcm(pi / 3, pi / 4, pi / 5, :ZYX);

julia> inv_rotation(D)
DCM{Float64}:
  0.353553  -0.492816  0.795068
  0.612372   0.764452  0.201527
 -0.707107   0.415627  0.572061

julia> ea = EulerAngleAxis(30 * pi / 180, [1, 0, 0]);

julia> inv_rotation(ea)
EulerAngleAxis{Float64}:
  Euler angle : 0.523599 rad  (30.0°)
  Euler axis  : [-1.0, -0.0, -0.0]

julia> Θ = EulerAngles(-pi / 3, -pi / 2, -pi, :YXZ);

julia> inv_rotation(Θ)
EulerAngles{Float64}:
  R(Z) :  3.14159 rad  ( 180.0°)
  R(X) :  1.5708  rad  ( 90.0°)
  R(Y) :  1.0472  rad  ( 60.0°)

julia> q = angle_to_quat(pi / 3, pi / 4, pi / 5, :ZYX);

julia> inv_rotation(q)
Quaternion{Float64}:
  + 0.820071 - 0.0652687⋅i - 0.45794⋅j - 0.336918⋅k
```
"""
@inline inv_rotation(D::DCM) = D'
@inline inv_rotation(ea::EulerAngleAxis) = inv(ea)
@inline inv_rotation(Θ::EulerAngles) = inv(Θ)
@inline inv_rotation(q::Quaternion) = conj(q)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule hidapi_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("hidapi")
JLLWrappers.@generate_main_file("hidapi", UUID("fc16ec7e-3cba-5ca4-a17a-1acd8165bbee"))
end  # module hidapi_jll

ext: jl


text:
using FileIO
using ImageIO
using Images: channelview, RGB
using LinearAlgebra
using ArgCheck
using Random
using JSON3
using JLD2
using ProgressLogging
using ResumableFunctions

using JSON3.StructTypes: StructType, Struct

const Vec3 = NTuple{3,Float32}

Base.@kwdef struct CameraView
    camera_direction::Vec3
    camera_origin::Vec3
    x_axis::Vec3
    y_axis::Vec3
    x_fov::Float32
    y_fov::Float32
end

function CameraView(path)
    data = open(path, "r") do io
        JSON3.read(io)
    end
    CameraView(
        camera_direction=Vec3(data["z"]),
        camera_origin=Vec3(data["origin"]),
        x_axis=Vec3(data["x"]),
        y_axis=Vec3(data["y"]),
        x_fov=Float32(data["x_fov"]),
        y_fov=Float32(data["y_fov"])
    )
end

"""
Get all of the rays in the view in raster scan order

Returns 3x2xB array of (origin, direction) pairs
"""
function bare_rays(cv::CameraView, width::Int, height::Int)
    z = collect(cv.camera_direction)
    @check size(z) == (3,)
    rgh = collect(LinRange(-1, 1, height))
    rgh = reshape(rgh, 1, 1, size(rgh)...)
    ys = tan(cv.y_fov / 2) .* rgh .* collect(cv.y_axis)
    @check size(ys) == (3, 1, height)
    rgw = collect(LinRange(-1, 1, width))
    rgw = reshape(rgw, 1, size(rgw)..., 1)
    xs = tan(cv.x_fov / 2) .* rgw .* collect(cv.x_axis)
    @check size(xs) == (3, width, 1)
    directions = reshape(xs .+ ys .+ z, 3, 1, :)
    directions ./= map(norm, eachslice(directions; dims=1))
    co = collect(cv.camera_origin)
    co = reshape(co, size(co)..., 1, 1)
    origins = reshape(repeat(co, 1, width, height), 3, 1, :)
    @check size(origins) == size(directions) == (3, 1, height*width)
    return Float32.(hcat(origins, directions))
end

struct FileNeRF
    imagepath::String
    cameraview::CameraView
end

function FileNeRF(imagepath)
    camerapath = first(splitext(imagepath)) * ".json"
    cv = CameraView(camerapath)
    return FileNeRF(imagepath, cv)
end

function image(fn::FileNeRF)
    return load(fn.imagepath)
end

function rays(fn::FileNeRF)
    # https://discourse.julialang.org/t/reading-png-rgb-channels-julia-vs-python/73599/12
    # Note: this is not the same as PIL.Image something about linear vs SRGB
    img = channelview(RGB.(image(fn)))
    bare = bare_rays(fn.cameraview, size(img)[begin+1:end]...) # TODO
    colors = Float32.(reshape(img, 3, 1, :))
    return hcat(bare, colors)
end

struct ModelMetadata
    bbox_min::Vec3
    bbox_max::Vec3
end

function ModelMetadata(path::String)
    metadata = open(path, "r") do io
        JSON3.read(io)
    end
    ModelMetadata(Vec3(metadata["min"]), Vec3(metadata["max"]))
end

struct NeRFDataset
    metadata::ModelMetadata
    views::Vector{FileNeRF}
end

function NeRFDataset(path::AbstractString)
    metadata = ModelMetadata(joinpath(path, "metadata.json"))
    img_paths = filter!(endswith("png"), readdir(path; join=true))
    views = map(FileNeRF, img_paths)
    return NeRFDataset(metadata, views)
end

Base.length(d::NeRFDataset) = length(d.views)
Base.getindex(d::NeRFDataset, idx::Integer) = rays(d.views[idx])

function iterate_batches(ds::NeRFDataset, savepath::AbstractString; rng::AbstractRNG=Random.GLOBAL_RNG, num_shards::Int=32, batch_size::Int=512, repeat::Bool=true)
    sd = ShardedNeRFDataset(ds, savepath; rng, num_shards)
    return iterate_batches(sd; batch_size, rng, repeat)
end

struct ShardedNeRFDataset{S<:AbstractString}
    files::Vector{S}
end

function ShardedNeRFDataset(path::AbstractString)
    files = filter!(endswith("jld2"), readdir(path; join=true))
    @check length(files) > 0
    return ShardedNeRFDataset(files)
end

Base.length(ds::ShardedNeRFDataset) = length(ds.files)

function ShardedNeRFDataset(ds::NeRFDataset, savepath::AbstractString; rng::AbstractRNG=Random.GLOBAL_RNG, num_shards::Int=32)
    if !isdir(savepath)
        mkpath(savepath)
    else
        files = filter!(endswith("jld2"), readdir(savepath; join=true))
        if length(files) == num_shards
            return ShardedNeRFDataset(files)
        end
    end
    shards = Vector{Array{Float32,3}}(undef, num_shards)
    @progress for view in ds.views
        rs = rays(view)
        assignments = rand(rng, 1:num_shards, (size(rs)[end],))
        for sid in 1:num_shards
            sub_batch = rs[:, :, assignments .== sid]
            if size(sub_batch, 3) > 0
                try
                    shards[sid] = cat(shards[sid], sub_batch; dims=3)
                catch e
                    if e isa UndefRefError
                        shards[sid] = sub_batch
                    end
                end
            end
        end
    end
    files = [joinpath(savepath, "$i.jld2") for i in 1:num_shards]
    for i in 1:num_shards
        jldopen(files[i], "w") do f
            f["rays"] = shards[i]
        end
    end
    return ShardedNeRFDataset(files)
end

@resumable function iterate_batches(ds::ShardedNeRFDataset; batch_size::Int, rng::AbstractRNG=Random.GLOBAL_RNG, repeat::Bool=false)
    shard_indices = randperm(rng, length(ds))
    current_batch = nothing
    while true
        for shard in shard_indices
            shard_rays = JLD2.load(ds.files[shard], "rays")
            if current_batch === nothing
                current_batch = shard_rays
            else
                current_batch = cat(current_batch, shard_rays; dims=3)
            end
            while size(current_batch)[end] >= batch_size
                @yield current_batch[:, :, 1:batch_size]
                current_batch = current_batch[:, :, batch_size+1:end]
            end
        end
        if !repeat
            break
        end
    end
    if size(current_batch)[end] > 0
        @yield current_batch
    end
end


text:
function magic_square(n, ::Val{:JuMP})
    N = n^2
    model = JuMP.Model(CBLS.Optimizer)
    magic_constant = n * (N + 1) / 2

    @variable(model, 1 ≤ X[1:n, 1:n] ≤ N, Int)
    @constraint(model, vec(X) in AllDifferent())

    for i in 1:n
        @constraint(model, X[i,:] in SumEqualParam(magic_constant))
        @constraint(model, X[:,i] in SumEqualParam(magic_constant))
    end
    @constraint(model, [X[i,i] for i in 1:n] in SumEqualParam(magic_constant))
    @constraint(model, [X[i,n + 1 - i] for i in 1:n] in SumEqualParam(magic_constant))

    return model, X
end

"""
    magic_square(n; modeler = :JuMP)

Create a model for the magic square problem of order `n`. The `modeler` argument accepts :JuMP (default), which refer to the solver the JuMP model.
"""
magic_square(n; modeler = :JuMP) = magic_square(n, Val(modeler))


text:
# Check ChangepointSampler works for different distributions

using Distributions

num_samples = 1000
λ = 50                 # Frequency of changes

println("Running ChangepointSampler tests...")

# Normal mean change sampler
mu = Normal(1.0, 1.0)
σ = 1.0
Y = ChangepointSampler(()->Normal(rand(mu), σ), λ)
rand(Y, num_samples)
@changepoint_sampler num_samples λ Normal(mu, σ)

# Normal variance change sampler
μ = 1.0
σ = Uniform(2.0, 15.0)
Y = ChangepointSampler(()->Normal(μ, rand(σ)), λ)
rand(Y, num_samples)
@changepoint_sampler num_samples λ Normal(μ, σ)

# Exponential change sampler
μ = Uniform(0.0, 10.0)
Y = ChangepointSampler(()->Exponential(rand(μ)), λ)
rand(Y, num_samples)
@changepoint_sampler num_samples λ Exponential(μ)

# Poisson change sampler
μ = Uniform(0.0, 10.0)
Y = ChangepointSampler(()->Poisson(rand(μ)), λ)
rand(Y, num_samples)
@changepoint_sampler num_samples λ Poisson(μ)


text:
using CEnum

struct libfive_interval
    lower::Cfloat
    upper::Cfloat
end

struct libfive_region2
    X::libfive_interval
    Y::libfive_interval
end

struct libfive_region3
    X::libfive_interval
    Y::libfive_interval
    Z::libfive_interval
end

struct libfive_vec2
    x::Cfloat
    y::Cfloat
end

struct libfive_vec3
    x::Cfloat
    y::Cfloat
    z::Cfloat
end

struct libfive_vec4
    x::Cfloat
    y::Cfloat
    z::Cfloat
    w::Cfloat
end

struct libfive_tri
    a::UInt32
    b::UInt32
    c::UInt32
end

struct libfive_contour
    pts::Ptr{libfive_vec2}
    count::UInt32
end

struct libfive_contours
    cs::Ptr{libfive_contour}
    count::UInt32
end

struct libfive_contour3
    pts::Ptr{libfive_vec3}
    count::UInt32
end

struct libfive_contours3
    cs::Ptr{libfive_contour3}
    count::UInt32
end

struct libfive_mesh
    verts::Ptr{libfive_vec3}
    tris::Ptr{libfive_tri}
    tri_count::UInt32
    vert_count::UInt32
end

struct libfive_mesh_coords
    verts::Ptr{libfive_vec3}
    vert_count::UInt32
    coord_indices::Ptr{Int32}
    coord_index_count::UInt32
end

struct libfive_pixels
    pixels::Ptr{Bool}
    width::UInt32
    height::UInt32
end

function libfive_contours_delete(cs)
    ccall((:libfive_contours_delete, libfive), Cvoid, (Ptr{libfive_contours},), cs)
end

function libfive_contours3_delete(cs)
    ccall((:libfive_contours3_delete, libfive), Cvoid, (Ptr{libfive_contours3},), cs)
end

function libfive_mesh_delete(m)
    ccall((:libfive_mesh_delete, libfive), Cvoid, (Ptr{libfive_mesh},), m)
end

function libfive_mesh_coords_delete(m)
    ccall((:libfive_mesh_coords_delete, libfive), Cvoid, (Ptr{libfive_mesh_coords},), m)
end

function libfive_pixels_delete(ps)
    ccall((:libfive_pixels_delete, libfive), Cvoid, (Ptr{libfive_pixels},), ps)
end

function libfive_opcode_enum(op)
    ccall((:libfive_opcode_enum, libfive), Cint, (Ptr{Cchar},), op)
end

function libfive_opcode_args(op)
    ccall((:libfive_opcode_args, libfive), Cint, (Cint,), op)
end

struct libfive_vars
    vars::Ptr{Ptr{Cvoid}}
    values::Ptr{Cfloat}
    size::UInt32
end

function libfive_vars_delete(j)
    ccall((:libfive_vars_delete, libfive), Cvoid, (Ptr{libfive_vars},), j)
end

mutable struct libfive_tree_ end

const libfive_tree = Ptr{libfive_tree_}

mutable struct libfive_evaluator_ end

const libfive_evaluator = Ptr{libfive_evaluator_}

# no prototype is found for this function at libfive.h:197:14, please use with caution
function libfive_tree_x()
    ccall((:libfive_tree_x, libfive), libfive_tree, ())
end

# no prototype is found for this function at libfive.h:202:14, please use with caution
function libfive_tree_y()
    ccall((:libfive_tree_y, libfive), libfive_tree, ())
end

# no prototype is found for this function at libfive.h:207:14, please use with caution
function libfive_tree_z()
    ccall((:libfive_tree_z, libfive), libfive_tree, ())
end

# no prototype is found for this function at libfive.h:212:14, please use with caution
function libfive_tree_var()
    ccall((:libfive_tree_var, libfive), libfive_tree, ())
end

function libfive_tree_is_var(t)
    ccall((:libfive_tree_is_var, libfive), Bool, (libfive_tree,), t)
end

function libfive_tree_const(f)
    ccall((:libfive_tree_const, libfive), libfive_tree, (Cfloat,), f)
end

function libfive_tree_get_const(t, success)
    ccall((:libfive_tree_get_const, libfive), Cfloat, (libfive_tree, Ptr{Bool}), t, success)
end

function libfive_tree_nullary(op)
    ccall((:libfive_tree_nullary, libfive), libfive_tree, (Cint,), op)
end

function libfive_tree_unary(op, a)
    ccall((:libfive_tree_unary, libfive), libfive_tree, (Cint, libfive_tree), op, a)
end

function libfive_tree_binary(op, a, b)
    ccall((:libfive_tree_binary, libfive), libfive_tree, (Cint, libfive_tree, libfive_tree), op, a, b)
end

function libfive_tree_id(t)
    ccall((:libfive_tree_id, libfive), Ptr{Cvoid}, (libfive_tree,), t)
end

function libfive_tree_eval_f(t, p)
    ccall((:libfive_tree_eval_f, libfive), Cfloat, (libfive_tree, libfive_vec3), t, p)
end

function libfive_tree_eval_r(t, r)
    ccall((:libfive_tree_eval_r, libfive), libfive_interval, (libfive_tree, libfive_region3), t, r)
end

function libfive_tree_eval_d(t, p)
    ccall((:libfive_tree_eval_d, libfive), libfive_vec3, (libfive_tree, libfive_vec3), t, p)
end

function libfive_tree_delete(ptr)
    ccall((:libfive_tree_delete, libfive), Cvoid, (libfive_tree,), ptr)
end

function libfive_tree_save(ptr, filename)
    ccall((:libfive_tree_save, libfive), Bool, (libfive_tree, Ptr{Cchar}), ptr, filename)
end

function libfive_tree_load(filename)
    ccall((:libfive_tree_load, libfive), libfive_tree, (Ptr{Cchar},), filename)
end

function libfive_tree_remap(p, x, y, z)
    ccall((:libfive_tree_remap, libfive), libfive_tree, (libfive_tree, libfive_tree, libfive_tree, libfive_tree), p, x, y, z)
end

function libfive_tree_optimized(t)
    ccall((:libfive_tree_optimized, libfive), libfive_tree, (libfive_tree,), t)
end

function libfive_tree_print(t)
    ccall((:libfive_tree_print, libfive), Ptr{Cchar}, (libfive_tree,), t)
end

function libfive_free_str(ptr)
    ccall((:libfive_free_str, libfive), Cvoid, (Ptr{Cchar},), ptr)
end

function libfive_tree_render_slice(tree, R, z, res)
    ccall((:libfive_tree_render_slice, libfive), Ptr{libfive_contours}, (libfive_tree, libfive_region2, Cfloat, Cfloat), tree, R, z, res)
end

function libfive_tree_render_slice3(tree, R, z, res)
    ccall((:libfive_tree_render_slice3, libfive), Ptr{libfive_contours3}, (libfive_tree, libfive_region2, Cfloat, Cfloat), tree, R, z, res)
end

function libfive_tree_save_slice(tree, R, z, res, f)
    ccall((:libfive_tree_save_slice, libfive), Cvoid, (libfive_tree, libfive_region2, Cfloat, Cfloat, Ptr{Cchar}), tree, R, z, res, f)
end

function libfive_tree_render_mesh(tree, R, res)
    ccall((:libfive_tree_render_mesh, libfive), Ptr{libfive_mesh}, (libfive_tree, libfive_region3, Cfloat), tree, R, res)
end

function libfive_tree_render_mesh_coords(tree, R, res)
    ccall((:libfive_tree_render_mesh_coords, libfive), Ptr{libfive_mesh_coords}, (libfive_tree, libfive_region3, Cfloat), tree, R, res)
end

function libfive_tree_save_mesh(tree, R, res, f)
    ccall((:libfive_tree_save_mesh, libfive), Bool, (libfive_tree, libfive_region3, Cfloat, Ptr{Cchar}), tree, R, res, f)
end

function libfive_evaluator_save_mesh(evaluator, R, f)
    ccall((:libfive_evaluator_save_mesh, libfive), Bool, (libfive_evaluator, libfive_region3, Ptr{Cchar}), evaluator, R, f)
end

function libfive_tree_save_meshes(trees, R, res, quality, f)
    ccall((:libfive_tree_save_meshes, libfive), Bool, (Ptr{libfive_tree}, libfive_region3, Cfloat, Cfloat, Ptr{Cchar}), trees, R, res, quality, f)
end

function libfive_tree_render_pixels(tree, R, z, res)
    ccall((:libfive_tree_render_pixels, libfive), Ptr{libfive_pixels}, (libfive_tree, libfive_region2, Cfloat, Cfloat), tree, R, z, res)
end

function libfive_tree_evaluator(tree, vars)
    ccall((:libfive_tree_evaluator, libfive), libfive_evaluator, (libfive_tree, libfive_vars), tree, vars)
end

function libfive_evaluator_update_vars(eval_tree, vars)
    ccall((:libfive_evaluator_update_vars, libfive), Bool, (libfive_evaluator, libfive_vars), eval_tree, vars)
end

function libfive_evaluator_delete(ptr)
    ccall((:libfive_evaluator_delete, libfive), Cvoid, (libfive_evaluator,), ptr)
end

function libfive_git_version()
    ccall((:libfive_git_version, libfive), Ptr{Cchar}, ())
end

function libfive_git_revision()
    ccall((:libfive_git_revision, libfive), Ptr{Cchar}, ())
end

function libfive_git_branch()
    ccall((:libfive_git_branch, libfive), Ptr{Cchar}, ())
end


text:
# ---
# title: 52. N-Queens II
# id: problem52
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Backtracking
# link: <https://leetcode.com/problems/n-queens-ii/description/>
# hidden: true
# ---
#
# The **n-queens** puzzle is the problem of placing `n` queens on an `n x n`
# chessboard such that no two queens attack each other.
#
# Given an integer `n`, return _the number of distinct solutions to the
# **n-queens puzzle**_.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/13/queens.jpg)
#
#
#
#     Input: n = 4
#     Output: 2
#     Explanation: There are two distinct solutions to the 4-queens puzzle as shown.
#
#
# **Example 2:**
#
#
#
#     Input: n = 1
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 9`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text: Dict{String,Any}("nest" => Dict{String,Any}("value" => Any[Dict{String,Any}("value" => Any[Dict{String,Any}("value" => "a","type" => "string")],"type" => "array"), Dict{String,Any}("value" => Any[Dict{String,Any}("value" => "b","type" => "string")],"type" => "array")],"type" => "array"))

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
"""
    BM_breadth!(S::BnBSolver,B::BnBModel,tL::Float64,tU::Float64,
                X1::Vector{Interval{Float64}},X2::Vector{Interval{Float64}},
                pos::Int64)

Stores two interval boxes `X1,X2` to the bottom of the stack along with their
respective lower, `tL` and upper bounds, `tU` and their position number in the
BnB tree. Also, assigns node numbers.
"""
function BM_breadth!(S::BnBSolver,B::BnBModel{T},tL::Float64,tU::Float64,
                             X1::Vector{T},
                             X2::Vector{T},pos::Int64) where {T}
  unshift!(B.box,X1,X2)
  unshift!(B.LBD,tL,tL)
  unshift!(B.UBD,tU,tU)
  unshift!(B.id,B.max_id+1,B.max_id+2)
  unshift!(B.pos,pos+1,pos+1)
  B.max_id += 2
end

"""
    BM_depth_best!(S::BnBSolver,B::BnBModel,tL::Float64,tU::Float64,
                X1::Vector{Interval{Float64}},X2::Vector{Interval{Float64}},
                pos::Int64)

Stores two interval boxes `X1,X2` to the top of the stack along with their
respective lower, `tL` and upper bounds, `tU` and their position number in the
BnB tree. Also, assigns node numbers.
"""
function BM_depth_best!(S::BnBSolver,B::BnBModel{T},tL::Float64,tU::Float64,
                             X1::Vector{T},
                             X2::Vector{T},pos::Int64) where {T}
  push!(B.box,X1,X2)
  push!(B.LBD,tL,tL)
  push!(B.UBD,tU,tU)
  push!(B.id,B.max_id+1,B.max_id+2)
  push!(B.pos,pos+1,pos+1)
  B.max_id += 2
end

"""
    BM_Single!(S::BnBSolver,B::BnBModel,tL::Float64,tU::Float64,
                             X::Vector{Interval{Float64}},pos::Int64)

Stores interval box `X` to the top of the stack along with their
respective lower, `tL` and upper bounds, `tU` and their position number in the
BnB tree.
"""
function BM_Single!(S::BnBSolver,B::BnBModel{T},tL::Float64,tU::Float64,
                             X::Vector{T},pos::Int64) where {T}
  push!(B.box,X)
  push!(B.LBD,tL)
  push!(B.UBD,tU)
  push!(B.id,B.max_id)
  push!(B.pos,pos)
end


text:
# The Shortest Path Faster Algorithm for single-source shortest path

###################################################################
#
#The type that capsulates the state of Shortest Path Faster Algorithm
#
###################################################################

using LightGraphs: nv, weights, outneighbors

"""
    struct SPFA <: ShortestPathAlgorithm

The structure used to configure and specify that [`shortest_paths`](@ref)
should use the [Shortest Path Faster Algorithm](https://en.wikipedia.org/wiki/Shortest_Path_Faster_Algorithm).
No additional configuration parameters are specified or required.

### Implementation Notes
`SPFA` supports the following shortest-path functionality:
- non-negative distance matrices / weights
- all destinations
"""
struct SPFA <: ShortestPathAlgorithm end
struct SPFAResult{T, U<:Integer} <: ShortestPathResult
    parents::Vector{U}
    dists::Vector{T}
end

function shortest_paths(g::AbstractGraph{U}, source::Integer, distmx::AbstractMatrix{T}, alg::SPFA) where {T, U<:Integer}

    nvg = nv(g)

    (source in 1:nvg) || throw(DomainError(source, "source should be in between 1 and $nvg"))
    dists = fill(typemax(T), nvg)
    parents = zeros(U, nvg)
    dists[source] = 0

    count = zeros(U, nvg)           # Vector to store the count of number of times a vertex goes in the queue.

    queue = Vector{U}()             # Vector used to implement queue
    inqueue = falses(nvg,1)         # BitArray to mark which vertices are in queue
    push!(queue, source)
    inqueue[source] = true

    @inbounds while !isempty(queue)
        v = popfirst!(queue)
        inqueue[v] = false

        @inbounds for v_neighbor in outneighbors(g, v)
            d = distmx[v,v_neighbor]
            if dists[v] + d < dists[v_neighbor]      # Relaxing edges
                dists[v_neighbor] = dists[v] + d
                parents[v_neighbor] = v

                if !inqueue[v_neighbor]
                    push!(queue,v_neighbor)
                    inqueue[v_neighbor] = true                   # Mark the vertex as inside queue.
                    count[v_neighbor] = count[v_neighbor]+1      # Increment the number of times a vertex enters a queue.

                    if count[v_neighbor] > nvg          # This step is just to check negative edge cycle. If this step is not there,
                        throw(NegativeCycleError())     # the algorithm will run infinitely in case of a negative weight cycle.
                                                        # If count[i]>nvg for any i belonging to [1,nvg], it means a negative edge
                                                        # cycle is present.
                    end
                end
            end
        end
    end

    return SPFAResult(parents, dists)
end

shortest_paths(g::AbstractGraph, s::Integer, alg::SPFA) = shortest_paths(g, s, weights(g), alg)

"""
Function which returns true if there is any negative weight cycle in the graph.
# Examples

```jldoctest
julia> g = complete_graph(3);

julia> d = [1 -3 1; -3 1 1; 1 1 1];

julia> has_negative_weight_cycle(g, d, SPFA())
true

julia> g = complete_graph(4);

julia> d = [1 1 -1 1; 1 1 -1 1; 1 1 1 1; 1 1 1 1];

julia> has_negative_weight_cycle(g, d, SPFA());
false
```
"""
function has_negative_weight_cycle(g::AbstractGraph, distmx::AbstractMatrix, alg::SPFA)
    try
        shortest_paths(g, 1, distmx, alg)
    catch e
        isa(e, ShortestPaths.NegativeCycleError) && return true
    end

    return false
end
has_negative_weight_cycle(g::AbstractGraph, ::SPFA) = false


text:
A = rand(3, 3)
B = ones(4, 4)
C = randn(5, 6)
K = A ⊗ B
X = collect(K)
v = rand(12)

K3 = A ⊗ B ⊗ C

@testset "multiply" begin
    @testset "Compatibility with other matrix ops" begin
        K = rand(3, 3) ⊗ randn(4, 4)
        K3 = rand(3, 3) ⊗ randn(4, 4) ⊗ randn(5, 6)
        KPow = kronecker(randn(3, 3), 3)

        for op in [adjoint, transpose, inv, pinv, conj]
            @test collect(op(K)) ≈ op(collect(K))
            if op != inv
                @test collect(op(K3)) ≈ op(collect(K3))
            else
                @test_throws DimensionMismatch op(K3)
            end
        end
    end

    @testset "Compatibility with Matrix Product" begin
        A_f = [randn(3, 4), randn(5, 6)]
        B_f = [randn(4, 2), randn(6, 7)]

        A = kronecker(A_f...)
        B = kronecker(B_f...)

        @test collect(A * B) ≈ collect(kronecker(map(*, A_f, B_f)...))

        C = randn(3, 3)
        D = randn(3, 2)
        C_pow = kronecker(C, 3)
        D_pow = kronecker(D, 3)
        @test collect(C_pow * D_pow) ≈ collect(kronecker(C*D, 3))

        @test_throws DimensionMismatch (D_pow * C_pow)
    end

    @testset "Vec trick" begin

        @test K * v ≈ X * v

        V = rand(4, 3)
        @test K * vec(V) ≈ X * vec(V)
        @test_throws DimensionMismatch K * V
        @test_throws DimensionMismatch K * reshape(V, 2, 6)

        v3 = randn(size(K3, 2))
        @test K3 * v3 ≈ collect(K3) * v3
        u = similar(v)
        @test mul!(u, K, v) ≈ X * v
    end

    @testset "Reshaped vec trick" begin
        @test K * v ≈ X * v

        V = sprand(4, 3, 1.0)
        @test K * vec(V) ≈ X * vec(V)
        @test_throws DimensionMismatch K * V
        @test_throws DimensionMismatch K * reshape(V, 2, 6)
        K3 = A ⊗ B ⊗ C
        v3 = sprand(size(K3, 2),1.0)
        @test K3 * vec(v3) ≈ collect(K3) * v3
        u = similar(v)
        @test mul!(u, K, v) ≈ X * v

    end

    @testset "AbstractKroneckerProduct * AbstractMatrix" begin
        rng = MersenneTwister(123456)
        a, b, x = randn(rng, 30, 20), randn(rng, 40, 50), randn(rng, 1000, 1100)

        @test kron(a, b) * x ≈ (a ⊗ b) * x
        @test kron(a, Eye(50)) * x ≈ (a ⊗ Eye(50)) * x
        @test kron(Eye(20), b) * x ≈ (Eye(20) ⊗ b) * x
    end

    @testset "sum" begin
        @test sum(K) ≈ sum(X)
        @test sum(K3) ≈ sum(collect(K3))

        @test sum(K, dims=1) ≈ sum(X, dims=1)
        @test sum(K3, dims=2) ≈ sum(collect(K3), dims=2)
        @test sum(K3, dims=2) isa AbstractKroneckerProduct

        @test sum(kronecker(A, 3)) ≈ sum(kron(A, A, A))
    end

    @inline function compare_against_kron(K, k, x, X, y, Y)
        @test collect(K) ≈ k

        res1 = (K * x)
        res2 = (K * X)
        @test res1 ≈ (k * x)
        @test res2 ≈ (k * X)

        if K isa AbstractKroneckerProduct
            res3 = (K \ y)
            res4 = (K \ Y)
            @test res3 ≈ (k \ y)
            @test res4 ≈ (k \ Y)
        end

        if size(K, 1) == size(x, 1)
            res1a = (x' * K)
            res1t = (transpose(x) * K)

            @test res1a ≈ (x' * k)
            @test res1t ≈ (transpose(x) * k)
            @test (K' * x) ≈ (k' * x) ≈ res1a'
            @test (transpose(K) * x) ≈ (transpose(k) * x) ≈ transpose(res1t)
        else
            @test_throws DimensionMismatch (x' * K)
            @test_throws DimensionMismatch (transpose(x) * K)
        end

        if size(K, 1) == size(X, 1)
            res2a = (X' * K)
            res2t = (transpose(X) * K)

            @test res2a ≈ (X' * k)
            @test res2t ≈ (transpose(X) * k)
            @test (K' * X) ≈ (k' * X) ≈ res2a'
            @test (transpose(K) * X) ≈ (transpose(k) * X) ≈ transpose(res2t)
        else
            @test_throws DimensionMismatch (X' * K)
            @test_throws DimensionMismatch (transpose(X) * K)
        end

    end

    @testset "2-factor square KroneckerProduct" begin
        A = randn(10, 10)
        B = randn(7, 7)
        x = randn(size(A, 2) * size(B, 2))
        X = randn(size(A, 2) * size(B, 2), 4)
        y = randn(size(A, 1) * size(B, 1))
        Y = randn(size(A, 1) * size(B, 1), 4)

        K1 = kronecker(A, B)
        K2 = kronecker(B, A)
        KS = kroneckersum(A, B)

        k1 = kron(A, B)
        k2 = kron(B, A)
        ks = kron(A, Diagonal(I, size(B, 1))) + kron(Diagonal(I, size(A, 1)), B)

        for (K, k) in [(K1, k1), (K2, k2), (KS, ks)]
            compare_against_kron(K, k, x, X, y, Y)
        end
    end

    @testset "2-factor rectangular KroneckerProduct" begin
        A = randn(10, 8)
        B = randn(7, 9)
        x = randn(size(A, 2) * size(B, 2))
        X = randn(size(A, 2) * size(B, 2), 4)
        y = randn(size(A, 1) * size(B, 1))
        Y = randn(size(A, 1) * size(B, 1), 4)

        K1 = kronecker(A, B)
        K2 = kronecker(B, A)
        k1 = kron(A, B)
        k2 = kron(B, A)

        for (K, k) in [(K1, k1), (K2, k2)]
            compare_against_kron(K, k, x, X, y, Y)
        end
    end

    @testset "3-factor rectangular KroneckerProduct" begin
        A = randn(10, 8)
        B = randn(7, 9)
        C = randn(6, 4)
        x = randn(size(A, 2) * size(B, 2) * size(C, 2))
        X = randn(size(A, 2) * size(B, 2) * size(C, 2), 4)
        y = randn(size(A, 1) * size(B, 1) * size(C, 1))
        Y = randn(size(A, 1) * size(B, 1) * size(C, 1), 4)

        K1 = kronecker(A, B, C)
        K2 = kronecker(B, A, C)
        k1 = kron(A, B, C)
        k2 = kron(B, A, C)

        for (K, k) in [(K1, k1), (K2, k2)]
            compare_against_kron(K, k, x, X, y, Y)
        end
    end

    @testset "10-factor square KroneckerProduct" begin
        matrices = [randn(2,2) for i in 1:10]
        v = randn(2^10)
        V = randn(2^10, 5)

        KP = kronecker(matrices...)
        KS = kroneckersum(matrices...)
        kp = kron(matrices...)
        ks = collect(KS)

        for (K, k) in [(KP, kp), (KS, ks)]
            compare_against_kron(K, k, v, V, v, V)
        end
    end

    @testset "10-factor square KroneckerPower" begin
        # reduce condition number of the matrix to avoid accidentally triggering a
        #  test failure
        A = randn(2,2) + 20I
        A /= opnorm(A)

        v = randn(2^10)
        V = randn(2^10, 5)

        K = kronecker(ntuple(_ -> A, 10)...)
        k = kron(ntuple(_ -> A, 10)...)

        compare_against_kron(K, k, v, V, v, V)
    end

    @testset "10-factor rectangular KroneckerProduct" begin
        matrices = [randn(3,2) for i in 1:10]
        x = randn(2^10)
        X = randn(2^10, 5)
        y = randn(3^10)
        Y = randn(3^10, 5)

        K = kronecker(matrices...)
        k = kron(matrices...)

        compare_against_kron(K, k, x, X, y, Y)
    end

    @testset "10-factor mixed KroneckerProduct" begin
        matrices1 = [randn(2,2) for i in 1:5]
        matrices2 = [randn(3,2) for i in 1:5]
        x = randn(2^10)
        X = randn(2^10, 5)
        y = randn(2^5 * 3^5)
        Y = randn(2^5 * 3^5, 5)

        K = kronecker(matrices1..., matrices2...)
        k = kron(matrices1..., matrices2...)

        @test collect(K) ≈ k

        compare_against_kron(K, k, x, X, y, Y)
    end

end


text:
"""
Ordinary differential equation in radial direction for the
hyperboloidally compactified Teukolsky equation.
"""
module RadialODE

export eig_vals_vecs_c

include("CustomTypes.jl")
include("Chebyshev.jl")
using .CustomTypes
import .Chebyshev as CH

#import IterativeSolvers: invpowm
using LinearAlgebra
using GenericSchur
using SparseArrays
using Printf: @printf

using ApproxFun

"""
    radial_discretized_eqn(
      nr::myI,
      s::myI,
      m::myI,
      a::myF,
      bhm::myF,
      om::myC,
      rmin::myF,
      rmax::myF
    )
"""
function radial_discretized_eqn(
      nr::myI,
      s::myI,
      m::myI,
      a::myF,
      bhm::myF,
      om::myC,
      rmin::myF,
      rmax::myF
)
   ## Using position space method
   #=
   D1 = CH.mat_D1(rmin,rmax,nr)
   D2 = D1*D1

   Id = sparse(I,nr,nr)

   X1 = CH.mat_X(rmin,rmax,nr)
   X2 = X1*X1
   X3 = X1*X2
   X4 = X1*X3
   =#

   ## Using Ultraspherical polynomial method
   X1 = ApproxFun.Fun(rmin..rmax);
   X2 = X1^2
   X3 = X1^3
   X4 = X1^4
   D1 = ApproxFun.Derivative()
   D2 = D1^2

   A = (
        tomyC(2*im)*om*I
        -
        tomyF(2)*(tomyF(1) + s)*X1
        +
        tomyF(2)*(
             im*om*((a^2) - tomyF(8)*(bhm^2))
             +
             im*m*a
             +
             (s + tomyF(3))*bhm
         )*X2
        +
        tomyF(4)*(tomyF(2)*im*om*bhm - tomyF(1))*(a^2)*X3
   )
   B = (
     (
      ((a^2) - tomyF(16)*(bhm^2))*(om^2)
      +
      tomyF(2)*(m*a + tomyF(2)*im*s*bhm)*om
     )*I
     +
     tomyF(2)*(
      tomyF(4)*((a^2) - tomyF(4)*(bhm^2))*bhm*(om^2)
      +
      (tomyF(4)*m*a*bhm - tomyF(4)*im*(s + tomyF(2))*(bhm^2) + im*(a^2))*om
      +
      im*m*a
      +
      (s + tomyF(1))*bhm
     )*X1
     +
     tomyF(2)*(
          tomyF(8)*(bhm^2)*(om^2)
          +
          tomyF(6)*im*bhm*om
          -
          tomyF(1)
     )*(a^2)*X2
   )

   return (
      -
      (X2 - tomyF(2)*bhm*X3 + (a^2)*X4)*D2
      +
      A*D1
      +
      B
   ) : ApproxFun.space(X1)
end

"""
    radial_discretized_eqn_p(
      np::myI,
      s::myI,
      m::myI,
      a::myF,
      bhm::myF,
      om::myC
    )
"""
function radial_discretized_eqn_p(
      np::myI,
      s::myI,
      m::myI,
      a::myF,
      bhm::myF,
      om::myC
)
   X = Vector{myI}()
   Y = Vector{myI}()
   V = Vector{myC}()

   A0 = tomyF(2)*im*om
   A1 = -tomyF(2)*(tomyF(1)+s)
   A2 = tomyF(2)*(im*om*((a^2) - tomyF(8)*(bhm^2)) + im*m*a + (s+tomyF(3))*bhm)
   A3 = tomyF(4)*(tomyF(2)*im*om*bhm - tomyF(1))*(a^2)

   B0 = ((a^2) - tomyF(16)*(bhm^2))*(om^2) + tomyF(2)*(m*a + tomyF(2)*im*s*bhm)*om
   B1 = tomyF(2)*(
           tomyF(4)*((a^2) - tomyF(4)*(bhm^2))*bhm*(om^2)
           +
           (tomyF(4)*m*a*bhm - tomyF(4)*im*(s+2)*(bhm^2) + im*(a^2))*om
           +
           im*m*a
           +
           (s + tomyF(1))*bhm
          )
   B2 = tomyF(2)*(tomyF(8)*((bhm*om)^2) + tomyF(6)*im*om*bhm - tomyF(1))*(a^2)

   push!(X,1)
   push!(Y,1)
   push!(V,B0)

   push!(X,1)
   push!(Y,2)
   push!(V,A0)

   push!(X,2)
   push!(Y,1)
   push!(V,B1)

   push!(X,2)
   push!(Y,2)
   push!(V,A1+B0)

   push!(X,2)
   push!(Y,3)
   push!(V,2*A0)

   for i=3:(np-1)
      n = i-1
      push!(X,i)
      push!(Y,i-2)
      push!(V,-(n-2)*(n-3)*(a^2) + (n-2)*A3 + B2)

      push!(X,i)
      push!(Y,i-1)
      push!(V,2*(n-1)*(n-2)*bhm + (n-1)*A2 + B1)

      push!(X,i)
      push!(Y,i)
      push!(V,-n*(n-1) + n*A1 + B0)

      push!(X,i)
      push!(Y,i+1)
      push!(V,(n+1)*A0)
   end

   n = np-1
   push!(X,np)
   push!(Y,np-2)
   push!(V,-(n-2)*(n-3)*(a^2) + (n-2)*A3 + B2)

   push!(X,np)
   push!(Y,np-1);
   push!(V,2*(n-1)*(n-2)*bhm + (n-1)*A2 + B1)

   push!(X,np)
   push!(Y,np)
   push!(V,-n*(n-1) + n*A1 + B0)

   mat = dropzeros(sparse(X,Y,V))
   for i=1:np
      mat[:,i] .*= 2^(i-1)
   end

   mat = dropzeros(sparse(X,Y,V))

   #for i=1:np
   #   mat[:,i] .*= 2^(i-1)
   #end

   return mat
end

"""
    eig_vals_vecs_c(
      nr::myI,
      s::myI,
      m::myI,
      a::myF,
      om::myC
    )

Compute eigenvectors and eigenvalues for the radial equation
using a pseudospectral Chebyshev polynomial method.
The black hole mass is always one.
"""
function eig_vals_vecs_c(
      nr::myI,
      s::myI,
      m::myI,
      a::myF,
      om::myC
   )
   bhm  = tomyF(1) ## always have unit black hole mass
   rmin = tomyF(0)
   rmax = tomyF(abs(a)>0 ? (bhm/(a^2))*(tomyF(1) - sqrt(tomyF(1)- ((a/bhm)^2))) : tomyF(0.5)/bhm)

   L = radial_discretized_eqn(nr,s,m,a,bhm,om,rmin,rmax)
   t = ApproxFun.eigs(L,nr,tolerance=1e-12)

   mini = argmin(abs.(t[1]))

   #nMat = radial_discretized_eqn_p(nr,s,m,a,bhm,om)
   #nt = eigen(Matrix(nMat), permute=true, scale=true, sortby=abs)
   #println(nMat)

   return -t[1][mini], t[2][mini], CH.cheb_pts(rmin,rmax,nr)
end

end # module


text:
"""
    curl(dx::AbstractVector, A::AbstractArray)

Calculate 2nd order cell-centered ∇×A where `A` is a 4D array of size (3, nx, ny, nz) and
`dx` is a vector of grid intervals in each dimension.
"""
function curl(dx::AbstractVector, A::AbstractArray{T,N}) where {T,N}
   @assert N == 4 && length(dx) == 3 "Input vector shall be indexed in 3D!"

   @views Ax, Ay, Az = A[1,:,:,:], A[2,:,:,:], A[3,:,:,:]

   B = zeros(T, size(A))
   @views Bx, By, Bz = B[1,:,:,:], B[2,:,:,:], B[3,:,:,:]

   invdx = inv.(dx)

   if any(==(1), size(A)) # 2D
      if size(A,4) == 1
         @inbounds for j in 2:size(A,3)-1, i in 2:size(A,2)-1
            ∂Ay∂x = (Ay[i-1,j,1] - 2Ay[i,j,1] + Ay[i+1,j,1]) * invdx[1]
            ∂Az∂x = (Az[i-1,j,1] - 2Az[i,j,1] + Az[i+1,j,1]) * invdx[1]
            ∂Ax∂y = (Ax[i-1,j,1] - 2Ax[i,j,1] + Ax[i+1,j,1]) * invdx[2]
            ∂Az∂y = (Az[i-1,j,1] - 2Az[i,j,1] + Az[i+1,j,1]) * invdx[2]
            ∂Ax∂z = (Ax[i-1,j,1] - 2Ax[i,j,1] + Ax[i+1,j,1]) * invdx[3]
            ∂Ay∂z = (Ay[i-1,j,1] - 2Ay[i,j,1] + Ay[i+1,j,1]) * invdx[3]

            Bx[i,j,1] = ∂Az∂y - ∂Ay∂z
            By[i,j,1] = ∂Ax∂z - ∂Az∂x
            Bz[i,j,1] = ∂Ay∂x - ∂Ax∂y
         end
      elseif size(A,3) == 1
         @inbounds for k in 2:size(A,4)-1, i in 2:size(A,2)-1
            ∂Ay∂x = (Ay[i-1,1,k] - 2Ay[i,1,k] + Ay[i+1,1,k]) * invdx[1]
            ∂Az∂x = (Az[i-1,1,k] - 2Az[i,1,k] + Az[i+1,1,k]) * invdx[1]
            ∂Ax∂y = (Ax[i-1,1,k] - 2Ax[i,1,k] + Ax[i+1,1,k]) * invdx[2]
            ∂Az∂y = (Az[i-1,1,k] - 2Az[i,1,k] + Az[i+1,1,k]) * invdx[2]
            ∂Ax∂z = (Ax[i-1,1,k] - 2Ax[i,1,k] + Ax[i+1,1,k]) * invdx[3]
            ∂Ay∂z = (Ay[i-1,1,k] - 2Ay[i,1,k] + Ay[i+1,1,k]) * invdx[3]

            Bx[i,1,k] = ∂Az∂y - ∂Ay∂z
            By[i,1,k] = ∂Ax∂z - ∂Az∂x
            Bz[i,1,k] = ∂Ay∂x - ∂Ax∂y
         end
      end
   else # 3D
      @inbounds for k in 2:size(A,4)-1, j in 2:size(A,3)-1, i in 2:size(A,2)-1
         ∂Ay∂x = (Ay[i-1,j,k] - 2Ay[i,j,k] + Ay[i+1,j,k]) * invdx[1]
         ∂Az∂x = (Az[i-1,j,k] - 2Az[i,j,k] + Az[i+1,j,k]) * invdx[1]
         ∂Ax∂y = (Ax[i-1,j,k] - 2Ax[i,j,k] + Ax[i+1,j,k]) * invdx[2]
         ∂Az∂y = (Az[i-1,j,k] - 2Az[i,j,k] + Az[i+1,j,k]) * invdx[2]
         ∂Ax∂z = (Ax[i-1,j,k] - 2Ax[i,j,k] + Ax[i+1,j,k]) * invdx[3]
         ∂Ay∂z = (Ay[i-1,j,k] - 2Ay[i,j,k] + Ay[i+1,j,k]) * invdx[3]

         Bx[i,j,k] = ∂Az∂y - ∂Ay∂z
         By[i,j,k] = ∂Ax∂z - ∂Az∂x
         Bz[i,j,k] = ∂Ay∂x - ∂Ax∂y
      end
   end
   B
end


text:
export CalculateComparer, fn, reset!

type CalculateComparer{C<:AbstractComparer} <: AbstractComparer
    factory::Function
    c::C
    _results::Dict{Any, Nullable{Any}}
end

CalculateComparer(factory::Function) = CalculateComparer{Comparer}(factory, Comparer(), Dict{Any, Nullable{Any}}())

CalculateComparer(factory::Function, fn::Function) = CalculateComparer{Comparer}(factory, Comparer(fn), Dict{Any, Nullable{Any}}())

CalculateComparer{C}(factory::Function, c::C) = CalculateComparer{C}(factory, c, Dict{Any, Nullable{Any}}())

fn(cc::CalculateComparer) = (a, b) ->
begin
    baseFn = fn(cc.c)
    aResult = create!(cc, a)
    bResult = create!(cc, b)
    baseFn(aResult, bResult)
end

function create!(cc::CalculateComparer, source)
    result = get(cc._results, source, Nullable{Any}())
    if !isnull(result)
        return get(result)
    end
    cc._results[source] = cc.factory(source)
end

function reset!(cc::CalculateComparer)
    reset!(cc.c)
    empty!(cc._results)
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
int_rules_5_1_3 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 5.1.3*(d+e*x^2)^p*(a+b*arcsin(c*x))^n =#
    #= Antiderivative(((~a')+(~b')*asin((~c')*(~x)))^(~n')/sqrt((~d)+(~e')*(~x)^2),~x) :=  1/c*Simp(sqrt(1-c^2*x^2)/sqrt(d+e*x^2))*Subst(Antiderivative((a+b*x)^n,x),x, asin(c*x)) <-- FreeQ([a,b,c,d,e,n],x) && EqQ(c^2*d+e,0) =#
    #= Antiderivative(((~a')+(~b')*acos((~c')*(~x)))^(~n')/sqrt((~d)+(~e')*(~x)^2),~x) :=  -1/c*Simp(sqrt(1-c^2*x^2)/sqrt(d+e*x^2))*Subst(Antiderivative((a+b*x)^n,x),x, acos(c*x)) <-- FreeQ([a,b,c,d,e,n],x) && EqQ(c^2*d+e,0) =#
    @apply_utils Antiderivative(1 / (sqrt(~d + ~(e') * (~x) ^ 2) * (~(a') + ~(b') * asin(~(c') * ~x))), ~x) => (1 / (~b * ~c)) * Simp(sqrt(1 - (~c) ^ 2 * (~x) ^ 2) / sqrt(~d + ~e * (~x) ^ 2)) * log(~a + ~b * asin(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && EqQ((~c) ^ 2 * ~d + ~e, 0)
    @apply_utils Antiderivative(1 / (sqrt(~d + ~(e') * (~x) ^ 2) * (~(a') + ~(b') * acos(~(c') * ~x))), ~x) => ((-1 / (~b * ~c)) * Simp(sqrt(1 - (~c) ^ 2 * (~x) ^ 2) / sqrt(~d + ~e * (~x) ^ 2)) * log(~a + ~b * acos(~c * ~x))) / (~b * ~c * sqrt(~d)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && EqQ((~c) ^ 2 * ~d + ~e, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n') / sqrt(~d + ~(e') * (~x) ^ 2), ~x) => (1 / (~b * ~c * (~n + 1))) * Simp(sqrt(1 - (~c) ^ 2 * (~x) ^ 2) / sqrt(~d + ~e * (~x) ^ 2)) * (~a + ~b * asin(~c * ~x)) ^ (~n + 1) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n') / sqrt(~d + ~(e') * (~x) ^ 2), ~x) => (-1 / (~b * ~c * (~n + 1))) * Simp(sqrt(1 - (~c) ^ 2 * (~x) ^ 2) / sqrt(~d + ~e * (~x) ^ 2)) * (~a + ~b * acos(~c * ~x)) ^ (~n + 1) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asin(~(c') * ~x)), ~x) => With([u = IntHide((~d + ~e * (~x) ^ 2) ^ ~p, ~x)], Dist(~a + ~b * asin(~c * ~x), u, ~x) - ~b * ~c * Antiderivative(SimplifyIntegrand(u / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && IGtQ(~p, 0))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acos(~(c') * ~x)), ~x) => With([u = IntHide((~d + ~e * (~x) ^ 2) ^ ~p, ~x)], Dist(~a + ~b * acos(~c * ~x), u, ~x) + ~b * ~c * Antiderivative(SimplifyIntegrand(u / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && IGtQ(~p, 0))
    @apply_utils Antiderivative(sqrt(~d + ~(e') * (~x) ^ 2) * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => ((~x * sqrt(~d + ~e * (~x) ^ 2) * (~a + ~b * asin(~c * ~x)) ^ ~n) / 2 - ((~b * ~c * ~n) / 2) * Simp(sqrt(~d + ~e * (~x) ^ 2) / sqrt(1 - (~c) ^ 2 * (~x) ^ 2)) * Antiderivative(~x * (~a + ~b * asin(~c * ~x)) ^ (~n - 1), ~x)) + (1 / 2) * Simp(sqrt(~d + ~e * (~x) ^ 2) / sqrt(1 - (~c) ^ 2 * (~x) ^ 2)) * Antiderivative((~a + ~b * asin(~c * ~x)) ^ ~n / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && GtQ(~n, 0))
    @apply_utils Antiderivative(sqrt(~d + ~(e') * (~x) ^ 2) * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => (~x * sqrt(~d + ~e * (~x) ^ 2) * (~a + ~b * acos(~c * ~x)) ^ ~n) / 2 + ((~b * ~c * ~n) / 2) * Simp(sqrt(~d + ~e * (~x) ^ 2) / sqrt(1 - (~c) ^ 2 * (~x) ^ 2)) * Antiderivative(~x * (~a + ~b * acos(~c * ~x)) ^ (~n - 1), ~x) + (1 / 2) * Simp(sqrt(~d + ~e * (~x) ^ 2) / sqrt(1 - (~c) ^ 2 * (~x) ^ 2)) * Antiderivative((~a + ~b * acos(~c * ~x)) ^ ~n / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && GtQ(~n, 0))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => ((~x * (~d + ~e * (~x) ^ 2) ^ ~p * (~a + ~b * asin(~c * ~x)) ^ ~n) / (2 * ~p + 1) + ((2 * ~d * ~p) / (2 * ~p + 1)) * Antiderivative((~d + ~e * (~x) ^ 2) ^ (~p - 1) * (~a + ~b * asin(~c * ~x)) ^ ~n, ~x)) - ((~b * ~c * ~n) / (2 * ~p + 1)) * Simp((~d + ~e * (~x) ^ 2) ^ ~p / (1 - (~c) ^ 2 * (~x) ^ 2) ^ ~p) * Antiderivative(~x * (1 - (~c) ^ 2 * (~x) ^ 2) ^ (~p - 1 / 2) * (~a + ~b * asin(~c * ~x)) ^ (~n - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (GtQ(~n, 0) && GtQ(~p, 0)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => (~x * (~d + ~e * (~x) ^ 2) ^ ~p * (~a + ~b * acos(~c * ~x)) ^ ~n) / (2 * ~p + 1) + ((2 * ~d * ~p) / (2 * ~p + 1)) * Antiderivative((~d + ~e * (~x) ^ 2) ^ (~p - 1) * (~a + ~b * acos(~c * ~x)) ^ ~n, ~x) + ((~b * ~c * ~n) / (2 * ~p + 1)) * Simp((~d + ~e * (~x) ^ 2) ^ ~p / (1 - (~c) ^ 2 * (~x) ^ 2) ^ ~p) * Antiderivative(~x * (1 - (~c) ^ 2 * (~x) ^ 2) ^ (~p - 1 / 2) * (~a + ~b * acos(~c * ~x)) ^ (~n - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (GtQ(~n, 0) && GtQ(~p, 0)))
    @apply_utils Antiderivative((~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n') / (~d + ~(e') * (~x) ^ 2) ^ (3 / 2), ~x) => (~x * (~a + ~b * asin(~c * ~x)) ^ ~n) / (~d * sqrt(~d + ~e * (~x) ^ 2)) - ((~b * ~c * ~n) / ~d) * Simp(sqrt(1 - (~c) ^ 2 * (~x) ^ 2) / sqrt(~d + ~e * (~x) ^ 2)) * Antiderivative((~x * (~a + ~b * asin(~c * ~x)) ^ (~n - 1)) / (1 - (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && GtQ(~n, 0))
    @apply_utils Antiderivative((~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n') / (~d + ~(e') * (~x) ^ 2) ^ (3 / 2), ~x) => (~x * (~a + ~b * acos(~c * ~x)) ^ ~n) / (~d * sqrt(~d + ~e * (~x) ^ 2)) + ((~b * ~c * ~n) / ~d) * Simp(sqrt(1 - (~c) ^ 2 * (~x) ^ 2) / sqrt(~d + ~e * (~x) ^ 2)) * Antiderivative((~x * (~a + ~b * acos(~c * ~x)) ^ (~n - 1)) / (1 - (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && GtQ(~n, 0))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => (-(~x) * (~d + ~e * (~x) ^ 2) ^ (~p + 1) * (~a + ~b * asin(~c * ~x)) ^ ~n) / (2 * ~d * (~p + 1)) + ((2 * ~p + 3) / (2 * ~d * (~p + 1))) * Antiderivative((~d + ~e * (~x) ^ 2) ^ (~p + 1) * (~a + ~b * asin(~c * ~x)) ^ ~n, ~x) + ((~b * ~c * ~n) / (2 * (~p + 1))) * Simp((~d + ~e * (~x) ^ 2) ^ ~p / (1 - (~c) ^ 2 * (~x) ^ 2) ^ ~p) * Antiderivative(~x * (1 - (~c) ^ 2 * (~x) ^ 2) ^ (~p + 1 / 2) * (~a + ~b * asin(~c * ~x)) ^ (~n - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (GtQ(~n, 0) && (LtQ(~p, -1) && NeQ(~p, -3 / 2))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => ((-(~x) * (~d + ~e * (~x) ^ 2) ^ (~p + 1) * (~a + ~b * acos(~c * ~x)) ^ ~n) / (2 * ~d * (~p + 1)) + ((2 * ~p + 3) / (2 * ~d * (~p + 1))) * Antiderivative((~d + ~e * (~x) ^ 2) ^ (~p + 1) * (~a + ~b * acos(~c * ~x)) ^ ~n, ~x)) - ((~b * ~c * ~n) / (2 * (~p + 1))) * Simp((~d + ~e * (~x) ^ 2) ^ ~p / (1 - (~c) ^ 2 * (~x) ^ 2) ^ ~p) * Antiderivative(~x * (1 - (~c) ^ 2 * (~x) ^ 2) ^ (~p + 1 / 2) * (~a + ~b * acos(~c * ~x)) ^ (~n - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (GtQ(~n, 0) && (LtQ(~p, -1) && NeQ(~p, -3 / 2))))
    @apply_utils Antiderivative((~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n') / (~d + ~(e') * (~x) ^ 2), ~x) => (1 / (~c * ~d)) * Subst(Antiderivative((~a + ~b * ~x) ^ ~n * sec(~x), ~x), ~x, asin(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && IGtQ(~n, 0))
    @apply_utils Antiderivative((~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n') / (~d + ~(e') * (~x) ^ 2), ~x) => (-1 / (~c * ~d)) * Subst(Antiderivative((~a + ~b * ~x) ^ ~n * csc(~x), ~x), ~x, acos(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && IGtQ(~n, 0))
    #= Antiderivative(((~d)+(~e')*(~x)^2)^(~p')*((~a')+(~b')*asin((~c')*(~x)))^(~n),~x) :=  d^p*(1-c^2*x^2)^(p+1/2)*(a+b*asin(c*x))^(n+1)/(b*c*(n+1)) + c*d^p*(2*p+1)/(b*(n+1))*Antiderivative(x*(1-c^2*x^2)^(p-1/2)*(a+b*asin(c*x))^ (n+1),x) <-- FreeQ([a,b,c,d,e,p],x) && EqQ(c^2*d+e,0) && LtQ(n,-1) &&  (IntegerQ(p) || GtQ(d,0)) =#
    #= Antiderivative(((~d)+(~e')*(~x)^2)^(~p')*((~a')+(~b')*acos((~c')*(~x)))^(~n),~x) :=  -d^p*(1-c^2*x^2)^(p+1/2)*(a+b*acos(c*x))^(n+1)/(b*c*(n+1)) - c*d^p*(2*p+1)/(b*(n+1))*Antiderivative(x*(1-c^2*x^2)^(p-1/2)*(a+b*acos(c*x))^ (n+1),x) <-- FreeQ([a,b,c,d,e,p],x) && EqQ(c^2*d+e,0) && LtQ(n,-1) &&  (IntegerQ(p) || GtQ(d,0)) =#
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~n, ~x) => (sqrt(1 - (~c) ^ 2 * (~x) ^ 2) * (~d + ~e * (~x) ^ 2) ^ ~p * (~a + ~b * asin(~c * ~x)) ^ (~n + 1)) / (~b * ~c * (~n + 1)) + ((~c * (2 * ~p + 1)) / (~b * (~n + 1))) * Simp((~d + ~e * (~x) ^ 2) ^ ~p / (1 - (~c) ^ 2 * (~x) ^ 2) ^ ~p) * Antiderivative(~x * (1 - (~c) ^ 2 * (~x) ^ 2) ^ (~p - 1 / 2) * (~a + ~b * asin(~c * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && LtQ(~n, -1))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~n, ~x) => (-(sqrt(1 - (~c) ^ 2 * (~x) ^ 2)) * (~d + ~e * (~x) ^ 2) ^ ~p * (~a + ~b * acos(~c * ~x)) ^ (~n + 1)) / (~b * ~c * (~n + 1)) - ((~c * (2 * ~p + 1)) / (~b * (~n + 1))) * Simp((~d + ~e * (~x) ^ 2) ^ ~p / (1 - (~c) ^ 2 * (~x) ^ 2) ^ ~p) * Antiderivative(~x * (1 - (~c) ^ 2 * (~x) ^ 2) ^ (~p - 1 / 2) * (~a + ~b * acos(~c * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && LtQ(~n, -1))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => (1 / (~b * ~c)) * Simp((~d + ~e * (~x) ^ 2) ^ ~p / (1 - (~c) ^ 2 * (~x) ^ 2) ^ ~p) * Subst(Antiderivative((~x) ^ ~n * cos(-(~a) / ~b + ~x / ~b) ^ (2 * ~p + 1), ~x), ~x, ~a + ~b * asin(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && IGtQ(2 * ~p, 0))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => (-1 / (~b * ~c)) * Simp((~d + ~e * (~x) ^ 2) ^ ~p / (1 - (~c) ^ 2 * (~x) ^ 2) ^ ~p) * Subst(Antiderivative((~x) ^ ~n * sin(-(~a) / ~b + ~x / ~b) ^ (2 * ~p + 1), ~x), ~x, ~a + ~b * acos(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && IGtQ(2 * ~p, 0))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asin(~(c') * ~x)), ~x) => With([u = IntHide((~d + ~e * (~x) ^ 2) ^ ~p, ~x)], Dist(~a + ~b * asin(~c * ~x), u, ~x) - ~b * ~c * Antiderivative(SimplifyIntegrand(u / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~c) ^ 2 * ~d + ~e, 0) && (IGtQ(~p, 0) || ILtQ(~p + 1 / 2, 0)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acos(~(c') * ~x)), ~x) => With([u = IntHide((~d + ~e * (~x) ^ 2) ^ ~p, ~x)], Dist(~a + ~b * acos(~c * ~x), u, ~x) + ~b * ~c * Antiderivative(SimplifyIntegrand(u / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~c) ^ 2 * ~d + ~e, 0) && (IGtQ(~p, 0) || ILtQ(~p + 1 / 2, 0)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~a + ~b * asin(~c * ~x)) ^ ~n, (~d + ~e * (~x) ^ 2) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (NeQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~p) && (GtQ(~p, 0) || IGtQ(~n, 0))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~a + ~b * acos(~c * ~x)) ^ ~n, (~d + ~e * (~x) ^ 2) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (NeQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~p) && (GtQ(~p, 0) || IGtQ(~n, 0))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => Unintegrable((~d + ~e * (~x) ^ 2) ^ ~p * (~a + ~b * asin(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p], ~x)
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => Unintegrable((~d + ~e * (~x) ^ 2) ^ ~p * (~a + ~b * acos(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p], ~x)
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~p * (~f + ~(g') * ~x) ^ ~q * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => ((-((~d) ^ 2) * ~g) / ~e) ^ ~q * Antiderivative((~d + ~e * ~x) ^ (~p - ~q) * (1 - (~c) ^ 2 * (~x) ^ 2) ^ ~q * (~a + ~b * asin(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && (EqQ(~e * ~f + ~d * ~g, 0) && (EqQ((~c) ^ 2 * (~d) ^ 2 - (~e) ^ 2, 0) && (HalfIntegerQ(~p, ~q) && (GeQ(~p - ~q, 0) && (GtQ(~d, 0) && LtQ(~g / ~e, 0))))))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~p * (~f + ~(g') * ~x) ^ ~q * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => ((-((~d) ^ 2) * ~g) / ~e) ^ ~q * Antiderivative((~d + ~e * ~x) ^ (~p - ~q) * (1 - (~c) ^ 2 * (~x) ^ 2) ^ ~q * (~a + ~b * acos(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && (EqQ(~e * ~f + ~d * ~g, 0) && (EqQ((~c) ^ 2 * (~d) ^ 2 - (~e) ^ 2, 0) && (HalfIntegerQ(~p, ~q) && (GeQ(~p - ~q, 0) && (GtQ(~d, 0) && LtQ(~g / ~e, 0))))))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~p * (~f + ~(g') * ~x) ^ ~q * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => (((~d + ~e * ~x) ^ ~q * (~f + ~g * ~x) ^ ~q) / (1 - (~c) ^ 2 * (~x) ^ 2) ^ ~q) * Antiderivative((~d + ~e * ~x) ^ (~p - ~q) * (1 - (~c) ^ 2 * (~x) ^ 2) ^ ~q * (~a + ~b * asin(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && (EqQ(~e * ~f + ~d * ~g, 0) && (EqQ((~c) ^ 2 * (~d) ^ 2 - (~e) ^ 2, 0) && (HalfIntegerQ(~p, ~q) && GeQ(~p - ~q, 0))))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~p * (~f + ~(g') * ~x) ^ ~q * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => (((~d + ~e * ~x) ^ ~q * (~f + ~g * ~x) ^ ~q) / (1 - (~c) ^ 2 * (~x) ^ 2) ^ ~q) * Antiderivative((~d + ~e * ~x) ^ (~p - ~q) * (1 - (~c) ^ 2 * (~x) ^ 2) ^ ~q * (~a + ~b * acos(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && (EqQ(~e * ~f + ~d * ~g, 0) && (EqQ((~c) ^ 2 * (~d) ^ 2 - (~e) ^ 2, 0) && (HalfIntegerQ(~p, ~q) && GeQ(~p - ~q, 0))))
end

max_stars_repo_path: src/rules/5 Inverse trig functions/5.1 Inverse sine/5.1.3 (d+e x^2)^p (a+b arcsin(c x))^n.jl


text:
using GSW
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
    as = 35.0; ct = 10.0; pr = 10.0
    println(GSW.gsw_rho(as, ct, pr))
end

# write your own tests here
#@test 1 == 2


text:
# First, document and link test case

export vdpol

vdpol = let
    tcname = :vdpol
    T = Float64 # the datatype used
    Tarr = Matrix
    dof = 2 # degrees of freedom
    dae=0
    # stiffness of system, one of the three constants
    stiffness = [nonstiff, mildlystiff, stiff][3]

    mu = 1e3
    eps = 1/mu^2 # rescaled parameter
    function fn!(t,y,dydt)
        # the ode function dydt = f(t,y)
        dydt[1] = y[2]
        prod =
        dydt[2] = ( (1.0 - y[1]^2)*y[2] - y[1]) / eps
        return nothing
    end
    fn!( ; T_::Type=T, dof_=dof) = zeros(T_,dof_)

    function jac!(t,y,dfdy)
        # the Jacobian of f
        #
        # Set to nothing if it is not known
        dfdy[1,1] = 0.0
        dfdy[1,2] = 1.0
        dfdy[2,1] = (-2.0 * y[1]*y[2]-1.0)/eps
        dfdy[2,2] = (1.0 - y[1]^2)/eps
        return nothing
    end
    jac!( ; T_::Type=T, dof_=dof) = zeros(T_,dof_,dof_)  # if the problem is large better return an appropriate sparse matrix

    mass! = nothing
    ic = T[2, 0] # vector of initial conditions
    tspan = T[0, 2] # integration interval
    refsol = T[0.1706167732170483e1, -0.8928097010247975e0] # reference solution at tspan[2]
    refsolinds = trues(dof)
    scd_absinds = falses(dof)
    tc = TestCaseExplicit{tcname, T, Tarr}(
                             stiffness,
                             dae,
                             dof,
                             fn!,
                             jac!,
                             mass!,
                             ic,
                             tspan,
                             refsol,
                             refsolinds,
                             scd_absinds)
    tc_all[tcname] = tc
    tc_stiff[tcname] = tc
end


text:
module StippleCharts

import Genie
import Stipple

using Stipple.Reexport

#===#

const COMPONENTS = [:apexchart => :VueApexCharts]

#===#

const assets_config = Genie.Assets.AssetsConfig(package = "StippleCharts.jl")

function deps() :: String
  if ! Genie.Assets.external_assets(Stipple.assets_config)

    Genie.Router.route(Genie.Assets.asset_path(assets_config, :js, file="apexcharts.min")) do
      Genie.Renderer.WebRenderable(
        Genie.Assets.embedded(Genie.Assets.asset_file(cwd=normpath(joinpath(@__DIR__, "..")), file="apexcharts.min.js")),
        :javascript) |> Genie.Renderer.respond
    end

    Genie.Router.route(Genie.Assets.asset_path(assets_config, :js, file="vue-apexcharts.min")) do
      Genie.Renderer.WebRenderable(
        Genie.Assets.embedded(Genie.Assets.asset_file(cwd=normpath(joinpath(@__DIR__, "..")), file="vue-apexcharts.min.js")),
        :javascript) |> Genie.Renderer.respond
    end

  end

  string(
    Genie.Renderer.Html.script(src="$(Genie.Assets.asset_path(assets_config, :js, file="apexcharts.min"))"),
    Genie.Renderer.Html.script(src="$(Genie.Assets.asset_path(assets_config, :js, file="vue-apexcharts.min"))")
  )
end

#===#

include("Charts.jl")
@reexport using .Charts

function __init__()
  Stipple.DEPS[@__MODULE__] = deps
end

end # module


text:
### Options

type Options
    maxiter::Int        # maximum number of iterations
    ftol::Float64       # function value change tolerance
    xtol::Float64       # solution change tolerance
    grtol::Float64      # gradient norm tolerance
    armijo::Float64     # Armijo coefficient for line search
    beta::Float64       # backtracking ratio
    verbosity::Symbol   # verbosity (:none | :final | :iter)
end

function Options(;maxiter::Integer=200,
                  ftol::Real=1.0e-6,
                  xtol::Real=1.0e-8,
                  grtol::Real=1.0e-8,
                  armijo::Real=0.5,
                  beta::Real=0.5,
                  verbosity::Symbol=:none)

     maxiter > 1 || error("maxiter must be an integer greater than 1.")
     ftol > 0 || error("ftol must be a positive real value.")
     xtol > 0 || error("xtol must be a positive real value.")
     grtol > 0 || error("grtol must be a positive real value.")
     0 < armijo < 1 || error("armijo must be a real value in (0, 1).")
     0 < beta < 1 || error("beta must be a real value in (0, 1).")
     (verbosity == :none || verbosity == :final || verbosity == :iter) ||
         error("verbosity must be either :none, :final, or :iter.")

     Options(convert(Int, maxiter),
             convert(Float64, ftol),
             convert(Float64, xtol),
             convert(Float64, grtol),
             convert(Float64, armijo),
             convert(Float64, beta),
             verbosity)
end

### Objective Function

abstract Functional{T}

immutable RiskFun{T,
                  XT<:StridedArray,
                  YT<:StridedArray,
                  RModel<:RiskModel} <: Functional{T}

    rmodel::RModel
    X::XT
    Y::YT
end

RiskFun{T<:AbstractFloat}(rmodel::RiskModel, X::StridedArray{T}, Y::StridedArray) =
    RiskFun{T, typeof(X), typeof(Y), typeof(rmodel)}(rmodel, X, Y)

value{T<:AbstractFloat}(f::RiskFun{T}, θ::StridedArray{T}) = value(f.rmodel, θ, f.X, f.Y)

value!{T<:AbstractFloat}(buffer, f::RiskFun{T}, θ::StridedArray{T}) =
    value!(buffer, f.rmodel, θ, f.X, f.Y)

value_and_grad!{T<:AbstractFloat}(f::RiskFun{T}, g::StridedArray{T}, θ::StridedArray{T}) =
    value_and_addgrad!(f.rmodel, zero(T), g, one(T), θ, f.X, f.Y)

value_and_grad!{T<:AbstractFloat}(buffer, f::RiskFun{T}, g::StridedArray{T}, θ::StridedArray{T}) =
    value_and_addgrad!(buffer, f.rmodel, zero(T), g, one(T), θ, f.X, f.Y)

immutable RegRiskFun{T,
                     XT<:StridedArray,
                     YT<:StridedArray,
                     RModel<:RiskModel,
                     Reg<:Regularizer} <: Functional{T}
    rmodel::RModel
    reg::Reg
    X::XT
    Y::YT
end

RegRiskFun{T<:AbstractFloat}(rmodel::RiskModel, reg::Regularizer, X::StridedArray{T}, Y::StridedArray) =
    RegRiskFun{T, typeof(X), typeof(Y), typeof(rmodel), typeof(reg)}(rmodel, reg, X, Y)

value{T<:AbstractFloat}(f::RegRiskFun{T}, θ::StridedArray{T}) =
    value(f.rmodel, θ, f.X, f.Y) + value(f.reg, θ)

value!{T<:AbstractFloat}(buffer, f::RegRiskFun{T}, θ::StridedArray{T}) =
    value!(buffer, f.rmodel, θ, f.X, f.Y) + value(f.reg, θ)

function value_and_grad!{T<:AbstractFloat}(f::RegRiskFun{T}, g::StridedArray{T}, θ::StridedArray{T})
    v_risk, _ = value_and_addgrad!(f.rmodel, zero(T), g, one(T), θ, f.X, f.Y)
    v_regr, _ = value_and_addgrad!(f.reg, one(T), g, one(T), θ)
    return (v_risk + v_regr, g)
end

function value_and_grad!{T<:AbstractFloat}(buffer, f::RegRiskFun{T}, g::StridedArray{T}, θ::StridedArray{T})
    v_risk, _ = value_and_addgrad!(buffer, f.rmodel, zero(T), g, one(T), θ, f.X, f.Y)
    v_regr, _ = value_and_addgrad!(f.reg, one(T), g, one(T), θ)
    return (v_risk + v_regr, g)
end

### Solution

immutable Solution{Sol<:StridedArray}
    sol::Sol
    fval::Float64
    niters::Int
    converged::Bool
end

function Base.show(io::IO, r::Solution)
    println(io, "RiskMinSolution:")
    println(io, "- sol:       $(size(r.sol)) $(typeof(r.sol))")
    println(io, "- fval:      $(r.fval)")
    println(io, "- niters:    $(r.niters)")
    println(io, "- converged: $(r.converged)")
end

### Solver

abstract Solver
abstract DescentSolver <: Solver

### Line search

function backtrack!{T<:AbstractFloat}(
    buffer,
    f::Functional{T},       # objective function
    θr::Array{T},           # destination solution
    θ::Array{T},            # last solution
    v::T,                   # objective value at θ
    g::Array{T},            # gradient at θ
    p::Array{T},            # descent direction
    α::T,                   # initial step size
    opts::Options)          # options

    armijo = convert(T, opts.armijo)
    β = convert(T, opts.beta)
    dv = dot(vec(p), vec(g))
    dv > zero(T) || error("The descent direction is invalid.")

    _xmcy!(θr, θ, α, p)   # θr <- θ - α p
    v2 = value!(buffer, f, θr)
    while v2 > v - armijo * α * dv
        α > eps(T) || error("Failed to find a proper step size.")
        α *= β
        _xmcy!(θr, θ, α, p)   # θr <- θ - α p
        v2 = value!(buffer, f, θr)
    end
    return α
end

function prox_backtrack!{T<:AbstractFloat}(
    f::Functional{T},       # smooth-part of the objective function
    reg::Regularizer,       # non-smoothed part (regularizer)
    θr::Array{T},           # destination solution
    θ::Array{T},            # last solution
    vf::T,                  # objective value at θ, w.r.t. f
    vr::T,                  # regularization value
    g::Array{T},            # gradient at θ, w.r.t. f
    p::Array{T},            # descent direction
    α::T,                   # initial step size
    opts::Options)          # options

    # debug
    # g_ = similar(g)
    # vf_, _ = value_and_grad!(f, g_, θ)
    # vr_ = value(reg, θ)
    # @assert vecnorm(g - g_) < 1.0e-12
    # @assert abs(vf - vf_) < 1.0e-12
    # @assert abs(vr - vr_) < 1.0e-12
    # @assert vecnorm(p - g) < 1.0e-12

    armijo = opts.armijo
    β = convert(T, opts.beta)
    dv = dot(vec(p), vec(g))
    sqr_pnorm = sumabs2(p)
    vr = value(reg, θ)
    # @show dv

    found = false
    while !found
        _xmcy!(θr, θ, α, p)  # θr <- θ - α p
        prox!(reg, θr, θr, α)
        vf2 = value(f, θr)
        vr2 = value(reg, θr)
        # τ = vf - α * dv + (α/2) * sqr_pnorm + vr
        if vf2 + vr2 < vf + vr
            found = true
        end
        if !found
            if α > eps(T)
                α *= β
            else
                error("Failed to find a proper step size.")
            end
        end
    end

    return α
end

### test convergence

test_convergence{T<:AbstractFloat}(θ::Array{T}, θpre::Array{T}, v::T, vpre::T, g::Array{T}, opt::Options) =
    abs(v - vpre) < convert(T, opt.ftol) ||
    vecnorm(g) < convert(T, opt.grtol) ||
    _l2diff(θ, θpre) < convert(T, opt.xtol)

test_convergence{T<:AbstractFloat}(θ::Array{T}, θpre::Array{T}, v::T, vpre::T, opt::Options) =
    abs(v - vpre) < convert(T, opt.ftol) ||
    _l2diff(θ, θpre) < convert(T, opt.xtol)

### auxiliary functions

function _sqrl2diff{T<:AbstractFloat}(x::StridedArray{T}, y::StridedArray{T})
    @assert length(x) == length(y)
    s = zero(T)
    @inbounds for i = 1:length(x)
        s += abs2(x[i] - y[i])
    end
    return s
end

_l2diff{T<:AbstractFloat}(x::StridedArray{T}, y::StridedArray{T}) =
    sqrt(_sqrl2diff(x, y))

function _xmcy!{T<:Real}(r::Array{T}, x::Array{T}, c::T, y::Array{T})
    @inbounds for i = 1:length(x)
        r[i] = x[i] - c * y[i]
    end
    r
end


text:
with(File(...), Paper(...), Ink(...), Pen(...)) do
    draw(Ink("red")) do
        move(...)
        line(...)
        line(...)
        arc(...)
    end
    paint(Ink("blue")) do
        line(...)
        line(...)
        start()
    end
end


text:
# 3D Grid examples
# ===============
#

# ## Quadrilateral
function quadrilateral()
    X=collect(0:0.25:1)
    Y=collect(0:0.2:1)
    Z=collect(0:0.1:1)
    simplexgrid(X,Y,Z)
end
# ![](quadrilateral.svg)


text:
# --- Calling Into MATLAB from embedded Julia --- #

# Call a matlab function specified by name
function call_matlab(nout::Integer, fn::String, args::Vector{MATLAB.MxArray})

    # initialize outputs
    outs = fill(C_NULL, nout)

    # call matlab
    ptr = ccall(mex_call_matlab_with_trap[], Ptr{Cvoid},
        (Int32, Ptr{Ptr{Cvoid}}, Int32, Ptr{Ptr{Cvoid}}, Ptr{UInt8}),
        nout, outs, length(args), args, fn)

    # check for errors
    if ptr != C_NULL
        throw(MatlabException(ptr))
    end

    # return dressed up outputs
    return MATLAB.MxArray.(outs, false)
end

# this version handles argument conversions to and from MATLAB
call_matlab(nout::Integer, fn::String, args...) = MATLAB.jvalue.(call_matlab(nout, fn, MATLAB.mxarray.(collect(args))))

# Make MxArray callable. Works for strings or function handles.
(mx::MATLAB.MxArray)(nout::Integer, args::Vector{MATLAB.MxArray}) = call_matlab(nout, "feval", vcat(mx, args))
(mx::MATLAB.MxArray)(args...) = MATLAB.jvalue(mx(1, MATLAB.mxarray.(args))[1])


text:
using NPZ
using CSV
using JSON
using Plots
using DataFrames
using Statistics
using CartesianGeneticProgramming

include("utils.jl")

struct Month
    topo0
    topo1
    tp
    hs
    direction
end

gene_json = JSON.parsefile("path/to/gene")
fitness = gene_json["fitness"]
functions = gene_json["fs"]
chromosome = convert(Array{Float64,1}, gene_json["chromosome"])
xs = gene_json["xs"]
ys = gene_json["ys"]

ind = CGPInd(get_config("topo.yaml"), chromosome)

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

dataset_path = get_dataset_path()
n_months = 236
months = []

for i in 1:n_months
    topo_path = string(dataset_path, '/', i, "_topo.npy")
    topo = npzread(topo_path)
    topo0 = norm_to_range(topo[1, :], -10.0, 10.0, -1, 1)
    topo1 = norm_to_range(topo[2, :], -10.0, 10.0, -1, 1)

    waves_path = string(dataset_path, '/', i, "_waves.csv")
    waves_csv = CSV.read(waves_path, DataFrame)

    tp_month = norm_to_range(waves_csv.tp, 0, 25, -1, 1)
    hs_month = norm_to_range(waves_csv.hs, 0, 10, -1, 1)
    direction_month = norm_to_range(waves_csv.direction, 0, 360, -1, 1)

    tp_weeks = []
    hs_weeks = []
    direction_weeks = []

    for i in range(0, stop=27)
        tp = tp_month[(i*8) + 1 : (i*8) + 8]
        hs = hs_month[(i*8) + 1 : (i*8) + 8]
        direction = direction_month[(i*8) + 1 : (i*8) + 8]
        push!(tp_weeks, tp)
        push!(hs_weeks, hs)
        push!(direction_weeks, direction)
    end
    month = Month(topo0, topo1, tp_weeks, hs_weeks, direction_weeks)
    global months
    push!(months, month)
end

month = months[1]

topo0 = deepcopy(month.topo0)
for w_i in 1:size(month.tp)[1]
    global topo0

    x_row = vcat(topo0, month.tp[w_i], month.hs[w_i], month.direction[w_i])
    topo0 = process(ind, x_row)
end
println("final plot:")

display(plot(norm_to_range(topo0, -1, 1, -10.0, 10.0)))
readline()

max_stars_repo_name: mahmoud-al-najar/cgp-morphcast

max_issues_repo_name: mahmoud-al-najar/cgp-morphcast

max_forks_repo_name: mahmoud-al-najar/cgp-morphcast


text:
export Weighted
struct Weighted{W,T}
    ℓ :: W
    val :: T
end

using Printf
function Base.show(io::IO, ℓx::Weighted)
    @printf io "Weighted(%g.4\n" (ℓx.ℓ)
    println(",", ℓx.val)
end


text:
using IntervalWavelets
using Base.Test

println("Testing interior scaling functions...")

#=
Orthogonality of integer translates of the internal scaling function
=#

EPS = 1e-3
R = 10
for N in 2:10
	C = ifilter(N)

	# Union of supports of all translated scaling functions
	S = support(C)
	extended_supp = DaubSupport( left(S), 2*right(S)-1 )
	x = dyadic_rationals( extended_supp, R )

	# Base scaling funtion
	y1 = zeros(x)
	phi = DaubScaling(C, R)
	@test_approx_eq_eps IntervalWavelets.l2norm(dyadic_rationals(S,R), phi) 1.0 EPS

	supp_index = [1:length(phi);]
	y1[supp_index] = phi

	# Integer translates
	y2 = similar(y1)
	for k in 1:right(S)-1
		fill!(y2, 0.0)
		y2[k*2^R-1+supp_index] = phi
		#@show IntervalWavelets.inner(x, y1, y2)
		@test_approx_eq_eps IntervalWavelets.inner(x, y1, y2) 0.0 EPS
	end
end


text:
export BufferStream, mem_usage, append_chunk
import Base: write, read, readbytes!, readavailable, unsafe_write, wait, close, eof, isopen, skip, length

"""
    BufferStream

Simple buffered stream that appends data to an internal chunk list and allows for easy,
simultaneous reading of that data.  Intended for shuffling data between tasks.  Best
performance is achieved by using `readavailable()`, which avoids copies on reads.
"""
mutable struct BufferStream <: IO
    chunks::Vector{AbstractVector{UInt8}}
    chunk_read_idx::Int
    read_cond::Threads.Condition
    write_cond::Threads.Condition
    is_open::Bool
    max_len::Int

    BufferStream(max_len::Int = 0) = new(AbstractVector{UInt8}[], 1, Threads.Condition(), Threads.Condition(), true, max_len)
end

# Close the stream, writing not allowed (but can still read until `eof()`)
function close(bs::BufferStream)
    bs.is_open = false
    lock(bs.read_cond) do
        notify(bs.read_cond)
    end
end

isopen(bs::BufferStream) = bs.is_open
eof(bs::BufferStream) = !isopen(bs) && isempty(bs.chunks)
function length(bs::BufferStream)
    lock(bs.write_cond) do
        if isempty(bs.chunks)
            return 0
        end
        len = length(bs.chunks[1]) - bs.chunk_read_idx + 1
        for idx in 2:length(bs.chunks)
            len += length(bs.chunks[idx])
        end
        return len
    end
end
function mem_usage(bs::BufferStream)
    lock(bs.write_cond) do
        return sum(Int[length(chunk) for chunk in bs.chunks])
    end
end

function append_chunk(bs::BufferStream, data::AbstractVector{UInt8})
    # Disallow writing if we're not open
    if !isopen(bs)
        throw(ArgumentError("Stream is closed"))
    end

    # Copy the data so that users can't clobber our internal list
    lock(bs.write_cond) do
        data_written = 0
        while data_written < length(data)
            if bs.max_len == 0
                space_available = length(data)
            else
                space_available = bs.max_len - length(bs)
            end
            if space_available == 0
                wait(bs.write_cond)
                continue
            end
            bytes_to_write = min(space_available, length(data) - data_written)
            push!(bs.chunks, data[data_written+1:data_written+bytes_to_write])
            # Notify someone who was waiting for some data
            lock(bs.read_cond) do
                notify(bs.read_cond; all=false)
            end
            data_written += bytes_to_write
        end
    end

    return length(data)
end

# Helper methods for read/write
write(bs::BufferStream, x::UInt8) = write(bs, [x])
function unsafe_write(bs::BufferStream, ref::Ptr{UInt8}, nbytes::UInt)
    return append_chunk(bs, copy(unsafe_wrap(Array, ref, (nbytes,))))
end

# Read a single byte.
function read(bs::BufferStream, ::Type{UInt8})
    data = UInt8[0]
    readbytes!(bs, data, 1)
    return data[1]
end

# Read everything
function read(bs::BufferStream)
    ret = UInt8[]
    lock(bs.read_cond) do
        while !eof(bs)
            append!(ret, readavailable(bs))
        end
    end
    return ret
end

# Completely consume the first chunk.
function readavailable(bs::BufferStream)
    lock(bs.read_cond) do
        if isempty(bs.chunks)
            if !isopen(bs)
                return UInt8[]
            end
            wait(bs.read_cond)

            # Handle cancellation explicitly
            if isempty(bs.chunks)
                return UInt8[]
            end
        end

        # We're gonna consume the rest of this chunk
        local chunk
        lock(bs.write_cond) do
            chunk = popfirst!(bs.chunks)
            notify(bs.write_cond; all=false)
        end
        chunkview = view(chunk, bs.chunk_read_idx:length(chunk))
        bs.chunk_read_idx = 1
        return chunkview
    end
end

# Read up to `maxlen` bytes, potentially partially consuming the first chunk.
function readbytes!(bs::BufferStream, data::Vector{UInt8}, maxlen::Int)
    lock(bs.read_cond) do
        if isempty(bs.chunks)
            if !isopen(bs)
                return 0
            end
            wait(bs.read_cond)

            # Handle cancellation explicitly
            if isempty(bs.chunks)
                return 0
            end
        end

        # Grab the first chunk that is available
        local chunk
        lock(bs.write_cond) do
            firstchunk = first(bs.chunks)

            # take view of only what we will read
            chunk = view(firstchunk, bs.chunk_read_idx:min(bs.chunk_read_idx + maxlen - 1,length(firstchunk)))

            # Update datastructure to consume the view we just took
            if maxlen > length(firstchunk) - bs.chunk_read_idx
                popfirst!(bs.chunks)
                bs.chunk_read_idx = 1
                notify(bs.write_cond; all=false)
            else
                bs.chunk_read_idx += maxlen
            end
        end

        for idx in 1:length(chunk)
            data[idx] = chunk[idx]
        end
        return length(chunk)
    end
end

function skip(bs::BufferStream, n::Int)
    data = Array{UInt8}(undef, min(2*1024*1024, n))
    while n > 0
        n -= readbytes!(bs, data, min(length(data), n))
        if eof(bs)
            break
        end
    end
    return n
end

max_stars_repo_name: staticfloat/SimpleBufferStream

max_issues_repo_name: staticfloat/SimpleBufferStream

max_forks_repo_name: staticfloat/SimpleBufferStream


text:
"""
    add_obj_fun!(
        ata_model::SoysterMaximinModel;
        psychometrics = false,
        items::Vector{Psychometrics.Item} = Psychometrics.Item[],
        items_file = "",
        kwargs...
    )

# Description

It adds the objective function as specified in the `settings_file`. It requires the [`start_ata`](#ATA.start_ata) build step.

Computes the IIFs at predefined ability points using the sampled item parameter estimates.
For each pair of item and ability point \$(i, θ_k)\$, the IIF\$(θ_k)_i\$ is computed as the mean minus the standard deviation
of the IIF\$(θ_k)_{ir}\$ computed on the R sampled item parameter estimates, where \$r = 1, \\ldots, R\$.

The vector of items can be passed either through the argument `items` (`Vector{Psychometrics.Item}`) or through the argument `items_file` (string file name).
In both cases, the items in the vector must have the same order of the items in the item bank and they must contain the R sampled item parameters in the field `parameters.chain`.

- **`ata_model::DeJongMaximinModel`** : Required.
"""
function add_obj_fun!(
    ata_model::SoysterMaximinModel;
    psychometrics = false,
    items::Vector{Psychometrics.Item} = Psychometrics.Item[],
    items_file = "",
    kwargs...,
)

    try
        if psychometrics
            soyster_load_parameters_chain!(
                ata_model;
                items_file = items_file,
                items = items,
                kwargs...,
            )
            success!(ata_model, "IIFs for all item parameters samples computed.")

        end
    catch e
        error!(ata_model, string(sprint(showerror, e)))

    end
    return nothing
end


text:
using Distributed

while nprocs() < Sys.CPU_THREADS
    addprocs(1; exeflags="--project")
    @info "adding worker"
end

@everywhere using HTTP, Sockets, Dates, PyCall
@everywhere @pyimport textblob
@everywhere @pyimport textblob.np_extractors as npx

@info "creating worker pool"
global wp = WorkerPool(workers())

@everywhere begin
    include("./src/solrplant_api.jl")
end
# SERVER FUNCTION
function process_text(;text="")
    @info "Processing Text"
    extractor = npx.ConllExtractor()
    f = remotecall(extract_plants, wp, text, extractor)
    @info "Result Text"
    @show f
    return fetch(f)
end

function build_server()
    headers = Dict{AbstractString,AbstractString}(
        "Server"            => "Julia/$VERSION",
        "Content-Type"      => "text/plain; charset=utf-8",
        "Content-Language"  => "en",
        "Date"              => Dates.format(now(Dates.UTC), Dates.RFC1123Format),
        # "Access-Control-Allow-Origin" => "https://bcbi.brown.edu",
        "Access-Control-Allow-Origin" => "*",
        "Access-Control-Allow-Methods" => "POST, GET" )

    # route handlers
    h_text = HTTP.Handlers.HandlerFunction( (req) -> begin

            @info "Entered Handler"
            @show headers
            @info "Request"
            @show req

            try
                @info "Received text"
                text = String(req.body)
                @show text
                response = HTTP.Response(200, HTTP.Headers(collect(headers)), body = process_text(text=text))
                return response
            catch
                @warn "Incorrect request"
                return HTTP.Response(400, HTTP.Headers(collect(headers)), body = "Incorrect request")
            end

        end)

    r = HTTP.Router()

    HTTP.register!(r, "POST", "/", h_text)

    return HTTP.Servers.Server(r, ratelimit=typemax(Int64)//1)

end

@info "building server"
s = build_server()

@info "running server"
port = 8081
HTTP.Servers.serve(s, "0.0.0.0", port)

max_stars_repo_name: bcbi/SolrPlantAPI

max_issues_repo_name: bcbi/SolrPlantAPI

max_forks_repo_name: bcbi/SolrPlantAPI


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ d1e8bbf9-c687-42da-88e4-95d731fc85ad
md"""# 熟悉Julia"""

# ╔═╡ a94472c3-8a03-4398-b8c0-1cf71ee4152d
md"""## 安装Julia 1.6.2

到[Julia的官方站点](https://julialang.org/)下载最新的Julia（目前版本是1.6.2），我们可以使用[适用于Windows的x64免安装版本](https://julialang-s3.julialang.org/bin/winnt/x64/1.6/julia-1.6.2-win64.zip)。也可以在国内的镜像站点下载，比如直接从[清华的镜像](https://mirrors.tuna.tsinghua.edu.cn/julia-releases/bin/winnt/x64/1.6/julia-1.6.2-win64.zip)下载。下载完后解压到合适的目录，比如 D:\soft_using\julia-1.6.2。

接下来，我们*需要把julia的可执行文件julia.exe的目录 D:\soft_using\julia-1.6.2\bin 加到环境变量PATH中去*。如果不清楚什么是环境变量或者什么是PATH，请[参阅此文](https://www.cnblogs.com/qtiger/p/13903609.html)。
"""

# ╔═╡ f5d320e6-e8bd-455d-9ce9-f2a9e37ffc29
md"""## 安装VSCode
Visual Studio Code(简称VSCode/VSC)是当今最流行的免费开源代码编辑器，由微软出品。它支持几乎所有主流的开发语言的语法高亮、智能代码补全、自定义热键、括号匹配、代码片段、代码对比Diff、Git等特性，支持插件扩展，并针对网页开发和云端应用开发做了优化。软件跨平台支持 Win、Mac 以及 Linux，运行流畅。

到[VSCode的官方站点](https://code.visualstudio.com/)下载最新的安装包，根据向导安装即可。可把vscode加到文件夹和文件的右键菜单中，便于使用。接下来给vscode装上Julia插件。关于VSCode的入门操作，可以参考[菜鸟教程上的VScode教程](https://www.runoob.com/w3cnote/vscode-tutorial.html)。"""

# ╔═╡ d4b9d09d-58d1-4d3b-8754-180b27874b6a
md"""## Julia安装Pluto包
首先参考[Julia中文文档](https://docs.juliacn.com/latest/manual/getting-started/)了解一些Julia的入门基础知识。

在julia提示符下，通过以下命令安装Pluto包：

```
using Pkg
Pkg.add("Pluto")
```
"""

# ╔═╡ 48e4cc4c-eec8-4768-a7b4-339607b323bc
md"""
在Julia提示符下通过如下命令启动Pluto：

```
using Pluto
Pluto.run()```

通常系统会自动调用默认的浏览器打开一个页面，然后就开始Pluto之旅吧。
"""

# ╔═╡ 52b492d0-85c6-4e21-8efd-7d9a25f06cde
md"""# 测试一下

我们这里简要测试一下Julia的命令，看看系统运行的样子。"""

# ╔═╡ f8240954-d993-4229-a438-96e00ef461ce
x=2.0

# ╔═╡ 6811adfd-6403-427c-97ea-47cb0b06ee33
sin(x)

# ╔═╡ 90608d89-2b8a-4efb-99f2-a3f0fa13b40d
@bind y html"<input type=range>"

# ╔═╡ 1d59eb51-24a2-49d2-8d11-fdd681e417bf
md"""我们通过以上这个命令，把y与拖拉杆绑定，下面我们使用一下这个变量y试试，并且计算一下sin（y）看看"""

# ╔═╡ 28d55a1e-7f63-4505-8ab7-cb475eef55f3
y

# ╔═╡ b105c8da-174b-4820-8bf1-55d855ba0edd
sin(y)

# ╔═╡ Cell order:
# ╟─d1e8bbf9-c687-42da-88e4-95d731fc85ad
# ╟─a94472c3-8a03-4398-b8c0-1cf71ee4152d
# ╟─f5d320e6-e8bd-455d-9ce9-f2a9e37ffc29
# ╟─d4b9d09d-58d1-4d3b-8754-180b27874b6a
# ╠═48e4cc4c-eec8-4768-a7b4-339607b323bc
# ╟─52b492d0-85c6-4e21-8efd-7d9a25f06cde
# ╠═f8240954-d993-4229-a438-96e00ef461ce
# ╠═6811adfd-6403-427c-97ea-47cb0b06ee33
# ╠═90608d89-2b8a-4efb-99f2-a3f0fa13b40d
# ╟─1d59eb51-24a2-49d2-8d11-fdd681e417bf
# ╠═28d55a1e-7f63-4505-8ab7-cb475eef55f3
# ╠═b105c8da-174b-4820-8bf1-55d855ba0edd


text:
using ModelingToolkit, SciMLBase, Serialization

@parameters t
@variables x(t)
D = Differential(t)

sys = ODESystem([D(x) ~ -0.5*x], defaults=Dict(x=>1.0))
for prob in [
    eval(ModelingToolkit.ODEProblem{false}(sys, nothing, nothing, SciMLBase.NullParameters())),
    eval(ModelingToolkit.ODEProblemExpr{false}(sys, nothing, nothing, SciMLBase.NullParameters()))
]
    _fn = tempname()

    open(_fn, "w") do f
        serialize(f, prob)
    end

    _cmd = "using ModelingToolkit, Serialization; deserialize(\"$_fn\")"

    run(`$(Base.julia_cmd()) -e $(_cmd)`)
end

include("../examples/rc_model.jl")
io = IOBuffer()
write(io, rc_model)
sys = include_string(@__MODULE__, String(take!(io)))
@test sys == flatten(rc_model)


text:
import ArrayInterface:
    ArrayInterface,
    StaticInt,
    contiguous_axis,
    contiguous_batch_size,
    parent_type,
    stride_rank,
    dense_dims

parent_type(::Type{<:PencilArray{T,N,A}}) where {T,N,A} = A

contiguous_axis(::Type{A}) where {A <: PencilArray} =
    _contiguous_axis(
        contiguous_axis(parent_type(A)),
        permutation(A),
    )

_contiguous_axis(x::Nothing, ::AbstractPermutation) = x
_contiguous_axis(x::StaticInt, ::NoPermutation) = x
@inline function _contiguous_axis(x::StaticInt{i}, p::Permutation) where {i}
    i == -1 && return x
    StaticInt(p[Val(i)])
end

contiguous_batch_size(::Type{A}) where {A <: PencilArray} =
    contiguous_batch_size(parent_type(A))

function stride_rank(::Type{A}) where {A <: PencilArray}
    rank = stride_rank(parent_type(A))
    rank === nothing && return nothing
    iperm = Tuple(inv(permutation(A)))
    iperm === nothing && return rank
    ArrayInterface.permute(rank, Val(iperm))
end

function dense_dims(::Type{A}) where {A <: PencilArray}
    dense = dense_dims(parent_type(A))
    dense === nothing && return nothing
    perm = Tuple(inv(permutation(A)))
    perm === nothing && return dense
    ArrayInterface.permute(dense, Val(perm))
end

ArrayInterface.size(A::PencilArray) =
    permutation(A) * ArrayInterface.size(parent(A))

ArrayInterface.strides(A::PencilArray) =
    permutation(A) * ArrayInterface.strides(parent(A))


text: isupper(str::String) = all( c -> isuppercase(c) || isspace(c), str )

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
#### Libs
using Statistics

#### Constants
const MODEL_FILENAME = "pls_model.jld" # jld filename for storing the model
const MODEL_ID       = "pls_model"     # if od the model in the filesystem jld data

#### An abstract pls model
abstract type PLSModel{T} end

#### PLS1 type
mutable struct PLS1Model{T<:AbstractFloat} <:PLSModel{T}
    W::Matrix{T}           # a set of vectors representing correlation weights of input data (X) with the target (Y)
    b::Matrix{T}           # a set of scalar values representing a latent value for dependent variables or target (Y)
    P::Matrix{T}           # a set of latent vetors for the input data (X)
    nfactors::Int          # a scalar value representing the number of latent variables
    mx::Matrix{T}          # mean stat after for z-scoring input data (X)
    my::T                  # mean stat after for z-scoring target data (Y)
    sx::Matrix{T}          # standard deviation stat after z-scoring input data (X)
    sy::T                  # standard deviation stat after z-scoring target data (X)
    nfeatures::Int         # number of input (X) features columns
    centralize::Bool       # store information of centralization of data. if true, tehn it is passed to transform function
end

## PLS1: constructor
function Model(X::Matrix{T},
               Y::Vector{T},
               nfactors::Int,
               centralize::Bool) where T<:AbstractFloat
    (nrows,ncols) = size(X)
    ## Allocation
    return PLS1Model(zeros(T,ncols,nfactors), ## W
            zeros(T,1,nfactors),       ## b
            zeros(T,ncols,nfactors), ## P
            nfactors,
            mean(X,dims=1),
            mean(Y),
            std(X,dims=1),
            std(Y),
            ncols,
            centralize)
end

########################################################################################
#### PLS2 type
mutable struct PLS2Model{T<:AbstractFloat} <:PLSModel{T}
    W::Matrix{T}           # a set of vectors representing correlation weights of input data (X) with the target (Y)
    Q::Matrix{T}           #
    b::Matrix{T}           # a set of scalar values representing a latent value for dependent variables or target (Y)
    P::Matrix{T}           # a set of latent vetors for the input data (X)
    nfactors::Int          # a scalar value representing the number of latent variables
    mx::Matrix{T}          # mean stat after for z-scoring input data (X)
    my::Matrix{T}          # mean stat after for z-scoring target data (Y)
    sx::Matrix{T}          # standard deviation stat after z-scoring input data (X)
    sy::Matrix{T}          # standard deviation stat after z-scoring target data (X)
    nfeatures::Int         # number of input (X) features columns
    ntargetcols::Int       # number of target (Y) columns
    centralize::Bool       # store information of centralization of data. if true, tehn it is passed to transform function
end

## PLS2: constructor
function Model(X::Matrix{T},
        Y::Matrix{T}, # this is the diference from PLS1 param constructor!
        nfactors::Int,
        centralize::Bool) where T<:AbstractFloat
    (nrows,ncols) = size(X)
    (n,m)         = size(Y)
    ## Allocation
    return PLS2Model(zeros(T,ncols,nfactors), ## W
            zeros(T,m,nfactors),       ## Q
            zeros(T,n,nfactors),       ## b
            zeros(T,ncols,nfactors),   ## P
            nfactors,
            mean(X,dims=1),
            mean(Y,dims=1),
            std(X,dims=1),
            std(Y,dims=1),
            ncols,
            m,
            centralize)
end

################################################################################
#### KPLS type
mutable struct KPLSModel{T<:AbstractFloat} <:PLSModel{T}
    X::Matrix{T}           # Training set
    K::Matrix{T}           # Kernel matrix
    B::Matrix{T}           # Regression matrix
    nfactors::Int          # a scalar value representing the number of latent variables
    mx::Matrix{T}          # mean stat after for z-scoring input data (X)
    my::AbstractArray{T}   # mean stat after for z-scoring target data (Y)
    sx::Matrix{T}          # standard deviation stat after z-scoring input data (X)
    sy::AbstractArray{T}   # standard deviation stat after z-scoring target data (X)
    nfeatures::Int         # number of input (X) features columns
    ntargetcols::Int       # number of target (Y) columns
    centralize::Bool       # store information of centralization of data. if true, tehn it is passed to transform function
    kernel::AbstractString
    width::Float64
end

## KPLS: constructor
function Model(X::Matrix{T},
            Y::AbstractArray{T}, # this is the diference from PLS1 param constructor!
            nfactors::Int,
            centralize::Bool,
            kernel::String,
            width::Float64) where T<:AbstractFloat
    (nrows,ncols) = size(X)
    (n,m)         = size(Y[:,:])
    ## Allocation
    return KPLSModel(zeros(T,nrows,ncols), ## X
            zeros(T,nrows,nrows),          ## K
            zeros(T,ncols,m),              ## B
            nfactors,
            mean(X,dims=1),
            mean(Y,dims=1),
            std(X,dims=1),
            std(Y,dims=1),
            ncols,
            m,
            centralize,
            kernel,
            width)
end

######################################################################################################
## Load and Store models (good for production)
function load(; filename::AbstractString = MODEL_FILENAME, modelname::AbstractString = MODEL_ID)
    local M
    jldopen(filename, "r") do file
        M = read(file, modelname)
    end
    M
end

function save(M::PLSModel; filename::AbstractString = MODEL_FILENAME, modelname::AbstractString = MODEL_ID)
    jldopen(filename, "w") do file
        write(file, modelname, M)
    end
end


text:
export omega, omega_unscaled

# Tools to compute the sparsity pattern

function omega(h, M, Y, X, Ny, Nx, ϵY, L::LinearTransform)
    Ne = size(X, 2)
    Omega = zeros(Nx+Ny, Nx+Ny)

    precisionY = (1/ϵY^2)
    sqprecisionY = precisionY^2

    # Fill block
    @showprogress for i=1:Ne
        Xi  = X[:, i]

        # Compute Jacobian of the forward model
        J = Jacobian(h, Xi, Ny, Nx)

        # Compute Hessian of the forward model
        H = Hessian(h, Xi, Ny, Nx)

        # Compute Hessian of the log pdf of the transport map
        Hlog = hess_x_log_pdf(M, X[:,i:i])

        # Compute Hessian of the log pdf of the transport map
        # Fill Ωxy block
        @inbounds for k=1:Ny
            for j=1:Nx
                Omega[Ny+j, k] += (J[k,j]*L.L.diag[k]*L.L.diag[Ny+j])^2
            end
        end

        # Fill Ωxx
        @inbounds for k=1:Nx
            for j=1:Nx
                Jj = view(J,:,j)
                Jk = view(J,:,k)
                Hjk = view(H,:,j,k)
                # @show L.L.diag[Ny+k]*L.L.diag[Ny+j]*precisionY*dot(Jj, Jk)
                # @show Hlog[1,j,k]
                # @show (L.L.diag[Ny+k]*L.L.diag[Ny+j])^2*dot(Hjk, Hjk)*precisionY
                Omega[Ny+j, Ny+k] += (-L.L.diag[Ny+k]*L.L.diag[Ny+j]*precisionY*dot(Jj, Jk) + Hlog[1,j,k])^2 +
                                  (L.L.diag[Ny+k]*L.L.diag[Ny+j])^2*dot(Hjk, Hjk)*precisionY
            end
        end
    end

    view(Omega, Ny+1:Ny+Nx, 1:Ny) .*= sqprecisionY
    view(Omega, 1:Ny, Ny+1:Ny+Nx) .= view(Omega, Ny+1:Ny+Nx, 1:Ny)'

    rmul!(Omega, 1.0/Ne)

    # Fill Ωyy with the identity
    @inbounds for i=1:Ny
        Omega[i,i] = 1.0
    end

    return Omega
end

function omega_unscaled(h, M, Y, X, Ny, Nx, ϵY, L::LinearTransform)
    Ne = size(X, 2)
    Ω = zeros(Nx+Ny, Nx+Ny)

    precisionY = (1/ϵY^2)
    sqprecisionY = precisionY^2

    # Fill block
    @showprogress for i=1:Ne
        Xi  = X[:, i]

        # Compute Jacobian of the forward model
        J = Jacobian(h, Xi, Ny, Nx)

        # Compute Hessian of the forward model
        H = Hessian(h, Xi, Ny, Nx)

        # Compute Hessian of the log pdf of the transport map
        Hlog = hess_x_log_pdf(M, X[:,i:i])

        # Compute Hessian of the log pdf of the transport map
        # Fill Ωxy block
        @inbounds for k=1:Ny
            for j=1:Nx
                Ω[Ny+j, k] += J[k,j]^2
            end
        end

        # Fill Ωxx
        @inbounds for k=1:Nx
            for j=1:Nx
                Jj = view(J,:,j)
                Jk = view(J,:,k)
                Hjk = view(H,j,k)
                Ω[Ny+j, Ny+k] += (-precisionY*dot(Jj, Jk) + 1/(L.L.diag[Ny+k]*L.L.diag[Ny+j])*Hlog[1,j,k])^2 +
                                 dot(Hjk, Hjk)*precisionY
            end
        end
    end

    view(Ω, Ny+1:Ny+Nx, 1:Ny) .*= sqprecisionY
    view(Ω, 1:Ny, Ny+1:Ny+Nx) .= view(Ω, Ny+1:Ny+Nx, 1:Ny)'

    rmul!(Ω, 1/Ne)

    # Fill Ωyy with the identity
    @inbounds for i=1:Ny
        Ω[i,i] = 1.0
    end

    return Ω
end


text:
import StateSpaceReconstruction:
	GroupSlices.groupslices,
    GroupSlices.firstinds,
    GroupSlices.groupinds

"""
    BinVisits

A data structure helping to organize information about bins that
get visited by the orbit.

## Fields

`first_visited_by`: Vector with the same length as there are visited bins.
`first_visited_by[i]` identifies which point of the orbit first visits the
i-th bin. If `first_visited_by[i] = 5`, then the i-th bin first got visited
by the 5th point of the orbit.

`visitors`. Vector with same length as there are visited bins.
`visitors[i]` will give a vector containing  the (column)
indices of the points visiting that bin. If `visitors[i] = [4, 12]`,
then the i-th bin was visited by points #4 and #12 of the embedding.

`visits_whichbin`: Vector with same length as there are number of points.
Each element of this vector indicate which bin the corresponding
point visits (expressed in terms of the unique bin labels).
If `visits_whichbin[i] = 5`, then the i-th point of the orbit
visits the 5th visited bin.

"""
struct BinVisits
    first_visited_by::Vector{Int}
    visitors::Vector{Vector{Int}}
    visits_whichbin::Vector{Int}
end

####################
# Pretty printing
####################
function summarise(bv::BinVisits)
    _type = typeof(bv)
    L = length(bv.visits_whichbin)
    L_unique = length(bv.visitors)
    "$_type for $L points that occupy $L_unique bins"
end
Base.show(io::IO, bv::BinVisits) = println(io, summarise(bv))

"""
    whichpoints_visit_whichbins(A, U) -> Vector{Int}

Given a sequence of unique elements from an array `U`, say
U = [a b c e], which are taken from the array `A`, say
A = [a b c a c b b a c e], we express the elements of the full
vector as their corresponding integer position in the array of
unique elements. The corresponding position vector will be
P = [1 2 3 1 3 2 2 1 3 4].

"""
function whichpoints_visit_whichbins(A, U)
    a = [A[:,i] for i in 1:size(A, 2)]
    u = [U[:,i] for i in 1:size(U, 2)]
    inds = Int[]
    npts = length(a)
    npts_unique = length(u)
    for j = 1:npts
        for i = 1:npts_unique
            if a[j] == u[i]
                push!(inds, i)
            end
        end
    end
    inds
end

"""
    get_binvisits(visited_bin_labels)

Organise bin labels by identifying bins that are visited multiple times,
and finding which points fall in which bins.

`visited_bin_labels` is an array where each column cᵢ (corresponding to
the point pᵢ) represents a tuple of integers identifying the position
of that bin, relative to some minimum value, given a stepsize `ϵ`.
See `assign_bin_labels` or `encode` for more information.
"""
function get_binvisits(visited_bin_labels::Array{Int, 2})
    slices = groupslices(visited_bin_labels, 2)
    first_visited_by = firstinds(slices)
    visitors = groupinds(slices)
    visits_whichbin = whichpoints_visit_whichbins(visited_bin_labels, unique(visited_bin_labels, dims = 2))

    BinVisits(first_visited_by, visitors, visits_whichbin)
end

get_binvisits(points::Vector{Vector{T}}) where T = get_binvisits(hcat(points...,))

export
BinVisits,
get_binvisits


text:
using Revise
using JLD2
using ViewerGL
using LarSurf

# created by surface_extraction_parallel_ircad01.jl

# mask_labels = ["portalvein", "liver"]
mask_labels = ["liver", "portalvein"]
views = []
view = ViewerGL.GLAxis(ViewerGL.Point3d(-1,-1,-1),ViewerGL.Point3d(1,1,1))
push!(views, view)

for mask_label in mask_labels
	@info "mask_label: $mask_label"
	@load "ircad_$(mask_label)_sm.jld2" Vs FVtri

    view = ViewerGL.GLGrid(Vs,FVtri,ViewerGL.Point4d(1,1,1,0.1))
	push!(views, view)
end

ViewerGL.VIEW(views)