text:
module Variance
using ValidatedNumerics;
using Distributions;
using Dierckx;

abstract BaseEstimator

function estimate(values, method::BaseEstimator)
    error("Unsupported method $method")
end

immutable byConfidenceProbability <: BaseEstimator
    confidence_probability :: Real
end

function estimate(values, method :: byConfidenceProbability)
    average = mean(values)
    count = length(values)
    gamma_down, gamma_up = getGammas(method.confidence_probability)
    quantile_down = getChiSquareQuantile(gamma_down, count)
    quantile_up = getChiSquareQuantile(gamma_up, count)

    reduced = mapreduce((x) -> (x - average)^2, +, values)

    @interval(reduced/quantile_down, reduced/quantile_up)
end

immutable byMeanAbsoluteDeviation <: BaseEstimator
    confidence_probability :: Real
end

function estimate(values, method :: byMeanAbsoluteDeviation)
    average = mean(values)
    count = length(values)
    quantile_down, quantile_up = getMeanAbsDeviationQuantiles(method.confidence_probability, count)

    mean_abs_deviation = mapreduce((x) -> abs(x - average), +, values)/count

    @interval(mean_abs_deviation/quantile_down, mean_abs_deviation/quantile_up)^2
end

immutable byPointVariance <: BaseEstimator
    confidence_probability :: Real
    variance :: Real
end

function estimate(values, method::byPointVariance)
    count = length(values)
    gamma_down, gamma_up = getGammas(method.confidence_probability)
    quantile_down = getChiSquareQuantile(gamma_down, count)
    quantile_up = getChiSquareQuantile(gamma_up, count)

    fixed_variance = method.variance * (1 + 0.254 / (count - 1))^2
    term = (count - 1) * fixed_variance

    @interval(term / quantile_down, term / quantile_up)
end

function getGammas(confidence_probability)
    (1 + confidence_probability) * 0.5, (1 - confidence_probability) * 0.5
end

function getChiSquareQuantile(value, count)
    d = Chisq(count - 1)
    quantile(d, value)
end

function getMeanAbsDeviationQuantiles(confidence_probability, count)
    if confidence_probability == 0.90
        const knots = [
        2 1.386 0.044;
        3 1.276 0.166;
        4 1.224 0.254;
        5 1.187 0.305;
        6 1.158 0.360;
        7 1.135 0.394;
        8 1.116 0.422;
        9 1.100 0.445;
        10 1.086 0.464;
        ]
    elseif confidence_probability == 0.95
        const knots = [
        2 1.585 0.022;
        3 1.417 0.116;
        4 1.344 0.199;
        5 1.292 0.260;
        6 1.253 0.306;
        7 1.222 0.342;
        8 1.196 0.372;
        9 1.175 0.396;
        10 1.156 0.417;
        ]
    elseif confidence_probability == 0.99
        const knots = [
        2 1.985 0.004;
        3 1.703 0.073;
        4 1.590 0.145;
        5 1.507 0.203;
        6 1.445 0.250;
        7 1.397 0.287;
        8 1.358 0.318;
        9 1.326 0.344;
        10 1.299 0.366;
        ]
    else
        error("Value=$confidence_probability is not allowed")
    end
    x_knots = knots[:, 1]
    if count > last(x_knots)
        error("Error $count too high")
    end

    y_down_knots = knots[:, 2]
    down_itp = Spline1D(x_knots, y_down_knots, k=1, bc="extrapolate")

    y_up_knots = knots[:, 3]
    up_itp = Spline1D(x_knots, y_up_knots, k=1, bc="extrapolate")
    evaluate(down_itp, count), evaluate(up_itp, count)
end

end


text:
module ParaxialBeamPropagation

using FFTW, Interpolations,
ClassicalOrthogonalPolynomials,SpecialFunctions, Plots,
BenchmarkTools, PlutoUI,Printf,LaTeXStrings

export get_fourier_transform,display_fourier_transform,
propagate_beam,display_beam,animate_beam,HG,LG,airy_beam,vertical_obstacle,horizontal_obstacle,
circular_obstacle,rectangular_obstacle

include("Propagation.jl")
include("InitialProfiles.jl")

end

ext: jl


text:
export PendulumEnv

struct PendulumEnvParams{T}
    max_speed::T
    max_torque::T
    g::T
    m::T
    l::T
    dt::T
    max_steps::Int
end

mutable struct PendulumEnv{A,T,R<:AbstractRNG} <: AbstractEnv
    params::PendulumEnvParams{T}
    action_space::A
    observation_space::Space{Vector{ClosedInterval{T}}}
    state::Vector{T}
    done::Bool
    t::Int
    rng::R
    reward::T
    n_actions::Int
end

"""
    PendulumEnv(;kwargs...)

# Keyword arguments

- `T = Float64`
- `max_speed = T(8)`
- `max_torque = T(2)`
- `g = T(10)`
- `m = T(1)`
- `l = T(1)`
- `dt = T(0.05)`
- `max_steps = 200`
- `continuous::Bool = true`
- `n_actions::Int = 3`
- `rng = Random.GLOBAL_RNG`
"""
function PendulumEnv(;
    T = Float64,
    max_speed = T(8),
    max_torque = T(2),
    g = T(10),
    m = T(1),
    l = T(1),
    dt = T(0.05),
    max_steps = 200,
    continuous::Bool = true,
    n_actions::Int = 3,
    rng = Random.GLOBAL_RNG,
)
    high = T.([1, 1, max_speed])
    action_space = continuous ? -2.0..2.0 : Base.OneTo(n_actions)
    env = PendulumEnv(
        PendulumEnvParams(max_speed, max_torque, g, m, l, dt, max_steps),
        action_space,
        Space(ClosedInterval{T}.(-high, high)),
        zeros(T, 2),
        false,
        0,
        rng,
        zero(T),
        n_actions,
    )
    reset!(env)
    env
end

Random.seed!(env::PendulumEnv, seed) = Random.seed!(env.rng, seed)

pendulum_observation(s) = [cos(s[1]), sin(s[1]), s[2]]
angle_normalize(x) = Base.mod((x + Base.π), (2 * Base.π)) - Base.π

RLBase.action_space(env::PendulumEnv) = env.action_space
RLBase.state_space(env::PendulumEnv) = env.observation_space
RLBase.reward(env::PendulumEnv) = env.reward
RLBase.is_terminated(env::PendulumEnv) = env.done
RLBase.state(env::PendulumEnv) = pendulum_observation(env.state)

function RLBase.reset!(env::PendulumEnv{A,T}) where {A,T}
    env.state[1] = 2 * π * (rand(env.rng, T) .- 1)
    env.state[2] = 2 * (rand(env.rng, T) .- 1)
    env.t = 0
    env.done = false
    env.reward = zero(T)
    nothing
end

function (env::PendulumEnv{<:ClosedInterval})(a::AbstractFloat)
    @assert a in env.action_space
    _step!(env, a)
end

function (env::PendulumEnv{<:Base.OneTo})(a::Int)
    @assert a in env.action_space
    float_a = (4 / (env.n_actions - 1)) * (a - (env.n_actions - 1) / 2 - 1)
    _step!(env, float_a)
end

function _step!(env::PendulumEnv, a)
    env.t += 1
    th, thdot = env.state
    a = clamp(a, -env.params.max_torque, env.params.max_torque)
    costs = angle_normalize(th)^2 + 0.1 * thdot^2 + 0.001 * a^2
    newthdot =
        thdot +
        (
            -3 * env.params.g / (2 * env.params.l) * sin(th + pi) +
            3 * a / (env.params.m * env.params.l^2)
        ) * env.params.dt
    th += newthdot * env.params.dt
    newthdot = clamp(newthdot, -env.params.max_speed, env.params.max_speed)
    env.state[1] = th
    env.state[2] = newthdot
    env.done = env.t >= env.params.max_steps
    env.reward = -costs
    nothing
end


text:
# ---
# title: 1552. Magnetic Force Between Two Balls
# id: problem1552
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Binary Search
# link: <https://leetcode.com/problems/magnetic-force-between-two-balls/description/>
# hidden: true
# ---
#
# In universe Earth C-137, Rick discovered a special form of magnetic force
# between two balls if they are put in his new invented basket. Rick has `n`
# empty baskets, the `ith` basket is at `position[i]`, Morty has `m` balls and
# needs to distribute the balls into the baskets such that the **minimum
# magnetic force** between any two balls is **maximum**.
#
# Rick stated that magnetic force between two different balls at positions `x`
# and `y` is `|x - y|`.
#
# Given the integer array `position` and the integer `m`. Return _the required
# force_.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/08/11/q3v1.jpg)
#
#
#
#     Input: position = [1,2,3,4,7], m = 3
#     Output: 3
#     Explanation: Distributing the 3 balls into baskets 1, 4 and 7 will make the magnetic force between ball pairs [3, 3, 6]. The minimum magnetic force is 3. We cannot achieve a larger minimum magnetic force than 3.
#
#
# **Example 2:**
#
#
#
#     Input: position = [5,4,3,2,1,1000000000], m = 2
#     Output: 999999999
#     Explanation: We can use baskets 1 and 1000000000.
#
#
#
#
# **Constraints:**
#
#   * `n == position.length`
#   * `2 <= n <= 10^5`
#   * `1 <= position[i] <= 10^9`
#   * All integers in `position` are **distinct**.
#   * `2 <= m <= position.length`
#
#
## @lc code=start
using LeetCode

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


text:
"""
Arrays with fixed size and circular indexing.
"""
module CircularArrays

export CircularArray, CircularVector

"""
    CircularArray{T, N, A} <: AbstractArray{T, N}

`N`-dimensional array backed by an `AbstractArray{T, N}` of type `A` with fixed size and circular indexing.

    array[index] == array[mod1(index, size)]
"""
struct CircularArray{T, N, A} <: AbstractArray{T, N}
    data::A
    CircularArray{T,N}(data::AbstractArray{T,N}) where {T,N} = new{T,N,typeof(data)}(data)
end

"""
    CircularVector{T,A} <: AbstractVector{T}

One-dimensional array backed by an `AbstractArray{T, 1}` of type `A` with fixed size and circular indexing.
Alias for [`CircularArray{T,1,A}`](@ref).

    array[index] == array[mod1(index, size)]
"""
const CircularVector{T} = CircularArray{T, 1}

@inline clamp_bounds(arr::CircularArray, I::Tuple{Vararg{Int}})::AbstractArray{Int, 1} = map(Base.splat(mod), zip(I, axes(arr.data)))

CircularArray(data::AbstractArray{T,N}) where {T,N} = CircularArray{T,N}(data)
CircularArray{T}(data::AbstractArray{T,N}) where {T,N} = CircularArray{T,N}(data)
CircularArray(def::T, size) where T = CircularArray(fill(def, size))

@inline Base.getindex(arr::CircularArray, i::Int) = @inbounds getindex(arr.data, mod(i, Base.axes1(arr.data)))
@inline Base.setindex!(arr::CircularArray, v, i::Int) = @inbounds setindex!(arr.data, v, mod(i, Base.axes1(arr.data)))
@inline Base.getindex(arr::CircularArray, I::Vararg{Int}) = @inbounds getindex(arr.data, clamp_bounds(arr, I)...)
@inline Base.setindex!(arr::CircularArray, v, I::Vararg{Int}) = @inbounds setindex!(arr.data, v, clamp_bounds(arr, I)...)
@inline Base.size(arr::CircularArray) = size(arr.data)
@inline Base.axes(arr::CircularArray) = axes(arr.data)

@inline Base.checkbounds(::CircularArray, _...) = nothing

@inline _similar(arr::CircularArray, ::Type{T}, dims) where T = CircularArray(similar(arr.data,T,dims))
@inline Base.similar(arr::CircularArray, ::Type{T}, dims::Tuple{Base.DimOrInd, Vararg{Base.DimOrInd}}) where T = _similar(arr,T,dims)
# Ambiguity resolution with Base
@inline Base.similar(arr::CircularArray, ::Type{T}, dims::Tuple{Int64,Vararg{Int64}}) where T = _similar(arr,T,dims)

CircularVector(data::AbstractArray{T, 1}) where T = CircularVector{T}(data)
CircularVector(def::T, size::Int) where T = CircularVector{T}(fill(def, size))

Base.IndexStyle(::Type{<:CircularVector}) = IndexLinear()

end


text:
"""
TODO: Add description of distance module.
"""
module distance
using utility
export dist, Distance

"""
The base type used by all distance metrics.
"""
abstract Distance

# Load each of the distance metrics.
include("distance/euclidean.jl")
include("distance/hamming.jl")
include("distance/levenshtein.jl")
end


text:
module GWKT
# export PropagatorKT, Weight, addChild
include("diagtree.jl")
export DiagTree

include("manual/manual.jl")
export Manual

end


text:
export init_mapf_problem

