text:
###Types to hold intermediate results of geometry calculations
abstract AbstractGeometry{O<:DerivativeOrder}

type BaseGeometry <: AbstractGeometry{ZeroOrder}
end

type GradientGeometry <: AbstractGeometry{FirstOrder}
end

type TensorGeometry <: AbstractGeometry{SecondOrder}
end

### Functions that, given a model and associated evaluation function, calculate the geometry of the parameter space

### Calculate the geometry over all samples
function geometry!(m::AbstractModel,s::AbstractSample)
    ###Calculation of logprior is the same for all types of models and samples
    logprior!(s.logprior,parameters(m),s.values)
    ###For all samaples falling inside the prior, calculate the geometry
    for i=1:numsamples(s)
        if isfinite(s.logprior[i])
            geometry!(m,s,i)
        else
            s.loglikelihood[i] = -Inf
        end
    end
    return s
end

function geometry!(m::AbstractModel,s::AbstractSample{ZeroOrder},i::Integer)
    r = evaluate!(m,s.values[:,i])
    s.loglikelihood[i] = loglikelihood(m,r)
end

# function geometry!(m::AbstractModel,s::AbstractSample{FirstOrder},i::Integer)
#     r = evaluate!(m,s.values[:,i])
#     s.loglikelihood[i] = loglikelihood(m,r)
#     gradlogprior!(m,s,r,i)
#     gradloglikelihood!(m,s,r,i)
# end

# function geometry!(m::AbstractModel,s::AbstractSample{SecondOrder},i::Integer)
#     r::AbstractArray = evaluate!(m,s.values[:,i])
#     s.loglikelihood[i] = loglikelihood(m,r)
#     gradlogprior!(m,s,r,i)
#     gradloglikelihood!(m,s,r,i)j
#     tensorlogprior!(m,s,r,i)
#     tensorloglikelihood!(m,s,r,i)
# end

### Calculate the likelihood of a set of data values given the measurement data
@inline loglikelihood(m::AbstractModel,data::AbstractArray) = loglikelihood(noisemodel(m),measurements(m),data)

# ###Calculation of the gradient from finite differences
# function gradient(m::AbstractModel,s::AbstractSample{GradientOrder},g::AbstractGeometry{GradientOrder})
#     #pre-allocate the array for finite difference results
#     #fds  = similar(evaldata,size(evaldata)...,nparas(s))
# end

#gradlogprior!(m::ODEModel,s::Sample{}) = (s.gradlogprior = map((p,v)->((logpdf(p,v+m.gradientepsilon)-logpdf(p,v))/m.gradientepsilon),m.parameters.priors,s.values))

###Calculate the loglikelihood function

###Calculate the gradient of the loglikelihood
#gradienthelper(m::ODEModel,data::Matrix{Float64},sol::Matrix{Float64},grad::Array{Float64,3}) = vec(sum((grad.*(data-sol)./m.variance),(1,2))) #.* does automatic broadcast(), summing result along first 2 dimensions
#gradloglikelihood(m::ODEModel,sol::Matrix{Float64},grad::Array{Float64,3}) = gradienthelper(m,m.measurements,sol,grad)

###Calculate the metric tensor
#tensorvalue(m::ODEModel,grad::Array{Float64,3},i::Int,j::Int) = sum(grad[:,:,i].*grad[:,:,j]./m.variance)

###Generate pseudodata for the approximate metric tensor calculations
#pseudodata(m::ODEModel,d::AbstractArray) = applynoise!(m.noise,d)
#tangentvector(m::ODEModel,sol::Matrix{Float64},grad::Array{Float64,3}) = gradienthelper(m,pseudodata(m,sol),sol,grad)

###Calculate the gradient of the loglikelihood
#gradienthelper(m::ODEModel,data::Matrix{Float64},sol::Matrix{Float64},grad::Array{Float64,3}) = vec(sum((grad.*(data-sol)./m.variance),(1,2))) #.* does automatic broadcast(), summing result along first 2 dimensions
#gradloglikelihood(m::ODEModel,sol::Matrix{Float64},grad::Array{Float64,3}) = gradienthelper(m,m.measurements,sol,grad)

#   #calculate finite differences
#   for i = 1:length(s.values)
#     np = copy(s.values)
#     np[i]  = np[i] + m.gradientepsilon
#     realepsilon = np[i] - s.values[i]
#     fd[:,:,i] = (evaluate(m,np)-sol)/realepsilon
#   end
#   s.gradloglikelihood = gradloglikelihood(m,sol,fd)

#   #return the finite differences solution
#   fd

# function tensor!(m::MCModel,t::TensorSample,sol::Matrix{Float64},grad::Array{Float64,3})
#   np = length(t.values)
#   for i=1:np
#     for j=i:np
#       t.tensorloglikelihood[i,j] = tensorvalue(m,grad,i,j)
#       t.tensorloglikelihood[j,i] = t.tensorloglikelihood[i,j]
#     end
#   end
# end

# function tensor!(m::MCModel,t::ApproximateTensorSample,sol::Matrix{Float64},grad::Array{Float64,3})
#   for i=1:size(t.tangentvectors,2)
#     t.tangentvectors[:,i] = tangentvector(m,sol,grad)
#   end
#   t.tensorloglikelihood = (t.tangentvectors*t.tangentvectors')/length(t.tangentvectors)
# end


text:
@doc raw"""
    UCME(k, testpredictions, testtargets)

Estimator of the unnormalized calibration mean embedding (UCME) with kernel `k` and sets of
`testpredictions` and `testtargets`.

Kernel `k` on the product space of predictions and targets has to be a `Kernel` from the
Julia package
[KernelFunctions.jl](https://github.com/JuliaGaussianProcesses/KernelFunctions.jl)
that can be evaluated for inputs that are tuples of predictions and targets.

The number of test predictions and test targets must be the same and at least one.

# Details

The estimator is biased and guaranteed to be non-negative. Its sample complexity
is ``O(mn)``, where ``m`` is the number of test locations and ``n`` is the total number of
samples.

Let ``(T_i)_{i=1,\ldots,m}`` be the set of test locations, i.e., test predictions and
corresponding targets, and let ``(P_{X_j}, Y_j)_{j=1,\ldots,n}`` be a data set of
predictions and corresponding targets. The plug-in estimator of ``\mathrm{UCME}_{k,m}^2``
is defined as
```math
m^{-1} \sum_{i=1}^{m} {\bigg(n^{-1} \sum_{j=1}^n k\big(T_i, (P_{X_j}, Y_j)\big)
- \mathbb{E}_{Z \sim P_{X_j}} k\big(T_i, (P_{X_j}, Z)\big)\bigg)}^2.
```

# References

Widmann, D., Lindsten, F., & Zachariah, D. (2021).
[Calibration tests beyond classification](https://openreview.net/forum?id=-bxf89v3Nx).
To be presented at *ICLR 2021*.
"""
struct UCME{K<:Kernel,TP,TT} <: CalibrationErrorEstimator
    """Kernel."""
    kernel::K
    """Test predictions."""
    testpredictions::TP
    """Test targets."""
    testtargets::TT

    function UCME{K,TP,TT}(kernel::K, testpredictions::TP, testtargets::TT) where {K,TP,TT}
        check_nsamples(testpredictions, testtargets)
        return new{K,TP,TT}(kernel, testpredictions, testtargets)
    end
end

function UCME(kernel::Kernel, testpredictions, testtargets)
    return UCME{typeof(kernel),typeof(testpredictions),typeof(testtargets)}(
        kernel, testpredictions, testtargets
    )
end

function (estimator::UCME)(predictions::AbstractVector, targets::AbstractVector)
    @unpack kernel, testpredictions, testtargets = estimator

    # obtain number of samples
    nsamples = check_nsamples(predictions, targets)

    # compute average over test locations
    estimate = mean(zip(testpredictions, testtargets)) do (tp, ty)
        unsafe_ucme_eval_testlocation(kernel, predictions, targets, tp, ty)
    end

    return estimate
end

function unsafe_ucme_eval_testlocation(
    kernel::Kernel,
    predictions::AbstractVector,
    targets::AbstractVector,
    testprediction,
    testtarget,
)
    # compute average over predictions and targets for the given test location
    estimate = mean(zip(predictions, targets)) do (p, y)
        unsafe_ucme_eval(kernel, p, y, testprediction, testtarget)
    end

    return estimate^2
end

function unsafe_ucme_eval(
    kernel::Kernel,
    p::AbstractVector{<:Real},
    y::Integer,
    testp::AbstractVector{<:Real},
    testy::Integer,
)
    res = sum(((z == y) - pz) * kernel((p, z), (testp, testy)) for (z, pz) in enumerate(p))
    return res
end
function unsafe_ucme_eval(kernel::Kernel, p::Real, y::Bool, testp::Real, testy::Bool)
    noty = !y
    return (y - p) * kernel((p, y), (testp, testy)) +
           (noty - p) * kernel((p, noty), (testp, testy))
end

function unsafe_ucme_eval(kernel::KernelTensorProduct, p, y, testp, testy)
    κpredictions, κtargets = kernel.kernels
    return unsafe_ucme_eval_targets(κtargets, p, y, testp, testy) * κpredictions(p, testp)
end
# resolve method ambiguity
function unsafe_ucme_eval(
    kernel::KernelTensorProduct,
    p::AbstractVector{<:Real},
    y::Integer,
    testp::AbstractVector{<:Real},
    testy::Integer,
)
    κpredictions, κtargets = kernel.kernels
    return unsafe_ucme_eval_targets(κtargets, p, y, testp, testy) * κpredictions(p, testp)
end
function unsafe_ucme_eval(
    kernel::KernelTensorProduct, p::Real, y::Bool, testp::Real, testy::Bool
)
    κpredictions, κtargets = kernel.kernels
    return unsafe_ucme_eval_targets(κtargets, p, y, testp, testy) * κpredictions(p, testp)
end

