text:
Base.maximum(f, x::AbstractArray{Union{T,Missing},1}; threads = false) where T <: Union{INTEGERS, FLOATS, TimeType} = threads ? hp_maximum(f, x) : stat_maximum(f, x)
Base.maximum(x::AbstractArray{Union{T,Missing},1}; threads = false) where T <: Union{INTEGERS, FLOATS, TimeType}  = threads ? hp_maximum(identity, x) : stat_maximum(identity, x)

Base.minimum(f, x::AbstractArray{Union{T,Missing},1}; threads = false) where T <: Union{INTEGERS, FLOATS, TimeType} = threads ? hp_minimum(f, x) : stat_minimum(f, x)
Base.minimum(x::AbstractArray{Union{T,Missing},1}; threads = false) where T <: Union{INTEGERS, FLOATS, TimeType}= threads ? hp_minimum(identity, x) : stat_minimum(identity, x)

# not optimised for simd
Base.maximum(f, x::AbstractVector{Union{Missing, T}}) where T <: AbstractString = mapreduce(f, _stat_max_fun, x)
Base.minimum(f, x::AbstractVector{Union{Missing, T}}) where T <: AbstractString = mapreduce(f, _stat_min_fun, x)
Base.maximum(x::AbstractVector{Union{Missing, T}}) where T <: AbstractString = maximum(identity, x)
Base.minimum(x::AbstractVector{Union{Missing, T}}) where T <: AbstractString = minimum(identity, x)

Base.sum(f, x::AbstractArray{Union{T,Missing},1}; threads = false) where T <: Union{INTEGERS, FLOATS} = threads ? hp_sum(f, x) : stat_sum(f, x)
Base.sum(x::AbstractArray{Union{T,Missing},1}; threads = false) where T <: Union{INTEGERS, FLOATS} = threads ? hp_sum(identity, x) : stat_sum(identity, x)

Statistics.mean(f, x::AbstractArray{Union{T,Missing},1}) where T <: Union{INTEGERS, FLOATS} = stat_mean(f, x)
Statistics.mean(x::AbstractArray{Union{T,Missing},1}) where T <: Union{INTEGERS, FLOATS} = stat_mean(x)

wsum(f, x::AbstractVector, w::AbstractVector) = stat_wsum(f, x, w)
wsum(x::AbstractVector, w::AbstractVector) = stat_wsum(identity, x, w)

wmean(f, x::AbstractVector, w::AbstractVector) = stat_wmean(f, x, w)
wmean(x::AbstractVector, w::AbstractVector) = stat_wmean(identity, x, w)

Statistics.var(f, x::AbstractArray{Union{T,Missing},1}, dof = true) where T <: Union{INTEGERS, FLOATS}= stat_var(f, x, dof)
Statistics.var(x::AbstractArray{Union{T,Missing},1}, dof = true) where T <: Union{INTEGERS, FLOATS}= stat_var(x, dof)

Statistics.std(f, x::AbstractArray{Union{T,Missing},1}, dof = true) where T <: Union{INTEGERS, FLOATS}= stat_std(f, x, dof)
Statistics.std(x::AbstractArray{Union{T,Missing},1}, dof = true) where T <: Union{INTEGERS, FLOATS}= stat_std(x, dof)

Statistics.median(x::AbstractArray{Union{T,Missing},1}) where T = stat_median(x)
Statistics.median!(x::AbstractArray{Union{T,Missing},1}) where T = stat_median!(x)

Base.extrema(f, x::AbstractArray{Union{T,Missing},1}; threads = false) where T <: Union{INTEGERS, FLOATS, TimeType} = threads ? (hp_minimum(f, x), hp_maximum(f, x)) : (stat_minimum(f, x), stat_maximum(f, x))
Base.extrema(x::AbstractArray{Union{T,Missing},1}; threads = false) where T  <: Union{INTEGERS, FLOATS, TimeType}= threads ? (hp_minimum(identity, x), hp_maximum(identity, x)) : (stat_minimum(identity, x), stat_maximum(identity, x))

Base.argmax(f, x::AbstractArray{Union{T,Missing},1}) where T <: Union{INTEGERS, FLOATS, TimeType, AbstractString} = stat_findmax(f, x)[2]
Base.argmax(x::AbstractArray{Union{T,Missing},1}) where T <: Union{INTEGERS, FLOATS, TimeType, AbstractString} = stat_findmax(x)[2]
Base.argmin(f, x::AbstractArray{Union{T,Missing},1}) where T <: Union{INTEGERS, FLOATS, TimeType, AbstractString} = stat_findmin(f, x)[2]
Base.argmin(x::AbstractArray{Union{T,Missing},1}) where T <: Union{INTEGERS, FLOATS, TimeType, AbstractString} = stat_findmin(x)[2]

Base.findmax(f, x::AbstractArray{Union{T,Missing},1}) where T <: Union{INTEGERS, FLOATS, TimeType, AbstractString} = stat_findmax(f, x)
Base.findmax(x::AbstractArray{Union{T,Missing},1}) where T <: Union{INTEGERS, FLOATS, TimeType, AbstractString} = stat_findmax(x)
Base.findmin(f, x::AbstractArray{Union{T,Missing},1}) where T <: Union{INTEGERS, FLOATS, TimeType, AbstractString} = stat_findmin(f, x)
Base.findmin(x::AbstractArray{Union{T,Missing},1}) where T <: Union{INTEGERS, FLOATS, TimeType, AbstractString} = stat_findmin(x)

function Base.cumsum(x::AbstractArray{Union{T,Missing},1}; missings = :ignore) where T <: Union{INTEGERS, FLOATS}
    if missings == :ignore
        stat_cumsum_ignore(x)
    elseif missings == :skip
        stat_cumsum_skip(x)
    else
        throw(ArgumentError("`missings` must be either `:ignore` or `:skip`"))
    end
end
function Base.cumsum!(outx::AbstractVector, x::AbstractArray{Union{T,Missing},1}; missings = :ignore) where T <: Union{INTEGERS, FLOATS}
    if missings == :ignore
        stat_cumsum!_ignore(outx, x)
    elseif missings == :skip
        stat_cumsum!_skip(outx, x)
    else
        throw(ArgumentError("`missings` must be either `:ignore` or `:skip`"))
    end
end
function Base.cumprod(x::AbstractArray{Union{T,Missing},1}; missings = :ignore) where T <: Union{INTEGERS, FLOATS}
    if missings == :ignore
        stat_cumprod_ignore(x)
    elseif missings == :skip
        stat_cumprod_skip(x)
    else
        throw(ArgumentError("`missings` must be either `:ignore` or `:skip`"))
    end
end
function Base.cumprod!(outx::AbstractVector, x::AbstractArray{Union{T,Missing},1}; missings = :ignore) where T <: Union{INTEGERS, FLOATS}
    if missings == :ignore
        stat_cumprod!_ignore(outx, x)
    elseif missings == :skip
        stat_cumprod!_skip(outx, x)
    else
        throw(ArgumentError("`missings` must be either `:ignore` or `:skip`"))
    end
end

function cummin(x::AbstractArray{Union{T,Missing},1}; missings = :ignore) where T <: Union{INTEGERS, FLOATS, TimeType}
    if missings == :ignore
        stat_cummin_ignore(x)
    elseif missings == :skip
        stat_cummin_skip(x)
    else
        throw(ArgumentError("`missings` must be either `:ignore` or `:skip`"))
    end
end
function cummin!(outx, x::AbstractArray{Union{T,Missing},1}; missings = :ignore) where T <: Union{INTEGERS, FLOATS, TimeType}
    if missings == :ignore
        stat_cummin!_ignore(outx, x)
    elseif missings == :skip
        stat_cummin!_skip(outx, x)
    else
        throw(ArgumentError("`missings` must be either `:ignore` or `:skip`"))
    end
end
function cummax(x::AbstractArray{Union{T,Missing},1}; missings = :ignore) where T <: Union{INTEGERS, FLOATS, TimeType}
    if missings == :ignore
        stat_cummax_ignore(x)
    elseif missings == :skip
        stat_cummax_skip(x)
    else
        throw(ArgumentError("`missings` must be either `:ignore` or `:skip`"))
    end
end
function cummax!(outx, x::AbstractArray{Union{T,Missing},1}; missings = :ignore) where T <: Union{INTEGERS, FLOATS, TimeType}
    if missings == :ignore
        stat_cummax!_ignore(outx, x)
    elseif missings == :skip
        stat_cummax!_skip(outx, x)
    else
        throw(ArgumentError("`missings` must be either `:ignore` or `:skip`"))
    end
end


text:
using LinearAlgebra
using TensorOperations
import Strided: UnsafeStridedView

@testset "Arrays" begin

    A = FermiMDrand(5,5)
    @test FermiMDArray(A) === A
    @test size(A) == (5,5)
    @test A[2,2] == A.data[2,2]
    A[2,2] = 0.75
    @test A[2,2] == 0.75
    @test ndims(A) == 2
    @test length(A) == 25
    At = similar(A)
    At .= transpose(A)
    @test permutedims(A, (2,1)) == At
    permutedims!(A, (2,1))
    @test A == At

    @test diag(A) == diag(A.data)

    # Make it Hermitian
    A = A + transpose(A)

    @test A^2 == A.data^2
    @test A^-1 ≈ A^-1.0

    Hvals, Hvec = diagonalize(Hermitian(A), hermitian=true)
    adjoint(Hvec) * A*Hvec ≈ diagm(Hvals)
    vals, vec = diagonalize(A, hermitian=false)
    vec^-1 * A*vec ≈ diagm(vals)

    b = rand(5,5)
    B = FermiMDArray(b)
    @test A + b == A + B == b + A
    @test A - b == A - B == -(b - A)
    @test 2*A == 2*A.data == A*2

    @test FermiMDArray(2.0) == 2.0

    C = FermiMDzeros(5,5)
    @test C == zeros(5,5)

    Bv = similar(B, (5,))
    Bv .= B[1,:]
    @test Bv == B.data[1,:]

    Ct = rand(Float32, (5,5))
    C .= Ct
    @test C.data == Base.convert(Array{Float64}, Ct)

    @test UnsafeStridedView(A) == UnsafeStridedView(A.data)

    idx = [ (1,1,1,1), (2,2,2,2), (3,3,3,3)]
    val = rand(3)
    A = FermiSparse(idx, val)
    @test ndims(A) == 4

    # x = @capture_out display(A)
    # println("")
    # @test occursin("Fermi Memory-held Dense Array", x)
end

ext: jl


text:
using Test, GitHubAPI

@testset "Licenses" begin
    @test isa(get_licenses(), Vector{String})
end
@testset "Search" begin
    @test parse_repos!("postgresql")
    @test parse_repos!("postgresql",
                       "2019-10-29T14:37:16+00:00..2019-08-15T00:00:00-04:00")
end
@testset "Repository" begin
    @test parse_commits!("uva-bi-sdad/GitHubAPI.jl")
    @test parse_commits!("uva-bi-sdad/GitHubAPI.jl", "2018-08-12 17:29:38+00")
end


text:
ROS.@include geometry_msgs: TransformStamped, PoseStamped, PointStamped, Vector3Stamped, Quaternion

ROS.init("testTF")
nh = ROS.NodeHandle()

b = ROS.tf.TransformBroadcaster()
sb = ROS.tf.StaticTransformBroadcaster()
bf = ROS.tf.Buffer()
tl = ROS.tf.TransformListener(bf)

ts = ROS.geometry_msgs_TransformStamped()
ts.header.stamp = ROS.now()
ts.header.frame_id = "base_link"
ts.child_frame_id = "gps_link"
ts.transform.translation.y = 1.6
ts.transform.rotation.w = 1
ts2 = ROS.geometry_msgs_TransformStamped()
ts2.header.stamp = ROS.now()
ts2.header.frame_id = "base_link"
ts2.child_frame_id = "camera_link"
ts2.transform.translation.x = 0.3
ts2.transform.rotation.w = 1
tss = [ts,ts2]

ROS.tf.sendTransform(b, ts)
ROS.tf.sendTransform(b, tss)
ROS.tf.sendTransform(sb, ts)
ROS.tf.sendTransform(sb, tss)

@test ROS.tf.canTransform(bf, "base_link", "camera_link", ROS.Time(0), ROS.Duration(1))
@test ROS.tf.canTransform(bf, "gps_link", ROS.Time(0), "camera_link", ROS.Time(0), "base_link", ROS.Duration(1))