function init_mapf_problem(graph,start_vtxs,goal_vtxs,
        cost_model=SumOfTravelTime(),
        heuristic = PerfectHeuristic(graph,start_vtxs,goal_vtxs),
        args...)
    env = CBSEnv.LowLevelEnv(
        graph = graph,
        cost_model = cost_model,
        heuristic = heuristic
    )
    starts = map(s->CBSEnv.State(vtx=s,t=0),start_vtxs)
    goals = map(s->CBSEnv.State(vtx=s,t=0),goal_vtxs)
    MAPF(env,starts,goals)
end

export
    init_mapf_1,
    init_mapf_2,
    init_mapf_3,
    init_mapf_4,
    init_mapf_5,
    init_mapf_6,
    init_mapf_7,
    init_mapf_8

function init_mapf_1(args...)
    vtx_grid = initialize_dense_vtx_grid(4,4)
    #  1   2   3   4
    #  5   6   7   8
    #  9  10  11  12
    # 13  14  15  16
    starts = [1,4]
    goals = [13,16]
    graph = construct_factory_env_from_vtx_grid(vtx_grid)
    init_mapf_problem(graph,starts,goals,args...)
end

"""
    switch places
"""
function init_mapf_2(args...)
    vtx_grid = initialize_dense_vtx_grid(4,4)
    #  1   2   3   4
    #  5   6   7   8
    #  9  10  11  12
    # 13  14  15  16
    starts = [1,4]
    goals = [4,1]
    graph = construct_factory_env_from_vtx_grid(vtx_grid)
    init_mapf_problem(graph,starts,goals,args...)
end

"""
    congested
"""
function init_mapf_3(args...)
    vtx_grid = initialize_dense_vtx_grid(4,4)
    #  1   2   3   4
    #  5   6   7   8
    #  9  10  11  12
    # 13  14  15  16
    starts = [1,2,3,4,5,6,7,8]
    goals = [13,14,15,16,9,10,11,12]
    graph = construct_factory_env_from_vtx_grid(vtx_grid)
    init_mapf_problem(graph,starts,goals,args...)
end

"""
    almost switch corners. With the fat path heuristic, the paths should be:
    - Robot 1: [1,5,9,13,14,15]
    - Robot 2: [16,12,8,4,3,2]
"""
function init_mapf_4(args...)
    vtx_grid = initialize_dense_vtx_grid(4,4)
    #  1   2   3   4
    #  5   6   7   8
    #  9  10  11  12
    # 13  14  15  16
    starts = [1,16]
    goals = [15,2]
    graph = construct_factory_env_from_vtx_grid(vtx_grid)
    init_mapf_problem(graph,starts,goals,args...)
end

"""
    init_mapf_5

PIBT demo from paper. CBS fails to solve it even with 80000 iterations.
"""
function init_mapf_5(args...)
    vtx_grid = initialize_dense_vtx_grid(2,4)
    #  1   2   3   4
    #  5   6   7   8
    starts = [2,5,6,3,4,7,8]
    goals = [4,8,5,7,8,4,3]
    graph = construct_factory_env_from_vtx_grid(vtx_grid)
    init_mapf_problem(graph,starts,goals,args...)
end

"""
    PIBT demo from paper
"""
function init_mapf_6(args...)
    vtx_grid = initialize_dense_vtx_grid(3,3)
    # 1  2  3
    # 4  5  6
    # 7  8  9
    starts = [2,4]
    goals = [8,6]
    graph = construct_factory_env_from_vtx_grid(vtx_grid)
    init_mapf_problem(graph,starts,goals,args...)
end

"""
    init_mapf_7

Robot 3 has to get past robots 1 and 2.
"""
function init_mapf_7(args...)
    vtx_grid = initialize_dense_vtx_grid(2,3)
    # 1  2  3
    # 4  5  6
    starts = [2,4,1]
    goals = [2,4,6]
    graph = construct_factory_env_from_vtx_grid(vtx_grid)
    init_mapf_problem(graph,starts,goals,args...)
end

"""
    init_mapf_8

Robots 1 and 2 have to get past robots 3 and 4. Requires about 660 iterations of
regular CBS.
"""
function init_mapf_8(args...)
    vtx_grid = initialize_dense_vtx_grid(2,3)
    # 1  2  3
    # 4  5  6
    starts = [1,4,2,5]
    goals = [6,3,2,5]
    graph = construct_factory_env_from_vtx_grid(vtx_grid)
    init_mapf_problem(graph,starts,goals,args...)
end

export mapf_test_problems

mapf_test_problems() = [
    init_mapf_1,
    init_mapf_2,
    init_mapf_3,
    init_mapf_4,
    init_mapf_5,
]


text:
#
# Address, Person types for deep match test
#

struct Address
    street::AbstractString
    city::AbstractString
    zip::AbstractString
end

struct Person
    firstname::AbstractString
    lastname::AbstractString
    address::Address
end

#
# Untyped lambda calculus definitions
#

abstract type Term end

struct Var <: Term
    name::AbstractString
end

struct Fun <: Term
    arg::AbstractString
    body::Term
end

struct App <: Term
    f::Term
    v::Term
end


text:
"""Abstract cell type"""
abstract type AbstractCell end

"""Dimension of the cell"""
dim(c::AbstractCell) = throw(MethodError(dim,(typeof(c),)))

"""Cell comparison"""
==(a::AbstractCell, b::AbstractCell) = throw(MethodError(==,(typeof(a), typeof(b))))

"""
    faces(cell)

Get a collection of `cell` faces.
"""
faces(c::AbstractCell) = throw(MethodError(faces,(typeof(c),)))

"""
    hash(cell)

Get a unique identifier of the `cell`
"""
hash(c::AbstractCell) = throw(MethodError(hash,(typeof(c),)))

"""
    union(u, v)

Create a new cell from a combination of vertices of cells `u` and `v`.
"""
union(u::C, v::C) where{C<:AbstractCell} = throw(MethodError(union,(C,C)))

"""Abstract simplex type"""
abstract type AbstractSimplex <: AbstractCell end

"""
    values(cell)

Get an array of `cell` values.
"""
values(c::AbstractSimplex) = throw(MethodError(values,(typeof(c),))) #TODO: consider to move to AbstractSimplex

"""
    vertecies(cell)

Get a collection of `cell` vertecies"""
vertecies(c::AbstractSimplex) = throw(MethodError(vertecies,(typeof(c),)))


text:
using Documenter, BayesianLinearRegressors

makedocs(;
    modules=[BayesianLinearRegressors],
    format=Documenter.HTML(),
    pages=["Home" => "index.md"],
    repo="https://github.com/JuliaGaussianProcesses/BayesianLinearRegressors.jl/blob/{commit}{path}#L{line}",
    sitename="BayesianLinearRegressors.jl",
    authors="Will Tebbutt <wt0881@my.bristol.ac.uk>",
    assets=String[],
)

deploydocs(; repo="github.com/JuliaGaussianProcesses/BayesianLinearRegressors.jl")


text:
module ChangePoint

include("../src/Jaynes.jl")
using .Jaynes
using Gen

# ------------ Change point model ------------ #

jfunc = @jaynes function change_point_model(N::Int)
    n ~ DiscreteUniform(N)
    λ₁ ~ Gamma(2, 2)
    λ₂ ~ Gamma(3, 2)
    x = [i >= n ? {:x => i} ~ Poisson(λ₂) : {:x => i} ~ Poisson(λ₁) for i in 1 : N]
    x
end

# Sample trace.
tr = simulate(jfunc, (50, ))
display(tr)

# ------------ Gibbs ------------ #

@kern function gibbs_kernel(trace)
    trace ~ mh(trace, select(:n))
    trace ~ mh(trace, select(:λ₁))
    trace ~ mh(trace, select(:λ₂))
end

# ------------ Inference ------------ #

infer = () -> begin
    obs = constrain([i <= 30 ? (:x => i, ) => rand(Poisson(3)) : (:x => i, ) => rand(Poisson(7)) for i in 1 : 50])

    # Sample trace.
    tr, w = generate(jfunc, (50, ), obs)

    trs = Trace[]
    for i in 1 : 5000
        tr, _ = gibbs_kernel(tr)
        i % 50 == 0 && push!(trs, tr)
    end

    # Posterior est.
    λ₁_est, λ₂_est, n_est = zip(map(trs) do tr
                                    (tr[:λ₁], tr[:λ₂], tr[:n])
                                end...)
    len = length(λ₁_est)
    println("est: $((sum(λ₁_est) / len,
                     sum(λ₂_est) / len,
                     sum(n_est) / len))")

end
infer()
@time infer()

end # module


text:
# Add dependencies
using Pkg
Pkg.add(Pkg.PackageSpec(;name="Git", version="1.2.1"))
Pkg.add(Pkg.PackageSpec(;name="TOML", version="1.0.0"))
using TOML
using Git

# Extract the version number to be updated
VERSION = ""
if length(ARGS) > 0
    VERSION = ARGS[1]
end

GITHUB_REPOSITORY = ENV["GITHUB_REPOSITORY"]
TOKEN = ""
URL = ""
if haskey(ENV, "GITHUB_TOKEN")
    TOKEN = ENV["GITHUB_TOKEN"]
    URL = "https://x-access-token:$(TOKEN)@github.com/$(GITHUB_REPOSITORY).git"
end

# Read the Project.toml file in the package
fname = "Project.toml";
dict_project = TOML.parsefile(fname)

# If no version is given, use the version number in the Project.toml
if isempty(VERSION)
    VERSION = dict_project["version"]
else
    # Santity
    VERSION = split(VERSION, "-")[1]
    VERSION = replace(VERSION, "v" => "")
    VERSION = replace(VERSION, "V" => "")

    # Update the version number in the Project.toml
    if dict_project["version"] != VERSION
        dict_project["version"] = VERSION
        open(fname, "w") do io
            TOML.print(io, dict_project)
        end

        # Commit the new Project.toml
        run(`$(git()) add Project.toml`)
        run(`$(git()) commit -m "Update version to v$(VERSION)"`)
        run(`$(git()) restore .`)
        if isempty(URL)
            run(`$(git()) push`)
        else
            run(`$(git()) push $(URL)`)
        end
    end
end

# Set output
run(`echo "::set-output name=version::$(VERSION)"`)


text:
data=readdlm("MIST_iso_1470351180.iso")

using CGS
rho_sun = MSUN/(4pi/3.*RSUN^3)
# Convert to density:
mass = vec(data[:,4])
density = mass./10.^(3.*data[:,12])*rho_sun
teff = vec(10.^data[:,11])

using PyPlot

fig,axes = subplots(2,1)
ax = axes[1,1]
ax[:plot](density,teff,"r-",alpha=0.6,label="t = 1.5 Gyr, [Fe/H]=0.0, MIST",linewidth=3)

rhoA=0.001167
srhoA=0.000052
tA = 4750.
stA = 250.
ax[:plot](rhoA+srhoA*[-1,1],tA*[1,1],"g-",linewidth=3,alpha=0.6)
ax[:plot](rhoA*[1,1],tA+stA*[-1,1],"g-",linewidth=3,alpha=0.6,label="Star A")

rhoB=0.0824
srhoB=0.0034
tB = 6500.
stB = 500.
ax[:semilogx](rhoB+srhoB*[-1,1],tB*[1,1],"b-",linewidth=3,alpha=0.6)
ax[:semilogx](rhoB*[1,1],tB+stB*[-1,1],"b-",linewidth=3,alpha=0.6,label="Star B")
ax[:axis]([1e-4,1e0,4e3,8e3])
#ax[:set_xlabel]("Density [g/cc]")
ax[:set_ylabel]("Teff [K]")
ax[:legend](loc = "upper left")

ax = axes[2,1]

ax[:semilogx](density,mass,"r-",alpha=0.6,label="t = 1.5 Gyr, [Fe/H]=0.0, MIST",linewidth=3)
mA = 1.9
smA = 0.27
ax[:semilogx](rhoA+srhoA*[-1,1],mA*[1,1],"g-",linewidth=3,alpha=0.6)
ax[:semilogx](rhoA*[1,1],mA+smA*[-1,1],"g-",linewidth=3,alpha=0.6,label="Star A")
ax[:semilogx](rhoB+srhoB*[-1,1,1,-1,-1],[0,0,10,10,0],"b-",linewidth=3,alpha=0.6,label="Star B")
ax[:axis]([1e-4,1e0,1.7,2.0])
ax[:set_xlabel]("Density [g/cc]")
ax[:set_ylabel]("Mass [Msun]")

PyPlot.savefig("isochrone.png",bbox_inches="tight")


text:
using MathMLRepr
using Mxma

#m = MExpr(:(f(x)=x^2+1))
#m = MExpr(:(a=1; b=x^2+2x+1))
#r = mcall1(m)
#println(r)

# Julia expressions => mstr
#mcall(:(b=x^2+2x+1))
#m1 = MExpr("sprint(b)\$")
#m1 = MExpr("sprint(b);")
#m1 = MExpr("a: 1\$ c: 2\$ print(b)\$")
#m1 = MExpr("sprint(b)")
#m1 = MExpr("a: matrix([1,2],[3,4])")
#mcall(MExpr(:(A = [1 2; 3 4])))
#mcall(MExpr(:(A = [1/2, 2/3, 3/4, 4/5])))
#mcall(MExpr(:(A = 34//87)))
#m1 = MExpr("sprint(A)\$")
#r = mcall(m1)
#println(r)