function unsafe_ucme_eval_targets(
    kernel::Kernel,
    p::AbstractVector{<:Real},
    y::Integer,
    testp::AbstractVector{<:Real},
    testy::Integer,
)
    res = sum(((z == y) - pz) * kernel(z, testy) for (z, pz) in enumerate(p))
    return res
end
function unsafe_ucme_eval_targets(
    kernel::Kernel, p::Real, y::Bool, testp::Real, testy::Bool
)
    noty = !y
    res = (y - p) * kernel(y, testy) + (noty - p) * kernel(noty, testy)
    return res
end

function unsafe_ucme_eval_targets(
    κtargets::WhiteKernel,
    p::AbstractVector{<:Real},
    y::Integer,
    testp::AbstractVector{<:Real},
    testy::Integer,
)
    @inbounds res = (y == testy) - p[testy]
    return res
end
function unsafe_ucme_eval_targets(
    kernel::WhiteKernel, p::Real, y::Bool, testp::Real, testy::Bool
)
    return (2 * testy - 1) * (y - p)
end

ext: jl


text:
module DiffusionModels

using Parameters, Distributions

import Base: rand, get

export get, fpt, rand, ddm_fpt, ddm_rand
export AbstractDrift, ConstDrift, VarDrift
export AbstractSigma, ConstSigma, VarSigma
export AbstractBounds
export AbstractConstBounds, ConstSymBounds, ConstAsymBounds
export AbstractVarBounds, VarSymBounds, VarAsymBounds
export DecayingBound, DecayingBounds
export AbstractNDT, NonDecisionTime
export DiffusionModel

include("drifttypes.jl")
include("sigmatypes.jl")
include("boundtypes.jl")
include("ndttypes.jl")
include("diffusiontypes.jl")
include("dm.jl")

end # module


text:
# The Numerical stability checks code in this file is based on the code from the SDDP.jl package,
# from the below mentioned commit and file.
# commit :8cd305188caffc50a1734913053fc81bba613778
# link to file :https://github.com/odow/SDDP.jl/blob/d353fe5a2903421e7fed6d609eb9377c35d715a1/src/print.jl#L190

########### JuMP model utils #########

mutable struct NumericalBounds
    min::Float64
    max::Float64
    min_index::Any
    max_index::Any
end

NumericalBounds() = NumericalBounds(Inf, -Inf, nothing, nothing)

set_min!(v::NumericalBounds, value::Real) = v.min = value
set_max!(v::NumericalBounds, value::Real) = v.max = value
set_min_index!(v::NumericalBounds, idx) = v.min_index = idx
set_max_index!(v::NumericalBounds, idx) = v.max_index = idx

mutable struct ConstraintBounds
    coefficient::NumericalBounds
    rhs::NumericalBounds
    function ConstraintBounds()
        return new(NumericalBounds(), NumericalBounds())
    end
end

function _update_coefficient_bounds(
    v::ConstraintBounds,
    constraint::JuMP.ScalarConstraint,
    idx,
)
    _update_numerical_bounds(v.coefficient, constraint.func, idx)
    return
end

function _update_rhs_bounds(v::ConstraintBounds, constraint::JuMP.ScalarConstraint, idx)
    _update_numerical_bounds(v.rhs, constraint.set, idx)
end

mutable struct VariableBounds
    bounds::NumericalBounds
    function VariableBounds()
        return new(NumericalBounds())
    end
end

function _update_variable_bounds(v::VariableBounds, variable::JuMP.VariableRef, idx)
    if JuMP.is_binary(variable)
        set_min!(v.bounds, 0.0)
        _update_numerical_bounds(v.bounds, 1.0, idx)
    else
        if JuMP.has_lower_bound(variable)
            _update_numerical_bounds(v.bounds, JuMP.lower_bound(variable), idx)
        end
        if JuMP.has_upper_bound(variable)
            _update_numerical_bounds(v.bounds, JuMP.upper_bound(variable), idx)
        end
    end
    return
end

function _update_numerical_bounds(v::NumericalBounds, value::Real, idx)
    if !isapprox(value, 0.0)
        if v.min > abs(value)
            set_min!(v, value)
            set_min_index!(v, idx)
        elseif v.max < abs(value)
            set_max!(v, value)
            set_max_index!(v, idx)
        end
    end
    return
end

function _update_numerical_bounds(bonuds::NumericalBounds, func::JuMP.GenericAffExpr, idx)
    for coefficient in values(func.terms)
        _update_numerical_bounds(bonuds, coefficient, idx)
    end
    return
end

function _update_numerical_bounds(bonuds::NumericalBounds, func::MOI.LessThan, idx)
    return _update_numerical_bounds(bonuds, func.upper, idx)
end

function _update_numerical_bounds(bonuds::NumericalBounds, func::MOI.GreaterThan, idx)
    return _update_numerical_bounds(bonuds, func.lower, idx)
end

function _update_numerical_bounds(bonuds::NumericalBounds, func::MOI.EqualTo, idx)
    return _update_numerical_bounds(bonuds, func.value, idx)
end

function _update_numerical_bounds(bonuds::NumericalBounds, func::MOI.Interval, idx)
    _update_numerical_bounds(bonuds, func.upper, idx)
    return _update_numerical_bounds(bonuds, func.lower, idx)
end

# Default fallback for unsupported constraints.
_update_numerical_bounds(range::NumericalBounds, func, idx) = nothing

function get_constraint_numerical_bounds(model::OperationsProblem; verbose = false)
    if verbose
        constraint_bounds = Dict()
        for (const_key, constriant_array) in
            get_constraints(get_optimization_container(model))
            bounds = ConstraintBounds()
            for idx in Iterators.product(constriant_array.axes...)
                con_obj = JuMP.constraint_object(constriant_array[idx...])
                _update_coefficient_bounds(bounds, con_obj, idx)
                _update_rhs_bounds(bounds, con_obj, idx)
            end
            constraint_bounds[const_key] = bounds
        end
        return constraint_bounds
    else
        bounds = ConstraintBounds()
        for (const_key, constriant_array) in
            get_constraints(get_optimization_container(model))
            for idx in Iterators.product(constriant_array.axes...)
                con_obj = JuMP.constraint_object(constriant_array[idx...])
                _update_coefficient_bounds(bounds, con_obj, (const_key, idx))
                _update_rhs_bounds(bounds, con_obj, (const_key, idx))
            end
        end
        return bounds
    end
end

function get_variable_numerical_bounds(model::OperationsProblem; verbose = false)
    if verbose
        variable_bounds = Dict()
        for (variable_key, variable_array) in
            get_variables(get_optimization_container(model))
            bounds = VariableBounds()
            for idx in Iterators.product(variable_array.axes...)
                var = variable_array[idx...]
                _update_variable_bounds(bounds, var, idx)
            end
            variable_bounds[variable_key] = bounds
        end
        return variable_bounds
    else
        bounds = VariableBounds()
        for (variable_key, variable_array) in
            get_variables(get_optimization_container(model))
            for idx in Iterators.product(variable_array.axes...)
                var = variable_array[idx...]
                _update_variable_bounds(bounds, var, (variable_key, idx))
            end
        end
        return bounds
    end
end


text:
using AutoGrad, Base.Test
# Uncomment these if you want lots of messages:
# import Base.Test: default_handler, Success, Failure, Error
# default_handler(r::Success) = info("$(r.expr)")
# default_handler(r::Failure) = warn("$(r.expr) FAILED")
# default_handler(r::Error)   = warn("$(r.expr): $(r.err)")

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

info("Test indexing...")
a1 = rand(3)                    # FAIL: some high-order (b1sum) tests with rand(2)
t1 = (a1...)
d1 = Dict(); for i=1:length(a1); d1[i]=a1[i]; end

s0(x)=x[1]^2+x[2]^2
s1 = grad(s0)
s1sum(x)=(y=s1(x);y[1]+y[2])
s2 = grad(s1sum)
@test check_grads(s0,a1)
@test check_grads(s0,t1)
@test check_grads(s0,d1)
@test check_grads(s1sum,a1)
@test check_grads(s1sum,t1)
@test check_grads(s1sum,d1)

using AutoGrad: sumvalues
f0(x)=(a=0;for i=1:length(x);a+=x[i]^2;end;a)
f1=grad(f0)
f1sum(x)=sumvalues(f1(x))
f2=grad(f1sum)
@test check_grads(f0,a1)
@test check_grads(f0,t1)
@test check_grads(f0,d1)
@test check_grads(f1sum,a1)
@test check_grads(f1sum,t1)
@test check_grads(f1sum,d1)

r0(x)=(s=0; for i=2:length(x); s+=(1-x[i-1])^2 + 100*(x[i]-x[i-1]^2)^2; end; s)
r1 = grad(r0)
r1sum(x)=sumvalues(r1(x))
r2 = grad(r1sum)
@test check_grads(r0,a1)
@test check_grads(r0,t1)
@test check_grads(r0,d1)
@test check_grads(r1sum,a1)
@test check_grads(r1sum,t1)
@test check_grads(r1sum,d1)

info("Test rosenbrock with map...")
b0(x) = sum(map((i, j) -> (1 - j)^2 + 100*(i - j^2)^2, x[2:end], x[1:end-1]))
b1 = grad(b0)
b1sum(x)=sumvalues(b1(x))
@test check_grads(b0,a1)
@test check_grads(b0,t1)
@test check_grads(b1sum,a1) # fail with size 2
@test check_grads(b1sum,t1) # fail with size 2
@time b1sum(rand(10000))

info("Test higher order gradients...")
g1 = grad(sin); @test g1(1)==cos(1)
g2 = grad(g1);  @test g2(1)==-sin(1)
g3 = grad(g2);  @test g3(1)==-cos(1)
g4 = grad(g3);  @test g4(1)==sin(1)
g5 = grad(g4);  @test g5(1)==cos(1)
g6 = grad(g5);  @test g6(1)==-sin(1)
g7 = grad(g6);  @test g7(1)==-cos(1)
g8 = grad(g7);  @test g8(1)==sin(1)
g9 = grad(g8);  @test g9(1)==cos(1)