ROS.tf.lookupTransform(bf, "base_link", "camera_link", ROS.Time(0), ROS.Duration(1))
ROS.tf.lookupTransform(bf, "gps_link", ROS.Time(0), "camera_link", ROS.Time(0), "base_link", ROS.Duration(1))

p = ROS.geometry_msgs_PoseStamped()
p.header.frame_id = "camera_link"
p.pose.position.x = 1
p.pose.position.y = 0.6
p.pose.orientation.w = 1

ROS.tf.transform(bf, p, "gps_link", ROS.Time(0), "base_link", ROS.Duration(2))
ROS.tf.transform(bf, p, "base_link", ROS.Duration(2))

p = ROS.geometry_msgs_PointStamped()
p.header.frame_id = "camera_link"
p.point.x = 1
p.point.y = 0.6

ROS.tf.transform(bf, p, "gps_link", ROS.Time(0), "base_link", ROS.Duration(2))
ROS.tf.transform(bf, p, "base_link", ROS.Duration(2))

p = ROS.geometry_msgs_Vector3Stamped()
p.header.frame_id = "camera_link"
p.vector.x = 1
p.vector.y = 0.6

ROS.tf.transform(bf, p, "gps_link", ROS.Time(0), "base_link", ROS.Duration(2))
ROS.tf.transform(bf, p, "base_link", ROS.Duration(2))

q1 = ROS.geometry_msgs_Quaternion()
q1.z = 1
q1.w = 1
q1 = ROS.tf.normalize(q1)
@test q1.z ≈ 0.707 atol=1e-03
q2 = ROS.tf.setRPY(float(π), 0, 0.00001)
@test q2.x ≈ 1 atol=0.01
q3 = q1*q2
q4 = q2+q1
@test (q1-q2).x ≈ -1 atol=0.01
@test (-q1).z == -q1.z
axis = ROS.tf.getAxis(q3)
@test axis[1] ≈ 0.707 atol=0.01
@test axis[2] ≈ 0.707 atol=0.01
@test axis[3] ≈ 0 atol=0.01

@test ROS.tf.angle(q1, q2) ≈ π/2 atol=0.01
@test ROS.tf.angleShortestPath(q1, q2) ≈ π atol=0.01
@test ROS.tf.dot(q1, q2) ≈ 0 atol=0.01
@test ROS.tf.getAngle(q1) ≈ π/2 atol=0.01
@test ROS.tf.getAngleShortestPath(q1) ≈ π/2 atol=0.01
@test ROS.tf.length(q2) ≈ 1 atol=0.01
@test ROS.tf.length2(q3) ≈ 1 atol=0.01
@test ROS.tf.inverse(q1).z == -q1.z

q5 = ROS.tf.setEuler(0, 0, float(π))
@test q5.x ≈ 0 atol=0.01
@test q5.y ≈ 0 atol=0.01
@test q5.z ≈ 1 atol=0.01
@test q5.w ≈ 0 atol=0.01

q6 = ROS.tf.slerp(q1, q5, 0.3)
@test q6.x ≈ 0 atol=0.01
@test q6.y ≈ 0 atol=0.01
@test q6.z ≈ 0.85 atol=0.01
@test q6.w ≈ 0.52 atol=0.01

q7 = ROS.tf.setRotation(q1, [1,0,0], float(π/2))
@test q7.x ≈ 0.707 atol=1e-03
@test q7.y ≈ 0 atol=0.01
@test q7.z ≈ 0 atol=0.01
@test q7.w ≈ 0.707 atol=1e-03

println("All $(basename(@__FILE__)) tests passed.")


text:
f(u,p,t) = p.*u

using DifferentialEquations, Plots
u0 = [10.0]
p = [-0.3]
tspan = (0.0,10.0)
prob = ODEProblem(f,u0,tspan,p)
sol = solve(prob)
plot(sol)
ylims!(0.0,10.0)

using Distributions
u0_dist = [Uniform(-10.0,10.0)]

prob_func(prob,i,repeat) = remake(prob, u0 = rand.(u0_dist))
ensemble_prob = EnsembleProblem(prob,prob_func=prob_func)

ensemblesol = solve(ensemble_prob,Tsit5(),EnsembleThreads(),trajectories=100000)

plot(ensemblesol, vars = (0,1), lw=1,alpha=0.1, label=nothing, idxs = 1:250)

g(sol) = sol(4.0)
mean([g(sol) for sol in ensemblesol])

using DiffEqUncertainty
expectation(g, prob, u0_dist, p, MonteCarlo(), Tsit5(); trajectories=100000)

expectation(g, prob, u0_dist, p, Koopman(), Tsit5())

exp(p[1]*4.0)*mean(u0_dist[1])

@time expectation(g, prob, u0_dist, p, MonteCarlo(), Tsit5(); trajectories=100000)

@time expectation(g, prob, u0_dist, p, Koopman(), Tsit5())

u0_dist = [Uniform(0.0,10.0)]
@time expectation(g, prob, u0_dist, p, MonteCarlo(), Tsit5(); trajectories=100_000)

@time expectation(g, prob, u0_dist, p, Koopman(), Tsit5())[1]

exp(p[1]*4.0)*mean(u0_dist[1])

u0_dist = [Normal(3.0,2.0)]
expectation(g, prob, u0_dist, p, Koopman(), Tsit5())

exp(p[1]*4.0)*mean(u0_dist[1])

u0_dist = [truncated(Normal(3.0,2.0),-1000,1000)]
expectation(g, prob, u0_dist, p, Koopman(), Tsit5())

u0_dist = [truncated(Normal(3.0,2.0),-5,11)]
expectation(g, prob, u0_dist, p, Koopman(), Tsit5())

u0_dist = [truncated(Normal(3.0,2.0),3-1000,3+1000)]
expectation(g, prob, u0_dist, p, Koopman(), Tsit5())

g(sol) = [sol(4.0)[1], sol(6.0)[1]]
expectation(g, prob, u0_dist, p, Koopman(), Tsit5(); nout = 2)

exp.(p.*[4.0,6.0])*mean(u0_dist[1])

saveat = tspan[1]:.5:tspan[2]
g(sol) = Matrix(sol)
mean_koop = expectation(g, prob, u0_dist, p, Koopman(), Tsit5(); nout = length(saveat), saveat=saveat)

plot(t->exp(p[1]*t)*mean(u0_dist[1]),tspan..., xlabel="t", label="analytical")
scatter!(collect(saveat),mean_koop.u[:],marker=:o, label=nothing)

function g(sol, power, counter)
    counter[power] = counter[power] + 1
    sol(4.0)[1]^power
end

counters = [0,0,0]
x_koop = expectation(s->g(s,1,counters), prob, u0_dist, p, Koopman(), Tsit5())
x2_koop = expectation(s->g(s,2,counters), prob, u0_dist, p, Koopman(), Tsit5())
x3_koop = expectation(s->g(s,3,counters), prob, u0_dist, p, Koopman(), Tsit5())
counters

function g(sol, counter)
    counter[1] = counter[1] + 1
    v = sol(4.0)[1]
    [v, v^2, v^3]
end

counter = [0]
expectation(s->g(s,counter), prob, u0_dist, p, Koopman(), Tsit5(); nout = 3)
counter

function g(sol)
    x = sol(4.0)[1]
    [x, x^2]
end

koop = expectation(g, prob, u0_dist, p, Koopman(), Tsit5(); nout = 2)
mean_koop = koop[1]
var_koop = koop[2] - mean_koop^2

exp(2*p[1]*4.0)*var(u0_dist[1])

saveat = tspan[1]:.5:tspan[2]
g(sol) = [Matrix(sol)'; (Matrix(sol).^2)']

koop = expectation(g, prob, u0_dist, p, Koopman(), Tsit5(); nout = length(saveat)*2, saveat=saveat)
μ = koop.u[1:length(saveat)]
σ = sqrt.(koop.u[length(saveat)+1:end] - μ.^2)

plot(t->exp(p[1]*t)*mean(u0_dist[1]),tspan..., ribbon = t->-sqrt(exp(2*p[1]*t)*var(u0_dist[1])), label="Analytical Mean, 1 std bounds")
scatter!(collect(saveat),μ,marker=:x, yerror = σ, c=:black, label = "Koopman Mean, 1 std bounds")

function g(sol)
    v = sol(4.0)[1]
    [v, v^2, v^3]
end

koop = expectation(g, prob, u0_dist, p, Koopman(), Tsit5(); nout = 3)
mean_koop = koop[1]
var_koop = koop[2] - mean_koop^2
(koop[3] - 3.0*mean_koop*var_koop - mean_koop^3) / var_koop^(3/2)

g(sol) = sol(4.0)[1]
centralmoment(5, g, prob, u0_dist, p, Koopman(), Tsit5(),
                ireltol = 1e-9, iabstol = 1e-9)

using Quadrature, Cuba

u0_dist = [truncated(Normal(3.0,2.0),-5,11)]
p_dist = [truncated(Normal(-.7, .1), -1,0)]

g(sol) = sol(6.0)[1]

expectation(g, prob, u0_dist, p_dist, Koopman(), Tsit5(), EnsembleThreads();
                quadalg = CubaSUAVE(), batch=10)[1]

using BenchmarkTools

@btime expectation(g, prob, u0_dist, p_dist, Koopman(), Tsit5();
                quadalg = CubaSUAVE())[1]

@btime expectation(g, prob, u0_dist, p_dist, Koopman(), Tsit5(), EnsembleThreads();
                quadalg = CubaSUAVE(), batch=10)[1]

using DiffEqGPU

function f(du, u,p,t)
    @inbounds begin
        du[1] = p[1]*u[1];
    end
    nothing
end

u0 = Float32[10.0]
p = Float32[-0.3]
tspan = (0.0f0,10.0f0)
prob = ODEProblem(f,u0,tspan,p)

g(sol) = sol(6.0)[1]

u0_dist = [truncated(Normal(3.0f0,2.0f0),-5f0,11f0)]
p_dist = [truncated(Normal(-.7f0, .1f0), -1f0,0f0)]

@btime expectation(g, prob, u0_dist, p_dist, Koopman(), Tsit5(), EnsembleGPUArray();
                   quadalg = CubaSUAVE(), batch=1000)[1]

using SciMLTutorials
SciMLTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
module TestGaussianProcesses

using MLJBase
using RDatasets
using Test
using Random:seed!
import CategoricalArrays

seed!(113355)

data = dataset("MASS", "crabs")
X = MLJBase.selectcols(data, [:FL, :RW, :CL, :CW, :BD])
y = MLJBase.selectcols(data, :Sp)

# load code to be tested:
import MLJModels
import GaussianProcesses
using MLJModels.GaussianProcesses_

baregp = GPClassifier()

# split the rows:
allrows = eachindex(y)
train, test = partition(allrows, 0.7, shuffle=true)

fitresult, cache, report =
    MLJBase.fit(baregp, 1, MLJBase.selectrows(X, train), y[train])
yhat = predict(baregp, fitresult, MLJBase.selectrows(X, test));

@test sum(yhat .== y[test]) / length(y[test]) >= 0.7 # around 0.7

fitresult, cache, report = MLJBase.fit(baregp, 1, X, y)
yhat2 = predict(baregp, fitresult, MLJBase.selectrows(X, test));

# gp = machine(baregp, X, y)
# fit!(gp)
# yhat2 = predict(gp, MLJBase.selectrows(X, test))

@test sum(yhat2 .== y[test]) / length(y[test]) >= 0.7

info_dict(baregp)

end # module
true


text:
# # Эксперимент №2. Удачный пример.

# ## Набор параметров
α       = 0.004;        # Параметр регуляризации
w       = 0.0001;       # Эмпирический параметр регуляризации
S       = 1000;         # Количество итераций
β       = 0.0001;       # Шаг минимизации
#
x_tp    = 0.02;         # Стартовое местоположение фронта
T_end   = 0.45;         # Регулируем конечное местоположение фронта
ε       = 10^1.5;       # Крутизна фронта
Nx      = 500;          # Число интервалов по ``X``
Mt      = 1000;        j# Число интервалов по ``T``

# -----------------------------------------------------------------------------
using NonLinearReactionAdvectionDiffusionWithFrontData
using NonLinearReactionAdvectionDiffusionWithFrontData: heterogeneity_map;
using Serialization;
using Plots; pyplot();
# -----------------------------------------------------------------------------