# Julia types => mstr
#A = rand(1:9, 2, 4)
#A = [1//2, 3//4, 5//6]
#A = rand(1:9, 3, 3); A = convert.(Rational, A)
#A = [:π :(x^2); :(x+1) :(2x)]

#s = @mxset(A)
# Some examples
#s = @mxc f = x^2+2x+1
#println(s)
#r = @mxc factor(f)
#println(r)
#r = @mxc expand((x+1)*(x+2))
#println(r)
#A = rand(1:9, 2, 4)
#r = @mxset A
#println(r)
#r = @mxc 2*[1,2,3,4]
#println(r)

# latex output
#r = @mxc f = x^2 + 3x + 2
#r = @mxc factor(f)
#io = IOBuffer()
#show(io, MIME("text/mathml+xml"), r)
#s = String(take!(io))
#println(s)

# Issue
#@mxc begin a=1;b=2 end
#r = @mx b=[1/2 2/3; 3/4 4/5]
#r = @mx b=[1 2; 3 4]
#r = @mx "cauchy_matrix([a,b,c])"
#io = IOBuffer()
#show(io, MIME("text/mathml+xml"), r)
#s = String(take!(io))
#println(s)
#println(typeof(r))
#println(length(Base.Multimedia.displays))
#display(r);println()
#display([r, r])
#display("text/plain", r); println()

# convert MExpr to LATEX
s = @mx fex=t^(1/2)*e^(-a*t/4)
t = convert(LATEX, s)
println(s)
println(t)


text:
"""
    mutable struct SalientPoleExponential <: Machine
        base_machine::SalientPoleMachine
        saturation_coeffs::Tuple{Float64, Float64}

3-states salient-pole synchronous machine with exponential saturation:
IEEE Std 1110 §5.3.2 (Model 2.1). GENSAE in PSSE and PSLF.

# Arguments:
- `base_machine::SalientPoleMachine`: Salient Pole Machine model.
- `saturation_coeffs::Tuple{Float64, Float64}``: Saturation coefficients for exponential model.
"""
mutable struct SalientPoleExponential <: Machine
    base_machine::SalientPoleMachine
    saturation_coeffs::Tuple{Float64, Float64}
end
IS.@forward((SalientPoleExponential, :base_machine), SalientPoleMachine)

function SalientPoleExponential(
    R::Float64,
    Td0_p::Float64,
    Td0_pp::Float64,
    Tq0_pp::Float64,
    Xd::Float64,
    Xq::Float64,
    Xd_p::Float64,
    Xd_pp::Float64,
    Xl::Float64,
    Se::Tuple{Float64, Float64},
)
    saturation_coeffs = get_exponential_saturation(Se)
    return SalientPoleExponential(
        SalientPoleMachine(R, Td0_p, Td0_pp, Tq0_pp, Xd, Xq, Xd_p, Xd_pp, Xl, Se),
        saturation_coeffs,
    )
end

function SalientPoleExponential(; R, Td0_p, Td0_pp, Tq0_pp, Xd, Xq, Xd_p, Xd_pp, Xl, Se)
    return SalientPoleExponential(R, Td0_p, Td0_pp, Tq0_pp, Xd, Xq, Xd_p, Xd_pp, Xl, Se)
end

function SalientPoleExponential(::Nothing)
    return SalientPoleExponential(;
        R = 0.0,
        Td0_p = 0.0,
        Td0_pp = 0.0,
        Tq0_pp = 0.0,
        Xd = 0.0,
        Xq = 0.0,
        Xd_p = 0.0,
        Xd_pp = 0.0,
        Xl = 0.0,
        Se = (0.0, 0.0),
    )
end

get_base_machine(value::SalientPoleExponential) = value.base_machine
get_saturation_coeffs(value::SalientPoleExponential) = value.saturation_coeffs

set_base_machine!(value::SalientPoleExponential, val::SalientPoleMachine) =
    value.base_machine = val
set_saturation_coeffs!(value::SalientPoleExponential, val::Tuple{Float64, Float64}) =
    value.saturation_coeffs = val

function IS.deserialize_struct(::Type{SalientPoleExponential}, data::Dict)
    vals = IS.deserialize_to_dict(SalientPoleExponential, data)
    return SalientPoleExponential(vals[:base_machine], vals[:saturation_coeffs])
end


text:
using Test
using DataDeps: fetch_default, fetch_base, fetch_http

ENV["DATADEPS_ALWAYS_ACCEPT"] = true

@testset "easy https url" begin
    url = "https://www.angio.net/pi/digits/10000.txt"
    # This is easy because the filename is in the URL
    # So it works with both `fetch_base` and `fetch_http`
    # HTTP.jl has tests for much more difficult cases, and fetch_http supports those
    @testset "$fetch_func" for fetch_func in (fetch_default, fetch_base, fetch_http)
        mktempdir() do localdir
            localpath = fetch_func(url, localdir)
            @test isfile(localpath)
            @test localpath == joinpath(localdir, "10000.txt")
            @test stat(localpath).size == 10_001
        end
    end
end


text:
using Documenter
using PiGPIO
using Literate

Literate.markdown(joinpath(@__DIR__, "..", "examples", "01_blink.jl"), joinpath(@__DIR__, "src", "examples"))
Literate.markdown(joinpath(@__DIR__, "..", "examples", "02_blink_twice.jl"),joinpath(@__DIR__, "src", "examples"))
Literate.markdown(joinpath(@__DIR__, "..", "examples", "03_rgb.jl"), joinpath(@__DIR__, "src", "examples"))

