text:
#-------------------------------------------------------------------------------
# Set of functions useful for the execution of the main algorithm.
#-------------------------------------------------------------------------------

"""

    predefined_sample_radius(a::Vector{Float64}, b::Vector{Float64})

Calculates the initial value of the radius of the sample region 'δ' based on the lower and upper limits 'a' and 'b'.

    - 'a': n-dimensional vector with the lower bounds.

    - 'b': n-dimensional vector with the upper bounds.

Returns the initial value of 'δ'.

"""
function predefined_sample_radius(
                                    a::Vector{Float64},
                                    b::Vector{Float64}
                                    )

    n = min( length(a), length(b) )
    aux = b[1] - a[1]
    aux2 = 0.0

    for i = 2:n

        aux2 = b[i] - a[i]

        if aux2 < aux

            aux = aux2

        end

    end

    return min( aux / 2.0, 1.0 )

end

"""

    verify_initial_room(n::Int64, δ::Float64, a::Vector{Float64},
                        b::Vector{Float64})

Checks whether the bounds satisfy the conditions 'b' >= 'a' + 2 * 'δ'.

    - 'n': dimension of the search space.

    - 'δ': radius of the sample set.

    - 'a': n-dimensional vector with the lower bounds.

    - 'b': n-dimensional vector with the upper bounds.

Returns a boolean value.

"""
function verify_initial_room(
                                n::Int64,
                                δ::Float64,
                                a::Vector{Float64},
                                b::Vector{Float64}
                                )

    for i = 1:n

        if ( b[i] - a[i] ) < ( 2.0 * δ )

            return false

        end

    end

    return true

end

"""

    correct_guess_bounds!(n::Int64, δ::Float64, a::Vector{Float64},
                            b::Vector{Float64}, x::Vector{Float64},
                            ao::Vector{Float64}, bo::Vector{Float64})

Fixes initial guess 'x' and sets the lower and upper bound vectors 'ao' and 'bo'
shifted to the origin 'x'.

    - 'n': dimension of the search space.

    - 'δ': radius of the sample set.

    - 'a': n-dimensional vector with the lower bounds.

    - 'b': n-dimensional vector with the upper bounds.

The function modifies the argument:

    - 'x': n-dimensional vector (initial guess).

    - 'ao': n-dimensional vector with the shifted lower bounds.

    - 'bo': n-dimensional vector with the shifted upper bounds.

"""
function correct_guess_bounds!(
                                n::Int64,
                                δ::Float64,
                                a::Vector{Float64},
                                b::Vector{Float64},
                                x::Vector{Float64},
                                ao::Vector{Float64},
                                bo::Vector{Float64}
                                )

    for i=1:n

        ao[i] = a[i] - x[i]
        bo[i] = b[i] - x[i]

        if ao[i] >= - δ

            if ao[i] >= 0.0

                x[i] = a[i]
                ao[i] = 0.0
                bo[i] = b[i] - a[i]

            else

                x[i] = a[i] + δ
                ao[i] = - δ
                bo[i] = max( b[i] - x[i], δ )

            end

        elseif bo[i] <= δ

            if bo[i] <= 0.0

                x[i] = b[i]
                ao[i] = a[i] - b[i]
                bo[i] = 0.0

            else

                x[i] = b[i] - δ
                ao[i] = min( a[i] - x[i], - δ )
                bo[i] = δ

            end

        end

    end

end

"""

    relative_reduction(model::AbstractModel, func_list::Array{Function, 1},
                        pred_red::Float64, y::Vector{Float64})

Computes the simplified relative reduction (using f_{i} instead of fmin) in the point 'y'.

    - 'model': model of LinearModel or QuadraticModel type.

    - 'func_list': list containing the functions that determine the objective
    function fmin.

    - 'pred_red': precalculated predicted reduction.

    - 'y': n-dimensional vector (point of interest).

Returns the value of the simplified relative reduction, the real reduction, the value of f_i(y) and the index 'model.imin'.

"""
function relative_reduction(
                            model::AbstractModel,
                            func_list::Array{Function, 1},
                            pred_red::Float64,
                            y::Vector{Float64}
                            )

    f_y = fi_eval(func_list, model.imin[], y)
    real_red = model.fval[model.kopt[] + 1] - f_y
    ρ = real_red / pred_red

    return ρ, real_red, f_y, model.imin[]

end

"""

    relative_reduction(model::AbstractModel, func_list::Array{Function, 1}, r::Int64,
                        pred_red::Float64, y::Vector{Float64}, imin_set::Vector{Bool})

Computes the relative reduction obtained in the point 'y'.

    - 'model': model of LinearModel or QuadraticModel type.

    - 'func_list': list containing the functions that determine the objective
    function fmin.

    - 'r': number of functions that make up the objective function fmin.

    - 'pred_red': precalculated predicted reduction.

    - 'y': n-dimensional vector (point of interest).

The function modifies the argument:

    - 'imin_set': boolean vector with the indexes belonging to the I_{min}(y) set.

Returns the value of the relative reduction, the real reduction, the value of fmin(y) and an index 'idx_y' belonging to Imin(y).

"""
function relative_reduction!(
                            model::AbstractModel,
                            func_list::Array{Function, 1},
                            r::Int64,
                            pred_red::Float64,
                            y::Vector{Float64},
                            imin_set::Vector{Bool}
                            )

    fmin_y, idx_y = fmin_eval!(func_list, r, y, imin_set)
    real_red = model.fval[model.kopt[] + 1] - fmin_y
    ρ = real_red / pred_red

    return ρ, real_red, fmin_y, idx_y

end

"""

    reconstruct_original_point!(model::AbstractModel, idx::Int64, a::Vector{Float64}, b::Vector{Float64},
                                ao::Vector{Float64}, bo::Vector{Float64}, x::Vector{Float64})

Constructs the model based in the given information.

    - 'model': model of LinearModel or QuadraticModel type.

    - 'idx': position of the point in 'Y' set.

    - 'a': n-dimensional vector with the lower bounds.

    - 'b': n-dimensional vector with the upper bounds.

    - 'ao': n-dimensional vector with the shifted lower bounds.

    - 'bo': n-dimensional vector with the shifted upper bounds.

The function modifies the argument:

    - 'x': n-dimensional vector (point of interest).

"""
function reconstruct_original_point!(
                                    model::AbstractModel,
                                    idx::Int64,
                                    a::Vector{Float64},
                                    b::Vector{Float64},
                                    ao::Vector{Float64},
                                    bo::Vector{Float64},
                                    x::Vector{Float64}
                                    )

    for i = 1:model.n

        if model.Y[idx, i] == ao[i]

            x[i] = a[i]

        elseif model.Y[idx, i] == bo[i]

            x[i] = b[i]

        else

            x[i] = min( max( a[i], model.xbase[i] + model.Y[idx, i] ), b[i] )

        end

    end

end

"""

    print_warning(flag::Int64)

Prints information about the exit flag.

    - 'flag': exit flag.

"""
function print_warning(
                        exit_flag::Symbol
                        )

    if exit_flag == :success

        printstyled("Success: ", bold=true, color=:light_green)
        println("The algorithm succeeded.")

    elseif exit_flag == :max_iterations

        printstyled("Warning: ", bold=true, color=:light_red)
        println("The maximum number of iterations has been reached.")

    elseif exit_flag == :max_evaluations

        printstyled("Warning: ", bold=true, color=:light_red)
        println("The maximum number of function evaluations has been reached.")

    elseif exit_flag == :nondescent

        printstyled("Warning: ", bold=true, color=:light_red)
        println("The calculated direction is not descent for the model.")

    elseif exit_flag == :better_point

        printstyled("Notice: ", bold=true, color=:light_yellow)
        println("A better value will be obtained when evaluating the function on 'xopt'.")

    else

        printstyled("Warning: ", bold=true, color=:light_red)
        println("Exit flag not specified.")

    end
    println("--------------------------------------------------------------------------------")

end

"""

    print_iteration(it_flag::Symbol, status_flag::Symbol, full_calc::Bool, nit::Int64,
                        nf::Int64, imin_idx::Int64, δ::Float64, Δ::Float64, π::Float64,
                        ρ::Float64, pred_red::Float64, real_red::Float64, fopt::Float64,
                        xopt::Vector{Float64}, d::Vector{Float64})

Prints information about the iterations.

    - 'it_flag': iteration flag.

    - 'status_flag': status flag.

    - 'full_calc': indicates if ρ was full calculated.

    - 'nit': number of iterations.

    - 'nf': number of function evaluations.

    - 'imin_idx': index corresponding to the available function value.

    - 'δ': sample set radius.

    - 'Δ': trust-region radius.

    - 'π': stationarity measure.

    - 'ρ': relative reduction.

    - 'pred_red': predicted reduction.

    - 'real_red': real reduction.

    - 'fopt': function value in 'xopt'

    - 'xopt': best point so far.

    - 'd': last computed direction.

"""
function print_iteration(
                            it_flag::Symbol,
                            status_flag::Symbol,
                            full_calc::Bool,
                            nit::Int64,
                            nf::Int64,
                            imin_idx::Int64,
                            δ::Float64,
                            Δ::Float64,
                            π::Float64,
                            ρ::Float64,
                            pred_red::Float64,
                            real_red::Float64,
                            fopt::Float64,
                            xopt::Vector{Float64},
                            d::Vector{Float64}
                            )

    if nit == 0

        println("--------------------------------------------------------------------------------")
        println("--------------------------------------------------------------------------------")

    end

    if ( nit == 0 && it_flag == :max_evaluations )

        println("Iteration     : $( nit )")
        println("Func. eval.   : $( nf )")
        println("Radius δ      : $( δ )")
        println("Radius Δ      : $( Δ )")
        println("I_min index   : $( imin_idx )")
        println("Best point    : $( xopt )")
        println("Func. val.    : $( fopt )")
        println("Iter. type    : $( it_flag )")
        println("--------------------------------------------------------------------------------")

    elseif it_flag == :nonspecified

        println("Iteration     : $( nit )")
        println("Func. eval.   : $( nf )")
        println("Radius δ      : $( δ )")
        println("Radius Δ      : $( Δ )")
        println("Stationarity π: $( π )")
        println("I_min index   : $( imin_idx )")
        println("Best point    : $( xopt )")
        println("Func. val.    : $( fopt )")
        println("Iter. type    : $( it_flag )")
        println("--------------------------------------------------------------------------------")

    elseif it_flag == :criticality

        println("Iteration     : $( nit )")
        println("Func. eval.   : $( nf )")
        println("Radius δ      : $( δ )")
        println("Radius Δ      : $( Δ )")
        println("Stationarity π: $( π )")
        println("I_min index   : $( imin_idx )")
        println("Best point    : $( xopt )")
        println("Func. val.    : $( fopt )")
        println("Iter. type    : $( it_flag )")
        println("Step status   : $( status_flag )")
        println("Direction d   : $( d )")
        println("--------------------------------------------------------------------------------")

    else

        println("Iteration     : $( nit )")
        println("Func. eval.   : $( nf )")
        println("Radius δ      : $( δ )")
        println("Radius Δ      : $( Δ )")
        println("Stationarity π: $( π )")
        println("I_min index   : $( imin_idx )")
        println("Best point    : $( xopt )")
        println("Func. val.    : $( fopt )")
        println("Iter. type    : $( it_flag )")
        println("Step status   : $( status_flag )")
        println("Direction d   : $( d )")
        println("Full ρ        : $( full_calc )")
        println("Rel. reduc. ρ : $( ρ )")
        println("Pred. reduc.  : $( pred_red )")
        println("Real reduc.   : $( real_red )")
        println("--------------------------------------------------------------------------------")

    end

end

"""

    print_info(model::AbstractModel, output::LOWDEROutput, exit_flag::Symbol, it_flag::Symbol,
                status_flag::Symbol, verbose::Int64, nit::Int64, nf::Int64, δ::Float64, Δ::Float64,
                π::Float64, full_calc::Bool, pred_red::Float64, real_red::Float64, ρ::Float64,
                d::Vector{Float64})

Prints information about the last iteration of LOWDER.

    - 'model': model of LinearModel or QuadraticModel type.

    - 'output': object of LOWDEROutput type.

    - 'it_flag': iteration flag.

    - 'status_flag': status flag.

    - 'full_calc': indicates if fmin was calculated at 'xopt'.

    - 'verbose': level of vebosity.

    - 'δ': sample set radius.

    - 'Δ': trust-region radius.

    - 'π': stationarity measure.

    - 'ρ': relative reduction.

    - 'pred_red': predicted reduction.

    - 'real_red': real reduction.

    - 'd': last computed direction.

"""
function print_info(
                    model::AbstractModel,
                    output::LOWDEROutput,
                    it_flag::Symbol,
                    status_flag::Symbol,
                    full_calc::Bool,
                    verbose::Int64,
                    δ::Float64,
                    Δ::Float64,
                    π::Float64,
                    ρ::Float64,
                    pred_red::Float64,
                    real_red::Float64,
                    d::Vector{Float64}
                    )

    if verbose != 0

        print_iteration(it_flag, status_flag, full_calc, output.iter, output.nf, output.index, δ, Δ, π, ρ, pred_red, real_red, output.f, output.solution, d)

        if verbose ≥ 2

            print_warning(output.status)

            if !( model.kopt[] == 0 || full_calc )

                print_warning(:better_point)

            end

            if verbose == 3

                show_output(output)

            end

        end

    end