# ## Решение на точных данных

# ### Решение прямой задачи для генерирования экспериментальной информации ----------------------------------------------------------------------------- a, b, t₀, T, N, M, ε, Xₙ, Tₘ, qₙ, ulₘ, urₘ, u₀ = dparams(x_tp = x_tp,
# -----------------------------------------------------------------------------
a, b, t₀, T, N, M, ε, Xₙ, Tₘ, qₙ, ulₘ, urₘ, u₀ = dparams(x_tp = x_tp,
                                                         ε = ε,
                                                         Nx = Nx,
                                                         Mt = Mt,
                                                         T_end = T_end);
u, XX, TP = solve(u₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ);
ϕl, ϕr, ϕ, f1_data, f2_data = generate_obs_data(u, Xₙ, N, Tₘ, M, qₙ, ulₘ, urₘ);
directP = draft(u, Xₙ, N, Tₘ, M, title = "Эскиз прямого решения")
savefig(directP, "direct1.png");
# -----------------------------------------------------------------------------

# ### Найдем начальное приближение
# -----------------------------------------------------------------------------
q_guess = initial_guess(f1_data, Xₙ, N, Tₘ, M, ulₘ, urₘ, 0.005);
initial_guessP = plot(Xₙ, qₙ, label="Истинное");
initial_guessP = plot!(Xₙ, q_guess, label="Найденное")
savefig("initial_guess.png");
# -----------------------------------------------------------------------------

# -----------------------------------------------------------------------------
# В следующей строке, мы записываем все параметры в latex формате,
# не забываем экранировать все спецсимволы.
# Дальше, мы отобразим все эти параметры на графике.
using Printf;
params = "\$w = $(w), \\varepsilon = $(@sprintf("%.4f", ε)), "*
"\\beta = $(β), f_1 \\in $(@sprintf("[%.2f, %.2f]", extrema(f1_data)...)), "*
"\\alpha = $(α)\$"
nothing #hide
# -----------------------------------------------------------------------------

# ### Старт с найденного приближения на точных данных
# -----------------------------------------------------------------------------
q₀ = q_guess;
@time qs, Js, Qs = minimize(q₀, u₀, ulₘ, urₘ, Xₙ, N, Tₘ, M, ε, f1_data, f2_data,
                            S = S, β = β, w = w, showProgress = true)
serialize("qs.jld", qs);
serialize("Js.jld", Js);
serialize("Qs.jld", Qs);

# ### Эскиз процесса минимизации
# -----------------------------------------------------------------------------
a, b, c = minimization_draft(qₙ, Qs, Xₙ, N, Js, zoom = true,
                              annotate_string = params, zoom_chunk = 17//18)
plot(a, size = (800, 800))
plot(b)
plot(c)
withguessP = plot(a, b, c, layout = (3,1), size = (800, 2400));
savefig(withguessP, "withguess1.png")
nothing; #hide
# -----------------------------------------------------------------------------


text:
struct BDay <: DatePeriod
    value::Int64
    calendar::Union{Symbol, String}
    BDay(v::Number, cal::Union{Symbol, String}) = new(v, cal)
end

Base.:(+)(dt::Date, z::BDay) = advancebdays(z.calendar, dt, z.value)
Base.:(-)(dt::Date, z::BDay) = advancebdays(z.calendar, dt, -1 * z.value)


text:
"""
    downsample(times, frequency)

Downsample time steps.  Keep time step at every other `days` or more.
"""
function downsample(times, days)
    # Uniform time steps at every other `days`
    down_times = collect(times[1]:days:times[end])

    # Ensure last value is included
    if down_times[end] != times[end]
        push!(down_times, times[end])
    end

    # Remove values that are not in times
    down_times = intersect(down_times, times)

    return down_times
end

"""
    find_integer_times(times)

Find time steps that are integer.  Return the indices and values found
as `(inds,vals)` where `inds` are the indices of `vals` in `times`.  Throw exception if there are duplicated integers in
`times`.  The unit of `times` in reservoir simulation is usually in days.
"""
function find_integer_times(times)
    integer_time_inds = findall(x->trunc(x) == x, times)
    integer_times = times[integer_time_inds]

    # Check for uniqueness
    if length(integer_times) != length(unique(integer_times))
        # This is caused by small time steps, and simulator
        # does not print all digits.
        throw(ErrorException("Non-unique integer times detected"))
    end

    return (integer_time_inds, integer_times)
end

"""
    find_indices(filter, values)

Find the indices of `filter` in `values`.  Throw exception if there is
duplication in `values`.
"""
function find_indices(filter, values)
    if length(values) != length(unique(values))
        # This is caused by small time steps, and simulator
        # does not print all digits.
        throw(ErrorException("Non-unique values detected"))
    end

    inds = findall(x-> x in filter, values)

    return inds
end

"""
    find_duplication(times)

Find time steps that are duplicated.  Duplication is a result of low-precision
output from simulator.  The `times` vector is sorted.

# Examples
```julia-repl
julia> inds, vals = find_duplication([1, 1.5, 2, 2, 2, 2, 3])
([4, 5, 6], [2.0, 2.0, 2.0])
```
"""
function find_duplication(times)
    duplicated_times = Vector{eltype(times)}() # == eltype(times)[]
    duplicated_time_inds = Int64[]
    unique_times = Set{eltype(times)}()

    for (ind,time) in enumerate(times)
        if time in unique_times
            push!(duplicated_time_inds, ind)
            push!(duplicated_times, time)
        else
            push!(unique_times, time)
        end
    end

    return duplicated_time_inds, duplicated_times
end

"""

Repair
`[1, 1.5, 2, 2, 2, 2, 3]` to [1, 1.5, 2, 2.25, 2.50, 2.75, 3].
Notice the last time step cannot be duplicated.
"""
function repair_times(times)
    # Last time step cannot be duplicated
    if length(times) > 1 && times[end-1] == times[end]
        throw(ErrorException("Last time step is duplicated"))
    end

    # eltype(times) -> Float64 for example
    ret = eltype(times)[]

    ind = 1
    while (ind <= length(times))
        push!(ret, times[ind])

        ind_end = find_last_duplication(ind, times)

        if ind_end == ind # no duplication
            ind += 1
            continue
        end

        # Duplication detected
        # Add delta to all the duplicated value except the first one

        # Example
        # [1, 1.5, 2, 2, 2, 2, 3]
        # ind     = 3
        # ind_end = 6
        # delta = (3-2)/(6+1-3) = 0.25
        delta = (times[ind_end+1] - times[ind])/(ind_end+1-ind)
        for i in ind:ind_end-1
            push!(ret, ret[i] + delta)
        end

        ind = ind_end + 1
    end

    return ret
end

"""

Find the index to the last duplication in the ordered list `values`,
starting from `values[ind]`.

Return `6` for `ind=3` and `values=[1, 1.5, 2, 2, 2, 2, 3]`

# Examples
```julia-repl
julia> find_last_duplication(3, [1, 1.5, 2, 2, 2, 2, 3])
6
```
"""
function find_last_duplication(ind, values)
    value = values[ind]
    for i = ind+1:length(values)
        if values[i] != value
            return i-1
        end
    end

    return ind # == length(values)
end

"""

Deprecated
Smooth production data
Backward integration
"""
function smooth_production(df::DataFrame, time_col_name, col_names)
    times = df[!,time_col_name]

    # Find times at integer
    # TODO: User array operation
    int_times = Int64[] # values of times where time is int
    ind_times = Int64[] # index into times where time is int
    for (ind,time) in enumerate(times)
       if trunc(time) == time
            push!(int_times, time)
            push!(ind_times, ind)
       end
    end

    count = length(ind_times)
    rdf = DataFrame(
        "TIME" => int_times
    )
    #@show ind_times
    for col_name in col_names # each column
        col = df[!,col_name]

        new_col = Float64[]
        push!(new_col, col[1])
        # For each integration interval, i.e., time = 0-1, 1-2, 2-3, 3-4 ...
        for (first_ind, last_ind, first_time, last_time) in zip(ind_times[1:end-1], ind_times[2:end], int_times[1:end-1], int_times[2:end])

            #@show first_ind, last_ind
            total = 0.0
            # index into the data within an interval,
            # backward integration
            for ind = first_ind+1:last_ind
                total += col[ind] * (times[ind] - times[ind-1])
            end
            total /= (last_time - first_time)
            push!(new_col, total)
        end
        rdf[:,col_name] = new_col
        #@show rdf[!,col_name]
    end

    return rdf
end
"""

# Arguments
- `selected_times`: Smoothed to these selected times
- `time_col`: Name of time column
- `cols`: Names of the columns to be smoothed

"""
function smooth_production(df::DataFrame, times, selected_times, cols)
    #times = df[!, time_col]
    #@show times
    selected_indices = find_indices(selected_times, times) #findall(x->x in selected_times, times)
    #@show selected_times
    if length(selected_times) != length(selected_indices)
        #@show length(selected_indices)
        #@show length(selected_times)
        #@show selected_indices
        for t in selected_indices

                #@show times[t]

        end
    end
    @assert length(selected_times) == length(selected_indices) "Not all selected times found"

    selected_count = length(selected_times) # No. of selected times
    # DataFrame for selected values
    sdf = DataFrame("TIME" => selected_times)

    for col_name in cols # column names
        col = df[!, col_name]
        new_col = Float64[]
        push!(new_col, col[1]) # backward time difference so 1st is the same
        for i in 1:selected_count-1 # Each interval between selected times
            ind1  = selected_indices[i]
            time1 =   selected_times[i]
            ind2  = selected_indices[i+1]
            time2 =   selected_times[i+1]

            # if time1 == time2
            #     @show time1
            # end
            # Cal average
            # This method is mass conservative.
            total = 0.0
            #@show time1
            for ind = ind1+1 : ind2 # Each interval between a pair of selected times
                # backward finite-difference in time
                total += col[ind] * (times[ind] - times[ind-1])
            end
            total /= (time2-time1)
            if isnan(total)
                @show total
            end
            push!(new_col, total)
        end

        sdf[:,col_name] = new_col
    end

    return sdf
end


text:
function Base.show(io::IO,hd::VTUHeader{T}) where {T<:Union{UInt32,UInt64}}
	print(io,"VTUHeader{$T}(")
	print(io,Int(hd.num_blocks),", ")
	print(io,Int(hd.blocksize),", ")
	print(io,Int(hd.last_blocksize),", ")
	print(io,"$T[")
	for (i,n) in enumerate(hd.compressed_blocksizes)
		print(io,Int(n))
		i==hd.num_blocks ? nothing : print(io,", ")
	end
	print(io,"])")
end
raw_headerlength(T::Union{Type{UInt32},Type{UInt64}},nblocks) = length(Base64.base64encode(rand(UInt8,3*sizeof(T)+nblocks*sizeof(T))))
headerlength(T::Union{Type{UInt32},Type{UInt64}},nblocks) = 3*sizeof(T)+nblocks*sizeof(T)
bytes(header::VTUHeader) = reinterpret(UInt8,vcat(header.num_blocks,header.blocksize,header.last_blocksize,header.compressed_blocksizes...))


text:
########################################################
## alternative convention: counter clockwise rotation ##
########################################################

## type declarations
##------------------

immutable Rot0CCW <: counterCWRotation
    description::String
    degrees::Int

    function Rot0CCW(descr::String, degr::Int)
        if degr != 0
            error("Wrong degrees for zero rotation!")
        end
        new(descr, degr)
    end

    Rot0CCW() = Rot0CW("No rotation", 0)
end

immutable Rot90CCW <: counterCWRotation
    description::String
    degrees::Int

    function Rot90CCW(descr::String, degr::Int)
        if degr != 90
            error("Wrong degrees for 90 degrees rotation!")
        end
        new(descr, degr)
    end

    Rot90CCW() = Rot90CCW("Counter clockwise rotation of 90 degrees",
                          90)
end

immutable Rot180CCW <: counterCWRotation
    description::String
    degrees::Int

    function Rot180CCW(descr::String, degr::Int)
        if (degr != 180)
            error("Wrong degrees for 180 degrees rotation!")
        end
        new(descr, degr)
    end

    errMsg = "Counter clockwise rotation of 180 degrees"
    Rot180CCW() =  Rot180CCW(errMsg, 180)

end

immutable Rot270CCW <: counterCWRotation
    description::String
    degrees::Int

    function Rot270CCW(descr::String, degr::Int)
        if (degr != 270)
            error("Wrong degrees for 270 degrees rotation!")
        end
        new(descr, degr)
    end

    errMsg = "Clockwise rotation of 270 degrees"
    Rot270CCW() = Rot270CCW(errMsg, 270)
end

#######################################################
## resolving pdf and cdf arguments for modifications ##
#######################################################

## counter clockwise rotation, array of points
##--------------------------------------------

function resolve_args(rot::Rot0CCW,
                      u1::FloatVec, u2::FloatVec)
    return (u1, u2)
end

function resolve_args(rot::Rot90CCW,
                      u1::FloatVec, u2::FloatVec)
    return (u2, 1-u1)
end

function resolve_args(rot::Rot180CCW,
                      u1::FloatVec, u2::FloatVec)
    return (1-u1, 1-u2)
end

function resolve_args(rot::Rot270CCW,
                      u1::FloatVec, u2::FloatVec)
    return (1-u2, u1)
end

##########################
## function evaluations ##
##########################

## pdf
##----

function pdf(cop::PairCop, rot::counterCWRotation,
             u1::FloatVec, u2::FloatVec)
    v1, v2 = resolve_args(rot, u1, u2)
    return pdf(cop, v1, v2)
end

## cdf
##----

function cdf(cop::PairCop, rot::Rot0CCW,
             u1::FloatVec, u2::FloatVec)
    return cdf(cop, u1, u2)
end

function cdf(cop::PairCop, rot::Rot90CCW,
             u1::FloatVec, u2::FloatVec)
    return u2 - cdf(cop, u2, 1-u1)
end

function cdf(cop::PairCop, rot::Rot180CCW,
             u1::FloatVec, u2::FloatVec)
    return u1 + u2 + cdf(cop, 1-u1, 1-u2) - 1
end

function cdf(cop::PairCop, rot::Rot270CCW,
             u1::FloatVec, u2::FloatVec)
    return u1 - cdf(cop, 1-u2, u1)
end

## h-functions
##------------

function hfun(cop::PairCop, rot::Rot0CCW,
              u1::FloatVec, u2::FloatVec)
    return hfun(cop, u1, u2)
end

function hfun(cop::PairCop, rot::Rot90CCW,
              u1::FloatVec, u2::FloatVec)
    return 1 - vfun(cop, 1-u1, u2)
end

function hfun(cop::PairCop, rot::Rot180CCW,
              u1::FloatVec, u2::FloatVec)
    return 1 - hfun(cop, 1-u1, 1-u2)
end

function hfun(cop::PairCop, rot::Rot270CCW,
              u1::FloatVec, u2::FloatVec)
    return vfun(cop, u1, 1-u2)
end

## v-functions
##------------

function vfun(cop::PairCop, rot::Rot0CCW,
              u2::FloatVec, u1::FloatVec)
    return vfun(cop, u2, u1)
end

function vfun(cop::PairCop, rot::Rot90CCW,
              u2::FloatVec, u1::FloatVec)
    return hfun(cop, u2, 1-u1)
end

function vfun(cop::PairCop, rot::Rot180CCW,
              u2::FloatVec, u1::FloatVec)
    return 1 - vfun(cop, 1-u2, 1-u1)
end

function vfun(cop::PairCop, rot::Rot270CCW,
              u2::FloatVec, u1::FloatVec)
    return 1 - hfun(cop, 1-u2, u1)
end

## inverse h-functions
##--------------------

function hinv(cop::PairCop, rot::Rot0CCW,
              u1::FloatVec, u2::FloatVec)
    return hinv(cop, u1, u2)
end

function hinv(cop::PairCop, rot::Rot90CCW,
              u1::FloatVec, u2::FloatVec)
    return 1 - vinv(cop, 1-u1, u2)
end

function hinv(cop::PairCop, rot::Rot180CCW,
              u1::FloatVec, u2::FloatVec)
    return 1 - hinv(cop, 1-u1, 1-u2)
end

function hinv(cop::PairCop, rot::Rot270CCW,
              u1::FloatVec, u2::FloatVec)
    return vinv(cop, u1, 1-u2)
end

## inverse v-functions
##--------------------

function vinv(cop::PairCop, rot::Rot0CCW,
              u2::FloatVec, u1::FloatVec)
    return vinv(cop, u2, u1)
end

function vinv(cop::PairCop, rot::Rot90CCW,
              u2::FloatVec, u1::FloatVec)
    return hinv(cop, u2, 1-u1)
end

function vinv(cop::PairCop, rot::Rot180CCW,
              u2::FloatVec, u1::FloatVec)
    return 1 - vinv(cop, 1-u2, 1-u1)
end

function vinv(cop::PairCop, rot::Rot270CCW,
              u2::FloatVec, u1::FloatVec)
    return 1 - hinv(cop, 1-u2, u1)
end


text:
import PyCall
import PyPlot

plt = PyPlot

"""
    plot!(fitter::Fitter; kwargs...)

    (fitter::Fitter) |> plot!

Plots the data and fitting functions associated with `fitter`. Updates the
settings of `fitter` from the values given in `kwargs` before fitting.
Returns the canvas.
"""
function plot!(fitter::Fitter; kwargs...)
    set!(fitter; kwargs...)

    figure_number = fitter._figure_number
    if figure_number < 0 || figure_number ∉ plt.get_fignums()
        figure_number = 0
        while figure_number ∈ plt.get_fignums()
            figure_number += 1
        end
        fitter._figure_number = figure_number
    end

    fig = plt.figure(figure_number)

    plot(fitter, fig)
end

"""
    plot(fitter::Fitter)

    (fitter::Fitter) |> plot

Plots the data and fitting functions associated with `fitter`. Returns
the canvas.

# Caveats

Using `plot!` is recommended unless minimizing side effects is required
by a specific application
"""
function plot(fitter::Fitter, fig)
    plt.figure(fig[:number])

    ax_main = plt.subplot2grid((4, 1), (0, 0), rowspan=3)
    plt.xscale(fitter[:xscale])
    plt.yscale(fitter[:yscale])

    ax_resid = plt.subplot2grid((4, 1), (3, 0), rowspan=1, sharex=ax_main)
    plt.xscale(fitter[:xscale])

    xmin, xmax = xlims(fitter)
    xextra = 0.1 * (xmax - xmin)

    xmin -= xextra
    xmax += xextra

    ax_main[:set_xlim](xmin, xmax)

    ax_resid[:set_xlabel](fitter[:xlabel])
    ax_resid[:set_ylabel]("Studentized residuals")
    ax_main[:set_ylabel](fitter[:ylabel])
    plt.setp(ax_main[:get_xticklabels](), visible=false)

    fit_mask = data_mask(fitter)

    xdata_included = xdata(fitter)[fit_mask]
    ydata_included = ydata(fitter)[fit_mask]
    eydata_included = eydata(fitter)[fit_mask]

    ax_main[:errorbar](
        xdata_included, ydata_included, eydata_included;
        fitter[:style_data]...)

    inverted_mask = (@__dot__ ~fit_mask)

    xdata_excluded = xdata(fitter)[inverted_mask]
    ydata_excluded = ydata(fitter)[inverted_mask]
    eydata_excluded = eydata(fitter)[inverted_mask]

    ax_main[:errorbar](
        xdata_excluded, ydata_excluded, eydata_excluded;
        fitter[:style_outliers]...)

    if fitter[:plot_curve] || fitter[:plot_guess]
        x_plot = linspace(xmin, xmax, fitter[:fpoints])
    end

    if fitter[:plot_guess]
        y_guess = apply_f(fitter, x_plot, guesses(fitter))
        ax_main[:plot](x_plot, y_guess; fitter[:style_guess]...)
    end

    try
        residuals = studentized_residuals(fitter)

        try
            ax_resid[:errorbar](
                xdata_included, residuals,
                ones(xdata_included); fitter[:style_data]...)
        catch e
            if isa(e, PyCall.PyError)
                msg = """
                    Length of xdata does not match length of residuals.
                    `fit!` must be called after calling `apply_mask!` or
                    `ignore_residuals!` in order to update residuals."""
                warn(msg)
            else
                rethrow(e)
            end
        end

        if fitter[:plot_curve]
            y_curve = apply_f(fitter, x_plot)
            ax_main[:plot](x_plot, y_curve; fitter[:style_fit]...)
            ax_resid[:plot]([xmin, xmax], [0, 0]; fitter[:style_fit]...)
        end
    catch e
        if isa(e, NoResultsException)
            warn(e.msg)
        else
            rethrow(e)
        end
    end

    fig
end

plot(fitter::Fitter) = plot(fitter, plt.figure())


text:
################################################################################
#
#             acb_root_ctx.jl: handling roots using arb
#
# This file is part of hecke.
#
# Copyright (c) 2015: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015 Tommy Hofmann
#
################################################################################
#
#  TODO:
#  - Documentation
#
################################################################################

function BigFloat(x::arb)
  n = ccall((:arb_bits, :libarb), Int, (Ref{arb}, ), x)
  p = precision(BigFloat)
  setprecision(BigFloat, max(n+10,p))
  a = BigFloat()
  b = BigFloat()
  ccall((:arb_get_interval_mpfr, :libarb), Nothing, (Ref{BigFloat}, Ref{BigFloat}, Ref{arb}), a, b, x)
  c = (a+b)/2
  setprecision(BigFloat, p)
  return c

end

################################################################################
#
#  Field access
#
################################################################################

prec(x::acb_root_ctx) = x.prec

poly(x::acb_root_ctx) = x.poly

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, x::acb_root_ctx)
  print(io, "acb_root_ctx of $(poly(x)) with precision $(prec(x))")