makedocs(
    sitename = "PiGPIO",
    format = Documenter.HTML(),
    modules = [PiGPIO],
    pages = [
        "index.md",
        "Tutorial" => "tutorial.md"
        "API Docs" => "api.md",
        "Examples" => [
            "Blink Once" => "examples/01_blink.md",
            "Blink Twice" => "examples/02_blink_twice.md",
            "Red-Green-Blue" => "examples/03_rgb.md"
        ]
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "<repository url>"
)=#


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/AbaqusReader.jl/blob/master/LICENSE

import Base.parse

# Define element type and number of nodes in element
element_has_nodes(::Type{Val{:C3D4}}) = 4
element_has_type( ::Type{Val{:C3D4}}) = :Tet4

element_has_nodes(::Type{Val{:C3D4H}}) = 4
element_has_type( ::Type{Val{:C3D4H}}) = :Tet4

element_has_nodes(::Type{Val{:C3D8}}) = 8
element_has_type( ::Type{Val{:C3D8}}) = :Hex8

element_has_nodes(::Type{Val{:C3D10}}) = 10
element_has_type(::Type{Val{:C3D10}}) = :Tet10

element_has_nodes(::Type{Val{:C3D10H}}) = 10
element_has_type(::Type{Val{:C3D10H}}) = :Tet10

element_has_nodes(::Type{Val{:C3D20}}) = 20
element_has_type(::Type{Val{:C3D20}}) = :Hex20

element_has_nodes(::Type{Val{:C3D20E}}) = 20

element_has_nodes(::Type{Val{:S3}}) = 3
element_has_type( ::Type{Val{:S3}}) = :Tri3

element_has_nodes(::Type{Val{:STRI65}}) = 6
element_has_type(::Type{Val{:STRI65}}) = :Tri6

element_has_nodes(::Type{Val{:CPS4}}) = 4
element_has_type(::Type{Val{:CPS4}}) = :Quad4

element_has_nodes(::Type{Val{:T2D2}}) = 2
element_has_type(::Type{Val{:T2D2}}) = :Seg2

element_has_nodes(::Type{Val{:T3D2}}) = 2
element_has_type(::Type{Val{:T3D2}}) = :Seg2

element_has_nodes(::Type{Val{:B33}}) = 2
element_has_type(::Type{Val{:B33}}) = :Seg2

"""Checks for a comment or empty line

Function return true, if line starts with comment character "**"
or has length of 0
"""
function empty_or_comment_line(line::T) where T<:AbstractString
    startswith(line, "**") || (length(line) == 0)
end

"""Match words from both sides of '=' character
"""
function matchset(definition)
    regexp = r"([\w\_\-]+[ ]*=[ ]*[\w\_\-]+)"
    collect(m.match for m = eachmatch(regexp, definition))
end

"""Parse string to get set type and name
"""
function parse_definition(definition)
    set_defs = Dict()
    set_definition = matchset(definition)
    set_definition == nothing && return nothing
    for x in set_definition
        name, vals = map(strip, split(x, "="))
        set_defs[lowercase(name)] = vals
    end
    set_defs
end

"""Parse all the numbers from string
"""
function parse_numbers(line, type_::Type{T})::Vector{T} where T
    regexp = r"[0-9]+"
    matches = collect((m.match for m = eachmatch(regexp, line)))
    map(x-> Base.parse(type_, x), matches)
end

"""Add set to model, if set exists
"""
function add_set!(model, definition, model_key, abaqus_key, ids)
    has_set_def = parse_definition(definition)
    if haskey(has_set_def, "elset")
        set_name = has_set_def[abaqus_key]
        @info("Adding $abaqus_key: $set_name")
        model[model_key][set_name] = ids
    end
end

"""Parse nodes from the lines
"""
function parse_section(model, lines, ::Symbol, idx_start, idx_end, ::Type{Val{:NODE}})
    nnodes = 0
    ids = Int[]
    definition = lines[idx_start]
    for line in lines[idx_start + 1: idx_end]
        if !(empty_or_comment_line(line))
            m = collect((m.match for m = eachmatch(r"[-0-9.eE+]+", line)))
            node_id = parse(Int, m[1])
            coords = parse.(Float64, m[2:end])
            model["nodes"][node_id] = coords
            push!(ids, node_id)
            nnodes += 1
        end
    end
    @info("$nnodes nodes found")
    add_set!(model, definition, "node_sets", "nset", ids)
end

"""Custon regex to find match from string. Index used if there are multiple matches
"""
function regex_match(regex_str, line, idx)
    return match(regex_str, line).captures[idx]
end

"""Custom list iterator

Simple iterator for comsuming element list. Depending
on the used element, connectivity nodes might be listed
in multiple lines, which is why iterator is used to handle
this problem.
"""
function consumeList(arr, start, stop)
    idx = start - 1
    function _it()
        idx += 1
        if idx > stop
            return nothing
        end
        arr[idx]
    end
    _it
end

"""Parse elements from input lines

Reads element ids and their connectivity nodes from input lines.
If elset definition exists, also adds the set to model.
"""
function parse_section(model, lines, ::Symbol, idx_start, idx_end, ::Type{Val{:ELEMENT}})
    ids = Int[]
    definition = lines[idx_start]
    regexp = r"TYPE=([\w\-\_]+)"i
    m = match(regexp, definition)
    m == nothing && error("Could not match regexp $regexp to line $definition")
    element_type = uppercase(m[1])
    eltype_sym = Symbol(element_type)
    eltype_nodes = element_has_nodes(Val{eltype_sym})
    element_type = element_has_type(Val{eltype_sym})
    @info("Parsing elements. Type: $(m[1]). Topology: $(element_type)")
    list_iterator = consumeList(lines, idx_start+1, idx_end)
    line = list_iterator()
    while line != nothing
        numbers = parse_numbers(line, Int)
        if !(empty_or_comment_line(line))
            id = numbers[1]
            push!(ids, id)
            connectivity = numbers[2:end]
            while length(connectivity) != eltype_nodes
                @assert length(connectivity) < eltype_nodes
                line = list_iterator()
                numbers = parse_numbers(line, Int)
                push!(connectivity, numbers...)
            end
            model["elements"][id] = connectivity
            model["element_types"][id] = element_type
        end
        line = list_iterator()
    end
    add_set!(model, definition, "element_sets", "elset", ids)
end

"""Parse node and elementset from input lines
"""
function parse_section(model, lines, key, idx_start, idx_end, ::Union{Type{Val{:NSET}},
        Type{Val{:ELSET}}})
    data = Int[]
    set_regex_string = Dict(:NSET  => r"NSET=([\w\-\_]+)"i,
                            :ELSET => r"ELSET=([\w\-\_]+)"i)
    selected_set = key == :NSET ? "node_sets" : "element_sets"
    definition = lines[idx_start]
    regex_string = set_regex_string[key]
    set_name = regex_match(regex_string, definition, 1)
    @info("Creating $(lowercase(string(key))) $set_name")

    if endswith(strip(uppercase(definition)), "GENERATE")
        line = lines[idx_start + 1]
        first_id, last_id, step_ = parse_numbers(line, Int)
        set_ids = collect(first_id:step_:last_id)
        push!(data, set_ids...)
    else
        for line in lines[idx_start + 1: idx_end]
            if !(empty_or_comment_line(line))
                set_ids = parse_numbers(line, Int)::Vector{Int}
                push!(data, set_ids...)
            end
        end
    end
    model[selected_set][set_name] = data
end

"""Parse SURFACE keyword
"""
function parse_section(model, lines, ::Symbol, idx_start, idx_end, ::Type{Val{:SURFACE}})
    data = Vector{Tuple{Int, Symbol}}()
    definition = lines[idx_start]

    has_set_def = parse_definition(definition)
    has_set_def != nothing || return
    set_type = get(has_set_def, "type", "UNKNOWN")
    set_name = has_set_def["name"]

    for line in lines[idx_start + 1: idx_end]
        empty_or_comment_line(line) && continue
        m = match(r"(?P<element_id>\d+),.*(?P<element_side>S\d+).*", line)
        element_id = parse(Int, m[:element_id])
        element_side = Symbol(m[:element_side])
        push!(data, (element_id, element_side))
    end
    model["surface_types"][set_name] = Symbol(set_type)
    model["surface_sets"][set_name] = data
    return
end

"""Find lines, which contain keywords, for example "*NODE"
"""
function find_keywords(lines)
    indexes = Int[]
    for (idx, line) in enumerate(lines)
        if startswith(line, "*") && !startswith(line, "**")
            push!(indexes, idx)
        end
    end
    return indexes
end

"""Main function for parsing Abaqus input file.

Function parses Abaqus input file and generates a dictionary of
all the available keywords.
"""
function parse_abaqus(fid::IOStream)
    model = Dict{String, Dict}()
    model["nodes"] = Dict{Int, Vector{Float64}}()
    model["node_sets"] = Dict{String, Vector{Int}}()
    model["elements"] = Dict{Int, Vector{Int}}()
    model["element_types"] = Dict{Int, Symbol}()
    model["element_sets"] = Dict{String, Vector{Int}}()
    model["surface_sets"] = Dict{String, Vector{Tuple{Int, Symbol}}}()
    model["surface_types"] = Dict{String, Symbol}()
    keyword_sym::Symbol = :none

    lines = readlines(fid)
    keyword_indexes = find_keywords(lines)
    nkeyword_indexes = length(keyword_indexes)
    push!(keyword_indexes, length(lines)+1)
    idx_start = keyword_indexes[1]

    for idx_end in keyword_indexes[2:end]
        keyword_line = strip(uppercase(lines[idx_start]))
        keyword = strip(regex_match(r"\s*([\w ]+)", keyword_line, 1))
        k_sym = Symbol(keyword)
        args = Tuple{Dict, Vector{Int}, Symbol, Int, Int, Type{Val{k_sym}}}
        if hasmethod(parse_section, args)
            parse_section(model, lines, k_sym, idx_start, idx_end-1, Val{k_sym})
        else
            @warn("Unknown section: '$(keyword)'")
        end
        idx_start = idx_end
    end
    return model
end

"""
    abaqus_read_mesh(fn::String)

Read ABAQUS mesh from file `fn`. Returns a dict with elements, nodes,
element sets, node sets and other topologically imporant things, but
not the actual model with boundary conditions, load steps and so on.
"""
function abaqus_read_mesh(fn::String)
    return open(parse_abaqus, fn)
end


text:
# sbf_load.jl
#
#   Load an SBF binary file, converting as much as possible in the process
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
# Julia 1.6.1
# 25.01.2018
# Last Edit: 22.07.2021

# A Labview string will have a four byte dimension array before the actual data
# s = String(read(fid,UInt8,5))
# a = Int64(read(fid,UInt16)) # Read little-endian UInt16 and convert it to Int64
# p_val = read(fid,Float64) # Read little-endian Float64

"""# sbf_load

Load SBF (Stoffflüsse Binary Format) files and convert the values if conversions are available and readable

---

### Examples

`header_info, data = sbf_load(\"K:\\\\Data\\\\Site_20171222_125500.sbf\")` Load the data from the file into a dataframe, data, converting the first column to DateTime. header_info is an array of strings representing the \'description\' of the file, the \'column names\', \'units\' for each column (m/s, K, °C, etc.), \'data types\' for each column (Float64, UInt16, etc.), and each column\'s \'conversion equations\'.\n\n
`header_info, data = sbf_load(\"K:\\\\Data\\\\Site_20171222_125500.sbf\")` Load the data from the file into a dataframe, data, converting the first column to DateTime. header_info is an array of strings representing the \'description\' of the file, the \'column names\', \'units\' for each column (m/s, K, °C, etc.), \'data types\' for each column (Float64, UInt16, etc.), and each column\'s \'conversion equations\'.\n\n
`header_info, data = sbf_load(\"Bavaria\",\"K:\\\\Data\",DateTime(2018,1,1),DateTime(2018,2,15,16))` Load all SBF files with the \"Bavaria\" prefix from the \"K:\\Data\" directory\n\n
`header_info, data = sbf_load(\"Bavaria\",\"K:\\\\Data\",DateTime(2018,1,1),DateTime(2018,2,15,16))` Load SBF files from the \"K:\\Data\" directory between January 1, 2018 and February 15, 2018 @ 16:00\n\n
`header_info, data = sbf_load(\"Bavaria\",\"K:\\\\Data\",DateTime(2018,1,1),DateTime(2018,2,15,16),verbose=true,veryverbose=true)` Load the time data from the previous example but show loading progress and column information

---

`sbf_load(F,verbose=false,veryverbose=false)`\n
* **F**::String = File name, including path
* **verbose**::Bool (optional) = Display information about the loading/processing, FALSE is default
* **veryverbose**::Bool (optional) = Display column names, units, and conversions for each file loaded, FALSE is default

`sbf_load(root_name,directory,mindate,maxdate;verbose=false,veryverbose=false)`\n
* **root_name**::String = Root name of the series of files to be loaded, \"Bavaria\" in the case of \"Bavaria_20180226_162200.sbf\"
* **directory**::String = Directory where the SBF files are locoated
* **mindate**::DateTime (optional) = Start of the time period to load, DateTime(0) is default
* **maxdate**::DateTime (optional) = End of the time period to load, DateTime(9999) is default
* **verbose**::Bool (optional) = Display information about the loading/processing, FALSE is default
* **veryverbose**::Bool (optional) = Display column names, units, and conversions for each file loaded, FALSE is default
"""
function sbf_load(root::String,src::String,mindate::DateTime=DateTime(0),maxdate::DateTime=DateTime(9999);verbose=false,veryverbose=false)
	##############
	##  Checks  ##
	##############
	!isdir(src) ? error(src * " is not a directory") : nothing

	#########################
	## List and Sort Files ##
	#########################
	verbose ? println("Listing SLT Files") : nothing
	(files,folders) = dirlist(src,regex=Regex(root * "_\\d{8}_\\d{6}\\.sbf\$"))

	# Parse File Dates
	verbose ? println("Sorting " * string(length(files)) * " Files") : nothing
	dfmt = Dates.DateFormat("yyyymmdd_HHMMSS")
	fdates = Array{DateTime}(undef,length(files))
	for i=1:1:length(files)
		fdates[i] = DateTime(files[i][end-18:end-4],dfmt)
	end

	# Remove files beyond the given time bounds
	f = findall(mindate .<= fdates .< maxdate)
	fdates = fdates[f]
	files = files[f]

	# Sort the Files By Date
	f = sortperm(fdates)
	fdates = fdates[f]
	files = files[f]

	# If no files are found return nothing
	if isempty(files)
		return Array{Any}(undef,0),DataFrame()
	end

	#################
	##  Load Data  ##
	#################
	verbose ? println("\t1: " * files[1]) : nothing
	prev_folder = dirname(files[1]) # Previous directory
	header_info, data = sbf_load(files[1],verbose=verbose,veryverbose=veryverbose)
	header_info = [[files[1];header_info]]
	for j=2:1:length(files)
		if verbose
			if dirname(files[j]) == prev_folder
				println("\t" * string(j) * ": " * basename(files[j]))
			else
				println("\t" * string(j) * ": " * files[j])
			end
		end
		header_temp, data_temp = sbf_load(files[j],verbose=verbose,veryverbose=veryverbose)
		header_temp = [files[j];header_temp]

		# Compare Headers
		if header_info[end][3] == header_temp[3] && header_info[end][4] == header_temp[4]
			data = [data;data_temp]

			if header_info[end][2] != header_temp[2] || header_info[end][5] == header_temp[5] || header_info[end][6] == header_temp[6]
				header_info = [header_info;[header_temp]]
			end
		else
			println("##  ERROR  ##")

			println(files[j])
			for k=2:1:length(header_info[end][3])
				println("\t" * header_info[end][3][k] * "(" * header_info[end][4][k] * ")")
			end

			println(files[j-1])
			for k=2:1:length(header_temp[3])
				println("\t" * header_temp[3][k] * "(" * header_temp[4][k] * ")")
			end

			error("File headers do not match. See above error information")
		end
	end

	return header_info, data
end # sbf_load(root::String,src::String;mindate::DateTime=DateTime(0),maxdate::DateTime=DateTime(9999),verbose::Bool=false)

function sbf_load(src::String;verbose::Bool=false,veryverbose::Bool=false)
	#################
	##  Constants  ##
	#################
	fid = open(src,"r")
	magicbytes = read!(fid,Array{Int8}(undef,2))
	if magicbytes != Int8[63,42]
		close(fid)
		error(src * " is not a Stoffflüsse binary (SBF) file")
	end

	########################
	##  Read Header Info  ##
	########################
	fid,description,cols,units,types,conversions,lsize = sbf_header(fid)
	types[1] = Float64 # Temporary until files are updated
	types0 = deepcopy(types) # Copy the original types for loading purposes

	#####################################
	##  Generate Conversion Functions  ##
	#####################################
	# Regular Expressions
	regex1 = r"^\w\/\-?\d+(\.\d+)?$" # Equations of the form x/1234 or x/123.423
	regex2 = r"^\-?\d+(\.\d+)?\w\+(\()\-?\d+(\.\d+)?(\))$" # Equations of the form 12.3456x+(-12.3456)
	regex2a = r"\-?\d+(\.\d+)?" # Find the number before the variable
	regex3 = r"^(\()\-?\d+(\.\d+)?\/\-?\d+(\.\d+)?(\))\w\+(\()\-?\d+(\.\d+)?(\))$" # Equations of the form (641.23523/-12.3456)x+(-12.3456)

	# Simple Functions
	dn = function doNothing(x);return x;end
	cv = function convert_to_Float64(x);return convert(Float64,x);end

	# Generated Functions
	F = fill!(Array{Function}(undef,length(conversions)),dn)
	for i=1:length(conversions)
		if occursin("timestamp",lowercase(cols[i])) || occursin("time",lowercase(cols[i]))
			if conversions[i] == "Seconds since January 1, 1904"
				veryverbose ? println("\t\t" * string(i) * ": Timestamp = Seconds since January 1, 1904") : nothing
				F[i] = function timestamp_aerodyne(t)
					DateTime(1904,1,1,1,0,0) + Dates.Second(floor(t)) + Dates.Millisecond(floor(1000(t - floor(t))))
				end
				types[i] = DateTime # For saving purposes
			end
		else
			if types[i] != Bool
				#println("F[" * string(i) * "] = function (y::" * string(types0[i]) * ")::Float64;x = convert(Float64,y);return " * conversions[i] * ";end")
				#F[i] = eval(parse("function (y::" * string(types0[i]) * ")::Float64;x = convert(Float64,y);return " * conversions[i] * ";end"))
				if occursin(regex1,replace(conversions[i]," " => ""))
					# Equations of the form x/1234 or x/123.423
					f1 = findall((in)("/"),conversions[i])
					y = Base.parse(Float64,conversions[i][f1[1] + 1:end])
					veryverbose ? println("\t\t" * string(i) * ": " * cols[i] * "(" * units[i] * ") = x/" * string(y)) : nothing
					if types[i] == Float64
						F[i] = function (x)
							return x/y
						end
					else
						F[i] = function (x)
							return convert(Float64,x)/y
						end
					end
				elseif occursin(regex2,replace(conversions[i]," " => ""))
					# Equations of the form 12.3456x+(-12.3456)
					f1 = findall((in)("("),conversions[i])
					f2 = findall((in)(")"),conversions[i])
					m = Meta.parse(match(regex2a,conversions[i]).match)
					b = Base.parse(Float64,conversions[i][f1[1] + 1:f2[1] - 1])
					veryverbose ? println("\t\t" * string(i) * ": " * cols[i] * "(" * units[i] * ") = " *string(m) * "x + (" * string(b) * ")") : nothing
					if types[i] == Float64
						F[i] = function (x)
							return m*x + b
						end
					else
						F[i] = function (x)
							return m*convert(Float64,x) + b
						end
					end
				elseif occursin(regex3,replace(conversions[i]," " => ""))
					# Equations of the form (641.23523/-12.3456)x+(-12.3456)
					f1 = findall((in)("("),conversions[i])
					f2 = findall((in)(")"),conversions[i])
					f3 = findall((in)("/"),conversions[i])
					m1 = Base.parse(Float64,conversions[i][f1[1] + 1:f3[1] - 1])
					m2 = Base.parse(Float64,conversions[i][f3[1] + 1:f2[1] - 1])
					b = Base.parse(Float64,conversions[i][f1[2] + 1:f2[2] - 1])
					veryverbose ? println("\t\t" * string(i) * ": " * cols[i] * "(" * units[i] * ") = (" *string(m1) * "/" * string(m2) * ")x + (" * string(b) * ")") : nothing
					if m2 != 0
						m = m1/m2
						if types[i] == Float64
							F[i] = function (x)
								return m*x + b
							end
						else
							F[i] = function (x)
								return m*convert(Float64,x) + b
							end
						end
					else
						verbose ? println("\t\t" * string(i) * ": " * cols[i] * "(" * units[i] * ") = x") : nothing
						F[i] = cv
					end
				else
					verbose ? println("\t\t" * string(i) * ": " * cols[i] * "(" * units[i] * ") = x") : nothing
					F[i] = cv
				end
				types[i] = Float64 # For saving purposes
			else
				F[i] = dn # Convert value to Float64
				#	F[i] = eval(parse("x -> " * conversions[i]))
			end
		end
	end

	############################
	##  Calculate line count  ##
	############################
	pos_datastart = position(fid)
	fsize = filesize(src)
	lcount = Int(floor((fsize - pos_datastart)/lsize))

	#################
	##  Load Data  ##
	#################
	#d = DataFrame(types,Symbol[Symbol(i) for i in cols],lcount,makeunique=true) # Preallocate dataframe
	d = DataFrame(Symbol[Symbol(i) for i in cols] .=> [Vector{T}(undef, lcount) for T in types],makeunique=true) # Preallocate dataframe

	colcount = length(cols)
	for i=1:lcount
		for j=1:colcount
			d[i,j] = F[j](read(fid,types0[j])) # Read and convert
		end
		read!(fid,Array{Int8}(undef,1)) # Skip the end of line character
	end

	# Completeness Check
	if !eof(fid)
		warn("\t" * src * " was closed early")
	end

	close(fid)

	return [description,cols,units,types,conversions], d
end


text:
module Surprise
# package code goes here

using Persa
using PyCall

const surprise = PyNULL()

function __init__()
    copy!(surprise, pyimport("surprise"))
end

abstract type SurpriseModel{T} <: Persa.Model{T} end

function Persa.train!(model::SurpriseModel, ds::Persa.Dataset)
    ds_surprise = transform(ds)

    model.object.fit(ds_surprise)

    return nothing
end

include("utils.jl")
include("memorybased.jl")
include("modelbased.jl")

end # module


text:
struct Firm{T <: NTuple}
  P::Function #inverse demand function
  ϕ::Function #production as a function of capital
  τ::T #tax
  α::Float64
  γ::Float64
  δ::Float64 #depreciation
  Θ::T #technology per state

  ψ::Vector{NTuple} #kroniker product of all possible state values
  π::Matrix{CartesianIndex} #policy states
  v::Matrix{Float64} #value states

  R::Array{Float64,3}
  M::Array{T,3}
end

function getlinearprice(ζ::NTuple{2,Float64})::Function
  local P::Function

  ζ[2] > 0.0 && @warn "Upward sloping demand curve! If not intentional, reparameterize ζ[1]"

  #set the function, allowing for a linear demand curve
  P(Y::Float64) = ζ[1] + ζ[2] * Y
  return P
end

#gets the cobb-douglas produciton function
function getCD(γ::Float64)::Function
  ϕ(k::Float64)::Float64 = k^γ

  return ϕ
end

function getyl(yh::Float64, k::Float64, ph::Float64, pl::Float64, F::Firm)
  kern::Float64 = (k^(F.γ*F.α)/pl) - (ph/pl)* (yh/F.Θ[1])^F.α
  yl::Float64 = (kern > 0.0) ? (kern)^(1/F.α)*F.Θ[2] : 0.0

  return yl
end

function getYlcheck(yh::Float64, k::Float64, ph::Float64, pl::Float64, F::Firm)
  Yl::Float64 = ((1/pl)*(k^(F.γ*F.α)-ph*(yh/F.Θ[1])^F.α))^(1/F.α)*F.Θ[2]
end

function getyh(yh::Float64, k::Float64, ph::Float64, pl::Float64, F::Firm)
  kern::Float64 = (k^(F.γ*F.α)/pl) - (ph/pl)* (yh/F.Θ[1])^F.α
  yl::Float64 = (kern > 0.0) ? (kern)^(1/F.α)*F.Θ[2] : 0.0

  return yl
end

function getyhcheck(yh::Float64, k::Float64, ph::Float64, pl::Float64, F::Firm)
  Yl::Float64 = ((1/pl)*(k^(F.γ*F.α)-ph*(yh/F.Θ[1])^F.α))^(1/F.α)*F.Θ[2]
end

function getkernel(pₜ::Float64, F::T where T<:Firm)::T
end

function Firm(;τ::T = DEF_τ,
  δ::Float64 = DEF_δ,
  Θ::T = DEF_Θ,
  α::Float64=DEF_α,
  γ::Float64 = DEF_γ,
  kstates::Vector{Float64} = DEF_KSTATES,
  yhstates::Vector{Float64} = DEF_YSTATES,
  ζ::T= DEF_ζ,
  p::Matrix{Float64} = DEF_p)::Firm{T} where T<:NTuple #for performance reasons, don't use vectors

  ψ::Vector{NTuple} = combine(kstates, yhstates)
  N::Int = length(ψ)
  S::Int = length(Θ)

  P::Function = getlinearprice(ζ)
  ϕ::Function = getCD(γ)

  π::Matrix{CartesianIndex} = Matrix{CartesianIndex}(undef, N, S)
  v::Matrix{Float64} = zeros(N,S)
  R::Array{Float64,3} = Array{Float64,3}(undef, N, N, S)
  M::Array{T,3} = fill(T(zeros(S)), N, N, S) #preallocate

  F::Firm{T} = Firm{T}(P, ϕ, τ, α, γ, δ, Θ, ψ, π, v, R, M)

  #println("getYl: ",getyl(1.4, 2.3, .3, .9, F)) #NOTE: Checks this formula
  #println("getYlcheck: ",getylcheck(1.4, 2.3, .3, .9, F)) #NOTE: Checks this formula

  #now build up R
  for s::Int ∈ 1:S
    for r ∈ 1:N, c ∈ 1:N
      (kₜ₁::Float64, yₜ₁s::Float64) =  ψ[r]
      (kₜ::Float64, yₜ::Float64) =  ψ[c]
      τₜ::Float64 = τ[s]

      #get the low state for the subsequent period
      yₜ₁l::Float64 = getyl(yₜ₁, kₜ₁, p[1,s], p[2,s], F::Firm)

      #need to get the current low state production if we are in the low state
      pₜ::Float64 = P(yₜ)

      dₜ::Float64 = pₜ*yₜ*(1-τₜ)-(kₜ₁-(1-δ)*kₜ)+τₜ*δ*kₜ

      F.R[r,c,s] = dₜ

    end
  end

  return F
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
@testset "pooling" begin

    # Test for agreement between CPU NNlib and CuDNN versions, across a variety of kwargs
    for num_spatial_dims in (1, 2, 3)
        # Initialize data we'll run our tests over
        C_in = 3
        batch_size = 1
        x = rand(Float64, fill(8, num_spatial_dims)..., C_in, batch_size)

        # Test that pooling is equivalent across GPU/CPU
        pdims = PoolDims(x, 2)
        y = maxpool(x, pdims)
        dy = ones(size(y))
        gputest(x -> maxpool(x, pdims), x)
        gputest((dy, y, x) -> ∇maxpool(dy, y, x, pdims), dy, y, x, checkgrad=false)
        gputest(x -> maxpool(x, pdims), x)
        gputest((dy, y, x) -> ∇maxpool(dy, y, x, pdims), dy, y, x, checkgrad=false)
    end
end


text:
using Documenter, APTDecoder

makedocs(modules = [APTDecoder], sitename = "APTDecoder.jl")

deploydocs(
    repo = "github.com/Alexander-Barth/APTDecoder.jl.git",
)

max_stars_repo_name: Alexander-Barth/APTDecoder

max_issues_repo_name: Alexander-Barth/APTDecoder

max_forks_repo_name: Alexander-Barth/APTDecoder


text:
# Wrapper around CounterValues so subtraction automatically handles wrapping
"""
    CounterValue(x::UInt)

A raw value returned from a performance counter.
This type will automatically detect and correct for counter roll-over.

```julia
julia> a = CounterTools.CounterValue(0x0)
CV(0)

julia> b = CounterTools.CounterValue(0x1)
CV(1)

julia> b - a
1

julia> UInt(a - b)
0x00007fffffffffff
```
"""
struct CounterValue
    value::UInt64
end
value(x::CounterValue) = x.value
Base.iszero(x::CounterValue) = iszero(x.value)

function Base.:-(x::CounterValue, y::CounterValue)
    # Test if overflow happend.
    start = (value(x) < value(y)) ? (UInt(1) << 48) : zero(UInt)
    if (value(x) < value(y))
        @show value(x)
        @show value(y)
    end
    return convert(Int, start + value(x) - value(y))
end

Base.show(io::IO, x::CounterValue) = print(io, "CV($(value(x)))")

#####
##### Bit maniuplation functions
#####
clearbit(x, i) = x & ~(one(x) << i)
setbit(x, i) = x | (one(x) << i)

isbitset(x, i) = !iszero(x & mask(i))

clearbits(x, i) = reduce(clearbit, i; init = x)
setbits(x, i) = reduce(setbit, i; init = x)

mask(lo::T, hi::T) where {T} = (one(lo) << (hi + 1)) - (1 << lo)
mask(i) = one(i) << i

hex(i::Integer) = string(i; base = 16)

#####
##### MSR Utilities
#####

msrpath(cpu::IndexZero) = "/dev/cpu/$(value(cpu))/msr"
msrpath(cpu::Integer) = msrpath(indexzero(cpu))
function readmsr(cpu, register)
    # Path to the kernel interface for MSRs
    path = msrpath(cpu)

    # We seek to the register and read a 64 bit int
    val = open(path) do f
        seek(f, value(indexzero(register)))
        return read(f, Int64)
    end
    return val
end

function writemsr(cpu::INDEX_TYPES, register, v)
    path = msrpath(cpu)
    open(path; write = true) do f
        seek(f, value(indexzero(register)))
        write(f, v)
    end
    return nothing
end

#####
##### Reading from Counters
#####

# Fixed counter values.
@enum FixedCounter::Int64 begin
    RetiredInstructions = 0
    ClockUnhaltedThread = 1
    ClockUnhaltedCore = 2
end

# This path goes through MSRs and is expected to be much much slower than the rdpmc
# based instructions.
readcounter(cpu, counter::INDEX_TYPES) = CounterValue(readmsr(cpu, PMC_MSRS[counter]))

"""
unsafe_rdpmc(i::Integer)

Read the contents of performance counter `i`.
If `i` is between 0 and 3, a programable counter is read.
If `i` is in 2^30, 2^30+1, or 2^30 + 2, a fixed function counter is read.

This function is unsafe since reading from an illegal value will recklessly segfault.
"""
function unsafe_rdpmc(i::Union{Integer,IndexZero,FixedCounter})
    high, low = unsafe_partial_rdpmc(i)
    return CounterValue((widen(high) << 32) | low)
end

unsafe_partial_rdpmc(i::Integer) = unsafe_partial_rdpmc(indexzero(i))
unsafe_partial_rdpmc(i::FixedCounter) = unsafe_partial_rdpmc(IndexZero(2^30 + Int64(i)))

function unsafe_partial_rdpmc(i::IndexZero)
    Base.@_inline_meta
    # This is reverse engineered from `ref.cpp` in the `ref/` directory and from Julia's
    # Tuple constructing syntax.
    #
    # The Tuple constructing syntax can be queried using
    # ```
    # f(a, b) = (a, b)
    # @code_llvm f(UInt32(0), UInt32(0))
    # ```
    val = Base.llvmcall(
        raw"""
        %val = call { i32, i32 } asm sideeffect "rdpmc", "={ax},={dx},{cx},~{dirflag},~{fpsr},~{flags}"(i32 %0)
        %low = extractvalue { i32, i32 } %val, 0
        %high = extractvalue { i32, i32 } %val, 1

        %.arr.0 = insertvalue [2 x i32] undef, i32 %high, 0
        %.arr.1 = insertvalue [2 x i32] %.arr.0, i32 %low, 1
        ret [2 x i32] %.arr.1
        """,
        Tuple{UInt32, UInt32},
        Tuple{Int32},
        convert(Int32, value(i))
    )
    return val
end


text:
function execute_query(db, query)
    if !MPI.Initialized()
        DBInterface.execute(db, query)
    else
        query = [c for c in query]
        MPI.send(query, 0, 0, MPI.COMM_WORLD)
    end
end

function load_instance_in_db!(db::SQLite.DB,
                        name::AbstractString, scenario::Int,
                        source_path::AbstractString, source_type::AbstractString,
                        date::DateTime)
    query = """
    INSERT INTO instances(name, scenario, source_path, source_type, date)
    VALUES('$name', $scenario, '$source_path', '$source_type', '$date');
    """
    execute_query(db, query)
end

function insert_decomposition!(db::SQLite.DB, origin_id, uuid,
                               origin_name, origin_scenario, extension_alg,
                               preprocess_path, date,
                               clique_path, cliquetree_path, graph_path; kwargs...)
    query = "INSERT INTO decompositions(uuid, origin_id, origin_name, origin_scenario, extension_alg, preprocess_path, date, clique_path, cliquetree_path, graph_path"
    features = [(k, v) for (k, v) in kwargs]
    for feature in features
        feature_name = feature[1]
        query *= ", $feature_name"
    end
    query *= ") VALUES('$uuid', $origin_id, '$origin_name', '$origin_scenario', '$extension_alg', '$preprocess_path', '$date', '$clique_path', '$cliquetree_path', '$graph_path'"
    for feature in features
        feature_value = feature[2]
        query *= ", $feature_value"
    end
    query *= ");"
    execute_query(db, query)
end

function insert_merge!(db::SQLite.DB, in_id::Int, out_id::Int,
                       heuristics::AbstractString, treshold_name::AbstractString,
                       treshold_percent::Float64, nb_added_edge::Int)
    query = """
    INSERT INTO merges(in_id, out_id, heuristics, treshold_name, treshold_percent, nb_added_edge)
    VALUES ($in_id, $out_id, '$heuristics', '$treshold_name', $treshold_percent, $nb_added_edge);
    """
    execute_query(db, query)
end

function insert_combination!(db::SQLite.DB, in_id1::Int, in_id2::Int, out_id::Int,
                             how::AbstractString, extension_alg::AbstractString)

    query = """
    INSERT INTO combinations(in_id1, in_id2, out_id, how, extension_alg)
    VALUES ($in_id1, $in_id2, $out_id, '$how', '$extension_alg');
    """
    execute_query(db, query)
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
struct AStarEpsilonHEntry{D <: Number}
    a_star_entry::AStarHEntry{D}
    focal_heuristic::D
end

struct CompareFocalHeap
end
compare(comp::CompareFocalHeap, e1::AStarEpsilonHEntry, e2::AStarEpsilonHEntry) = (e1.focal_heuristic, e1.a_star_entry.fvalue, -e1.a_star_entry.gvalue) < (e2.focal_heuristic, e2.a_star_entry.fvalue, -e2.a_star_entry.gvalue)

@with_kw mutable struct AStarEpsilonStates{D<:Number}
    a_star_states::AStarStates{D}                                           = AStarStates{D}()
    # Subset of heap
    focal_heap::MutableBinaryHeap{AStarEpsilonHEntry{D},CompareFocalHeap}   = MutableBinaryHeap{AStarEpsilonHEntry{D},CompareFocalHeap}()
    focal_hmap::Dict{Int64,Int64}                                           = Dict{Int64,Int64}()
    best_fvalue::D                                                          = zero(D)
    focal_heur_value::Dict{Int64,D}                                         = Dict{Int64,D}()
end

"""
Execute expand operation on the chosen node, while implicitly generating its neighbors based on a
visitor method, and populating `neighbors` with the neighboring vertices.
"""
function process_neighbors_implicit!(
    state::AStarEpsilonStates{D},
    graph::AbstractGraph{V},
    edge_wt_fn::Function,
    neighbors::Vector{Int64},
    parent_entry::AStarEpsilonHEntry{D},
    visitor::AbstractDijkstraVisitor,
    eps_weight::Float64,
    admissible_heuristic::Function,
    focal_state_heuristic::Function,
    focal_transition_heuristic::Function) where {V, D <: Number}

    dv = zero(D)
    u = parent_entry.a_star_entry.v_idx
    du = parent_entry.a_star_entry.gvalue

    for iv in neighbors

        # Default color 0
        v_color = get(state.a_star_states.colormap, iv, 0)

        if v_color == 0

            # Inserting for the first time
            state.a_star_states.dists[iv] = dv = du + edge_wt_fn(graph.vertices[u], graph.vertices[iv])
            state.a_star_states.parent_indices[iv] = u
            state.a_star_states.colormap[iv] = 1
            Graphs.discover_vertex!(visitor, graph.vertices[u], graph.vertices[iv], dv)

            new_fvalue = dv + admissible_heuristic(graph.vertices[iv])
            new_focal_heur = parent_entry.focal_heuristic + focal_state_heuristic(graph.vertices[iv]) +
                                                            focal_transition_heuristic(graph.vertices[u], graph.vertices[iv])

            new_entry = AStarEpsilonHEntry(AStarHEntry(iv, dv, new_fvalue), new_focal_heur)

            state.a_star_states.hmap[iv] = push!(state.a_star_states.heap, new_entry.a_star_entry)
            state.focal_heur_value[iv] = new_focal_heur

            # Only insert in focal list if condition satisfied
            if new_fvalue <= eps_weight * state.best_fvalue
                state.focal_hmap[iv] = push!(state.focal_heap, new_entry)
            end

        elseif v_color == 1

            dv = du + edge_wt_fn(graph.vertices[u], graph.vertices[iv])

            # Update cost-to-come if cheaper from current parent
            if dv < state.a_star_states.dists[iv]

                state.a_star_states.dists[iv] = dv
                state.a_star_states.parent_indices[iv] = u

                Graphs.update_vertex!(visitor, graph.vertices[u], graph.vertices[iv], dv)

                updated_fvalue = dv + admissible_heuristic(graph.vertices[iv])
                updated_entry = AStarEpsilonHEntry(AStarHEntry(iv, dv, updated_fvalue), parent_entry.focal_heuristic)

                old_fvalue = state.a_star_states.heap.nodes[state.a_star_states.heap.node_map[state.a_star_states.hmap[iv]]].value.a_star_entry.fvalue
                # @show old_fvalue

                # @show updated_entry
                update!(state.a_star_states.heap, state.a_star_states.hmap[iv], updated_entry.a_star_entry)

                # Enter into focal list if new fvalue is good enough
                # but it was not before
                if updated_fvalue <= eps_weight * state.best_fvalue &&
                    old_fvalue > eps_weight * state.best_fvalue
                    state.focal_hmap[iv] = push!(state.focal_heap, updated_entry)
                end

            end
        end
    end
end

"""
Runs A*-epsilon or focal search algorithm on the given list-of-indices graph.
As with all A* versions (except a_star_spath), the graph is implicit and the graph visitor must generate the neighbors of the expanded
vertex just-in-time.

Arguments:
    - `graph::AbstractGraph{V} where {V}`
    - `edge_wt_fn::F1 where {F1 <: Function}` Maps (u,v) to the edge weight of the u -> v edge
    - `source::Int64`
    - `visitor::AbstractDijkstraVisitor` The graph visitor that implements include_vertex!, which is called when
                                         a vertex is expanded
    - `eps_weight::Float64` w sub-optimality factor for focal search
    - `admissible_heuristic::F2 where {F2 <: Function}` Maps u to h(u); an admissible heuristic for the cost-to-go
    - `focal_state_heuristic::F3 where {F3 <: Function}`
    - `focal_transition_heuristic::F4 where {F4 <: Function}`
    - `::Type{D} = Float64` The type of the edge weight value

Returns:
    - `::AStarEpsilonStates` The result of the search
"""
function a_star_implicit_epsilon_path!(
    graph::AbstractGraph{V},                # the graph
    edge_wt_fn::F1, # distances associated with edges
    source::Int64,             # the source index
    visitor::AbstractDijkstraVisitor,# visitor object\
    eps_weight::Float64,
    admissible_heuristic::F2,      # Heuristic function for vertices
    focal_state_heuristic::F3,
    focal_transition_heuristic::F4,
    ::Type{D} = Float64) where {V, D <: Number, F1 <: Function, F2 <: Function, F3 <: Function, F4 <: Function}

    state = AStarEpsilonStates{D}()
    set_source!(state.a_star_states, source)
    source_heur = admissible_heuristic(graph.vertices[source])
    state.best_fvalue = source_heur

    # Initialize heap
    source_entry = AStarEpsilonHEntry(AStarHEntry(source, zero(D), source_heur), focal_state_heuristic(graph.vertices[source]))
    source_handle = push!(state.a_star_states.heap, source_entry.a_star_entry)
    state.a_star_states.hmap[source] = source_handle

    # Initialize focal heap
    focal_handle = push!(state.focal_heap, source_entry)
    state.focal_hmap[focal_handle] = source_handle

    while ~(isempty(state.a_star_states.heap))

        # Enter open-set items that are now valid into focal list
        old_best_fvalue = state.best_fvalue
        state.best_fvalue = top(state.a_star_states.heap).fvalue

        if state.best_fvalue > old_best_fvalue

            # Iterate over open set  in increasing order of fvalue and insert in focal list if valid
            for node in sort(state.a_star_states.heap.nodes, by = x->x.value.fvalue)
                fvalue = node.value.fvalue

                if fvalue > eps_weight * old_best_fvalue && fvalue <= eps_weight * state.best_fvalue
                    new_focal_entry = AStarEpsilonHEntry(node.value, state.focal_heur_value[node.value.v_idx])
                    state.focal_hmap[node.value.v_idx] = push!(state.focal_heap, new_focal_entry)
                end

                if fvalue > eps_weight * state.best_fvalue
                    break
                end
            end
        end

        # pick next vertex to include
        focal_entry, focal_handle = top_with_handle(state.focal_heap)
        heap_handle = state.a_star_states.hmap[focal_entry.a_star_entry.v_idx]

        ui = focal_entry.a_star_entry.v_idx
        du = focal_entry.a_star_entry.gvalue
        state.a_star_states.colormap[ui] = 2

        # Will be populated by include_vertex!
        nbrs = Vector{Int64}(undef, 0)

        if Graphs.include_vertex!(visitor, graph.vertices[state.a_star_states.parent_indices[ui]], graph.vertices[ui], du, nbrs) == false
            return state
        end
        # @show [graph.vertices[n] for n in nbrs]

        # Delete from open list before considering neighbors
        # TODO: Check this!!!
        pop!(state.focal_heap)
        # delete!(state.focal_heap, focal_handle)
        delete!(state.focal_hmap, focal_entry.a_star_entry.v_idx)
        delete!(state.a_star_states.heap, heap_handle)
        delete!(state.a_star_states.hmap, focal_entry.a_star_entry.v_idx)

        # process u's neighbors
        process_neighbors_implicit!(state, graph, edge_wt_fn, nbrs, focal_entry, visitor,
                                    eps_weight, admissible_heuristic, focal_state_heuristic, focal_transition_heuristic)
        Graphs.close_vertex!(visitor, graph.vertices[ui])
    end

    return state
end


text:
#=
using Pkg
Pkg.generate("jAnalogAlarmClock")
Pkg.activate("jAnalogAlarmClock")
Pkg.add(PackageSpec(url="https://github.com/JuliaGraphics/Gtk.jl", rev="master"))
Pkg.add(PackageSpec(url="https://github.com/JuliaGraphics/Luxor.jl", rev="master"))
Pkg.add(PackageSpec(url="https://github.com/JuliaGraphics/ColorSchemes.jl", rev="master"))
=#

#=
##########################
# WITH PackageCompiler.jl WORK only with modied C code
##########################
using PackageCompiler

build_executable("./jAnalogAlarmClock/src/alarmClock.jl",
"alarmClock","./jAnalogAlarmClock/src/program.c";
builddir = "./jAnalogAlarmClock/bin",
release = true, Release = true)

#############################
# WITH ApplicationBuilder.jl  BUG
/home/ygo/.julia/packages/PackageCompiler/oT98U/examples/program.c:43:5: internal compiler error: Erreur du bus
     for (i = 1; i < argc; i++) {
#############################

OLD Pkg.add(PackageSpec(url="https://github.com/NHDaly/ApplicationBuilder.jl", rev="master"))
Pkg.add(PackageSpec(url="https://github.com/JuliaLang/PackageCompiler.jl/", rev="master"))

using ApplicationBuilder
build_app_bundle("$(homedir())/julia/jAnalogAlarmClock/src/alarmClock.jl";
                 builddir="$(homedir())/julia/jAnalogAlarmClock/BINjClock" ,
                 appname="jClock", verbose=true,
                 create_installer=true)

=#
module jAnalogAlarmClock
    using Colors, Cairo, Compat, FileIO
    using Gtk
    using Luxor
    global L=Luxor
    include("drawclock.jl")
    global winx = 260
    global winy = 260
    global curcolor = "white"

    function callClock(tt)
            drawclock(60)
            Gtk.draw(c)
    end
    function mydraw()
        L.background(curcolor) # hide
        drawclock(60)
    end
    function cbresize(w)
            wdth, hght = screen_size(w)
            #wzize = get_gtk_property(w, :size_request)
            println("Windows changed $wdth $hght")
    end
    function mainwin()
        # main win
        win = GtkWindow("alarmclock",winx,winy)
        # screen = get_gtk_property(win, :screen,GdkScreen)
        # get_gtk_property(screen, :rgba_colormap, GdkColormap)
        # set_gtk_property!(win, :visual , 0.9)
        vbox = GtkBox(:v)
        #gtk canvas
        global c = Gtk.Canvas(winx,winy)
        # create luxor drawing
        global currentdrawing =  L.Drawing(winx,winy, "alarmClock.png")
        global luxctx = currentdrawing.cr

        # Define the popup menu
        popupmenu = Gtk.Menu()
        setAlarm = Gtk.MenuItem("Set Alarm")
        switchDeco = Gtk.MenuItem("Switch decoration")
        push!(popupmenu, setAlarm)
        push!(popupmenu, switchDeco)

        c.mouse.button3press = (widget,event) -> popup(popupmenu, event)

        signal_connect(setAlarm, :activate) do widget
            println("open set alarm dialog")
        end

        signal_connect(switchDeco, :activate) do widget
            if get_gtk_property(win, :decorated, Bool)
                set_gtk_property!(win, :decorated ,false)
            else
                set_gtk_property!(win, :decorated ,true)
            end
        end

        signal_connect(win, :configure_event ) do widget
            cbresize(win)
        end
        @guarded Gtk.draw(c) do widget
            ctx = Gtk.getgc(c)
            mydraw()
            Cairo.set_source_surface(ctx, currentdrawing.surface, 0, 0)
            Cairo.paint(ctx)
            Gtk.fill(ctx)
            Gtk.reveal(c)
        end

        #set_gtk_property!(win, :decorated ,false)
        set_gtk_property!(win, :opacity , 0.9)
        # set_gtk_property!(vbox, :opacity , 0.9)
        set_gtk_property!(c, :opacity , 0.9)

        push!(win, vbox)
        push!(vbox, c)

        showall(win)
        # This next line is crucial: otherwise your popup menu shows as a thin bar
        Gtk.showall(popupmenu)

        # create and set timer
        global tt = Timer(callClock, 1, interval = 1.0)

        win
    end

    # function main for static compiler
    Base.@ccallable function julia_main(ARGS::Vector{String})::Cint
        win = mainwin()
    	if !isinteractive()
    	    c = Condition()
    	    signal_connect(win, :destroy) do widget
    		notify(c)
    	    end
    	    wait(c)
    	end
        return 0
    end

# if interactive kill timer when destroy win
    if isinteractive()
        win = mainwin()
        signal_connect(win, :destroy) do widget
                close(alarmclock.tt)
        end
    end
end

max_stars_repo_name: nodrygo/jAnalogAlarmClock

max_issues_repo_name: nodrygo/jAnalogAlarmClock

max_forks_repo_name: nodrygo/jAnalogAlarmClock


text:
##############################################################################################
using Dates
buffSvr= Channel(100);

function SVrReadFileTsk(fileName )
    fi1 =open(fileName,"r")
    seek( fi1,512)
    cfgdata =readCfg(fi1)
    while(!eof( fi1))
        dataBuf = Array{UInt16,1}(undef,8192 )
        DT = readdDT(fi1)
        a= readdataType(fi1)
        # senddataType(sock)
        readindatas(fi1,1,dataBuf)
        # println(" read 1 row data ")
        put!(buffSvr, (cfgdata, DT,a,dataBuf ) )
        yield()
        # sleep(0.3) #seconds
    end
    close(fi1 )
    sleep(1000)
    print("exit out ")
end

function sendTestBufferFrame(sock,typeCode::UInt16,cfgdata::Array{UInt16,1} ,dataBuf::Array{UInt16,1}, DT,a)
    # println("---", typeCode ," ",DT ," ",length(dataBuf)," ",a)
    # println()
    SendHeadID(sock)
    sendTypeCode(sock,typeCode)
    sendcheckCode(sock)
    a= sendmachineName(sock, glbMachineName)#"TEST-567890")
    # DT = readdDT(fi1)
    a= sendDT(sock,DT)
    yield()
    sendparas(sock,cfgdata)
    yield()
    senddataType(sock)

    if typeCode ==1 || typeCode == 2
        sendDataPayLoad(sock,UInt(512),dataBuf)
    else
        sendDataPayLoad(sock,UInt(1),dataBuf)
    end
    yield()
    sendtailCheckCode(sock)
    yield()
end

function  SendoutBuffData( sock )
    # iscontinue= true
    cnt = 0
    while true
        if isready(buffSvr)
            # iscontinue = true
            cnt=0
            cfgdata, DT,a,dataBuf =take!( buffSvr )
            typeCode  = UInt16( glbtypecode)
            sendTestBufferFrame(sock,typeCode,cfgdata ,dataBuf, DT,a)
            println("$(now()),$typeCode ")
            sleep(0.2)
        else
            println("empty send Queue ")
            cnt+=1;
            sleep( 0.1)
            if cnt >100
                return
            end
        end
    end
end
function SVrSendSocketTskEntry(Port ,IPAddr)
    # println("SVrSendSocketTskEntry start ")
    if IPAddr ==""
        server = listen(IPv4(0),Port)
    else
        server = listen(IPv4(IPAddr),Port)
    end
    println(server)
    try
        while true # 1
            sock = accept(server)
            println("Client up")
            SendoutBuffData(sock)
            close(sock)
        end
    catch ex
        println( ex)
        exit(1)
    finally
        println("Except happen")
    end
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
# models.jl
# Specification of the 1D models

"""
    nondispersive_wavepacket(x0, λ)

Generates an N-unit 1D non-dispersive wavepacket
Shamelessly copied from Wikipedia:
https://en.wikipedia.org/wiki/Wave_packet

x0 - center of wavepacket
λ  - width of wavepacket

returns
ψ - complex wavefunction
"""
function nondispersive_wavepacket(x0, λ)
    ψ=[ exp(-(x-x0)^2)*(cos(2π*(x-x0)/λ) - im * sin(2π*(x-x0)/λ) ) for x=1:N ]
    return ψ
end

"""
    planewave(λ)

Generates a N-unit 1D planewave with wavelength λ.

λ - wavelength (r*a0)

returns
ψ - wavefunction
"""
function planewave(λ)
    k=2π/λ
    ψ=[ exp(-im* k*x) for x=1:N ]
    return ψ
end


text:
# Render documentation.

using Documenter
using MatrixOptim

makedocs(sitename="MatrixOptim.jl")


text:
decompose_trace(trace::Optim.OptimizationTrace) = last(trace)

function __map_optimizer_args(prob::OptimizationProblem, opt::Union{Optim.AbstractOptimizer, Optim.Fminbox,Optim.SAMIN, Optim.ConstrainedOptimizer};
    cb=nothing,
    maxiters::Union{Number, Nothing}=nothing,
    maxtime::Union{Number, Nothing}=nothing,
    abstol::Union{Number, Nothing}=nothing,
    reltol::Union{Number, Nothing}=nothing,
    kwargs...)
    if !isnothing(abstol)
        @warn "common abstol is currently not used by $(opt)"
    end

    mapped_args = (;extended_trace=true, kwargs...)

    if !isnothing(cb)
        mapped_args = (; mapped_args...,  callback = cb)
    end

    if !isnothing(maxiters)
        mapped_args = (; mapped_args..., iterations=maxiters)
    end

    if !isnothing(maxtime)
        mapped_args = (; mapped_args..., time_limit=maxtime)
    end

    if !isnothing(reltol)
        mapped_args = (; mapped_args..., f_tol=reltol)
    end

    return Optim.Options(;mapped_args...)
end

function __solve(prob::OptimizationProblem, opt::Optim.AbstractOptimizer, data = DEFAULT_DATA;
    kwargs...)
    if !isnothing(prob.lb) | !isnothing(prob.ub)
        if !(opt isa Union{Optim.Fminbox, Optim.SAMIN, Optim.AbstractConstrainedOptimizer})
            if opt isa Optim.ParticleSwarm
                opt = Optim.ParticleSwarm(;lower=prob.lb, upper=prob.ub, n_particles=opt.n_particles)
            elseif opt isa Optim.SimulatedAnnealing
                @warn "$(opt) can currently not be wrapped in Fminbox(). The lower and upper bounds thus will be ignored. Consider using a different optimizer or open an issue with Optim.jl"
            else
                opt = Optim.Fminbox(opt)
            end
        end
    end

    return ___solve(prob, opt, data; kwargs...)
end

function ___solve(prob::OptimizationProblem, opt::Optim.AbstractOptimizer,
                 data = DEFAULT_DATA;
                 cb = (args...) -> (false),
                 maxiters::Union{Number, Nothing} = nothing,
                 maxtime::Union{Number, Nothing} = nothing,
                 abstol::Union{Number, Nothing}=nothing,
                 reltol::Union{Number, Nothing}=nothing,
                 progress = false,
                 kwargs...)
    local x, cur, state

    if data != DEFAULT_DATA
        maxiters = length(data)
    end

    cur, state = iterate(data)

    function _cb(trace)
        cb_call = opt == Optim.NelderMead() ? cb(decompose_trace(trace).metadata["centroid"],x...) : cb(decompose_trace(trace).metadata["x"],x...)
        if !(typeof(cb_call) <: Bool)
            error("The callback should return a boolean `halt` for whether to stop the optimization process.")
        end
        cur, state = iterate(data, state)
        cb_call
    end

    maxiters = _check_and_convert_maxiters(maxiters)
    maxtime = _check_and_convert_maxtime(maxtime)

    f = instantiate_function(prob.f,prob.u0,prob.f.adtype,prob.p)

    !(opt isa Optim.ZerothOrderOptimizer) && f.grad === nothing && error("Use OptimizationFunction to pass the derivatives or automatically generate them with one of the autodiff backends")

    _loss = function(θ)
        x = f.f(θ, prob.p, cur...)
        __x = first(x)
        return prob.sense === MaxSense ? -__x : __x
    end

    fg! = function (G,θ)
        if G !== nothing
            f.grad(G, θ, cur...)
            if prob.sense === MaxSense
                G .*= false
            end
        end
        return _loss(θ)
    end

    if opt isa Optim.KrylovTrustRegion
        hv = function (H,θ,v)
            f.hv(H,θ,v,cur...)
            if prob.sense === MaxSense
                H .*= false
            end
        end
        optim_f = Optim.TwiceDifferentiableHV(_loss, fg!, hv, prob.u0)
    else
        gg = function (G, θ)
            f.grad(G, θ, cur...)
            if prob.sense === MaxSense
                G .*= false
            end
        end

        hh = function (H,θ)
            f.hess(H,θ,cur...)
            if prob.sense === MaxSense
                H .*= false
            end
        end
        optim_f = Optim.TwiceDifferentiable(_loss, gg, fg!, hh, prob.u0)
    end

    opt_args = _map_optimizer_args(prob, opt, cb=_cb, maxiters=maxiters, maxtime=maxtime,abstol=abstol, reltol=reltol; kwargs...)

    t0 = time()
    opt_res = Optim.optimize(optim_f, prob.u0, opt, opt_args)
    t1 = time()
    opt_ret = Symbol(Optim.converged(opt_res))

    SciMLBase.build_solution(prob, opt, opt_res.minimizer, prob.sense === MaxSense ? -opt_res.minimum : opt_res.minimum; original=opt_res, retcode=opt_ret)
end

function ___solve(prob::OptimizationProblem, opt::Union{Optim.Fminbox,Optim.SAMIN},
                 data = DEFAULT_DATA;
                 cb = (args...) -> (false),
                 maxiters::Union{Number, Nothing} = nothing,
                 maxtime::Union{Number, Nothing} = nothing,
                 abstol::Union{Number, Nothing}=nothing,
                 reltol::Union{Number, Nothing}=nothing,
                 progress = false,
                 kwargs...)

    local x, cur, state

    if data != DEFAULT_DATA
        maxiters = length(data)
    end

    cur, state = iterate(data)

    function _cb(trace)
        cb_call = !(opt isa Optim.SAMIN) && opt.method == Optim.NelderMead() ? cb(decompose_trace(trace).metadata["centroid"],x...) : cb(decompose_trace(trace).metadata["x"],x...)
        if !(typeof(cb_call) <: Bool)
            error("The callback should return a boolean `halt` for whether to stop the optimization process.")
        end
        cur, state = iterate(data, state)
        cb_call
    end

    maxiters = _check_and_convert_maxiters(maxiters)
    maxtime = _check_and_convert_maxtime(maxtime)

    f = instantiate_function(prob.f,prob.u0,prob.f.adtype,prob.p)

    !(opt isa Optim.ZerothOrderOptimizer) && f.grad === nothing && error("Use OptimizationFunction to pass the derivatives or automatically generate them with one of the autodiff backends")

    _loss = function(θ)
        x = f.f(θ, prob.p, cur...)
        __x = first(x)
        return prob.sense === MaxSense ? -__x : __x
    end
    fg! = function (G,θ)
        if G !== nothing
            f.grad(G, θ, cur...)
            if prob.sense === MaxSense
                G .*= false
            end
        end
        return _loss(θ)
    end

    gg = function (G, θ)
        f.grad(G, θ, cur...)
        if prob.sense === MaxSense
            G .*= false
        end
    end
    optim_f = Optim.OnceDifferentiable(_loss, gg, fg!, prob.u0)

    opt_args = _map_optimizer_args(prob, opt, cb=_cb, maxiters=maxiters, maxtime=maxtime,abstol=abstol, reltol=reltol; kwargs...)

    t0 = time()
    opt_res = Optim.optimize(optim_f, prob.lb, prob.ub, prob.u0, opt, opt_args)
    t1 = time()
    opt_ret = Symbol(Optim.converged(opt_res))

    SciMLBase.build_solution(prob, opt, opt_res.minimizer, opt_res.minimum; original=opt_res, retcode=opt_ret)
end

function ___solve(prob::OptimizationProblem, opt::Optim.ConstrainedOptimizer,
                 data = DEFAULT_DATA;
                 cb = (args...) -> (false),
                 maxiters::Union{Number, Nothing} = nothing,
                 maxtime::Union{Number, Nothing} = nothing,
                 abstol::Union{Number, Nothing}=nothing,
                 reltol::Union{Number, Nothing}=nothing,
                 progress = false,
                 kwargs...)

    local x, cur, state

    if data != DEFAULT_DATA
        maxiters = length(data)
    end

    cur, state = iterate(data)

      function _cb(trace)
      cb_call = cb(decompose_trace(trace).metadata["x"],x...)
      if !(typeof(cb_call) <: Bool)
          error("The callback should return a boolean `halt` for whether to stop the optimization process.")
      end
      cur, state = iterate(data, state)
      cb_call
    end

    maxiters = _check_and_convert_maxiters(maxiters)
    maxtime = _check_and_convert_maxtime(maxtime)

    f = instantiate_function(prob.f,prob.u0,prob.f.adtype,prob.p,prob.ucons === nothing ? 0 : length(prob.ucons))

    f.cons_j ===nothing && error("This optimizer requires derivative definitions for nonlinear constraints. If the problem does not have nonlinear constraints, choose a different optimizer. Otherwise define the derivative for cons using OptimizationFunction either directly or automatically generate them with one of the autodiff backends")

    _loss = function(θ)
        x = f.f(θ, prob.p, cur...)
        __x = first(x)
        return prob.sense === MaxSense ? -__x : __x
    end
    fg! = function (G,θ)
        if G !== nothing
            f.grad(G, θ, cur...)
            if prob.sense === MaxSense
                G .*= false
            end
        end
        return _loss(θ)
    end
    gg = function (G, θ)
        f.grad(G, θ, cur...)
        if prob.sense === MaxSense
            G .*= false
        end
    end

    hh = function (H,θ)
        f.hess(H,θ,cur...)
        if prob.sense === MaxSense
            H .*= false
        end
    end
    optim_f = Optim.TwiceDifferentiable(_loss, gg, fg!, hh, prob.u0)

    cons! = (res, θ) -> res .= f.cons(θ);

    cons_j! = function(J, x)
        f.cons_j(J, x)
    end

    cons_hl! = function (h, θ, λ)
        res = [similar(h) for i in 1:length(λ)]
        f.cons_h(res, θ)
        for i in 1:length(λ)
            h .+= λ[i]*res[i]
        end
    end

    lb = prob.lb === nothing ? [] : prob.lb
    ub = prob.ub === nothing ? [] : prob.ub
    optim_fc = Optim.TwiceDifferentiableConstraints(cons!, cons_j!, cons_hl!, lb, ub, prob.lcons, prob.ucons)

    opt_args = _map_optimizer_args(prob, opt, cb=_cb, maxiters=maxiters, maxtime=maxtime,abstol=abstol, reltol=reltol; kwargs...)

    t0 = time()
    opt_res = Optim.optimize(optim_f, optim_fc, prob.u0, opt, opt_args)
    t1 = time()
    opt_ret = Symbol(Optim.converged(opt_res))

    SciMLBase.build_solution(prob, opt, opt_res.minimizer, opt_res.minimum; original=opt_res, retcode=opt_ret)
end


text:
"""
    SphereSphereCollision

    collision between two spheres

    origin_parent: position of contact on parent body relative to center of mass
    origin_child: position of contact on parent body relative to center of mass
    radius_parent: radius of contact for parent body
    radius_child: radius of contact for child body
"""
mutable struct SphereSphereCollision{T,O,I,OI} <: Collision{T,O,I,OI}
    origin_parent::SVector{I,T}
    origin_child::SVector{I,T}
    radius_parent::T
    radius_child::T
end

# function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, collision::SphereSphereCollision)
#     summary(io, collision)
#     println(io, "")
#     println(io, "origin_parent: "*string(collision.origin_parent))
#     println(io, "origin_child:  "*string(collision.origin_child))
#     println(io, "radius_parent:  "*string(collision.radius_parent))
#     println(io, "radius_child:  "*string(collision.radius_child))
# end

# distance
function distance(collision::SphereSphereCollision, xp, qp, xc, qc)
    # contact origin points
    cop = contact_point_origin(xp, qp, collision.origin_parent)
    coc = contact_point_origin(xc, qc, collision.origin_child)

    # distance between contact origins
    d = norm(cop - coc, 2)

    # minimum distance between spheres
    return d - (collision.radius_parent + collision.radius_child)
end

function ∂distance∂x(gradient::Symbol, collision::SphereSphereCollision, xp, qp, xc, qc)
    # contact origin points
    cop = contact_point_origin(xp, qp, collision.origin_parent)
    coc = contact_point_origin(xc, qc, collision.origin_child)

    # distance between contact origins
    d = norm(cop - coc, 2)
    ∂norm∂d = ∂norm∂x(cop - coc)

    if gradient == :parent
        D = ∂norm∂d *  1.0 * ∂contact_point_origin∂x(xp, qp, collision.origin_parent)
    elseif gradient == :child
        D = ∂norm∂d * -1.0 * ∂contact_point_origin∂x(xc, qc, collision.origin_child)
    end

    if gradient == :parent
        FD = FiniteDiff.finite_difference_jacobian(x -> distance(collision, x, qp, xc, qc), xp)
    elseif gradient == :child
        FD = FiniteDiff.finite_difference_jacobian(x -> distance(collision, xp, qp, x, qc), xc)
    end

    return FD

    # @assert norm(D - FD, Inf) < 1.0e-5

    return D
end

function ∂distance∂q(gradient::Symbol, collision::SphereSphereCollision, xp, qp, xc, qc)
    # contact origin points
    cop = contact_point_origin(xp, qp, collision.origin_parent)
    coc = contact_point_origin(xc, qc, collision.origin_child)

    # distance between contact origins
    d = norm(cop - coc, 2)
    ∂norm∂d = ∂norm∂x(cop - coc)

    if gradient == :parent
        D = ∂norm∂d *  1.0 * ∂contact_point_origin∂q(xp, qp, collision.origin_parent)
    elseif gradient == :child
        D = ∂norm∂d * -1.0 * ∂contact_point_origin∂q(xc, qc, collision.origin_child)
    end

    if gradient == :parent
        FD = FiniteDiff.finite_difference_jacobian(q -> distance(collision, xp, Quaternion(q..., false), xc, qc), vector(qp))
    elseif gradient == :child
        FD = FiniteDiff.finite_difference_jacobian(q -> distance(collision, xp, qp, xc, Quaternion(q..., false)), vector(qc))
    end

    return FD
    # @assert norm(D - FD, Inf) < 1.0e-5

    return D
end

# contact point in world frame
function contact_point(relative::Symbol, collision::SphereSphereCollision, xp, qp, xc, qc)
    # contact origin points
    cop = contact_point_origin(xp, qp, collision.origin_parent)
    coc = contact_point_origin(xc, qc, collision.origin_child)

    # direction of minimum distance (child to parent)
    d = cop - coc
    dir = normalize(d)

    # contact point
    if relative == :parent
        return cop - collision.radius_parent * dir
    elseif relative == :child
        return coc + collision.radius_child * dir
    end
end

function ∂contact_point∂x(relative::Symbol, jacobian::Symbol, collision::SphereSphereCollision, xp, qp, xc, qc)
    # contact origin points
    cop = contact_point_origin(xp, qp, collision.origin_parent)
    coc = contact_point_origin(xc, qc, collision.origin_child)

    # direction of minimum distance (child to parent)
    d = cop - coc
    dir = normalize(d)

    if relative == :parent
        # cop - collision.radius_parent * dir
        if jacobian == :parent
            ∂c∂x = ∂contact_point_origin∂x(xp, qp, collision.origin_parent)
            X = ∂c∂x
            X -= collision.radius_parent * ∂normalize∂x(d) * ∂c∂x
        elseif jacobian == :child
            X = -1.0 * collision.radius_parent * ∂normalize∂x(d) * -1.0 * ∂contact_point_origin∂x(xc, qc, collision.origin_child)
        end
    elseif relative == :child
        # coc + collision.radius_child * dir
        if jacobian == :parent
            X = 1.0 * collision.radius_child * ∂normalize∂x(d) * ∂contact_point_origin∂x(xp, qp, collision.origin_parent)
        elseif jacobian == :child
            ∂c∂x = ∂contact_point_origin∂x(xc, qc, collision.origin_child)
            X = ∂c∂x
            X += collision.radius_child * ∂normalize∂x(d) * -1.0 * ∂c∂x
        end
    end

    if jacobian == :parent
        FD =  FiniteDiff.finite_difference_jacobian(x -> contact_point(relative, collision, x, qp, xc, qc), xp)
    elseif jacobian == :child
        FD = FiniteDiff.finite_difference_jacobian(x -> contact_point(relative, collision, xp, qp, x, qc), xc)
    end

    return FD
    # @assert norm(X - FD, Inf) < 1.0e-5

    return X
end

function ∂contact_point∂q(relative::Symbol, jacobian::Symbol, collision::SphereSphereCollision, xp, qp, xc, qc)
    # contact origin points
    cop = contact_point_origin(xp, qp, collision.origin_parent)
    coc = contact_point_origin(xc, qc, collision.origin_child)

    # direction of minimum distance (child to parent)
    d = cop - coc
    dir = normalize(d)

    if relative == :parent
        # cop - collision.radius_parent * dir
        if jacobian == :parent
            ∂c∂q = ∂contact_point_origin∂q(xp, qp, collision.origin_parent)
            Q = ∂c∂q
            Q -= collision.radius_parent * ∂normalize∂x(d) * ∂c∂q
        elseif jacobian == :child
            Q = -1.0 * collision.radius_parent * ∂normalize∂x(d) * -1.0 * ∂contact_point_origin∂q(xc, qc, collision.origin_child)
        end
    elseif relative == :child
        # coc + collision.radius_child * dir
        if jacobian == :parent
            Q = 1.0 * collision.radius_child * ∂normalize∂x(d) * ∂contact_point_origin∂q(xp, qp, collision.origin_parent)
        elseif jacobian == :child
            ∂c∂q = ∂contact_point_origin∂q(xc, qc, collision.origin_child)
            Q = ∂c∂q
            Q += collision.radius_child * ∂normalize∂x(d) * -1.0 * ∂c∂q
        end
    end

    if jacobian == :parent
        FD = FiniteDiff.finite_difference_jacobian(q -> contact_point(relative, collision, xp, Quaternion(q..., false), xc, qc), vector(qp))
    elseif jacobian == :child
        FD = FiniteDiff.finite_difference_jacobian(q -> contact_point(relative, collision, xp, qp, xc, Quaternion(q..., false)), vector(qc))
    end

    return FD

    # @assert norm(Q - FD, Inf) < 1.0e-5

    return Q
end


text:
#### Types

"""
    AbstractParameterSet

The top-level super-type parameter set.
"""
abstract type AbstractParameterSet end
const APS = AbstractParameterSet

"""
    AbstractEarthParameterSet <: AbstractParameterSet

An earth parameter set, specific to planet Earth.
"""
abstract type AbstractEarthParameterSet <: AbstractParameterSet end
const AEPS = AbstractEarthParameterSet


text:
"""
$(TYPEDEF)
$(TYPEDFIELDS)

"""
struct QuantLogUniform <: LogUniformQuantDist
    low::NestedFloat
    high::NestedFloat
    q::NestedFloat
end

function quantloguniform(low::Float64, high::Float64, q::Float64)::Float64
    return round(exp(uniform(low, high)) / q) * q
end

loguniformquant(x::QuantLogUniform, low, high, q) = quantloguniform(low, high, q)