info("Test neural net...")
n0(w,x,y)=sum(((w[3]*max(0,w[1]*x.+w[2]).+w[4])-y).^2)
n1 = grad(n0)
n1sum(w,x,y)=sum(map(sum,n1(w,x,y)))
n1sumd(w,x,y)=sum(map(sum,values(n1(w,x,y))))
wa = Any[rand(2,3),rand(2),rand(2,2),rand(2)]
wt = (wa...)
wd = Dict(); for i=1:length(wa); wd[i]=wa[i]; end
@test check_grads(n0, wa, rand(3,10), rand(2,10))
@test check_grads(n0, wt, rand(3,10), rand(2,10))
@test check_grads(n0, wd, rand(3,10), rand(2,10))
@test check_grads(n1sum, wa, rand(3,10), rand(2,10))
@test check_grads(n1sum, wt, rand(3,10), rand(2,10))
# This needs more work:
# @test check_grads(n1sumd, wd, rand(3,10), rand(2,10))  # FAIL

info("Test primitives...")
using AutoGrad: runtests
runtests()

### UTILS...

"Find out where different methods are."
function where(k)
    f = eval(k)
    a = (rand(), rand(2), rand(2,2), rand(2,2,2))
    for x in a
        try println(@which f(x)); end
        for y in a
            try println(@which f(x,y)); end
        end
    end
end

"See which scalar/array args f accepts."
function testargs(k)
    print("$k:")
    f = eval(k)
    try f(rand()); print(" (N,)"); end
    try f(rand(2)); print(" (V,)"); end
    try f(rand(2,2)); print(" (M,)"); end
    try f(rand(2,2,2)); print(" (T,)"); end
    try f(rand(),rand()); print(" (N,N)"); end
    try f(rand(),rand(2)); print(" (N,V)"); end
    try f(rand(),rand(2,2)); print(" (N,M)"); end
    try f(rand(),rand(2,2,2)); print(" (N,T)"); end
    try f(rand(2),rand()); print(" (V,N)"); end
    try f(rand(2),rand(2)); print(" (V,V)"); end
    try f(rand(2),rand(2,2)); print(" (V,M)"); end
    try f(rand(2),rand(2,2,2)); print(" (V,T)"); end
    try f(rand(2,2),rand()); print(" (M,N)"); end
    try f(rand(2,2),rand(2)); print(" (M,V)"); end
    try f(rand(2,2),rand(2,2)); print(" (M,M)"); end
    try f(rand(2,2),rand(2,2,2)); print(" (M,T)"); end
    try f(rand(2,2,2),rand()); print(" (T,N)"); end
    try f(rand(2,2,2),rand(2)); print(" (T,V)"); end
    try f(rand(2,2,2),rand(2,2)); print(" (T,M)"); end
    try f(rand(2,2,2),rand(2,2,2)); print(" (T,T)"); end
    println()
end

function test1arg(x; o...)
    for k in keys(math1arg)
        fx = eval(k)
        f = isa(x,Number) ? fx : a->sum(fx(a))
        (f===fx) || name(f,(:sum,k))
        @test check_grads(f, x; o...)
    end
end

function test2arg(x1,x2,flist=keys(math2arg); o...)
    for k in flist
        fx = eval(k)
        f = isa(x1,Number) && isa(x2,Number) ? fx : (a,b)->sum(fx(a,b))
        (f===fx) || name(f,(:sum,k))
        @test check_grads(f, x1, x2; o...)
    end
end

function test2arg_old(x1,x2;eps=1e-5)
    # First, second, or both arguments can be arrays
    # Grad only works for scalar output, so we'll sum the results
    psum(a)=(isa(a,Number) ? a : (s=0;for i=1:length(a);s+=a[i];end;s))
    for k in keys(math2arg)
        fx = eval(k)
        f = (x1,x2)->psum(fx(x1,x2))
        y = f(x1,x2)
        dbg((k,:x,x1,x2,:y,y))
        g1,g2 = grad(f,1),grad(f,2)
        dx = [g1(x1,x2),g2(x1,x2)]
        # Need to separate here based on array or not, maybe write a gcheck?
        Dx = [(f(x1+eps,x2)-f(x1-eps,x2))/2eps, Dx2 = (f(x1,x2+eps)-f(x1,x2-eps))/2eps]
        println((k,:diff,abs(dx-Dx),:dx,dx,:Dx,Dx))
        isapprox2(dx, Dx) || error("$k derivative is wrong: $((dx,Dx))")
    end
    return true
end

function test2arg_scalar(eps=1e-5)
    for k in keys(math2arg)
        f = eval(k)
        g1,g2 = grad(f,1),grad(f,2)
        x1,x2 = rand(2)
        y = f(x1,x2)
        dbg((k,:x,x1,x2,:y,y))
        dx = [g1(x1,x2),g2(x1,x2)]
        Dx = [(f(x1+eps,x2)-f(x1-eps,x2))/2eps, Dx2 = (f(x1,x2+eps)-f(x1,x2-eps))/2eps]
        println((k,:diff,abs(dx-Dx),:dx,dx,:Dx,Dx))
        # @test_approx_eq_eps dx Dx 1e-6
        isapprox2(dx, Dx) || error("$k derivative is wrong: $((dx,Dx))")
    end
    return true
end

function test1arg_scalar(eps=1e-5)
    for k in keys(math1arg)
        f = eval(k)
        g = grad(f)
        x = rand()
        y = f(x)
        dbg((k,:x,x,:y,y))
        dx = g(x)
        Dx = (f(x+eps)-f(x-eps))/2eps
        println((k,:diff,abs(dx-Dx),:dx,dx,:Dx,Dx))
        # @test_approx_eq_eps dx Dx 1e-6
        isapprox(dx, Dx) || error("$k derivative is wrong: $((dx,Dx))")
    end
    return true
end

function test1arg_array(eps=1e-5)
    for k in keys(math1arg)
        fx = eval(k)
        f = x->(a=fx(x);a[1]+a[2])
        g = grad(f)
        x = rand(2,3)
        y = f(x)
        dbg((k,:x,x,:y,y))
        dx = g(x)
        Dx = similar(x)
        for i=1:length(x)
            xi=x[i]
            x[i]=xi+eps; f1=f(x)
            x[i]=xi-eps; f2=f(x)
            x[i]=xi
            Dx[i]=(f1-f2)/2eps
        end
        println((k,:maxdiff,maximum(abs(dx-Dx)),:avgdx,mean(abs(dx)))) # :dx,dx,:Dx,Dx))
        isapprox2(dx, Dx) || error("$k derivative is wrong: $((dx,Dx))")
    end
    return true
end

function isapprox2(x, y;
                   maxeps::Real = max(eps(eltype(x)), eps(eltype(y))),
                   rtol::Real=maxeps^(1/3), atol::Real=maxeps^(1/2))
    size(x) == size(y) || (warn("isapprox: $(size(x))!=$(size(y))"); return false)
    x = convert(Array, x)
    y = convert(Array, y)
    d = abs(x-y)
    s = abs(x)+abs(y)
    maximum(d - rtol * s) <= atol
end

# include("math1arg.jl")
# @test test1arg_scalar()
# @test test1arg_array()

function test1(x=(1.,2.))
    # foo(x)=sin(x[1])+cos(x[2])
    foo(x) = x[1]+x[2]
    goo = grad(foo)
    @show goo(x)
    # This does not work because goo is not scalar valued!
    # hoo = grad(goo)
    # @show hoo(x)
end

function test2()
    gsin = grad(sin)
    hsin = grad(gsin)
    #@show sin(1.0)
    #@show gsin(1.0)
    @show hsin(1.0)
end

function test3()
    foo2(x,y)=sin(x)+cos(y)
    goo2 = grad(foo2)
    goo22 = grad(foo2, 2)
    @show goo2(1,2)
    @show goo22(1,2)
end


text:
using ShirleyRayTracer
using Test

function single_scancol()
    scene = Scene(Camera(Point3(13.,2.,3.), zero(Point3), Vec3(0,1,0), 20, 16/9, 0.1, 10.0))
    add!(scene, Sphere(Point3(0, 1, 0), 1.0, Dielectric(1.5)))
	add!(scene, Sphere(Point3(-4, 1, 0), 1.0, Lambertian(0.4,0.2,0.1)))
	add!(scene, Sphere(Point3(4, 1, 0), 1.0, Metal(0.7,0.6,0.5, 0.0)))
    width = 1200
    height = 100
    nsamples = 1
    max_depth = 5
    x = 600
	image = Vector{RGB{N0f8}}(undef, height)
	for y in 1:height
        image[y] = render_pixel(x, height-y+1, scene, height, width, nsamples, max_depth)
    end

    image
end

@testset "ShirleyRayTracer.jl" begin
    @test single_scancol()
end


text:
using Test, LinearAlgebra, Convex, Suppressor
using CVChannel

# Check whether a Mosek license exists in the test ENV
const has_MOSEK_LICENSE = hasMOSEKLicense()

@testset "./src/optimizer_interface.jl" begin

@testset "qsolve!()" begin
    function _state_optimization()
        Π = [0.5 0.5;0.5 0.5]
        ρ = HermitianSemidefinite(2)

        objective = real(tr(Π*ρ))
        constraints = [
            real(tr(ρ)) == 1,
            imag(tr(ρ)) == 0
        ]

        return ρ, objective, constraints
    end

    # NOTE: These test should really be testing which optimizer is running, but
    # Suppressor is not able to capture the outputs for some unknown reason.
    @testset "runs with SCS backend" begin
        (ρ, objective, constraints) = _state_optimization()
        problem = maximize(objective, constraints)

        # unable to capture output from SCS with suppressor
        # testing in default quiet mode
        qsolve!(problem)

        @test problem.optval ≈ 1
        @test ρ.value ≈ [0.5 0.5;0.5 0.5]
    end

    if has_MOSEK_LICENSE
        @testset "runs with mosek backend" begin
            (ρ, objective, constraints) = _state_optimization()
            problem = maximize(objective, constraints)

            mosek_out = @capture_out qsolve!(problem, quiet=false, use_mosek=true)

            @test occursin(r"Problem\n\s*Name\s*:\s*\n\s*Objective sense\s*:\s*max\s*\n", mosek_out)
            @test problem.optval ≈ 1
            @test isapprox(ρ.value, [0.5 0.5;0.5 0.5], atol=1e-5)
        end
    else
        @warn "MOSEK tests skipped because a license was not found."
    end