end

################################################################################
#
#  Refining
#
################################################################################

# This will refine the roots to given target_prec
# If none is given, double the precision
function refine(x::acb_root_ctx, target_prec::Int = 2*prec(x))

  if target_prec < prec(x)
    return nothing
  end

  x.roots = _refine_roots!(x.poly, x._roots, target_prec, target_prec)
  x.prec = target_prec

  r, s = x.signature

  for i in 1:r
    x.real_roots[i] = real(x.roots[i])
  end

  j = 0
  for i in r+1:degree(x.poly)
    if ispositive(imag(x.roots[i]))
      j += 1
      x.complex_roots[j] = x.roots[i]
    end
  end
  @assert j == s
  nothing
end

function _evaluate(x::fmpq_poly, y::acb)
  z = parent(y)(0)
  for i in 0:degree(x)
    z = z + coeff(x,i)*y^i
  end
  return z
end

################################################################################
#
#  Complex root isolation for fmpz_poly/fmpq_poly
#
################################################################################

# Assume that roots points to an array of deg many acb's.
# This function checks if the radius of the real and imaginary parts
# are smaller then 2^(-abs_tol).
function _validate_size_of_zeros(roots::Ptr{acb_struct}, deg::Int, abs_tol::Int)
  ok = true
  if abs_tol < 0
    return true
  end

  for i = 0 : deg-1
    re = ccall((:acb_real_ptr, :libarb), Ptr{Nemo.arb_struct},
          (Ptr{acb}, ), roots + i * sizeof(acb_struct))
    im = ccall((:acb_imag_ptr, :libarb), Ptr{Nemo.arb_struct},
          (Ptr{acb}, ), roots + i * sizeof(acb_struct))
    t = ccall((:arb_rad_ptr, :libarb), Ptr{Nemo.mag_struct}, (Ptr{arb}, ), re)
    u = ccall((:arb_rad_ptr, :libarb), Ptr{Nemo.mag_struct}, (Ptr{arb}, ), im)
    ok = ok && (ccall((:mag_cmp_2exp_si, :libarb), Cint,
        (Ptr{Nemo.mag_struct}, Int), t, -abs_tol) <= 0)
    ok = ok && (ccall((:mag_cmp_2exp_si, :libarb), Cint,
        (Ptr{Nemo.mag_struct}, Int), u, -abs_tol) <= 0)
    if !ok
      break
    end
  end
  if !ok
    return false
  end
  return true
end

# Return the roots of x with radii < 2^(-abs_tol) as Vector{acb}.
# It is assumed that x is squarefree
function _roots(x::Union{fmpq_poly, fmpz_poly}, abs_tol::Int = 32, initial_prec::Int = 0, max_iter = 0::Int)
  d = degree(x)
  roots = acb_vec(d)

  wp = _roots!(roots, x, abs_tol, initial_prec, max_iter, false)

  res = array(AcbField(wp, false), roots, d)
  acb_vec_clear(roots, d)
  return res
