text:
## Script to estimate the rotational transform between camera and IMU

using SensorFeatureTracking, TransformUtils, FixedPointNumbers, Colors, Plots
using ProfileView

projdir = joinpath(dirname(@__FILE__), "..")
# srcdir = joinpath(projdir,"src")
datadir = joinpath(projdir,"Data")

## Extract example data
# Extract example imu data
cd(dirname(@__FILE__))

imudata = Vector{IMU_DATA}()

imu_data,imu_fields = readcsv("../Data/testSequence_IMU.csv",header=true)

for i = 1:size(imu_data,1)
	push!(imudata,IMU_DATA(Int64(imu_data[i,1]),imu_data[i,9:11],imu_data[i,3:5]))
end

# read image log file and get timestaps
timestamps = readcsv("../Data/testSequence_CAM.csv", Int64, skipstart=4)[:,2]
N = size(timestamps,1)

## Set up camera model
focald = 524.040
cu = 319.254
cv = 251.227
cam = CameraModel(640,480,[focald, focald],[cv, cu], 0., [0])

##
# go through frames and calculate delta rotations and save as Quaternion
index = PInt64(1)
ΔQ_s = Vector{Quaternion}(N)
trackAX = Vector{Vector{Float64}}(N)
gyroAX = Vector{Vector{Float64}}(N)

for i = 1:(N)#length(timestamps)
	(valid, R) = integrateGyroBetweenFrames!(index,timestamps[i],imudata)
	ΔQ_s[i] = convert(Quaternion,SO3(R))

	gyroAX[i]  = ΔQ_s[i].v ./ norm(ΔQ_s[i].v)
	trackAX[i] = [NaN, NaN, NaN]
end

## Track features and save delta rotations between frames, Block trakcer for now
startFrame = 50
numFrames = 50

im1 = load(joinpath(datadir,"testSequence/image_$(startFrame).jpg"))
im2 = load(joinpath(datadir,"testSequence/image_$(startFrame+1).jpg"))

(ro,co) = size(im1)
feats = getApproxBestShiTomasi(im1,nfeatures=100, stepguess=0.95)

flow = BlockTracker(feats, matchFunction = compute_ssd)#flow_value_threshold = 50000

image1 = deepcopy(im1)
oldfeats = deepcopy(flow.features)

block_tracker!(flow, im1, im2)

# blank image for feature overlay
blankImg = zeros(Gray{N0f8},ro,co)
map((ft1,ft2)->drawfeatureLine!(blankImg, Feature(ft1),Feature(ft2)),oldfeats,flow.features)

# there must be at least 4 valid keypoints to track
q = estimateRotationFromKeypoints(oldfeats,flow.features, cam)

trackAX[startFrame+1] = q.v./norm(q.v)

#
for frame_idx = startFrame+2:startFrame+numFrames -1
    println("Processing Frame: ", frame_idx)
    im1 = deepcopy(im2)
    im2 = load(joinpath(datadir,"testSequence/image_$(frame_idx).jpg"))

    # im2 = im2[range...]
    block_tracker!(flow, im1, im2)

    map((ft1,ft2)->drawfeatureLine!(blankImg, Feature(ft1),Feature(ft2)),flow.features, oldfeats)

	q = estimateRotationFromKeypoints(oldfeats,flow.features, cam)
	trackAX[frame_idx] = q.v./norm(q.v)

    oldfeats = deepcopy(flow.features)
end

# overlay tracking traces with image
# o_img = (reinterpret(UInt8,blankImg)) .| (reinterpret(UInt8,im2))
# b_img = reinterpret(Gray{N0f8},o_img)
o_img = (reinterpret(UInt8,blankImg)) .| (reinterpret(UInt8,im2))
r_img = reinterpret(Gray{N0f8},o_img)
n_img = .~(reinterpret(UInt8,blankImg)) .& (reinterpret(UInt8,im2))
n_img = reinterpret(Gray{N0f8},n_img)
#
# colorview(RGB, n_img, n_img, b_img)
# colorview(RGB, n_img, g_img, n_img)
colorview(RGB, r_img, n_img, n_img)

## finally compute angle between camera frame and imu frame using HornAbsoluteOrientation

# Compute best fit alignment between all these vectors
oQc = HornAbsoluteOrientation(gyroAX, trackAX)

## Rotate trackAX into alignment with IMU Axis to test
# TODO: does not seem to be working
# iQc = Quaternion()
# 	iQc.s = 0.697815224446998;
# 	iQc.v = [0.715977169639074, 0.014522897934957, 0.014821960773685]

trackAX_b = TransformUtils.rotate.(oQc, trackAX)

# repack in matrix for plotting TODO: I'm sure Vectors can somehow be plotted?
trackAX_b_m = zeros(N,3)
gyroAX_m = zeros(N,3)
trackAX_m = zeros(N,3)
for i = 1:N
	trackAX_b_m[i,:] = trackAX_b[i]
	gyroAX_m[i,:] = gyroAX[i]
	trackAX_m[i,:] = trackAX[i]
end
##
p1 = plot(gyroAX_m[startFrame:startFrame+numFrames,:] - trackAX_b_m[startFrame:startFrame+numFrames,:])
##
# p2 = plot(trackAX_b_m[startFrame:startFrame+numFrames,:])
##
p2 = plot(gyroAX_m[startFrame:startFrame+numFrames,:] - trackAX_m[startFrame:startFrame+numFrames,:])
##
plot(p1,p2,layout=(2,1))

max_stars_repo_name: dehann/SensorFeatureTracking

max_issues_repo_name: dehann/SensorFeatureTracking

max_forks_repo_name: dehann/SensorFeatureTracking


text:
module Smooth

using ..PythonOT: PythonOT
using ..PyCall: PyCall

export smooth_ot_dual

"""
    smooth_ot_dual(μ, ν, C, ε; reg_type="l2", kwargs...)

Compute the optimal transport plan for a regularized optimal transport problem
with source and target marginals `μ` and `ν`, cost matrix `C` of size
`(length(μ), length(ν))`, and regularization parameter `ε`.

The optimal transport map `γ` is of the same size as `C` and solves
```math
\\inf_{\\gamma \\in \\Pi(\\mu, \\nu)} \\langle \\gamma, C \\rangle
+ \\varepsilon \\Omega(\\gamma),
```
where ``\\Omega(\\gamma)`` is the L2-regularization term
``\\Omega(\\gamma) = \\|\\gamma\\|_F^2/2`` if `reg_type="l2"` (the default) or
the entropic regularization term
``\\Omega(\\gamma) = \\sum_{i,j} \\gamma_{i,j} \\log \\gamma_{i,j}`` if `reg_type="kl"`.

The function solves the dual formulation[^BSR2018]
```math
\\max_{\\alpha, \\beta} \\mu^{\\mathsf{T}} \\alpha + \\nu^{\\mathsf{T}} \\beta −
\\sum_{j} \\delta_{\\Omega}(\\alpha + \\beta_j - C_j),
```
where ``C_j`` is the ``j``th column of the cost matrix and ``\\delta_{\\Omega}`` is the
conjugate of the regularization term ``\\Omega``.

This function is a wrapper of the function
[`smooth_ot_dual`](https://pythonot.github.io/gen_modules/ot.smooth.html#ot.smooth.smooth_ot_dual)
in the Python Optimal Transport package. Keyword arguments are listed in the documentation
of the Python function.

# Examples

```jldoctest; setup=:(using PythonOT.Smooth)
julia> μ = [0.5, 0.2, 0.3];

julia> ν = [0.0, 1.0];

julia> C = [0.0 1.0; 2.0 0.0; 0.5 1.5];

julia> smooth_ot_dual(μ, ν, C, 0.01)
3×2 Matrix{Float64}:
 0.0  0.5
 0.0  0.2
 0.0  0.300001
```

[^BSR2018]: Blondel, M., Seguy, V., & Rolet, A. (2018). Smooth and Sparse Optimal Transport. In *Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics (AISTATS)*.
"""
function smooth_ot_dual(μ, ν, C, ε; kwargs...)
    return PythonOT.pot.smooth.smooth_ot_dual(μ, ν, C, ε; kwargs...)
end
end

ext: jl


text:
module ObjCWrap

include("primitives.jl")
include("methods.jl")
include("call.jl")
include("syntax.jl")

include("foundation.jl")
include("classes.jl")

end


text:
@doc """
High pass filter applied in forward and reverse direction

Simply a wrapper for the DSP.jl functions

#### Arguments

* `signals`: Signal data in the format samples x channels
* `cutOff`: Cut off frequency in Hz
* `fs`: Sampling rate
* `order`: Filter orde

#### Returns

* filtered signal
* filter used on signal
""" ->
function highpass_filter{T <: AbstractFloat}(signals::Array{T}, cutOff::Number, fs::Number, order::Int)
    debug("Highpass filtering $(size(signals)[end]) channels.  Pass band > $(cutOff) Hz")
    Wn = cutOff/(fs/2)
    highpass_filter(signals, Wn, order)
end

function highpass_filter{T <: AbstractFloat}(signals::Array{T}, Wn::Number, order::Int)
    debug("Filter order = $order, Wn = $Wn")
    f = digitalfilter(Highpass(Wn), Butterworth(order))
    signals = filtfilt(f, signals)
    return signals, f
end

@doc """
Low pass filter applied in forward and reverse direction

Simply a wrapper for the DSP.jl functions

#### Input

* `signals`: Signal data in the format samples x channels
* `cutOff`: Cut off frequency in Hz
* `fs`: Sampling rate
* `order`: Filter orde

#### Output

* filtered signal
* filter used on signal
""" ->
function lowpass_filter{T <: AbstractFloat}(signals::Array{T}, cutOff::Number, fs::Number, order::Int)
    debug("Lowpass filtering $(size(signals)[end]) channels.  Pass band < $(cutOff) Hz")
    Wn = cutOff/(fs/2)
    lowpass_filter(signals, Wn, order)
end

function lowpass_filter{T <: AbstractFloat}(signals::Array{T}, Wn::Number, order::Int)
    debug("Filter order = $order, Wn = $Wn")
    f = digitalfilter(Lowpass(Wn), Butterworth(order))
    signals = filtfilt(f, signals)
    return signals, f
end

@doc """
Band pass filter
""" ->
function bandpass_filter(signals::Array, lower::Number, upper::Number, fs::Number, n::Int, rp::Number)
    # Type 1 Chebychev filter
    # TODO filtfilt does not work. Why not?

    signals = convert(Array{Float64}, signals)

    f = digitalfilter(Bandpass(lower, upper, fs=fs), Chebyshev1(n, rp))

    Logging.info("Bandpass filtering $(size(signals)[end]) channels.     $lower < Hz < $upper")
    debug("Filter order = $n, fs = $fs")

    signals = filt(f, signals)
    signals = filt(f, flipdim(signals, 1))
    signals = flipdim(signals, 1)

    return signals, f
end

#######################################
#
# Filter compensation
#
#######################################

function compensate_for_filter(d::Dict, spectrum::AbstractArray, fs::Real)
    frequencies = linspace(0, 1, Int(size(spectrum, 1))) * fs / 2

    key_name = "filter"
    key_numb = 1
    key      = string(key_name, key_numb)

    while haskey(d, key)

        spectrum = compensate_for_filter(d[key], spectrum, frequencies, fs)

        debug("Accounted for $key response in spectrum estimation")

        key_numb += 1
        key = string(key_name, key_numb)
    end

    return spectrum
end

@doc """
Recover the spectrum of signal by compensating for filtering done.

#### Arguments

* `filter`: The filter used on the spectrum
* `spectrum`: Spectrum of signal
* `frequencies`: Array of frequencies you want to apply the compensation to
* `fs`: Sampling rate

#### Returns

Spectrum of the signal after comensating for the filter

#### TODO

Extend this to arbitrary number of dimensions rather than the hard coded 3
""" ->
function compensate_for_filter(filter::FilterCoefficients, spectrum::AbstractArray, frequencies::AbstractArray, fs::Real)
    filter_response     = freqz(filter, frequencies, fs)

    for f = 1:length(filter_response)
        spectrum[f, :, :] = spectrum[f, :, :] ./ abs(filter_response[f])^2
    end

    return spectrum
end


text:
"""
    Realized firm profits taking the interest rate R as given
"""
function firm_profit(a, K; par::Pars)
    @unpack_Pars par
    Y^(1/θ)*exp(a)^((θ - 1)/θ)*K^α_tilde - R*K
end

function Eπ(s, K; par::Pars, bpar::BPars)
    @unpack_Pars par
    ea = EA(s, par, bpar)
    Y^(1/θ)*ea*K^α_tilde - R*K
end

"""
    Constructor for the Bayesian Parameters (weights on private signal and prior)
"""
function BPars(τ, par::Pars)
    @unpack_Pars par
    τ   = τ
    σ²  = 1/(τ + σ²ₐ^-1) # posterior uncertainty
    w_a = τ*σ²           # weight on private signal
    w_p = σ²ₐ^-1*σ²      # weight on prior

    BPars(τ   = τ,
          w_a = w_a,
          w_p = w_p,
          σ²  = σ²)
end

signal(a, eps, bpar::BPars) = a + eps/√bpar.τ
signal(a, eps, τ) = a + eps/√τ
Tau_of_K_old(K) = K

function EA(a, par::Pars, bpar::BPars)
    @unpack_Pars par
    @unpack_BPars bpar
    exp(((θ - 1)/θ)*w_a*a + ((θ - 1)/θ)^2*σ²/2)
end

function capital_demand(a, par::Pars, bpar::BPars)
    @unpack_Pars par
    @unpack_BPars bpar
    (Y^(1/θ)*EA(a, par, bpar)/R)^(1/(1-α_tilde))
end

extractK(pos; K_grid) = K_grid[pos[2]]

function solve_model(par::Pars; V_guess = nothing, conv = 1e-8)
    @unpack_Pars par

    V_guess    === nothing ? V_guess = similar(K_grid) : nothing
    K_sol      = zeros(a_n^2,K_n)
    V_new      = zeros(K_n)
    V_next_big = repeat(V_guess, inner = a_n^2)
    nodes, w   = qnwnorm(a_n, [par.a_t, 0], [par.σ²ₐ, 1])
    a          = nodes[:,1]
    ε          = nodes[:,2]
    diff       = 1e8
    while abs(diff) > conv
        V_new   = iterate(V_guess, K_sol, V_new, a, ε, w, V_next_big, par)
        diff    = sum((V_guess - V_new).^2)
        @show diff
        V_guess = copy(V_new)
    end
    return V_new