end

@testset "hasMOSEKLicense()" begin
    test_dir = joinpath(@__DIR__,"files")

    # WARNING: Modifying ENV is a good way to shoot yourself in the foot.
    # making a copy of the environment variables
    env_copy = copy(ENV)

    # restores the environment variables to original values
    _restoreENV() = begin
        ENV["MOSEKLM_LICENSE_FILE"] = haskey(env_copy, "MOSEKLM_LICENSE_FILE") ? env_copy["MOSEKLM_LICENSE_FILE"] : ""
        ENV["HOME"] = haskey(env_copy, "HOME") ? env_copy["HOME"] : ""
        ENV["PROFILE"] = haskey(env_copy, "PROFILE") ? env_copy["PROFILE"] : ""
    end

    @testset "modifying ENV to be different license states" begin
        ENV["MOSEKLM_LICENSE_FILE"] = ""
        ENV["HOME"] = ""
        ENV["PROFILE"] = ""
        @test !hasMOSEKLicense()

        ENV["MOSEKLM_LICENSE_FILE"] = joinpath(test_dir,"mosek","mosek.lic")
        @test hasMOSEKLicense()

        ENV["MOSEKLM_LICENSE_FILE"] = ""
        ENV["HOME"] = test_dir
        @test hasMOSEKLicense()

        ENV["HOME"] = ""
        ENV["PROFILE"] = test_dir
        @test hasMOSEKLicense()

        _restoreENV()
    end

    # NOTE: This must run last!
    # In case any tests failed, making sure environment variables are properly restored
    _restoreENV()
end

end

max_forks_repo_name: ChitambarLab/cv-channel


text:
#=P.1.Tricks & Tips in Numerical Computing
JuliaCon 2018 - M\cr|NA
Credits: Professor Mr. Nick Higham @nhigham
# today the lecturing Professor @nhigham of Numerical Computing
University of Manchesterh

http://maths.manchester.ac.uk/~higham

nickhigham.wordpress.com

=#

"""P.2.Introduction

1. Exploiting complex Artihmetic
2. Understanding multivalued functions
3. Exploiting associativity of matrix products
4. Randomization to avoid pathological cases
5. Making the most of Low precision arithmetic

http://bit.ly/tricks18
#redirects to a seemingly empty site

A trick used threetimes becomes a standard technique  -George Polya -Hungarian mathematician [ETH Zürich] (Palo Alto), CA, U.S.
#---

"""
x = 10^5

#f(x) =
#f = f(x)

#to view only (secant method )
# f'(x) = (f(x + h) - f(x)) / h

#=p.3.Differentiation With(out) a Difference

≈ abs(f(x))

Error O(h) (especialy)

Note:
this is a differential function
=#
range = 1:10
size(range)[1]

#=
calculates b on different subset of ranges (for each n)
```inputs:

a: initial position
n: steps required (see: nLookup)
h: stepsize
```

Example |--- |...|--- |
        a1  b1   bn-1  bn
say a = 1, b3 =4 , & h 1 (pre-defined constant) , let ni be n1, n2,...ni =s.t.= 1,2,3,...  for i = 1,2,3,.. , i∈ N+ Positive Natural numbers,
then: # for i =1,2,3,..
b[i] = a*(n[i]*h)
e.g.

b[1] = a[1] + (n[1] * h) = 1 + (1 * 1 )  = 2  #infer: sub-range [1,2]  (b=4 !isa b[1]=2  (now a[2] = b[1]=2 for next op ) a starts at 2
b[2] = a[2] + (n[2] * h) = 2 + (2 * 1)  =  4  #infer: sub-rage [2,4]   (b4=4 isa b[2]=4 HALT!)
b is reached , with 2 steps
Recursing back from  b[i]= b = 4 : if a[1] !=1 yet (2 steps )
if b
b[4]

=#
using Test,  BenchmarkTools, DataStructures

"""
Returns the number of steps

    Example:

    for range [a= 1 , b=3]
    size of range = b - a = 3 - 1 = 2  range length
    if h is 1 then stepsrequired = size(range) * h = 2 * 1 = 2 steps
"""
function nLookup(a, b, h)

    _size = nothing  # b -a
    _size = b - a # range length
    stepsrequired = ceil( _size * h )  #steps #whatever h is , it returns ceiled whole number
    return stepsrequired
end
    """
    f_prime

    ```inputs
    f: the function operator
    x: the single variable
    a: lower range bound
    b: upper range bound
    h: the stepsize : initialized to 1
    ```

    example:

   for range [a= 1 , b=3]
    size of range = b - a = 3 - 1 = 2  range length
    if h is 1 then stepsrequired = size(range) * h = 2 * 1 = 2 steps
    a , a+h, a+2h, ...,b
    a=1 , 1+1, 1+1+1 = b
    """
function f_prime!(f,x,a,b,;h=1)

    result =nothing
    n = nLookup(a,b,h) #returns `number of steps` required to do  n  :: # is it always a whole number #expectation : whole number
    sum = 0
    if i > n  && i < 1 #Upon fallback
        #return result

    else #otherwise # do the work
        #which variables are changing
        #nth  step:

      #  f_prime!()

    end

    if b >= a
        result = f(x=b) - f(x=a) / (b - a)

    elseif b <= a
        result = f(x = a) - f(x=b) / (a - b)

    return result

end

function f_double_prime(f, x, a = 1, b = 10)

    for i in enumerate(2)

    f = f_prime!(f, x, a, b) #f(b) - f(a)
    δfδx = deepcopy(f)
end

function h(u, x, f, f_double_prime) # f(x) is not a valid function Argument #Solution: remove arguments
    #pythagorean theorem
    return h ≈ sqrt(u * abs(f(x) / f_double_prime(x)))
end