end

# Assume that roots points to approximations of the roots of x.
# This function updates roots inplace to radii <= 2^(-abs_tol) and returns
# the roots as Vector{acb}.
# It is assumed that x is squarefree
function _refine_roots!(x::Union{fmpq_poly, fmpz_poly}, roots::Ptr{acb_struct},
                                                        abs_tol::Int = 32,
                                                        initial_prec::Int = 0,
                                                        max_iter::Int = 0)
  wp = _roots!(roots, x, abs_tol, initial_prec, max_iter, true)
  res = array(AcbField(wp, false), roots, degree(x))
  return res
end

# This is the workhorse.
# It computes the roots of x with with radii <= 2^(-abs_tol)
# The result will be stored in roots
# If have_approx = true, it is assumed that roots contains approximations
# to the roots.
function _roots!(roots::Ptr{acb_struct}, x::Union{fmpq_poly, fmpz_poly},
                                         abs_tol::Int = 32,
                                         initial_prec::Int = 0,
                                         max_iter::Int = 0,
                                         have_approx::Bool = false)
  deg = degree(x)

  initial_prec = (initial_prec >= 2) ? initial_prec : 32

  wp = initial_prec

  while true
    in_roots = roots
    step_max_iter = (max_iter >= 1) ? max_iter : min(max(deg, 32), wp)
    y = acb_poly(x, wp)

    if have_approx
      isolated = ccall((:acb_poly_find_roots, :libarb), Int,
              (Ptr{acb_struct}, Ref{acb_poly}, Ptr{acb_struct}, Int, Int),
              roots, y, in_roots, step_max_iter, wp)
    else
      isolated = ccall((:acb_poly_find_roots, :libarb), Int,
              (Ptr{acb_struct}, Ref{acb_poly}, Ptr{acb_struct}, Int, Int),
              roots, y, C_NULL, step_max_iter, wp)
    end

    have_approx = true

    if isolated == deg
      ok = _validate_size_of_zeros(roots, deg, abs_tol)
      real_ok = ccall((:acb_poly_validate_real_roots, :libarb),
          Bool, (Ptr{acb_struct}, Ref{acb_poly}, Int), roots, y, wp)

      if !real_ok
          ok = false
      else
        for i = 0 : deg - 1
            im = ccall((:acb_imag_ptr, :libarb), Ptr{Nemo.arb_struct},
                (Ptr{acb}, ), roots + i * sizeof(acb_struct))
            if ccall((:arb_contains_zero, :libarb), Bool, (Ptr{Nemo.arb_struct}, ), im)
                ccall((:arb_zero, :libarb), Nothing, (Ptr{Nemo.arb_struct}, ), im)
            end
        end
      end

      if ok
        break
      end
    end

    wp = wp * 2
    if wp > 2^22
      error("Aborting since required precision is > 2^22")
    end
  end

  ccall((:_acb_vec_sort_pretty, :libarb), Nothing,
        (Ptr{acb_struct}, Int), roots, deg)

  return wp
end

function radiuslttwopower(x::arb, e::Int)
  t = ccall((:arb_rad_ptr, :libarb), Ptr{Nemo.mag_struct}, (Ref{arb}, ), x)
  b = ccall((:mag_cmp_2exp_si, :libarb), Cint,
          (Ptr{Nemo.mag_struct}, Int), t, e) <= 0
  return b
end

function radiuslttwopower(x::acb, e::Int)
  re = ccall((:acb_real_ptr, :libarb), Ptr{Nemo.arb_struct},
          (Ref{acb}, ), x)
  im = ccall((:acb_imag_ptr, :libarb), Ptr{Nemo.arb_struct},
          (Ref{acb}, ), x)
  t = ccall((:arb_rad_ptr, :libarb), Ptr{Nemo.mag_struct}, (Ptr{arb}, ), re)
  u = ccall((:arb_rad_ptr, :libarb), Ptr{Nemo.mag_struct}, (Ptr{arb}, ), im)
  ok = (ccall((:mag_cmp_2exp_si, :libarb), Cint,
              (Ptr{Nemo.mag_struct}, Int), t, e) <= 0)

  ok = ok && (ccall((:mag_cmp_2exp_si, :libarb), Cint,
              (Ptr{Nemo.mag_struct}, Int), u, e) <= 0)
  return ok
end

function round(x::arb, p::Int)
  z = ArbField(p, false)()
  ccall((:arb_set_round, :libarb), Nothing, (Ref{Nemo.arb}, Ref{Nemo.arb}, Int), z, x, p)
  return z
end

function round(x::acb, p::Int)
  z = AcbField(p, false)()
  ccall((:acb_set_round, :libarb), Nothing, (Ref{Nemo.acb}, Ref{Nemo.acb}, Int), z, x, p)
  return z
end

function arb_trim(x::arb)
  z = arb()
  ccall((:arb_trim, :libarb), Nothing, (Ref{Nemo.arb}, Ref{Nemo.arb}), z, x)
  z.parent = ArbField(arb_bits(z), false)
  return z
end

function round!(z::arb, x::arb, p::Int)
  ccall((:arb_set_round, :libarb), Nothing, (Ref{Nemo.arb}, Ref{Nemo.arb}, Int), z, x, p)
  z.parent = ArbField(p, false)
  return z
end

function round!(z::acb, x::acb, p::Int)
  ccall((:acb_set_round, :libarb), Nothing, (Ref{Nemo.acb}, Ref{Nemo.acb}, Int), z, x, p)
  z.parent = AcbField(p, false)
  return z
end

function bits(x::arb)
  return ccall((:arb_bits, :libarb), Int, (Ref{Nemo.arb}, ), x)
end

function rel_error_bits(x::arb)
  return ccall((:arb_rel_error_bits, :libarb), Int, (Ref{Nemo.arb}, ), x)
end

function rel_accuracy(x::arb)
  return ccall((:arb_rel_accuracy_bits, :libarb), Int, (Ref{Nemo.arb}, ), x)
end

function set!(z::arb, x::arb)
  ccall((:arb_set, :libarb), Nothing, (Ref{Nemo.arb}, Ref{Nemo.arb}), z, x)
  return z
end

function bits(x::acb)
  return ccall((:acb_bits, :libarb), Int, (Ref{Nemo.acb}, ), x)
end

function rel_error_bits(x::acb)
  return ccall((:acb_rel_error_bits, :libarb), Int, (Ref{Nemo.acb}, ), x)
end

function rel_accuracy(x::acb)
  return ccall((:acb_rel_accuracy_bits, :libarb), Int, (Ref{Nemo.acb}, ), x)
end

function set!(z::acb, x::acb)
  ccall((:acb_set, :libarb), Nothing, (Ref{Nemo.acb}, Ref{Nemo.acb}), z, x)
  return z
end

function expand!(x::Union{arb, acb}, max_radius_2exp::Int)
  if rel_accuracy(x) < 0
    # Radius has less precision then the midpoint
    # Think of 0.100001 +/- 10
    # Reducing the precision of the midpoint won't help.
    return x
  end
  if isexact(x)
    return x
  end
  z = deepcopy(x)
  p = bits(x)
  q = div(p, 2)
  if q < 2
    return x
  end
  y = round(x, q)
  while radiuslttwopower(y, max_radius_2exp) && q > 4
    q = div(q, 2)
    set!(x, y)
    y = round!(y, y, q)
  end
  x.parent = parent_type(typeof(x))(max(2, bits(x)))
  @assert radiuslttwopower(x, max_radius_2exp)
  return x
end

function expand(x::Union{arb, acb}, max_radius_2exp::Int)
  z = deepcopy(x)
  expand!(z, max_radius_2exp)
  return z
end


text:
# ---
# title: 1092. Shortest Common Supersequence
# id: problem1092
# author: Indigo
# date: 2022-03-08
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/shortest-common-supersequence/description/>
# hidden: true
# ---
#
# Given two strings `str1` and `str2`, return the shortest string that has both
# `str1` and `str2` as subsequences.  If multiple answers exist, you may return
# any of them.
#
# _(A string S is a subsequence of string T if deleting some number of
# characters from T (possibly 0, and the characters are chosen _anywhere_ from
# T) results in the string S.)_
#
#
#
# **Example 1:**
#
#
#
#     Input: str1 = "abac", str2 = "cab"
#     Output: "cabac"
#     Explanation:
#     str1 = "abac" is a subsequence of "cabac" because we can delete the first "c".
#     str2 = "cab" is a subsequence of "cabac" because we can delete the last "ac".
#     The answer provided is the shortest such string that satisfies these properties.
#
#
#
#
# **Note:**
#
#   1. `1 <= str1.length, str2.length <= 1000`
#   2. `str1` and `str2` consist of lowercase English letters.
#
#
## @lc code=start
using LeetCode

function shortest_common_supersequence(str1::String, str2::String)
    m, n = length(str1), length(str2)
    dp = OffsetArray(fill(0, m + 1, n + 1), -1, -1)
    for i in 1:m, j in 1:n
        dp[i, j] = max(
            dp[i - 1, j], dp[i, j - 1], dp[i - 1, j - 1] + Int(str1[i] == str2[j])
        )
    end
    res = ""
    coord = CartesianIndex(m, n)
    while coord != CartesianIndex(0, 0)
        c1 = coord - CartesianIndex(1, 0)
        c2 = coord - CartesianIndex(0, 1)
        if c1 ∈ CartesianIndices(dp) && dp[coord] == dp[c1]
            res *= str1[coord[1]]
            coord -= CartesianIndex(1, 0)
        elseif c2 ∈ CartesianIndices(dp) && dp[coord] == dp[c2]
            res *= str2[coord[2]]
            coord -= CartesianIndex(0, 1)
        else
            res *= str1[coord[1]]
            coord -= CartesianIndex(1, 1)
        end
    end
    return reverse(res)
end
## @lc code=end


text:
export LinearRegressor, RidgeRegressor, LassoRegressor, ElasticNetRegressor,
       RobustRegressor, HuberRegressor, QuantileRegressor, LADRegressor,
       LogisticClassifier, MultinomialClassifier

const SymStr = Union{Symbol,String}

include("regressors.jl")
include("classifiers.jl")

const REG_MODELS = (LinearRegressor, RidgeRegressor, LassoRegressor,
                    ElasticNetRegressor, RobustRegressor, HuberRegressor,
                    QuantileRegressor, LADRegressor)
const CLF_MODELS = (LogisticClassifier, MultinomialClassifier)
const ALL_MODELS = (REG_MODELS..., CLF_MODELS...)

#= ==========
   REGRESSORS
   ========== =#

function MMI.fit(m::Union{REG_MODELS...}, verb::Int, X, y)
    Xmatrix = MMI.matrix(X)
    reg     = glr(m)
    solver  = m.solver === nothing ? _solver(reg, size(Xmatrix)) : m.solver
    # get the parameters
    θ = fit(reg, Xmatrix, y; solver=solver)
    # return
    return θ, nothing, NamedTuple{}()
end

MMI.predict(m::Union{REG_MODELS...}, θ, Xnew) = apply_X(MMI.matrix(Xnew), θ)

function MMI.fitted_params(m::Union{REG_MODELS...}, θ)
    m.fit_intercept && return (coefs = θ[1:end-1], intercept = θ[end])
    return (coefs = θ, intercept = nothing)
end

#= ===========
   CLASSIFIERS
   =========== =#

function MMI.fit(m::Union{CLF_MODELS...}, verb::Int, X, y)
    Xmatrix  = MMI.matrix(X)
    yplain   = convert.(Int, MMI.int(y))
    classes  = MMI.classes(y[1])
    nclasses = length(classes)
    if nclasses == 2
        # recode
        yplain[yplain .== 1] .= -1
        yplain[yplain .== 2] .= 1
        c = 1
    else
        c = nclasses
    end
    clf    = glr(m)
    solver = m.solver === nothing ? _solver(clf, size(Xmatrix)) : m.solver
    # get the parameters
    θ = fit(clf, Xmatrix, yplain, solver=solver)
    # return
    return (θ, c, classes), nothing, NamedTuple{}()
end