end

"""

    projection_active_set!(v::Vector{Float64}, active_set::Vector{Bool},
                            proj_v::Vector{Float64}; sym::Bool=false)

Computes the projection of the vector 'v' on the set of active constraints

    - 'v': n-dimensional vector (point of interest).

    - 'active_set': boolean n-dimensional vector with the active constraints.

The function modifies the argument:

    - 'proj_v': n-dimensional vector (projection of 'v').

The optional argument is:

    - 'sym': indicates if it is calculated based on the symmetric direction of 'v'. Set to 'false' by default.

"""
function projection_active_set!(
                                v::Vector{Float64},
                                active_set::Vector{Bool},
                                proj_v::Vector{Float64};
                                sym::Bool=false
                                )

    for i=1:length(v)

        if active_set[i]

            proj_v[i] = 0.0

        else

            if sym

                proj_v[i] = - v[i]

            else

                proj_v[i] = v[i]

            end

        end

    end

end

"""

    step_projection!(model::AbstractModel, a::Vector{Float64},
                        b::Vector{Float64}, x::Vector{Float64},
                        d::Vector{Float64})

Projects 'x' := 'xopt' + 'd' over the box constraints and
modifies 'd' to conform to the new point 'x'.

    - 'model': model of LinearModel or QuadraticModel type.

    - 'a': n-dimensional vector with the lower bounds.

    - 'b': n-dimensional vector with the upper bounds.

The function modifies the argument:

    - 'x': n-dimensional vector (point of interest).

    - 'd': n-dimensional vector (direction).

"""
function step_projection!(
                            model::AbstractModel,
                            a::Vector{Float64},
                            b::Vector{Float64},
                            x::Vector{Float64},
                            d::Vector{Float64}
                            )

    for i = 1:model.n

        x[i] = max( a[i], min( model.xopt[i] + d[i], b[i] ) )
        d[i] = x[i] - model.xopt[i]

    end

end

"""

    compute_alpha_linear(Δ::Float64, a::Vector{Float64}, b::Vector{Float64},
                            x::Vector{Float64}, d::Vector{Float64}, s::Vector{Float64})

Calculates the step length α for 'trsbox!' routine if the model is linear.

    - 'Δ': trust-region radius.

    - 'a': n-dimensional vector with the lower bounds.

    - 'b': n-dimensional vector with the upper bounds.

    - 'x': n-dimensional vector (current point).

    - 'd': n-dimensional vector (first direction).

    - 's': n-dimensional vector (second direction).

Returns the value of α, a boolean indicating whether α_Δ is chosen,
and an index indicating which box constraint is met in the case where α_B is chosen.

"""
function compute_alpha_linear(
                                Δ::Float64,
                                a::Vector{Float64},
                                b::Vector{Float64},
                                x::Vector{Float64},
                                d::Vector{Float64},
                                s::Vector{Float64}
                                )

    roots = []
    α_j = Inf
    α_B = Inf
    idx = 0

    # Computes α_B, the largest number such that 'a' ≤ 'x' + 'd' + 'α'*'s' ≤ 'b'.
    for j=1:length(x)

        if s[j] < 0.0

            α_j = ( a[j] - x[j] - d[j] ) / s[j]

        elseif s[j] > 0.0

            α_j = ( b[j] - x[j] - d[j] ) / s[j]

        end

        if α_j < α_B

            α_B = α_j
            idx = j

        end

    end

    # Computes α_Δ, the largest number such that ||'d' + 'α'*'s'|| ≤ 'Δ'.

    a_2 = dot(s, s)
    a_1 = 2.0 * dot(d, s)
    a_0 = dot(d, d) - Δ ^ 2.0

    solve_quadratic!(a_2, a_1, a_0, roots)

    α_Δ = maximum(roots)

    if ( α_Δ == NaN ) || ( α_Δ < 0.0 )

        return α_B, false, idx

    else

        α = min(α_B, α_Δ)

        if α == α_Δ

            return α, true, 0

        else

            return α, false, idx

        end

    end

end

"""

    new_search_direction!(pdTpd::Float64, pdTpg::Float64, pgTpg::Float64,
                            proj_d::Vector{Float64}, proj_grad::Vector{Float64})

Calculates a new search direction for the 'trsbox!' routine.

    - 'pdTpd': product of the direction 'd' by itself.

    - 'pdTpg': product of the projection of the direction 'd' by the gradient vector
                of the model evaluated in 'xopt'.

    - 'pgTpg': product of the projection of the gradient vector of the model evaluated in 'xopt' by itself.

    - 'proj_d': projection of the search direction 'd'

The function modifies the argument:

    - 'proj_grad': projection of the gradient of the model in 'xopt'.

"""
function new_search_direction!(
                                pdTpd::Float64,
                                pdTpg::Float64,
                                pgTpg::Float64,
                                proj_d::Vector{Float64},
                                proj_grad::Vector{Float64}
                                )

    roots = []
    α = 0.0
    β = 0.0

    if pdTpg ≈ 0.0

        # If P_{I}(d)'P_{I}(∇model(xk+d)) = 0, then the new direction is a multiple of P_{I}(∇Q(xk+d)).

        solve_quadratic!( pgTpg, 0.0, pdTpd, roots)

        β = minimum(roots)

        if ( β == NaN ) || ( β > 0.0 )

            @. proj_grad = 0.0

        else

            @. proj_grad *= β

        end

    else

        # Otherwise, the new search direction is a linear combination of P_{I}(d) and P_{I}(∇model(xk+d)).

        aux = ( pdTpd ^ 2.0 * pgTpg ) / pdTpg ^ 2.0 - pdTpd

        solve_quadratic!( aux, 0.0, - pdTpd, roots )

        α = maximum(roots)

        if ( α == NaN ) || ( α * pdTpd ≈ 0.0 )

            @. proj_grad = 0.0

        else

            β = - pgTpg / ( α * pdTpd )

            if ( α * pdTpg + β * pgTpg ) < 0.0

                @. proj_grad = α * proj_d + β * proj_grad

            else

                α = minimum(roots)

                if ( α == NaN ) || ( α * pdTpd ≈ 0.0 )

                    @. proj_grad = 0.0

                else

                    β = - pgTpg / ( α * pdTpd )

                    if ( α * pdTpg + β * pgTpg ) < 0.0

                        @. proj_grad = α * proj_d + β * proj_grad

                    else

                        @. proj_grad = 0.0

                    end

                end

            end

        end

    end

end

"""

    solve_quadratic!(a::Float64, b::Float64, c::Float64,
                        roots::Vector{Any})

Calculates the roots of the quadratic equation a * x ^ 2 + b * x + c = 0.

    - 'a': quadratic term.

    - 'b': linear term.

    - 'c': constant.

The function modifies the argument:

    - 'roots': vector with the roots of the quadratic equation.

"""
function solve_quadratic!(
                            a::Float64,
                            b::Float64,
                            c::Float64,
                            roots::Vector{Any}
                            )

    if a ≈ 0.0

        root = - c / b

        if !( isinf(root) )

            append!(roots, root)

        elseif ( b ≈ 0.0 ) && ( c ≈ 0.0 )

            append!(roots, 0.0)

        else

            append!(roots, NaN)

        end

    else

        b = b / a
        c = c / a

        if c ≈ 0.0

            append!(roots, 0.0)
            append!(roots, -b)

        else

            Δ = b ^ 2.0 - 4.0 * c

            if isinf(Δ)

                root_1 = -b
                root_2 = c / root_1

                if isinf(root_2)

                    append!(roots, root_1)

                else

                    append!(roots, root_1)
                    append!(roots, root_2)

                end

            else

                if Δ < 0.0

                    append!(roots, NaN)

                elseif Δ == 0.0

                    append!(roots, - 0.5 * b)

                else

                    if b < 0.0

                        root_1 = 0.5 * ( -b + sqrt(Δ) )

                    else

                        root_1 = 0.5 * ( - b - sqrt(Δ) )

                    end

                    root_2 = c / root_1

                    if isinf(root_2)

                        append!(roots, root_1)

                    else

                        append!(roots, root_1)
                        append!(roots, root_2)

                    end

                end

            end

        end

    end

end

"""

    binary_search(lower_value::Float64, upper_value::Float64,
                    stop_condition::Function, ε::Float64)

Performs a binary search on the range [low_value, high_value] to
satisfy 'stop_condition' with tolerance 'ϵ'.

    - 'lower_value': lower value.

    - 'upper_value': upper value.

    - 'stop_condition': stop condition.

    - 'ϵ': tolerance.

Returns a estimate that satisfy the 'stop_condition'.

"""
function binary_search(
                        lower_value::Float64,
                        upper_value::Float64,
                        stop_condition::Function,
                        ε::Float64
                        )

    if !( stop_condition(lower_value) )

        return NaN

    elseif stop_condition(upper_value)

        return upper_value

    else

        while true

            mean_value = (lower_value + upper_value) / 2.0

            if stop_condition(mean_value)

                if (upper_value - mean_value) ≤ ε

                    return mean_value

                else

                    lower_value = mean_value

                end

            else

                if (mean_value - lower_value) ≤ ε

                    return lower_value

                else

                    upper_value = mean_value

                end

            end

        end

    end

end

"""

    cond_θB(θ::Float64, a::Vector{Float64}, b::Vector{Float64}, xopt::Vector{Float64},
        d::Vector{Float64}, proj_d::Vector{Float64}, s::Vector{Float64})

Verifies in 'Θ' satisfies the ΘB condition for 'trsbox!' routine.

    - 'Θ': value that will be tested.

    - 'a': n-dimensional vector with the lower bounds.

    - 'b': n-dimensional vector with the upper bounds.

    - 'xopt': best sample point so far.

    - 'd': first search direction 'd'.

    - 'proj_d': projection of the direction 'd'.

    - 's': second search direction 's'.

Returns a boolean.

"""
function cond_θB(
                    θ::Float64,
                    a::Vector{Float64},
                    b::Vector{Float64},
                    xopt::Vector{Float64},
                    d::Vector{Float64},
                    proj_d::Vector{Float64},
                    s::Vector{Float64}
                    )

    for i=1:length(a)

        xnewi = xopt[i] + d[i] - proj_d[i] + cos(θ) * proj_d[i] + sin(θ) * s[i]

        if ( xnewi < a[i] ) || ( xnewi > b[i] )

            return false

        end

    end

    return true

end

"""

    cond_θQ_linear(θ::Float64, gTd::Float64, gTpd::Float64, gTs::Float64)

Verifies in 'Θ' satisfies the ΘQ condition for 'trsbox!' routine, in the liner case.

    - 'Θ': value that will be tested.

    - 'gTd': product of the direction 'd' by the gradient vector of the model evaluated in 'xopt'.

    - 'gTd': product of the projection of the direction 'd' by the gradient vector of the
    model evaluated in 'xopt'.

    - 'gTs': product of the direction 's' by the gradient vector of the model evaluated in 'xopt'.

Returns a boolean.

"""
function cond_θQ_linear(
                        θ::Float64,
                        gTd::Float64,
                        gTpd::Float64,
                        gTs::Float64
                        )

    if ( gTd - gTpd + cos(θ) * gTpd + sin(θ) * gTs ) < 0.0

        return true

    else

        return false

    end

end

"""

    update_active_set!(a::Vector{Float64}, b::Vector{Float64},
                        xopt::Vector{Float64}, d::Vector{Float64},
                        active_set::Vector{Bool})

Updates the active set.

    - 'a': n-dimensional vector with the lower bounds.

    - 'b': n-dimensional vector with the upper bounds.

    - 'xopt': best sample point so far.

    - 'd': search direction 'd'.

The function modifies the argument:

    - 'active_set': boolean n-dimensional vector with the active constraints.

"""
function update_active_set!(
                            a::Vector{Float64},
                            b::Vector{Float64},
                            xopt::Vector{Float64},
                            d::Vector{Float64},
                            active_set::Vector{Bool}
                            )

    for i=1:length(a)

        if !( active_set[i] )

            xnewi = xopt[i] + d[i]

            if ( a[i] == xnewi ) || ( xnewi == b[i] )

                active_set[i] == true

            end

        end

    end

end