# Q,which order of complexity is this algorithm?
#=
function (model=sqrt(v[m, n]), approx=x_minus_one(v[m, n])
    for i in enumerate(m) #rows
        #tmp = sum(v, 1) ; #rowErrors[]
        # rowErrors = [rowErrors,tmp]
        #print(tmp)
        for j in enumerate(n) # cols #entering j's context

            error[m, k] = abs(sqrt(v[m, n]) - x_minus_one(v[m, n]))
            k += 1

            #REVEALED: Sqare Array is a must
            #DimensionMismatch (matrix is not square dimensions are 10x1 )# either sqrt or x_minus_one requires rectangular dimensions
        end
    end
end
=#

#Example: the square-root sqrt: most irrational rational

a = sqrt(x);
b = x_minus_one(x) = x^-1; #  Matrix{Float64} <: DenseMatrix{Float64} <: AbstractMatrix{Float64} <: Any
m = 10;
n = 10;
v = rand(m, n) # sqrt only accepted a Squared matrix
sqrt(v) ≈ x_minus_one(v) #false they are not equal!!!!!!!!! # not even a false # Q.what is ϵ here (I have never framed ϵ) just working with the Maestro for not
#once matrix is rectangular i.e. 10x10
sqrt(v) == x_minus_one(v)
#for all item vij
sumError = nothing;
error = zeros(m, n) #nothing ;

rowErrors = nothing;
colErrors = nothing;
tmp = nothing;

k = 1

for i in enumerate(m) #rows
    #tmp = sum(v, 1) ; #rowErrors[]
    # rowErrors = [rowErrors,tmp]
    #print(tmp)
    for j in enumerate(n) # cols #entering j's context
        #tmp = sum(v, 2)
        #   colErrors = [colErrors, tmp] # col error Array
        #  print(tmp)
        error[m, k] = abs(sqrt(v[m, n]) - x_minus_one(v[m, n]))
        k += 1
        #sumError += error #sum error of each time
        #REVEALED: hidden Issue
        #DimensionMismatch (matrix is not square dimensions are 10x1 )# either sqrt or x_minus_one requires rectangular dimensions
    end
end
#mean = ∑
print(error)
#print(colErrors)
#print(rowErrors)
using BenchmarkTools
function dispersion()
colDispersion = _maxCol - _minCol

end

rowErrors = cumsum(v, dims = 1)
#=10×10 Matrix{Float64}:
 0.755532  0.590352  0.111428  0.887575  0.0507037  0.903207  0.358227  0.186171  0.761759  0.883808
 1.10254   0.596795  0.428979  0.958866  0.156567   1.09377   0.999251  0.762358  1.12141   1.36038
 ⋮                                                  ⋮
 4.3395    3.67298   3.29483   6.04047   3.35655    4.74542   2.36259   4.56519   3.20877   5.0396
 4.87931   4.34269   4.2687    **6.44228**   3.38897    5.0163    *3.0492*    4.72759   4.09149   5.76818
 largest error (column independent) col 4:  6.44228 (>5.76818 last cumuli)
 min error            col(7): *3.0492*
 error range = [*3.0492*, 6.44228]
=#
6.44228 - 3.0492
size(v)
colErrors = cumsum(v, dims = 2)

_minCol = minimum(colErrors)
_maxCol = maximum(colErrors)
colDispersion = _maxCol - _minCol

minimum(A[:, 4])
#=
 0.755532   1.34588   1.45731   2.34489   2.39559   3.2988   3.65703  3.8432   4.60496  **5.48876**
 0.347007   0.35345   0.671001  0.742292  0.848155  1.03871  1.67974  2.25592  2.61558  3.09215
 ⋮                                                  ⋮
 0.0863792  0.132913  0.530735  1.41574   1.61419   1.62391  2.04341  2.63593  2.69061  3.63327
 0.539811   1.20953   2.18339   2.5852    2.61763   2.88851  3.57511  3.73751  4.62024  5.34882

 Max:largest (row independent) error row(1) 10,1: 5.48876
 Min: (row independent) error
 min ( 0.0863792) error  0.0863792
 =#

@btime rowErrormat = cumsum(error, dims = 1) 143.349 ns (1allocation:896bytes)
_maxRow = maximum(rowErrormat)
_minRow = minimum(rowErrormat)
rowDispersion = _maxRow - _minRow

colDispersion < rowDispersion # false positive
colDispersion > rowDispersion  #true positive
colDispersion > rowDispersion === true #it should be true,the only statement under question# it displays false!
ratio = max(colDispersion, rowDispersion) / min(colDispersion, rowDispersion)
ratio = rowDispersion / colDispersion
#---below is an overkill
xor(colDispersion < rowDispersion, colDispersion > rowDispersion) #always true [true,false]
and = &
and(colDispersion < rowDispersion, colDispersion > rowDispersion) # hiccup with and # always false

#= there is an error in estimating the sqrt function, perturbuted into
 0.0       0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0       0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 ⋮                             ⋮
 0.0       0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.518967  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

 only 1 error created  a big problem (of not a Neglegible size)
=#
colErrormat = cumsum(error, dims = 2)
end
end
#=
 0.0       0.0       0.0       0.0       0.0       0.0       0.0       0.0       0.0       0.0
 0.0       0.0       0.0       0.0       0.0       0.0       0.0       0.0       0.0       0.0
 ⋮                                                 ⋮
 0.0       0.0       0.0       0.0       0.0       0.0       0.0       0.0       0.0       0.0
 0.518967  0.518967  0.518967  0.518967  0.518967  0.518967  0.518967  0.518967  0.518967  0.518967

 #
=#

#=
for f(x) = exp(x) with x = 1
Optimal

h≈ (u*abs(f(x)/f''(x)))^1/2
Prof. Higham:
and you could see the error

"error touches the red line then error stars to blow up"
me: approx. at (10^8) - 10^7
this is a well known phenomenon - v-shaped curve

the errors in the evaluation starting to
Dominate the truncation errors
mathematically going to zero

what is the point where moment is reached
approx. the value shown here
-unit round-off: 16 *mod f(x) / mod (f''(x))

you're limited to how much accuracy you're going to get

Q.can we do better?
=#

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
@testset "Test apparent wind angle" begin
    @test isapprox(0.666980, sail_route.awa(60, 3.086, 5.144), rtol=3)
end


text:
# This is just an experiment with linear only tasks #

### General Utility Functions ###

# Data Read: Reads the structured data from the given instance file
# Input:  - filename: A string containing the file name with its relative path.
# Output: - n: An Integer containing the number of tasks;
#         - tasks_times: An (n) Array of Integers with each task's time;
#         - prec_relations: A (n,n) 2-dimensional Binary Array representing
#                           each task's Precedence Relation.
#         - times_sum: An Integer representing the sum of all the tasks' times.
function data_read(filename::String)
    # Gets every line in the instance file for processing
    lines = []
    open(filename,"r") do f
        lines = readlines(filename)
    end

    # The number n of tasks is given in the first line
    n = parse(Int64, lines[1])

    # Gets each task's time in the next n lines and the sum of all times
    task_times = Int64[]
    times_sum = 0
    for i = 2:(n+1)
        curr_time = parse(Int64, lines[i])
        append!(task_times, curr_time)
        times_sum += curr_time
    end

    # Gets the i,j Direct Precedence Relations until the "-1,-1" end mark
    # and transforms it into a Binary Matrix of the Relations
    prec_relations = zeros(Int64, (n,n))
    curr_line = n+2
    # Julia has no "Do While" so a first verification is necessary to avoid a break
    rs = split(lines[curr_line], ",")
    r = [parse(Int64, rs[1]), parse(Int64, rs[2])]
    while (r[1] > 0)
        prec_relations[r[1], r[2]] = 1

        curr_line += 1
        rs = split(lines[curr_line], ",")
        r = [parse(Int64, rs[1]), parse(Int64, rs[2])]
    end

    println("Data successfully read from file \"", filename, "\" containing ", n, " tasks.")
    println("")

    # Returns the structured data
    return n, task_times, prec_relations, times_sum
end

# Solution Write: Writes the final solution and time in a file with stdout format
# Input:  - filename: A string with the name for the file that will be written;
#         - solution: The final best solution encountered during execution;
#         - time: The time it took to find the given solution.
# Output: None.
function solution_write(filename::String, solution, time)
    open(filename, "w") do f
        if (time == 0)
            println(f, "No feasible solution found.\n")

            println("- FINAL RESULT: No feasible solution found.")
        else
            println(f, "Best solution: ", solution[2], "\r\n")
            println(f, "With score: ", solution[1], "\r\n")
            println(f, "Found in: ", time, " seconds.\r\n")

            println("- FINAL RESULT: ")
            println("- Best solution: ", solution[2])
            println("- With score: ", solution[1])
            println("- Found in: ", time, " seconds.")
        end
    end

    println("Solution saved in file: ", filename)
    println("")
end

# Args Read: Reads any (or provides default) Argments passed through command line
# Input:  - args: An Array of Strings containing the argments from the command line.
# Output: - file: A string with the name of the file for the instance given;
#         - stations: An Integer with the number of stations available, defaults to 1;
#         - population: An Integer with the size of the population, defaults to 100.
function args_read(args = ["instances/HAHN.IN2", "1", "500"])
    if (length(args) == 0)
        println("WARNING: No argument provided, using default instance Hahn and stations and population sizes of 1 and 100, respectively")
        file = "instances/HAHN.IN2"
        stations = 1
        population = 500
    elseif (length(args) == 1)
        println("WARNING: Only one argument provided, using default stations and population sizes of 1 and 100, respectively")
        file = args[1]
        stations = 1
        population = 500
    elseif (length(args) == 2)
        println("WARNING: Population size not provided, using default of 500")
        file = args[1]
        stations = parse(Int64, args[2])
        population = 500
    else
        file = args[1]
        stations = parse(Int64, args[2])
        population = parse(Int64, args[3])
    end

    println("Arguments successfully read. Using these parameters:")
    println("File Path and Name: \"", file, "\".")
    println("Number of Stations: ", stations, ".")
    println("Size of Population: ", population, ".")
    println("")

    return file, stations, population
end

### Genetic Algorithm Functions ###

# Rand Gene: Generates a random gene containing a solution for the problem
# Input:  - s: An Integer representing the given number of stations available;
#         - n: An Integer representing the number of tasks to be assigned.
# Output: - gene: A (s,n) 2-dimensional Integer Array representing the sequence
#                 of tasks n assigned to each station s.
function rand_gene(s::Int64, n::Int64)
    gene = Int[]
    task_counter = ones(Int64, s)

    if (s > 1)
        cut = rand(1:n)
        for i = 1:(s-1)
            append!(gene, cut)
            cut = rand(cut:n)
        end
    end

    return gene
end

# Rand Population:
# Input:  - size: An Integer representing the size of the generated population;
#         - s, n, task_times: Needed from gene and fitness functions.
# Output: - population: A random population containing "size"s [fitness, gene] pairs.
function rand_population(size::Int64, s::Int64, n::Int64, task_times)
    population = []
    for i = 1:size
        gene = rand_gene(s, n)
        fit = fitness(s, n, gene, task_times)
        curr_individual = [fit, gene]
        append!(population, [curr_individual])
    end

    return population
end

# Crossover: Creates a new Gene by doing a crossover between two parents
# Input:  - parent_1: A [fitness, gene] individual representing the first parent;
#         - parent_2: A [fitness, gene] individual representing the second parent.
#         - s, n, task_times: Needed from fitness function.
# Output: - new_individual: A new [fitness, gene] individual.
function crossover(s, n, parent1, parent2, task_times)
    gene1 = parent1[2]
    gene2 = parent2[2]

    new_gene = gene1
    rnd = rand(1:2)

    if (rnd == 1)
        for i = 1:(s-1)
            new_gene[i] = min(new_gene[i], gene2[i])
        end
    else
        for i = 1:(s-1)
            it = s - i
            new_gene[it] = max(new_gene[it], gene2[it])
        end
    end

    fit = fitness(s, n, new_gene, task_times)
    new_individual = [fit, new_gene]

    return new_individual
end

# Fitness: Calculates the fitness of the given gene
#          In SALBP, it's the longest cycle among all stations
#          A cycle of a station is the sum the times of all the tasks given to it
# Input:  - gene: A (s,n) 2-dimensional Integer Array representing a gene;
#         - task_times: An (n) Array of Integers with each task's time.
#         - s, n: Needed from gene function.
# Output: - max_cycle: An Integer representing the fitness of the given gene.
function fitness(s, n, gene, task_times)
    max_cycle = 0

    if (s == 1)
        max_cycle = sum(task_times)
    elseif (s == 2)
        max_cycle = max(task_times[1:gene[1]], task_times[(gene[1]+1):n])
    else
        curr_cut = 0
        for i = 1:(s-1)
            if(curr_cut != gene[i])
                curr_cut += 1
                curr_cycle = sum(task_times[curr_cut:gene[i]])
                curr_cut = gene[i]
                max_cycle = max(max_cycle, curr_cycle)
            end
        end
        curr_cut += 1
        if (curr_cut < n)
            curr_cycle = sum(task_times[curr_cut:n])
            max_cycle = max(max_cycle, curr_cycle)
        end
    end

    return max_cycle
end

# Genetic Algorithm:
# Input:  - init_pop: An array containing the initial [fitness, gene] pairs population
#                     of the previously defined size;
#         - pop_size, s, n, task_times, times_sum, prec_relations: Needed from gene, population
#                                                                  and fitness functions;
#         - max_gen: An integer containing the stop criteria of the GA, max generations.
# Output: - best_solution: The [fitness, gene] pair of the best solution found;
#         - time_taken: The time it took to find the best solution found;
function genetic_algorithm(init_pop, pop_size::Int64, s::Int64, n::Int64, task_times, times_sum::Int64, max_gen::Int64, prec_relations)
    tic()
    population = init_pop
    best_found_fit = times_sum + 1
    best_solution = [0, 0]
    time_taken = 0

    println("Genetic Algorithm started.")

    # Until the end criteria is hit
    for i = 1:max_gen
        population = sort(population, lt=lexless)
        curr_best = population[1][1]

        if (curr_best <= times_sum)
            if (curr_best < best_found_fit)
                best_found_fit = curr_best
                best_solution = population[1]
                time_taken += toc()
                tic()
                println("IT [", i, "] - ", "Found a better solution with score: ", best_solution[1])
            end
        end

        # Now to create a new population for the next generation
        class_A = ceil(Int64, pop_size * 0.2)
        class_B = ceil(Int64, pop_size * 0.3)
        class_C = ceil(Int64, pop_size * 0.5)

        # Class A stays the same, changes start in Class B
        for i = class_A:class_C-1
            # Class B crossover
            curr_parents = i - class_A + 1
            population[i] = crossover(s, n, population[curr_parents], population[curr_parents + 1], task_times)
        end
        # Class C just replaces the worst half of the population for new random genes
        population[(class_C+1):end,:] = rand_population(class_C, s, n, task_times)
    end

    return best_solution, time_taken
end

### Main Function ###
function main()
    max_gen = 500
    out_filename = "main_test.txt"

    println("- Reading arguments...")
    println("")
    file_name, num_stations, population_size = args_read(ARGS)

    println("- Setting up data from given file...")
    println("")
    num_tasks, task_times, prec_relations, times_sum = data_read(file_name)

    println("- Starting Genetic Algorithm execution...")
    println("")
    initial_population = rand_population(population_size, num_stations, num_tasks, task_times)
    best_solution, time_taken = genetic_algorithm(initial_population, population_size, num_stations, num_tasks, task_times, times_sum, max_gen, prec_relations)

    solution_write(out_filename, best_solution, time_taken)

    println("- Execution ended.")
    println("")
end

main()

max_stars_repo_name: AlexandreSSJr/SALBP-GA

max_issues_repo_name: AlexandreSSJr/SALBP-GA

max_forks_repo_name: AlexandreSSJr/SALBP-GA


text:
import DelayEmbeddings.Dataset

import ..InvariantMeasures:
    TransferOperatorEstimatorRectangularBinning,
    TransferOperatorRectangularBinning
import CausalityToolsBase: encode, get_minima_and_edgelengths, RectangularBinning

function invariantmeasure(data, ϵ::RectangularBinning,
        estimator = :TransferOperatorEstimatorRectangularBinning;
        kwargs...) where T

    if estimator == :TransferOperatorEstimatorRectangularBinning
        # Encode the points (identify which bin, given the
        # rectangular binning scheme ϵ, each point falls in)
        mini, edgelengths = get_minima_and_edgelengths(data, ϵ)
        encoded_pts = encode(data, mini, edgelengths)

        # Which are the visited bins, which points
        # visits which bin, repetitions, etc...
        binvisits = get_binvisits(encoded_pts)

        # Use that information to estimate transfer operator
        TO = estimate_transferoperator_from_binvisits(binvisits)

        return invariantmeasure(TO; kwargs...)
    end
end

function invariantmeasure(data::StateSpaceReconstruction.AbstractEmbedding,
        ϵ::RectangularBinning,
        estimator = :TransferOperatorEstimatorRectangularBinning;
        kwargs...)

    invariantmeasure(data.points, ϵ, estimator, kwargs...)
end

function invariantmeasure(data::Dataset,
        ϵ::RectangularBinning,
        estimator = :TransferOperatorEstimatorRectangularBinning;
        kwargs...)

    invariantmeasure(transpose(Matrix(data)), ϵ, estimator, kwargs...)
end

"""
    invariantmeasure(data,
        ϵ::Union{Int, Float64, Vector{Int}, Vector{Float64}},
        estimator = TransferOperatorEstimatorRectangularBinning;
        kwargs...)

Estimate the invariant measure from a rectangular partition of `data` using
the estimator provided with the `estimator` argument. The default
is `estimator = :TransferOperatorEstimatorRectangularBinning`.

This is done by discretizing the state space into rectangular bins with
edge lengths dictated by the binning scheme `ϵ`. We then approximate the
transfer operator over the discretized state space, and compute the
invariant measure over the bins from the transfer operator.

Returns an [`InvariantDistribution`](@ref) instance.

## Arguments
- **`data`**: The data from which to estimate the invariant measure. The following
    data types are currently accepted:
    - `Dataset` instances from `DynamicalSystems.jl`/`DelayEmbeddings.jl`.
    - `AbstractEmbedding` subtypes from `StateSpaceReconstruction.jl`.
    - `AbstractArray{T, 2}` instances where each column represents a point.
- **`ϵ`**: A valid `RectangularBinning` binning scheme.
"""
function invariantmeasure(data,
        ϵ::RectangularBinning,
        estimator::Symbol; kwargs...)
    invariantmeasure(data, ϵ, estimator, kwargs...)
end

export invariantmeasure

"""
    invariantmeasure(to::AbstractTransferOperator;
            N::Int = 200,
            tolerance::Float64 = 1e-8,
            delta::Float64 = 1e-8)

Compute the invariant probability distribution from a transfer operator.

## Computing an invariant probability distribution
The distribution is taken as a left eigenvector of the transfer matrix,
obtained by repeated application of the transfer operator on a randomly
initialised distribution until the probability distribution converges.

## Keyword arguments
- `N`: the maximum number of iterations.
- `tolerance` and `delta`: decides when convergence is achieved.
"""
function invariantmeasure(to::AbstractTransferOperator;
			N::Int = 200,
            tolerance::Float64 = 1e-8,
            delta::Float64 = 1e-8)
    #=
    # Start with a random distribution `Ρ` (big rho). Normalise it so that it
    # sums to 1 and forms a true probability distribution over the simplices.
    =#
    Ρ = rand(Float64, 1, size(to.transfermatrix, 1))
    Ρ = Ρ ./ sum(Ρ, dims = 2)

    #=
    # Start estimating the invariant distribution. We could either do this by
    # finding the left-eigenvector of M, or by repeated application of M on Ρ
    # until the distribution converges. Here, we use the latter approach,
    # meaning that we iterate until Ρ doesn't change substantially between
    # iterations.
    =#
    distribution = Ρ * to.transfermatrix

    distance = norm(distribution - Ρ) / norm(Ρ)

    check = floor(Int, 1 / delta)
    check_pts = floor.(Int, transpose(collect(1:N)) ./ check) .* transpose(collect(1:N))
    check_pts = check_pts[check_pts .> 0]
    num_checkpts = size(check_pts, 1)
    check_pts_counter = 1

    counter = 1
    while counter <= N && distance >= tolerance
        counter += 1
        Ρ = distribution

        # Apply the Markov matrix to the current state of the distribution
        distribution = Ρ * to.transfermatrix

        if (check_pts_counter <= num_checkpts &&
           counter == check_pts[check_pts_counter])

            check_pts_counter += 1
            colsum_distribution = sum(distribution, dims = 2)[1]
            if abs(colsum_distribution - 1) > delta
                distribution = distribution ./ colsum_distribution
            end
        end

        distance = norm(distribution - Ρ) / norm(Ρ)
    end
    distribution = dropdims(distribution, dims = 1)

    # Do the last normalisation and check
    colsum_distribution = sum(distribution)

    if abs(colsum_distribution - 1) > delta
        distribution = distribution ./ colsum_distribution
    end

    # Find partition elements with strictly positive measure.
    δ = tolerance/size(to.transfermatrix, 1)
    simplex_inds_nonzero = findall(distribution .> δ)

    # Extract the elements of the invariant measure corresponding to these indices
    return PerronFrobenius.InvariantDistribution(distribution,simplex_inds_nonzero)
end

export invariantmeasure

"""
    invariantmeasure(to::AbstractTransferOperator;
            N::Int = 200,
            tolerance::Float64 = 1e-8,
            delta::Float64 = 1e-8)

Compute the invariant probability distribution from a transfer operator.

## Computing an invariant probability distribution
The distribution is taken as a left eigenvector of the transfer matrix,
obtained by repeated application of the transfer operator on a randomly
initialised distribution until the probability distribution converges.

## Keyword arguments
- `N`: the maximum number of iterations.
- `tolerance` and `delta`: decides when convergence is achieved.
"""
function invariantmeasure(to::AbstractTransferOperator;
			N::Int = 200,
            tolerance::Float64 = 1e-8,
            delta::Float64 = 1e-8)
    #=
    # Start with a random distribution `Ρ` (big rho). Normalise it so that it
    # sums to 1 and forms a true probability distribution over the simplices.
    =#
    Ρ = rand(Float64, 1, size(to.transfermatrix, 1))
    Ρ = Ρ ./ sum(Ρ, dims = 2)

    #=
    # Start estimating the invariant distribution. We could either do this by
    # finding the left-eigenvector of M, or by repeated application of M on Ρ
    # until the distribution converges. Here, we use the latter approach,
    # meaning that we iterate until Ρ doesn't change substantially between
    # iterations.
    =#
    distribution = Ρ * to.transfermatrix

    distance = norm(distribution - Ρ) / norm(Ρ)

    check = floor(Int, 1 / delta)
    check_pts = floor.(Int, transpose(collect(1:N)) ./ check) .* transpose(collect(1:N))
    check_pts = check_pts[check_pts .> 0]
    num_checkpts = size(check_pts, 1)
    check_pts_counter = 1

    counter = 1
    while counter <= N && distance >= tolerance
        counter += 1
        Ρ = distribution

        # Apply the Markov matrix to the current state of the distribution
        distribution = Ρ * to.transfermatrix

        if (check_pts_counter <= num_checkpts &&
           counter == check_pts[check_pts_counter])

            check_pts_counter += 1
            colsum_distribution = sum(distribution, dims = 2)[1]
            if abs(colsum_distribution - 1) > delta
                distribution = distribution ./ colsum_distribution
            end
        end

        distance = norm(distribution - Ρ) / norm(Ρ)
    end
    distribution = dropdims(distribution, dims = 1)

    # Do the last normalisation and check
    colsum_distribution = sum(distribution)

    if abs(colsum_distribution - 1) > delta
        distribution = distribution ./ colsum_distribution
    end

    # Find partition elements with strictly positive measure.
    δ = tolerance/size(to.transfermatrix, 1)
    simplex_inds_nonzero = findall(distribution .> δ)

    # Extract the elements of the invariant measure corresponding to these indices
    return PerronFrobenius.InvariantDistribution(distribution, simplex_inds_nonzero)
end


text:
using netgengen

#defaults
Rref = 1.0
min_radius=1
max_radius=2
padding = max_radius*2
geofilename = "randomcereal.geo"
siffilename = "randomcereal.sif"

a = 8
b = 17
p = [0.0,0.0,0.0]
q = [1.0,0.0,0.0]

num_coil = 10
num_tries = 10000

type Coil
  a::Float64
  b::Float64
  radius::Float64
  center::Array{Float64,1}
  normal::Array{Float64,1}
  csg_geom::torus
  Coil(name::String,p,q,R::Float64;modifiers...) =
  (x = new();
   x.a = a;
   x.b = b;
   x.center=p;
   x.normal=q;
   x.radius = R;
   x.csg_geom = torus(name, x.center, x.normal, x.b*x.radius/(x.a+x.b), x.a*x.radius/(x.a+x.b);modifiers...);
   return x)
end

function coil_intersects(coil_a::Coil, coil_b::Coil, sep=0.0)
  coil_distance(coil_a, coil_b) = norm(coil_a.center-coil_b.center)
  if coil_distance(coil_a, coil_b) > coil_a.radius + coil_b.radius + sep
    return false
  end
  return true
end

function cram_coil!(coils::Array{Coil,1}, coil)
  for x in coils
    if coil_intersects(x, coil)
      return false
    end
  end
  push!(coils, coil)
  return true
end

coil_array = Array{Coil,1}()
n = 1
m = 1
sizes=linspace(min_radius,max_radius,100)

# make the coils
while n <= num_coil && m <= num_tries
  normal = randn(3)
  normal = normal/norm(normal)
  origin = randn(3)*Rref
  radius = rand(sizes)
  if cram_coil!(coil_array, Coil("coil_$(n)", origin, normal, radius))
    n = n + 1
  else
    m = m + 1
  end
end

# find enclosing box + padding
P = coil_array[1].center-coil_array[1].radius
Q = coil_array[1].center+coil_array[1].radius

for c in coil_array
  for k = 1:3
    if c.center[k]-c.radius<=P[k]+padding
      P[k] = c.center[k]-c.radius-padding
    end
  end
  for k = 1:3
    if c.center[k]+c.radius>=Q[k]-padding
      Q[k] = c.center[k]+c.radius+padding
    end
  end
end

not_tori=Array{CSGObject,1}()
for x in coil_array
  push!(not_tori, not(x.csg_geom))
end
airbrick=brick("airbrick",P,Q,bc=1,transparent=true)
push!(not_tori, airbrick)
air = intersection("air",not_tori)

geo_io=open("randomcereal.geo","w")
print(geo_io, "algebraic3d\n")
tlo(air,geo_io)
for x in coil_array
  tlo(x.csg_geom,geo_io)
end
close(geo_io)

sifinfo = open("randomcereal.sif","w")
for n = 1:length(coil_array)
  c = coil_array[n]
  co = c.center
  cn = c.normal
  print(sifinfo,"!Coil $(n)------------------------\n")
  print(sifinfo,"coil normal(3) = real $(cn[1]) $(cn[2]) $(cn[3])\n")
  print(sifinfo,"coil origin(3) = real $(co[1]) $(co[2]) $(co[3])\n\n")
end
close(sifinfo)

max_stars_repo_name: juhanikataja/randomcereal

max_issues_repo_name: juhanikataja/randomcereal

max_forks_repo_name: juhanikataja/randomcereal


text:
# Autogenerated wrapper script for Dex_jll for x86_64-linux-musl
export dex, webtemplates

JLLWrappers.@generate_wrapper_header("Dex")
JLLWrappers.@declare_executable_product(dex)
JLLWrappers.@declare_file_product(webtemplates)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        dex,
        "bin/dex",
    )

    JLLWrappers.@init_file_product(
        webtemplates,
        "share/webtemplates.tar.gz",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using StanSample, Test

binom_model = "
// Inferring a Rate
data {
  int<lower=1> n;
  int<lower=0> k;
}
parameters {
  real<lower=0,upper=1> theta;
  real<lower=0,upper=1> thetaprior;
}
model {
  // Prior Distribution for Rate Theta
  theta ~ beta(1, 1);
  thetaprior ~ beta(1, 1);

  // Observed Counts
  k ~ binomial(n, theta);
}
generated quantities {
  int<lower=0> postpredk;
  int<lower=0> priorpredk;

  postpredk <- binomial_rng(n, theta);
  priorpredk <- binomial_rng(n, thetaprior);
}
"

sm = SampleModel("binomial", binom_model)

binom_data = Dict("n" => 10, "k" => 5)

(sample_file, log_file) = stan_sample(sm, data=binom_data)

if !(sample_file == nothing)
  chn = read_samples(sm)
  describe(chn)

  # Ceate a ChainDataFrame
  summary_df = read_summary(sm)
  @test summary_df[:theta, :mean][1] ≈ 0.24 atol=0.8

end


text:
using ApproxFun, SingularIntegralEquations, Base.Test

println("Chebyshev addition test")

d = Interval(-1.,1.)
fK(x,y) = exp(abs2(y-x))
FLR = LowRankFun(fK,Chebyshev(d),Chebyshev(d))
FPF = ProductFun(fK,Chebyshev(d),Chebyshev(d))
G = convolutionProductFun(fK,Chebyshev(d),Chebyshev(d))
G1 = convolutionProductFun(fK,Ultraspherical(1,d),Chebyshev(d))
G2 = convolutionProductFun(fK,Ultraspherical(1,d),ChebyshevDirichlet{1,1}(d))

@test norm(pad(coefficients(FLR),size(coefficients(G),1),size(coefficients(G),2))-coefficients(G))<100000eps()
@test norm(pad(coefficients(FPF),size(coefficients(G),1),size(coefficients(G),2))-coefficients(G))<100000eps()
@test norm(fK(.123,.456)-G(.123,.456))≤10000eps()
@test norm(G(.123,.456)-G1(.123,.456))≤10000eps()
@test norm(G(.123,.456)-G2(.123,.456))≤10000eps()

println("Fourier on PeriodicInterval tests")

f2 = Fun(θ->exp(sin(θ))+sin(cos(θ)),Fourier())
FLR = LowRankFun((θ,ϕ)->f2(ϕ-θ),Fourier(),Fourier())
FPF = ProductFun((θ,ϕ)->f2(ϕ-θ),Fourier(),Fourier())
G = convolutionProductFun(f2,Fourier(),Fourier())

@test norm(pad(coefficients(FLR),size(coefficients(G),1),size(coefficients(G),2))-coefficients(G))<400eps()
@test norm(pad(coefficients(FPF),size(coefficients(G),1),size(coefficients(G),2))-coefficients(G))<400eps()
@test norm(f2(.456-.123)-G(.123,.456))≤400eps()

f2 = Fun(θ->1+sin(cos(θ)),CosSpace())
G = convolutionProductFun(f2,Fourier(),Fourier())
@test norm(f2(.456-.123)-G(.123,.456))≤400eps()

f2 = Fun(θ->sin(sin(θ)),SinSpace())
G = convolutionProductFun(f2,Fourier(),Fourier())
@test norm(f2(.456-.123)-G(.123,.456))≤400eps()

println("Laurent on PeriodicInterval tests")

f2 = Fun(θ->exp(sin(θ))+sin(cos(θ)),Laurent())
G = convolutionProductFun(f2,Laurent(),Laurent())
@test norm(f2(.456-.123)-G(.123,.456))≤400eps()

f2 = Fun(θ->π+e*exp(im*θ)+sqrt(2)*exp(im*2θ)+catalan*exp(im*3θ)+γ*exp(im*4θ),Taylor(PeriodicInterval()))
G = convolutionProductFun(f2,Laurent(),Laurent())
@test norm(f2(.456-.123)-G(.123,.456))≤10eps()

f2 = Fun(θ->e*exp(-im*θ)+sqrt(2)*exp(-im*2θ)+catalan*exp(-im*3θ)+γ*exp(-im*4θ),Hardy{false}(PeriodicInterval()))
G = convolutionProductFun(f2,Laurent(),Laurent())
@test norm(f2(.456-.123)-G(.123,.456))≤10eps()

println("Timing tests: ")

d = Interval(-2.5,-.5)
fK2(x,y) = besselj0(100(y-x))
convolutionProductFun(fK2,Chebyshev(d),Chebyshev(d))
@time G = convolutionProductFun(fK2,Chebyshev(d),Chebyshev(d))
convolutionProductFun(fK2,Ultraspherical(1,d),ChebyshevDirichlet{1,1}(d))
@time G = convolutionProductFun(fK2,Ultraspherical(1,d),ChebyshevDirichlet{1,1}(d))
println("Chebyshev addition: Time should be ~0.01 seconds.")

f2 = Fun(θ->besselj0(500*abs(2sin(θ/2))),CosSpace())
convolutionProductFun(f2,Fourier(),Fourier())
@time G = convolutionProductFun(f2,Fourier(),Fourier())
println("CosSpace addition: Time should be ~0.05 seconds.")

f2 = Fun(θ->besselj0(500*abs(2sin(θ/2))),Laurent())
convolutionProductFun(f2,Laurent(),Laurent())
@time G = convolutionProductFun(f2,Laurent(),Laurent())
println("Laurent addition: Time should be ~0.14 seconds.")


text:
test("Add path to clipper") do
    path = Vector{IntPoint}()

    push!(path, IntPoint(0, 0))
    push!(path, IntPoint(0, 1))

    c = Clip()

    @test add_path!(c, path, PolyTypeSubject, true) == false

    push!(path, IntPoint(1, 1))

    @test add_path!(c, path, PolyTypeSubject, true) == true
end

test("Union") do
  path1 = Vector{IntPoint}()
  push!(path1, IntPoint(0, 0))
  push!(path1, IntPoint(0, 1))
  push!(path1, IntPoint(1, 1))
  push!(path1, IntPoint(1, 0))

  path2 = Vector{IntPoint}()
  push!(path2, IntPoint(1, 0))
  push!(path2, IntPoint(1, 1))
  push!(path2, IntPoint(2, 1))
  push!(path2, IntPoint(2, 0))

  c = Clip()
  add_path!(c, path1, PolyTypeSubject, true)
  add_path!(c, path2, PolyTypeSubject, true)

  result, polys = execute(c, ClipTypeUnion, PolyFillTypeEvenOdd, PolyFillTypeEvenOdd)

  @test result == true
  @test string(polys) == "Array{Clipper.IntPoint,1}[Clipper.IntPoint[[0,0],[2,0],[2,1],[0,1]]]"
end

test("Difference") do
  path1 = Vector{IntPoint}()
  push!(path1, IntPoint(0, 0))
  push!(path1, IntPoint(0, 10))
  push!(path1, IntPoint(10, 10))
  push!(path1, IntPoint(10, 0))

  path2 = Vector{IntPoint}()
  push!(path2, IntPoint(4, 0))
  push!(path2, IntPoint(4, 10))
  push!(path2, IntPoint(6, 10))
  push!(path2, IntPoint(6, 0))

  c = Clip()
  add_path!(c, path1, PolyTypeSubject, true)
  add_path!(c, path2, PolyTypeClip, true)

  result, polys = execute(c, ClipTypeDifference, PolyFillTypeEvenOdd, PolyFillTypeEvenOdd)

  @test result == true
  @test string(polys) == "Array{Clipper.IntPoint,1}[Clipper.IntPoint[[10,10],[6,10],[6,0],[10,0]],Clipper.IntPoint[[0,10],[0,0],[4,0],[4,10]]]"
end

test("GetBounds") do
  path1 = Vector{IntPoint}()
  push!(path1, IntPoint(0, 0))
  push!(path1, IntPoint(0, 1))
  push!(path1, IntPoint(1, 1))
  push!(path1, IntPoint(1, 0))

  path2 = Vector{IntPoint}()
  push!(path2, IntPoint(1, 0))
  push!(path2, IntPoint(1, 1))
  push!(path2, IntPoint(2, 1))
  push!(path2, IntPoint(2, 0))

  c = Clip()
  add_path!(c, path1, PolyTypeSubject, true)
  add_path!(c, path2, PolyTypeSubject, true)

  rect = get_bounds(c)

  @test rect.left == 0
  @test rect.top == 0
  @test rect.right == 2
  @test rect.bottom == 1
end

test("Clear") do
  path1 = Vector{IntPoint}()
  push!(path1, IntPoint(0, 0))
  push!(path1, IntPoint(0, 10))
  push!(path1, IntPoint(10, 10))
  push!(path1, IntPoint(10, 0))

  path2 = Vector{IntPoint}()
  push!(path2, IntPoint(4, 0))
  push!(path2, IntPoint(4, 10))
  push!(path2, IntPoint(6, 10))
  push!(path2, IntPoint(6, 0))

  c = Clip()
  add_path!(c, path1, PolyTypeSubject, true)
  add_path!(c, path2, PolyTypeSubject, true)

  Clipper.clear!(c)

  rect = get_bounds(c)

  @test rect.left == 0
  @test rect.top == 0
  @test rect.right == 0
  @test rect.bottom == 0
end

test("AddPaths") do
  path1 = Vector{IntPoint}()
  push!(path1, IntPoint(0, 0))
  push!(path1, IntPoint(0, 1))
  push!(path1, IntPoint(1, 1))
  push!(path1, IntPoint(1, 0))

  path2 = Vector{IntPoint}()
  push!(path2, IntPoint(1, 0))
  push!(path2, IntPoint(1, 1))
  push!(path2, IntPoint(2, 1))
  push!(path2, IntPoint(2, 0))

  paths = Vector{IntPoint}[path1, path2]

  c = Clip()
  add_paths!(c, paths, PolyTypeSubject, true)

  result, polys = execute(c, ClipTypeUnion, PolyFillTypeEvenOdd, PolyFillTypeEvenOdd)

  @test result == true
  @test string(polys) == "Array{Clipper.IntPoint,1}[Clipper.IntPoint[[0,0],[2,0],[2,1],[0,1]]]"
end

test("Orientation") do
  path1 = Vector{IntPoint}()
  push!(path1, IntPoint(0, 0))
  push!(path1, IntPoint(0, 1))
  push!(path1, IntPoint(1, 1))
  push!(path1, IntPoint(1, 0))

  @test orientation(path1) == false

  path1 = Vector{IntPoint}()
  push!(path1, IntPoint(1, 0))
  push!(path1, IntPoint(1, 1))
  push!(path1, IntPoint(0, 1))
  push!(path1, IntPoint(0, 0))

  @test orientation(path1) == true
end

test("Area") do
  path1 = Vector{IntPoint}()
  push!(path1, IntPoint(0, 0))
  push!(path1, IntPoint(0, 1))
  push!(path1, IntPoint(1, 1))

  @test area(path1) == -0.5
end


text:
"""
    init_dataset(stat_list, time_vec, save_folder; path_runoff)

Create an initial dataset.
"""
function init_dataset(stat_list, time_vec, save_folder;
    path_runoff = "/home/jmg/flood_forecasting/runoff_observed")

    # Read metadata from excel sheet

    df_meta = read_metadata()

    df_meta = subset_metadata(df_meta, stat_list)

    # Dictonary that links drainage basin keys to senorge indices

    dbk_ind = read_dbk_ind()

    dbk_ind = subset_dbk_ind(dbk_ind, df_meta)

    # Dictonary that links drainage basin keys to dataframes with
    # geographic information (landuse, elevation...)

    df_geo = watershed_info(dbk_ind)

    # Process metadata for the individual elevation bands

    metadata_elevbands(save_folder, df_geo, df_meta)

    # Read meteorological data

    met_var = "tm"

    tair_data, tair_time = read_senorge_data(time_vec, met_var, df_geo; scale_func = x -> (x .- 2732.0)./10.0)

    met_var = "rr"

    prec_data, prec_time = read_senorge_data(time_vec, met_var, df_geo; scale_func = x -> x/10.0)

    # Write data to files

    write_met_data(save_folder, "tair", time_vec, tair_data, df_meta)

    write_met_data(save_folder, "prec", time_vec, prec_data, df_meta)

    # Process runoff data

    process_runoff_data(stat_list, time_vec, df_meta, path_runoff, save_folder)

end

max_forks_repo_name: jmgnve/NveData


text:
function solution()
  i, map, n, q = 0, Dict(), 10000, 1

  while true
    q += 1
    p = pop!(map, q, 0)

    if p == 0
      map[q * q] = q

      if i == n
        return q
      else
        i += 1
      end
    else
      x = p + q

      while get(map, x, 0) != 0
        x += p
      end

      map[x] = p
    end
  end
end

if length(ARGS) == 1 && ARGS[1] == "-a"
  println(solution())
  exit(0)
end

for line in eachline(STDIN)
  iters = int(line)

  start = time_ns()
  for i in 1:iters
    solution()
  end
  end_ = time_ns()

  println(end_ - start)
end


text:
mutable struct MSQNode{T}
    @atomic next::Union{MSQNode{T},Nothing}
    @atomic value::T

    MSQNode{T}() where {T} = new{T}(nothing)
    MSQNode{T}(next::Union{MSQNode{T},Nothing}, value::T) where {T} = new{T}(next, value)
end

# TODO: rename it to ConcurrentLinkedQueue?
mutable struct ConcurrentQueue{T}
    @atomic head::MSQNode{T}
    @atomic tail::MSQNode{T}
end

ConcurrentQueue() = ConcurrentQueue{Any}()
function ConcurrentQueue{T}() where {T}
    n = MSQNode{T}()
    return ConcurrentQueue{T}(n, n)
end

Base.eltype(::Type{ConcurrentQueue{T}}) where {T} = T

function Base.push!(queue::ConcurrentQueue{T}, v) where {T}
    v = convert(T, v)
    node = MSQNode{T}(nothing, v)

    tail = @atomic queue.tail
    while true
        next = @atomic tail.next
        if next === nothing
            next, ok = @atomicreplace(tail.next, next => node)
            if ok
                @atomicreplace(queue.tail, tail => node)
                return queue
            end
            tail = @atomic queue.tail
        else
            tail, ok = @atomicreplace(queue.tail, tail => next)
            if ok
                tail = next
            end
        end
    end
end

function ConcurrentCollections.maybepopfirst!(queue::ConcurrentQueue)
    head = @atomic queue.head
    tail = @atomic queue.tail
    while true
        node = @atomic head.next
        head′ = @atomic queue.head
        if head === head′
            if head === tail
                if node === nothing
                    return nothing
                end
                tail, _ = @atomicreplace(queue.tail, tail => node)
                continue  # not need to reload tail
            else
                node::MSQNode
                head, ok = @atomicreplace(queue.head, head => node)
                if ok
                    return Some(node.value)
                end
            end
        else
            head = head′
        end
        tail = @atomic queue.tail
    end
end

function Base.popfirst!(queue::ConcurrentQueue)
    r = maybepopfirst!(queue)
    if r === nothing
        error("queue is empty")
    else
        return something(r)
    end
end