function MMI.predict(m::Union{CLF_MODELS...}, (θ, c, classes), Xnew)
    Xmatrix = MMI.matrix(Xnew)
    preds   = apply_X(Xmatrix, θ, c)
    # binary classification
    if c == 1
        preds  .= sigmoid.(preds)
        preds   = hcat(1.0 .- preds, preds) # scores for -1 and 1
        return [MMI.UnivariateFinite(classes, preds[i, :]) for i in 1:size(Xmatrix,1)]
    end
    # multiclass
    preds .= softmax(preds)
    return [MMI.UnivariateFinite(classes, preds[i, :]) for i in 1:size(Xmatrix,1)]
end

function MMI.fitted_params(m::Union{CLF_MODELS...}, (θ, c, classes))
    if c > 1
        if m.fit_intercept
            W = reshape(θ, div(length(θ), c), c)
            return (coefs = W, intercept = nothing)
        end
        W = reshape(θ, p+1, c)
        return (coefs = W[1:p, :], intercept = W[end, :])
    end
    # single class
    m.fit_intercept && return (coefs = θ[1:end-1], intercept = θ[end])
    return (coefs = θ, intercept = nothing)
end

#= =======================
   METADATA FOR ALL MODELS
   ======================= =#

MMI.metadata_pkg.(ALL_MODELS,
    name="MLJLinearModels",
    uuid="6ee0df7b-362f-4a72-a706-9e79364fb692",
    url="https://github.com/alan-turing-institute/MLJLinearModels.jl",
    julia=true,
    license="MIT",
    is_wrapper=false)

descr_(M) = descr(M) *
    "\n→ based on [MLJLinearModels](https://github.com/alan-turing-institute/MLJLinearModels.jl)" *
    "\n→ do `@load $(MMI.name(M)) pkg=\"MLJLinearModels\" to use the model.`" *
    "\n→ do `?$(MMI.name(M))` for documentation."
lp_(M) = "MLJLinearModels.$(MMI.name(M))"

for M in REG_MODELS
    MMI.metadata_model(M,
        input=MMI.Table(MMI.Continuous),
        target=AbstractVector{MMI.Continuous},
        weights=false,
        descr=descr_(M), path=lp_(M))
end
for M in CLF_MODELS
    MMI.metadata_model(M,
        input=MMI.Table(MMI.Continuous),
        target=AbstractVector{<:MMI.Finite},
        weights=false,
        descr=descr_(M), path=lp_(M))
end


text:
function getvar()
    println("Predefined vars that can be calculated for each cell/particle:")
    println("----------------------------------------------------------------")
    println("=============================[gas]:=============================")
    println("       -all the non derived hydro vars-")
    println(":cpu, :level, :rho, :cx, :cy, :cz, :vx, :vy, :vz, :p, var6,...")
    println()
    println("              -derived hydro vars-")
    println(":x, :y, :z")
    println(":mass, :cellsize, :volume, :freefall_time")
    println(":cs, :mach, :jeanslength, :jeansnumber")
    println()
    println("==========================[particles]:==========================")
    println("        all the non derived  vars:")
    println(":cpu, :level, :id, :family, :tag ")
    println(":x, :y, :z, :vx, :vy, :vz, :mass, :birth, :metal....")
    println()
    println("              -derived particle vars-")
    println(":age")
    println()
    println("===========================[clumps]:===========================")
    println(":peak_x or :x, :peak_y or :y, :peak_z or :z")
    println(":v, :ekin,...")
    println()
    #println("===========================[sinks]:============================")
    println("=====================[gas or particles]:=======================")
    println(":v, :ekin")
    println()
    println("related to a given center:")
    println("---------------------------")
    println(":r_cylinder, :r_sphere (radial components)")
    #println(":ϕ, :θ")
    println(":vr_cylinder") #, vr_sphere (radial components)")
    println(":vϕ") #", :vθ")
    #println(":l, :lx, :ly, :lz :lr, :lϕ, :lθ (angular momentum)")
    println("----------------------------------------------------------------")
    return
end

"""
#### Get variables or derived quantities from the dataset:
- overview the list of predefined quantities with: getinfo()
- select variable(s) and their unit(s)
- give the spatial center (with units) of the data within the box (relevant e.g. for radius dependency)
- relate the coordinates to a direction (x,y,z)
- pass a modified database
- pass a mask to exclude elements (cells/particles/...) from the calculation

```julia
getvar(   dataobject::DataSetType, var::Symbol;
        filtered_db::JuliaDB.AbstractIndexedTable=JuliaDB.table([1]),
        center::Array{<:Any,1}=[0.,0.,0.],
        center_unit::Symbol=:standard,
        direction::Symbol=:z,
        unit::Symbol=:standard,
        mask::MaskType=[false],
        ref_time::Real=dataobject.info.time)

return Array{Float64,1}

```

#### Arguments
##### Required:
- **`dataobject`:** needs to be of type: "DataSetType"
- **`var(s)`:** select a variable from the database or a predefined quantity (see field: info, function getvar(), dataobject.data)
##### Predefined/Optional Keywords:
- **`filtered_db`:** pass a filtered or manipulated database together with the corresponding DataSetType object (required argument)
- **`center`:** in units given by argument `center_unit`; by default [0., 0., 0.]; the box-center can be selected by e.g. [:bc], [:boxcenter], [value, :bc, :bc], etc..
- **`center_unit`:** :standard (code units), :Mpc, :kpc, :pc, :mpc, :ly, :au , :km, :cm (of typye Symbol) ..etc. ; see for defined length-scales viewfields(info.scale)
- **`direction`:** todo
- **`unit(s)`:** return the variable in given units
- **`mask`:** needs to be of type MaskType which is a supertype of Array{Bool,1} or BitArray{1} with the length of the database (rows)
- **`ref_time`:** reference zero-time for particle age calculation

### Defined Methods - function defined for different arguments
- getvar(   dataobject::DataSetType, var::Symbol; ...) # one given variable -> returns 1d array
- getvar(   dataobject::DataSetType, var::Symbol, unit::Symbol; ...) # one given variable with its unit -> returns 1d array
- getvar(   dataobject::DataSetType, vars::Array{Symbol,1}; ...) # several given variables -> array needed -> returns dictionary with 1d arrays
- getvar(   dataobject::DataSetType, vars::Array{Symbol,1}, units::Array{Symbol,1}; ...) # several given variables and their corresponding units -> both arrays -> returns dictionary with 1d arrays
- getvar(   dataobject::DataSetType, vars::Array{Symbol,1}, unit::Symbol; ...) # several given variables that have the same unit -> array for the variables and a single Symbol for the unit -> returns dictionary with 1d arrays

#### Examples
```julia
# read simulation information
julia> info = getinfo(420)
julia> gas = gethydro(info)

# Example 1: get the mass for each cell of the hydro data (1dim array)
mass1 = getvar(gas, :mass)  # in [code units]
mass = getvar(gas, :mass) * gas.scale.Msol # scale the result from code units to solar masses
mass = getvar(gas, :mass, unit=:Msol) # unit calculation, provided by a keyword argument
mass = getvar(gas, :mass, :Msol) # unit calculation provided by an argument

# Example 2: get the mass and |v| (several variables) for each cell of the hydro data
quantities = getvar(gas, [:mass, :v]) # in [code units]
returns: Dict{Any,Any} with 2 entries:
  :mass => [8.9407e-7, 8.9407e-7, 8.9407e-7, 8.9407e-7, 8.9407e-7, 8.9407e-7, 8…
  :v => [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  2.28274e-7, 2.…

quantities = getvar(gas, [:mass, :v], units=[:Msol, :km_s]) # unit calculation, provided by a keyword argument
quantities = getvar(gas, [:mass, :v], [:Msol, :km_s]) # unit calculation provided by an argument

# Example 3: get several variables in the same units by providing a single argument
quantities = getvar(gas, [:vx, :vy, :vz], :km_s)
...
"""
function getvar(   dataobject::DataSetType, var::Symbol;
                    filtered_db::JuliaDB.AbstractIndexedTable=JuliaDB.table([1]),
                    center::Array{<:Any,1}=[0.,0.,0.],
                    center_unit::Symbol=:standard,
                    direction::Symbol=:z,
                    unit::Symbol=:standard,
                    mask::MaskType=[false],
                    ref_time::Real=dataobject.info.time)

    center = center_in_standardnotation(dataobject.info, center, center_unit)

    # construct corresponding DataSetType from filtered database to use the calculations below
    if typeof(filtered_db) != IndexedTable{StructArrays.StructArray{Tuple{Int64},1,Tuple{Array{Int64,1}},Int64}}
        dataobject = construct_datatype(filtered_db, dataobject);
    end

    return get_data(dataobject, [var], [unit], direction, center, mask, ref_time )
end

function getvar(   dataobject::DataSetType, var::Symbol, unit::Symbol;
                    filtered_db::JuliaDB.AbstractIndexedTable=JuliaDB.table([1]),
                    center::Array{<:Any,1}=[0.,0.,0.],
                    center_unit::Symbol=:standard,
                    direction::Symbol=:z,
                    mask::MaskType=[false],
                    ref_time::Real=dataobject.info.time)

    center = center_in_standardnotation(dataobject.info, center, center_unit)

    # construct corresponding DataSetType from filtered database to use the calculations below
    if typeof(filtered_db) != IndexedTable{StructArrays.StructArray{Tuple{Int64},1,Tuple{Array{Int64,1}},Int64}}
        dataobject = construct_datatype(filtered_db, dataobject);
    end

    return get_data(dataobject, [var], [unit], direction, center, mask, ref_time )
end

function getvar(   dataobject::DataSetType, vars::Array{Symbol,1}, units::Array{Symbol,1};
                    filtered_db::JuliaDB.AbstractIndexedTable=JuliaDB.table([1]),
                    center::Array{<:Any,1}=[0.,0.,0.],
                    center_unit::Symbol=:standard,
                    direction::Symbol=:z,
                    mask::MaskType=[false],
                    ref_time::Real=dataobject.info.time)

    center = center_in_standardnotation(dataobject.info, center, center_unit)

    # construct corresponding DataSetType from filtered database to use the calculations below
    if typeof(filtered_db) != IndexedTable{StructArrays.StructArray{Tuple{Int64},1,Tuple{Array{Int64,1}},Int64}}
        dataobject = construct_datatype(filtered_db, dataobject);
    end

    return get_data(dataobject, vars, units, direction, center, mask, ref_time )
end

function getvar(   dataobject::DataSetType, vars::Array{Symbol,1}, unit::Symbol;
                    filtered_db::JuliaDB.AbstractIndexedTable=JuliaDB.table([1]),
                    center::Array{<:Any,1}=[0.,0.,0.],
                    center_unit::Symbol=:standard,
                    direction::Symbol=:z,
                    mask::MaskType=[false],
                    ref_time::Real=dataobject.info.time)

    center = center_in_standardnotation(dataobject.info, center, center_unit)
    units = fill(unit, length(vars)) # use given unit for all variables

    # construct corresponding DataSetType from filtered database to use the calculations below
    if typeof(filtered_db) != IndexedTable{StructArrays.StructArray{Tuple{Int64},1,Tuple{Array{Int64,1}},Int64}}
        dataobject = construct_datatype(filtered_db, dataobject);
    end

    return get_data(dataobject, vars, units, direction, center, mask, ref_time )
end

function getvar(   dataobject::DataSetType, vars::Array{Symbol,1};
                    filtered_db::JuliaDB.AbstractIndexedTable=JuliaDB.table([1]),
                    center::Array{<:Any,1}=[0.,0.,0.],
                    center_unit::Symbol=:standard,
                    direction::Symbol=:z,
                    units::Array{Symbol,1}=[:standard],
                    mask::MaskType=[false],
                    ref_time::Real=dataobject.info.time)

    center = center_in_standardnotation(dataobject.info, center, center_unit)

    # construct corresponding DataSetType from filtered database to use the calculations below
    if typeof(filtered_db) != IndexedTable{StructArrays.StructArray{Tuple{Int64},1,Tuple{Array{Int64,1}},Int64}}
        dataobject = construct_datatype(filtered_db, dataobject);
    end

    #vars = unique(vars)
    return get_data(dataobject, vars, units, direction, center, mask, ref_time )
end

function center_in_standardnotation(dataobject::InfoType, center::Array{<:Any,1}, center_unit::Symbol)

    # check for :bc, :boxcenter
    Ncenter = length(center)
    if Ncenter  == 1
        if in(:bc, center) || in(:boxcenter, center)
            bc = 0.5
            center = [bc, bc, bc]
        end
    else
        for i = 1:Ncenter
            if center[i] == :bc || center[i] == :boxcenter
                bc = 0.5
                center[i] = bc
            elseif center_unit != :standard
                center[i] = center[i] / dataobject.boxlen .* getunit(dataobject, center_unit)
            end
        end
    end
    return center