end

function iterate(V_guess, K_sol, V_new, a, ε, w, V_next_big, par::Pars)
    @unpack_Pars par
    for i in 1:K_n
        K_old          = K_grid[i]
        τ              = Tau_of_K_old(K_old)
        s              = signal.(a, ε, τ)
        grid           = gridmake(s, K_grid)
        W_big          = Eπ.(grid[:,1], grid[:,2]; par, bpar = BPars(τ, par)) .+ (1-δ).*K_old - γ/2*grid[:,2].^2/K_old + R^-1*V_next_big
        W_big_reshaped = reshape(W_big, (a_n^2, K_n))
        W, pos         = findmax(W_big_reshaped, dims = 2)
        K_sol[:,i]     = extractK.(pos; K_grid)
        #V_new_big[:,i] = sum(w.*(firm_profit.(a, K_sol[:,i]; par)[:] + par.R^-1*extractK.(pos; K_grid = V_guess)[:]))
        V_new[i]       = sum(w.*(firm_profit.(a, K_sol[:,i]; par)[:] .+ (1-δ).*K_old - γ/2*K_sol[:,i].^2/K_old + R^-1*extractK.(pos; K_grid = V_guess)[:]))
    end
    return V_new
end

function K_solve(s, V, K_old; par::Pars)
    @unpack_Pars par
    W_big  = Eπ.(s, K_grid; par, bpar = BPars(Tau_of_K_old(K_old), par)) .+ (1-δ).*K_old - γ/2*K_grid.^2/K_old + R^-1*V
    W, pos = findmax(W_big)

    return K_grid[pos]
end


text:
module testliquid

using OnlineStats, Qwt, Distributions, StatsBase, OnlineAI

# create data

nin = 1
nout = 1
T = 2000
# x = collect(linspace(-15., 15., T))

# x is a simple differenced AR(1)... y is the future val
x = randn(T)
for i in 2:T
  x[i] += x[i-1] * 0.95
end
# dx = diff(x)
lookahead = 100
y = x
# y = x[lookahead:end]
# x = x[1:end-lookahead]
T = length(x)

# create liquid state machine
params = LiquidParams(λ = 1.0,
                      w = 8,
                      h = 4,
                      decayRateDist = Uniform(0.99, 1.0),
                      pctInput = 0.2,
                      pctOutput = 0.8,
                      readout = FireWindowReadout(30),
                      baseThreshold = 1.0)
                      # readout = StateReadout())
lsm = LiquidStateMachine(params, nin, nout)
# liquid = lsm.liquid

# set up visualization
viz = visualize(lsm)

# fit the lsm
plotiter = 10
# anim = animation(viz.window, "/Users/tom/Pictures/gifs")
for (i,t) = enumerate(lookahead:T)
  fit!(lsm, vec(x[t-lookahead+1,:]), vec(y[t,:]))

  if i % plotiter == 0
    fit!(viz, vec(y[t,:]))
    # saveframe(anim)
  end
end
# makegif(anim)

# # plot the grf functions
# grfs = GaussianReceptiveField[input.grf for input in lsm.input.inputs[1,:]]
# plt3 = plot(grfs, -15.:.1:15.)

end

tl = testliquid
# lsm = tl.lsm


text:
"""
    Dataset

structure holding inner representation of input data samples and provides replacement of missing values so they can be used by universal solver

# Fields
- `passengers`: number of data samples present
- `ids`: id numbers assigned to individual samples to distinguish them easily
- `X`: matrix of real numbers where one column represents features of a single data sample
"""
struct Dataset
    passengers::Integer
    ids::Vector{<:Integer}
    X::Matrix{<:Real}

    """
        Dataset(data::DataFrame; normalize=true)

    constucts Dataset object from `DataFrame`

    when keyword argument `normalize` is set to ``true``, the data are normalized for better function of some classifiers
    """
    function Dataset(data::DataFrame; normalize=true)
        ids = data[!, :PassengerId]
        passengers = length(ids)
        X = ones(8, passengers)

        sex = ones(Real, passengers)
        sex[data[!, :Sex] .== "female"] .= -1
        X[2, :] = sex

        X[3, :] = data[!, :Pclass]

        age = data[!, :Age]
        avgage = round(sum(skipmissing(age)) / (passengers - sum(ismissing.(age))); digits=2)
        age[ismissing.(age)] .= avgage
        X[4, :] = age

        X[5, :] = data[!, :SibSp]
        X[6, :] = data[!, :Parch]
        fare = data[!, :Fare]
        avgfare = round(sum(skipmissing(fare)) / (passengers - sum(ismissing.(fare))); digits=2)
        fare[ismissing.(fare)] .= avgfare
        X[7, :] = data[!, :Fare]

        embarked = data[!, :Embarked]
        embarkedprime = skipmissing(embarked)
        miss = argmax([sum(embarkedprime .== "C"), sum(embarkedprime .== "Q"), sum(embarkedprime .== "S")])
        embarked[ismissing.(embarked)] .= "M"
        X[8, embarked .== "M"] .= miss
        X[8, embarked .== "C"] .= 1
        X[8, embarked .== "Q"] .= 2
        X[8, embarked .== "S"] .= 3

        if normalize
            normalize!(X)
            X[1, :] .= 1
        end

        return new(passengers, ids, X)
    end

    """
        Dataset(data::Matrix{<:Real})

    construct Dataset object when the data are already prepared in matrix `data`

    skips normalization and filling of missing values
    """
    function Dataset(data::Matrix{<:Real})
        _, n = size(data)
        return new(n, collect(1:n), data)
    end
end

"""
    normalize!(X::Matrix{<:Real})

normalize data in `X` in place by subtracting the mean and dividing by standard deviation of the samples
"""
function normalize!(X::Matrix{<:Real})
    X .-= mean(X; dims=2)
    X ./= std(X; dims=2)
end

"""
    (dataset::Dataset)() = dataset.X

functor to return data matrix from `dataset` without the need to access the field directly
"""
(dataset::Dataset)() = dataset.X

"""
    Labels

structure holding inner representation of labels assigned to given dataset (it can be both training or testing dataset)

# Fields
- `ids`: id numbers assigned to samples; needed for output in Kaggle specified format
- `survived`: labels of given samples (-1 for negative class, 1 for positive class)
"""
struct Labels
    ids::Vector{<:Integer}
    survived::Vector{<:Integer}

    """
        Labels(data::DataFrame)

    extract correct labels from given `data` DataFrame
    """
    function Labels(data::DataFrame)
        ids = data[!, :PassengerId]
        survived = data[!, :Survived]
        labels = ones(Integer, length(survived))
        labels[survived .== 0] .= -1
        return new(ids, labels)
    end

    """
        Labels(data::DataFrame, survived::Vector{<:Integer})

    extract `ids` from DataFrame and assign them classes from `survived` vector
    """
    function Labels(data::Dataset, survived::Vector{<:Integer})
        return new(data.ids, survived)
    end

    """
        Labels(survived::Vector{<:Integer})

    wrap classes from `survived` vector into Labels object

    `ids` are set to increasing sequence starting from 1
    """
    function Labels(survived::Vector{<:Integer})
        n = length(survived)
        return new(collect(1:n), survived)
    end
end

"""
    (labels::Labels)() = labels.survived

functor to return labels vector from `Labels` without the need to access the field directly
"""
(labels::Labels)() = labels.survived

"""
    classificationerror(correct::Labels, prediction::Labels)

compute ratio of misclassified data samples from `prediction` when given `correct` labels
"""
function classificationerror(correct::Labels, prediction::Labels)
    return sum(correct.survived .!= prediction.survived) / length(correct.survived)
end

"""
    out(labels::Labels)

transform input representation of labels into required one by Kaggle (i.e. transofrm -1 to 0 label)
"""
function out(labels::Labels)
    survived = copy(labels.survived)
    survived[survived .== -1] .= 0
    return labels.ids, survived
end


text:
module Utilities

using DynamicPPL: AbstractSampler, Sampler
using Distributions, Bijectors
using StatsFuns, SpecialFunctions
using MCMCChains: Chains, setinfo
import Distributions: sample

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

include("robustinit.jl")
include("helper.jl")

end # module


text:
"""
    struct Mesh{T}

* `x`: X-coordinate mesh in meters
* `y`: Y-coordinate mesh in meters
* `z`: Z-coordinate mesh in meters
* `connections`: defines which points are connected to eah other

"""

struct Mesh{T}
    x::Vector{T}
    y::Vector{T}
    z::Vector{T}
    connections::Vector{Vector{Int64}}
end

function mesh(p::AbstractSurfacePrimitive{T}, n_arc::Int64)::Mesh{T} where {T}
    vs = vertices(p, n_arc)
    x, y, z = broadcast(i -> getindex.(vs, i), (1,2,3))
    c = connections(p, n_arc)
    Mesh{T}(x,y,z,c)
end

function mesh(p::AbstractSurfacePrimitive{T}, n_arc::Int64, n_vert_lines::Int64)::Mesh{T} where {T}
    vs = vertices(p, n_arc)
    x, y, z = broadcast(i -> getindex.(vs, i), (1,2,3))
    c = connections(p, n_arc, n_vert_lines)
    Mesh{T}(x,y,z,c)
end

@recipe function f(m::Mesh{T}) where {T}
    seriestype := :mesh3d
    xguide --> "x"
    yguide --> "y"
    zguide --> "z"
    unitformat --> :slash
    label --> ""
    if occursin("GRBackend", string(typeof(plotattributes[:plot_object].backend)))
        linewidth --> 2
    else
        linewidth --> 0.2
    end
    linecolor --> :black
    seriescolor --> 1
    fillalpha --> 0.5
    connections := m.connections
    m.x*internal_length_unit, m.y*internal_length_unit, m.z*internal_length_unit
end


text:
# This solution is quite messy - should revisit sometime

function single_anagram(subject, candidate)
    lc_sub = lowercase(subject)
    lc_cand = lowercase(candidate)
    sub_sort = sort(collect(lc_sub))
    cand_sort = sort(collect(lc_cand))

    return (sub_sort == cand_sort) & (lc_sub != lc_cand)

end

function detect_anagrams(subject, candidates)

    output =  [cand for cand in candidates if single_anagram(subject, lowercase(cand))]

end

max_stars_repo_name: jfcliang/jl-exercism-solutions

max_issues_repo_name: jfcliang/jl-exercism-solutions

max_forks_repo_name: jfcliang/jl-exercism-solutions


text:
module Constants

"""
Gravitational constant in Msun kpc (km / s)^2
"""
const G = 4.30091727e-06

"""
Critical density in cosmologically-invariant units of Msun h^2 / kpc3
"""
const rho_crit_h2 = 2.774848e+02

"""
radians per arcsecond
"""
const rad_per_arcsec = 4.84813681109536e-06

end


text:
## NullDict
"""
$(TYPEDEF)

An immutable associative type that signifies an empty dictionary and does not
allocate any memory.
"""
struct NullDict{K, V} <: AbstractDict{K, V}
end
Base.haskey(::NullDict, k) = false
Base.length(::NullDict) = 0
Base.iterate(::NullDict) = nothing


text:
include("complex_roots.jl");
include("spectrum.jl");
R = 2.5; rectangle = [-R,-0.001,-2*R,2*R];
maxevals = 10^4;

β = 2.5
U(x)=1/β*((1+x^2)^(β/2) - 1);
dUdx(x)=x*(1 + x^2)^(β/2-1);

(ψ,dψdz)=construct_ψ(U,dUdx);

#using QZ40
integrand_plus(z::Complex) = - dψdz(z)/(1-ψ(z));
integrand_minus(z::Complex) = dψdz(z)/(1+ψ(z));
roots_plus = qz40(integrand_plus,rectangle,true,maxevals)
roots_minus = qz40(integrand_minus,rectangle,true,maxevals)

using Plots;
using LaTeXStrings;

plot([0], [0], seriestype=:scatter,c=:white,markershape= :diamond,label="")
plot!(real(roots_minus),imag(roots_minus),seriestype=:scatter,c=:black,label=L"L^-")
plot!(real(roots_plus),imag(roots_plus),seriestype=:scatter,c=:white,markershape= :diamond,label=L"L^+")
epsilon = 0.1
perturbation_minus = zeros(Complex, size(roots_minus))
for i in 1:length(roots_minus)
    perturbation_minus[i] = epsilon * refresh_perturbation_symmetric(U,dUdx, ψ, roots_minus[i], -1)
end
perturbation_plus = zeros(Complex, size(roots_plus))
for i in 1:length(roots_plus)
    perturbation_plus[i] = epsilon * refresh_perturbation_symmetric(U,dUdx, ψ, roots_plus[i], +1)
end

quiver!(real(roots_minus),imag(roots_minus), quiver=(real(perturbation_minus),imag(perturbation_minus)),c=:black)

quiver!(real(roots_plus),imag(roots_plus), quiver=(real(perturbation_plus),imag(perturbation_plus)),c=:black,linestyle=:dot)
savefig("perturbation-beta-$(β).pdf")

max_stars_repo_name: jbierkens/spectral

max_issues_repo_name: jbierkens/spectral

max_forks_repo_name: jbierkens/spectral


text:
using Spectra
using Base.Test

# Dummy data
x = collect(0.:1.:100.)
y = collect(0.:1.:100.)
x2 = collect(-100.:0.1:100.)
y_poly4 = 1.0+ x2 + 0.25*x2.^2 + 0.0001.*x2.^3 + 1.e-6.*x2.^4
roi = [0.0 100.0]

# testing baselines
y_fit_gcvspl,bas_fit_gcvspl = baseline(x[:],y[:],roi,"gcvspline",p=1.0)
@test_approx_eq_eps(y_fit_gcvspl,(y[:]-y[:]),1e-10)
@test_approx_eq_eps(bas_fit_gcvspl,y[:],1e-10)

y_fit_dierckx,bas_fit_dierckx = baseline(x[:],y[:],roi,"Dspline",p=1.0)
@test_approx_eq_eps(y_fit_dierckx,(y[:]-y[:]),1e-10)
@test_approx_eq_eps(bas_fit_dierckx,y[:],1e-10)