"""

    save_info!(model::AbstractModel, it_flag::Symbol, status_flag::Symbol,
                full_calc::Bool, nit::Int64, nf::Int64, δ::Float64,
                Δ::Float64, π::Float64, ρ::Float64, pred_red::Float64,
                real_red::Float64, d::Vector{Float64}, file::IO)

Saves information about the iterations.

    - 'model': model of LinearModel or QuadraticModel type.

    - 'it_flag': iteration flag.

    - 'status_flag': status flag.

    - 'full_calc': indicates if ρ was full calculated.

    - 'nit': number of iterations.

    - 'nf': number of function evaluations.

    - 'δ': sample set radius.

    - 'Δ': trust-region radius.

    - 'π': stationarity measure.

    - 'ρ': relative reduction.

    - 'pred_red': predicted reduction.

    - 'real_red': real reduction.

    - 'd': last computed direction.

    - 'file': IO.

"""
function save_info!(
                    model::AbstractModel,
                    it_flag::Symbol,
                    status_flag::Symbol,
                    full_calc::Bool,
                    nit::Int64,
                    nf::Int64,
                    δ::Float64,
                    Δ::Float64,
                    π::Float64,
                    ρ::Float64,
                    pred_red::Float64,
                    real_red::Float64,
                    d::Vector{Float64},
                    file::IO
                    )

    text = "$( nit ) $( nf ) $( model.imin[] ) $( model.kopt[] ) $( δ ) $( Δ ) $( π ) $( it_flag ) $( status_flag ) $( full_calc ) $( ρ ) $( pred_red ) $( real_red ) $( model.c[] )"

    print( file, text )

    for i = 1:model.n

        print( file, " $( model.g[i] )" )

    end

    for i = 1:model.n

        print( file, " $( model.xbase[i] )" )

    end

    for i = 1:model.n

        print( file, " $( model.xopt[i] )" )

    end

    for i = 1:model.m

        print( file, " $( model.fval[i] )" )

    end

    for i = 1:(model.m - 1)

        for j = 1:model.n

            print( file, " $( model.Y[i, j] )" )

        end

    end

    for i = 1:(model.n - 1)

        print( file, " $( d[i] )" )

    end

    println( file, " $( d[model.n] )" )

end

"""

    choose_index_altmov(model::LinearModel)

Chooses the point that must leave the sampling set for 'altmov' steps.

    - 'model': model of LinearModel or QuadraticModel type.

Returns an index 't' which indicates the point that must leave the sample set.

"""
function choose_index_altmov(
                                model::AbstractModel
                            )

    ( val, idx_t ) = findmax( model.dst )

    if idx_t == model.kopt[]

        idx_t = 0

    end

    return idx_t

end

max_stars_repo_name: aschwertner/LOWDER

max_issues_repo_name: aschwertner/LOWDER

max_forks_repo_name: aschwertner/LOWDER


text:
"""
    default_PARF(grid, ΔT, iterations)
Generate default hourly surface PAR.
"""
function default_PARF(grid, ΔT, iterations)
    PAR = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.3871666666666666, 87.10258333333333, 475.78150000000016, 929.2737916666669,
           1232.3633333333337, 1638.918916666667, 1823.7921666666664, 1906.2769583333336, 1776.0280416666667,
           1678.5026249999999, 1410.216666666667, 815.4129583333336, 525.104, 135.993, 2.9493750000000003, 0.0, 0.0, 0.0,]
    total_days = (ΔT * iterations) ÷ 86400 + 1
    PAR_day = zeros(grid.Nx, grid.Ny, 24)
    for i in 1:24
        PAR_day[:,:,i] .= PAR[i]
    end
    PAR_domain = repeat(PAR_day, outer = (1,1,total_days))
    return PAR_domain
end

"""
    default_temperature(grid, ΔT, iterations)
Generate default hourly temperature.
"""
function default_temperature(grid, ΔT, iterations)
    temp = [26.646446609406727, 26.56698729810778, 26.517037086855467, 26.5, 26.517037086855467, 26.56698729810778,
            26.646446609406727, 26.75, 26.87059047744874, 27.0, 27.12940952255126, 27.25, 27.353553390593273,
            27.43301270189222, 27.482962913144533, 27.5, 27.482962913144533, 27.43301270189222, 27.353553390593273,
            27.25, 27.12940952255126, 27.0, 26.87059047744874, 26.75]
    total_days = (ΔT * iterations) ÷ 86400 + 1
    temp_day = zeros(grid.Nx, grid.Ny, grid.Nz, 24)
    for i in 1:24
        temp_day[:,:,end,i] .= temp[i]
    end
    # vertical temperature gradient
    for j in grid.Nz-1:-1:1
        CUDA.@allowscalar temp_day[:,:,j,:] .= temp_day[:,:,j+1,:] .+ 0.04 * grid.zC[j]
    end
    temp_domain = repeat(temp_day, outer = (1,1,1,total_days))
    return temp_domain
end

"""
    bgc_params_default()
Generate default biogeochemical parameter values
"""
function bgc_params_default()
    params = Dict(
        "kw"       => 0.046,              # PAR attenuation (/m)
        "kc"       => 0.04,               # PAR attenuation (m²/mgChl)
        "kDOC"     => 1/30/86400,         # Remineralization rate for DOC, turn over time: a month (per second)
        "Nit"      => 1/30/86400,         # Nitrification rate for NH4
        "kDON"     => 1/30/86400,         # Remineralization rate for DON, turn over time: a month (per second)
        "kDOP"     => 1/30/86400,         # Remineralization rate for DON, turn over time: a month (per second)
        "kPOC"     => 1/30/86400,         # Remineralization rate for POC, turn over time: a month (per second)
        "kPON"     => 1/30/86400,         # Remineralization rate for PON, turn over time: a month (per second)
        "kPOP"     => 1/30/86400,         # Remineralization rate for PON, turn over time: a month (per second)
        "κh"       => 0.0e-6,             # Horizontal diffusion
        "κv"       => 0.0e-6,             # Vertical diffusion
        "κhP"      => 0.0e-6,             # Horizontal diffusion for individuals
        "κvP"      => 0.0e-6,             # Vertical diffusion for individuals
    )
    return params
end

"""
    phyt_params_default(N::Int64, mode::AbstractMode)
Generate default phytoplankton parameter values based on `AbstractMode` and species number `N`.
"""
#=
CH includes cabohydrate and lipids
┌─────────┬────────────────┬────────────────────┬───────────────────────────┬──────────────────────────┐
│         │ Micromonas sp. │ Ostreococcus tauri │ Thalassiosira weissflogii │ Thalassiosira pseudonana │
├─────────┼────────────────┼────────────────────┼───────────────────────────┼──────────────────────────┤
│ PRO:DNA │          34.83 │              22.75 │                     67.99 │                   135.36 │
│ RNA:DNA │           1.73 │               1.71 │                      3.70 │                     6.49 │
│ CHL:DNA │           3.45 │               2.16 │                      9.83 │                    18.05 │
│  CH:DNA │          23.32 │              19.01 │                     79.26 │                   121.48 │
│   k_pro │       7.06e-05 │           1.50e-04 │                  7.64e-05 │                 1.17e-04 │
│   k_dna │       1.62e-07 │           5.44e-07 │                  1.16e-07 │                 9.61e-08 │
│   k_rna │       3.24e-07 │           6.71e-07 │                  6.60e-07 │                 6.60e-07 │
│k_pro_sat│       4.50e-13 │           4.82e-14 │                  3.21e-11 │                 2.73e-12 │
│k_dna_sat│       8.29e-16 │           2.07e-16 │                  1.45e-12 │                 6.01e-14 │
│k_rna_sat│       1.00e-12 │           2.40e-14 │                  1.76e-10 │                 1.74e-11 │
└─────────┴────────────────┴────────────────────┴───────────────────────────┴──────────────────────────┘
=#
function phyt_params_default(N::Int64, mode::MacroMolecularMode)
    params=Dict(
        "Nsuper"   => [1],       # Number of phyto cells each super individual represents
        "C_DNA"    => [1.8e-13], # DNA C quota of phyto cells (mmolC/cell)
        "mean"     => [1.2],     # Mean of the normal distribution of initial phyto individuals
        "var"      => [0.3],     # Variance of the normal distribution of initial phyto individuals
        "RNA2DNA"  => [1.73],    # Initial RNA:DNA ratio in phytoplankton (mmol C/mmolC) from Micromonas sp.
        "PRO2DNA"  => [34.83],   # Initial protein:DNA ratio in phytoplankton (mmol C/mmolC) from Micromonas sp.
        "CH2DNA"   => [23.3],    # Initial carbohydrate+lipid:DNA ratio in phytoplankton (mmol C/mmolC) from Micromonas sp.
        "Chl2DNA"  => [3.45],    # Initial Chla:DNA ratio in phytoplankton (mmolC/mmolC) from Micromonas sp.
        "α"        => [2.0e-2],  # Irradiance absorption coeff (m²/mgChl)
        "Φ"        => [4.0e-5],  # Maximum quantum yield (mmolC/μmol photon)
        "T⁺"       => [308.0],   # Maximal temperature for growth (K)
        "Ea"       => [5.3e4],   # Free energy
        "PCmax"    => [6.2e-5],  # Maximum primary production rate (per second)
        "VDOCmax"  => [0.0],     # Maximum DOC uptake rate (mmol C/mmol C/second)
        "VNH4max"  => [6.9e-6],  # Maximum N uptake rate (mmol N/mmol C/second)
        "VNO3max"  => [6.9e-6],  # Maximum N uptake rate (mmol N/mmol C/second)
        "VPO4max"  => [1.2e-6],  # Maximum P uptake rate (mmol P/mmol C/second)
        "KsatNH4"  => [0.005],   # Half-saturation coeff (mmol N/m³)
        "KsatNO3"  => [0.010],   # Half-saturation coeff (mmol N/m³)
        "KsatPO4"  => [0.003],   # Half-saturation coeff (mmol P/m³)
        "KsatDOC"  => [0.0],     # Half-saturation coeff (mmol C/m³)
        "NSTmax"   => [0.12],    # Maximum N reserve in total N (mmol N/mmol N)
        "NSTmin"   => [0.01],    # Minimum N reserve in total N (mmol N/mmol N)
        "PSTmax"   => [0.70],    # Maximum P reserve in total P (mmol P/mmol P)
        "PSTmin"   => [0.40],    # Minimum P reserve in total P (mmol P/mmol P)
        "CHmax"    => [0.4],     # Maximum Carbohydrate in cell (mmol C/mmol C)
        "CHmin"    => [0.1],     # Minimum Carbohydrate in cell (mmol C/mmol C)
        "respir"   => [1.2e-6],  # Respiration rate(per second)
        "k_pro"    => [6.0e-5],  # Protein synthesis rate (mmol C/mmol C/second)
        "k_sat_pro"=> [4.5e-13], # Hafl saturation constent for protein synthesis (mmol C/cell)
        "k_dna"    => [1.6e-7],  # DNA synthesis rate (mmol C/mmol C/second)
        "k_sat_dna"=> [1.0e-15], # Hafl saturation constent for DNA synthesis (mmol C/cell)
        "k_rna"    => [3.0e-7],  # RNA synthesis rate (mmol C/mmol C/second)
        "k_sat_rna"=> [1.0e-12], # Hafl saturation constent for RNA synthesis (mmol C/cell)
        "Chl2N"    => [3.0],     # Maximum Chla:N ratio in phytoplankton
        "R_NC_PRO" => [1/4.5],   # N:C ratio in protein (from Inomura et al 2020.)
        "R_NC_DNA" => [1/2.9],   # N:C ratio in DNA (from Inomura et al 2020.)
        "R_PC_DNA" => [1/11.1],  # P:C ratio in DNA
        "R_NC_RNA" => [1/2.8],   # N:C ratio in RNA (from Inomura et al 2020.)
        "R_PC_RNA" => [1/10.7],  # P:C ratio in RNA
        "dvid_P"   => [1.0e-5],  # Division probability per second
        "grz_P"    => [0.0],     # Grazing probability per second
        "mort_P"   => [5e-5],    # Probability of cell natural death per second
        "mort_reg" => [0.5],     # Regulation of cell natural death
        "grazFracC"=> [0.7],     # Fraction goes into dissolved organic pool
        "grazFracN"=> [0.7],     # Fraction goes into dissolved organic pool
        "grazFracP"=> [0.7],     # Fraction goes into dissolved organic pool
        "mortFracC"=> [0.5],     # Fraction goes into dissolved organic pool
        "mortFracN"=> [0.5],     # Fraction goes into dissolved organic pool
        "mortFracP"=> [0.5],     # Fraction goes into dissolved organic pool
    )

    if N == 1
        return params
    else
        return generate_n_species_params(N, params)
    end
end

"""
    phyt_params_default(N::Int64, mode::AbstractMode)
Generate default phytoplankton parameter values based on `AbstractMode` and species number `N`.
"""
function phyt_params_default(N::Int64, mode::QuotaMode)
    params=Dict(
        "Nsuper"   => [1],       # Number of phyto cells each super individual represents
        "Cquota"   => [1.8e-11], # C quota of phyto cells at size = 1.0
        "mean"     => [1.2],     # Mean of the normal distribution of initial phyto individuals
        "var"      => [0.3],     # Variance of the normal distribution of initial phyto individuals
        "Chl2Cint" => [0.10],    # Initial Chla:C ratio in phytoplankton (mgChl/mmolC)
        "α"        => [2.0e-2],  # Irradiance absorption coeff (m²/mgChl)
        "Φ"        => [4.0e-5],  # Maximum quantum yield (mmolC/μmol photon)
        "T⁺"       => [305.0],   # Maximal temperature for growth
        "Ea"       => [5.3e4],   # Free energy
        "PCmax"    => [4.2e-5],  # Maximum primary production rate (per second)
        "VDOCmax"  => [0.0],     # Maximum DOC uptake rate (mmol C/mmol C/second)
        "VNH4max"  => [6.9e-6],  # Maximum N uptake rate (mmol N/mmol C/second)
        "VNO3max"  => [6.9e-6],  # Maximum N uptake rate (mmol N/mmol C/second)
        "VPO4max"  => [1.2e-6],  # Maximum P uptake rate (mmol P/mmol C/second)
        "PC_b"     => [0.6],     # Shape parameter for size
        "VDOC_b"   => [0.6],     # Shape parameter for size
        "VN_b"     => [0.6],     # Shape parameter for size
        "VP_b"     => [0.6],     # Shape parameter for size
        "KsatNH4"  => [0.005],   # Half-saturation coeff (mmol N/m³)
        "KsatNO3"  => [0.010],   # Half-saturation coeff (mmol N/m³)
        "KsatPO4"  => [0.003],   # Half-saturation coeff (mmol P/m³)
        "KsatDOC"  => [0.0],     # Half-saturation coeff (mmol C/m³)
        "Nqmax"    => [0.12],    # Maximum N quota in cell (mmol N/mmol C)
        "Nqmin"    => [0.05],    # Minimum N quota in cell (mmol N/mmol C)
        "Pqmax"    => [0.01],    # Maximum P quota in cell (mmol P/mmol C)
        "Pqmin"    => [0.004],   # Minimum P quota in cell (mmol P/mmol C)
        "Cqmax"    => [0.4],     # Maximum C quota in cell (mmol C/mmol C)
        "Cqmin"    => [0.1],     # Minimum C quota in cell (mmol C/mmol C)
        "k_mtb"    => [3.5e-5],  # Metabolic rate (per second)
        "k_mtb_b"  => [0.25],    # Metabolic rate
        "respir_a" => [1.2e-6],  # Respiration rate(per second)
        "respir_b" => [0.6],     # Shape parameter for size
        "Chl2N"    => [3.0],     # Maximum Chla:N ratio in phytoplankton
        "R_NC"     => [16/106],  # N:C ratio in cell biomass
        "R_PC"     => [1/106],   # N:C ratio in cell biomass
        "grz_P"    => [0.0],     # Grazing probability per second
        "dvid_P"   => [5e-5],    # Probability of cell division per second.
        "dvid_type"=> [1],       # The type of cell division, 1:sizer, 2:adder.
        "dvid_stp" => [6.0],     # Steepness of sigmoidal function
        "dvid_reg" => [1.9],     # Regulations of cell division (cell size)
        "dvid_stp2"=> [2.0],     # Steepness of sigmoidal function
        "dvid_reg2"=> [12.0],    # Regulations of cell division (clock time)
        "mort_P"   => [5e-5],    # Probability of cell natural death per second
        "mort_reg" => [0.5],     # Regulation of cell natural death
        "grazFracC"=> [0.7],     # Fraction goes into dissolved organic pool
        "grazFracN"=> [0.7],     # Fraction goes into dissolved organic pool
        "grazFracP"=> [0.7],     # Fraction goes into dissolved organic pool
        "mortFracC"=> [0.5],     # Fraction goes into dissolved organic pool
        "mortFracN"=> [0.5],     # Fraction goes into dissolved organic pool
        "mortFracP"=> [0.5],     # Fraction goes into dissolved organic pool
        "ther_mort"=> [0],       # thermal mortality, 1 for on, 0 for off
    )

    if N == 1
        return params
    else
        return generate_n_species_params(N, params)
    end
end

function phyt_params_default(N::Int64, mode::CarbonMode)
    params=Dict(
        "Nsuper"   => [1],       # Number of phyto cells each super individual represents
        "Cquota"   => [1.8e-11], # C quota of phyto cells at size = 1.0
        "mean"     => [1.2],     # Mean of the normal distribution of initial phyto individuals
        "var"      => [0.3],     # Variance of the normal distribution of initial phyto individuals
        "Chl2C"    => [0.10],    # Chla:C ratio in phytoplankton (mgChl/mmolC)
        "α"        => [2.0e-2],  # Irradiance absorption coeff (m²/mgChl)
        "Φ"        => [4.0e-5],  # Maximum quantum yield (mmolC/μmol photon)
        "T⁺"       => [298.0],   # Maximal temperature for growth
        "Ea"       => [5.3e4],   # Free energy
        "PCmax"    => [4.2e-5],  # Maximum primary production rate (per second)
        "PC_b"     => [0.6],     # Shape parameter for size
        "respir_a" => [1.2e-6],  # Respiration rate(per second)
        "respir_b" => [0.6],     # Shape parameter for size
        "grz_P"    => [0.0],     # Grazing probability per second
        "dvid_P"   => [5e-5],    # Probability of cell division per second.
        "dvid_type"=> [1],       # The type of cell division, 1:sizer, 2:adder.
        "dvid_stp" => [6.0],     # Steepness of sigmoidal function
        "dvid_reg" => [1.9],     # Regulations of cell division (sizer)
        "dvid_stp2"=> [2.0],     # Steepness of sigmoidal function
        "dvid_reg2"=> [12.0],    # Regulations of cell division (sizer)
        "mort_P"   => [5e-5],    # Probability of cell natural death per second
        "mort_reg" => [0.5],     # Regulation of cell natural death
        "grazFracC"=> [0.7],     # Fraction goes into dissolved organic pool
        "mortFracC"=> [0.5],     # Fraction goes into dissolved organic pool
        "ther_mort"=> [0],       # thermal mortality, 1 for on, 0 for off
    )

    if N == 1
        return params
    else
        return generate_n_species_params(N, params)
    end
end

function generate_n_species_params(N, params)
    p = []
    for key in keys(params)
        pa = (key, fill(params[key][1],N))
        push!(p, pa)
    end
    return Dict(p)
end

ext: jl


text:
# # Custom Affine Transformations
#
# This page explains how to implement your own affine transformation operations by way of
# example.
# You can add your own custom affine transformations into Stheno using the same
# mechanism as all of the existing transformations (addition, multiplication, composition,
# etc).
# First, load up the relevant packages.

using AbstractGPs
using LinearAlgebra
using Stheno

# ## The Affine Transformation

# Suppose that, for some reason, you wish to implement the affine transformation of a single
# process `f` given by `(Af)(x) = f(x) + f(x + 3) - 2`.
# In order to define this transformation, first create a function which accepts `f` and
# returns a `CompositeGP`:
using Stheno: AbstractGP, CompositeGP, SthenoAbstractGP

A(f::SthenoAbstractGP) = CompositeGP((A, f), f.gpc)

# The first argument to `CompositeGP` contains `A` itself and any data needed to fully
# specify the process results from this transformation. In this case the only piece of
# information required is `f`, but really any data can be put in this argument.
# For example, if we wished to replace the translation of `-3` by a parameter, we could do
# so, and make it a part of this first argument.
#
# The second argument is the book-keeping object that needs to be passed around in order to
# know how to compute covariances properly.
# It is because of this argument that we restrict the accepted GPs to be `SthenoAbstractGP`,
# as we can safely assume that these have a `gpc` field.

# We'll now define a type alias in order to simplify some methods later on:
const A_args = Tuple{typeof(A), SthenoAbstractGP};

# ## Most Important Methods

# We must now define methods of three functions on
# `A_args`: `mean`, `cov`, and `var`.
# First the `mean` -- this method should accept both an `A_args` and an `AbstracVector`, and
# return the mean vector of `A(f)` at `x`.
# Some textbook calculations reveal that this is

Stheno.mean((A, f)::A_args, x::AbstractVector) = mean(f, x) .+ mean(f, x .+ 3) .- 2

# The first argument here is _always_ going to be precisely the tuple of arguments passed
# into the `CompositeGP` constructor above.
# You can assume that you can compute any statistics of `f` that the AbstractGPs API
# provides.

# We now turn our attention to `cov`. The first method we consider is
# `cov(args::A_args, x::AbstractVector, y::AbstractVector)`, which should return the
# cross-covariance matrix between all pairs of points in `x` and `y` under the transformed
# process, `A(f)`.
# Again, some standard manipulations reveal that this covariance is given by

function Stheno.cov((A, f)::A_args, x::AbstractVector, y::AbstractVector)
    return cov(f, x, y) + cov(f, x, y .+ 3) + cov(f, x .+ 3, y) + cov(f, x .+ 3, y .+ 3)
end

# The last substantially new method to implement is
# `cov(args::A_args, g::AbstractGP, x::AbstractVector, y::AbstractVector)`, which should
# return the cross-covariance matrix between `A(f)` at `x` and `g` at `y`.
# When implementing this method, you can assume you have access to functions like
# `cov(f, g, x, y)` etc:

function Stheno.cov((A, f)::A_args, g::AbstractGP, x::AbstractVector, y::AbstractVector)
    return cov(f, g, x, y) + cov(f, g, x .+ 3, y)
end

# ## Additional (Required) Methods

# There are a number of other methods that you should implement.
# These are all just special cases or slight modifications of the three methods above, and
# should be straightforward to implement given that you've implemented the above methods.

# First, lets build a GPPP containing an instance of our transformation so that some
# properties can be verified.
# The definition of the methods being implemented is demonstrated by checking an
# equality after defining each method.

gppp = @gppp let
    f = GP(SEKernel())
    Af = A(f)
end

# Also create some input vectors.

x_f = GPPPInput(:f, randn(3))
y_f = GPPPInput(:f, randn(6))
x_Af = GPPPInput(:Af, randn(3))
y_Af = GPPPInput(:Af, randn(6))
z_Af = GPPPInput(:Af, randn(3))

# The covariance matrix at a single pair of inputs:

function Stheno.cov((A, f)::A_args, x::AbstractVector)
    return cov(f, x) + cov(f, x, x .+ 3) + cov(f, x .+ 3, x) + cov(f, x .+ 3)
end

cov(gppp, x_Af, x_Af) ≈ cov(gppp, x_Af)

# The diagonal of the covariance matrix at a single pair of inputs:

function Stheno.var((A, f)::A_args, x::AbstractVector)
    return var(f, x) + var(f, x .+ 3) + var(f, x, x .+ 3) + var(f, x .+ 3, x)
end

var(gppp, x_Af) ≈ diag(cov(gppp, x_Af))

# The diagonal of the cross-covariance matrix for equal-length inputs:

function Stheno.var((A, f)::A_args, x::AbstractVector, y::AbstractVector)
    return var(f, x, y) + var(f, x, y .+ 3) + var(f, x .+ 3, y) + var(f, x .+ 3, y .+ 3)
end

var(gppp, x_Af, z_Af) ≈ diag(cov(gppp, x_Af, z_Af))

# The diagonal of the cross-covariance between different processes for equal-length inputs:

function Stheno.var((A, f)::A_args, g::AbstractGP, x::AbstractVector, y::AbstractVector)
    return var(f, g, x, y) + var(f, g, x .+ 3, y)
end

var(gppp, x_Af, x_f) ≈ diag(cov(gppp, x_Af, x_f))

# `cov` and `var` between processes when `Af`'s arguments are the second argument, rather
# than the first:

function Stheno.cov(g::AbstractGP, (A, f)::A_args, x::AbstractVector, y::AbstractVector)
    return cov(g, f, x, y) + cov(g, f, x, y .+ 3)
end

cov(gppp, x_f, x_Af) ≈ cov(gppp, x_Af, x_f)'

function Stheno.var(g::AbstractGP, (A, f)::A_args, x::AbstractVector, y::AbstractVector)
    return var(g, f, x, y) + var(g, f, x, y .+ 3)
end

var(gppp, x_f, x_Af) ≈ var(gppp, x_Af, x_f)

# ## Checking Your Implementation

# Given the numerous methods above, it's a really good idea to utilise the functionality
# provided by AbstractGPs.jl to check that you've implemented them all consistently with one
# another.

using AbstractGPs.TestUtils: test_internal_abstractgps_interface
using Random

rng = MersenneTwister(123456);
test_internal_abstractgps_interface(rng, gppp, x_Af, y_Af);
test_internal_abstractgps_interface(rng, gppp, x_Af, y_f);
test_internal_abstractgps_interface(rng, gppp, x_f, y_Af);

# Roughly speaking, provided that you've implemented the first three methods correctly, this
# test ought to catch any glaring problems if you've made a mistake in the rest.
# If course, it won't check that your implementations of the first three methods correctly
# implement the desired affine transformation, so you should write whatever tests you need
# in order to convince yourself of that.


text:
function documentationproject(
        m::Module;
        inlineincludes = false,
        executecode = true,
        refmodules = (m,),
        watchpackage = true,
        assetdir = p"assets",
    )
    dir = Path(pkgdir(m))
    doctree = Pollen.loaddoctree(joinpath(dir, "toc.md"))

    rewriters = Rewriter[]

    watchpackage && push!(rewriters, PackageWatcher([m]))

    push!(rewriters, DocumentFolder(dir))

    # Change each document :body into an :article
    push!(rewriters, Replacer(x -> Pollen.withtag(x, :article), SelectTag(:body)))

    # Add slugified IDs to all :h1, :h2, :h3, and :h4 tags
    push!(rewriters, AddID(SEL_H1234))

    # Run cells marked with "cell" attributes and include results in doc, as in Publish.jl
    executecode && push!(rewriters, ExecuteCode())

    # Find references to defined symbols in package and create reference pages for them
    push!(rewriters, Referencer(refmodules))

    # Insert each doc into a HTML template, and either link to CSS and JS or inline them
    push!(rewriters,
        HTMLTemplater(
            joinpath(ASSETDIR, "hugobook.html"),
            [p"hugobook.css", p"ansi.css",
            p"loadhighlightjs.js", p"julia.min.js", p"highlight.min.js"];
            assetdir = ASSETDIR,
            inlineincludes = inlineincludes,
            insertpos = NthChild(2, SelectAttrEq(:class, "book-page"))
        ))

    isdir(joinpath(dir, assetdir)) && push!(rewriters, Assets(joinpath(dir, assetdir)))

    # Insert table of contents, a sidebar with document tree, and project title into
    # the template
    push!(rewriters, Inserter([
        (Pollen.toccreator(hierarchysels = SelectTag.((:h1, :h2, :h3, :h4))), FirstChild(SelectAttrEq(:id, "toc"))),
        ((p, x) -> doctree, FirstChild(SelectAttrEq(:id, "sidebar"))),
        ((p, x) -> XLeaf("$m.jl"), FirstChild(SelectAttrEq(:id, "title"))),
        (createtitle, FirstChild(SelectTag(:head))),
    ]))

    # Add "html" extensions to all links, e.g. ".md" -> ".md.html"
    push!(rewriters, ChangeLinkExtension("html"))

    # Replace absolute links by relative ones so that they are properly resolved, even
    # if output files are served under a subresource like "/Package.jl/dev/..."
    push!(rewriters, RelativeLinks())
    push!(rewriters, RelativeLinks(:script, :src))
    push!(rewriters, RelativeLinks(:link, :href))

    # Change non-html tags like :toc to a :div[class="toc"]
    push!(rewriters, HTMLify())

    # Include fonts
    push!(rewriters, Assets(joinpath(ASSETDIR, p"fonts"), dstdir = p"fonts/"))

    push!(rewriters, HTMLRedirect(p"README.md"))

    project = Project(rewriters)
end

function serve(m::Module; kwargs...)
    project = documentationproject(m)
    serve(project; kwargs...)
end

"""
    serve(project)
    serve(module)
"""
function serve(project::Project; lazy = true)
    builder = FileBuilder(HTML(), Path(mktempdir()))
    server = Server(project, builder)
    mode = lazy ? ServeFilesLazy() : ServeFiles()
    runserver(server, mode)
end

const ASSETDIR = joinpath(Path(pkgdir(Pollen)), p"static")


text:
using Pkg

Pkg.add("StatsBase")
Pkg.add("HDF5")
Pkg.add("DataStructures")
Pkg.add("NearestNeighbors")
Pkg.add("JSON")
Pkg.add("Serialization")
Pkg.add("CSV")
Pkg.add("DataFrames")

max_stars_repo_name: lzzppp/DERT

max_issues_repo_name: lzzppp/DERT

max_forks_repo_name: lzzppp/DERT


text:
# Configuration for running simple neutral copy error simulation
#=
Recommended command line to run:
>  julia neutral.jl examples/sn_example1
=#
export simtype
@everywhere simtype = 3
@everywhere const N = 10        # population size
const mutstddev = 0.04
const ngens = 1001
const initial_value = 1.0
const num_trials=1000
const record_interval = 200
const log_error = false
const wright_fisher_copy = false
const conformist_probability=0.3


text:
dropdown_options = [
    Dict("label" => "Blobs", "value" => "make_blobs"),
    Dict("label" => "Gaussian Quantiles", "value" => "make_gaussian_quantiles"),
    Dict("label" => "Moons", "value" => "make_moons"),
    Dict("label" => "Two Spirals", "value" => "make_twospirals"),
    Dict("label" => "Circles", "value" => "make_circles"),
]

app.layout = html_div(id="page-index") do
    html_div(id="container") do
        html_div(id="header") do
            html_a(href="http://www.ufrrj.br/", target="_blank") do
                html_img(id="logo", src="/assets/logo.png", alt="logo")
            end,
            html_a(href="https://atislabs.com.br/", target="_blank") do
                html_h2("ATISLabs")
            end
        end,
        html_hr(),
        html_div(id="main") do
            html_h1("Synthetic Datasets Dashboard"),
            html_div(id = "link") do
                dcc_dropdown(id="drop", options = dropdown_options, value = "make_blobs")
            end,
            dcc_location(id="url", refresh = false),
            html_div(id="page-content")
        end,
        html_hr(),
        html_div(id="footer") do
            html_a(href="http://www.dcc.ufrrj.br", target="_blank") do
                html_p("© DCC/UFRRJ 2020")
            end
        end
    end
end

callback!(app,Output("page-content", "children"),[Input("drop","value")]) do drop
    dropload(drop)
end


text:
## TYPES TO BE SUBTYPED

abstract type TuningStrategy <: MLJBase.MLJType end
MLJBase.show_as_constructed(::Type{<:TuningStrategy}) = true

## METHODS TO BE IMPLEMENTED

# for initialization of state (compulsory)
setup(tuning::TuningStrategy, model, range, verbosity) = range

# for adding extra user-inspectable information to the history:
extras(tuning::TuningStrategy, history, state, E) = NamedTuple()

# for generating batches of new models and updating the state (but not
# history):
function models end

# for adding to the default report:
tuning_report(tuning::TuningStrategy, history, state) = NamedTuple()

# for declaring the default number of models to evaluate:
default_n(tuning::TuningStrategy, range) = DEFAULT_N

# for encoding the selection_heuristics supported by given strategy:
supports_heuristic(heuristic::Any, strategy::Any) = false


text:
##                                                                            ##
##              Problem set #1 for Structural modelling course                ##
##                               Nicolás Urdaneta                             ##
##                                                                            ##

## Install & load packages

#=
Pkg.add("JLD")
Pkg.add("JLD2")
Pkg.add("Random")
Pkg.add("LinearAlgebra")
Pkg.add("Statistics")
Pkg.add("CSV")
Pkg.add("FreqTables")
Pkg.add("DataFrames")
Pkg.add("Distributions")
=#
#Pkg.add("HDF5")
using JLD
using JLD2
using Random
using LinearAlgebra
using Statistics
using CSV
using FreqTables
using DataFrames
using Distributions

## 0. Github account
    ## Nurdaneta
    ## Forked

## 1. Initializing variables and basic matrix operations

    # a) Matrices of random numbers

Random.seed!(1234)

function q1()
    # i.-iv. create matrices
    A = rand(Uniform(-5,10), 10,7)
    # Or A = 15*rand(10,7).-5
    B = rand(Normal(-2,15), 10,7)
    # Or B = 15*randn(10,7).-2
    C = [A[1:5,1:5] B[1:5,6:7]]
    # Or C = cat(A[1:5,1:5],B[1:5,end-1:end]; dims=2)
    D = copy(A)
    D[D.>0] .= 0
    # Or D = A.*(A.<=0)

    # b) List number elements of A
    println(length(A))

    # c) List unique elements of D
    println(length(unique(D)))

    # d) reshape function, create E: vec operator to B
    E = reshape(B, (70,1))
    E = vec(B)
        # Or E = reshape(B,length(B))
        # Or E = B[:]

    # e) 3 dimensional matrix
    F = cat(A,B;dims=3)

    # f) permute
    F = permutedims(F,[3 1 2])

    # g) kronecker product
    G = kron(B,C)
    # Will give error: G = kron(C,F)

    # h) saving
    save("matrixpractice.jld","A",A,"B",B,"C",C,"D",D,"E",E,"F",F,"G",G)

    # i) saving
    save("firstmatrix.jld","A",A,"B",B,"C",C,"D",D)

    # j) export to csv
    CSV.write("Cmatrix.csv",DataFrame(C))

    # export to DAT
    CSV.write("Dmatrix.dat",DataFrame(D); delim="\t")

    # Function around code
    return A,B,C,D
end

    ## 2. Loops and comprehensions

    function q2(A,B,C)
        # a)
    AB = Array{Float64}(undef, 10, 7)
    for i in 1:10
        for j in 1:7
            AB[i,j] = A[i,j]*B[i,j]
        end
    end
    # OR:     AB = [A[i,j]*B[i,j] for i=1:size(A,1),j=1:size(A,2)]

    AB2 = A.*B

        # b)
    Cprime = [C[i] for i in eachindex(C) if (C[i]<=5 && C[i]>=-5)]
    Cprime2 = C[(C.>=-5) .& (C.<=5)]
        # c)
    X = Array{Float64}(undef, 15169, 6, 5)
    X[1:15169,1,1:5] .= 1

    for i in 1:5
        X[1:15169,2,i] = rand(Uniform(0,1), 15169,1)
        X[1:15169,2,i]  = (X[1:15169,2,i] .< 0.75*(6-i)/5)

        X[1:15169,3,i] = rand(Normal(15+i-1,5(i-1)), 15169,1)

        X[1:15169,4,i] = rand(Normal(pi*(6-i)/3,1/exp(1)), 15169,1)

        X[1:15169,5,i] = rand(Binomial(20,0.6),15169)

        X[1:15169,6,i] = rand(Binomial(20,0.5),15169)
    end

    #= OR:
    N = 15_169
    K = 6
    T = 5
    X = cat([cat([ones(N,1) rand(N,1).<=(0.75*(6-t)/5) (15+t-1).+(5*(t-1)).*randn(N,1) (π*(6-t)/3).+(1/exp(1)).*randn(N,1) rand(Binomial(20,0.6),N) rand(Binomial(20,0.5),N)];dims=3) for t=1:T]...;dims=3) # discrete_normal binomial
    =#

        # d)
    β =  Array{Float64}(undef, 6, 5)
    β[1,1:5] = [(1+0.25*(i-1)) for i in 1:5 ]
    β[2,1:5] = [log(i) for i in 1:5]
    β[3,1:5] = [-sqrt(i) for i in 1:5]
    β[4,1:5] = [exp(i)-exp(i+1) for i in 1:5]
    β[5,1:5] = [i for i in 1:5]
    β[6,1:5] = [i/3 for i in 1:5]

    #= OR
    β = vcat([cat([1+0.25*(t-1) log(t) -sqrt(t) exp(t)-exp(t+1) t t/3];dims=1) for t=1:T]...)'
    =#

        # e)
    Y = hcat([cat(X[:,:,t]*β[:,t] + .36*randn(N,1);dims=2) for t=2:T]...)

    return nothing
end

## 3. Reading in Data and calculating summary statistics

function q3()

        # a)
    using CSV; data = CSV.read("nlsw88.csv")
    #save("nlsw88.jld","nlsw88",nlsw88)
    @save "nlsw88.jld2" data
    #using JLD; save("nlsw88.jld", "data", data)

        # b)
    mean(data.never_married)
    mean(data.collgrad)

        # c)
    freqtable(data[!,"race"])
    freqtable(data, :race)

        # d)
        summarystats = describe(data)

        # e)
    freqtable(data[!,"industry"],data[!,"occupation"])
    freqtable(data, :industry, :occupation)

        # f)
    wageonly = data[:,[:industry, :occupation,:wage]]
    grouper = groupby(wageonly, [:industry,:occupation])
    combine(grouper, valuecols(grouper) .=> mean)

    return nothing
end

## 4. Practice with functions

funciton q4()

    # a)
    mats = load("firstmatrix.jld")
    A = mats["A"]
    B = mats["B"]
    C = mats["C"]
    D = mats["D"]

    # b)
    function matrixops(m1,m2)

        if size(m1)!=size(m2)
            error("inputs must have the same size.")
        end

        r1 = m1.*m2
        r2 = m1'*m2
        r3 = sum(m1)+sum(m2)

        return r1, r2, r3

    end

    # d)
    matrixops(A,B)

    # e)
    matrixops(C,D)

    # f)
    matrixops(convert(Array, data.ttl_exp), convert(Array, data.wage))

    return nothing
end


text:
export lanczosBidiag

function lanczosBidiag(A::SparseMatrixCSC{T1,Int},b::Array{T2,1}, k::Int) where {T1,T2}
	T  = promote_type(T1,T2)
	x1 = zeros(T,size(A,1))
	x2 = zeros(T,size(A,2))

	Af(x,flag) = (flag=='F') ? mul!(x1,A,x,1.0,0.0) : mul!(x2,transpose(A),x,1.0,0.0)
	return lanczosBidiag(Af,b,k)
end

lanczosBidiag(A,b,k) = lanczosBidiag((x,flag) -> ((flag=='F') ? A*x : A'*x),b,k)

"""
U, B, V =  lanczosBidiag(A,p::Vector,k::Int)

Lanczos bidiagonalization of matrix A.

Input:

  A       - function computing A*x = A(x,'F') and A'*x = A(x,'T')
  p       - starting vector
  k       - dimension of subspace

Output:

  U,B,V   - Lanczos vectors
"""
function lanczosBidiag(A::Function,p::Vector,k::Int)

	m      = length(p)
	beta   = norm(p)
	u      = p/beta
	ATu    = A(u,'T')

	n      = length(ATu)
	v      = zeros(n)

	# alocate space for Lanczos vectors
	U = zeros(m,k+1)
	V = zeros(n,k)
	B = zeros(k,2)
	# Prepare for Lanczos iteration.
	U[:,1] = u

	for i=1:k # perform Lanczos bidiagonalization with reorthogonalization.
		if i>1;
			ATu = A(u,'T')
		end
		r = ATu - beta*v

		for j=1:i-1 # reorthogonalization of V
			r -= dot(V[:,j],r)*V[:,j]
		end

		alpha = norm(r)
		v = r/alpha
		B[i,2] = alpha
		V[:,i] = v

		Av = A(v,'F')
		p = Av - alpha*u

		for j=1:i; p = p - dot(U[:,j],p)*U[:,j]; end
		beta = norm(p)
		u    = p/beta

		B[i,1]   = beta
		U[:,i+1] = u
	end
	II,JJ,VV = SparseArrays.spdiagm_internal(-1=>B[:,1], 0=>B[:,2])
	B = sparse(II,JJ,VV,k+1,k)
	# Bt = spdiagm((B[:,1],B[:,2]),[-1,0],k+1,k)
	return U, B, V
end


text:
type Grid
	obj::Ptr{Void}
end

function Grid()
    obj = Ref{Ptr{Void}}(C_NULL)
        err = ccall(("ElDefaultGrid", libEl), Cuint,
            (Ref{Ptr{Void}},), obj)
    err == 0 || throw(ElError(err))
    return Grid(obj[])
end


text:
module MLToolkit

using LinearAlgebra, Random

export Σ

const Σ = sum

# declare Matrix class and constructors
include("Matrix.jl")
# declare SupervisedLearner abstract type
include("SupervisedLearner.jl")
# get arg parsing stuff
include("ParseArgs.jl")
# BaselineLearner class
include("BaselineLearner.jl")

learners = Dict("baseline" => BaselineLearner)

function getlearner(model::AbstractString, args::Dict{Symbol,AbstractString})::SupervisedLearner
	if haskey(learners, model)
		# pass extra command line args to the learner constructor as keyword arguments
		learners[model](;args...)
	else
		error("Unrecognized model: $model. Supported models are $(join(keys(learners), ", ", ", and "))")
	end
end

function run()
	# If you want to use a specific seed for the RNG, use `--seed <seed>` on the command line
	Random.seed!() # this seeds the global RNG. If you pass in seed it gets reseeded in parseargs
	args = parseargs(ARGS)
	learner = getlearner(args.learner, args.other)
	evalmode = args.evaluation

	matrix = loadarff(args.arff)
	if args.normalize
		@info "Using normalized data"
		extrema = normalize(matrix)
	end
	# load (and normalize) test data as well if the evalmode is Static
	if evalmode isa Static
		evalmode.data = loadarff(evalmode.filename)
		if args.normalize
			normalize(evalmode.data, extrema)
		end
	end

	@info "Dataset" matrix.datasetname rows(matrix) columns(matrix) args.learner evalmode

	# precompile so that the compilation time doesn't skew the first runtime, if we're timing anything
	precompile(train, (typeof(learner), Matrix, Matrix))
	precompile(measureaccuracy, (typeof(learner), Matrix, Matrix))
	precompile(measureaccuracy, (typeof(learner), Matrix, Matrix, Union{Matrix,Nothing}))
	trainandtest(learner, matrix, evalmode)
end

function trainandtest(learner::SupervisedLearner, data::Matrix, ::Training)
	@info "Calculating accuracy on training set"
	features = copymatrix(data, 1:rows(data), 1:columns(data)-1)
	labels = copymatrix(data, 1:rows(data), columns(data))
	confusion = initconfusionmatrix(labels)
	elapsedtime = @elapsed train(learner, features, labels)
	accuracy = measureaccuracy(learner, features, labels, confusion)
	@info "Results" timetotrain=elapsedtime accuracy
	@debug "Confusion matrix: (Row=target value, Col=predicted value)" confusion
end

function trainandtest(learner::SupervisedLearner, data::Matrix, evalmode::Static)
	testdata = evalmode.data
	@info "Calculating accuracy on separate test set" evalmode.filename rows(testdata)
	features = copymatrix(data, 1:rows(data), 1:columns(data)-1)
	labels = copymatrix(data, 1:rows(data), columns(data))
	elapsedtime = @elapsed train(learner, features, labels)
	trainaccuracy = measureaccuracy(learner, features, labels)
	testfeatures = copymatrix(testdata, 1:rows(testdata), 1:columns(testdata)-1)
	testlabels = copymatrix(testdata, 1:rows(testdata), columns(testdata))
	confusion = initconfusionmatrix(labels)
	testaccuracy = measureaccuracy(learner, testfeatures, testlabels, confusion)
	@info "Results" timetotrain=elapsedtime trainaccuracy testaccuracy
	@debug "Confusion matrix: (Row=target value, Col=predicted value)" confusion
end

function trainandtest(learner::SupervisedLearner, data::Matrix, evalmode::Rando)
	trainpercent = evalmode.percenttest
	if trainpercent < 0 || trainpercent > 1
		error("Percentage for random evaluation must be between 0 and 1")
	end
	@info "Calculating accuracy on a random hold-out set" trainpercent testpercent=1-trainpercent
	shuffle!(data)
	trainsize = trunc(Int, trainpercent * rows(data))
	trainfeatures = copymatrix(data, 1:trainsize, 1:columns(data)-1)
	trainlabels = copymatrix(data, 1:trainsize, columns(data))
	testfeatures = copymatrix(data, trainsize+1:rows(data), 1:columns(data)-1)
	testlabels = copymatrix(data, trainsize+1:rows(data), columns(data))
	elapsedtime = @elapsed train(learner, trainfeatures, trainlabels)
	trainaccuracy = measureaccuracy(learner, trainfeatures, trainlabels)
	confusion = initconfusionmatrix(testlabels)
	testaccuracy = measureaccuracy(learner, testfeatures, testlabels, confusion)
	@info "Results" timetotrain=elapsedtime trainaccuracy testaccuracy
	@debug "Confusion matrix: (Row=target value, Col=predicted value)" confusion
end

function trainandtest(learner::SupervisedLearner, data::Matrix, evalmode::Cross)
	folds = evalmode.numfolds
	folds <= 0 && error("Number of folds must be greater than 0")
	@info "Calculating accuracy using cross-validation" folds
	sumaccuracy = 0.0
	elapsedtime = 0.0
	shuffle!(data)
	for i in 0:folds-1
		begin_i = i * div(rows(data), folds)
		end_i = (i + 1) * div(rows(data), folds)
		trainfeatures = copymatrix(data, [1:begin_i; end_i+1:rows(data)], 1:columns(data)-1)
		trainlabels = copymatrix(data, [1:begin_i; end_i+1:rows(data)], columns(data))
		testfeatures = copymatrix(data, begin_i+1:end_i, 1:columns(data)-1)
		testlabels = copymatrix(data, begin_i+1:end_i, columns(data))
		elapsedtime += @elapsed train(learner, trainfeatures, trainlabels)
		accuracy = measureaccuracy(learner, testfeatures, testlabels)
		sumaccuracy += accuracy
		@info "Fold $(i+1)" testset=begin_i+1:end_i accuracy
	end
	@info "Average results" timetotrain=elapsedtime/folds accuracy=sumaccuracy/folds
end

function initconfusionmatrix(labelmatrix::Matrix)
	dims = valuecount(labelmatrix, 1)
	confusion = Matrix(dims, dims)
	for i in 1:dims
		setattributename(confusion, i, attributevalue(labelmatrix, 1, i-1))
	end
	confusion
end

end # module

max_stars_repo_name: coreywoodfield/mltoolkit-julia

max_issues_repo_name: coreywoodfield/mltoolkit-julia

max_forks_repo_name: coreywoodfield/mltoolkit-julia


text:
# julia1_4.2.jl Vector and Matrix Operations
#
# Most of the math functions work equally well on Vectors
# and matrices.
# The same goes for operations but some need a preceding
# dot (.) to be corresponding element by element ops
# Those include .* .+ .- .% ./ .\ and .^ or any op with otherwise
#  different normal vector/matrix behavior. .+ and .- were enforced
#  recently for clarity of intent.
#  library functions that are operating on element-wise also require a dot
# The following types have a variety of elementary
# operations and functions with optimized methods:
# Type	            Description
# Diagonal	        Diagonal matrix
# Hermitian	        Hermitian matrix - Square, self conjugate
# UpperTriangular	Upper right triangular matrix
# LowerTriangular	Lower left triangular matrix
# Tridiagonal	    Tridiagonal matrix
# SymTridiagonal	Symmetric tridiagonal matrix
# Bidiagonal	    (true)Upper/(false)Lower bidiagonal matrix

# In addition, The LAPACK library is available.
# and a large variety of matrix factorizations

# Examples:
using LinearAlgebra
A = [2.5 2 -4; 3 -1 -6; -10 2.3 4.4]
println(eigvals(A))
println(eigvecs(A))
# 3 ways to get an inverse
Ainv=inv(A)
A^-1
A\I  # I is the identity/uniform-scaling operator
# checking out the inverse
A*Ainv

Am1=A^2 .+ exp.(A) .- 1

B = [2.5 2 -4; 2 -1 -3; -4 -3 4.4]
sB=Symmetric(B)  # less computation
# solve for x in sB*x=v linear system of equations
v=[1.0,2,3]
x=sB\v
println(x)  # Linear Solutions
#
###
#
# more than 12 different factorizations are available
# Single Valued Decomposition SVD
#
(U,S,V)=svd(B)  # such that U*diagm(S)*V'
println(U)      # Matrix
println(V)      # Matrix
println(S)      # Vector diagonal of matrix
println(U*Diagonal(S)*V') # V' transpose of matrix V
#
# LU decomposition L=lower triangular part,
# U=upper triangular part, p= permutation vector
(L,U,p)=lu(B)
println(L)  # lower triangle
println(U)  # upper triangle
println(p)  # permutation vector for rows
println((L*U)[p,:]) # rearrange rows by p

# Make a random Matrix
R=rand(3,3) .+ Diagonal([1.0, 1, 1])  # random matrix R + unit diag. matrix
Ri=inv(R)            # now invert it

# Sparce matrices - type=SparseMatrixCSC Compressed Sparce Cols
#        colVec       rowVec       Data
using SparseArrays
#          row index    col index    values
S=sparse([1, 2, 3, 3],[1, 3, 1, 2],[3, 3, 2, 0])
S=dropzeros(S)
R=S*I
R*rand(3,3)


text:
using Cat

const Tensor = Array{Float64}
@category TensorNetwork

@arrow TensorNetwork TSum :: (Tensor, Tensor) --> Tensor
@arrow TensorNetwork TProd :: (Tensor, Tensor) --> Tensor
@arrow TensorNetwork Contract :: (Tensor, Int, Int) --> Tensor

function example()
    g = tensor(2) # metric
    R = tensor(4) # Riemann curvature tensor
    T = tensor(2) # Stress-energy tensor

    Ric = Contract(R, 1, 3) # Ricci tensor
    Rsc = Contract(Ric, 1, 2) # Ricci scalar
    EFE = Ric - 0.5*Rsc*g - (8*π)*T
end

@alias TensorNetwork TSum Base.:+ 2
@alias TensorNetwork TProd Base.:* 2
Base.:*(a::Float64, b::TensorNetwork.Arrow) = TProd(a*ones(), b)
Base.:*(a::TensorNetwork.Arrow, b::Float64) = TProd(b*ones(), a)

tensor(nd::Int64) = TensorNetwork.Identity{Tensor{nd}}()

@interpretation Einsum (=>) TensorNetwork begin
    letters::Dict{TensorNetwork.Arrow, Char}
end


text:
using Weave
using InteractiveUtils

const main_dir = joinpath(@__DIR__, "..")
const src_dir = joinpath(main_dir, "script")
const notebook_dir = joinpath(main_dir, "notebook")
const test_dir = joinpath(main_dir, "test")

"""
Generates the notebook and lightened source files for test
"""
function generate_output()
    for folder in readdir(src_dir)
        for file in readdir(joinpath(src_dir, folder))
            if endswith(file, ".jl")
                @info("Building $(folder)/$(file)")
                filename = split(file, ".")[1]
                testpath = joinpath(test_dir, folder, "$file")
                notebookpath = joinpath(notebook_dir, folder)
                srcpath = joinpath(src_dir, folder, file)
                try
                    Weave.notebook(srcpath, out_path=notebookpath, timeout=-1, nbconvert_options="--allow-errors")
                    Weave.tangle(srcpath, out_path=testpath)
                catch e
                    @warn("Error weaving $folder/$file:\n$e")
                end
                println("")
            end
        end
    end
end


text:
using SimpleTraits
using Base.Test

trait = SimpleTraits.trait

# @test_throws MethodError trait(4)
@test_throws ErrorException istrait(4)

# definition & adding
@traitdef Tr1{X}
@test trait(Tr1{Int})==Not{Tr1{Int}}
@test !istrait(Tr1{Int})
@traitimpl Tr1{Integer}
@test trait(Tr1{Int})==Tr1{Int}
@test istrait(Tr1{Int})
@test trait(Tr1{Bool})==Tr1{Bool}
@test trait(Tr1{AbstractString})==Not{Tr1{AbstractString}}
@test !istrait(Tr1{AbstractString})

# Logic.  trait(Tr) returns the same trait Tr if it is fulfilled and
# Not{Tr} otherwise.  This is a bit confusing.
@test trait(Tr1{AbstractString})==Not{Tr1{AbstractString}}
@test istrait(Tr1{AbstractString})==false
@test trait(Not{Tr1{AbstractString}})==Not{Tr1{AbstractString}}
@test istrait(Not{Tr1{AbstractString}})==true
@test trait(Not{Not{Tr1{AbstractString}}})==Not{Tr1{AbstractString}}
@test istrait(Not{Not{Tr1{AbstractString}}})==false
@test trait(Not{Not{Not{Tr1{AbstractString}}}})==Not{Tr1{AbstractString}}
@test istrait(Not{Not{Not{Tr1{AbstractString}}}})==true

@test trait(Not{Tr1{Integer}})==Tr1{Integer}
@test istrait(Not{Tr1{Integer}})==false
@test trait(Not{Not{Tr1{Integer}}})==Tr1{Integer}
@test istrait(Not{Not{Tr1{Integer}}})==true
@test trait(Not{Not{Not{Tr1{Integer}}}})==Tr1{Integer}
@test istrait(Not{Not{Not{Tr1{Integer}}}})==false

@traitdef Tr2{X,Y}
@test trait(Tr2{Int,AbstractFloat})==Not{Tr2{Int,AbstractFloat}}
@traitimpl Tr2{Integer, Float64}
@test trait(Tr2{Int, Float64})==Tr2{Int, Float64}
@test trait(Tr2{Int, Float32})==Not{Tr2{Int, Float32}}

# trait functions
@traitfn f{X; Tr1{X}}(x::X) = 1  # def 1
@traitfn f{X; !Tr1{X}}(x::X) = 2
@test f(5)==1
@test f(5.)==2

@traitfn f{X,Y; Tr2{X,Y}}(x::X,y::Y,z) = 1
@test f(5,5., "a")==1
@test_throws MethodError f(5,5, "a")==2
@traitfn f{X,Y; !Tr2{X,Y}}(x::X,y::Y,z) = 2
@test f(5,5, "a")==2

# This will overwrite the definition def1 above
@traitfn f{X; !Tr2{X,X}}(x::X) = 10
@traitfn f{X; Tr2{X,X}}(x::X) = 100
@test f(5)==10
@test f(5.)==10
@traitimpl Tr2{Integer, Integer}
@test f(5.)==10
@test !(f(5)==100)
# need to update method cache:
@traitfn f{X; Tr2{X,X}}(x::X) = 100
@test f(5)==100
@test f(5.)==10

# VarArg
@traitfn vara{X; Tr1{X}}(x::X, y...) = y
@test vara(5, 7, 8)==(7,8)
# @test vara(5.0, 7, 8)==((7,8),) # hangs in lowering because of https://github.com/JuliaLang/julia/issues/13183
@traitfn vara2{X; Tr1{X}}(x::X...) = x
@test vara2(5, 7, 8)==(5, 7, 8)
@test_throws MethodError vara2(5, 7, 8.0)

@traitfn vara3{X; Tr1{X}}(::X...) = X
@test vara3(5, 7, 8)==Int
@test_throws MethodError vara3(5, 7, 8.0)

# with macro
@traitfn @inbounds gg{X; Tr1{X}}(x::X) = x
@test gg(5)==5
@traitfn @generated ggg{X; Tr1{X}}(x::X) = X<:AbstractArray ? :(x+1) : :(x)
@test ggg(5)==5
@traitimpl Tr1{AbstractArray}
@test ggg([5])==[6]

# traitfn with Type
@traitfn ggt{X; Tr1{X}}(::Type{X}, y) = (X,y)
@test ggt(Array, 5)==(Array, 5)

# traitfn with ::X
@traitfn gg27{X; Tr1{X}}(::X) = X
@test gg27([1])==Array{Int,1}

######
# Other tests
#####
include("base-traits.jl")

if VERSION >= v"0.5.0-dev"
    include("backtraces.jl")
end


text:
module Phonon

using AbInitioSoftwareBase: load
using Comonicon: @cast
using Express.PhononWorkflow.Recipes: buildworkflow
using QuantumESPRESSOExpress
using SimpleWorkflows: run!
# using ..ExpressCommands: @load_plugin

"""
Run a `config` file, better with absolute path.

# Arguments
- `config`: the file to be run. Available extensions are `.json`, `.yaml`, `.yml` or `.toml`.
"""
@cast function run(config)
    # @load_plugin
    workflow = buildworkflow(config)
    dict = load(config)
    run!(workflow; saveas = dict["save"]["status"])
    return workflow
end

end


text:
function calculator()
    while true
        calculatInterface()
        println("press c to continue or q to quit")
        response = readline()
        if cmp(response, "c") == 0
            continue
        end
        if cmp(response, "q") == 0
            break
        end
    end
end

function calculatInterface()
    result = 0
    println("Enter the first number")
    a = getEnteredNumber()
    println("Enter The operation to perform")
    operation = readline()
    println("Enter last number")
    b = getEnteredNumber()
    if cmp(operation, "+") == 0
        result = a + b
    end
    if cmp(operation, "-") == 0
        result = a - b
    end
    if cmp(operation, "*") == 0
        result = a * b
    end
    if cmp(operation, "/") == 0
        result = a / b
    end
    if cmp(operation, "^") == 0
        result = a ^ b
    end

    println("resultat: ", result)
end

function getEnteredNumber()
    x = readline()
    x = parse(Float64, x)
    return x
end

calculator()


text:
struct SquaredEuclidean{T1<:AbstractArray,T2<:AbstractArray} <: DistanceTransform
    dt::T1
    v::T2
    z::T1
end

"""
	SquaredEuclidean(
		f::AbstractArray,
		dt = zeros(Float32, size(f)),
		v = ones(Int64, size(f)),
		z = zeros(Float32, size(f) .+ 1)
	)

Prepares an array to be `transform`ed using the squared euclidean algorithm
laid out in 'Distance Transforms of Sampled Functions' [Felzenszwalb and
Huttenlocher] (DOI: 10.4086/toc.2012.v008a019)

# Arguments
- f/img: 1D, 2D, or 3D to be transformed based on location
    to the nearest background (0) pixel
- dt: Empty array that is the size of `f` or `img`. Will be filled
    with the distance transform values of each element in `f` or `img`
- v: `ones(Int64, length(f))` or
    `ones(Int64, size(img))`
- z: `zeros(Float32, length(f) + 1)` or
    `zeros(Float32, size(img) .+ 1)`
"""
function SquaredEuclidean(
    f::AbstractArray,
    dt=zeros(Float32, size(f)),
    v=ones(Int64, size(f)),
    z=zeros(Float32, size(f) .+ 1),
)
    return SquaredEuclidean(dt, v, z)
end

"""
	transform(f::AbstractVector{T}, tfm::SquaredEuclidean)
    transform(img::AbstractMatrix{T}, tfm::SquaredEuclidean)
    transform!(img::AbstractMatrix{T}, tfm::SquaredEuclidean, nthreads)
	transform!(img::CuArray{T,2}, tfm::SquaredEuclidean)

Applies a squared euclidean distance transform to an input image.
Returns an array with spatial information embedded in the array
elements.

# Arguments
- f/img: 1D, 2D, or 3D to be transformed based on location
    to the nearest background (0) pixel
- tfm: `SquaredEuclidean` type
- nthreads: The number of threads on the computer `Threads.nthreads()`.
    Allows you to use a parallelized `transform`
    function if you have access to multiple threads.

# Citation
'Distance Transforms of Sampled Functions' [Felzenszwalb and
Huttenlocher](DOI: 10.4086/toc.2012.v008a019)
"""
function transform(f::AbstractVector{T}, tfm::SquaredEuclidean) where {T}
    dt = tfm.dt
    v = tfm.v
    z = tfm.z
    n = length(f)
    k = 1
    z[1] = -Inf32
    z[2] = Inf32

    # Lower envelope operation
    for q in 2:n
        while true
            s = ((f[q] + q^2) - (f[v[k]] + v[k]^2)) / (2 * q - 2 * v[k])
            if s ≤ z[k]
                k -= 1
            else
                k += 1
                v[k] = q
                z[k] = s
                z[k + 1] = Inf32
                break
            end
        end
    end

    # Distance transform operation
    k = 1
    for q in 1:n
        while z[k + 1] < q
            k = k + 1
        end
        dt[q] = (q - v[k])^2 + f[v[k]]
    end
    return dt
end

# This function is called to for multi-dimensional transforms
function _transform(f, tfm::SquaredEuclidean, dt, v, z) where {T}
    n = length(f)
    k = 1
    z[1] = -Inf32
    z[2] = Inf32

    # Lower envelope operation
    for q in 2:n
        while true
            s = ((f[q] + q^2) - (f[v[k]] + v[k]^2)) / (2 * q - 2 * v[k])
            if s ≤ z[k]
                k -= 1
            else
                k += 1
                v[k] = q
                z[k] = s
                z[k + 1] = Inf32
                break
            end
        end
    end

    # Distance transform operation
    k = 1
    for q in 1:n
        while z[k + 1] < q
            k = k + 1
        end
        dt[q] = (q - v[k])^2 + f[v[k]]
    end
    return dt
end

# This function is called for the GPU version
function _transform(f::AbstractArray{T,1}, dt, v, z) where {T}
    n = length(f)
    k = 1
    z[1] = -1.0f12
    z[2] = 1.0f12

    # Lower envelope operation
    for q in 2:n
        while true
            s = ((f[q] + q^2) - (f[v[k]] + v[k]^2)) / (2 * q - 2 * v[k])
            if s ≤ z[k]
                k -= 1
            else
                k += 1
                v[k] = q
                z[k] = s
                z[k + 1] = 1.0f12
                break
            end
        end
    end

    # Distance transform operation
    k = 1
    for q in 1:n
        while z[k + 1] < q
            k = k + 1
        end
        dt[q] = (q - v[k])^2 + f[v[k]]
    end
    return dt
end

function transform(img::AbstractMatrix{T}, tfm::SquaredEuclidean) where {T}
    rows, columns = size(img)
    dt = tfm.dt
    v = tfm.v
    z = tfm.z
    for x in 1:rows
        dt[x, :] = _transform(img[x, :], tfm, dt[x, :], v[x, :], z[x, :])
    end

    for y in 1:columns
        dt[:, y] = _transform(img[:, y], tfm, dt[:, y], v[:, y], z[:, y])
    end

    return dt
end

function _transform(img::AbstractMatrix{T}, tfm, dt, v, z) where {T}
    rows, columns = size(img)
    for x in 1:rows
        dt[x, :] = _transform(img[x, :], tfm, dt[x, :], v[x, :], z[x, :])
    end

    for y in 1:columns
        dt[:, y] = _transform(img[:, y], tfm, dt[:, y], v[:, y], z[:, y])
    end

    return dt
end

function transform(img::AbstractArray{T, 3}, tfm::SquaredEuclidean) where {T}
    rows, columns, depth = size(img)
    dt = tfm.dt
    v = tfm.v
    z = tfm.z

    for d in 1:depth
        for x in 1:rows
            dt[x, :, d] = _transform(img[x, :, d], tfm, dt[x, :, d], v[x, :, d], z[x, :, d])
        end

        for y in 1:columns
            dt[:, y, d] = _transform(img[:, y, d], tfm, dt[:, y, d], v[:, y, d], z[:, y, d])
        end
    end
    return dt
end

function transform!(img::AbstractMatrix{T}, tfm::SquaredEuclidean, nthreads) where {T}
    @assert nthreads > 1
    dt = tfm.dt
    v = tfm.v
    z = tfm.z
    rows, columns = size(img)
    Threads.@threads for x in 1:rows
        @views _transform(img[x, :], tfm, dt[x, :], fill!(v[x, :], 1), fill!(z[x, :], 0))
    end
    Threads.@threads for y in 1:columns
        @views _transform(img[:, y], tfm, dt[:, y], fill!(v[:, y], 1), fill!(z[:, y], 0))
    end
    return dt
end

function transform!(img::CuArray{T,2}, tfm::SquaredEuclidean) where {T}
    dt = tfm.dt
    v = tfm.v
    z = tfm.z
    rows, columns = size(img)
    @floop CUDAEx() for x in 1:rows
        @views _transform(img[x, :], dt[x, :], fill!(v[x, :], 1), fill!(z[x, :], 0))
    end

    @floop CUDAEx() for y in 1:columns
        @views _transform(img[:, y], dt[:, y], fill!(v[:, y], 1), fill!(z[:, y], 0))
    end
    return dt
end

function transform!(img::CuArray{T,3}, tfm::SquaredEuclidean) where {T}
    dt = tfm.dt
    v = tfm.v
    z = tfm.z
    rows, columns, slices = size(img)
    @floop CUDAEx() for x in 1:rows
        @views _transform(img[x, :, :], dt[x, :, :], fill!(v[x, :, :], 1), fill!(z[x, :, :], 0))
    end

    @floop CUDAEx() for y in 1:columns
        @views _transform(img[:, y, :], dt[:, y, :], fill!(v[:, y, :], 1), fill!(z[:, y, :], 0))
    end

    @floop CUDAEx() for s in 1:slices
        @views _transform(img[:, :, s], dt[:, :, s], fill!(v[:, :, s], 1), fill!(z[:, :, s], 0))
    end
    return dt
end


text:
# Sparsity-promoting least squares migration with source subsampling using the linearized Bregman algorithm
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: December 2016
#

export linbreg_lsrtm

function linbreg_lsrtm(J::TimeModeling.joJacobian, dD::TimeModeling.joData; iterations=10, lambda=[], shots_per_iteration=1, snapshot=[])

# Sparsity promoting LSRTM using the linearized Bregman method

	# Setup Curvelet transform
	C = joCurvelet2D(J.model.n[1], J.model.n[2]; zero_finest=true, DDT=Float32, RDT=Float64)

	# Main loop
	x = zeros(Float32,prod(J.model.n))
	z = zeros(Float32,prod(J.model.n))
	for j=1:iterations
		println("Iteration ", j, " of ", iterations)

		# Subsample shots
		shotIdx = randperm(J.info.nsrc)[1:shots_per_iteration]

		# Subsample observed data
		dDsub = subsample(dD, shotIdx)

		# Subsample Jacobian
		Jsub = subsample(J, shotIdx)

		# Residual
		dDpred = Jsub*x
		residual = dDpred - dDsub

		# Gradient
		gradient = Jsub'*residual

		# Step size
		step = norm(residual)^2/norm(gradient)^2

		# Automatic thresholding
		if isempty(lambda) && j==1
			lambda = convert(Float32,0.1*maximum(abs.(step*C*gradient)))
		end

		# Update dual variable
		z = z - step*gradient;

		# Update primal variable
		x = C'*softThresholding(C*z, lambda)

		# Save snapshots
		if ~isempty(snapshot)
			save(join([snapshot, "_", string(j), ".jld"]), "z", reshape(z,J.model.n), "x", reshape(x,J.model.n))
		end
	end

	return x
end

max_stars_repo_name: SINBADconsortium/opesciSLIM

max_issues_repo_name: SINBADconsortium/opesciSLIM

max_forks_repo_name: SINBADconsortium/opesciSLIM


text:
using Clang
using Clang.Generators

cd(@__DIR__)
# current version checked in is v2.2.9
include_dir = normpath(joinpath(@__DIR__, "RadeonProRenderSDK", "RadeonProRender", "inc"))
# LIBCLANG_HEADERS are those headers to be wrapped.
headers = joinpath.(include_dir, [
    "RadeonProRender_v2.h",
    # "RadeonProRender_GL.h",
    # "RadeonProRender_VK.h"
])

# wrapper generator options
options = load_options(joinpath(@__DIR__, "rpr.toml"))

# add compiler flags, e.g. "-DXXXXXXXXX"
args = get_default_args()
push!(args, "-I$include_dir")

ctx = create_context(headers, args, options)

# run generator
build!(ctx, BUILDSTAGE_NO_PRINTING)

function rewrite!(e::Expr)
    if Meta.isexpr(e, :function)
        func_args = e.args[1].args
        last_arg = func_args[end]
        body = e.args[2]
        ccall_expr = body.args[end]
        type_tuple = ccall_expr.args[4]
        new_body = if startswith(string(last_arg), "out_")
            last_type_ptr = type_tuple.args[end]
            @assert Meta.isexpr(last_type_ptr, :curly)
            @assert last_type_ptr.args[1] == :Ptr
            last_type = last_type_ptr.args[2]
            pop!(func_args) # remove from input args
            :($(last_arg) = Ref{$(last_type)}(); check_error($(ccall_expr)); return $(last_arg)[])
        elseif ccall_expr.args[3] == :rpr_status
            :(check_error($(ccall_expr)))
        else
            ccall_expr
        end
        e.args[2] = new_body
    end
    return e
end

function rewrite!(dag::ExprDAG)
    for node in get_nodes(dag)
        for expr in get_exprs(node)
            rewrite!(expr)
        end
    end
end

rewrite!(ctx.dag)

quote
    function check_error(error_code)
        error_code == RPR_SUCCESS && return
        return error("Error code returned: $(error_code)")
    end
end

cd(@__DIR__)
build!(ctx, BUILDSTAGE_PRINTING_ONLY)


text:
endpoint(::BitbucketAPI, ::typeof(is_member), workspace::String, user::String) =
    Endpoint(:GET, "/workspaces/$workspace/members";
             allow_404=true,
             query=Dict(
                 :q => "user.nickname=\"$user\""
             )
             )
endpoint(::BitbucketAPI, ::typeof(is_member), workspace::String, user::UUID) =
    Endpoint(:GET, "/workspaces/$workspace/members";
             allow_404=true,
             query=Dict(
                 :q => "user.nickname=\"{$user}\""
             )
             )
postprocessor(::BitbucketAPI, ::typeof(is_member)) = DoSomething(ismemberorcollaborator)
into(::BitbucketAPI, ::typeof(is_member)) = Dict

@not_implemented(::BitbucketAPI, ::typeof(groups), )


text:
# Only supports 2D geometries for now

# pt is [x,y] and ring is [[x,y], [x,y],..]
function inring(pt::Vector{Float64}, ring::Vector{Vector{Float64}})
    intersect(i::Vector{Float64},j::Vector{Float64}) =
        (i[2] >= pt[2]) != (j[2] >= pt[2]) && (pt[1] <= (j[1] - i[1]) * (pt[2] - i[2]) / (j[2] - i[2]) + i[1])
    isinside = intersect(ring[1], ring[end])
    for k=2:length(ring)
        isinside = intersect(ring[k], ring[k-1]) ? !isinside : isinside
    end
    isinside
end

# pt is [x,y] and polygon is [ring, [ring, ring, ...]]
function inpolygon(pt::Vector{Float64}, polygon::Vector{Vector{Vector{Float64}}})
    if !inring(pt, polygon[1]) # check if it is in the outer ring first
        return false
    end
    for poly in polygon[2:end] # check for the point in any of the holes
        if inring(pt,poly)
            return false
        end
    end
    true
end

# GeoJSON.Point

inside(point::GeoJSON.Point, polygon::GeoJSON.Polygon) = inpolygon(point.coordinates, polygon.coordinates)

function inside(point::GeoJSON.Point, polygon::GeoJSON.MultiPolygon)
    for poly in polygon.coordinates
        if inpolygon(point.coordinates, poly)
            return true
        end
    end
    false
end

inside(point::GeoJSON.Point, polygon::GeoJSON.Feature) = inside(point, polygon.geometry)

function inside(point::GeoJSON.Point, collection::GeoJSON.GeometryCollection)
    for geom in collection.geometries
        if inside(point, geom)
            return true
        end
    end
    false
end

function inside(point::GeoJSON.Point, collection::GeoJSON.FeatureCollection)
    for feature in collection.features
        if inside(point, feature)
            return true
        end
    end
    false
end

# GeoJSON.MultiPoint
function inside(point::GeoJSON.MultiPoint, polygon::GeoJSON.Polygon)
    poly = polygon.coordinates
    for pt in point.coordinates
        if !inpolygon(pt, poly)
            return false
        end
    end
    true
end

function inside(point::GeoJSON.MultiPoint, polygon::GeoJSON.MultiPolygon)
    for pt in point.coordinates
        for poly in polygon.coordinates
            if !inpolygon(pt, poly)
                return false
            end
        end
    end
    true
end

inside(point::GeoJSON.MultiPoint, polygon::GeoJSON.Feature) = inside(point, polygon.geometry)

function inside(point::GeoJSON.MultiPoint, collection::GeoJSON.GeometryCollection)
    for geom in collection.geometries
        if inside(point, geom)
            return true
        end
    end
    false
end

function inside(point::GeoJSON.MultiPoint, collection::GeoJSON.FeatureCollection)
    for feature in collection.features
        if inside(point, feature)
            return true
        end
    end
    false
end

# GeoJSON.Feature
inside(point::GeoJSON.Feature, polygon::GeoJSON.Polygon) = inside(point.geometry, polygon)
inside(point::GeoJSON.Feature, polygon::GeoJSON.MultiPolygon) = inside(point.geometry, polygon)
inside(point::GeoJSON.Feature, polygon::GeoJSON.Feature) = inside(point.geometry, polygon)
inside(point::GeoJSON.Feature, collection::GeoJSON.GeometryCollection) = inside(point.geometry, collection)
inside(point::GeoJSON.Feature, collection::GeoJSON.FeatureCollection) = inside(point.geometry, collection)

for geojson_type in (GeoJSON.Point, GeoJSON.MultiPoint, GeoJSON.Feature)
    @eval begin
        within(point::$geojson_type, polygon::GeoJSON.Polygon) = inside(point, polygon) ? point : nothing
        within(point::$geojson_type, polygon::GeoJSON.MultiPolygon) = inside(point, polygon) ? point : nothing
        within(point::$geojson_type, polygon::GeoJSON.Feature) = inside(point, polygon) ? point : nothing
        within(point::$geojson_type, collection::GeoJSON.GeometryCollection) = inside(point, collection) ? point : nothing
        within(point::$geojson_type, collection::GeoJSON.FeatureCollection) = inside(point, collection) ? point : nothing
    end
end

# GeoJSON.GeometryCollection
within(collection::GeoJSON.GeometryCollection, polygon::GeoJSON.Polygon) = GeoJSON.GeometryCollection(Base.filter((geom,)->inside(geom,polygon), collection.geometries))
within(collection::GeoJSON.GeometryCollection, polygon::GeoJSON.MultiPolygon) = GeoJSON.GeometryCollection(Base.filter((geom,)->inside(geom,polygon), collection.geometries))
within(collection::GeoJSON.GeometryCollection, polygon::GeoJSON.Feature) = GeoJSON.GeometryCollection(Base.filter((geom,)->inside(geom,polygon), collection.geometries))
within(collection::GeoJSON.GeometryCollection, polygon::GeoJSON.GeometryCollection) = GeoJSON.GeometryCollection(Base.filter((geom,)->inside(geom,polygon), collection.geometries))
within(collection::GeoJSON.GeometryCollection, polygon::GeoJSON.FeatureCollection) = GeoJSON.GeometryCollection(Base.filter((geom,)->inside(geom,polygon), collection.geometries))

# GeoJSON.FeatureCollection
within(collection::GeoJSON.FeatureCollection, polygon::GeoJSON.Polygon) = GeoJSON.FeatureCollection(Base.filter((feature,)->inside(feature,polygon), collection.features))
within(collection::GeoJSON.FeatureCollection, polygon::GeoJSON.MultiPolygon) = GeoJSON.FeatureCollection(Base.filter((feature,)->inside(feature,polygon), collection.features))
within(collection::GeoJSON.FeatureCollection, polygon::GeoJSON.Feature) = GeoJSON.FeatureCollection(Base.filter((feature,)->inside(feature,polygon), collection.features))
within(collection::GeoJSON.FeatureCollection, polygon::GeoJSON.GeometryCollection) = GeoJSON.FeatureCollection(Base.filter((feature,)->inside(feature,polygon), collection.features))
within(collection::GeoJSON.FeatureCollection, polygon::GeoJSON.FeatureCollection) = GeoJSON.FeatureCollection(Base.filter((feature,)->inside(feature,polygon), collection.features))

# end

# function tag()
function tag(points::GeoJSON.FeatureCollection, polygons::GeoJSON.FeatureCollection, field::AbstractString, outfield::AbstractString)
    tagpoints = deepcopy(points)
    for pt in tagpoints.features
        if pt.properties == nothing
            pt.properties = Dict{AbstractString,Any}()
        end
        for poly in polygons.features
            if !haskey(pt.properties, outfield)
                if inside(pt, poly)
                    pt.properties[outfield] = poly.properties[field]
                end
            end
        end
    end
    tagpoints
end


text:
using JLD2
using PyPlot

function quick_and_dirty_plot(var, var_name)
    PyPlot.pcolormesh(var)
    PyPlot.title(var_name)
    PyPlot.colorbar()
    PyPlot.savefig("$var_name.png")
    PyPlot.close("all")
    @info "Saved $var_name.png"
end

face = 1 # just plot one face for now

jldopen("cubed_sphere_96_grid.jld2", "r") do file
    for key in keys(file["face$face"])
        quick_and_dirty_plot(file["face$face/$key"], key)
    end
end


text:
version = "6.5.1"
bfpkg_url = "https://downloads.openmicroscopy.org/bio-formats/$version/artifacts/bioformats_package.jar"
@info "Downloading version $version of bioformats_package.jar from $bfpkg_url"
download(bfpkg_url, joinpath(pwd(), "bioformats_package.jar"))

xsd_version = "2016-06"
xsd_url = "https://www.openmicroscopy.org/Schemas/OME/$xsd_version/ome.xsd"
@info "Downloading version $xsd_version of ome.xsd from $xsd_url"
download(xsd_url, joinpath(pwd(), "ome.xsd"))