end

"""
#### Get mass-array from the dataset (cells/particles/clumps/...):
```julia
getmass(dataobject::HydroDataType)
getmass(dataobject::PartDataType)
getmass(dataobject::ClumpDataType)

return Array{Float64,1}
```
"""
function getmass(dataobject::HydroDataType;)

    lmax = dataobject.lmax
    boxlen = dataobject.boxlen
    isamr = checkuniformgrid(dataobject, lmax)

    #return select(dataobject.data, :rho) .* (dataobject.boxlen ./ 2. ^(select(dataobject.data, :level))).^3
    if isamr
        return select( dataobject.data, (:rho, :level)=>p->p.rho * (boxlen / 2^p.level)^3 )
    else # if uniform grid
        return select( dataobject.data, :rho=>p->p * (boxlen / 2^lmax)^3 )
    end
end

function getmass(dataobject::PartDataType;)
    return getvar(dataobject, :mass)
end

function getmass(dataobject::ClumpDataType;)
    return getvar(dataobject, :mass)
end

"""
#### Get the x,y,z positions from the dataset (cells/particles/clumps/...):
```julia
getpositions( dataobject::DataSetType;
                    unit::Symbol=:standard,
                    direction::Symbol=:z,
                    center::Array{<:Any,1}=[0., 0., 0.],
                    center_unit::Symbol=:standard)

return x, y, z
```

#### Arguments
##### Required:
- **`dataobject`:** needs to be of type: "DataSetType"
##### Predefined/Optional Keywords:
- **`center`:** in unit given by argument `center_unit`; by default [0., 0., 0.]; the box-center can be selected by e.g. [:bc], [:boxcenter], [value, :bc, :bc], etc..
- **`center_unit`:** :standard (code units), :Mpc, :kpc, :pc, :mpc, :ly, :au , :km, :cm (of typye Symbol) ..etc. ; see for defined length-scales viewfields(info.scale)
- **`direction`:** todo
- **`unit`:** return the variables in given unit
- **`mask`:** needs to be of type MaskType which is a supertype of Array{Bool,1} or BitArray{1} with the length of the database (rows)

### Defined Methods - function defined for different arguments

- getpositions( dataobject::DataSetType; ...) # one given variable
- getpositions( dataobject::DataSetType, unit::Symbol; ...) # one given variable with its unit

"""
function getpositions( dataobject::DataSetType, unit::Symbol;
                        direction::Symbol=:z,
                        center::Array{<:Any,1}=[0., 0., 0.],
                        center_unit::Symbol=:standard,
                        mask::MaskType=[false])

    positions = getvar(dataobject, [:x, :y, :z],
                        center=center,
                        center_unit=center_unit,
                        direction=direction,
                        units=[unit, unit, unit],
                        mask=mask)

    return positions[:x], positions[:y], positions[:z]
end

function getpositions( dataobject::DataSetType;
                        unit::Symbol=:standard,
                        direction::Symbol=:z,
                        center::Array{<:Any,1}=[0., 0., 0.],
                        center_unit::Symbol=:standard,
                        mask::MaskType=[false])

    positions = getvar(dataobject, [:x, :y, :z],
                        center=center,
                        center_unit=center_unit,
                        direction=direction,
                        units=[unit, unit, unit],
                        mask=mask)

    return positions[:x], positions[:y], positions[:z]
end

"""
#### Get the extent of the dataset-domain:
```julia
function getextent( dataobject::DataSetType;
                     unit::Symbol=:standard,
                     center::Array{<:Any,1}=[0., 0., 0.],
                     center_unit::Symbol=:standard,
                     direction::Symbol=:z)

return (xmin, xmax), (ymin ,ymax ), (zmin ,zmax )
```

#### Arguments
##### Required:
- **`dataobject`:** needs to be of type: "DataSetType"
##### Predefined/Optional Keywords:
- **`center`:** in unit given by argument `center_unit`; by default [0., 0., 0.]; the box-center can be selected by e.g. [:bc], [:boxcenter], [value, :bc, :bc], etc..
- **`center_unit`:** :standard (code units), :Mpc, :kpc, :pc, :mpc, :ly, :au , :km, :cm (of typye Symbol) ..etc. ; see for defined length-scales viewfields(info.scale)
- **`direction`:** todo
- **`unit`:** return the variables in given unit

### Defined Methods - function defined for different arguments
- getextent( dataobject::DataSetType; # one given variable
- getextent( dataobject::DataSetType, unit::Symbol; ...) # one given variable with its unit

"""
function getextent( dataobject::DataSetType, unit::Symbol;
                     center::Array{<:Any,1}=[0., 0., 0.],
                     center_unit::Symbol=:standard,
                     direction::Symbol=:z)

    return  getextent( dataobject,
                         center=center,
                         center_unit=center_unit,
                         direction=direction,
                         unit=unit)
end

function getextent( dataobject::DataSetType;
                     unit::Symbol=:standard,
                     center::Array{<:Any,1}=[0., 0., 0.],
                     center_unit::Symbol=:standard,
                     direction::Symbol=:z)

    range = dataobject.ranges
    boxlen = dataobject.boxlen

    center = prepboxcenter(dataobject.info, center_unit, center) # code units
    center = center ./ dataobject.boxlen
    #selected_unit = 1.
    # if center_unit != :standard
    #     selected_unit = getunit(dataobject.info, center_unit)
    #     center = center ./ dataobject.boxlen .* selected_unit
    # end

    selected_unit = getunit(dataobject.info, unit)
    xmin = ( range[1] - center[1] ) * boxlen * selected_unit
    xmax = ( range[2] - center[1] ) * boxlen * selected_unit
    ymin = ( range[3] - center[2] ) * boxlen * selected_unit
    ymax = ( range[4] - center[2] ) * boxlen * selected_unit
    zmin = ( range[5] - center[3] ) * boxlen * selected_unit
    zmax = ( range[6] - center[3] ) * boxlen * selected_unit

    if direction == :y
        xmin_buffer = xmin
        xmax_buffer = xmax

        xmin= zmin
        xmax= zmax
        zmin= ymin
        zmax= ymax
        ymin= xmin_buffer
        ymax= xmax_buffer

    elseif direction == :x
        xmin_buffer = xmin
        xmax_buffer = xmax

        xmin= zmin
        xmax= zmax
        zmin= xmin_buffer
        zmax= xmax_buffer

    end

    return (xmin, xmax), (ymin ,ymax ), (zmin ,zmax )
end


text:
using NIDAQ, Base.Test

@testset "installation" begin
@test typeof(getproperties()) == Dict{String,Tuple{Any,Bool}}
@test haskey(getproperties(), "NIDAQMajorVersion")
global dev = devices()
end

@testset "device" begin
if length(dev) == 0
    info("no data acquisition devices found")
    exit()
else
    info("found $(dev).  testing with $(dev[1])")
    dev=dev[1]
    @test typeof(getproperties(dev)) == Dict{String,Tuple{Any,Bool}}
    global props = getproperties(dev)

    @test typeof(getproperties(dev)) == Dict{String,Tuple{Any,Bool}}
    @test typeof(analog_input_channels()) == Vector{String}
    @test typeof(analog_output_channels()) == Vector{String}
    @test typeof(digital_input_channels()) == Vector{String}
    @test typeof(digital_output_channels()) == Vector{String}
    @test typeof(counter_input_channels()) == Vector{String}
    @test typeof(counter_output_channels()) == Vector{String}
    @test typeof(analog_input_ranges()) == Matrix{Float64}
    @test typeof(analog_output_ranges()) == Matrix{Float64}
end
end

@testset "analog input" begin
if length(props["AIPhysicalChans"][1]) == 0
    info("$dev does not support AIPhysicalChans")
else
    t = analog_input(dev*"/ai0")
    @test typeof(getproperties(t)) == Dict{String,Tuple{Any,Bool}}
    @test typeof(getproperties(t,dev*"/ai0")) == Dict{String,Tuple{Any,Bool}}
    @test typeof(t) == NIDAQ.AITask
    @test start(t) == nothing
    @test length(read(t, Float64, 3)) == 3
    @test stop(t) == nothing
    @test analog_input(t, dev*"/ai1") == nothing
    @test start(t) == nothing
    @test length(read(t, UInt32, 6)) == 12
    @test stop(t) == nothing
    @test NIDAQ.CfgSampClkTiming(t.th, convert(Ref{UInt8},b""), 100.0, NIDAQ.Val_Rising,
              NIDAQ.Val_FiniteSamps, UInt64(10)) == 0
    @test start(t) == nothing
    @test length(read(t, Float64)) == 20
    @test stop(t) == nothing
    @test clear(t) == nothing
end
end

@testset "analog output" begin
if length(props["AOPhysicalChans"][1]) == 0
    info("$dev does not support AOPhysicalChans")
else
    t = analog_output(dev*"/ao0")
    @test typeof(t) == NIDAQ.AOTask
    @test start(t) == nothing
    @test write(t, rand(3)) == 3
    @test stop(t) == nothing
    @test analog_output(t, dev*"/ao1") == nothing
    @test start(t) == nothing
    @test write(t, rand(UInt32,6,2)) == 6
    @test stop(t) == nothing
    @test NIDAQ.CfgSampClkTiming(t.th, convert(Ref{UInt8},b""), 100.0, NIDAQ.Val_Rising,
              NIDAQ.Val_FiniteSamps, UInt64(10)) == 0
    @test write(t, rand(UInt32,10,2)) == 10
    @test start(t) == nothing
    @test NIDAQ.WaitUntilTaskDone(t.th,10.0) == 0
    @test stop(t) == nothing
    @test clear(t) == nothing
end
end

@testset "digital input" begin
if length(props["DILines"][1]) == 0
    info("$dev does not support DILines")
else
    t = digital_input(dev*"/Port0/Line0")
    @test typeof(t) == NIDAQ.DITask
    @test start(t) == nothing
    @test length(read(t, UInt32, 3)) == 3
    @test stop(t) == nothing
    @test digital_input(t, dev*"/Port0/Line1") == nothing
    @test start(t) == nothing
    @test length(read(t, UInt8, 6)) == 12
    @test stop(t) == nothing
    rslt = Ref{UInt32}(0)
    NIDAQ.DAQmxGetBufInputOnbrdBufSize(t.th, rslt)
    if rslt[] != 0 #If the device supports buffered digital input
        @test NIDAQ.CfgSampClkTiming(t.th, convert(Ref{UInt8},b""), 100.0, NIDAQ.Val_Rising,
              NIDAQ.Val_FiniteSamps, UInt64(10)) == 0
	@test start(t) == nothing
	@test length(read(t, UInt32)) == 20
	@test stop(t) == nothing
	@test clear(t) == nothing
    else
        info("Device does not support clocked (buffered) digital input")
    end
end
end

@testset "digital output" begin
if length(props["DOLines"][1]) == 0
    info("$dev does not support DOLines")
else
    t = digital_output(dev*"/Port0/Line0")
    @test typeof(t) == NIDAQ.DOTask
    @test start(t) == nothing
    @test write(t, round.(UInt32, [1,0,1,0,1,0])) == 6
    @test stop(t) == nothing
    @test digital_output(t, dev*"/Port0/Line1") == nothing
    @test start(t) == nothing
    @test write(t, round.(UInt8, [1 0; 0 0; 1 0; 0 1; 1 1; 0 1])) == 6
    @test stop(t) == nothing
    rslt = Ref{UInt32}(0)
    NIDAQ.DAQmxGetBufOutputOnbrdBufSize(t.th, rslt)
    if rslt[] != 0 #If the device supports buffered digital output
        @test NIDAQ.CfgSampClkTiming(t.th, convert(Ref{UInt8},b""), 100.0, NIDAQ.Val_Rising,
              NIDAQ.Val_FiniteSamps, UInt64(10)) == 0
        @test write(t, rand(UInt32,10,2)) == 10
        @test start(t) == nothing
        @test NIDAQ.WaitUntilTaskDone(t.th,10.0) == 0
        @test stop(t) == nothing
        @test clear(t) == nothing
    else
        info("Device does not support clocked (buffered) digital output")
    end