y_fit_linear,bas_fit_linear = baseline(x[:],y[:],roi,"poly",p=2.0)
@test_approx_eq_eps(y_fit_linear,(y[:]-y[:]),1e-10)
@test_approx_eq_eps(bas_fit_linear,y[:],1e-10)

y_fit_poly,bas_fit_poly = baseline(x2[:],y_poly4[:],roi,"poly",p=4.0)
@test_approx_eq_eps(y_fit_poly,(y_poly4[:]-y_poly4[:]),1e-3)
@test_approx_eq_eps(bas_fit_poly,y_poly4[:],1e-3)

#Tolerances for ML algo are high, it is just to test that they are running well...
y_fit_kr,bas_fit_kr = baseline(x[:],y[:],roi,"KRregression")
@test_approx_eq_eps(y_fit_kr,(y[:]-y[:]),1e-1)
@test_approx_eq_eps(bas_fit_kr,y[:],1e-1)

y_fit_svm,bas_fit_svm = baseline(x[:],y[:],roi,"SVMregression")
@test_approx_eq_eps(y_fit_svm,(y[:]-y[:]),10.)
@test_approx_eq_eps(bas_fit_svm,y[:],10.)

# testing long correction function
x_for_long = [20.,21.,22.,23.,24.,25.]
y_for_long = [1.0,1.0,1.0,1.0,1.0,1.0]
h_plank = 6.626070040e-34 #Plank constant
k_bolt = 1.38064852e-23  #Boltzman constant
c_light = 299792458. #speed of light in vaccum
nu0 = 1.0./(514.532).*1e9 #laser wavenumber at 514.532
nu = 100.0.*x_for_long # cm-1 to m-1
T = 23.0+273.15 # the temperature in K

x_long,long_res,eselong = tlcorrection([x_for_long y_for_long],23.0,514.532) # using the function
t0 = nu0.^3.*nu./((nu0-nu).^4)
t1= 1 - exp(-h_plank.*c_light.*nu./(k_bolt.*T)) # c in m/s  : t1 dimensionless
long_calc= y_for_long.*t0.*t1 # pour les y
long_calc = long_calc./trapz(x_for_long,long_calc) # area normalisation

@test_approx_eq(long_res,long_calc)
@test_approx_eq(x_for_long,x_long)

#figure()
#plot(x2,y_poly4,"o")
#plot(x2,bas_fit_poly,"-")


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

function cut!(cutmask::AbstractArray{Bool,N},
              simdev::AbstractArray{T,N}, TIdev::AbstractArray{T,N},
              tileind::CartesianIndex{N}, pasted::Set{CartesianIndex{N}},
              geoconfig::NamedTuple) where {N,T<:Real}
  tilesize = geoconfig.tilesize
  ovlsize  = geoconfig.ovlsize
  spacing  = geoconfig.spacing

  cutmask .= false
  for d=1:N
    # Cartesian index of previous and next tiles along dimension
    prev = CartesianIndex(ntuple(i -> i == d ? (tileind[d]-1) : tileind[i], N))
    next = CartesianIndex(ntuple(i -> i == d ? (tileind[d]+1) : tileind[i], N))

    # compute mask with previous tile
    if ovlsize[d] > 1 && prev ∈ pasted
      oslice = ntuple(i -> i == d ? (1:ovlsize[d]) : (1:tilesize[i]), N)
      inds = CartesianIndices(oslice)
      A = view(simdev, inds); B = view(TIdev, inds)
      cutmask[inds] .|= boykov_kolmogorov_cut(A, B, d)
    end

    # compute mask with next tile
    if ovlsize[d] > 1 && next ∈ pasted
      oslice = ntuple(i -> i == d ? (spacing[d]+1:tilesize[d]) : (1:tilesize[i]), N)
      inds = CartesianIndices(oslice)
      A = view(simdev, inds); B = view(TIdev, inds)
      cutmask[inds] .|= reverse(boykov_kolmogorov_cut(reverse(A, dims=d), reverse(B, dims=d), d), dims=d)
    end
  end
end

function boykov_kolmogorov_cut(A::AbstractArray, B::AbstractArray, dim::Integer)
  # permute dimensions so that the algorithm is
  # the same for cuts in x, y and z directions
  if dim == 1
    A = permutedims(A, [1,2,3])
    B = permutedims(B, [1,2,3])
  elseif dim == 2
    A = permutedims(A, [2,1,3])
    B = permutedims(B, [2,1,3])
  elseif dim == 3
    A = permutedims(A, [3,2,1])
    B = permutedims(B, [3,2,1])
  end

  E = abs.(A - B)

  mx, my, mz = size(E)
  nvox = mx*my*mz

  # compute gradients
  ∇xₐ = similar(E); ∇yₐ = similar(E); ∇zₐ = similar(E)
  ∇xᵦ = similar(E); ∇yᵦ = similar(E); ∇zᵦ = similar(E)
  for k=1:mz, j=1:my, i=1:mx-1
    ∇xₐ[i,j,k] = A[i+1,j,k] - A[i,j,k]
    ∇xᵦ[i,j,k] = B[i+1,j,k] - B[i,j,k]
  end
  for k=1:mz, j=1:my-1, i=1:mx
    ∇yₐ[i,j,k] = A[i,j+1,k] - A[i,j,k]
    ∇yᵦ[i,j,k] = B[i,j+1,k] - B[i,j,k]
  end
  for k=1:mz-1, j=1:my, i=1:mx
    ∇zₐ[i,j,k] = A[i,j,k+1] - A[i,j,k]
    ∇zᵦ[i,j,k] = B[i,j,k+1] - B[i,j,k]
  end
  if mx > 1
    ∇xₐ[mx,:,:] = ∇xₐ[mx-1,:,:]
    ∇xᵦ[mx,:,:] = ∇xᵦ[mx-1,:,:]
  end
  if my > 1
    ∇yₐ[:,my,:] = ∇yₐ[:,my-1,:]
    ∇yᵦ[:,my,:] = ∇yᵦ[:,my-1,:]
  end
  if mz > 1
    ∇zₐ[:,:,mz] = ∇zₐ[:,:,mz-1]
    ∇zᵦ[:,:,mz] = ∇zᵦ[:,:,mz-1]
  end
  map!(abs, ∇xₐ, ∇xₐ); map!(abs, ∇yₐ, ∇yₐ); map!(abs, ∇zₐ, ∇zₐ)
  map!(abs, ∇xᵦ, ∇xᵦ); map!(abs, ∇yᵦ, ∇yᵦ); map!(abs, ∇zᵦ, ∇zᵦ)

  # add source and sink terminals
  s = nvox + 1; t = nvox + 2

  # construct graph and capacity matrix
  G = DiGraph(nvox+2)
  C = spzeros(nvox+2, nvox+2)
  for k=1:mz, j=1:my, i=1:mx-1
    c = cart2lin((mx,my,mz), CartesianIndex(  i, j, k))
    d = cart2lin((mx,my,mz), CartesianIndex(i+1, j, k))
    add_edge!(G, c, d)
    add_edge!(G, d, c)
    C[c,d] = C[d,c] = (E[c] + E[d]) / (∇xₐ[c] + ∇xₐ[d] + ∇xᵦ[c] + ∇xᵦ[d])
  end
  for k=1:mz, j=1:my-1, i=1:mx
    c = cart2lin((mx,my,mz), CartesianIndex(i,   j, k))
    r = cart2lin((mx,my,mz), CartesianIndex(i, j+1, k))
    add_edge!(G, c, r)
    add_edge!(G, r, c)
    C[c,r] = C[r,c] = (E[c] + E[r]) / (∇yₐ[c] + ∇yₐ[r] + ∇yᵦ[c] + ∇yᵦ[r])
  end
  for k=1:mz-1, j=1:my, i=1:mx
    c = cart2lin((mx,my,mz), CartesianIndex(i, j,   k))
    o = cart2lin((mx,my,mz), CartesianIndex(i, j, k+1))
    add_edge!(G, c, o)
    add_edge!(G, o, c)
    C[c,o] = C[o,c] = (E[c] + E[o]) / (∇zₐ[c] + ∇zₐ[o] + ∇zᵦ[c] + ∇zᵦ[o])
  end
  for k=1:mz, j=1:my
    u = cart2lin((mx,my,mz), CartesianIndex( 1, j, k))
    v = cart2lin((mx,my,mz), CartesianIndex(mx, j, k))
    add_edge!(G, s, u)
    add_edge!(G, v, t)
    C[s,u] = C[v,t] = Inf
  end

  # Boykov-Kolmogorov max-flow/min-cut
  _, __, labels = maximum_flow(G, s, t, C, algorithm=BoykovKolmogorovAlgorithm())

  # remove source and sink terminals
  labels = labels[1:end-2]

  # cut mask
  M = falses(size(E))
  M[labels .== 1] .= true
  M[labels .== 0] .= true

  # permute back to original shape
  if dim == 1
    M = permutedims(M, [1,2,3])
  elseif dim == 2
    M = permutedims(M, [2,1,3])
  elseif dim == 3
    M = permutedims(M, [3,2,1])
  end

  M
end


text:
module InferenceAlgorithmTest

using Test
using ForneyLab

@testset "InferenceAlgorithm" begin
    algo = InferenceAlgorithm()
    @test algo.posterior_factorization == currentPosteriorFactorization()
    @test currentInferenceAlgorithm() === algo
end

end # module


text:
f(x) = x + 100

macro g()         # According to Julia docs, ...
  :(f(x) + 5)     # ... f is global, x is local, right?
end               # if so, f should refer to the f above?

(function main()
  local x = 3
  f(x) = x - 100  # f in the call environment subtracts 100
#  println(@g())   # So why does this do -92?
end)()

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
module renderables

import ..measure: Measure
import ..segment: Segment
import ..console: console_width
import ..style: get_style_codes, MarkupStyle, apply_style
import Term: split_lines, reshape_text, fillin, join_lines

export AbstractRenderable, Renderable, RenderableText

# ------------------------------- abstract type ------------------------------ #

"""
    AbstractRenderable
"""
abstract type AbstractRenderable end

Measure(renderable::AbstractRenderable) = renderable.measure

"""
    Base.string(r::AbstractRenderable)::String

Creates a string representation of a renderable
"""
function Base.string(r::AbstractRenderable)::String
    lines = [seg.text for seg in r.segments]
    return join(lines, "\n")
end

function Base.show(io::IO, renderable::AbstractRenderable)
    if io isa IOContext{Base.TTY}
        print(
            io,
            "$(typeof(renderable)) <: AbstractRenderable \e[2msize: $(renderable.measure)\e[0m"
        )
    else
        map(s -> println(io, s), renderable.segments)
    end
end

# ------------------------- generic renderable object ------------------------ #

"""
    Renderable

Generic `Renderable` object.
"""

struct Renderable <: AbstractRenderable
    segments::Vector{Segment}
    measure::Measure
end

"""
    Renderable(
        str::String; width::Union{Nothing,Int} = nothing
    )

Convenience method to construct a RenderableText
"""
function Renderable(
    str::AbstractString; width::Union{Nothing,Int} = nothing
)
    return RenderableText(str; width = width)
end

Renderable(ren::AbstractRenderable) = ren
Renderable() = Renderable(Vector{Segment}[], Measure(0, 0))

# ---------------------------------------------------------------------------- #
#                                TEXT RENDERABLE                               #
# ---------------------------------------------------------------------------- #

"""
    RenderableText

`Renderable` represnting a text.

See also [`Renderable`](@ref), [`TextBox`](@ref)
"""

struct RenderableText <: AbstractRenderable
    segments::Vector
    measure::Measure
    style::Union{Nothing, String}
end

"""
    RenderableText(text::String; width::Union{Nothing, Int, Symbol}=nothing)

Construct a `RenderableText` out of a string.

If a `width` is passed the text is resized to match the width.
"""
function RenderableText(text::AbstractString; style::Union{Nothing, String}=nothing, width::Union{Nothing,Int} = nothing)
    text = apply_style(text)

    # reshape text
    if !isnothing(width)
        width = min(console_width(stdout)-1, width)
        text = reshape_text(text, width)
    end
    text = fillin(text)

    # create renderable
    if isnothing(style)
        segments = Segment.(split_lines(text))
    else
        style_init, style_finish = get_style_codes(MarkupStyle(style))
        segments = map(ln -> Segment(style_init * ln * style_finish), split_lines(text))
    end

    return RenderableText(segments, Measure(segments), style)
end

"""
    RenderableText(rt::RenderableText; width::Union{Nothing,Int} = nothing)

Construct a RenderableText by possibly re-shaping a RenderableText
"""
function RenderableText(rt::RenderableText; style::Union{Nothing, String}=nothing, width::Union{Nothing,Int} = nothing)
    if rt.style == style && rt.measure.w == width
        return rt
    else
        text = join_lines([seg.text for seg in rt.segments])
        return RenderableText(text; style=style, width=width)
    end
end

"""
    RenderableText(text::Vector; width::Union{Nothing,Int} = nothing)

Construct a RenderableText out a vector of objects.
"""
function RenderableText(text::Vector; style::Union{Nothing, String}=nothing, width::Union{Nothing,Int} = nothing)
    return RenderableText(join(string(text), "\n"); style=style, width = width)
end

# -------------------------------- union type -------------------------------- #
RenderablesUnion = Union{AbstractString,AbstractRenderable}

end


text:
struct FreshVarHeuristic end

struct HeuristicScore
    # represents utilization factor for symbolic lower bounds in the layers
    lowers::Vector{Float64}
    # represents utilization factor for symbolic upper bounds in the layers
    uppers::Vector{Float64}
end

function backward_linear(solver::FreshVarHeuristic, L::Layer, input::HeuristicScore)
    W⁺ = max.(0., L.weights')
    W⁻ = abs.(min.(0., L.weights'))
    lowers = W⁺ * input.lowers + W⁻ * input.uppers
    uppers = W⁺ * input.uppers + W⁻ * input.lowers
    return HeuristicScore(lowers, uppers)
end

function backward_act(solver::FreshVarHeuristic, L::Layer{ReLU}, input, lbs, ubs)
    λ_l = relaxed_relu_gradient_lower.(lbs, ubs)
    λ_u = relaxed_relu_gradient.(lbs, ubs)

    lowers = λ_l .* input.lowers
    uppers = λ_u .* input.uppers
    return HeuristicScore(lowers, uppers)
end

function backward_act(solver::FreshVarHeuristic, L::Layer{Id}, input, lbs, ubs)
    return input
end

"""
Calculates factors for usage of node's lower and upper bound later on in the network.

input (HeuristicScore) - lower and upper bound usage factors of the output nodes

pre_activation (bool) - whether to return factors for the pre-activation values of a neuron
post_activation (bool) - whether to return factors for the post-activation values of a neuron
"""
function backward_network(solver::FreshVarHeuristic, net, lbs, ubs, input; pre_activation=true, post_activation=true)
    scores = []

    Z = backward_linear(solver, net.layers[end], input)
    post_activation && push!(scores, Z)

    for i in reverse(2:length(net.layers)-1)
        layer = net.layers[i]
        Ẑ = backward_act(solver, layer, Z, lbs[i], ubs[i])
        Z = backward_linear(solver, layer, Ẑ)

        pre_activation && push!(scores, Ẑ)
        post_activation && push!(scores, Z)
    end

    # we iterated backwards through the NN, but want scores forward
    return reverse(scores)
end

function cancellation_score(l, u)
    M = max(l, u)
    m = min(l, u)

    return M == 0. ? 0. : m / M
end

# Methods for calculating amount of overapproximation of nodes

abstract type OverapproxCalculator end

struct AreaCalculator end
struct RangeCalculator end
struct OneSideBound end
struct SplitAreaCalculator end

function calc_overapprox(method::AreaCalculator, l, u)
    l >= 0 && return 0.
    u <= 0 && return 0.
    upper_area = 0.5*(-l)*u
    if u > -l
        lower_area = 0.5*l*l
    else
        lower_area = 0.5*u*u
    end

    return upper_area + lower_area
end

function calc_overapprox(method::SplitAreaCalculator, l, u)
    l >= 0 && return 0., 0.
    u <= 0 && return 0., 0.
    upper_area = 0.5*(-l)*u
    if u > -l
        lower_area = 0.5*l*l
    else
        lower_area = 0.5*u*u
    end

    return upper_area, lower_area
end

function calc_overapprox(method::RangeCalculator, l, u)
    l >= 0 && return 0
    u <= 0 && return 0
    return u - l
end

function calc_overapprox(method::OneSideBound, l, u)
    l >= 0 && return 0
    u <= 0 && return 0
    return min(-l, u)
end

######## Local Layer introduction Heuristics

function fresh_var_earliest(max_vars, current_n_vars, los, his, var_frac)
    n_node = length(los)
    n_vars = min(max_vars - current_n_vars, floor(Int, var_frac * n_node))

    return fresh_var_first(los, his, n_vars)
end

function fresh_var_range_non_zero(max_vars, current_n_vars, los, his, var_frac)
    n_node = length(los)
    n_fixed_zero = sum(his .<= 0)
    n_node = n_node - n_fixed_zero  # fixed zero nodes are basically non-existent in the nn
    n_vars = min(max_vars - current_n_vars, ceil(Int, var_frac * n_node))
    return fresh_var_largest_range(los, his, n_vars)
end

function fresh_var_n_vars_non_zero(max_vars, current_n_vars, los, his, var_frac)
    n_node = length(los)
    n_fixed_zero = sum(his .<= 0)
    n_node = n_node - n_fixed_zero  # fixed zero nodes are basically non-existent in the nn

    λ = n_node / (current_n_vars + n_node)
    n_vars = min(max_vars - current_n_vars, ceil(Int, var_frac * λ * n_node))

    ranges = his - los

    p = sortperm(-ranges)
    p = p[(los[p] .< 0) .& (his[p] .> 0)] # only crossing relus
    return p[1:min(length(p), n_vars)]
end


text:
"""
For localized modes, the L₁-norm is more reliable, but, if the result is extended,
the fit procedure is unstable, and it more meaninfull to apply fit with L₂-norm instead.
"""
function get_single_ξ_and_R1(x, y)
    x_use, y_use = select_points(x, y)
    A_L₁, ξ_L₁, R1, y_fit_L₁ = get_coeffs_L1(x_use, y_use)

    if is_R1_localized(R1)
        return ξ_L₁, R1
    else
        A_L₂, ξ_L₂, R2, y_fit_L₂ = get_coeffs_L2(x_use, y_use)
        return ξ_L₂, R2
    end
end
function get_coeffs_L1(x_use, y_use)
    A_L₁, ξ_L₁ = get_exp_fit_coeffs_with_L1(x_use, y_use)
    y_fit_L₁ = model_exp(x_use, [A_L₁, ξ_L₁])
    R1 = compute_R1(y_use, y_fit_L₁)
    return A_L₁, ξ_L₁, R1, y_fit_L₁
end
function get_coeffs_L2(x_use, y_use)
    A_L₂, ξ_L₂ = get_exp_fit_coeffs_with_L2(x_use, y_use)
    y_fit_L₂ = model_exp(x_use, [A_L₂, ξ_L₂])
    R2 = compute_R2(y_use, y_fit_L₂)
    return A_L₂, ξ_L₂, R2, y_fit_L₂
end
function is_R1_localized(R1)
    # is an arbitrary decision to say that a mode is  "localized" if R1 > "somthing"
    if R1 > R1_threshold
        return true
    else
        return false
    end
end

"""
    x_use, y_use = select_points(DCM::Vector, ψ²::Vector)

computes how many points since the center of mass are important based upon the `probability_threshold` (constant defined as 0.9999)
"""
function select_points(DCM::Vector, ψ²::Vector)
    if length(ψ²) < 10
        error("Very few data. Need at least 10 points")
    end

    cumulative_probability = cumsum(ψ²)
    threshold = sum(cumulative_probability .< probability_threshold)

    x_use = DCM[1:threshold]
    y_use = ψ²[1:threshold]

    return x_use, y_use
end

######### Now, we actually compute the coefficients
"""
    A,ξ = get_exp_fit_coeffs_with_L1(x, y)

Assumes that (x,y) have a exponential decay of form `y(x) = A*exp(-x/ξ)`
"""
function get_exp_fit_coeffs_with_L1(x, y)
    fit_result =
        optimize(xx -> model_exp_minimize(xx, x, y), rand(2), Options(iterations = 30))
    A = minimizer(fit_result)[1]
    ξ = minimizer(fit_result)[2]
    return A, ξ
end
model_exp_minimize(x0, x, y) = sum(abs.(x0[1] * exp.(-x ./ x0[2]) .- y))
@. model_exp(x, p) = p[1] * exp.(-x ./ p[2])

"""
    A,ξ = get_exp_fit_coeffs_with_L2(x, y)

Assumes that (x,y) have a exponential decay in linear form, that is, `log(y(x)) = log(A) - x/ξ`
"""
function get_exp_fit_coeffs_with_L2(x, y)
    fit = curve_fit(model_linear, x, log.(y), [maximum(log.(y)), rand()])
    A = exp(coef(fit)[1])
    ξ = coef(fit)[2]
    return A, ξ
end
@. model_linear(x, p) = p[1] - x / p[2]

######### Quality of the Fit

"""
    compute_R2(y, y_predict)

returns the R2 (adjusted) for a data set "y", and its modeled values "y_predict"

both inputs should have the same size

example:
actual = [2, 4, 6, 7]
predict = [2.601,3.83,5.059,7.517]

compute_R2(actual, predict)
0.8430978644067797

Ref:
Coefficient of Determination, R-squared, ASK Academic Skills Kit, Newcastle University
Coefficient of Determination (R-Squared), Matlab Documentation
"""
compute_R2(y, y_predict) = generic_R(log.(y), log.(y_predict), SSE_R2, SST_R2)
compute_R1(y, y_predict) = generic_R(y, y_predict, SSE_R1, SST_R1)

function generic_R(y, y_predict, SSE_fnc::Function, SST_fnc::Function)
    if size(y) ≠ size(y_predict)
        @error("Inputs does not have similar dimension")
    end

    average_y = sum(y) / length(y)

    SSE = SSE_fnc(y, y_predict)# the sum of squared error
    SST = SST_fnc(y, average_y) # sum of squared total

    n = length(y)
    p = 2 # number of regression coefficients, for linear, p=2
    R2_adjusted = 1 - ((n - 1) / (n - p)) * (SSE / SST)

    if R2_adjusted < 0
        R2_adjusted = 1e-16
    end
    return R2_adjusted
end

SSE_R1(y, y_predict) = sum(abs(y[i] - y_predict[i]) for i in eachindex(y)) # the sum of squared error
SST_R1(y, average_y) = sum(abs(y[i] - average_y) for i in eachindex(y)) # sum of squared total

SSE_R2(y, y_predict) = sum((y[i] - y_predict[i])^2 for i in eachindex(y)) # the sum of squared error
SST_R2(y, average_y) = sum((y[i] - average_y)^2 for i in eachindex(y)) # sum of squared total


text:
module XGBoost_

#> export the new models you're going to define (and nothing else):
export XGBoostRegressor, XGBoostClassifier

#> for all Supervised models:
import MLJ
import MLJ: CanWeightTarget, CanRankFeatures
import MLJ: Nominal, Numeric, NA, Probababilistic, Multivariate,  Multiclass

#> for all classifiers:
using CategoricalArrays

#> import package:
import XGBoost

mutable struct XGBoostRegressor{Integer} <:MLJ.Supervised{Integer}  #check the parametric type here..
    num_round::Integer
    booster::String
    silent::Union{Int,Bool}
    disable_default_eval_metric::Real
    eta::Real
    gamma::Real
    max_depth::Real
    min_child_weight::Real
    max_delta_step::Real
    subsample::Real
    colsample_bytree::Real
    colsample_bylevel::Real
    lambda::Real
    alpha::Real
    tree_method::String
    sketch_eps::Real
    scale_pos_weight::Real
    updater::String
    refresh_leaf::Union{Int,Bool}
    process_type::String
    grow_policy::String
    max_leaves::Int
    max_bin::Int
    predictor::String
    sample_type::String
    normalize_type::String
    rate_drop::Real
    one_drop
    skip_drop::Real
    feature_selector::String
    top_k::Real
    tweedie_variance_power::Real
    objective
    base_score::Real
    eval_metric
    seed::Integer
    watchlist
    num_class::Integer
end

"""
# constructor:
A full list of the kwargs accepted, and their value ranges, consult
https://xgboost.readthedocs.io/en/latest/parameter.html.
The only required kwarg is num_round.
"""
function XGBoostRegressor(
    ;num_round=1
    ,booster="gbtree"
    ,silent=0  #> might be redundant due to verbosity
    ,disable_default_eval_metric=0
    ,eta=0.3
    ,gamma=0
    ,max_depth=6
    ,min_child_weight=1
    ,max_delta_step=0
    ,subsample=1
    ,colsample_bytree=1
    ,colsample_bylevel=1
    ,lambda=1
    ,alpha=0
    ,tree_method="auto"
    ,sketch_eps=0.03
    ,scale_pos_weight=1
    ,updater="grow_colmaker"
    ,refresh_leaf=1
    ,process_type="default"
    ,grow_policy="depthwise"
    ,max_leaves=0
    ,max_bin=256
    ,predictor="cpu_predictor" #> gpu version not currently working with Julia, maybe remove completely?
    ,sample_type="uniform"
    ,normalize_type="tree"
    ,rate_drop=0.0
    ,one_drop=0
    ,skip_drop=0.0
    ,feature_selector="cyclic"
    ,top_k=0
    ,tweedie_variance_power=1.5
    ,objective="reg:linear"
    ,base_score=0.5
    ,eval_metric="rmse"
    ,seed=0
    ,watchlist=[]
    ,num_class=1)

    model = XGBoostRegressor{Integer}(
    num_round
    ,booster
    ,silent  #> might be redundant due to verbosity
    ,disable_default_eval_metric
    ,eta
    ,gamma
    ,max_depth
    ,min_child_weight
    ,max_delta_step
    ,subsample
    ,colsample_bytree
    ,colsample_bylevel
    ,lambda
    ,alpha
    ,tree_method
    ,sketch_eps
    ,scale_pos_weight
    ,updater
    ,refresh_leaf
    ,process_type
    ,grow_policy
    ,max_leaves
    ,max_bin
    ,predictor #> gpu version not currently working with Julia, maybe remove completely?
    ,sample_type
    ,normalize_type
    ,rate_drop
    ,one_drop
    ,skip_drop
    ,feature_selector
    ,top_k
    ,tweedie_variance_power
    ,objective
    ,base_score
    ,eval_metric
    ,seed
    ,watchlist
    ,num_class)

     message = MLJ.clean!(model)           #> future proof by including these
     isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

mutable struct XGBoostClassifier{Integer} <:MLJ.Supervised{Integer}  #check the parametric type here..
    num_round::Integer
    booster::String
    silent::Union{Int,Bool}
    disable_default_eval_metric::Real
    eta::Real
    gamma::Real
    max_depth::Real
    min_child_weight::Real
    max_delta_step::Real
    subsample::Real
    colsample_bytree::Real
    colsample_bylevel::Real
    lambda::Real
    alpha::Real
    tree_method::String
    sketch_eps::Real
    scale_pos_weight::Real
    updater::String
    refresh_leaf::Union{Int,Bool}
    process_type::String
    grow_policy::String
    max_leaves::Int
    max_bin::Int
    predictor::String
    sample_type::String
    normalize_type::String
    rate_drop::Real
    one_drop
    skip_drop::Real
    feature_selector::String
    top_k::Real
    tweedie_variance_power::Real
    objective
    base_score::Real
    eval_metric
    seed::Integer
    watchlist
    num_class::Integer
end

"""
# constructor:
A full list of the kwargs accepted, and their value ranges, consult
https://xgboost.readthedocs.io/en/latest/parameter.html.
The only required kwarg is num_round.
"""
function XGBoostClassifier(
    ;num_round=1
    ,booster="gbtree"
    ,silent=0  #> might be redundant due to verbosity
    ,disable_default_eval_metric=0
    ,eta=0.3
    ,gamma=0
    ,max_depth=6
    ,min_child_weight=1
    ,max_delta_step=0
    ,subsample=1
    ,colsample_bytree=1
    ,colsample_bylevel=1
    ,lambda=1
    ,alpha=0
    ,tree_method="auto"
    ,sketch_eps=0.03
    ,scale_pos_weight=1
    ,updater="grow_colmaker"
    ,refresh_leaf=1
    ,process_type="default"
    ,grow_policy="depthwise"
    ,max_leaves=0
    ,max_bin=256
    ,predictor="cpu_predictor" #> gpu version not currently working with Julia, maybe remove completely?
    ,sample_type="uniform"
    ,normalize_type="tree"
    ,rate_drop=0.0
    ,one_drop=0
    ,skip_drop=0.0
    ,feature_selector="cyclic"
    ,top_k=0
    ,tweedie_variance_power=1.5
    ,objective="binary:logistic"
    ,base_score=0.5
    ,eval_metric="rmse"
    ,seed=0
    ,watchlist=[]
    ,num_class=1)

    model = XGBoostClassifier{Integer}(
    num_round
    ,booster
    ,silent  #> might be redundant due to verbosity
    ,disable_default_eval_metric
    ,eta
    ,gamma
    ,max_depth
    ,min_child_weight
    ,max_delta_step
    ,subsample
    ,colsample_bytree
    ,colsample_bylevel
    ,lambda
    ,alpha
    ,tree_method
    ,sketch_eps
    ,scale_pos_weight
    ,updater
    ,refresh_leaf
    ,process_type
    ,grow_policy
    ,max_leaves
    ,max_bin
    ,predictor #> gpu version not currently working with Julia, maybe remove completely?
    ,sample_type
    ,normalize_type
    ,rate_drop
    ,one_drop
    ,skip_drop
    ,feature_selector
    ,top_k
    ,tweedie_variance_power
    ,objective
    ,base_score
    ,eval_metric
    ,seed
    ,watchlist
    ,num_class)

     message = MLJ.clean!(model)           #> future proof by including these
     isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

function MLJ.clean!(model::XGBoostRegressor)
    warning = ""
    if(model.booster=="gblinear" &&(model.updater!=("shotgun") && model.updater!=("coord_descent")))
        model.updater="shotgun"
        warning *= "updater has been changed to shotgun, the default option for booster=\"gblinear\""
    end
    if(model.objective in ["multi:softmax","multi:softprob","binary:logistic","binary:logitraw","binary:hinge"])
            warning *="\n objective function is more suited to XGBoostClassifier"
    end
    if(model.objective in ["multi:softmax","multi:softprob"])
        if model.num_class==1
            model.num_class=2
            warning *= "\n num_class has been changed to 2"
        end
        if model.eval_metric=="rmse"
            model.eval_metric="mlogloss"
            warning *= "\n eval_metric has been changed to mlogloss"
        end
    end
    return warning
end

function MLJ.clean!(model::XGBoostClassifier)
    warning = ""
    if(model.booster=="gblinear" &&(model.updater!=("shotgun") && model.updater!=("coord_descent")))
        model.updater="shotgun"
        warning *= "updater has been changed to shotgun, the default option for booster=\"gblinear\""
    end
    if(!(model.objective in ["multi:softmax","multi:softprob","binary:logistic","binary:logitraw","binary:hinge"]))
            warning *="\n objective function is more suited to XGBoostClassifier"
    end
    if(model.objective in ["multi:softmax","multi:softprob"])
        if model.num_class==1
            model.num_class=2
            warning *= "\n num_class has been changed to 2"
        end
        if model.eval_metric=="rmse"
            model.eval_metric="mlogloss"
            warning *= "\n eval_metric has been changed to mlogloss"
        end
    end
    return warning
end

#> The following optional method (the fallback does nothing, returns
#> empty warning) is called by the constructor above but also by the
#> fit methods below:

#> A required `fit` method returns `fitresult, cache, report`. (Return
#> `cache=nothing` unless you are overloading `update`)

function MLJ.fit(model::Union{XGBoostRegressor,XGBoostClassifier}
             , verbosity            #> must be here even if unsupported in pkg
             , X::Array{<:Real,2}
             , y::Vector)

    dm = XGBoost.DMatrix(X,label=y)
    fitresult = XGBoost.xgboost( dm
                               , model.num_round
                               , booster = model.booster
                               , silent = verbosity
                               , disable_default_eval_metric = model.disable_default_eval_metric
                               , eta = model.eta
                               , gamma = model.gamma
                               , max_depth = model.max_depth
                               , min_child_weight = model.min_child_weight
                               , max_delta_step = model.max_delta_step
                               , subsample = model.subsample
                               , colsample_bytree = model.colsample_bytree
                               , colsample_bylevel = model.colsample_bylevel
                               , lambda = model.lambda
                               , alpha = model.alpha
                               , tree_method = model.tree_method
                               , sketch_eps = model.sketch_eps
                               , scale_pos_weight = model.scale_pos_weight
                               , updater = model.updater
                               , refresh_leaf = model.refresh_leaf
                               , process_type = model.process_type
                               , grow_policy = model.grow_policy
                               , max_leaves = model.max_leaves
                               , max_bin = model.max_bin
                               , predictor = model.predictor
                               , sample_type = model.sample_type
                               , normalize_type = model.normalize_type
                               , rate_drop = model.rate_drop
                               , one_drop = model.one_drop
                               , skip_drop = model.skip_drop
                               , feature_selector = model.feature_selector
                               , top_k = model.top_k
                               , tweedie_variance_power = model.tweedie_variance_power
                               , objective = model.objective
                               , base_score = model.base_score
                               , eval_metric=model.eval_metric
                               , seed = model.seed
                               , watchlist=model.watchlist
                               , num_class=model.num_class)

    #> return package-specific statistics (eg, feature rankings,
    #> internal estimates of generalization error) in `report`, which
    #> should be `nothing` or a dictionary keyed on symbols.

    cache = nothing
    report = nothing

    return fitresult, cache, report

end

function MLJ.fit(model::Union{XGBoostRegressor,XGBoostClassifier}
             , verbosity            #> must be here even if unsupported in pkg
             , dm::XGBoost.DMatrix)

    fitresult = XGBoost.xgboost( dm
                               , model.num_round
                               , booster = model.booster
                               , silent = verbosity
                               , disable_default_eval_metric = model.disable_default_eval_metric
                               , eta = model.eta
                               , gamma = model.gamma
                               , max_depth = model.max_depth
                               , min_child_weight = model.min_child_weight
                               , max_delta_step = model.max_delta_step
                               , subsample = model.subsample
                               , colsample_bytree = model.colsample_bytree
                               , colsample_bylevel = model.colsample_bylevel
                               , lambda = model.lambda
                               , alpha = model.alpha
                               , tree_method = model.tree_method
                               , sketch_eps = model.sketch_eps
                               , scale_pos_weight = model.scale_pos_weight
                               , updater = model.updater
                               , refresh_leaf = model.refresh_leaf
                               , process_type = model.process_type
                               , grow_policy = model.grow_policy
                               , max_leaves = model.max_leaves
                               , max_bin = model.max_bin
                               , predictor = model.predictor
                               , sample_type = model.sample_type
                               , normalize_type = model.normalize_type
                               , rate_drop = model.rate_drop
                               , one_drop = model.one_drop
                               , skip_drop = model.skip_drop
                               , feature_selector = model.feature_selector
                               , top_k = model.top_k
                               , tweedie_variance_power = model.tweedie_variance_power
                               , objective = model.objective
                               , base_score = model.base_score
                               , eval_metric=model.eval_metric
                               , seed = model.seed
                               , watchlist=model.watchlist
                               , num_class = model.num_class)

    #> return package-specific statistics (eg, feature rankings,
    #> internal estimates of generalization error) in `report`, which
    #> should be `nothing` or a dictionary keyed on symbols.

    cache = nothing
    report = nothing

    return fitresult, cache, report

end

#Order of predict matters, see https://discourse.julialang.org/t/fun-with-kwargs-methods/15711/5

function MLJ.predict(model::Union{XGBoostRegressor,XGBoostClassifier}
        , fitresult::XGBoost.Booster
        , Xnew::Union{XGBoost.DMatrix,Array{<:Real,2}}
        ; ntree_limit::Integer)

    return XGBoost.predict(fitresult, Xnew,ntree_limit=ntree_limit)
end

function MLJ.predict(model::Union{XGBoostRegressor,XGBoostClassifier}
        , fitresult::XGBoost.Booster
        , Xnew::Union{XGBoost.DMatrix,Array{<:Real,2}})

    return XGBoost.predict(fitresult, Xnew)
end

#module end
end

using .XGBoost_
export XGBoostClassifier,XGBoostRegressor


text:
function main()
end

max_stars_repo_name: IvanJericevich/aeron-julia

max_issues_repo_name: IvanJericevich/aeron-julia

max_forks_repo_name: IvanJericevich/aeron-julia


text:
module DashBioUtils

using Glob, GZip

using HTTP, StringEncodings, FASTX
using StatsBase

include("ngl_parser.jl")
include("protein_reader.jl")
include("xyz_reader.jl")

end # module


text:
using DelimitedFiles

function ReadData(fname::AbstractString)
    return readdlm(fname, '\n', Int)
end

function CountIncreases(data::Array)
    count = 0
    for i = 1:(length(data)-1)
        count += data[i+1]>data[i]
    end

    return count
end

function CountWindowIncreases(data::Array)
    count = 0
    for i = 1:(length(data)-3)
        count += sum(data[i+1:i+3]) > sum(data[i:i+2])
    end

    return count
end

# test data
#data = [199, 200, 208, 210, 200, 207, 240, 269, 260, 263];

# real data
data = ReadData("../data/day1.txt")

# first half
println(CountIncreases(data))

# second half
println(CountWindowIncreases(data))


text:
"""
    FelsensteinFunction(tree_postorder::Vector{N}, pi_::T, rates::Vector{Float64}, data::Array{Float64,3}, n_c::Int64, blv::Vector{Float64}) where {T<:Real, N<:GeneralNode}

This function calculates the log-likelihood of an evolutiuonary model using the
Felsensteins pruning algorithm.
The function is written such that it is differentiable by Zygote 0.5.3.

Returns log-likelihood as a Real number.

* `tree_postorder` : Vector of Nodes.

* `pi` : Real number used in calculation.

* `rates` : Vector of Floats, used in calculation.

* `data` : Array of Floats, currently a placeholder.

* `n_c` : Int64, nchar value derived from tree and dataset.

* `blv` : Vector of Floats, branchlength vector derived from tree.
"""
function FelsensteinFunction(tree_postorder::Vector{N}, pi_::Array{Float64}, rates::Float64,
                             data::Array{Float64,3}, n_c::Int64, blv::Array{Float64,1})::Float64 where {N<:GeneralNode}

    mu::Float64 =  1.0 / (2.0 * prod(pi_))
    mutationArray::Vector{Array{Float64, 2}} = calc_trans.(blv, pi_[1], mu, rates)
    root_node::N = last(tree_postorder)
    rns::Array{Float64} = zeros(Float64, 1, n_c)

    @views for node in tree_postorder
        if node.nchild > 0
            res = node_loop(node, mutationArray, data)
            if !node.root
                scaler::Array{Float64} = maximum(res, dims=1)
                rns = rns .+ log.(scaler)
                res = res ./ scaler
            end #if

            data[:,:,node.num] = res
        end #if
    end # for
    return likelihood_root(data[:,:,root_node.num], pi_, rns)
end # function

@inline function likelihood_root(root::Array{Float64,2}, pi_::Array{Float64}, rns::Array{Float64,2})::Float64
    sum(log.(sum(root .* pi_, dims=1))+rns)
end

function node_loop(node::N, mutationArray::Vector{Array{Float64, 2}}, data)::Array{Float64,2} where {N<:GeneralNode, S<:Real}
    # creating a new array is necessary because Zygote can not differentiate otherwise
    out = ones(size(data[:,:,node.num]))
    @inbounds @views for child in node.children
            out = out .* BLAS.gemm('N','N',mutationArray[child.num], data[:,:,child.num])
    end
    out
end

function calc_trans(time::Float64, pi_::S, mu::Float64, r::Float64)::Array{Float64,2} where {S<:Real}

    v = exp(-r*time*mu)
    v1 = pi_ - (pi_ - 1)*v
    v2 = pi_ - pi_* v
    v3 = (pi_ - 1)*(v - 1)
    v4 = pi_*(v - 1) + 1
    mm::Array{Float64}=[v1 v3;
                        v2 v4]
    return mm
end

function calc_trans(time::Float64, pi_::Array{Float64}, mu::Float64, r::Float64)::Array{Float64,2}
    res::Array{Float64} = repeat(reshape(pi_, (1,2)), 2)
    di = Diagonal(ones(2))
    v = exp(-r*time*mu)
    res = res .*(1.0-v)
    res .+ (di .* v)
end

"""
    FelsensteinFunction(tree_postorder::Vector{N}, pi_::Array{Float64}, rates::Array{Float64},
                         U::Array{Float64,2}, D::Array{Float64}, Uinv::Array{Float64,2},
                         data::Array{Float64,4}, c_grad::Bool = true) where {N<:GeneralNode}

This function calculates the log-likelihood of an evolutiuonary model using the
Felsensteins pruning algorithm. If `c_grad` equals `true` (default) the analytic gradient
regarding the branch lengths of the tree gets computed as well.
"""
function FelsensteinFunction(tree_postorder::Vector{N}, pi_::Array{Float64}, rates::Float64,
                      U::A, D::V, Uinv::A, mu::Float64,
                     data::Array{Float64,3}, substitutionModel::Function,  c_grad::Bool = true) where {N<:GeneralNode, M<:Number, A<:AbstractArray, V<:AbstractVector}
    Nbases, Nsites, Nnodes = size(data)

    grv::Vector{Float64} = Vector{Float64}(undef, Nnodes-1)
    Down::Array{Float64,3} = similar(data)
    pre_order_partial::Array{Float64,3} = similar(data)
    ll = fels_ll(tree_postorder, data, D, U, Uinv, rates, mu, Down, pi_, substitutionModel)

    if c_grad
        grv = fels_grad(tree_postorder, data, D, U, Uinv, rates, mu,
                        Down, pi_, pre_order_partial, grv, substitutionModel)
    end

    return ll, grv, data, Down, pre_order_partial

end

function fels_grad(tree_postorder::Vector{N}, data::Array{Float64,3},
         D::V, U::A, Uinv::A, rate::Float64, mu::Float64, Down::Array{Float64,3},
         pi_::Array{Float64}, pre_order_partial::Array{Float64,3}, grv::Vector{Float64}, substitutionModel::Function)::Vector{Float64} where {N <: GeneralNode, M <: Number,
                                                                                A<:AbstractArray, V<:AbstractVector}

    root_node = last(tree_postorder)
    pre_order_partial[:, :, root_node.num] .= pi_

    @inbounds @views for node in reverse(tree_postorder)[2:end]
           mother = get_mother(node)
            pre_order_partial[:, :, node.num] .= pre_order_partial[:, :, mother.num]
            @simd for c in 1:mother.nchild
                child = mother.children[c]
                if child.num != node.num
                    pre_order_partial[:, :, node.num] .*= Down[:, :, child.num]
                end
            end

            ptg = U * ((D .* (rate * mu)) .* diagm(exp.(D .* (rate*mu*node.inc_length)))) * Uinv
            gradi = sum(pre_order_partial[:, :,  node.num] .* (ptg * data[:, :, node.num]),dims=1)
            tr = calculate_transition(substitutionModel, rate, mu, node.inc_length, U, Uinv, D, pi_)
            pre_order_partial[:, :,  node.num] .= transpose(tr) * pre_order_partial[:, :, node.num]
            grv[node.num] = sum(gradi ./ sum(pre_order_partial[:, :,  node.num] .* data[:, :,  node.num], dims=1))

            if node.nchild > 0
                pre_order_partial[:, :,  node.num] ./= maximum(pre_order_partial[:, :,  node.num], dims = 1)
            end
    end # for
    grv
end

function fels_ll(tree_postorder::Vector{N}, data::Array{Float64,3},
          D::V, U::A, Uinv::A,
         rate::Float64, mu::Float64,Down::Array{Float64,3},
         pi_::Array{Float64}, substitutionModel::Function)::Float64 where {N <: GeneralNode, M<:Number, A<:AbstractArray, V<:AbstractVector}

    ll::Float64 = 0.0
    @inbounds @views for node in tree_postorder
        if node.nchild > 0
            data[:, :, node.num] .= 1.0
            for child in node.children
                mul!(Down[:, :, child.num], calculate_transition(substitutionModel, rate, mu, child.inc_length, U, Uinv, D, pi_), data[:, :, child.num])
                data[:, :, node.num] .*= Down[:, :, child.num]
            end

            if !node.root
                ll += sum(log.(maximum(data[:, :, node.num], dims=1)))
                data[:, :, node.num] ./= maximum(data[:, :, node.num], dims=1)
                #ll += sum(log.(scaler))
            else
                ll += sum(log.(sum(data[:, :, last(tree_postorder).num] .* pi_, dims=1)))
            end
        end #if
    end # for
    ll
end

function recursive_ll(node::N, data::Array{Float64,3},
                        D::V, U::A, Uinv::A,
                    rate::Float64, mu::Float64,Down::Array{Float64,3},
                    pi_::Array{Float64}, substitutionModel::Function)::Float64 where {N <: GeneralNode, M<:Number, A<:AbstractArray, V<:AbstractVector}
    ll = 0.0
    if node.nchild > 0
        data[:, :, node.num] .= 1.0
        for child in node.children
            ll += recursive_ll(child, data, D, U, Uinv, rate, mu, Down, pi_, substitutionModel)

            tr = calculate_transition(substitutionModel, rate, mu, child.inc_length, U, Uinv, D, pi_)
            Down[:, :, child.num] .= tr * data[:, :, child.num]
            data[:, :, node.num] .*= Down[:, :, child.num]
        end

        if !node.root
            scaler = maximum(data[:, :, node.num], dims=1)
            data[:, :, node.num] ./= scaler
            ll += sum(log.(scaler))
        else
            ll += sum(log.(sum(data[:, :, node.num] .* pi_, dims=1)))
        end
    end
    return ll
end

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


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

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

# ╔═╡ 27bdb510-ad38-11eb-3854-4dc8f7ca7154
begin
	using Plots,PlutoUI
	plotlyjs()
end

# ╔═╡ 7a8e9c8e-1471-41f5-81e1-5c7d92994521
md"""# What the ϵ-δ!!!

En este cuaderno vamos a explorar el concepto de continuidad de una funcion utilizando la definicion de ϵ-δ.
"""

# ╔═╡ 30983170-71fa-4b4d-8142-39f8f60768a0
md"Primero definamos la funcion ``f(x)``"

# ╔═╡ 8ca53e3b-219c-4d87-bd32-28aa2db75cca
f(x) = x

# ╔═╡ e8cee677-00b9-4511-a56d-fff53282719e
md"Y eligamos un punto ``x_{0} \in Dom(f)`` en el cual vamos a analizar la funcion"

# ╔═╡ a0aa85d3-2c69-4617-931a-b0c060626630
x₀ = 1

# ╔═╡ 4ff601e9-085c-494d-b405-22f56889a475
md"""
Usen los siguientes controles para responder la siguientes preguntas

#### A)

Si me dan un```\epsilon > 0``, cuanto tiene que valer ``δ`` para que la imagen ``f(x)`` de todos los puntos ``x`` que estan a una distancia menor a ``δ`` de ``x_0`` (``| x - x_0| < \delta ``) esten a una distancia menor a ``\epsilon`` de ``f(x_0)`` (``|f(x) - f(x_0) | < \epsilon ``)?

#### B)

Puedo hacer ``\epsilon`` tan chico como yo quiera y encontrar un ``\delta`` que cumpla la condicion de arriba?
"""

# ╔═╡ e5d0e972-6009-4a0d-849f-6a43103e99a0
# δ₁,δ₂,ϵ₀,E=0.1,0.2,1.0,1.0
md"""
ϵ₀ = $(@bind ϵ₀ Slider(0.01:0.01:1,default=1,show_value=true))
E = $(@bind E Slider(-5:1,default=0,show_value=true))

δ₁ = $(@bind δ₁ Slider(0.01:0.01:1,default=1,show_value=true))
δ₂ = $(@bind δ₂ Slider(0.01:0.01:1,default=1,show_value=true))
"""

# ╔═╡ d99e9b14-4f51-468f-8d13-9db1fd29f5b8
md"ϵ = ``ϵ₀*10^E`` = $(ϵ = ϵ₀*10^E)  |  δ = max(δ₁,δ₂) = $(δ = max(δ₁,δ₂))"

# ╔═╡ b18fa593-cea1-446c-a56e-4b22cb2aee4c
begin
x0 = floor(x₀)
md"""
``X`` = $(@bind x_i Scrubbable(x0-5:0.1:x0,default=x0-3)):
$(@bind x_s Scrubbable(0.001:0.05:1,default=1,format=".2g")):
$(@bind x_f Scrubbable(x0:0.1:x0+5,default=x0+3))

El dominio de puntos para graficar *inicio:espaciado:final*
"""
end

# ╔═╡ 17bb2999-2bf3-4eef-94d0-7f610f67cd97
begin
	xs = x_i:x_s:x_f
	y₀ = f(x₀)
	plot(xs,f.(xs), label=false)
	scatter!([(x₀,y₀)], label=false,color=:red)
	hline!([y₀-ϵ,y₀+ϵ], label=false)
	vline!([x₀ - δ₁,x₀ + δ₂], label=false)
end

# ╔═╡ f3789190-904d-4d1b-afa2-e2e952c84b4d
md"""
**Ejemplos de funciones que pueden probar**

1) f(x) un polinomio de grado igual o mayor a 2.

2) f(x) una funcion homografica.

3) La siguiente funcion definida en partes

```math
f(x) = \begin{cases}
	  x & x \leq x_{0} \\
	  x - 1 & x \gt x _{0}
   \end{cases}
```
Pueden copiar el siguiente codigo
> f(x) =  x <= x₀ ? x : x-1

"""

# ╔═╡ Cell order:
# ╟─27bdb510-ad38-11eb-3854-4dc8f7ca7154
# ╟─7a8e9c8e-1471-41f5-81e1-5c7d92994521
# ╟─30983170-71fa-4b4d-8142-39f8f60768a0
# ╠═8ca53e3b-219c-4d87-bd32-28aa2db75cca
# ╟─e8cee677-00b9-4511-a56d-fff53282719e
# ╠═a0aa85d3-2c69-4617-931a-b0c060626630
# ╟─4ff601e9-085c-494d-b405-22f56889a475
# ╟─e5d0e972-6009-4a0d-849f-6a43103e99a0
# ╟─d99e9b14-4f51-468f-8d13-9db1fd29f5b8
# ╟─17bb2999-2bf3-4eef-94d0-7f610f67cd97
# ╟─b18fa593-cea1-446c-a56e-4b22cb2aee4c
# ╟─f3789190-904d-4d1b-afa2-e2e952c84b4d

max_stars_repo_name: akielbowicz/cbc-analisis-ingenieria

max_issues_repo_name: akielbowicz/cbc-analisis-ingenieria

max_forks_repo_name: akielbowicz/cbc-analisis-ingenieria


text:
import PipeLayout: linecoefs, pwl_ineqs, pwl_inverse

@testset "linear interpolation of 2d points" begin
    # simple line through origin
    s, o = linecoefs(0, 0, 1, 2)
    @test s ≈ 2
    @test o ≈ 0

    # order of points doesn't matter
    s, o = linecoefs(1, 2, 0, 0)
    @test s ≈ 2
    @test o ≈ 0

    # neg. slope with offset
    s, o = linecoefs(2, 2, 4, 1)
    @test s ≈ -0.5
    @test o ≈ 3
end

@testset "inequalities for convex PWL function" begin
    # symmetric quadratic, with flat segment in the middle
    xs = [-5, -3, -1, 1, 3, 5]
    ineqs = pwl_ineqs(xs, xs.^2)
    @test size(ineqs) == (5, 3)
    @test ineqs[1,:] ./ ineqs[end,:]   ≈ [-1, 1, 1]
    @test ineqs[2,:] ./ ineqs[end-1,:] ≈ [-1, 1, 1]
    @test ineqs[3,:] ./ ineqs[3,2]     ≈ [ 0, 1, 1]
end

@testset "inverses from monotoneous PWL functions" begin
    @testset "monotoneous increasing" begin
        xs = Float64[1, 2, 3, 4]
        ys = xs.^2
        @test_throws AssertionError pwl_inverse(xs, ys, 0.5)
        @test pwl_inverse(xs, ys, 1.0) ≈ 1.0
        @test pwl_inverse(xs, ys, 2.5) ≈ 1.5
        @test pwl_inverse(xs, ys, 16.0) ≈ 4.0
        @test_throws AssertionError pwl_inverse(xs, ys, 16.5)
    end

    @testset "monotoneous decreasing" begin
        xs = Float64[1, 2, 3, 4]
        ys = -xs.^2
        @test_throws AssertionError pwl_inverse(xs, ys, -0.5)
        @test pwl_inverse(xs, ys, -1.0) ≈ 1.0
        @test pwl_inverse(xs, ys, -2.5) ≈ 1.5
        @test pwl_inverse(xs, ys, -16.0) ≈ 4.0
        @test_throws AssertionError pwl_inverse(xs, ys, -16.5)
    end

    @testset "with constant pieces" begin
        xs = Float64[1, 2, 3]
        ys = Float64[1, 1, 5]
        @test_throws AssertionError pwl_inverse(xs, ys, 0.5)
        @test pwl_inverse(xs, ys, 1) ≈ 1.0
        @test pwl_inverse(xs, ys, 3) ≈ 2.5
        @test pwl_inverse(xs, ys, 5) ≈ 3.0
        @test_throws AssertionError pwl_inverse(xs, ys, 5.5)
    end
end


text:
module problems

    include("str.jl")
    include("arr.jl")
    include("matrix.jl")
    include("linkedlist.jl")

end # module

max_stars_repo_name: bmoretz/Daily-Coding-Problem

max_issues_repo_name: bmoretz/Daily-Coding-Problem

max_forks_repo_name: bmoretz/Daily-Coding-Problem


text:
using Distributions
using ConjugatePriors

import ConjugatePriors: posterior, posterior_rand, posterior_mode, posterior_randmodel, fit_map

n = 100
w = rand(100)

# auxiliary tools

function ccount(K, x)
	r = zeros(K)
	for i = 1:length(x)
		r[x[i]] += 1.0
	end
	r
end

function ccount(K, x, w)
	r = zeros(K)
	for i = 1:length(x)
		r[x[i]] += w[i]
	end
	r
end

@testset "Beta - Bernoulli" begin

    pri = Beta(1.0, 2.0)

    x = rand(Bernoulli(0.3), n)
    p = posterior(pri, Bernoulli, x)
    @test isa(p, Beta)
    @test p.α ≈ pri.α + sum(x)
    @test p.β ≈ pri.β + (n - sum(x))

    f = fit_map(pri, Bernoulli, x)
    @test isa(f, Bernoulli)
    @test succprob(f) ≈ mode(p)

    p = posterior(pri, Bernoulli, x, w)
    @test isa(p, Beta)
    @test p.α ≈ pri.α + sum(x .* w)
    @test p.β ≈ pri.β + (sum(w) - sum(x .* w))

    f = fit_map(pri, Bernoulli, x, w)
    @test isa(f, Bernoulli)
    @test succprob(f) ≈ mode(p)

    @testset "posterior_rand & posterior_randmodel" begin

        pri = Beta(1.0, 2.0)
        x = rand(Bernoulli(0.3), n)
        post = posterior(pri, Bernoulli, x)

        pv = posterior_rand(pri, Bernoulli, x)
        @test isa(pv, Float64)
        @test 0. <= pv <= 1.

        pv = posterior_rand(pri, Bernoulli, x, w)
        @test isa(pv, Float64)
        @test 0. <= pv <= 1.

        pm = posterior_randmodel(pri, Bernoulli, x)
        @test isa(pm, Bernoulli)
        @test 0. <= succprob(pm) <= 1.

        pm = posterior_randmodel(pri, Bernoulli, x, w)
        @test isa(pm, Bernoulli)
        @test 0. <= succprob(pm) <= 1.
    end

end

@testset "Beta - Binomial" begin

    pri = Beta(1.0, 2.0)

    x = rand(Binomial(10, 0.3), n)
    p = posterior(pri, Binomial, (10, x))
    @test isa(p, Beta)
    @test p.α ≈ pri.α + sum(x)
    @test p.β ≈ pri.β + (10n - sum(x))

    f = fit_map(pri, Binomial, (10, x))
    @test isa(f, Binomial)
    @test ntrials(f) == 10
    @test succprob(f) ≈ mode(p)

    p = posterior(pri, Binomial, (10, x), w)
    @test isa(p, Beta)
    @test p.α ≈ pri.α + sum(x .* w)
    @test p.β ≈ pri.β + (10 * sum(w) - sum(x .* w))

    f = fit_map(pri, Binomial, (10, x), w)
    @test isa(f, Binomial)
    @test ntrials(f) == 10
    @test succprob(f) ≈ mode(p)

end

@testset "Dirichlet - Categorical" begin

    pri = Dirichlet([1., 2., 3.])

    x = rand(Categorical([0.2, 0.3, 0.5]), n)
    p = posterior(pri, Categorical, x)
    @test isa(p, Dirichlet)
    @test p.alpha ≈ pri.alpha + ccount(3, x)

    f = fit_map(pri, Categorical, x)
    @test isa(f, Categorical)
    @test probs(f) ≈ mode(p)

    p = posterior(pri, Categorical, x, w)
    @test isa(p, Dirichlet)
    @test p.alpha ≈ pri.alpha + ccount(3, x, w)

    f = fit_map(pri, Categorical, x, w)
    @test isa(f, Categorical)
    @test probs(f) ≈ mode(p)

end

@testset "Dirichlet - Multinomial" begin

    pri = Dirichlet([1., 2., 3.])

    x = rand(Multinomial(100, [0.2, 0.3, 0.5]), 1)
    p = posterior(pri, Multinomial, x)
    @test isa(p, Dirichlet)
    @test p.alpha ≈ pri.alpha + x

    r = posterior_mode(pri, Multinomial, x)
    @test r ≈ mode(p)

    x = rand(Multinomial(10, [0.2, 0.3, 0.5]), n)
    p = posterior(pri, Multinomial, x)
    @test isa(p, Dirichlet)
    @test p.alpha ≈ pri.alpha + vec(sum(x, dims=2))

    r = posterior_mode(pri, Multinomial, x)
    @test r ≈ mode(p)

    p = posterior(pri, Multinomial, x, w)
    @test isa(p, Dirichlet)
    @test p.alpha ≈ pri.alpha + vec(x * w)

    r = posterior_mode(pri, Multinomial, x, w)
    @test r ≈ mode(p)

end

@testset "Gamma - Exponential" begin

    pri = Gamma(1.5, 2.0)

    x = rand(Exponential(2.0), n)
    p = posterior(pri, Exponential, x)
    @test isa(p, Gamma)
    @test shape(p) ≈ shape(pri) + n
    @test rate(p) ≈ rate(pri) + sum(x)

    f = fit_map(pri, Exponential, x)
    @test isa(f, Exponential)
    @test rate(f) ≈ mode(p)

    p = posterior(pri, Exponential, x, w)
    @test isa(p, Gamma)
    @test shape(p) ≈ shape(pri) + sum(w)
    @test rate(p) ≈ rate(pri) + sum(x .* w)

    f = fit_map(pri, Exponential, x, w)
    @test isa(f, Exponential)
    @test rate(f) ≈ mode(p)

end


text:
#===============================================================================
# Installs the required packages for the project
#
# Author: Carl Cortright
# Date: 1/3/2017
#
===============================================================================#
using Conda

Pkg.add("PyCall")
Pkg.add("ScikitLearn")
Pkg.add("DataFrames")
Conda.add("scikit-learn")
Pkg.update() # Making sure everything is up-to-date

max_stars_repo_name: carlcortright/Numerai

max_issues_repo_name: carlcortright/Numerai

max_forks_repo_name: carlcortright/Numerai


text:
@static if is_apple()
    immutable Cpasswd
        pw_name::Cstring
        pw_passwd::Cstring
        pw_uid::Cint
        pw_gid::Cint
        pw_change::Cint
        pw_class::Cstring
        pw_gecos::Cstring
        pw_dir::Cstring
        pw_shell::Cstring
        pw_expire::Cint
        pw_fields::Cint
    end
elseif is_linux()
    immutable Cpasswd
       pw_name::Cstring
       pw_passwd::Cstring
       pw_uid::Cint
       pw_gid::Cint
       pw_gecos::Cstring
       pw_dir::Cstring
       pw_shell::Cstring
    end
else
    immutable Cpasswd
        pw_name::Cstring
        pw_uid::Cint
        pw_gid::Cint
        pw_dir::Cstring
        pw_shell::Cstring
    end

    Cpasswd() = Cpasswd(pointer("NA"), 0, 0, pointer("NA"), pointer("NA"))
end

immutable Cgroup
    gr_name::Cstring
    gr_passwd::Cstring
    gr_gid::Cint
end

Cgroup() = Cgroup(pointer("NA"), pointer("NA"), 0)

immutable User
    name::String
    uid::UInt64
    gid::UInt64
    dir::String
    shell::String

    function User(ps::Cpasswd)
        new(
            unsafe_string(ps.pw_name),
            UInt64(ps.pw_uid),
            UInt64(ps.pw_gid),
            unsafe_string(ps.pw_dir),
            unsafe_string(ps.pw_shell)
        )
    end
end

User(passwd::Ptr{Cpasswd}) = User(unsafe_load(passwd))

function Base.show(io::IO, user::User)
    print(io, "$(user.uid) ($(user.name))")
end

function User(name::String)
    ps = @static if is_unix()
        ccall((:getpwnam, "libc"), Ptr{Cpasswd}, (Ptr{UInt8},), name)
    else
        Cpasswd()
    end

    User(ps)
end

function User(uid::UInt)
    ps = @static if is_unix()
        ccall((:getpwuid, "libc"), Ptr{Cpasswd}, (UInt64,), uid)
    else
        Cpasswd()
    end

    User(ps)
end

function User()
    uid = @static is_unix() ? ccall((:geteuid, "libc"), Cint, ()) : 0
    User(UInt64(uid))
end

immutable Group
    name::String
    gid::UInt64

    Group(gr::Cgroup) = new(unsafe_string(gr.gr_name), UInt64(gr.gr_gid))
end

Group(group::Ptr{Cgroup}) = Group(unsafe_load(group))

function Base.show(io::IO, group::Group)
    print(io, "$(group.gid) ($(group.name))")
end

function Group(name::String)
    ps = @static if is_unix()
        ccall((:getgrnam, "libc"), Ptr{Cgroup}, (Ptr{UInt8},), name)
    else
        Cgroup()
    end

    Group(ps)
end

function Group(gid::UInt)
    gr = @static if is_unix()
        ccall((:getgrgid, "libc"), Ptr{Cgroup}, (UInt64,), gid)
    else
        Cgroup()
    end

    Group(gr)
end

function Group()
    gid = @static is_unix() ? ccall((:getegid, "libc"), Cint, ()) : 0
    Group(UInt64(gid))
end


text:
# ---
# title: 1457. Pseudo-Palindromic Paths in a Binary Tree
# id: problem1457
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Bit Manipulation, Tree, Depth-first Search
# link: <https://leetcode.com/problems/pseudo-palindromic-paths-in-a-binary-tree/description/>
# hidden: true
# ---
#
# Given a binary tree where node values are digits from 1 to 9. A path in the
# binary tree is said to be **pseudo-palindromic** if at least one permutation
# of the node values in the path is a palindrome.
#
# _Return the number of **pseudo-palindromic** paths going from the root node to
# leaf nodes._
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/05/06/palindromic_paths_1.png)
#
#
#
#     Input: root = [2,3,1,3,1,null,1]
#     Output: 2
#     Explanation: The figure above represents the given binary tree. There are three paths going from the root node to leaf nodes: the red path [2,3,3], the green path [2,1,1], and the path [2,3,1]. Among these paths only red path and green path are pseudo-palindromic paths since the red path [2,3,3] can be rearranged in [3,2,3] (palindrome) and the green path [2,1,1] can be rearranged in [1,2,1] (palindrome).
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/05/07/palindromic_paths_2.png)**
#
#
#
#     Input: root = [2,1,1,1,3,null,null,null,null,null,1]
#     Output: 1
#     Explanation: The figure above represents the given binary tree. There are three paths going from the root node to leaf nodes: the green path [2,1,1], the path [2,1,3,1], and the path [2,1]. Among these paths only the green path is pseudo-palindromic since [2,1,1] can be rearranged in [1,2,1] (palindrome).
#
#
# **Example 3:**
#
#
#
#     Input: root = [9]
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * The given binary tree will have between `1` and `10^5` nodes.
#   * Node values are digits from `1` to `9`.
#
#
## @lc code=start
using LeetCode

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


text:
module DeclarativeFormatting

export FormatRequirement, IncompatibleFormatException, allowable, best,
       takenonzero, conflict, reconcile, FormatSpecification

#= Format Specification rules =#
abstract type FormatRequirement end

mutable struct IncompatibleFormatException <: Exception
    msg::String
end

#= Format helpers =#
function allowable(x, y)
    if isempty(x)
        return y
    elseif isempty(y)
        return x
    end
    merged = typeof(x)()
    for k in keys(x)
        if haskey(y, k)
            merged[k] = x[k] + y[k]
        end
    end
    merged
end

function best(x)
    top = maximum(values(x))
    for (k, v) in x
        if v == top
            return k
        end
    end
end

function takenonzero(a, b, zero)
    if a == zero
        b
    elseif b == zero
        a
    elseif a == b
        a
    else throw(IncompatibleFormatException("Can't reconcile $a and $b")) end
end

"""
Return true if two `FormatRequirement`s conflict.
"""
conflict(x, y) = typeof(x) ≡ typeof(y)

"""
Reconcile two conflicting `FormatRequirement`s if possible.
"""
function reconcile end

"""
A collection of `FormatRequirement`s.
"""
struct FormatSpecification
    reqs::Vector{FormatRequirement}
end

function Base.push!(spec::FormatSpecification, nextreq::FormatRequirement)
    for (i, req) in enumerate(spec.reqs)
        if conflict(nextreq, req)
            # zap the conflicting requirement
            deleteat!(spec.reqs, i)

            # reconcile the two conflicting requirements
            # then add the mutally agreed requirements
            altreqs = reconcile(nextreq, req)
            append!(spec, altreqs)
            return spec
        end
    end
    # no conflict, return with new requirement
    push!(spec.reqs, nextreq)
    spec
end
function Base.append!(spec::FormatSpecification, reqs)
    for req in reqs
        push!(spec, req)
    end
    spec
end

function Base.union(x::FormatSpecification, y::FormatSpecification)
    result = FormatSpecification(copy(x.reqs))
    for req in y.reqs
        push!(result, req)
    end
    result
end

function Base.get(spec::FormatSpecification, ::Type{T}, d) where T<:FormatRequirement
    for req in spec.reqs
        if isa(req, T)
            return req
        end
    end
    d
end

end  # module DeclarativeFormat


text:
module BMatching
using Random, Statistics

using UnPack
using Printf
using Erdos

export run_bp

getref(v::Vector, i::Integer) = pointer(v, i)

include("bmatching_bp.jl")

end #module


text:
instructions = open("navigation_data.txt") do f
    map(r -> (r[1], parse(Int, r[2:end])), readlines(f))
end

function rotate(v::Vector{Int64}, a::Int64)::Vector{Int64}
    b = deg2rad(a)
    x1, y1 = v
    x2 = cos(b)*x1 - sin(b)*y1
    y2 = sin(b)*x1 + cos(b)*y1

    return round.(Int, [x2, y2])
end

ship = [0, 0]
direction = 0
for (cmd, val) in instructions
    global ship, direction

    if cmd == 'N'
        ship += [0, val]
    elseif cmd == 'E'
        ship += [val, 0]
    elseif cmd == 'S'
        ship += [0, -val]
    elseif cmd == 'W'
        ship += [-val, 0]
    elseif cmd == 'F'
        ship += val * round.(Int, [cos(deg2rad(direction)), sin(deg2rad(direction))])
    elseif cmd == 'L'
        direction = (direction + val) % 360
    elseif cmd == 'R'
        direction = (direction - val) % 360
    end
end

println(sum(abs.(ship)))

waypoint = [10, 1] # relative to ship
ship = [0, 0]
for (cmd, val) in instructions
    global ship, waypoint

    if cmd == 'N'
        waypoint += [0, val]
    elseif cmd == 'E'
        waypoint += [val, 0]
    elseif cmd == 'S'
        waypoint += [0, -val]
    elseif cmd == 'W'
        waypoint += [-val, 0]
    elseif cmd == 'F'
        ship += waypoint * val
    elseif cmd == 'L'
        waypoint = rotate(waypoint, val)
    elseif cmd == 'R'
        waypoint = rotate(waypoint, -val)
    end
    # println("$(i): ship: $(ship), wp: $(waypoint)")
end

println(sum(abs.(ship)))


text:
# Regex to validate a Roman numeral
const VALID_ROMAN_PATTERN =
    r"""
    ^\s*                   # Skip leading whitespace
    (
    M*                     # Thousands
    (C{0,9}|CD|DC{0,4}|CM) # Hundreds
    (X{0,9}|XL|LX{0,4}|XC) # Tens
    (I{0,9}|IV|VI{0,4}|IX) # Ones
    )
    \s*$                   # Skip trailing whitespace
    """ix # Be case-insensitive and verbose

const OVERBAR_CHAR = Char(0x0305)
_place_overbar(c::AbstractChar) = string(c, OVERBAR_CHAR)
_place_overbar(s::AbstractString) = join(_place_overbar(c) for c in s)

const NUMERAL_MAP = [
    (1_000_000, _place_overbar("M"))
    (500_000, _place_overbar("D"))
    (100_000, _place_overbar("C"))
    (50_000, _place_overbar("L"))
    (10_000, _place_overbar("X"))
    (9000, _place_overbar("IX"))
    (8000, _place_overbar("VIII"))
    (7000, _place_overbar("VII"))
    (6000, _place_overbar("VI"))
    (5000, _place_overbar("V"))
    (4000, _place_overbar("IV"))
    (1000, "M")
    (900,  "CM")
    (500,  "D")
    (400,  "CD")
    (100,  "C")
    (90,   "XC")
    (50,   "L")
    (40,   "XL")
    (10,   "X")
    (9,    "IX")
    (5,    "V")
    (4,    "IV")
    (1,    "I")
]

function fromroman(str::AbstractString)
    m = match(VALID_ROMAN_PATTERN, str)
    m ≡ nothing && throw(Meta.ParseError(str * " is not a valid roman numeral"))
    # Strip whitespace and make uppercase
    str = uppercase(m.captures[1])
    i = 1
    val = 0
    strlen = length(str)
    for (num_val, numeral) in NUMERAL_MAP
        numlen = length(numeral)
        while i+numlen-1 <= strlen && str[i:i+numlen-1] == numeral
            val += num_val
            i += numlen
        end
    end
    return val
end

import Base: parse
parse(::Type{RomanNumeral}, s::AbstractString) = RomanNumeral(fromroman(s))

# using Compat: @warn
function toroman(val::Integer)
    val <= 0 && throw(DomainError(val, "in ancient Rome there were only strictly positive numbers"))
    val > maximum(first(t) for t in NUMERAL_MAP) && @warn("Roman numerals do not handle large numbers well") # instead if maximum, use NUMERAL_MAP[1][1]?

    str = IOBuffer()
    for (num_val, numeral) in NUMERAL_MAP
        i = div(val, num_val)
        # Never concatenate an empty string to `str`
        i == 0 && continue
        print(str, repeat(numeral,i))
        val -= i*num_val
        # Stop when ready
        val == 0 && break
    end
    String(take!(str))
end


text:
macro _checkStatus(stat,msg)
   return quote
      if ($(esc(stat)) == 0)
         throw(CALCEPHException($(esc(msg))))
      end
   end
end

macro _checkPointer(ptr,msg)
   return quote
      if ($(esc(ptr)) == C_NULL)
         throw(CALCEPHException($(esc(msg))))
      end
   end
end

macro _checkOrder(order)
   return quote
      local or = $(esc(order))
      if (or<0) || (or>3)
         throw(CALCEPHException("Order must be between 0 and 3."))
      end
   end
end

"""
    Ephem

  Ephemeris descriptor. Create with:

    eph = Ephem(filename)
    eph = Ephem([filename1,filename2...])

  The ephemeris descriptor will be used to access the ephemeris and related
  data stored in the specified files.

  Because, Julia GC is lazy, you may want to free the memory managed by eph
  before you get rid of the reference to eph with:

    finalize(eph)

  or after by forcing the GC to run:

    gc()

"""
mutable struct Ephem
   data :: Ptr{Cvoid}
   function Ephem(files::Vector{<:AbstractString})
      ptr = ccall((:calceph_open_array, libcalceph), Ptr{Cvoid},
                  (Int, Ptr{Ptr{UInt8}}), length(files), files)
      @_checkPointer ptr "Unable to open ephemeris file(s)!"
      obj = new(ptr)
      finalizer(_ephemDestructor,obj) # register object destructor
      return obj
   end
end

# to be called by gc when cleaning up
# not in the exposed interface but can be called with finalize(e)
function _ephemDestructor(eph::Ephem)
   if (eph.data == C_NULL)
      return
   end
   ccall((:calceph_close, libcalceph), Cvoid, (Ptr{Cvoid},), eph.data)
   eph.data = C_NULL
   return
end

Ephem(file::AbstractString) = Ephem([file])

"""
    prefetch(eph)

  This function prefetches to the main memory all files associated to the ephemeris descriptor eph.

"""
function prefetch(eph::Ephem)
    @_checkPointer eph.data "Ephemeris is not properly initialized!"
    stat = ccall((:calceph_prefetch, libcalceph), Int, (Ptr{Cvoid},), eph.data)
    @_checkStatus stat "Unable to prefetch  ephemeris!"
    return
end


text:
module CorrelationTrackers
using Reexport: @reexport
@reexport using AnnealingAPI
using CorrelationFunctions
using Base.Iterators: zip, countfrom, takewhile, take
using CircularArrays: CircularArray
using LinearAlgebra: norm
using Images: imgradients, KernelFactors
using Interpolations: interpolate, extrapolate, Gridded, Linear, Line

include("slices.jl")
include("tracker.jl")
include("extrapolation.jl")
include("helpers.jl")

export
    CorrelationTracker,
    ExtrapolatedData,
    default_trackers,
    # Extra trackers
    CCTracker, cross_correlation
end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ e3b54b66-d8cd-430f-af1a-5258b8cdc4ba
using Base: SimpleVector, show_can_elide, isgensym, unwrap_unionall

# ╔═╡ c6955676-657f-4413-81cc-e0e2fc8dc6d8
function stripmodules(s::Symbol)
	split(string(s),'.')[end]  |> Symbol
end

# ╔═╡ 2ae5b44c-678a-42da-934f-7ff057704381
md"""
# \_toexpr
"""

# ╔═╡ 69c2a190-0a1b-4b87-a143-ab9b890a065b
md"""
The function `_toexpr` is used to process the components of a method signature to reconstruct an expression that could be used to bring the method into scope from the parent module (loaded with `@plutoinclude`) to the current module
"""

# ╔═╡ 4fcb74bc-f9ed-4969-b9cb-44c0991a2788
# This function is basicaly copied and adapted from https://github.com/JuliaLang/julia/blob/743a37898d447d047002efcc19ce59825ef63cc1/base/show.jl#L604-L648
function _toexpr(v::Val, env::SimpleVector, orig::SimpleVector, wheres::Vector; to::Module, from::Module, importedlist::Vector{Symbol}, fromname::Symbol)
	ex = Expr(:curly)
	n = length(env)
    elide = length(wheres)
    function egal_var(p::TypeVar, @nospecialize o)
        return o isa TypeVar &&
            ccall(:jl_types_egal, Cint, (Any, Any), p.ub, o.ub) != 0 &&
            ccall(:jl_types_egal, Cint, (Any, Any), p.lb, o.lb) != 0
    end
	println("n = $n")
    for i = n:-1:1
        p = env[i]
		println("p = $p")
        if p isa TypeVar
            if i == n && egal_var(p, orig[i]) && show_can_elide(p, wheres, elide, env, i)
                n -= 1
                elide -= 1
            elseif p.lb === Union{} && isgensym(p.name) && show_can_elide(p, wheres, elide, env, i)
                elide -= 1
            elseif p.ub === Any && isgensym(p.name) && show_can_elide(p, wheres, elide, env, i)
                elide -= 1
            end
        end
    end
	if n > 0
        for i = 1:n
            p = env[i]
            if p isa TypeVar
                if p.lb === Union{} && something(findfirst(@nospecialize(w) -> w === p, wheres), 0) > elide
                    push!(ex.args, Expr(:(<:), _toexpr(v, p.ub;to, from, importedlist, fromname)))
                elseif p.ub === Any && something(findfirst(@nospecialize(w) -> w === p, wheres), 0) > elide
                    push!(ex.args, Expr(:(>:), _toexpr(v, p.lb; to, from, importedlist, fromname)))
                else
                    push!(ex.args, _toexpr(v, p; to, from, importedlist, fromname))
                end
            else
               push!(ex.args, _toexpr(v, p; to, from, importedlist, fromname))
            end
        end
    end
    resize!(wheres, elide)
    ex
end

# ╔═╡ bc123a56-31b7-4599-8acb-19cdd2934f30
# function _toexpr(v::Val, x::Type; to::Module, from::Module, importedlist::Vector{Symbol}, fromname::Symbol)
# 	alias = Base.make_typealias(x)
# 	if alias === nothing
# 		# This is not an alias
# 	else
# 		wheres = Base.make_wheres()
# 	end
# end

# ╔═╡ 6d20a40c-7407-4bbe-ba97-b05a541efb46
# function _toexpr(v::Val, x::DataType, wheres::Vector = TypeVar[]; to::Module, from::Module, importedlist::Vector{Symbol}, fromname::Symbol)
# 	parameters = x.parameters::SimpleVector
#     name = x.name.wrapper |> Symbol |> stripmodules
# 	# println("name = $name, stripped = $(stripmodules(name))")
# 	val = (isdefined(to,name) || name ∈ importedlist) ? name : :($fromname.$name)
# 	if isempty(parameters)
# 		return val
# 	end
# 	orig = if v isa Val{:wheres}
# 		unwrap_unionall(x.name.wrapper).parameters
# 	elseif v isa Val{:types}
# 		parameters
# 	else
# 		error("Unsupported Val direction")
# 	end
# 	ex = _toexpr(v, parameters, orig, wheres; to, from, importedlist, fromname)
# 	if isempty(ex.args)
# 		return val
# 	else
# 		pushfirst!(ex.args, val)
# 	end
# 	return ex
# end

# ╔═╡ 3b18e1ce-9b6a-4c81-b78c-d1618a2fd4cb
# function _toexpr(v::Val, x::UnionAll; to::Module, from::Module, importedlist::Vector{Symbol}, fromname::Symbol)
# 	wheres = TypeVar[]
# 	while x isa UnionAll
# 		push!(wheres,x.var)
# 		x = x.body
# 	end
# 	ex = _toexpr(v, x, wheres; to, from, importedlist, fromname)
# end

# ╔═╡ 94b3e64c-fd50-410e-9619-9122e6ef35cd
function _toexpr(v::Val, x; to::Module, from::Module, importedlist::Vector{Symbol}, fromname::Symbol)
	return x
end

# ╔═╡ 72481e47-60cc-4e86-8cbd-9bbb7576c995
function _toexpr(v::Val{:types},x::TypeVar; to::Module, from::Module, importedlist::Vector{Symbol}, fromname::Symbol)
	# If this has a name, return the name
	isgensym(x.name) || return x.name
end

# ╔═╡ 99336a92-3510-467d-a2cb-5cd4d3758406
function _toexpr(v::Val{:wheres},x::TypeVar; to::Module, from::Module, importedlist::Vector{Symbol}, fromname::Symbol)
	if isgensym(x.name)
		if x.lb === Union{}
			ex = Expr(:(<:), _toexpr(v,x.ub; to, from, importedlist, fromname))
		elseif x.ub === Any
			ex = Expr(:(:>), _toexpr(v,x.lb; to, from, importedlist, fromname))
		else
			ex = ()
		end
	else
		if x.lb === Union{} && x.ub === Any
			ex = x.name
		elseif x.lb === Union{}
			ex = Expr(:(<:), x.name, _toexpr(v,x.ub; to, from, importedlist, fromname))
		elseif x.ub === Any
			ex = Expr(:(:>), x.name, _toexpr(v,x.lb; to, from, importedlist, fromname))
		else
			ex = Expr(:comparison,_toexpr(v,x.lb; to, from, importedlist, fromname),:(<:),x.name,:(<:),_toexpr(v,x.ub; to, from, importedlist, fromname))
		end
	end
	return ex
end

# ╔═╡ 782a1a5d-b04f-43d5-94ac-71e4b91b0222
function _toexpr(v::Val,u::Union; to::Module, from::Module, importedlist::Vector{Symbol}, fromname::Symbol)
	ex = Expr(:curly)
	push!(ex.args,:Union)
	push!(ex.args, _toexpr(v,u.a; to, from, importedlist, fromname))
	push!(ex.args, _toexpr(v,u.b; to, from, importedlist, fromname))
	ex
end

# ╔═╡ 6c270dd6-4f9d-49cb-8b09-ff30d1dc4bd6
md"""
# Other Helper Funcs
"""

# ╔═╡ d895c070-d923-4f01-bfd8-e9129b5f290d
function _method_expr(mtd::Method; to::Module, from::Module, importedlist::Vector{Symbol}, fromname::Symbol)
	s = mtd.name
	lhs = Expr(:call)
	# Add the method name
	# push!(lhs.args,mtd.name)
	push!(lhs.args,s)
	nms = map(Base.method_argnames(mtd)[2:end]) do nm
		nm === Symbol("#unused#") ? gensym() : nm
	end
	tv = Any[]
    sig = mtd.sig
    while isa(sig, UnionAll)
        push!(tv, sig.var)
        sig = sig.body
    end
	# Get the argument types, stripped from TypeVars
	tps = sig.parameters[2:end]
	for (nm,sig) ∈ zip(nms,tps)
		push!(lhs.args, Expr(:(::),_toexpr(Val(:types), nm; to, from, importedlist, fromname),_toexpr(Val(:types), sig; to, from, importedlist, fromname)))
	end
	if !isempty(tv)
		lhs = Expr(:where,lhs,map(x -> _toexpr(Val(:wheres),x; to, from, importedlist, fromname),tv)...)
		# lhs = Expr(:where,lhs,tv...)
	end
	lhs
	# Add the function call
	rhs = :($fromname.$s())
	# Push the variables
	for (nm,sig) ∈ zip(nms,tps)
		if sig isa Core.TypeofVararg
			push!(rhs.args, Expr(:(...),nm))
		else
			push!(rhs.args, nm)
		end
	end
	rhs = Expr(:block, rhs)
	Expr(:(=), lhs, rhs)
end

# ╔═╡ 3381acd4-4d69-45d2-95c9-b555339161c0
function _copymethods!(ex::Expr, s::Symbol; to::Module, from::Module, importedlist::Vector{Symbol}, fromname::Symbol)
	f = getfield(from,s)
	ml = methods(f,from)
	for mtd ∈ ml
		push!(ex.args, _method_expr(mtd; to, from, importedlist, fromname))
	end
	ex
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0-rc2"
manifest_format = "2.0"

[deps]
"""

# ╔═╡ Cell order:
# ╠═e3b54b66-d8cd-430f-af1a-5258b8cdc4ba
# ╠═c6955676-657f-4413-81cc-e0e2fc8dc6d8
# ╠═2ae5b44c-678a-42da-934f-7ff057704381
# ╠═69c2a190-0a1b-4b87-a143-ab9b890a065b
# ╠═4fcb74bc-f9ed-4969-b9cb-44c0991a2788
# ╠═bc123a56-31b7-4599-8acb-19cdd2934f30
# ╠═6d20a40c-7407-4bbe-ba97-b05a541efb46
# ╠═3b18e1ce-9b6a-4c81-b78c-d1618a2fd4cb
# ╠═94b3e64c-fd50-410e-9619-9122e6ef35cd
# ╠═72481e47-60cc-4e86-8cbd-9bbb7576c995
# ╠═99336a92-3510-467d-a2cb-5cd4d3758406
# ╠═782a1a5d-b04f-43d5-94ac-71e4b91b0222
# ╠═6c270dd6-4f9d-49cb-8b09-ff30d1dc4bd6
# ╠═3381acd4-4d69-45d2-95c9-b555339161c0
# ╠═d895c070-d923-4f01-bfd8-e9129b5f290d
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002