text:
using LinearAlgebra
using Distributions

abstract type SurrogateOptimizationAlgorithm end

struct SRBF <: SurrogateOptimizationAlgorithm end
struct LCBS <: SurrogateOptimizationAlgorithm end
struct EI <: SurrogateOptimizationAlgorithm end
struct DYCORS <: SurrogateOptimizationAlgorithm end
struct SOP{P} <: SurrogateOptimizationAlgorithm
    p::P
end

function merit_function(point,w,surr::AbstractSurrogate,s_max,s_min,d_max,d_min,box_size)
    if length(point)==1
        D_x = box_size+1
        for i = 1:length(surr.x)
            distance = norm(surr.x[i]-point)
            if distance < D_x
                D_x = distance
            end
        end
        return w*(surr(point) - s_min)/(s_max-s_min) + (1-w)*((d_max - D_x)/(d_max - d_min))
    else
        D_x = norm(box_size)+1
        for i = 1:length(surr.x)
            distance = norm(surr.x[i] .- point)
            if distance < D_x
                D_x = distance
            end
        end
        return w*(surr(point) - s_min)/(s_max-s_min) + (1-w)*((d_max - D_x)/(d_max - d_min))
    end
end

"""
The main idea is to pick the new evaluations from a set of candidate points where each candidate point is generated as an N(0, sigma^2)
distributed perturbation from the current best solution.
The value of sigma is modified based on progress and follows the same logic as
in many trust region methods: we increase sigma if we make a lot of progress
(the surrogate is accurate) and decrease sigma when we aren’t able to make progress
(the surrogate model is inaccurate).
More details about how sigma is updated is given in the original papers.

After generating the candidate points, we predict their objective function value
and compute the minimum distance to the previously evaluated point.
Let the candidate points be denoted by C and let the function value predictions
be s(x\\_i) and the distance values be d(x\\_i), both rescaled through a
linear transformation to the interval [0,1]. This is done to put the values on
the same scale.
The next point selected for evaluation is the candidate point x that minimizes
the weighted-distance merit function:

``merit(x) = ws(x) + (1-w)(1-d(x))``

where `` 0 \\leq w \\leq 1 ``.
That is, we want a small function value prediction and a large minimum distance
from the previously evaluated points.
The weight w is commonly cycled between
a few values to achieve both exploitation and exploration.
When w is close to zero, we do pure exploration, while w close to 1 corresponds to exploitation.
"""
function surrogate_optimize(obj::Function,::SRBF,lb,ub,surr::AbstractSurrogate,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=100)
    scale = 0.2
    success = 0
    failure = 0
    w_range = [0.3,0.5,0.7,0.95]

    #Vector containing size in each direction
    box_size = lb-ub
    success = 0
    failures = 0
    dtol = 1e-3*norm(ub-lb)
    d = length(surr.x)
    num_of_iterations = 0
    for w in Iterators.cycle(w_range)
        num_of_iterations += 1
        if num_of_iterations == maxiters
            index = argmin(surr.y)
            return (surr.x[index],surr.y[index])
        end
        for k = 1:maxiters
            incumbent_value = minimum(surr.y)
            incumbent_x = surr.x[argmin(surr.y)]

            new_lb = incumbent_x .- 3*scale*norm(incumbent_x .-lb)
            new_ub = incumbent_x .+ 3*scale*norm(incumbent_x .-ub)

            @inbounds for i = 1:length(new_lb)
                if new_lb[i] < lb[i]
                    new_lb = collect(new_lb)
                    new_lb[i] = lb[i]
                end
                if new_ub[i] > ub[i]
                    new_ub = collect(new_ub)
                    new_ub[i] = ub[i]
                end
            end

            new_sample = sample(num_new_samples,new_lb,new_ub,sample_type)
            s = zeros(eltype(surr.x[1]),num_new_samples)
            for j = 1:num_new_samples
                s[j] = surr(new_sample[j])
            end
            s_max = maximum(s)
            s_min = minimum(s)

            d_min = norm(box_size .+ 1)
            d_max = 0.0
            for r = 1:length(surr.x)
                for c = 1:num_new_samples
                    distance_rc = norm(surr.x[r] .- new_sample[c])
                    if distance_rc > d_max
                        d_max = distance_rc
                    end
                    if distance_rc < d_min
                        d_min = distance_rc
                    end
                end
            end

            #3)Evaluate merit function in the sampled points

            evaluation_of_merit_function = zeros(float(eltype(surr.x[1])),num_new_samples)
            @inbounds for r = 1:num_new_samples
                evaluation_of_merit_function[r] = merit_function(new_sample[r],w,surr,s_max,s_min,d_max,d_min,box_size)
            end
            new_addition = false
            adaptive_point_x = Tuple{}
            diff_x = zeros(eltype(surr.x[1]),d)
            while new_addition == false
                #find minimum
                new_min_y = minimum(evaluation_of_merit_function)
                min_index = argmin(evaluation_of_merit_function)
                new_min_x = new_sample[min_index]
                for l = 1:d
                    diff_x[l] = norm(surr.x[l] .- new_min_x)
                end
                bit_x = diff_x .> dtol
                #new_min_x has to have some distance from krig.x
                if false in bit_x
                    #The new_point is not actually that new, discard it!

                    deleteat!(evaluation_of_merit_function,min_index[1])
                    deleteat!(new_sample,min_index)

                    if length(new_sample) == 0
                        println("Out of sampling points")
                        index = argmin(surr.y)
                        return (surr.x[index],surr.y[index])
                    end
                else
                    new_addition = true
                    adaptive_point_x = Tuple(new_min_x)
                end
            end

            #4) Evaluate objective function at adaptive point
            adaptive_point_y = obj(adaptive_point_x)

            #5) Update surrogate with (adaptive_point,objective(adaptive_point)
            add_point!(surr,adaptive_point_x,adaptive_point_y)

            #6) How to go on?
            if surr(adaptive_point_x) < incumbent_value
                #success
                incumbent_x = adaptive_point_x
                incumbent_value = adaptive_point_y
                if failure == 0
                    success +=1
                else
                    failure = 0
                    success += 1
                end
            else
                #failure
                if success == 0
                    failure += 1
                else
                    success = 0
                    failure += 1
                end
            end

            if success == 3
                scale = scale*2
                if scale > 0.8*norm(ub-lb)
                    println("Exiting, scale too big")
                    index = argmin(surr.y)
                    return (surr.x[index],surr.y[index])
                end
                success = 0
                failure = 0
            end

            if failure == 5
                scale = scale/2
                #check bounds and go on only if > 1e-5*interval
                if scale < 1e-5
                    println("Exiting, too narrow")
                    index = argmin(surr.y)
                    return (surr.x[index],surr.y[index])
                end
                success = 0
                failure = 0
            end
        end
    end
end

"""
SRBF 1D:
surrogate_optimize(obj::Function,::SRBF,lb::Number,ub::Number,surr::AbstractSurrogate,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=100)
"""
function surrogate_optimize(obj::Function,::SRBF,lb::Number,ub::Number,surr::AbstractSurrogate,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=100)
#Suggested by:
#https://www.mathworks.com/help/gads/surrogate-optimization-algorithm.html
    scale = 0.2
    success = 0
    failure = 0
    w_range = [0.3,0.5,0.7,0.95]
    box_size = lb-ub
    success = 0
    failures = 0
    dtol = 1e-3*norm(ub-lb)
    num_of_iterations = 0
    for w in Iterators.cycle(w_range)
        num_of_iterations += 1
        if num_of_iterations == maxiters
            index = argmin(surr.y)
            return (surr.x[index],surr.y[index])
        end
        for k = 1:maxiters
            #1) Sample near incumbent (the 2 fraction is arbitrary here)
            incumbent_value = minimum(surr.y)
            incumbent_x = surr.x[argmin(surr.y)]

            new_lb = incumbent_x-scale*norm(incumbent_x-lb)
            new_ub = incumbent_x+scale*norm(incumbent_x-ub)
            if new_lb < lb
                new_lb = lb
            end
            if new_ub > ub
                new_ub = ub
            end
            new_sample = sample(num_new_samples,new_lb,new_ub,sample_type)

            #2) Create  merit function
            s = zeros(eltype(surr.x[1]),num_new_samples)
            for j = 1:num_new_samples
                s[j] = surr(new_sample[j])
            end
            s_max = maximum(s)
            s_min = minimum(s)

            d_min = box_size + 1
            d_max = 0.0
            for r = 1:length(surr.x)
                for c = 1:num_new_samples
                    distance_rc = norm(surr.x[r]-new_sample[c])
                    if distance_rc > d_max
                        d_max = distance_rc
                    end
                    if distance_rc < d_min
                        d_min = distance_rc
                    end
                end
            end
            #3) Evaluate merit function at the sampled points
            evaluation_of_merit_function = merit_function.(new_sample,w,surr,s_max,s_min,d_max,d_min,box_size)

            new_addition = false
            adaptive_point_x = zero(eltype(new_sample[1]));
            while new_addition == false
                #find minimum
                new_min_y = minimum(evaluation_of_merit_function)
                min_index = argmin(evaluation_of_merit_function)
                new_min_x = new_sample[min_index]

                diff_x = abs.(surr.x .- new_min_x)
                bit_x = diff_x .> dtol
                #new_min_x has to have some distance from krig.x
                if false in bit_x
                    #The new_point is not actually that new, discard it!
                    deleteat!(evaluation_of_merit_function,min_index)
                    deleteat!(new_sample,min_index)
                    if length(new_sample) == 0
                        println("Out of sampling points")
                        index = argmin(surr.y)
                        return (surr.x[index],surr.y[index])
                    end
                else
                new_addition = true
                adaptive_point_x = new_min_x
                end
            end
            #4) Evaluate objective function at adaptive point
            adaptive_point_y = obj(adaptive_point_x)

            #5) Update surrogate with (adaptive_point,objective(adaptive_point)
            add_point!(surr,adaptive_point_x,adaptive_point_y)

            #6) How to go on?
            if surr(adaptive_point_x) < incumbent_value
                #success
                incumbent_x = adaptive_point_x
                incumbent_value = adaptive_point_y
                if failure == 0
                    success +=1
                else
                    failure = 0
                    success += 1
                end
            else
                #failure
                if success == 0
                    failure += 1
                else
                    success = 0
                    failure += 1
                end
            end

            if success == 3
                scale = scale*2
                #check bounds cannot go more than [a,b]
                if scale > 0.8*norm(ub-lb)
                    println("Exiting, scale too big")
                    index = argmin(surr.y)
                    return (surr.x[index],surr.y[index])
                end
                success = 0
                failure = 0
            end

            if failure == 5
                scale = scale/2
                #check bounds and go on only if > 1e-5*interval
                if scale < 1e-5
                    println("Exiting, too narrow")
                    index = argmin(surr.y)
                    return (surr.x[index],surr.y[index])
                end
                sucess = 0
                failure = 0
            end
        end
    end
end

"""
This is an implementation of Lower Confidence Bound (LCB),
a popular acquisition function in Bayesian optimization.
Under a Gaussian process (GP) prior, the goal is to minimize:
``LCB(x) := E[x] - k * \\sqrt{(V[x])}``
default value ``k = 2``.
"""
function surrogate_optimize(obj::Function,::LCBS,lb::Number,ub::Number,krig,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=100)
    #Default value
    k = 2.0
    dtol = 1e-3 * norm(ub-lb)
    for i = 1:maxiters
        new_sample = sample(num_new_samples,lb,ub,sample_type)
        evaluations = zeros(eltype(krig.x[1]), num_new_samples)
        for j = 1:num_new_samples
            evaluations[j] = krig(new_sample[j]) +
                             k*sqrt(std_error_at_point(krig,new_sample[j]))
        end

        new_addition = false
        min_add_x = zero(eltype(new_sample[1]));
        min_add_y = zero(eltype(krig.y[1]));
        while new_addition == false
            #find minimum
            new_min_y = minimum(evaluations)
            min_index = argmin(evaluations)
            new_min_x = new_sample[min_index]

            diff_x = abs.(krig.x .- new_min_x)
            bit_x = diff_x .> dtol
            #new_min_x has to have some distance from krig.x
            if false in bit_x
                #The new_point is not actually that new, discard it!
                deleteat!(evaluations,min_index)
                deleteat!(new_sample,min_index)

                if length(new_sample) == 0
                    println("Out of sampling points")
                    index = argmin(krig.y)
                    return (krig.x[index],krig.y[index])
                end
             else
                new_addition = true
                min_add_x = new_min_x
                min_add_y = new_min_y
            end
        end
        if min_add_y < 1e-6*(maximum(krig.y) - minimum(krig.y))
            return
        else
            if (abs(min_add_y) == Inf || min_add_y == NaN)
                println("New point being added is +Inf or NaN, skipping.\n")
            else
                add_point!(krig,min_add_x,min_add_y)
            end
        end
    end
end

"""
This is an implementation of Lower Confidence Bound (LCB),
a popular acquisition function in Bayesian optimization.
Under a Gaussian process (GP) prior, the goal is to minimize:

``LCB(x) := E[x] - k * \\sqrt{(V[x])}``

default value ``k = 2``.
"""
function surrogate_optimize(obj::Function,::LCBS,lb,ub,krig,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=100)
    #Default value
    k = 2.0
    dtol = 1e-3 * norm(ub-lb)
    for i = 1:maxiters
        d = length(krig.x)
        new_sample = sample(num_new_samples,lb,ub,sample_type)
        evaluations = zeros(eltype(krig.x[1]),num_new_samples)
        for j = 1:num_new_samples
            evaluations[j] = krig(new_sample[j]) +
                             k*sqrt(std_error_at_point(krig,new_sample[j]))
        end

        new_addition = false
        min_add_x = Tuple{}
        min_add_y = zero(eltype(krig.y[1]));
        diff_x = zeros(eltype(krig.x[1]),d)
        while new_addition == false
            #find minimum
            new_min_y = minimum(evaluations)
            min_index = argmin(evaluations)
            new_min_x = new_sample[min_index]
            for l = 1:d
                diff_x[l] = norm(krig.x[l] .- new_min_x)
            end
            bit_x = diff_x .> dtol
            #new_min_x has to have some distance from krig.x
            if false in bit_x
                #The new_point is not actually that new, discard it!
                deleteat!(evaluations,min_index)
                deleteat!(new_sample,min_index)

                if length(new_sample) == 0
                    println("Out of sampling points")
                    index = argmin(krig.y)
                    return (krig.x[index],krig.y[index])
                end
             else
                new_addition = true
                min_add_x = new_min_x
                min_add_y = new_min_y
            end
        end
        if min_add_y < 1e-6*(maximum(krig.y) - minimum(krig.y))
            index = argmin(krig.y)
            return (krig.x[index],krig.y[index])
        else
            min_add_y = obj(min_add_x) # I actually add the objc function at that point
            if (abs(min_add_y) == Inf || min_add_y == NaN)
                println("New point being added is +Inf or NaN, skipping.\n")
            else
                add_point!(krig,Tuple(min_add_x),min_add_y)
            end
        end
    end
end

"""
Expected improvement method 1D
"""
function surrogate_optimize(obj::Function,::EI,lb::Number,ub::Number,krig,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=100)
        dtol = 1e-3*norm(ub-lb)
        eps = 0.01
        for i = 1:maxiters
            new_sample = sample(num_new_samples,lb,ub,sample_type)
            f_max = maximum(krig.y)
            evaluations = zeros(eltype(krig.x[1]),num_new_samples)
            point_found = false
            new_x_max = zero(eltype(krig.x[1]))
            new_y_max = zero(eltype(krig.x[1]))
            while point_found == false
                for j = 1:length(new_sample)
                    std = std_error_at_point(krig,new_sample[j])
                    u = krig(new_sample[j])
                    if abs(std) > 1e-6
                        z = (u - f_max - eps)/std
                    else
                        z = 0
                    end
                    evaluations[j] = (u-f_max-eps)*cdf(Normal(),z) + std*pdf(Normal(),z)
                end
                index_max = argmax(evaluations)
                x_new = new_sample[index_max]
                y_new = maximum(evaluations)
                diff_x = abs.(krig.x .- x_new)
                bit_x = diff_x .> dtol
                #new_min_x has to have some distance from krig.x
                if false in bit_x
                    #The new_point is not actually that new, discard it!
                    deleteat!(evaluations,index_max)
                    deleteat!(new_sample,index_max)
                    if length(new_sample) == 0
                        println("Out of sampling points")
                        index = argmin(krig.y)
                        return (krig.x[index],krig.y[index])
                    end
                 else
                    point_found = true
                    new_x_max = x_new
                    new_y_max = y_new
                end
            end
            if new_y_max < 1e-6*norm(maximum(krig.y)-minimum(krig.y))
                index = argmin(krig.y)
                return (krig.x[index],krig.y[index])
            end
            add_point!(krig,new_x_max,obj(new_x_max))
        end
        println("Completed maximum number of iterations")
end

"""
This is an implementation of Expected Improvement (EI),
arguably the most popular acquisition function in Bayesian optimization.
Under a Gaussian process (GP) prior, the goal is to
maximize expected improvement:

``EI(x) := E[max(f_{best}-f(x),0)``

"""
function surrogate_optimize(obj::Function,::EI,lb,ub,krig,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=100)
        dtol = 1e-3*norm(ub-lb)
        eps = 0.01
        for i = 1:maxiters
            d = length(krig.x)
            new_sample = sample(num_new_samples,lb,ub,sample_type)
            f_max = maximum(krig.y)
            evaluations = zeros(eltype(krig.x[1]),num_new_samples)
            point_found = false
            new_x_max = zero(eltype(krig.x[1]))
            new_y_max = zero(eltype(krig.x[1]))
            diff_x = zeros(eltype(krig.x[1]),d)
            while point_found == false
                for j = 1:length(new_sample)
                    std = std_error_at_point(krig,new_sample[j])
                    u = krig(new_sample[j])
                    if abs(std) > 1e-6
                        z = (u - f_max - eps)/std
                    else
                        z = 0
                    end
                    evaluations[j] = (u-f_max-eps)*cdf(Normal(),z) + std*pdf(Normal(),z)
                end
                index_max = argmax(evaluations)
                x_new = new_sample[index_max]
                y_new = maximum(evaluations)
                for l = 1:d
                    diff_x[l] = norm(krig.x[l] .- x_new)
                end
                bit_x = diff_x .> dtol
                #new_min_x has to have some distance from krig.x
                if false in bit_x
                    #The new_point is not actually that new, discard it!
                    deleteat!(evaluations,index_max)
                    deleteat!(new_sample,index_max)
                    if length(new_sample) == 0
                        println("Out of sampling points")
                        index = argmin(krig.y)
                        return (krig.x[index],krig.y[index])
                    end
                 else
                    point_found = true
                    new_x_max = x_new
                    new_y_max = y_new
                end
            end
            if new_y_max < 1e-6*norm(maximum(krig.y)-minimum(krig.y))
                index = argmin(krig.y)
                return (krig.x[index],krig.y[index])
            end
            add_point!(krig,Tuple(new_x_max),obj(new_x_max))
        end
        println("Completed maximum number of iterations")
end

function adjust_step_size(sigma_n,sigma_min,C_success,t_success,C_fail,t_fail)
    if C_success >= t_success
        sigma_n = 2*sigma_n
        C_success = 0
    end
    if C_fail >= t_fail
        sigma_n = max(sigma_n/2,sigma_min)
        C_fail = 0
    end
    return sigma_n,C_success,C_fail
end

function select_evaluation_point_1D(new_points1,surr1::AbstractSurrogate,numb_iters,maxiters)
    v = [0.3,0.5,0.8,0.95]
    k = 4
    n = length(surr1.x)
    if mod(maxiters-1,4) != 0
        w_nR = v[mod(maxiters-1,4)]
    else
        w_nR = v[4]
    end
    w_nD = 1 - w_nR

    l = length(new_points1)
    evaluations1 = zeros(eltype(surr1.y[1]),l)

    for i = 1:l
        evaluations1[i] = surr1(new_points1[i])
    end
    s_max = maximum(evaluations1)
    s_min = minimum(evaluations1)
    V_nR = zeros(eltype(surr1.y[1]),l)
    for i = 1:l
        if abs(s_max-s_min) <= 10e-6
            V_nR[i] = 1.0
        else
            V_nR[i] = (evaluations1[i] - s_min)/(s_max-s_min)
        end
    end

    #Compute score V_nD
    V_nD = zeros(eltype(surr1.y[1]),l)
    delta_n_x = zeros(eltype(surr1.x[1]),l)
    delta = zeros(eltype(surr1.x[1]),n)
    for j = 1:l
        for i = 1:n
            delta[i] = norm(new_points1[j]-surr1.x[i])
        end
        delta_n_x[j] = minimum(delta)
    end
    delta_n_max = maximum(delta_n_x)
    delta_n_min = minimum(delta_n_x)
    for i = 1:l
        if abs(delta_n_max-delta_n_min) <= 10e-6
            V_nD[i] = 1.0
        else
            V_nD[i] = (delta_n_max - delta_n_x[i])/(delta_n_max-delta_n_min)
        end
    end

    #Compute weighted score
    W_n = w_nR*V_nR + w_nD*V_nD
    return new_points1[argmin(W_n)]
end

"""
surrogate_optimize(obj::Function,::DYCORS,lb::Number,ub::Number,surr1::AbstractSurrogate,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=100)

DYCORS optimization method in 1D, following closely: Combining radial basis function
surrogates and dynamic coordinate search in high-dimensional expensive black-box optimization".

"""
function surrogate_optimize(obj::Function,::DYCORS,lb::Number,ub::Number,surr1::AbstractSurrogate,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=100)
    x_best = argmin(surr1.y)
    y_best = minimum(surr1.y)
    sigma_n = 0.2*norm(ub-lb)
    d = length(lb)
    sigma_min = 0.2*(0.5)^6*norm(ub-lb)
    t_success = 3
    t_fail = max(d,5)
    C_success = 0
    C_fail = 0
    for k = 1:maxiters
        p_select = min(20/d,1)*(1-log(k))/log(maxiters-1)
        # In 1D I_perturb is always equal to one, no need to sample
        d = 1
        I_perturb = d
        new_points = zeros(eltype(surr1.x[1]),num_new_samples)
        for i = 1:num_new_samples
            new_points[i] = x_best + rand(Normal(0,sigma_n))
            while new_points[i] < lb || new_points[i] > ub
                if new_points[i] > ub
                    #reflection
                    new_points[i] = maximum(surr1.x) - norm(new_points[i] - maximum(surr1.x))
                end
                if new_points[i] < lb
                    #reflection
                    new_points[i] = minimum(surr1.x) + norm(new_points[i]-minimum(surr1.x))
                end
            end
        end

        x_new = select_evaluation_point_1D(new_points,surr1,k,maxiters)
        f_new = obj(x_new)

        if f_new < y_best
            C_success = C_success + 1
            C_fail = 0
        else
            C_fail = C_fail + 1
            C_success = 0
        end

        sigma_n,C_success,C_fail = adjust_step_size(sigma_n,sigma_min,C_success,t_success,C_fail,t_fail)

        if f_new < y_best
            x_best = x_new
            y_best = f_new
            add_point!(surr1,x_best,y_best)
        end
    end
    index = argmin(surr1.y)
    return (surr1.x[index],surr1.y[index])
end

function select_evaluation_point_ND(new_points,surrn::AbstractSurrogate,numb_iters,maxiters)
    v = [0.3,0.5,0.8,0.95]
    k = 4
    n = size(surrn.x,1)
    d = size(surrn.x,2)
    if mod(maxiters-1,4) != 0
        w_nR = v[mod(maxiters-1,4)]
    else
        w_nR = v[4]
    end
    w_nD = 1 - w_nR

    l = size(new_points,1)
    evaluations = zeros(eltype(surrn.y[1]),l)
    for i = 1:l
        evaluations[i] = surrn(Tuple(new_points[i,:]))
    end
    s_max = maximum(evaluations)
    s_min = minimum(evaluations)
    V_nR = zeros(eltype(surrn.y[1]),l)
    for i = 1:l
        if abs(s_max-s_min) <= 10e-6
            V_nR[i] = 1.0
        else
            V_nR[i] = (evaluations[i] - s_min)/(s_max-s_min)
        end
    end

    #Compute score V_nD
    V_nD = zeros(eltype(surrn.y[1]),l)
    delta_n_x = zeros(eltype(surrn.x[1]),l)
    delta = zeros(eltype(surrn.x[1]),n)
    for j = 1:l
        for i = 1:n
            delta[i] = norm(new_points[j,:]-collect(surrn.x[i]))
        end
        delta_n_x[j] = minimum(delta)
    end
    delta_n_max = maximum(delta_n_x)
    delta_n_min = minimum(delta_n_x)
    for i = 1:l
        if abs(delta_n_max-delta_n_min) <= 10e-6
            V_nD[i] = 1.0
        else
            V_nD[i] = (delta_n_max - delta_n_x[i])/(delta_n_max-delta_n_min)
        end
    end
    #Compute weighted score
    W_n = w_nR*V_nR + w_nD*V_nD
    return new_points[argmin(W_n),:]
end

"""
surrogate_optimize(obj::Function,::DYCORS,lb::Number,ub::Number,surr1::AbstractSurrogate,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=100)

This is an implementation of the DYCORS strategy by Regis and Shoemaker:

Rommel G Regis and Christine A Shoemaker.
Combining radial basis function surrogates and dynamic coordinate search in high-dimensional expensive black-box optimization.
Engineering Optimization, 45(5): 529–555, 2013.

This is an extension of the SRBF strategy that changes how the
candidate points are generated. The main idea is that many objective
functions depend only on a few directions so it may be advantageous to
perturb only a few directions. In particular, we use a perturbation probability
to perturb a given coordinate and decrease this probability after each function
evaluation so fewer coordinates are perturbed later in the optimization.
"""
function surrogate_optimize(obj::Function,::DYCORS,lb,ub,surrn::AbstractSurrogate,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=100)
    x_best = collect(surrn.x[argmin(surrn.y)])
    y_best = minimum(surrn.y)
    sigma_n = 0.2*norm(ub-lb)
    d = length(lb)
    sigma_min = 0.2*(0.5)^6*norm(ub-lb)
    t_success = 3
    t_fail = max(d,5)
    C_success = 0
    C_fail = 0
    for k = 1:maxiters
        p_select = min(20/d,1)*(1-log(k))/log(maxiters-1)
        new_points = zeros(eltype(surrn.x[1]),num_new_samples,d)
        for j = 1:num_new_samples
            w = sample(d,0,1,sample_type)
            I_perturb = w .< p_select
            if ~(true in I_perturb)
                val = rand(1:d)
                I_perturb = vcat(zeros(Int,val-1),1,zeros(Int,d-val))
            end
            I_perturb = Int.(I_perturb)
            for i = 1:d
                if I_perturb[i] == 1
                    new_points[j,i] = x_best[i] + rand(Normal(0,sigma_n))
                else
                    new_points[j,i] = x_best[i]
                end
            end
        end

        for i = 1:num_new_samples
            for j = 1:d
                while new_points[i,j] < lb[j] || new_points[i,j] > ub[j]
                    if new_points[i,j] > ub[j]
                        new_points[i,j] = maximum(surrn.x)[j] - norm(new_points[i,j] - maximum(surrn.x)[j])
                    end
                    if new_points[i,j] < lb[j]
                        new_points[i,j] = minimum(surrn.x)[j] + norm(new_points[i]-minimum(surrn.x)[j])
                    end
                end
            end
        end

        #ND version
        x_new = select_evaluation_point_ND(new_points,surrn,k,maxiters)
        f_new = obj(x_new)

        if f_new < y_best
            C_success = C_success + 1
            C_fail = 0
        else
            C_fail = C_fail + 1
            C_success = 0
        end

        sigma_n,C_success,C_fail = adjust_step_size(sigma_n,sigma_min,C_success,t_success,C_fail,t_fail)

        if f_new < y_best
            x_best = x_new
            y_best = f_new
            add_point!(surrn,Tuple(x_best),y_best)
        end
    end
    index = argmin(surrn.y)
    return (surrn.x[index],surrn.y[index])
end

function obj2_1D(value,points)
    min = +Inf
    my_p = filter(x->abs(x-value)>10^-6,points)
    for i = 1:length(my_p)
        new_val = norm(my_p[i]-value)
        if new_val < min
            min = new_val
        end
    end
    return min
end

function I_tier_ranking_1D(P,surrSOP::AbstractSurrogate)
    #obj1 = objective_function
    #obj2 = obj2_1D
    Fronts = Dict{Int,Array{eltype(surrSOP.x[1]),1}}()
    i = 1
    while true
        F = []
        j = 1
        for p in P
            n_p = 0
            k = 1
            for q in P
                #I use equality with floats because p and q are in surrSOP.x
                #for sure at this stage
                p_index = j
                q_index = k
                val1_p = surrSOP.y[p_index]
                val2_p = obj2_1D(p,P)
                val1_q = surrSOP.y[q_index]
                val2_q = obj2_1D(q,P)
                p_dominates_q = (val1_p < val1_q || abs(val1_p-val1_q) <= 10^-5) &&
                                (val2_p < val2_q || abs(val2_p-val2_q) <= 10^-5) &&
                                ((val1_p < val1_q) || (val2_p < val2_q))

                q_dominates_p = (val1_p < val1_q || abs(val1_p-val1_q) < 10^-5) &&
                                (val2_p < val2_q || abs(val2_p-val2_q) < 10^-5) &&
                                ((val1_p < val1_q) || (val2_p < val2_q))
                if q_dominates_p
                    n_p += 1
                end
                k = k + 1
            end
        if n_p == 0
            # no individual dominates p
            push!(F,p)
        end
        j = j + 1
        end
        if length(F) > 0
            Fronts[i] = F
            P = setdiff(P,F)
            i = i + 1
        else
            return Fronts
        end
    end
    return F
end

function II_tier_ranking_1D(D::Dict,srg::AbstractSurrogate)
    for i = 1:length(D)
        pos = []
        yn = []
        for j = 1:length(D[i])
            push!(pos,findall(e->e==D[i][j],srg.x))
            push!(yn,srg.y[pos[j]])
        end
        D[i] = D[i][sortperm(D[i])]
    end
    return D
end

function Hypervolume_Pareto_improving(f1_new,f2_new,Pareto_set)
    if size(Pareto_set,1) == 1
        area_before = zero(eltype(f1_new))
    else
        my_p = Pareto_set
        #Area before
        v_ref = [maximum(Pareto_set[:,1]),maximum(Pareto_set[:,2])]
        my_p = vcat(my_p,v_ref)
        v = sortperm(my_p[:,2])
        my_p[:,1] = my_p[:,1][v]
        my_p[:,2] = my_p[:,2][v]
        area_before = zero(eltype(f1_new))
        for j = 1:length(v)-1
            area_before += (my_p[j+1,2]-my_p[j,2])*(v_ref[1]-my_p[j])
        end
    end
    #Area after
    Pareto_set = vcat(Pareto_set,[f1_new f2_new])
    v_ref = [maximum(Pareto_set[:,1]) maximum(Pareto_set[:,2])]
    Pareto_set = vcat(Pareto_set,v_ref)
    v = sortperm(Pareto_set[:,2])
    Pareto_set[:,1] = Pareto_set[:,1][v]
    Pareto_set[:,2] = Pareto_set[:,2][v]
    area_after = zero(eltype(f1_new))
    for j = 1:length(v)-1
        area_after += (Pareto_set[j+1,2]-Pareto_set[j,2])*(v_ref[1]-Pareto_set[j])
    end
    return area_after - area_before
end

"""
surrogate_optimize(obj::Function,::SOP,lb::Number,ub::Number,surr::AbstractSurrogate,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=100)

SOP Surrogate optimization method, following closely the following papers:

    - SOP: parallel surrogate global optimization with Pareto center selection for computationally expensive single objective problems by Tipaluck Krityakierne
    - Multiobjective Optimization Using Evolutionary Algorithms by Kalyan Deb
#Suggested number of new_samples = min(500*d,5000)
"""
function surrogate_optimize(obj::Function,sop1::SOP,lb::Number,ub::Number,surrSOP::AbstractSurrogate,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=100)
    d = length(lb)
    N_fail = 3
    N_tenure = 5
    tau = 10^-5
    num_P = sop1.p
    centers_global = surrSOP.x
    r_centers_global = 0.2*norm(ub-lb)*ones(length(surrSOP.x))
    N_failures_global = zeros(length(surrSOP.x))
    tabu = []
    N_tenures_tabu = []
    for k = 1:maxiters
        N_tenures_tabu .+= 1
        #deleting points that have been in tabu for too long
        del = N_tenures_tabu .> N_tenure

        if length(del) > 0
            for i = length(del)
                if del[i]
                    del[i] = i
                end
            end
            deleteat!(N_tenures_tabu,del)
            deleteat!(tabu,del)
        end

        ##### P CENTERS ######
        C = []

        #S(x) set of points already evaluated
        #Rank points in S with:
        #1) Non dominated sorting
        Fronts_I = I_tier_ranking_1D(centers_global,surrSOP)
        #2) Second tier ranking
        Fronts = II_tier_ranking_1D(Fronts_I,surrSOP)
        ranked_list = []
        for i = 1:length(Fronts)
            for j = 1:length(Fronts[i])
                push!(ranked_list,Fronts[i][j])
            end
        end
        ranked_list = eltype(surrSOP.x[1]).(ranked_list)

        centers_full = 0
        i = 1
        while i <= length(ranked_list) && centers_full == 0
            flag = 0
            for j = 1:length(ranked_list)
                for m = 1:length(tabu)
                    if abs(ranked_list[j]-tabu[m]) < tau
                        flag = 1
                    end
                end
                for l = 1:length(centers_global)
                    if abs(ranked_list[j]-centers_global[l]) < tau
                        flag = 1
                    end
                end
            end
            if flag == 1
                skip
            else
                push!(C,ranked_list[i])
                if length(C) == num_P
                    centers_full = 1
                end
            end
            i = i + 1
        end

        #I examined all the points in the ranked list but num_selected < num_p
        #I just iterate again using only radius rule
        if length(C) < num_P
            i = 1
            while i <= length(ranked_list) && centers_full == 0
                flag = 0
                for j = 1:length(ranked_list)
                    for m = 1:length(centers_global)
                        if abs(centers_global[j] - ranked_list[m]) < tau
                            flag = 1
                        end
                    end
                end
                if flag == 1
                    skip
                else
                    push!(C,ranked_list[i])
                    if length(C) == num_P
                        centers_full = 1
                    end
                end
                i = i + 1
            end
        end

        #If I still have num_selected < num_P, I double down on some centers iteratively
        if length(C) < num_P
            i = 1
            while i <= length(ranked_list)
                push!(C,ranked_list[i])
                if length(C) == num_P
                    centers_full = 1
                end
                i = i + 1
            end
        end

        #Here I have selected C = [] containing the centers
        r_centers = 0.2*norm(ub-lb)*ones(num_P)
        N_failures = zeros(num_P)
        #2.3 Candidate search
        new_points = zeros(eltype(surrSOP.x[1]),num_P,2)
        for i = 1:num_P
            N_candidates = zeros(eltype(surrSOP.x[1]),num_new_samples)
            #Using phi(n) just like DYCORS, merit function = surrogate
            #Like in DYCORS, I_perturb = 1 always
            evaluations = zeros(eltype(surrSOP.y[1]),num_new_samples)
            for j = 1:num_new_samples
                a = lb - C[i]
                b = ub - C[i]
                N_candidates[j] = C[i] + rand(truncated(Normal(0,r_centers[i]), a,b))
                evaluations[j] = surrSOP(N_candidates[j])
            end
            x_best = N_candidates[argmin(evaluations)]
            y_best = minimum(evaluations)
            new_points[i,1] = x_best
            new_points[i,2] = y_best
        end

        #new_points[i] now contains:
        #[x_1,y_1; x_2,y_2,...,x_{num_new_samples},y_{num_new_samples}]

        #2.4 Adaptive learning and tabu archive
        for i=1:num_P
            if new_points[i,1] in centers_global
                r_centers[i] = r_centers_global[i]
                N_failures[i] = N_failures_global[i]
            end

            f_1 = new_points[i,1]
            f_2 = obj2_1D(f_1,surrSOP.x)

            l = length(Fronts[1])
            Pareto_set = zeros(eltype(surrSOP.x[1]),l,2)

            for j = 1:l
                val = obj2_1D(Fronts[1][j],surrSOP.x)
                Pareto_set[j,1] = Fronts[1][j]
                Pareto_set[j,2] = val
            end
            if (Hypervolume_Pareto_improving(f_1,f_2,Pareto_set)<tau)
                #failure
                r_centers[i] = r_centers[i]/2
                N_failures[i] += 1
                if N_failures[i] > N_fail
                    push!(tabu,C[i])
                    push!(N_tenures_tabu,0)
                end
            else
                #P_i is success
                #Adaptive_learning
                add_point!(surrSOP,new_points[i,1],new_points[i,2])
                push!(r_centers_global,r_centers[i])
                push!(N_failures_global,N_failures[i])
            end
        end
    end
    index = argmin(surrSOP.y)
    return (surrSOP.x[index],surrSOP.y[index])
end


text:
m = Float64[
    1 3 8 4;
    2 5 3 1;
    6 3 2 5;
]

f = Float64[
    2 1;
    1 3
]

expected = Float64[
    22 28 26;
    24 22 24
]

@test forward(ConvolutionLayer(2, 1, f), m) == expected

m = Float64[
    1 3 8;
    2 5 3;
    6 3 2;
]

f = Float64[
    2 1;
    1 3;
    1 0;
]
expected = Float64[
    28 31
]

output = forward(ConvolutionLayer(size(f, 1), size(f, 2), 1, 1, f), m)
@test output == expected

ext: jl

max_stars_repo_name: IshitaTakeshi/NeuralNetwork

max_issues_repo_name: IshitaTakeshi/NeuralNetwork

max_forks_repo_name: IshitaTakeshi/NeuralNetwork


text:
using Documenter, BangBang

makedocs(
    ;
    modules = [BangBang],
    format = Documenter.HTML(),
    pages = [
        "Home" => "index.md",
        hide("internals.md"),
    ],
    repo = "https://github.com/JuliaFolds/BangBang.jl/blob/{commit}{path}#L{line}",
    sitename = "BangBang.jl",
    authors = "Takafumi Arakaki <aka.tkf@gmail.com>",
    doctest = false,  # run via tests
    strict = [
        :autodocs_block,
        :cross_references,
        :docs_block,
        :doctest,
        :eval_block,
        :example_block,
        :footnote,
        :linkcheck,
        :meta_block,
        # Ignoring `:missing_docs` since Documenter does not handle docstrings
        # defined outside the sub-module nicely:
        # :missing_docs,
        :parse_error,
        :setup_block,
    ],
    # Ref:
    # https://juliadocs.github.io/Documenter.jl/stable/lib/public/#Documenter.makedocs
)

deploydocs(; repo = "github.com/JuliaFolds/BangBang.jl", push_preview = true)


text:
using BinDeps

@BinDeps.setup

@linux_only begin
    libnames = ["libMagickWand"]
    suffixes = ["", "-Q16", "-6.Q16", "-Q8"]
    options = ["", "HDRI"]
    aliases = vec(libnames.*transpose(suffixes).*reshape(options,(1,1,length(options))))
    libwand = library_dependency("libMagickWand", aliases=aliases)
    provides(AptGet, "libmagickwand4", libwand)
    provides(Yum, "ImageMagick", libwand)
end

@windows_only begin
    libwand = library_dependency("CORE_RL_wand_")

    const OS_ARCH = WORD_SIZE == 64 ? "x86_64" : "x86"

    if WORD_SIZE == 32
        provides(Binaries,URI("http://www.imagemagick.org/download/binaries/ImageMagick-6.8.8-6-Q16-windows-dll.exe"),libwand,os = :Windows)
    else
        provides(Binaries,URI("http://www.imagemagick.org/download/binaries/ImageMagick-6.8.8-6-Q16-windows-x64-dll.exe"),libwand,os = :Windows)
    end
end

@osx_only begin
    if Pkg.installed("Homebrew") === nothing
            error("Homebrew package not installed, please run Pkg.add(\"Homebrew\")")
    end
    using Homebrew
    libwand = library_dependency("libMagickWand-6.Q16")
    provides( Homebrew.HB, "imagemagick", libwand, os = :Darwin )
end

@BinDeps.install


text:
using CSV
using DataFrames
using StatsBase
using StatsModels
using Statistics
using Lathe
using GLM

data = DataFrame(CSV.File("motor_cena.csv"))

#display(describe(data))
#display(countmap(data[!,:Boja]))
#display(countmap(data[!,:Stanje]))
#display(countmap(data[!,:Ostecenje]))
#display(describe(data))

select!(data, Not(:Boja))

data[ismissing.(data[!,:Stanje]), :Stanje] .= mode(skipmissing(data[!,:Stanje]))
data[ismissing.(data[!,:Ostecenje]), :Ostecenje] .= mode(skipmissing(data[!,:Ostecenje]))
data[ismissing.(data[!,:Kilometraza]), :Kilometraza] .= trunc(Int64, mode(skipmissing(data[!,:Kilometraza])))
dropmissing!(data,[:Cena])
filter!(row -> row.Cena != 0, data)

filter!(row -> row.Kilometraza <= 500000, data)
filter!(row -> row.Godiste <= 2021 && row.Godiste >= 1900, data)
filter!(row -> row.BrojCilindara > 0 && row.BrojCilindara < 200, data)

covKwKs = cov(data.kW, data.KS)

if covKwKs > 0.6
    select!(data, Not(:kW))
end

fm = @formula(Cena ~ Stanje+Tip+Godiste+Kilometraza+Kubikaza+KS+BrojCilindara+Ostecenje)
dataTrain, dataTest = Lathe.preprocess.TrainTestSplit(data, .75)
linearRegressor = lm(fm, dataTrain)

predictTest = predict(linearRegressor, dataTest)
errorTest = dataTest.Cena - predictTest
rmseTest = sqrt(mean(errorTest.*errorTest))

predictTrain = predict(linearRegressor, dataTrain)
errorTrain = dataTrain.Cena - predictTrain
rmseTrain = sqrt(mean(errorTrain.*errorTrain))

println("RMSE trening skupa je $rmseTrain")
println("RMSE test skupa je: $rmseTest")

if r2(linearRegressor) > 0.5
    println("Model je dobar za predvidjanje!")
end

#plotlyjs()

#plotGodisteCena = scatter(data.Godiste, data.Cena, title = "Godiste - Cena", xlabel="Godiste", ylabel = "Cena")

#display(plotGodisteCena)
#savefig(plotGodisteCena, "godisteCena.html")

#plotKilometrazaCena = scatter(data.Kilometraza, data.Cena, title = "Kilometraza - Cena", xlabel="Kilometraza", ylabel = "Cena")

#display(plotKilometrazaCena)
#savefig(plotKilometrazaCena, "kilometrazaCena.html")

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Julia 1.1.0
#     language: julia
#     name: julia-1.1
# ---

# # This notebook
#
# - [x] read gridded output + `float_trajectory*data` => `uvetc` dictionnary.
# - [x] recompute `u,v` from gridded output
# - [x] compare with `u,v` from `float_traj*data`
# - [x] solve for float trajectory with `DifferentialEquations.jl`
#
# _Notes:_ For documentation see <https://gaelforget.github.io/MeshArrays.jl/stable/>, <https://docs.juliadiffeq.org/latest/solvers/ode_solve.html> and <https://en.wikipedia.org/wiki/Displacement_(vector)>

# ## 1. import software

using IndividualDisplacements, MeshArrays, DifferentialEquations, Plots
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"PlotIndDisp.jl"))

# ## 2. reload trajectories from `MITgcm/pkg/flt`

dirIn="flt_example/"
prec=Float32
df=IndividualDisplacements.ReadDisplacements(dirIn,prec)
PyPlot.figure(); PlotBasic(df,300,100000.0)

# ## 3. Read gridded variables via `MeshArrays.jl`
#
# Put grid variables in a dictionary.
#
# _Note:_ `myread` function deals with tiled files from `flt_example/`.

# +
import IndividualDisplacements: myread

mygrid=gcmgrid("flt_example/","ll",1,[(80,42)], [80 42], Float32, read, write)
nr=8

GridVariables=Dict("XC" => myread(mygrid.path*"XC",MeshArray(mygrid,Float32)),
"YC" => myread(mygrid.path*"YC",MeshArray(mygrid,Float32)),
"XG" => myread(mygrid.path*"XG",MeshArray(mygrid,Float32)),
"YG" => myread(mygrid.path*"YG",MeshArray(mygrid,Float32)),
"dx" => 5000.0);
# -

# Put velocity fields in in another dictionary and merge the two dictionaries.

# +
t0=0.0 #approximation / simplification
t1=18001.0*3600.0

u0=myread(mygrid.path*"U.0000000001",MeshArray(mygrid,Float32,nr))
u1=myread(mygrid.path*"U.0000018001",MeshArray(mygrid,Float32,nr))
v0=myread(mygrid.path*"V.0000000001",MeshArray(mygrid,Float32,nr))
v1=myread(mygrid.path*"V.0000018001",MeshArray(mygrid,Float32,nr))

kk=3 #3 to match -1406.25 in pkg/flt output
u0=u0[:,kk]; u1=u1[:,kk];
v0=v0[:,kk]; v1=v1[:,kk];

u0=u0./GridVariables["dx"]
u1=u1./GridVariables["dx"]
v0=v0./GridVariables["dx"]
v1=v1./GridVariables["dx"]

uvt = Dict("u0" => u0, "u1" => u1, "v0" => v0, "v1" => v1, "t0" => t0, "t1" => t1)

uvetc=merge(uvt,GridVariables);
# -

# ## 4. Visualize velocity fields

# +
mskW=myread(mygrid.path*"hFacW",MeshArray(mygrid,Float32,nr))
mskW=1.0 .+ 0.0 * mask(mskW[:,kk],NaN,0.0)
mskS=myread(mygrid.path*"hFacS",MeshArray(mygrid,Float32,nr))
mskS=1.0 .+ 0.0 * mask(mskS[:,kk],NaN,0.0)

msk=Dict("mskW" => mskW, "mskS" => mskS)

uvetc=merge(uvetc,msk);
# -

heatmap(mskW[1,1].*u0[1,1],title="U at the start")

heatmap(mskS[1,1].*v0[1,1],title="V at the start")

heatmap(mskW[1,1].*u1[1,1]-u0[1,1],title="U end - U start")

# ## 5. Visualize trajectories from `MITgcm/pkg/flt`
#
# Select one trajectory

tmp=df[df.ID .== 200, :]
tmp[1:4,:]

# Super-impose trajectory over velocity field (first for u ...)

PyPlot.contourf(GridVariables["XG"].f[1], GridVariables["YC"].f[1], mskW.f[1].*u0.f[1])
PyPlot.plot(tmp[:,:lon],tmp[:,:lat],color="r")
colorbar()

# Super-impose trajectory over velocity field (... then for v)

PyPlot.contourf(GridVariables["XG"].f[1], GridVariables["YC"].f[1], mskS.f[1].*v0.f[1])
PyPlot.plot(tmp[:,:lon],tmp[:,:lat],color="r")
colorbar()

# ## 6. Recompute displacements from gridded flow fields

# +
comp_vel=IndividualDisplacements.VelComp
get_vel=IndividualDisplacements.VelCopy

uInit=[tmp[1,:lon];tmp[1,:lat]]./uvetc["dx"]
nSteps=Int32(tmp[end,:time]/3600)-2
du=fill(0.0,2);
# -

# Visualize and compare with actual grid point values -- jumps on the tangential component are expected with linear scheme:

tmpu=fill(0.0,100)
tmpv=fill(0.0,100)
tmpx=fill(0.0,100)
for i=1:100
    tmpx[i]=500.0 *i./uvetc["dx"]
    comp_vel(du,[tmpx[i];0.499./uvetc["dx"]],uvetc,0.0)
    tmpu[i]=du[1]
    tmpv[i]=du[2]
end
Plots.plot(tmpx,tmpu)
Plots.plot!(uvetc["XG"].f[1][1:10,1]./uvetc["dx"],uvetc["u0"].f[1][1:10,1],marker=".")
Plots.plot!(tmpx,tmpv)
Plots.plot!(uvetc["XG"].f[1][1:10,1]./uvetc["dx"],uvetc["v0"].f[1][1:10,1],marker=".")

tmpu=fill(0.0,100)
tmpv=fill(0.0,100)
tmpy=fill(0.0,100)
for i=1:100
    tmpy[i]=500.0 *i./uvetc["dx"]
    comp_vel(du,[0.499./uvetc["dx"];tmpy[i]],uvetc,0.0)
    tmpu[i]=du[1]
    tmpv[i]=du[2]
end
Plots.plot(tmpx,tmpu)
Plots.plot!(uvetc["YG"].f[1][1,1:10]./uvetc["dx"],uvetc["u0"].f[1][1,1:10],marker=".")
Plots.plot!(tmpx,tmpv)
Plots.plot!(uvetc["YG"].f[1][1,1:10]./uvetc["dx"],uvetc["v0"].f[1][1,1:10],marker=".")

# Compare recomputed velocities with those from `pkg/flt`

nSteps=2998
tmpu=fill(0.0,nSteps); tmpv=fill(0.0,nSteps);
tmpx=fill(0.0,nSteps); tmpy=fill(0.0,nSteps);
refu=fill(0.0,nSteps); refv=fill(0.0,nSteps);
for i=1:nSteps
    get_vel(du,[tmp[i,:lon],tmp[i,:lat]],tmp,tmp[i,:time])
    refu[i]=du[1]./uvetc["dx"]
    refv[i]=du[2]./uvetc["dx"]
    comp_vel(du,[tmp[i,:lon],tmp[i,:lat]]./uvetc["dx"],uvetc,tmp[i,:time])
    tmpu[i]=du[1]
    tmpv[i]=du[2]
end
#
Plots.plot(tmpu)
Plots.plot!(tmpv)
Plots.plot!(refu)
Plots.plot!(refv)

# ## 6. Recompute trajectories from gridded flow fields
#
# Solve through time using `DifferentialEquations.jl`

using DifferentialEquations
tspan = (0.0,nSteps*3600.0)
#prob = ODEProblem(get_vel,uInit,tspan,tmp)
prob = ODEProblem(comp_vel,uInit,tspan,uvetc)
sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
sol[1:4]

# Compare recomputed trajectories with those from `pkg/flt`

# +
ref=transpose([tmp[1:nSteps,:lon] tmp[1:nSteps,:lat]])
maxLon=80*5.e3
maxLat=42*5.e3
show(size(ref))
for i=1:nSteps-1
    ref[1,i+1]-ref[1,i]>maxLon/2 ? ref[1,i+1:end]-=fill(maxLon,(nSteps-i)) : nothing
    ref[1,i+1]-ref[1,i]<-maxLon/2 ? ref[1,i+1:end]+=fill(maxLon,(nSteps-i)) : nothing
    ref[2,i+1]-ref[2,i]>maxLat/2 ? ref[2,i+1:end]-=fill(maxLat,(nSteps-i)) : nothing
    ref[2,i+1]-ref[2,i]<-maxLat/2 ? ref[2,i+1:end]+=fill(maxLat,(nSteps-i)) : nothing
end
ref=ref./uvetc["dx"]

using Plots
Plots.plot(sol[1,:],sol[2,:],linewidth=5,title="Using Recomputed Velocities",
     xaxis="lon",yaxis="lat",label="Julia Solution") # legend=false
Plots.plot!(ref[1,:],ref[2,:],lw=3,ls=:dash,label="MITgcm Solution")
# -


text:
function bracket(::Type{UnivariateSolverUtils}, arg0::RealFieldUnivariateFunction, arg1::RealFieldElement, arg2::RealFieldElement, arg3::RealFieldElement)
    return jcall(UnivariateSolverUtils, "bracket", Vector{RealFieldElement}, (RealFieldUnivariateFunction, RealFieldElement, RealFieldElement, RealFieldElement), arg0, arg1, arg2, arg3)
end

function bracket(::Type{UnivariateSolverUtils}, arg0::RealFieldUnivariateFunction, arg1::RealFieldElement, arg2::RealFieldElement, arg3::RealFieldElement, arg4::RealFieldElement, arg5::RealFieldElement, arg6::jint)
    return jcall(UnivariateSolverUtils, "bracket", Vector{RealFieldElement}, (RealFieldUnivariateFunction, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, jint), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function bracket(::Type{UnivariateSolverUtils}, arg0::RealFieldUnivariateFunction, arg1::RealFieldElement, arg2::RealFieldElement, arg3::RealFieldElement, arg4::jint)
    return jcall(UnivariateSolverUtils, "bracket", Vector{RealFieldElement}, (RealFieldUnivariateFunction, RealFieldElement, RealFieldElement, RealFieldElement, jint), arg0, arg1, arg2, arg3, arg4)
end

function bracket(::Type{UnivariateSolverUtils}, arg0::UnivariateFunction, arg1::jdouble, arg2::jdouble, arg3::jdouble)
    return jcall(UnivariateSolverUtils, "bracket", Vector{jdouble}, (UnivariateFunction, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function bracket(::Type{UnivariateSolverUtils}, arg0::UnivariateFunction, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::jint)
    return jcall(UnivariateSolverUtils, "bracket", Vector{jdouble}, (UnivariateFunction, jdouble, jdouble, jdouble, jdouble, jdouble, jint), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function bracket(::Type{UnivariateSolverUtils}, arg0::UnivariateFunction, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::jint)
    return jcall(UnivariateSolverUtils, "bracket", Vector{jdouble}, (UnivariateFunction, jdouble, jdouble, jdouble, jint), arg0, arg1, arg2, arg3, arg4)
end

function force_side(::Type{UnivariateSolverUtils}, arg0::jint, arg1::UnivariateFunction, arg2::BracketedUnivariateSolver, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::AllowedSolution)
    return jcall(UnivariateSolverUtils, "forceSide", jdouble, (jint, UnivariateFunction, BracketedUnivariateSolver, jdouble, jdouble, jdouble, AllowedSolution), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function is_bracketing(::Type{UnivariateSolverUtils}, arg0::UnivariateFunction, arg1::jdouble, arg2::jdouble)
    return jcall(UnivariateSolverUtils, "isBracketing", jboolean, (UnivariateFunction, jdouble, jdouble), arg0, arg1, arg2)
end

function is_sequence(::Type{UnivariateSolverUtils}, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(UnivariateSolverUtils, "isSequence", jboolean, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function midpoint(::Type{UnivariateSolverUtils}, arg0::jdouble, arg1::jdouble)
    return jcall(UnivariateSolverUtils, "midpoint", jdouble, (jdouble, jdouble), arg0, arg1)
end

function solve(::Type{UnivariateSolverUtils}, arg0::UnivariateFunction, arg1::jdouble, arg2::jdouble)
    return jcall(UnivariateSolverUtils, "solve", jdouble, (UnivariateFunction, jdouble, jdouble), arg0, arg1, arg2)
end

function solve(::Type{UnivariateSolverUtils}, arg0::UnivariateFunction, arg1::jdouble, arg2::jdouble, arg3::jdouble)
    return jcall(UnivariateSolverUtils, "solve", jdouble, (UnivariateFunction, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function verify_bracketing(::Type{UnivariateSolverUtils}, arg0::UnivariateFunction, arg1::jdouble, arg2::jdouble)
    return jcall(UnivariateSolverUtils, "verifyBracketing", void, (UnivariateFunction, jdouble, jdouble), arg0, arg1, arg2)
end

function verify_interval(::Type{UnivariateSolverUtils}, arg0::jdouble, arg1::jdouble)
    return jcall(UnivariateSolverUtils, "verifyInterval", void, (jdouble, jdouble), arg0, arg1)
end

function verify_sequence(::Type{UnivariateSolverUtils}, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(UnivariateSolverUtils, "verifySequence", void, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end


text:
using Pluto

host = get(ENV, "HOST", "127.0.0.1")
notebook = "./src/notebook.jl"
if !isempty(host)
    Pluto.run(notebook = notebook, host = host)
else
    Pluto.run(notebook = notebook)
end


text:
using Test
using Distances: Euclidean, CosineDist
using SparseArrays
using LinearAlgebra
using UMAP
using UMAP: fuzzy_simplicial_set, compute_membership_strengths, smooth_knn_dists, smooth_knn_dist, spectral_layout, optimize_embedding, knn_search, combine_fuzzy_sets, fit_ab, SMOOTH_K_TOLERANCE

include("utils_tests.jl")
include("umap_tests.jl")


text:
# # Writing TIFFs

# This page is a tutorial for saving TIFFs using TiffImages.jl and covers some common
# use cases

#md # ```@contents
#md # Pages = ["writing.md"]
#md # Depth = 5
#md # ```

# You might want to write TIFFs to disk too. Now this can be done quite simply
# with TiffImages.jl. Say you have some AbstractArray type that you want to save, here
# we'll call it `data`:

using Random
using Images # for nice inline images

Random.seed!(123)
data = rand(RGB{N0f8}, 10, 10)

#md # !!! note
#md #     TiffImages.jl only works with AbstractArrays with `eltype`s of `<:Colorant` because
#md #     the writer needs to know how to represent the image data on disk. Make sure to
#md #     convert your `AbstractArrays` using before passing them. See the
#md #     [common strategies](#Strategies-for-saving-common-types) section
#md #     below for tips.

# ## Converting to `TiffImages.jl`'s TIFF type
# Next lets convert `data` to a TIFF type

using TiffImages
img = TiffImages.DenseTaggedImage(data)

# Wait nothing happened! Hang with me, lets take a closer look at our new object
# using the `dump` command. We can see that there's now new information
# associated with our data! TiffImages.jl usually represents TIFF images as simply the
# data and associated tags that describe the data

dump(img; maxdepth=1)

# The tags are organized as a vector of what are called [Image File
# Directories](https://www.awaresystems.be/imaging/tiff/faq.html) (IFDs). For a
# simple 2D image like what we have, the IFDs will be stored a vector of
# length=1. For 3D images, the length of the IFDs vector will equal the length
# of the image in the third dimension.

# Lets take a look at what tags there are:

ifd = first(img.ifds) # since our data is 2D
ifd

# ## Manipulating TIFF Tags
#
# These are some of the most basic tags that are required by the TIFF spec. We
# can even update it to add our own custom tags

ifd[TiffImages.IMAGEDESCRIPTION] = "This is very important data"
ifd

# We can even add tags that aren't in the standard set in
# [`TiffImages.TiffTag`](@ref) as long as they are a `UInt16`

ifd[UInt16(34735)] = UInt16[1, 2, 3]
ifd

# We can also delete tags if we decide we don't want them:

delete!(ifd, TiffImages.IMAGEDESCRIPTION)
ifd

#md # !!! warning
#md #     Careful with `delete!`, if any of core tags are deleted, TiffImages.jl and
#md #     other readers might fail to read the file

# ## Saving to disk
#
# Once you're happy with your TIFF object, you can write it to disk as follows:

open("test.tif", "w") do io
    write(io, img)
end

# And to just double check, we can load it right back in

TiffImages.load("test.tif")

# ## Strategies for saving common types

# The general strategy for saving arrays will differ a bit depending on the
# type. The key step is the convert or reinterpret the arrays so that the
# elements are subtypes of `Colors.Colorant`

# #### Unsigned Integers

# Say you want to save a 3D array of small integers as grayscale values.

data2 = rand(UInt8.(1:255), 5, 10)
eltype(data2)

# You can't directly save the `data2` since TiffImages.jl needs some color information
# to properly save the file. You can use
# [`reinterpret`](https://docs.julialang.org/en/v1/base/arrays/#Base.reinterpret)
# to accomplish this:

grays = reinterpret(Gray{N0f8}, data2)
img2 = TiffImages.DenseTaggedImage(grays)

# Here the data are first reinterpreted as `N0f8`s, which is a
# [`FixedPointNumber`](https://github.com/JuliaMath/FixedPointNumbers.jl) then
# wrapped with a Gray type that marks this as a grayscale image. TiffImages.jl uses
# this information to update the TIFF tags

# #### Floating point numbers

# With RGB we can reinterpret the first dimension of a 3D array as the 3
# different color components (red, green, and blue):

data = rand(Float64, 3, 5, 10);
colors = dropdims(reinterpret(RGB{eltype(data)}, data), dims=1) # drop first dimension
img3 = TiffImages.DenseTaggedImage(colors)

# Here we dropped the first dimension since it was collapsed into the RGB type
# when we ran the `reinterpret` command.

# #### Signed integers

# Say you want to save data that has negative integer values. In that case, you
# can't use `N0f8`, etc because those only worked for unsigned integers. You
# have to instead use `Q0f63`, etc, which is a different kind of fixed point
# number that uses one bit for the sign info (that's why it's `Q0f63`, not `Q0f64`!)

data = rand(-100:100, 5, 5)
#--------------------------
img4 = TiffImages.DenseTaggedImage(reinterpret(Gray{Q0f63}, data))
println(img4.ifds[1])

# As you can see the `SAMPLEFORMATS` and `BITSPERSAMPLE` tags correctly updated
# to show that this TIFF contains signed integers and 64-bit data, respectively.

#md # !!! warning
#md #     Currently, several of the display libraries struggle with showing
#md #     `Colorant`s backed by a signed type so you might run into errors, but
#md #     the data will still save properly


text:
using IterativeSolvers
using FactCheck

facts("rsvd") do

context("Small wide rectangular") do
    A = [1. 2 3 4; 5 6 7 8]
    S1 = svdfact(A)
    S2 = rsvd(A, 2, 0)

    @fact vecnorm(abs(S1[:U]) - abs(S2[:U])) --> less_than( √(eps()))
    @fact vecnorm(abs(S1[:Vt]) - abs(S2[:Vt])) --> less_than(√(eps()))
    @fact norm(S1[:S] - S2[:S]) --> less_than(√(eps()))
end

context("Small tall rectangular") do
    A = [1. 2; 3 4; 5 6; 7 8]
    S1 = svdfact(A)
    S2 = rsvd(A, 4, 0)

    @fact vecnorm(abs(S1[:U]) - abs(S2[:U])) --> less_than(√(eps()))
    @fact vecnorm(abs(S1[:Vt]) - abs(S2[:Vt])) --> less_than(√(eps()))
    @fact norm(S1[:S] - S2[:S]) --> less_than(√(eps()))
end

context("Small square") do
    A = [1. 2 3; 4 5 6; 7 8 9]
    S1 = svdfact(A)
    S2 = rsvd(A, 3, 0)

    @fact vecnorm(abs(S1[:U]) - abs(S2[:U]))  --> less_than(√(eps()))
    @fact vecnorm(abs(S1[:Vt]) - abs(S2[:Vt]))  --> less_than(√(eps()))
    @fact norm(S1[:S] - S2[:S])  --> less_than(√(eps()))
end

context("rrange_adaptive") do
    A = [1. 2 3; 4 5 6; 7 8 9]
    @fact size(IterativeSolvers.rrange_adaptive(A, 3, 1e-3)) --> (3,2)
end

end


text:
include("test_SimulationExamples.jl")
include("test_withoutMacros.jl")
include("test_withoutMacros_withoutVariables1.jl")
include("test_withoutMacros_withoutVariables2.jl")


text:
using Gadfly, Distributions, Random

set_default_plot_size(6inch, 16inch)

Random.seed!(1234)

x = rand(Normal(), 100)
y = rand(Normal(10), 100)

xd = Normal()
yd = Normal(10)

# two numeric vectors
pl1 = plot(x=x, y=y, Stat.qq, Geom.point)

# one numeric and one Distribution
pl2 = plot(x=x, y=yd, Stat.qq, Geom.point)
pl3 = plot(x=xd, y=y, Stat.qq, Geom.point)

# Apply different scales to x and y
pl4 = plot(x=x, y=exp.(y), Stat.qq, Geom.point, Scale.y_log10)
pl5= plot(x=exp.(x), y=y, Stat.qq, Geom.point, Scale.x_log10)

# Apply scales to Distributions
z = rand(Exponential(), 100)
pl5 = plot(x=z, y=Exponential(), Stat.qq, Geom.point)
pl6 = plot(x=log.(z), y=Exponential(), Stat.qq, Geom.point)
pl7 = plot(x=log.(z), y=Exponential(), Stat.qq, Geom.point, Scale.y_log)
pl8 = plot(x=z, y=Exponential(), Stat.qq, Geom.point, Scale.x_log, Scale.y_log)

# by analogy with Stat.func, computed (.y) aes should be scaled:
#plot(xmin=[1], xmax=[15], y=[exp], Stat.func, Geom.line)
# (y value is log(exp(x), or y=x)
#plot(xmin=[1], xmax=[15], y=[exp], Stat.func, Geom.line, Scale.y_log)

vstack(pl1, pl2, pl3, pl4, pl5, pl6, pl7, pl8)


text:
@testset "umap tests" begin
    A = sprand(10000, 10000, 0.001)
    B = dropzeros(A + A' - A .* A')

    @testset "constructor" begin
        @testset "argument validation tests" begin
            data = rand(5, 10)
            @test_throws ArgumentError UMAP_([1. 1.]; n_neighbors=0) # n_neighbors error
            @test_throws ArgumentError UMAP_([1. 1.], 0; n_neighbors=1) # n_comps error
            @test_throws ArgumentError UMAP_([1. 1.], 2; n_neighbors=1) # n_comps error
            @test_throws ArgumentError UMAP_([1. 1.; 1. 1.; 1. 1.];
                    n_neighbors=1, min_dist = 0.) # min_dist error
        end
    end

    @testset "input type stability tests" begin
        data = rand(5, 100)
        umap_ = UMAP_(data)
        @test umap_ isa UMAP_{Float64}
        @test size(umap_.graph) == (100, 100)
        @test size(umap_.embedding) == (2, 100)

        data = rand(Float32, 5, 100)
        @test UMAP_(data) isa UMAP_{Float32}
    end

    @testset "fuzzy_simpl_set" begin
        data = rand(20, 500)
        k = 5
        umap_graph = fuzzy_simplicial_set(data, k, Euclidean(), 1, 1.)
        @test issymmetric(umap_graph)
        @test all(0. .<= umap_graph .<= 1.)
        data = rand(Float32, 20, 500)
        umap_graph = fuzzy_simplicial_set(data, k, Euclidean(), 1, 1.f0)
        @test issymmetric(umap_graph)
        @test eltype(umap_graph) == Float32

        data = 2 .* rand(20, 1000) .- 1
        umap_graph = fuzzy_simplicial_set(data, k, CosineDist(), 1, 1.)
        @test issymmetric(umap_graph)
        @test all(0. .<= umap_graph .<= 1.)
    end

    @testset "smooth_knn_dists" begin
        dists = [0., 1., 2., 3., 4., 5.]
        rho = 1
        k = 6
        local_connectivity = 1
        bandwidth = 1.
        niter = 64
        sigma = smooth_knn_dist(dists, rho, k, bandwidth, niter)
        psum(ds, r, s) = sum(exp.(-max.(ds .- r, 0.) ./ s))
        @test psum(dists, rho, sigma) - log2(k)*bandwidth < SMOOTH_K_TOLERANCE

        knn_dists = [0. 0. 0.;
                     1. 2. 3.;
                     2. 4. 5.;
                     3. 4. 5.;
                     4. 6. 6.;
                     5. 6. 10.]
        rhos, sigmas = smooth_knn_dists(knn_dists, k, local_connectivity)
        @test rhos == [1., 2., 3.]
        diffs = [psum(knn_dists[:,i], rhos[i], sigmas[i]) for i in 1:3] .- log2(6)
        @test all(diffs .< SMOOTH_K_TOLERANCE)

        knn_dists = [0. 0. 0.;
                     0. 1. 2.;
                     0. 2. 3.]
        rhos, sigmas = smooth_knn_dists(knn_dists, 2, 1)
        @test rhos == [0., 1., 2.]

        rhos, sigmas = smooth_knn_dists(knn_dists, 2, 1.5)
        @test rhos == [0., 1.5, 2.5]
    end

    @testset "compute_membership_strengths" begin
        knns = [1 2 3; 2 1 2]
        dists = [0. 0. 0.; 2. 2. 3.]
        rhos = [2., 1., 4.]
        sigmas = [1., 1., 1.]
        true_rows = [1, 2, 2, 1, 3, 2]
        true_cols = [1, 1, 2, 2, 3, 3]
        true_vals = [0., 1., 0., exp(-1.), 0., 1.]
        rows, cols, vals = compute_membership_strengths(knns, dists, rhos, sigmas)
        @test rows == true_rows
        @test cols == true_cols
        @test vals == true_vals
    end

    @testset "optimize_embedding" begin
        layout = spectral_layout(B, 5)
        n_epochs = 1
        initial_alpha = 1.
        min_dist = 1.
        spread = 1.
        gamma = 1.
        neg_sample_rate = 5
        embedding = optimize_embedding(B, layout, n_epochs, initial_alpha, min_dist, spread, gamma, neg_sample_rate)
        @test embedding isa Array{Float64, 2}
    end

    @testset "spectral_layout" begin
        layout = spectral_layout(B, 5)
        @test layout isa Array{Float64, 2}
        @inferred spectral_layout(B, 5)
        layout32 = spectral_layout(convert(SparseMatrixCSC{Float32}, B), 5)
        @test layout32 isa Array{Float32, 2}
    end

end


text:
import TimeSeries

include("../data/data_5bus.jl")
include("../data/data_14bus.jl")

@testset "Test System constructors" begin
    tPowerSystem = System()

    sys5 = System(nodes5, generators5, loads5_DA, nothing, nothing,  1000.0)
    sys5 = System(nodes5, generators5, loads5_DA, branches5, nothing,  1000.0)

    battery5=[GenericBattery(name="Bat",
                             status=true,
                             bus=nodes5[2],
                             activepower=10.0,
                             energy=5.0,
                             capacity=(min=0.0, max=0.0),
                             inputactivepowerlimits=(min=0.0, max=50.0),
                             outputactivepowerlimits=(min=0.0, max=50.0),
                             efficiency=(in=0.90, out=0.80),
                            )]

    sys5b = System(nodes5, generators5, loads5_DA, nothing, battery5,  1000.0)

    generators_hg5 = [
        HydroFix("HydroFix", true, nodes5[2],
                 TechHydro(60.0, 15.0, (min=0.0, max=60.0), nothing, nothing, nothing,
                           nothing),
                 TimeSeries.TimeArray(DayAhead, solar_ts_DA)
        ),
        HydroCurtailment("HydroCurtailment", true, nodes5[3],
                         TechHydro(60.0, 10.0, (min=0.0, max=60.0), nothing, nothing,
                                   (up=10.0, down=10.0), nothing),
                         1000.0, TimeSeries.TimeArray(DayAhead, wind_ts_DA))
    ]

    sys5bh = System(nodes5, append!(generators5, generators_hg5), loads5_DA, branches5,
                         battery5,  1000.0)

     #Test Data for 14 Bus

    sys14 = System(nodes14, generators14, loads14, nothing, nothing, 1000.0)
    sys14 = System(nodes14, generators14, loads14, branches14, nothing, 1000.0)

    battery14 = [GenericBattery(name="Bat",
                                status=true,
                                bus=nodes14[2],
                                activepower=10.0,
                                energy=5.0,
                                capacity=(min=0.0, max=0.0),
                                inputactivepowerlimits=(min=0.0, max=50.0),
                                outputactivepowerlimits=(min=0.0, max=50.0),
                                efficiency=(in=0.90, out=0.80),
                               )]

    sys14b = System(nodes14, generators14, loads14, nothing, battery14, 1000.0)
    sys14b = System(nodes14, generators14, loads14, branches14, battery14, 1000.0)

    ps_dict = PowerSystems.parsestandardfiles(joinpath(MATPOWER_DIR, "case5_re.m"))
    sys = PowerSystems.System(ps_dict)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 5a14baba-4b8a-4917-b354-246775bb8c1d
#hide
using PlutoUI

# ╔═╡ 00f8bf47-47d2-4a99-86ba-d7cd0b7f268b
md""" ## First steps

Create a directory in your local system where you would want to start your work. For example, in your local system

$(LocalResource("../docs/images/testfolder.png"))

`cd` into this directory now. This will be the place where you are going to place your project. Now open a command prompt from inside this directory and type `julia`, this will launch a Julia shell. Here type `]` to launch the "pkg" prompt and then you can type `activate .` in order to create a restricted environment where you can restrict further installations or else you can skip typing `activate .` and install Nbdev by with the following command --> `add Nbdev`

$(LocalResource("../docs/images/addnbdev.png"))

Now go back to the Julia prompt and type `import Nbdev followed by Nbdev.new()`. This will create a new directory with "newproj" as the default name at the root with the following structure.

$(LocalResource("../docs/images/folderstruct.png"))

If you want have a different name for your project directory then give a name while executing the `new` command. For example `Nbdev.new("hellowproject")`.

"""

# ╔═╡ fbb6eadc-80a1-4cd9-ad27-bebcdd6e5981
md"""
!!! note

*Right now nbdev is highly opinionated about the folder structure like presence of "nbs", "src" and "docs" directory.*"""

# ╔═╡ b0c1fe8f-e796-4835-95f8-17990c7421a7
md"## First notebook"

# ╔═╡ e0694cd4-7d07-494d-a010-f705d1cbeb2c
md"""`cd` into the project directory that you created in previous step. this is where you will create your package and put all the related code.

Once you are inside this directory, open a terminal and start Pluto the usual way as you use to do and create a notebook inside the "nbs" directory.

At this point you might have installed Nbdev in your global environment. If this is not the case then follow any of the below methods."""

# ╔═╡ 15a5c4fd-301e-4f7c-9ad3-28ae2d214224
md"""

### Method 1 - Keeping Nbbev installation isolated for a specific Pluto notebook

* After you `cd` into your project directory (created with the command `Nbdev.new()`), open a terminal and type in the command `julia`.

* Then type in `]` to opne the pkg teminal.

* Activate the package environment by typing `activate .` . This will activate an isolated environment for you current package.

* In the Pkg terminal type `add Pluto`. This will install pluto at project level.

* After installation is done, move back to the Julia terminal and type in the following 👇

```
import Pluto
Pluto.run()
```

* Copy the notebook link which pops up in the terminal and paste it in the browser of choice. This will launch the Pluto GUI. Now, create a new notebook.

* In the first cell of your notebook type in the following code 👇

```
using Nbdev
```

* This will install Nbdev in our notebook environment. If you have followed this step then it's obvious that you want your developement environment to be isolated within a notebook. If this is the case then you can skip "method 2" and proceed with the next steps.
"""

# ╔═╡ 305c50dd-78bd-43d9-8f4e-7074bb046a59
md"""
!!! note
Do remember that if you keep your Nbdev installation isolated within a specific notebook then you may not be able to use some of the Nbdev commands at project level (outside a notebook i.e. from the terminal of your project). For example if you go with method 1 then you can run `Nbdev.build()` command from a terminal in your project root since Nbdev is not present in the Project.toml in your specific project. To do this you would need to install Nbdev at your project or global level anyways.

If you want to be able to use all the commands in Nbdev then use method 2.
"""

# ╔═╡ 2ac2c281-ad4a-4364-8f62-5c8202d4e9c1
md"### Method 2- Installing Nbdev at the project enviroment level and sharing the environment with all notebooks

* After you `cd` into your project directory (created with the command `Nbdev.new()`), open a terminal and type in the command `julia`.

* Then type in `]` to opne the pkg teminal.

* Activate the package environment by typing `activate .` . This will activate an isolated environment for you current package.

* Now,type add Nbdev. This wil install Nbdev at your project level.

* To share this installation of Nbdev (or any other dependencies which you had installed at your project level) use the following steps

  * Open a notebook inside the nbs directory.
  * In the first cell of your notebook use the following code
```
     begin
     import Pkg
     Pkg.activate(Base.current_project())
     end
```"

# ╔═╡ 9258500a-0a33-41cd-b32a-e09e820f84ee
md"## Start building your project"

# ╔═╡ dedef829-e43c-41c9-8238-b9d631dccb35
md"There are some anchors which tell Nbdev about the portion of your notebook that should be exported as code or should be part of documentation or should not be touched."

# ╔═╡ 99d43364-5ac6-43bb-a70c-20a79a463fdf
md"### #export"

# ╔═╡ fb11c61e-2b26-4b2a-9a34-835163ac7f87
md"when you put the comment #export at the begining of a cell, it tells nbdev to export that cell as a code into the source file."

# ╔═╡ c9f87989-522c-4adc-9a50-af95e3fa9cfa
md"
```
#export
function somefunc(x)
	x+1
end
```"

# ╔═╡ 0d4fb029-bfc8-4208-a6f4-bd5ea981ff82
md"Nbdev will now put the above cell as part of the source file (with the same name as the current notebook) inside the src directory."

# ╔═╡ ac65836e-6f4f-4212-b7f4-a648de9c761b
md"### #hide"

# ╔═╡ 789ebdb8-0017-4662-9120-a780584f8368
md"Suppose that you have included soem cells as par of experment while developing your package but you want that the cell should neither be exported as a code nor that should be exported as part of documentation.

In such a case use the `#hide` anchor at the begining of the cell."

# ╔═╡ 809102dd-b91a-458f-b5c7-2a49c7536d3b
md"
```
#hide
myexperiment = 1+1
```"

# ╔═╡ f8391a08-edf6-499c-8e8d-8f731dad8fd4
md"such a cell won't be exported."

# ╔═╡ 4d12cf7f-54d7-4a06-9749-a97246e7564e
md"Any other cell which neither has the `#hide` nor the `#export` anchor, is exported as an example block. For such a block the cell is executed and the output is captured.

The output and the code of an example cell appears in the final document."

# ╔═╡ c40eeb9c-dbda-4dc0-9fc3-e91be447ac25
md"### #noop"

# ╔═╡ 929d68db-fac5-4507-8316-a0dcf71df403
md"Nbdev captures the output for cells which are not marked with `#export` by default. but there might be times when you have an example for which you just want to show the code and not the output. Such cases can be handled by the `#noop` anchor. Just start such cells with `#noop` and nbdev will not capture output for such cells."

# ╔═╡ fac3e16d-76ff-476e-a2c4-ff34f7c3f1cc
md"""
```
begin
#noop
"
some experiment
"
function myexperiment
    1+1
end
myexperiment
end
```"""

# ╔═╡ 9fcf1320-0ddf-471c-a584-092586122347
md"In the above case even though Pluto displays the docstring as an output, that output won't be displayed in your document.

As another example let's consider something like the code below"

# ╔═╡ a3ea3f66-a66c-41d1-8c3c-c1cd3a879871
md"
```
#noop
add = 1+1
add
````
"

# ╔═╡ 77b75b82-c684-4377-882a-ff0a4283c7f5
md"The output for such code cell won't be captured in the documentation. Only the code will be displayed once you build the document."

# ╔═╡ a7714f66-3eff-481a-af7f-ae26c2b8e8e2
md"## Docstrings"

# ╔═╡ 926b30a2-1a86-4952-a991-d75b5ce6e545
md"Define a docstring for functions, structs etc. and use the special function `showDoc` to document that function or struct"

# ╔═╡ 5bd889ce-9522-462b-a843-508061926306
"""
> This is a function docstring
"""
function somefunc()
	1+1
end

# ╔═╡ 8c9f1b35-8f88-44ba-8a37-57b92b42980b
md"`showdoc(somefunc)` would show the function like this 👇"

# ╔═╡ 52e3492a-c27d-476f-9010-7d009fd0eeaf
md"`showdoc(somefunc)`"

# ╔═╡ 306239c9-7fd6-442e-8b97-db674398b307
md"## Exporting code"

# ╔═╡ a8955cbd-2530-405f-b876-18bb35e0b008
md"""Once you are done with your notebook and you can run `notebook2script()`. This reads the current notebook in the "nbs" folder, scrubs out all the cells which are do not start with `#export` and creates a module file in the src folder.

The module file in the src folder will have the same name as the current file with the first letter in uppercase.

The exported code file automatically has all the code inside the module with the same name as the file.

For example, the exported source code file would have the following structure 👇

```
module Example

function somefunc(x)
    x+1
end

end
```"""

# ╔═╡ 84f7040c-0ca9-4d61-abfe-789fc446ebb2
md"""
Once you are happy with the code you can export the source code by calling

`notebook2script(nbsdir, srcdir)`

Here "nbs_dir" is the directory where your notebooks are stored and the "src_dir" is the directory where you want your source code files to be exported. This is usually the `src` directory."""

# ╔═╡ 7b1fdf81-870a-4977-8a9d-c867ee2861ff
md"## Tests"

# ╔═╡ 6ace0489-21ae-488f-8ea5-d3087651394c
md"""Nbdev comes with support for [PlutoTests(in alpha)](https://github.com/JuliaPluto/PlutoTest.jl) out of the box.

What this means is that you can write reactive tests within the same notebook where you are writing your code. when you run that particular notebook, the affected tests are executed. So, you can write tests as you write your code.

In addition to reactivity you can also take advantage of other goodies in PlutoTests like "time travel" to see the different stages of your code (follow the link to PlutoTests for further details)."""

# ╔═╡ 50ef31b0-b5f2-4a8e-b123-1681645a031f
md"!!! note
Currently PlutoTests is in **alpha** so use it with caution. However, you can even use the native test suite in Julia to write your tests but won't get features like time travel while using the native test suite.

Automatic execution of tests using github actions or Travis CI is planned for future versions of PlutoTest"

# ╔═╡ bc79b045-8969-4014-b7ac-01282d969c8c
md"## Building documentation"

# ╔═╡ eb26aca4-4b12-435d-941f-edb18bda0c3c
md"""
Do not forget to edit the "index" notebook available in the "nbs" folder. This notebook would be used as your homepage and the contents of this notebook are also used in the README file."""

# ╔═╡ 547f7e70-71bf-44a6-ba62-0e35edf688fc
md"!!! note

*You will get the index notebook in the nbs directory when you create a new project using the nbdev template.*
"

# ╔═╡ b1e0024c-aa32-46a3-baea-534b1187fb7a
md"""To build documentation, goto the project root and launch the Julia prompt and then import Nbdev with `import Nbdev` and then call the following function `Nbdev.build(<your notebook directory)`. Your code at the Julia prompt should look like this

```
import Nbdev
Nbdev.build("nbs")
```

This will run all your notebooks in the notebook directory (/nbs usually) and pick up the cells which are not marked as `#export` or #hide and then create markdown files in the `/docs` directory."""

# ╔═╡ fa7161ef-219f-47b5-b4a0-b18e0d3b5b54
md"## Creating the Table of contents"

# ╔═╡ 54817dac-18b6-40cc-9ef6-b51370164839
md"""The default configuration in the "mkdocs.yml" file (available to you when you create your repository from the template) contains all the necessary default configuration required to create you project page.

In order to test how your website looks like, go to the project root and type `mkdocs serve` . This would serve the project documentation to view locally (localhost).

The table of contents in the default case would be created by mkdocs on the basis of the documents available under the "docs" directory.

If you want to give your own spin to the table of contents then add a config similar to the following in the mkdocs.yml file-->

```
nav:
    - Home: index.md
    - Docs:
      - Nav1: nav1.md
      - Nav2: nav2.md
```
"""

# ╔═╡ 2b4992b3-f1a4-4587-b0d6-cb30efb1dfa1
md"""
!!! note
*Make sure that the "Home" navigation should always be the index.md file*"""

# ╔═╡ 394597e3-42f8-4bd3-8f10-206b5ea5f225
md"For more details about setting up navigation refer [this](https://www.mkdocs.org/user-guide/writing-your-docs/) and [this](https://squidfunk.github.io/mkdocs-material/setup/extensions/python-markdown/?h=table+con+of#table-of-contents) "

# ╔═╡ 094e214e-6681-4de8-b9a4-e73c4b23ff1c
md"## Other required edits"

# ╔═╡ 23ac38d5-c44a-4b8b-86c4-0bc3833add14
md"""In the mkdocs.yml file, edit the "site_name" config to the name that you want to give to your site. The "repo_url" and "repo_name" should be your git repository url and the name of your git repository respectively."""

# ╔═╡ 33eaee48-fc83-4c71-949d-1c9b9a2b32aa
md"## Optional edits"

# ╔═╡ 952bc125-e3b6-474c-90c4-786fdf300dcc
md"There are loads f other optional edits which you can make to the config file in order to customize your project website. Refer the documentation [here](https://squidfunk.github.io/mkdocs-material/getting-started/) to know about the other options."

# ╔═╡ 7429a0be-4f40-4d36-89cd-df56377c4502
md"## Publishing your project website"

# ╔═╡ 8fa0648a-599b-4ed5-870f-6b3d74e25fcd
md"""The first step is to checkin all your project contents to your repoitory. As soon as you push your changes to the remote branch, the github action (placed under "/.github/workflows" creates a new "gh-pages" branch under your project repository, builds the project webpages and pushes the created pages to the ""gh-pages" branch.

If you go to the "Actions" tab in your github repository then you can view the Action under the "All workflows" section. Here, a green tick mark before the action name suggests that all went fine.

Now, go to the Settings tab in your repository, click on Pages menu item from the left menu and in the "Source" drop-down select "gh-pages" as the source and root as the source folder, Click on save.

A detailed step by step description about configuring github pages can be read [here](https://docs.github.com/en/pages/getting-started-with-github-pages/configuring-a-publishing-source-for-your-github-pages-site).

You project website would be published at the following url --> "<username>.github.io/<project reporsitory name>"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.38"
"""

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

[[AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.3"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "670e559e5c8e191ded66fa9ea89c97f10376bb4c"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.38"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═5a14baba-4b8a-4917-b354-246775bb8c1d
# ╠═00f8bf47-47d2-4a99-86ba-d7cd0b7f268b
# ╠═fbb6eadc-80a1-4cd9-ad27-bebcdd6e5981
# ╠═b0c1fe8f-e796-4835-95f8-17990c7421a7
# ╠═e0694cd4-7d07-494d-a010-f705d1cbeb2c
# ╠═15a5c4fd-301e-4f7c-9ad3-28ae2d214224
# ╠═305c50dd-78bd-43d9-8f4e-7074bb046a59
# ╠═2ac2c281-ad4a-4364-8f62-5c8202d4e9c1
# ╠═9258500a-0a33-41cd-b32a-e09e820f84ee
# ╠═dedef829-e43c-41c9-8238-b9d631dccb35
# ╠═99d43364-5ac6-43bb-a70c-20a79a463fdf
# ╠═fb11c61e-2b26-4b2a-9a34-835163ac7f87
# ╠═c9f87989-522c-4adc-9a50-af95e3fa9cfa
# ╠═0d4fb029-bfc8-4208-a6f4-bd5ea981ff82
# ╠═ac65836e-6f4f-4212-b7f4-a648de9c761b
# ╠═789ebdb8-0017-4662-9120-a780584f8368
# ╠═809102dd-b91a-458f-b5c7-2a49c7536d3b
# ╠═f8391a08-edf6-499c-8e8d-8f731dad8fd4
# ╠═4d12cf7f-54d7-4a06-9749-a97246e7564e
# ╠═c40eeb9c-dbda-4dc0-9fc3-e91be447ac25
# ╠═929d68db-fac5-4507-8316-a0dcf71df403
# ╠═fac3e16d-76ff-476e-a2c4-ff34f7c3f1cc
# ╠═9fcf1320-0ddf-471c-a584-092586122347
# ╠═a3ea3f66-a66c-41d1-8c3c-c1cd3a879871
# ╠═77b75b82-c684-4377-882a-ff0a4283c7f5
# ╠═a7714f66-3eff-481a-af7f-ae26c2b8e8e2
# ╠═926b30a2-1a86-4952-a991-d75b5ce6e545
# ╠═5bd889ce-9522-462b-a843-508061926306
# ╠═8c9f1b35-8f88-44ba-8a37-57b92b42980b
# ╠═52e3492a-c27d-476f-9010-7d009fd0eeaf
# ╠═306239c9-7fd6-442e-8b97-db674398b307
# ╠═a8955cbd-2530-405f-b876-18bb35e0b008
# ╠═84f7040c-0ca9-4d61-abfe-789fc446ebb2
# ╠═7b1fdf81-870a-4977-8a9d-c867ee2861ff
# ╠═6ace0489-21ae-488f-8ea5-d3087651394c
# ╠═50ef31b0-b5f2-4a8e-b123-1681645a031f
# ╠═bc79b045-8969-4014-b7ac-01282d969c8c
# ╠═eb26aca4-4b12-435d-941f-edb18bda0c3c
# ╠═547f7e70-71bf-44a6-ba62-0e35edf688fc
# ╠═b1e0024c-aa32-46a3-baea-534b1187fb7a
# ╠═fa7161ef-219f-47b5-b4a0-b18e0d3b5b54
# ╠═54817dac-18b6-40cc-9ef6-b51370164839
# ╠═2b4992b3-f1a4-4587-b0d6-cb30efb1dfa1
# ╠═394597e3-42f8-4bd3-8f10-206b5ea5f225
# ╠═094e214e-6681-4de8-b9a4-e73c4b23ff1c
# ╠═23ac38d5-c44a-4b8b-86c4-0bc3833add14
# ╠═33eaee48-fc83-4c71-949d-1c9b9a2b32aa
# ╠═952bc125-e3b6-474c-90c4-786fdf300dcc
# ╠═7429a0be-4f40-4d36-89cd-df56377c4502
# ╠═8fa0648a-599b-4ed5-870f-6b3d74e25fcd
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
#!/usr/bin/env julia

# Lets add workers programatically, and see if they distribute to different cores.

p = nworkers()
println("Workers allocated ", p)
println("Total workers ", workers())

println("Total processes ", nprocs())

@everywhere function print_cpus()
    run(`ps -eo pid,psr,pcpu,cmd` |> `grep julia`)
end

println("Before map")
print_cpus()

@everywhere function pcpus(mat)
    println("In map")
    print_cpus()
    svd(mat)
end

M = {rand(1000,1000) for i=1:p}

pmap(pcpus, M)

println("After map.")
print_cpus()


text:
using BusinessDays, Dates

d0 = Date(2015, 06, 29) ; d1 = Date(2100, 12, 20)

cal = BusinessDays.Brazil()
@time BusinessDays.initcache(cal)
bdays(cal, d0, d1) # force JIT compilation
@time bdays(cal, d0, d1)
@time for i in 1:1000000 bdays(cal, d0, d1) end


text:
@testset "Character Tables" begin
    for i in allofem
        n, cl, ct = Symmetry.pg_to_chartab(i)
        @test n == eval(Meta.parse(i*"cn"))
        @test isapprox(ct,eval(Meta.parse(i*"ct")))
    end
end


text:
#############################################################################
#############################################################################
#
# A script that runs all unit tests in the project.
#
#############################################################################
#############################################################################

include("../poly_factorization_project.jl")

####
# Execute unit tests for integers
###
include("mod_integers_test.jl")
mod_test_euclid_ints()
mod_test_ext_euclid_ints()

####
# Execute unit tests for polynomials
####
include("mod_polynomials_test.jl")
for i in [3,11,17,19,67,101]
    println("")
    println("testing p = $i")
    mod_sum_test_poly(i)
    mod_prod_test_poly(i)
    mod_prod_derivative_test_poly(i)
    mod_ext_euclid_test_poly(i)
    mod_division_test_poly(i)

end

####
# Execute unit tests for polynomial factorization
####
include("mod_factorization_test.jl")
mod_factor_test_poly()

max_stars_repo_path: test/mod_runtests .jl

max_issues_repo_path: test/mod_runtests .jl

max_forks_repo_path: test/mod_runtests .jl


text:
"""
$(SIGNATURES)

Calculate reference viscosity with variable hard sphere (VHS) model
"""
ref_vhs_vis(Kn, alpha, omega) =
    5.0 * (alpha + 1.0) * (alpha + 2.0) * √π /
    (4.0 * alpha * (5.0 - 2.0 * omega) * (7.0 - 2.0 * omega)) * Kn

"""
$(SIGNATURES)

Calculate collision time with variable hard sphere (VHS) model
"""
vhs_collision_time(ρ, λ, μᵣ, ω) = μᵣ * 2.0 * λ^(1.0 - ω) / ρ

"""
$(SIGNATURES)

For Rykov model prim[end] should be λₜ
"""
vhs_collision_time(prim::AV, muRef, omega) = muRef * 2.0 * prim[end]^(1.0 - omega) / prim[1] # for rykov model prim[end] should be λₜ

"""
$(SIGNATURES)

Calculate mixture collision time from AAP model
"""
function aap_hs_collision_time(prim::AM, mi, ni, me, ne, kn)
    ν = similar(prim, 2)

    ν[1] =
        prim[1, 1] / (mi * (ni + ne)) * 4.0 * sqrt(π) / 3.0 *
        sqrt(1.0 / prim[end, 1] + 1.0 / prim[end, 1]) / (sqrt(2.0) * π * kn) +
        prim[1, 2] / (me * (ni + ne)) * 4.0 * sqrt(π) / 3.0 *
        sqrt(1.0 / prim[end, 1] + 1.0 / prim[end, 2]) / (sqrt(2.0) * π * kn)
    ν[2] =
        prim[1, 1] / (mi * (ni + ne)) * 4.0 * sqrt(π) / 3.0 *
        sqrt(1.0 / prim[end, 1] + 1.0 / prim[end, 2]) / (sqrt(2.0) * π * kn) +
        prim[1, 2] / (me * (ni + ne)) * 4.0 * sqrt(π) / 3.0 *
        sqrt(1.0 / prim[end, 2] + 1.0 / prim[end, 2]) / (sqrt(2.0) * π * kn)

    return 1.0 ./ ν
end

"""
$(SIGNATURES)

Source term of AAP model in DifferentialEquations.jl
"""
function aap_hs_diffeq!(du, u, p, t)
    if length(u) == 6
        I₁, I₂, I₃, E₁, E₂, E₃ = u
        w = [
            I₁ E₁
            I₂ E₂
            I₃ E₃
        ]
    elseif length(u) == 8
        I₁, I₂, I₃, I₄, E₁, E₂, E₃, E₄ = u
        w = [
            I₁ E₁
            I₂ E₂
            I₃ E₃
            I₄ E₄
        ]
    elseif length(u) == 10
        I₁, I₂, I₃, I₄, I₅, E₁, E₂, E₃, E₄, E₅ = u
        w = [
            I₁ E₁
            I₂ E₂
            I₃ E₃
            I₄ E₄
            I₅ E₅
        ]
    else
    end

    τᵢ, τₑ, mi, ni, me, ne, kn, γ = p
    τ = [τᵢ, τₑ]

    # modified variables
    prim = mixture_conserve_prim(w, γ)
    mixprim = aap_hs_prim(prim, τ, mi, ni, me, ne, kn)
    mixw = mixture_conserve_prim(mixprim, γ)

    if length(u) == 6
        du[1] = (mixw[1, 1] - I₁) / τᵢ
        du[2] = (mixw[2, 1] - I₂) / τᵢ
        du[3] = (mixw[3, 1] - I₃) / τᵢ
        du[4] = (mixw[1, 2] - E₁) / τₑ
        du[5] = (mixw[2, 2] - E₂) / τₑ
        du[6] = (mixw[3, 2] - E₃) / τₑ
    elseif length(u) == 8
        du[1] = (mixw[1, 1] - I₁) / τᵢ
        du[2] = (mixw[2, 1] - I₂) / τᵢ
        du[3] = (mixw[3, 1] - I₃) / τᵢ
        du[4] = (mixw[4, 1] - I₄) / τᵢ
        du[5] = (mixw[1, 2] - E₁) / τₑ
        du[6] = (mixw[2, 2] - E₂) / τₑ
        du[7] = (mixw[3, 2] - E₃) / τₑ
        du[8] = (mixw[4, 2] - E₄) / τₑ
    elseif length(u) == 10
        du[1] = (mixw[1, 1] - I₁) / τᵢ
        du[2] = (mixw[2, 1] - I₂) / τᵢ
        du[3] = (mixw[3, 1] - I₃) / τᵢ
        du[4] = (mixw[4, 1] - I₄) / τᵢ
        du[5] = (mixw[5, 1] - I₅) / τᵢ
        du[6] = (mixw[1, 2] - E₁) / τₑ
        du[7] = (mixw[2, 2] - E₂) / τₑ
        du[8] = (mixw[3, 2] - E₃) / τₑ
        du[9] = (mixw[4, 2] - E₄) / τₑ
        du[10] = (mixw[5, 2] - E₅) / τₑ
    end

    return nothing
end

"""
$(SIGNATURES)

RHS-ODE of Boltzmann equation
"""
function boltzmann_ode!(df, f::AA{T,3}, p, t) where {T<:Real}
    Kn, M, phi, psi, phipsi = p
    df .= boltzmann_fft(f, Kn, M, phi, psi, phipsi)
end

"""
$(SIGNATURES)

RHS-ODE of Boltzmann equation with non-uniform velocity
"""
function boltzmann_nuode!(df, f::AA{T,3}, p, t) where {T<:Real}
    Kn, M, phi, psi, phipsi, u, v, w, vnu, u1, v1, w1, vuni = p

    nu = length(u)
    nv = length(v)
    nw = length(w)
    nu1 = length(u1)
    nv1 = length(v1)
    nw1 = length(w1)

    curve = itp.RegularGridInterpolator((u, v, w), f)
    _f = reshape(curve(vuni), nu1, nv1, nw1)
    _df = boltzmann_fft(_f, Kn, M, phi, psi, phipsi)

    curve1 = itp.RegularGridInterpolator((u1, v1, w1), _df)
    df .= reshape(curve1(vnu), nu, nv, nw)
end

"""
$(SIGNATURES)

RHS-ODE of BGK equation
"""
function bgk_ode!(df, f::AA{T}, p, t) where {T<:Real}
    g, τ = p
    df .= (g .- f) ./ τ
end


text:
module Perturbo
export ElecHam, bands, hamiltonian!

include("./types.jl")
include("./WignerSeitzCell.jl")
include("./load_data.jl")
include("./electronic_structure.jl")

end # module


text:
module IWSLT
using Fetch
using DataDeps

using ..Datasets: Dataset
import ..Datasets: testfile, devfile, trainfile
import ..Datasets: token_freq, get_vocab

export IWSLT2016

function __init__()
    iwslt2016_init()
end

include("./iwslt2016.jl")

end