end
end

@testset "counter output" begin
if isempty(props["COPhysicalChans"][1][1])
    info("$dev does not support COPhysicalChans")
else
    t = generate_pulses(dev*"/ctr0")
    @test typeof(t) == NIDAQ.COTask
    @test NIDAQ.CfgImplicitTiming(t.th, NIDAQ.Val_FiniteSamps, UInt64(10)) == 0
    @test start(t) == nothing
    @test stop(t) == nothing
    @test clear(t) == nothing
end
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

struct ConceptManager{T <: AbstractCoreTransaction} <: AbstractConceptManager
    transaction::T
end

struct Root
end

# Get root something type. The hard coded strings came from GraqlToken
Base.get(cm::ConceptManager, ::Type{EntityType}, ::Root) = get(cm, ThingType, "entity")
Base.get(cm::ConceptManager, ::Type{RelationType}, ::Root) = get(cm, ThingType, "relation")
Base.get(cm::ConceptManager, ::Type{AttributeType}, ::Root) = get(cm, ThingType, "attribute")
Base.get(cm::ConceptManager, ::Type{ThingType}, ::Root) = get(cm, ThingType, "thing")

function Base.get(cm::ConceptManager, ::Type{EntityType}, label::AbstractString)
    t = get(cm, ThingType, label)
    t !== nothing && t isa EntityType && return t
    return nothing
end

function put(cm::ConceptManager, ::Type{EntityType}, label::AbstractString)
    res = execute(cm, ConceptManagerRequestBuilder.put_entity_type_req(label))
    if which_oneof(res, :res) == :put_entity_type_res
        return instantiate(res.put_entity_type_res.entity_type)
    end
    return nothing
end

function Base.get(cm::ConceptManager, ::Type{RelationType}, label::AbstractString)
    t = get(cm, ThingType, label)
    t !== nothing && t isa RelationType && return t
    return nothing
end

function put(cm::ConceptManager, ::Type{RelationType}, label::AbstractString)
    res = execute(cm, ConceptManagerRequestBuilder.put_relation_type_req(label))
    if which_oneof(res, :res) == :put_relation_type_res
        return instantiate(res.put_relation_type_res.relation_type)
    end
    return nothing
end

function Base.get(cm::ConceptManager, ::Type{AttributeType}, label::AbstractString)
    t = get(cm, ThingType, label)
    t !== nothing && t isa AttributeType && return t
    return nothing
end

function put(cm::ConceptManager, ::Type{AttributeType}, label::AbstractString, value_type::EnumType)
    res = execute(cm, ConceptManagerRequestBuilder.put_attribute_type_req(label, value_type))
    if which_oneof(res, :res) == :put_attribute_type_res
        return instantiate(res.put_attribute_type_res.attribute_type)
    end
    return nothing
end

function Base.get(cm::ConceptManager, ::Type{ThingType}, label::AbstractString)
    req = ConceptManagerRequestBuilder.get_thing_type_req(label)
    res = execute(cm, req)
    if which_oneof(res, :res) == :get_thing_type_res
        if hasproperty(res.get_thing_type_res, :thing_type)
            return instantiate(res.get_thing_type_res.thing_type)
        else
            return nothing
        end
    end
    return nothing
end

function Base.get(cm::ConceptManager, iid::AbstractString)
    res = get_proto_thing(cm, iid)
    res !== nothing && return instantiate(res)
    return nothing
end

function Base.get(cm::ConceptManager, thing::AbstractThing)
    return get(cm, thing.iid)
end

function execute(cm::ConceptManager, req::Proto.Transaction_Req)
    result = execute(cm.transaction, req)
    return result !== nothing ? result.concept_manager_res : result
end

# get_proto_thing is supposed to be a vehicle to get things in its proto form to get
# a result for setting up set_has etc.
function get_proto_thing(cm::ConceptManager, iid::AbstractString)
    req = ConceptManagerRequestBuilder.get_thing_req(iid)
    res = execute(cm, req)
    if which_oneof(res, :res) == :get_thing_res
        if hasproperty(res.get_thing_res, :thing)
            return res.get_thing_res.thing
        else
            return nothing
        end
    end
    return nothing
end


text:
# getindex/view/setindex! ======================================================

#### getindex/view ####

# Integer returns a single value, but not for view
@propagate_inbounds Base.getindex(A::AbstractDimArray, i1::Integer, i2::Integer, I::Integer...) =
    Base.getindex(parent(A), i1, i2, I...)
# No indices. These just prevent stack overflows
@propagate_inbounds Base.getindex(A::AbstractDimArray) = Base.getindex(parent(A))
@propagate_inbounds Base.view(A::AbstractDimArray) = rebuild(A, Base.view(parent(A)), ())

for f in (:getindex, :view, :dotview)
    @eval begin
        #### Array getindex/view ###
        @propagate_inbounds Base.$f(A::AbstractDimArray, i::Integer) = Base.$f(parent(A), i)
        @propagate_inbounds Base.$f(A::AbstractDimArray, I...) =
            Base.$f(A, dims2indices(A, I)...)
        # Linear indexing returns parent type
        @propagate_inbounds Base.$f(A::AbstractDimArray, i::Union{Colon,AbstractVector{<:Integer}}) =
            Base.$f(parent(A), i)
        @propagate_inbounds Base.$f(A::AbstractDimArray, i::AbstractArray{<:Bool}) =
            Base.$f(parent(A), i)
        # Except 1D DimArrays
        @propagate_inbounds Base.$f(A::AbstractDimArray{<:Any, 1}, i::Union{Colon,AbstractVector{<:Integer}}) =
            rebuildsliced(Base.$f, A, Base.$f(parent(A), i), (i,))
        # Dimension indexing. Allows indexing with A[somedim=At(25.0)] for Dim{:somedim}
        @propagate_inbounds Base.$f(A::AbstractDimArray, args::Dimension...; kw...) =
            Base.$f(A, dims2indices(A, (args..., kwdims(values(kw))...))...)
        # Standard indices
        @propagate_inbounds Base.$f(A::AbstractDimArray, i1::StandardIndices, i2::StandardIndices, I::StandardIndices...) =
            rebuildsliced(Base.$f, A, Base.$f(parent(A), i1, i2, I...), (i1, i2, I...))
        @propagate_inbounds Base.$f(A::AbstractDimArray, I::CartesianIndex) = Base.$f(parent(A), I)
    end
end

#### setindex ####

@propagate_inbounds Base.setindex!(A::AbstractDimArray, x) = setindex!(parent(A), x)
@propagate_inbounds Base.setindex!(A::AbstractDimArray, x, args::Dimension...; kw...) =
    setindex!(A, x, dims2indices(A, (args..., kwdims(values(kw))...))...)
@propagate_inbounds Base.setindex!(A::AbstractDimArray, x, i, I...) =
    setindex!(A, x, dims2indices(A, (i, I...))...)
@propagate_inbounds Base.setindex!(A::AbstractDimArray, x, i1::StandardIndices, I::StandardIndices...) =
    setindex!(parent(A), x, i1, I...)
@propagate_inbounds Base.setindex!(A::AbstractDimArray, x, I::CartesianIndex) =
    setindex!(parent(A), x, I)


text:
using graphtheoryufrj
using DataStructures

"""
Requirement (PT-BR): Calcule a distancia e o menor caminho a partir dos
vertices 10, 20, 30, 40, 50 para todos os outros vertices do grafo. Obtenha
o tempo de execucao para calcular todas estas distancias. Reporte a distancia
e um caminho minimo entre os vertices de origem e o vertice 1.
"""
function test_one(infile_name::ASCIIString)
	println("- Starting the Test One")
	try
		initial_vertex = [10 20 30 40 50]

		println("\n-- Performing 5 Dijkstras on the Adjacency List: ")
		graph_list = simpleGraph(infile_name)

		# Do the Dijkstra to find the shortest path
		@time for i in initial_vertex
			@time dist,_ = dijkstra(graph_list.Graph, i)
			println(dist[1])
		end
	catch
		error("Error in \"test_one\" function.")
	end
	return "TEST_ONE_SUCESS"
end

"""
Requirement (PT-BR): Obtenha uma arvore geradora minima, informando seu peso.
Obtenha o tempo de execucao para resolver este problema.
"""
function test_two(infile_name::ASCIIString, algorithm::Function, initial_vertex)
	println("- Starting the Test Two")
	try
		println("\n-- Performing the algorithm to obtain the MST")
		graph_list = simpleGraph(infile_name)
		@time dist,_ = algorithm(graph_list.Graph, initial_vertex)
		println("MST weight: ", sum(dist))
	catch
		error("Error in \"test_two\" function.")
	end
	return "TEST_TWO_SUCESS"
end

"""
Requirement (PT-BR): Calcule a distancia media do grafo. Obtenha o tempo de
execucao para resolver este problema.
"""
function test_three(infile_name::ASCIIString)
	println("- Starting the Test Three")
	#try
		println("\n-- Performing the mean_distance algorithm")
		graph_list = simpleGraph(infile_name)
		@time mean = mean_distance_dijkstra(graph_list.Graph)
		println("Mean distance: ", mean)
	#catch
	#	error("Error in \"test_three\" function.")
	#end
	return "TEST_THREE_SUCESS"
end

"""
Requirement (PT-BR): Calcule a distancia e o caminho minimo entre Edsger W.
Dijkstra (o pesquisador) e os seguintes pesquisadores na rede de colaboracao:
Alan M. Turing, J. B. Kruskal, Jon M. Kleinberg, Eva Tardos, Daniel R. Figueiredo.
"""
function test_colab_one()
"""
2722 - Dijkstra
11365 - Turing
471365 - Kruskal
5709 - Kleinberg
11386 - Eva Tardos
343930 - Ratton
"""
	println("- Starting the Test One of the colaboration network")
	initial_vertex = 2722

	try
		println("\n-- Performing the algorithm to obtain the MST")
		graph_list = simpleGraph("../assets/rede_colaboracao.txt")
		@time dist, prev = dijkstra(graph_list.Graph, initial_vertex)
		println("Dijkstra to Turing: ", dist[11365])
		println("Dijkstra to Kruskal: ", dist[471365])
		println("Dijkstra to Kleinberg: ", dist[5709])
		println("Dijkstra to Eva Tardos: ", dist[11386])
		println("Dijkstra to Ratton: ", dist[343930])

		# Parents of Ratton
		u = 343930
		while(u != 2722)
			println(u)
			u = prev[u]
		end

		# 3 biggest degrees
		# @TODO: Improve this shit - Too many workaround
		for i in 1:length(dist)
			if(dist[i] == 100000000)
				dist[i] = 0
			end
		end

		i = findmax(dist)
		println("3 biggest degrees: ")
		println(i)
		dist[i[2]] = 0
		i = findmax(dist)
		println(i)
		dist[i[2]] = 0
		i = findmax(dist)
		println(i)

	catch
		error("Error in \"test_colab_one\" function.")
	end

	return "TEST_COLAB_ONE_SUCESS"
end

"""
Requirement (PT-BR): Obtenha uma arvore geradora minima, e responda as
seguintes perguntas:
(a) Determine os tres vertices de maior grau na MST.
(b) Determine os vizinhos de Edsger W. Dijkstra e de Daniel R. Figueiredo na MST.
"""
function test_colab_two()
	println("- Starting the Test Two of the colaboration network")
	try

	catch
		error("Error in \"test_colab_two\" function.")
	end
	return "TEST_COLAB_TWO_SUCESS"
end


text:
# rectified linear: f(x)=max(x,0)
rectlin(x)=x.*(x.>0)

Frectlin(x::Array{Float64,2})=(max(x,0),[]);
Frectlin_inplace(value,auxvalue,x::Array{Float64,2})=copy!(value,max(x,0))

Drectlin(derivativeIDX,f_c,faux_c,grad_c,grad_n,x)=axpy!(1.0,grad_c.*(x.>0),grad_n)

if PROC=="GPU"
    Frectlin(x::CudaArray)=(rectlin(x),[])
    Frectlin_inplace(value,aux,x::CudaArray)=rectlin!(x,value)

    Drectlin(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::CudaArray)=A_emult_Bg0!(grad_c,x,grad_n)#; t.*(x.>0)
end

Derivative[Frectlin]=Drectlin
Inplace[Frectlin]=Frectlin_inplace

rectlin(A::ADnode)=ADnode(Frectlin,A)
export rectlin

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff