text:
function get_build_number(pr_number::AbstractString;
                          organization_slug::AbstractString, pipeline_slug::AbstractString)
    _pr_number = strip(pr_number)
    !isempty(_pr_number)                  || throw(ArgumentError("You must provide a valid pull request number"))
    occursin(r"^[\d][\d]*?$", _pr_number) || throw(ArgumentError("You must provide a valid pull request number"))
    builds_per_page = 100
    max_pages = 100
    for page_number = 1:max_pages
        url = string(
            "https://api.buildkite.com/v2/",
            "organizations/$(organization_slug)/",
            "pipelines/$(pipeline_slug)/",
            "builds?per_page=$(builds_per_page)&page=$(page_number)",
        )
        cmd_string = string(
            "curl -X GET -H \"Authorization: Bearer $(ENV["BUILDKITE_API_TOKEN"])\" \"$(url)\"",
            "| ",
            "jq '[.[] | {pr: .pull_request.id, id: .id, number: .number}]' ",
            "| ",
            "jq '.[] | select(.pr == \"$(_pr_number)\") | .number'",
        )
        cmd = `bash -c $(cmd_string)`
        sleep(0.1) # helps us stay under the Buildkite API rate limits
        str = read(cmd, String)
        lines = strip.(strip.(strip.(split(strip(str), '\n')), Ref('"')))
        filter!(x -> !isempty(x), lines)
        if !isempty(lines)
            build_number = convert(String, strip(lines[1]))::String
            return build_number
        end # if
    end # for
    msg = string(
        "I tried $(max_pages) pages ",
        "(with $(builds_per_page) builds per page), ",
        "but I could not find any Buildkite builds for ",
        "GitHub pull request $(pr_number).",
    )
    throw(ErrorException(msg))
end

function get_failed_job_ids(build_number::AbstractString;
                            organization_slug::AbstractString, pipeline_slug::AbstractString)
    url = string(
        "https://api.buildkite.com/v2/",
        "organizations/$(organization_slug)/",
        "pipelines/$(pipeline_slug)/",
        "builds/$(build_number)",
    )
    cmd_string_1 = string(
        "curl -X GET -H \"Authorization: Bearer $(ENV["BUILDKITE_API_TOKEN"])\" \"$(url)\"",
        "| ",
        "jq '[.jobs | . | .[] | {id: .id, state: .state}]' ",
        "| ",
        "jq '[.[] | select((.state == \"failed\") or (.state == \"errored\"))]'",
    )
    cmd_string_2 = string(cmd_string_1, "| ", "jq '.[] .id'")
    cmd_1 = `bash -c $(cmd_string_1)`
    cmd_2 = `bash -c $(cmd_string_2)`
    sleep(0.1) # helps us stay under the Buildkite API rate limits
    run(cmd_1) # print the output to the log, for debugging purposes
    sleep(0.1) # helps us stay under the Buildkite API rate limits
    str = read(cmd_2, String)
    lines = strip.(strip.(strip.(split(strip(str), '\n')), Ref('"')))
    filter!(x -> !isempty(x), lines)
    failed_job_ids = convert(Vector{String}, lines)
    return failed_job_ids
end

function retry_job(build_number::AbstractString, job_id::AbstractString;
                   organization_slug::AbstractString, pipeline_slug::AbstractString)
    url = string(
        "https://api.buildkite.com/v2/",
        "organizations/$(organization_slug)/",
        "pipelines/$(pipeline_slug)/",
        "builds/$(build_number)/",
        "jobs/$(job_id)/retry",
    )
    cmd_string = string(
        "curl -X PUT -H \"Authorization: Bearer $(ENV["BUILDKITE_API_TOKEN"])\" \"$(url)\"",
    )
    cmd = `bash -c $(cmd_string)`
    sleep(0.1) # helps us stay under the Buildkite API rate limits
    run(cmd)
    return nothing
end

function retry_jobs(build_number::AbstractString, job_ids::AbstractVector{<:AbstractString};
                    organization_slug::AbstractString, pipeline_slug::AbstractString,)
    if isempty(job_ids)
        @info "There are no jobs to retry."
    end
    num_jobs = length(job_ids)
    for (i, job_id) in enumerate(job_ids)
        @info "$(i) of $(num_jobs). Attempting to retry job: $(job_id)"
        retry_job(build_number, job_id; organization_slug, pipeline_slug)
    end
    return nothing
end

function main(pr_number::AbstractString;
              organization_slug::AbstractString, pipeline_slug::AbstractString)
    @info "The Buildkite organization is $(organization_slug)"
    @info "The Buildkite pipeline is $(pipeline_slug)"
    @info "The GitHub pull request number is $(pr_number)"
    build_number = get_build_number(pr_number; organization_slug, pipeline_slug)
    @info "The build number is $(build_number)"
    failed_job_ids = get_failed_job_ids(build_number; organization_slug, pipeline_slug)
    @info "There are $(length(failed_job_ids)) failed jobs." failed_job_ids
    retry_jobs(build_number, failed_job_ids; organization_slug, pipeline_slug)
    return nothing
end

const pr_number         = ENV["PR_NUMBER"]
const organization_slug = ENV["BUILDKITE_ORGANIZATION_SLUG"]
const pipeline_slug     = ENV["BUILDKITE_PIPELINE_SLUG"]

main(pr_number; organization_slug, pipeline_slug)

max_stars_repo_name: JuliaLang/retry-buildkite

max_issues_repo_name: JuliaLang/retry-buildkite

max_forks_repo_name: JuliaLang/retry-buildkite


text:
function _bench_replica(
    num_samples,
    algorithm_obj::ReplicaAlgorithm{1}
)
    # num_samples = 2_048
    num_ants = 1
    num_correlators = 3
    enable_gpu = Val(true)
    system = GPSL1(use_gpu = Val(true));
    # system_h = GPSL1(use_gpu = Val(false));
    codes = system.codes
    # codes_text_mem_simple = CuTexture(
        # CuTextureArray(codes)
    # )
    # codes_text_mem = CuTexture(
    #     CuTextureArray(codes),
    #     address_mode = CUDA.ADDRESS_MODE_WRAP,
    #     interpolation = CUDA.NearestNeighbour(),
    #     normalized_coordinates = true
    # )
    code_frequency = get_code_frequency(system)
    code_length = get_code_length(system)
    start_code_phase = 0.0f0
    carrier_phase = 0.0f0
    carrier_frequency = 1500Hz
    prn = 1

    # Generate the signal;
    _, sampling_frequency = gen_signal(system, prn, carrier_frequency, num_samples, num_ants = NumAnts(num_ants), start_code_phase = start_code_phase, start_carrier_phase = carrier_phase)

    # Generate correlator;
    correlator = EarlyPromptLateCorrelator(NumAnts(num_ants), NumAccumulators(num_correlators))
    correlator_sample_shifts = get_correlator_sample_shifts(system, correlator, sampling_frequency, 0.5)
    num_of_shifts = correlator_sample_shifts[end] - correlator_sample_shifts[1]

    # Generate blank code and carrier replica, and downconverted signal;
    # code_replica_cpu = zeros(Float32, num_samples + num_of_shifts)
    code_replica = CUDA.zeros(Float32, num_samples + num_of_shifts)

    # Generate CUDA kernel tuning parameters;
    threads_per_block = 768
    blocks_per_grid = cld.(num_samples, threads_per_block)

    # @cuda threads=threads_per_block blocks=blocks_per_grid gen_code_replica_kernel!(
    #     code_replica,
    #     codes_text_mem, # texture memory codes
    #     code_frequency,
    #     sampling_frequency,
    #     start_code_phase,
    #     prn,
    #     num_samples + num_of_shifts,
    #     correlator_sample_shifts[1],
    #     code_length
    # )
    @benchmark CUDA.@sync @cuda threads=$threads_per_block blocks=$blocks_per_grid $gen_code_replica_kernel!(
        $code_replica,
        $codes, # texture memory codes
        $code_frequency,
        $sampling_frequency,
        $start_code_phase,
        $prn,
        $num_samples + $num_of_shifts,
        $correlator_sample_shifts[1],
        $code_length
    )
end

function _bench_replica(
    num_samples,
    algorithm_obj::ReplicaAlgorithm{2}
)
    # num_samples = 2_048
    num_ants = 1
    num_correlators = 3
    enable_gpu = Val(true)
    system = GPSL1(use_gpu = Val(true));
    # system_h = GPSL1(use_gpu = Val(false));
    codes = system.codes
    # codes_text_mem_simple = CuTexture(
        # CuTextureArray(codes)
    # )
    codes_text_mem = CuTexture(
        CuTextureArray(codes),
        address_mode = CUDA.ADDRESS_MODE_WRAP,
        interpolation = CUDA.NearestNeighbour(),
        normalized_coordinates = true
    )
    code_frequency = get_code_frequency(system)
    code_length = get_code_length(system)
    start_code_phase = 0.0f0
    carrier_phase = 0.0f0
    carrier_frequency = 1500Hz
    prn = 1

    # Generate the signal;
    _, sampling_frequency = gen_signal(system, prn, carrier_frequency, num_samples, num_ants = NumAnts(num_ants), start_code_phase = start_code_phase, start_carrier_phase = carrier_phase)

    # Generate correlator;
    correlator = EarlyPromptLateCorrelator(NumAnts(num_ants), NumAccumulators(num_correlators))
    correlator_sample_shifts = get_correlator_sample_shifts(system, correlator, sampling_frequency, 0.5)
    num_of_shifts = correlator_sample_shifts[end] - correlator_sample_shifts[1]

    # Generate blank code and carrier replica, and downconverted signal;
    # code_replica_cpu = zeros(Float32, num_samples + num_of_shifts)
    code_replica = CUDA.zeros(Float32, num_samples + num_of_shifts)

    # Generate CUDA kernel tuning parameters;
    threads_per_block = 768
    blocks_per_grid = cld.(num_samples, threads_per_block)

    # @cuda threads=threads_per_block blocks=blocks_per_grid gen_code_replica_texture_mem_kernel!(
    #     code_replica,
    #     codes_text_mem, # texture memory codes
    #     code_frequency,
    #     sampling_frequency,
    #     start_code_phase,
    #     prn,
    #     num_samples + num_of_shifts,
    #     correlator_sample_shifts[1],
    #     code_length
    # )
    @benchmark CUDA.@sync @cuda threads=$threads_per_block blocks=$blocks_per_grid $gen_code_replica_texture_mem_kernel!(
        $code_replica,
        $codes_text_mem, # texture memory codes
        $code_frequency,
        $sampling_frequency,
        $start_code_phase,
        $prn,
        $num_samples + $num_of_shifts,
        $correlator_sample_shifts[1],
        $code_length
    )
end

function run_replica_benchmark(benchmark_params::Dict)
    @unpack num_samples, algorithm = benchmark_params
    algorithm_obj = MEMDICT[algorithm]
    benchmark_results = _bench_replica(num_samples, algorithm_obj)
    benchmark_results_w_params = copy(benchmark_params)
    benchmark_results_w_params["Minimum"] = minimum(benchmark_results).time
    benchmark_results_w_params["Mean"] = mean(benchmark_results).time
    benchmark_results_w_params["Median"] = median(benchmark_results).time
    benchmark_results_w_params["Std"] = std(benchmark_results).time
    return benchmark_results_w_params
end

ext: jl

max_stars_repo_name: mfkiwl/GPUAcceleratedTracking

max_issues_repo_name: mfkiwl/GPUAcceleratedTracking

max_forks_repo_name: mfkiwl/GPUAcceleratedTracking


text:
function pancakesort!(arr::Vector{<:Real})
    len = length(arr)
    if len < 2 return arr end
    for i in len:-1:2
        j = indmax(arr[1:i])
        if i == j continue end
        arr[1:j] = reverse(arr[1:j])
        arr[1:i] = reverse(arr[1:i])
    end
    return arr
end

v = rand(-10:10, 10)
println("# unordered: $v\n -> ordered: ", pancakesort!(v))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Copyright (c) 2022 Code Komali
#
# This software is released under the MIT License.
# https://opensource.org/licenses/MIT
using Test
include("StudentData.jl")
using .StudentData

function swap!(arr,i,j)
    temp = arr[i]
    arr[i] = arr[j]
    arr[j] = temp
end

function partition(arr;
    startindex=1,
    endindex=lastindex(arr),
    desc = true,
    compareFn = desc ? (x,y)->x>y : (x,y) -> x<y)
    i = startindex + 1
    j = i
    pivot = arr[startindex]
    while j <= endindex
        if compareFn(arr[j], pivot)
            swap!(arr,i,j)
            i+=1
        end
        j+=1
    end
    swap!(arr,startindex,i-1)
    return i-1
end

function quicksort(arr;
    startindex=1,
    endindex=lastindex(arr),
    desc = true,
    compareFn = desc ? (x,y)->x>y : (x,y) -> x<y)
    startindex >= endindex && return arr
    j = partition(arr,startindex=startindex,endindex=endindex, desc=desc, compareFn = compareFn)
    quicksort(arr, startindex=startindex, endindex = j-1, desc=desc, compareFn = compareFn)
    quicksort(arr, startindex=j+1, endindex=endindex, desc=desc, compareFn = compareFn)
    return arr
end

# tests for quicksort
@testset "Quicksort tests" begin
    @test quicksort([11, 9, 7, 10, 8, 6]) == [11, 10, 9, 8, 7, 6]
    @test quicksort([11, 9, 7, 10, 8, 6], desc = false) == [6, 7, 8, 9, 10, 11]
    randArray = rand(Int8, 10)
    @test quicksort(randArray) == sort(randArray, rev = true) # Note: our quicksort implementation by default prefers descending order
    @test quicksort(randArray, desc = false) == sort(randArray)
    @test quicksort(students, compareFn = studentCompareFn) == sortedStudents
end

max_stars_repo_name: codekomali/JuliaPractice

max_issues_repo_name: codekomali/JuliaPractice

max_forks_repo_name: codekomali/JuliaPractice


text:
using JuliaWebAPI
using Test

include("srvr.jl")
include("clnt.jl")

function run_httprpctests()
    println("starting httpserver in async mode...")
    run_httprpcsrvr(JuliaWebAPI.JSONMsgFormat(), JuliaWebAPI.ZMQTransport(SRVR_ADDR, ZMQ.REP, true), true)
    wait_for_httpsrvr()
    println("starting client...")
    run_httpclnt()
end

!isempty(ARGS) && (ARGS[1] == "--runhttptests") && run_httprpctests()


text:
# interfacing with other packages

using LinearAlgebra
using LinearAlgebra: BlasFloat

Base.:(\)(A::Union{UpperTriangular{T, S},LowerTriangular{T, S}}, B::CuMatrix{T}) where {T<:BlasFloat, S<:AbstractCuSparseMatrix{T}}       = sm('N',A,B,'O')
Base.:(\)(transA::Transpose{T, UpperTriangular{T, S}}, B::CuMatrix{T}) where {T<:BlasFloat, S<:AbstractCuSparseMatrix{T}} = sm('T',parent(transA),B,'O')
Base.:(\)(transA::Transpose{T, LowerTriangular{T, S}}, B::CuMatrix{T}) where {T<:BlasFloat, S<:AbstractCuSparseMatrix{T}} = sm('T',parent(transA),B,'O')
Base.:(\)(adjA::Adjoint{T, UpperTriangular{T, S}},B::CuMatrix{T}) where {T<:BlasFloat, S<:AbstractCuSparseMatrix{T}} = sm('C',parent(adjA),B,'O')
Base.:(\)(adjA::Adjoint{T, LowerTriangular{T, S}},B::CuMatrix{T}) where {T<:BlasFloat, S<:AbstractCuSparseMatrix{T}} = sm('C',parent(adjA),B,'O')

LinearAlgebra.mul!(C::CuVector{T},A::CuSparseMatrix,B::CuVector) where {T} = mv!('N',one(T),A,B,zero(T),C,'O')
LinearAlgebra.mul!(C::CuVector{T},transA::Transpose{<:Any,<:CuSparseMatrix},B::CuVector) where {T} = mv!('T',one(T),parent(transA),B,zero(T),C,'O')
LinearAlgebra.mul!(C::CuVector{T},adjA::Adjoint{<:Any,<:CuSparseMatrix},B::CuVector) where {T} = mv!('C',one(T),parent(adjA),B,zero(T),C,'O')
LinearAlgebra.mul!(C::CuVector{T},A::HermOrSym{T,<:CuSparseMatrix{T}},B::CuVector{T}) where T = mv!('N',one(T),A,B,zero(T),C,'O')
LinearAlgebra.mul!(C::CuVector{T},transA::Transpose{<:Any, <:HermOrSym{T,<:CuSparseMatrix{T}}},B::CuVector{T}) where {T} = mv!('T',one(T),parent(transA),B,zero(T),C,'O')
LinearAlgebra.mul!(C::CuVector{T},adjA::Adjoint{<:Any, <:HermOrSym{T,<:CuSparseMatrix{T}}},B::CuVector{T}) where {T} = mv!('C',one(T),parent(adjA),B,zero(T),C,'O')

LinearAlgebra.mul!(C::CuMatrix{T},A::CuSparseMatrix{T},B::CuMatrix{T}) where {T} = mm2!('N','N',one(T),A,B,zero(T),C,'O')
LinearAlgebra.mul!(C::CuMatrix{T},A::CuSparseMatrix{T},transB::Transpose{<:Any, <:CuMatrix{T}})  where {T} = mm2!('N','T',one(T),A,parent(transB),zero(T),C,'O')
LinearAlgebra.mul!(C::CuMatrix{T},transA::Transpose{<:Any, <:CuSparseMatrix{T}},B::CuMatrix{T})  where {T} = mm2!('T','N',one(T),parent(transA),B,zero(T),C,'O')
LinearAlgebra.mul!(C::CuMatrix{T},transA::Transpose{<:Any, <:CuSparseMatrix{T}},transB::Transpose{<:Any, <:CuMatrix{T}}) where {T} = mm2!('T','T',one(T),parent(transA),parent(transB),zero(T),C,'O')
LinearAlgebra.mul!(C::CuMatrix{T},adjA::Adjoint{<:Any, <:CuSparseMatrix{T}},B::CuMatrix{T})  where {T} = mm2!('C','N',one(T),parent(adjA),B,zero(T),C,'O')

LinearAlgebra.mul!(C::CuMatrix{T},A::HermOrSym{<:Number, <:CuSparseMatrix},B::CuMatrix) where {T} = mm!('N',one(T),A,B,zero(T),C,'O')
LinearAlgebra.mul!(C::CuMatrix{T},transA::Transpose{<:Any, <:HermOrSym{<:Number, <:CuSparseMatrix}},B::CuMatrix) where {T} = mm!('T',one(T),parent(transA),B,zero(T),C,'O')
LinearAlgebra.mul!(C::CuMatrix{T},adjA::Adjoint{<:Any, <:HermOrSym{<:Number, <:CuSparseMatrix}},B::CuMatrix) where {T} = mm!('C',one(T),parent(adjA),B,zero(T),C,'O')

Base.:(\)(A::Union{UpperTriangular{T, S},LowerTriangular{T, S}}, B::CuVector{T}) where {T<:BlasFloat, S<:AbstractCuSparseMatrix{T}}       = sv2('N',A,B,'O')
Base.:(\)(transA::Transpose{T, UpperTriangular{T, S}},B::CuVector{T}) where {T<:BlasFloat, S<:AbstractCuSparseMatrix{T}} = sv2('T',parent(transA),B,'O')
Base.:(\)(transA::Transpose{T, LowerTriangular{T, S}},B::CuVector{T}) where {T<:BlasFloat, S<:AbstractCuSparseMatrix{T}} = sv2('T',parent(transA),B,'O')
Base.:(\)(adjA::Adjoint{T, UpperTriangular{T, S}},B::CuVector{T}) where {T<:BlasFloat, S<:AbstractCuSparseMatrix{T}}  = sv2('C',parent(adjA),B,'O')
Base.:(\)(adjA::Adjoint{T, LowerTriangular{T, S}},B::CuVector{T}) where {T<:BlasFloat, S<:AbstractCuSparseMatrix{T}}  = sv2('C',parent(adjA),B,'O')
Base.:(\)(A::AbstractTriangular{T,CuSparseMatrixHYB{T}},B::CuVector{T})       where T = sv('N',A,B,'O')
Base.:(\)(transA::Transpose{T, UpperTriangular{T, CuSparseMatrixHYB{T}}},B::CuVector{T}) where {T<:BlasFloat} = sv('T',parent(transA),B,'O')
Base.:(\)(transA::Transpose{T, LowerTriangular{T, CuSparseMatrixHYB{T}}},B::CuVector{T}) where {T<:BlasFloat} = sv('T',parent(transA),B,'O')
Base.:(\)(adjA::Adjoint{T, UpperTriangular{T, CuSparseMatrixHYB{T}}},B::CuVector{T}) where {T<:BlasFloat} = sv('C',parent(adjA),B,'O')
Base.:(\)(adjA::Adjoint{T, LowerTriangular{T, CuSparseMatrixHYB{T}}},B::CuVector{T}) where {T<:BlasFloat} = sv('C',parent(adjA),B,'O')

Base.:(+)(A::Union{CuSparseMatrixCSR,CuSparseMatrixCSC},B::Union{CuSparseMatrixCSR,CuSparseMatrixCSC}) = geam(A,B,'O','O','O')
Base.:(-)(A::Union{CuSparseMatrixCSR,CuSparseMatrixCSC},B::Union{CuSparseMatrixCSR,CuSparseMatrixCSC}) = geam(A,-one(eltype(A)),B,'O','O','O')


text:
function crash(mdp::PedCarMDP, s::PedCarMDPState)
    return crash(mdp, s.ego, s.ped, s.car)
end

function crash(mdp::PedCarMDP, ego::VehicleState, ped::VehicleState, car::VehicleState)
    return collision_checker(ego, car, mdp.ego_type, mdp.car_type) || collision_checker(ego, ped, mdp.ego_type, mdp.ped_type)
end

function POMDPs.discount(mdp::PedCarMDP)
    return mdp.γ
end

function POMDPs.isterminal(mdp::PedCarMDP, s::PedCarMDPState)
    if s.crash
        return true
    elseif s.ego.posF.s >= get_end(mdp.env.roadway[mdp.ego_goal]) &&
       get_lane(mdp.env.roadway, s.ego).tag == mdp.ego_goal
       return true
   end
   return false
end

function POMDPs.convert_s(::Type{V}, s::PedCarMDPState, mdp::PedCarMDP) where V<:AbstractArray
    n_routes = 4
    n_features = 4
    z = zeros(n_features*3 + n_routes + 1)
    z[1] = s.ego.posG.x / abs(mdp.env.params.x_min)
    z[2] = s.ego.posG.y / abs(mdp.env.params.y_min)
    z[3] = s.ego.posG.θ / π
    z[4] = s.ego.v / abs(mdp.env.params.speed_limit)
    z[5] = s.ped.posG.x / abs(mdp.env.params.x_min)
    z[6] = s.ped.posG.y / abs(mdp.env.params.y_min)
    z[7] = s.ped.posG.θ / π
    z[8] = s.ped.v / abs(mdp.env.params.speed_limit)
    z[9] = s.car.posG.x / abs(mdp.env.params.x_min)
    z[10] = s.car.posG.y / abs(mdp.env.params.y_min)
    z[11] = s.car.posG.θ / π
    z[12] = s.car.v / abs(mdp.env.params.speed_limit)
    # one hot encoding for the route
    routes = get_car_routes(mdp.env)
    for (i, r) in enumerate(routes)
        if r[1] == s.route[1] && r[end] == s.route[end]
            z[12+i] = 1.
        end
    end
    z[17] = float(s.crash)
    return z
end

function POMDPs.convert_s(::Type{PedCarMDPState}, z::V, mdp::PedCarMDP) where V<:AbstractArray{Float64}
    n_routes = 4
    n_features = 4
    @assert length(z) == n_features*3 + n_routes + 1
    ego_x = z[1]*abs(mdp.env.params.x_min)
    ego_y = z[2]*abs(mdp.env.params.y_min)
    ego_θ = z[3]*π
    ego_v = z[4]*abs(mdp.env.params.speed_limit)
    ego = VehicleState(VecSE2(ego_x, ego_y, ego_θ), mdp.env.roadway, ego_v)
    ped_x = z[5]*abs(mdp.env.params.x_min)
    ped_y = z[6]*abs(mdp.env.params.y_min)
    ped_θ = z[7]*π
    ped_v = z[8]*abs(mdp.env.params.speed_limit)
    ped = VehicleState(VecSE2(ped_x, ped_y, ped_θ), mdp.env.ped_roadway, ped_v)
    car_x = z[9]*abs(mdp.env.params.x_min)
    car_y = z[10]*abs(mdp.env.params.y_min)
    car_θ = z[11]*π
    car_v = z[12]*abs(mdp.env.params.speed_limit)
    car = VehicleState(VecSE2(car_x, car_y, car_θ), mdp.env.roadway, car_v)
    # one hot encoding for the route
    routes = get_car_routes(mdp.env)
    route = SVector{2, LaneTag}(LaneTag(0, 0), LaneTag(0, 0))
    for (i, r) in enumerate(routes)
        if z[12+i] == 1.
            route = SVector{2, LaneTag}(r[1], r[end])
        end
    end
    collision = Bool(z[17])
    return PedCarMDPState(collision, ego, ped, car, route)
end


text:
using CompressedSensing
using PyPlot
#This is a recreation of figure 1 from http://arxiv.org/pdf/1203.1548.pdf
#WARNING: this simulation takes a few minutes to complete!

#define settings used
sources=200
sensors=50
measurements= 10

noise=0
srand(0)

#how many monte carlo runs to do
max_monte=100

#at what sparsity should the algorithm be tested
S=10:5:40

mean_cor=zeros(length(S))

#For each sparsity level
for sparsity = 1:length(S)
    #run max_monte simulations and record the results
    for monte=1:max_monte
        #Create a signal that is sparsity sparse
        signal = zeros(sources,measurements)
        index=[] #Find "sparsity" random indecies
        while length(index)<S[sparsity]
            index = [index,unique(int(rand(1)*(sources-1)+1))]
        end
        #Set those random indecies to random gaussian values
        signal[index,:]=randn(S[sparsity],measurements)

        #Create an IID random measurement matrix
        measurement_Matrix = randn(sensors,sources)

        #"Measure" the signal, and add noise if desired
        measured = measurement_Matrix*signal+randn(sensors,measurements)*noise;

        #Attempt to rebuild the signal from the measurement and the measurement_Matrix
        re_signal = ZAP(measured,measurement_Matrix,verbose=false)

        #Record how often the answer is correct (we define correct to be >.99 correlation)
        mean_cor[sparsity]+=mean(sum(cor(signal, (re_signal)).*eye(measurements),1).>.99)./max_monte
    end
end

#plot the results
figure(figsize=(3.5,3))
xlabel("K Sparsity")
ylabel("Recovery Rate (%)")
ylim(0,1)
plot(S,mean_cor)
savefig("Example_5_Fig_1.png", bbox_inches="tight",pad_inches=.2)

max_stars_repo_path: examples/Example 5.jl

max_stars_repo_name: dahlend/CompressedSensing

max_issues_repo_path: examples/Example 5.jl

max_forks_repo_path: examples/Example 5.jl


text:
using Documenter, TopOpt
using DocumenterCitations

# Load packages to avoid precompilation output in the docs
# import ...

# Generate examples
include("generate.jl")

GENERATED_EXAMPLES = [
    joinpath("examples", f) for
    f in ("simp.md", "beso.md", "geso.md", "csimp.md", "global_stress.md")
]

bib = CitationBibliography(joinpath(@__DIR__, "biblio", "ref.bib"))
makedocs(
    bib;
    sitename="TopOpt.jl",
    format=Documenter.HTML(; prettyurls=get(ENV, "CI", nothing) == "true"),
    # doctest = false,
    pages=[
        "Home" => "index.md",
        "Examples" => GENERATED_EXAMPLES,
        "API Reference" => ["reference/TopOptProblems.md", "reference/Algorithms.md"],
        "Bibliography" => "bibliography.md",
    ],
)

# # make sure there are no *.vtu files left around from the build
# cd(joinpath(@__DIR__, "build", "examples")) do
#     foreach(file -> endswith(file, ".vtu") && rm(file), readdir())
# end

if get(ENV, "CI", nothing) == "true"
    deploydocs(; repo="github.com/JuliaTopOpt/TopOpt.jl.git", push_preview=true)
end


text:
# HIRES (chemical reaction, mildly stiff)
#
# http://www.dm.uniba.it/~testset/problems/hires.php
# http://www.unige.ch/~hairer/testset/testset.html
#
## Copyright (C) 2014, Mauro Werder <mauro_lc@runbox.com>
#
# Adapted from OdePkg, original author:
## Copyright (C) 2007-2012, Thomas Treichl <treichl@users.sourceforge.net>
## OdePkg - A package for solving ordinary differential equations and more
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; If not, see <http://www.gnu.org/licenses/>.

export hires

hires = let
    tcname = :hires
    T = Float64
    Tarr = Matrix
    dae = 0
    dof = 8

    # Computes the results for the HIRES problem
    function fn!(t, y, out)
        out[1] = -1.71 * y[1] + 0.43 * y[2] + 8.32 * y[3] + 0.0007
        out[2] =  1.71 * y[1] - 8.75 * y[2]
        out[3] = -10.03 * y[3] + 0.43 * y[4] + 0.035 * y[5]
        out[4] =  8.32 * y[2] + 1.71 * y[3] - 1.12 * y[4]
        out[5] = -1.745 * y[5] + 0.43 * (y[6] + y[7])
        out[6] = -280 * y[6] * y[8] + 0.69 * y[4] + 1.71 * y[5] - 0.43 * y[6] + 0.69 * y[7]
        out[7] =  280 * y[6] * y[8] - 1.81 * y[7]
        out[8] = -out[7]
        return nothing
    end
    # initializes storage for y:
    fn!(;T_::Type=T, dof_=dof) = zeros(T_,dof_)

    # Computes the JACOBIAN matrix for the HIRES problem
    function jac!(t, y, dfdy)
        dfdy[1,1] = -1.71
        dfdy[1,2] =  0.43
        dfdy[1,3] =  8.32
        dfdy[2,1] =  1.71
        dfdy[2,2] = -8.75
        dfdy[3,3] = -10.03
        dfdy[3,4] =  0.43
        dfdy[3,5] =  0.035
        dfdy[4,2] =  8.32
        dfdy[4,3] =  1.71
        dfdy[4,4] = -1.12
        dfdy[5,5] = -1.745
        dfdy[5,6] =  0.43
        dfdy[5,7] =  0.43
        dfdy[6,4] =  0.69
        dfdy[6,5] =  1.71
        dfdy[6,6] = -280 * y[8] - 0.43
        dfdy[6,7] =  0.69
        dfdy[6,8] = -280 * y[6]
        dfdy[7,6] =  280 * y[8]
        dfdy[7,7] = -1.81
        dfdy[7,8] =  280 * y[6]
        dfdy[8,6] = -280 * y[8]
        dfdy[8,7] =  1.81
        dfdy[8,8] = -280 * y[6]
        return nothing
    end
    # Returns a matrix which can hold the Jacobian of selected
    # type, initialized to zero.  Can also be used to make a
    # Matlab-style JPattern matrix with jac!(Bool)
    jac!(;T_::Type=T, dof_=dof) = zeros(T_,dof_,dof_)

    mass! = nothing

    ic = T[1, 0, 0, 0, 0, 0, 0, 0.0057]
    tspan = T[0.0, 321.8122]

    refsol = zeros(T, dof)
    refsol[1] = 0.73713125733256e-3
    refsol[2] = 0.14424857263161e-3
    refsol[3] = 0.58887297409675e-4
    refsol[4] = 0.11756513432831e-2
    refsol[5] = 0.23863561988313e-2
    refsol[6] = 0.62389682527427e-2
    refsol[7] = 0.28499983951857e-2
    refsol[8] = 0.28500016048142e-2

    refsolinds = trues(dof)

    scd_absinds = Int[]

    tc = TestCaseExplicit{tcname, T, Tarr}(
                             mildlystiff,
                             dae,
                             dof,
                             fn!,
                             jac!,
                             mass!,
                             ic,
                             tspan,
                             refsol,
                             refsolinds,
                             scd_absinds)

    tc_all[tcname] = tc
    tc_stiff[tcname] = tc
end


text:
# Script to test spglinecurvy

using GenSPGL
using MAT
fid = matopen("../compare/spgLine_vars.mat")

f = read(fid, "f")
dx = vec(read(fid, "dx"))
gtd = read(fid, "gtd")
lastFv = vec(read(fid, "lastFv"))
funForward = NLfunForward
funPenalty = funLS
b = vec(read(fid, "b"))
x = vec(read(fid, "x"))
feasSrchIt = convert(Int64, read(fid,"feasSrchIt"))
timeProject = zero(Float64)
linear_in = read(fid,"linear")

if linear_in == 0
    linear = false
else
    linear = true
end
A = funForward
opts = spgOptions(  optTol = 1e-5,
                    bpTol = 1e-5,
                    decTol = 1e-5,
                    project = TraceNorm_project,
                    primal_norm = TraceNorm_primal,
                    dual_norm = TraceNorm_dual,
                    proxy = true,
                    ignorePErr = true,
                    iterations = 150,
                    verbosity = 2)

# Avoid anon func
afunT(x) = reshape(x,355,709)
params = Dict{String, Any}([("nr", 40)
                                ("Ind", vec(b) .== 0)
                                ("numr", 355)
                                ("numc", 709)
                                ("funForward", NLfunForward)
                                ("afunT", afunT)
                                ("afun", afun)
                                ("mode", 1)
                                ("ls", 1)
                                ("logical", 0)
                                ("funPenalty", funLS)])

fNew, step, rNew, iter, err, localProdA = spgline(A,
        f,
        dx,
        gtd,
        x,
        maximum(lastFv),
        funForward,
        funPenalty,
        params,
        b,
        feasSrchIt,
        linear,
        opts,
        timeProject)


text:
using Soss

mstep = @model pars,state begin
    # Parameters
    α = pars.α        # Daily transmission rate
    β = pars.β        # Daily recovery rate
    γ = pars.γ        # Daily case fatality rate

    # Starting counts
    s0 = state.s      # Susceptible
    i0 = state.i      # Infected
    r0 = state.r      # Recovered
    d0 = state.d      # Deceased
    n = s0 + i0 + r0  # Population

    # Transitions between states
    si ~ Binomial(s0, α * i0 / n)
    ir ~ Binomial(i0, β)
    id ~ Binomial(i0 - ir, γ)

    # Updated counts
    return ( s = s0 - si
           , i = i0 + si - ir - id
           , r = r0 + ir
           , d = d0 + id
           )
end;

m = @model s0 begin
    α ~ Uniform()
    β ~ Uniform()
    γ ~ Uniform()
    pars = (α=α, β=β, γ=γ)
    x ~ Chain(s0) do s
            mstep(pars, s)
        end
end

using CSV
covid = CSV.read("/home/chad/git/covid-19/data/countries-aggregated.csv");

using DataFramesMeta

us = @where(covid, :Country .== "US");

# Some simple data cleaning - each column should be non-decreasing
# Fix decreasing values with geometric mean
for j in [:Confirmed, :Recovered, :Deaths]
    for i in 2:(size(us,1) - 1)
        if us[i, j] < us[i-1, j]
            us[i,j] = round(Int, sqrt(us[i-1,j] * us[i+1,j]))
        end
    end
end

using DataFrames
df = DataFrame(
    :i => us.Confirmed
    , :r => us.Recovered
    , :d => us.Deaths
    )
n = 331000000

df[!, :s] .= n .- df.i .- df.r .- df.d;
df[!, :si] .= 0;
df[:si][2:end] = .-diff(df.s);
df[!, :ir] .= 0;
df[:ir][2:end] = diff(df.r);
df[!, :id] .= 0;
df[:id][2:end] = diff(df.d);

# just the last 30 days
df = df[(end-29):end,:];

using NamedTupleTools
import NamedTupleTools
NamedTupleTools.namedtuple(dfrow::DataFrameRow) = namedtuple(pairs(dfrow).iter.is...)

nts = namedtuple(names(df)).(eachrow(df));

s0 = namedtuple(df[1,:]);
post = dynamicHMC(m(s0=s0),(x=namedtuple.(eachrow(df)),));
ppost=particles(post)

# R₀
ppost.α / (ppost.β + ppost.γ)

# Case fatality rate
ppost.γ / (ppost.β + ppost.γ)

# Implied infection duration
1/(ppost.β + ppost.γ)


text:
# Autogenerated wrapper script for SoapyUHD_jll for aarch64-linux-musl-cxx03
export libuhdSupport

using soapysdr_jll
using boost_jll
using USRPHardwareDriver_jll
JLLWrappers.@generate_wrapper_header("SoapyUHD")
JLLWrappers.@declare_library_product(libuhdSupport, "libuhdSupport.so")
function __init__()
    JLLWrappers.@generate_init_header(soapysdr_jll, boost_jll, USRPHardwareDriver_jll)
    JLLWrappers.@init_library_product(
        libuhdSupport,
        "lib/SoapySDR/modules0.8/libuhdSupport.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# NOTE: This file contains many parts that are copied from the file runtests.jl from the Package MPI.jl.
push!(LOAD_PATH, "../src") # FIXME: to be removed everywhere?

import ImplicitGlobalGrid # Precompile it.

excludedfiles = [ "test_excluded.jl"];

function runtests()
    exename   = joinpath(Sys.BINDIR, Base.julia_exename())
    testdir   = pwd()
    istest(f) = endswith(f, ".jl") && startswith(f, "test_")
    testfiles = sort(filter(istest, readdir(testdir)))

    nfail = 0
    printstyled("Testing package ImplicitGlobalGrid.jl\n"; bold=true, color=:white)
    for f in testfiles
        println("")
        if f ∈ excludedfiles
            println("Test Skip:")
            println("$f")
            continue
        end
        try
            run(`$exename -O3 --check-bounds=no $(joinpath(testdir, f))`)
        catch ex
            nfail += 1
        end
    end
    return nfail
end
exit(runtests())


text:
module Aref_LD_grad_function
export gradM
export final_T
export LG_Helper

using Roots
using DifferentialEquations
using Printf
include("leap_frog_definitions.jl")

barrier=5
condition_escape(u, t, integrator) =  maximum([abs(u[1]),abs(u[2]),abs(u[3]),abs(u[4])]) > barrier
affect_stop!(integrator) = terminate!(integrator)
callback_escape = DiscreteCallback(condition_escape, affect_stop!)
function condition_hits_PSS_q1(u, t, integrator) # Event when event_f(u,t) == 0
    u[1]
end
cb=ContinuousCallback(condition_hits_PSS_q1, affect_stop!, rootfind=true)

function gradM(mesh, H, t_end)
    delta=1e-7
    Q=mesh[2]
    P=mesh[1]

    mXb=[P-delta,Q]
    mXf=[P+delta,Q]
    mYb=[P,Q-delta]
    mYf=[P,Q+delta]

    ldXb=LD_Helper(mXb,H,t_end)
    ldXf=LD_Helper(mXf,H,t_end)
    # ldYb=LD_Helper(mYb,H,t_end)
    # ldYf=LD_Helper(mYf,H,t_end)
    dLDx=(ldXf-ldXb)/(2*delta)
    # dLDy=(ldYf-ldYb)/(2*delta)
    # M=(dLDx)^2+(dLDy)^2
    return maximum([log10(abs(dLDx)),4])

end

function LD_Helper(mesh, H,  t_end)
    Q=mesh[2]
    P=mesh[1]
     #for a given Q,P,H with X=0
    Y=Yfind_Aref(Q,P,H)
    if ~isempty(Y)
        u0=zeros(5)
        u0[1]=0 #X
        u0[2]=P #P
        u0[3]=Q #Q
        u0[4]=Y #Y
        u0[5]=0 # counter variable

        #constructor for ODE

        prob_b = ODEProblem(Eq_of_M_Lagrangian_Descriptors_Backwards,u0,(0., t_end))
        prob_f = ODEProblem(Eq_of_M_Lagrangian_Descriptors,u0,(0., t_end))

        sol_f=solve(prob_f, Tsit5(),maxiters=1e20,reltol=1e-10,abstol=1e-13,save_idxs = [5],save_every_step=false, save_end=true, dense=false)#,abstol=1e-9)
        # sol_b=solve(prob_b, Tsit5(),maxiters=1e20,reltol=1e-9,abstol=1e-12,save_idxs = [5],save_every_step=false, save_end=true, dense=false)#,abstol=1e-9)

        # sol_f=solve(prob_f,  Vern9(),maxiters=1e20,reltol=1e-14,abstol=1e-14,save_idxs = [5],save_every_step=false, save_end=true, dense=false)#,abstol=1e-9)
        # sol_b=solve(prob_b,  Vern9(),maxiters=1e20,reltol=1e-14,abstol=1e-14,save_idxs = [5],save_every_step=false, save_end=true, dense=false)#,abstol
        # sol_b=solve(prob_b, RK4(),maxiters=1e20,reltol=1e-6,abstol=1e-8)

        # sol_f=solve(prob_f, RK4(),maxiters=1e20,reltol=1e-8,abstol=1e-10)
        # sol_b=solve(prob_b, RK4(),maxiters=1e20,reltol=1e-8,abstol=1e-10)

        # sol_b=solve(prob_b, RK4(),maxiters=1e20)
        # sol_f=solve(prob_f, RK4(),maxiters=1e20)

        # sol_b=solve(prob_b, RK4(),maxiters=1e20,callback=callback_escape)
        # sol_f=solve(prob_f, RK4(),maxiters=1e20,callback=callback_escape)

        # sol_f=solve(prob_f, Vern9(),maxiters=1e20,reltol=1e-13,abstol=1e-15)
        # sol_b=solve(prob_b, Vern9(),maxiters=1e20,reltol=1e-13,abstol=1e-15)

        #find final value to check if energy was conserved
        # uf=zeros(5)
        # uf[1]=sol[1,end] #X
        # uf[2]=sol[2,end] #P
        # uf[3]=sol[3,end] #Q
        # uf[4]=sol[4,end] #Y
        # dH=abs(H_test(uf)-H)

        LD=sol_f[1,end]
        # LD=sol_b[1,end]
        # if maximum([abs(sol_f[1,end]),abs(sol_f[2,end]),abs(sol_f[3,)nd]),abs(sol_f[4,end]),abs(sol_b[1,end]),abs(sol_b[2,end]),abs(sol_b[3,end]),abs(sol_b[4,end])])<barrier
        return LD
    else
        return NaN
    end
end
function gradient_matrix(LD_Matrix,dx,dy)
    N,N=size(LD_Matrix)
    grad_X_Matrix=zeros(N-4,N-4)
    grad_Y_Matrix=zeros(N-4,N-4)
    for i=1:N-4
        for j=1:N-4
            grad_X_Matrix[i,j]=-LD_Matrix[i+4,j+1]+8*LD_Matrix[i+3,j+1]-8*LD_Matrix[i+1,j+1]+LD_Matrix[i,j+1]
            grad_Y_Matrix[i,j]=-LD_Matrix[i+1,j+4]+8*LD_Matrix[i+1,j+3]-8*LD_Matrix[i+1,j+1]+LD_Matrix[i+1,j]
        end
    end
    X= (grad_X_Matrix/(12*dx) ).^2
    Y= (grad_Y_Matrix/(12*dy) ).^2
    return X,Y
end

function final_T(H)
    Q=1e-8
    P=1e-8
    Y=Yfind_Aref(Q,P,H)
    if ~isempty(Y)
        u0=zeros(5)
        u0[1]=0 #X
        u0[2]=P #P
        u0[3]=Q #Q
        u0[4]=Y #Y
        u0[5]=0 # counter variable

        #constructor for ODE

        prob = ODEProblem(Eq_of_M_Lagrangian_Descriptors_Backwards,u0,(0., 20))
        sol=solve(prob, Vern9(),maxiters=1e20,reltol=1e-13,abstol=1e-15,callback=cb)

        return sol.t[end]
    end
end

function gradient_matrix_2(LD_Matrix,dx,dy)
    N,N=size(LD_Matrix)
    grad_X_Matrix=zeros(N-2,N-2)
    grad_Y_Matrix=zeros(N-2,N-2)
    for i=1:N-2
        for j=1:N-2
            grad_X_Matrix[i,j]=LD_Matrix[i+2,j+1]-LD_Matrix[i,j+1]
            grad_Y_Matrix[i,j]=LD_Matrix[i+1,j+2]-LD_Matrix[i+1,j]
        end
    end
    grad_M=( grad_X_Matrix/(2*dx) ).^2+( grad_Y_Matrix/(2*dy) ).^2
    return sqrt.(grad_M)
end

function gradient_matrix_4(LD_Matrix,dx,dy)
    N,N=size(LD_Matrix)
    grad_X_Matrix=zeros(N-4,N-4)
    grad_Y_Matrix=zeros(N-4,N-4)
    for i=1:N-4
        for j=1:N-4
            grad_X_Matrix[i,j]=-LD_Matrix[i+4,j+1]+8*LD_Matrix[i+3,j+1]-8*LD_Matrix[i+1,j+1]+LD_Matrix[i,j+1]
            grad_Y_Matrix[i,j]=-LD_Matrix[i+1,j+4]+8*LD_Matrix[i+1,j+3]-8*LD_Matrix[i+1,j+1]+LD_Matrix[i+1,j]
        end
    end
    grad_M=( grad_X_Matrix/(12*dx) ).^2+( grad_Y_Matrix/(12*dy) ).^2
    return sqrt.(grad_M)
end

function D_X(LD_Matrix,dx)
    N=size(LD_Matrix)
    N=N-6
    grad_X_Matrix=zeros(N,N)
    for i=1:N
        for j=1:N
            grad_X_Matrix[i,j]=-LD_Matrix[i+4,j+1]+8*LD_Matrix[i+3,j+1]-8*LD_Matrix[i+1,j+1]+LD_Matrix[i,j+1]
        end
    end
    grad_M=abs(grad_X_Matrix/(12*dx) )
    return grad_M
end

end


text:
using SparseArrays

function comp_merit(it::Class_iterate)
    return LinearAlgebra.norm(comp(it),Inf)^3 / get_mu(it)^2
end

#function comp_merit_predicted(it::Class_iterate, dir::Class_point, step_size::Float64)
#    return comp_predicted(it, dir, step_size)^3
#end

function comp(it::Class_iterate)
    return it.point.s .* it.point.y .- it.point.mu
end

function comp_ratio_max(it::Class_iterate)
    return max(maximum(it.point.s .* it.point.y ./ it.point.mu), maximum(it.point.mu ./ (it.point.s .* it.point.y)))
end

mutable struct Eval_NaN_error <: Exception
   num::Float64
   point::Array{Float64,1}
   place::String

   function Eval_NaN_error(num::Float64,point::Array{Float64,1},place::String)
     return new(num::Float64,point::Array{Float64,1},place::String)
   end
end

Base.showerror(io::IO, e::Eval_NaN_error) = print(io, "$(e.num) by $(e.place) not expected!");

function isbad(num::Float64)
    if isnan(num) || isinf(num)
      return true
      #my_warn("$num not expected!")
      #Cutest_NaN_error("$num not expected!", num, point, place)
    else
      return false
    end
end

function isbad(vec)
  for i = 1:length(vec)
    if isbad(vec[i])
      return true
    end
  end

  return false
end

function eval_r(it::Class_iterate)
    x = it.point.x
    a_norm_pen = it.a_norm_penalty_par * sum(get_cons(it))
    return a_norm_pen
end

function eval_grad_r(it::Class_iterate)
    x = it.point.x
    a_norm_grad = it.a_norm_penalty_par * eval_jac_T_prod(it, ones(length(get_cons(it))))
    return a_norm_grad
end

function get_fval(it::Class_iterate)
    return it.cache.fval
end

function get_grad(it::Class_iterate)
    return it.cache.grad #+
end

function get_cons(it::Class_iterate)
    return it.cache.cons
end

function get_jac(it::Class_iterate)
    return it.cache.J
end

function get_jac_T(it::Class_iterate)
    return it.cache.J_T
end

function eval_J_T_J(it::Class_iterate, diag_vals::Vector)
    return it.cache.J_T * sparse(Diagonal(diag_vals)) * it.cache.J
end

function eval_diag_J_T_J(it::Class_iterate,diag_vals::Vector)
    n = dim(it)
    di = zeros(n)
    for i = 1:n
        a = it.cache.J[:,i]
        for j in a.nzind
            di[i] += a[j]^2 * diag_vals[j]
        end
        #dot(a .* diag_vals,a)
    end
    return di
end

function eval_jac_prod(it::Class_iterate, x::Vector)
    return it.cache.J * x
end

function eval_jac_T_prod(it::Class_iterate, y::Vector)
    return it.cache.J_T * y
end

function get_primal_res(it::Class_iterate)
    return it.cache.cons .- it.point.s
end

function get_max_vio(it::Class_iterate)
    return -min(0.0,minimum(it.cache.cons))
end

function eval_phi(it::Class_iterate, mu::Float64)
    if all(it.point.s .> 0.0)
        return get_fval(it) - mu * sum( log.( it.point.s ) ) + mu * eval_r(it)
    else
        return Inf
    end
end

function eval_grad_phi(it::Class_iterate, mu::Float64)
    y_tilde = it.point.mu ./ it.point.s #
    return eval_grad_lag(it, mu, y_tilde)
end

function get_lag_hess(it::Class_iterate)
    return it.cache.H # remember this is a triangular matrix!!!!
end

function eval_grad_lag(it::Class_iterate, mu::Float64)
    return eval_grad_lag(it, mu, it.point.y)
end

function eval_grad_lag(it::Class_iterate, mu::Float64, y::Array{Float64,1})
    return get_grad(it) - eval_jac_T_prod(it, y) + mu * eval_grad_r(it)
end

function dynamic_eval_Jt_prod(it::Class_iterate)
    return eval_Jt_prod(it.nlp, it.point.x, it.point.y)
end

function eval_farkas(it::Class_iterate, y::Array{Float64,1})
    return max(farkas_certificate(it, y), stationary_infeasible_measure(it,y))
end

function farkas_certificate(it::Class_iterate, y::Array{Float64,1})
    c = it.cache.cons
    feas_obj = -sum(c .* y)
    if feas_obj > 0.0
      return LinearAlgebra.norm(it.cache.J_T * y, 1) / feas_obj
    else
      return Inf
    end
end

function stationary_infeasible_measure(it::Class_iterate, y::Array{Float64,1})
    return (LinearAlgebra.norm(it.cache.J_T * y, 1) + dot(it.point.s, it.point.y)) / LinearAlgebra.norm(y,1)
end

function eval_farkas(it::Class_iterate)
    return eval_farkas(it, it.point.y)
end

function eval_merit_function(it::Class_iterate, pars::Class_parameters)
    if is_feasible(it, pars.ls.comp_feas)
        if length(it.point.s) > 0
          comp_penalty = LinearAlgebra.norm(comp(it), Inf)^3 / (it.point.mu)^2
        else
          comp_penalty = 0.0
        end

        return eval_phi(it, it.point.mu) + comp_penalty
    else
        return Inf
    end
end

function eval_phi_diff(it::Class_iterate, candidate::Class_iterate, mu::Float64)
    fdiff = get_fval(candidate) - get_fval(it)
    rdiff = mu * eval_r(candidate) - mu * eval_r(it)
    logdiff = - mu * sum( log.( candidate.point.s ./ it.point.s) )

    return fdiff + rdiff + logdiff
end

function eval_merit_function_difference(it::Class_iterate, candidate::Class_iterate, pars::Class_parameters)
    if is_feasible(candidate, pars.ls.comp_feas)
        if length(it.point.s) > 0
          comp_penalty = LinearAlgebra.norm(comp(candidate), Inf)^3 / (it.point.mu)^2 - LinearAlgebra.norm(comp(it), Inf)^3 / (it.point.mu)^2
        else
          comp_penalty = 0.0
        end

        #phi_diff = eval_phi(candidate, candidate.point.mu) - eval_phi(it, it.point.mu)
        phi_diff = eval_phi_diff(it, candidate, candidate.point.mu)
        return phi_diff + comp_penalty
    else
        return Inf
    end
end

function is_lower_triangular(mat)
    for i = 1:size(mat,1)
        for j = (i+1):size(mat,2)
            if mat[i,j] != 0.0
                return false
            end
        end
    end
    return true
end

function vector_product(LowerTriangularMat::SparseMatrixCSC{Float64,Int64}, vector::Array{Float64,1})
    DEBUG_MODE = false
    if DEBUG_MODE
        warn("checking lower triangular property")
        @assert is_lower_triangular(LowerTriangularMat) == true
    end
    v1 = LowerTriangularMat * vector
    v2 = LowerTriangularMat' * vector
    return v1 + v2 - diag(LowerTriangularMat) .* vector
end

function hess_product(it::Class_iterate, vector::Array{Float64,1})
    return vector_product(it.cache.H, vector)
end

function phi_predicted_reduction_primal_dual(it::Class_iterate, dir::Class_point, step_size::Float64)
    H = it.cache.H
    J = it.cache.J

    v = J * dir.x
    J_gain = dot(v.^2, it.point.y ./ it.point.s)

    #@show H

    h_prod = vector_product(H, dir.x) # the hessian is lower triangular!

    return step_size * dot(dir.x, eval_grad_phi(it, it.point.mu)) + step_size^2 * 0.5 * (dot(dir.x, h_prod) + J_gain)
end

function comp_predicted(it::Class_iterate, dir::Class_point, step_size::Float64)
    y = it.point.y;
    s = it.point.s;
    return s .* y + dir.y .* s * step_size + dir.s .* y * step_size .- (it.point.mu + dir.mu * step_size)
end

function merit_function_predicted_reduction(it::Class_iterate, dir::Class_point, step_size::Float64)
    C_k = LinearAlgebra.norm(comp(it),Inf)
    P_k = LinearAlgebra.norm(comp_predicted(it, dir, step_size), Inf)
    #@show C_k, P_k
    if length(it.point.s) > 0
      comp_penalty = (P_k^3 - C_k^3) / (it.point.mu)^2
    else
      comp_penalty = 0.0
    end
    #@show comp_penalty

    phi_red = phi_predicted_reduction_primal_dual(it, dir, step_size)
    #@show predict_red

    #@show comp_penalty

    return phi_red + comp_penalty
end

function eval_kkt_err(it::Class_iterate, par::Class_parameters)
    return scaled_dual_feas(it, par) + LinearAlgebra.norm(comp(it), Inf)
end


text:
###############################################################################
#
# Root finding method
#
###############################################################################
"""
Abstract type of the ConstrainedRootSolvers methods

$(TYPEDEF)

"""
abstract type AbstractCRSMethod{FT<:AbstractFloat} end

"""
Bisection method for 1D root solvers

$(TYPEDEF)

# Fields

$(TYPEDFIELDS)

"""
Base.@kwdef mutable struct BisectionMethod{FT<:AbstractFloat} <:
                           AbstractCRSMethod{FT}
    "lower bound"
    x_min::FT = 0
    "upper bound"
    x_max::FT = 0
    "matrix that stores x and y"
    xy::Matrix{FT} = FT[x_min 0; (x_min+x_max)/2 0; x_max 0]

    # history Vector
    "history of all simulations"
    history::Vector = Vector{FT}[]
end

"""
Nelder-Mead method for 2D and above solvers

$(TYPEDEF)

# Fields

$(TYPEDFIELDS)

"""
Base.@kwdef mutable struct NelderMeadMethod{FT<:AbstractFloat} <:
                           AbstractCRSMethod{FT}
    "Number of parameters to optimize"
    N::Int = 2
    "Initial values"
    x_inis::Array{FT,1} = zeros(FT,N+1)
    "Simplex array of array with dimension (N+1) * (N+1)"
    simplex::Array{Array{FT,1},1} = [zeros(FT,N+1) for i in 1:(N+1)]

    # temporary containers
    "Centroid"
    cen_x::Array{FT,1} = deepcopy(x_inis)
    "Reflection"
    ref_x::Array{FT,1} = deepcopy(x_inis)
    "Expansion"
    exp_x::Array{FT,1} = deepcopy(x_inis)
    "Contraction"
    con_x::Array{FT,1} = deepcopy(x_inis)
end

"""
Newton's method constrained by mininum and maximum ranges for 1D root solver

$(TYPEDEF)

# Fields

$(TYPEDFIELDS)

"""
Base.@kwdef mutable struct NewtonBisectionMethod{FT<:AbstractFloat} <:
                           AbstractCRSMethod{FT}
    "Lower bound"
    x_min::FT = 0
    "Upper bound"
    x_max::FT = 1
    "Initial guess"
    x_ini::FT = (x_min + x_max) / 2

    # history Vector
    "history of all simulations"
    history::Vector = Vector{FT}[]
end

"""
Newton raphson method for 1D root solver

$(TYPEDEF)

# Fields

$(TYPEDFIELDS)

"""
Base.@kwdef mutable struct NewtonRaphsonMethod{FT<:AbstractFloat} <:
                           AbstractCRSMethod{FT}
    "Initial guess"
    x_ini::FT

    # history Vector
    "history of all simulations"
    history::Vector = Vector{FT}[]
end

"""
Reduce step method for 1D root solver. This method increases or decreases from
    initial guess until no improvement is found. Then the incremantal step
    decreases, and then the root solver continues.

$(TYPEDEF)

# Fields

$(TYPEDFIELDS)

"""
Base.@kwdef mutable struct ReduceStepMethod{FT<:AbstractFloat} <:
                           AbstractCRSMethod{FT}
    "Lower bound"
    x_min::FT = 0
    "Upper bound"
    x_max::FT = 1
    "Initial guess"
    x_ini::FT = (x_min + x_max) / 2
    "Initial step"
    Δ_ini::FT = (x_min + x_max) / 10

    # history Vector
    "history of all simulations"
    history::Vector = Vector{FT}[]
end

"""
Reduce step method for 2D and above root solver. This method increases or
    decreases each variable in the initial guess until no improvement is found.
    Then the incremental steps decreases, and then the root solver continues.

$(TYPEDEF)

# Fields

$(TYPEDFIELDS)

"""
Base.@kwdef mutable struct ReduceStepMethodND{FT<:AbstractFloat} <:
                           AbstractCRSMethod{FT}
    "Lower bound"
    x_mins::Array{FT,1} = zeros(FT,2)
    "Upper bound"
    x_maxs::Array{FT,1} =  ones(FT,2)
    "Initial guess"
    x_inis::Array{FT,1} = FT[0.5, 0.5]
    "Target x"
    x_targ::Array{FT,1} = deepcopy(x_inis)
    "Temporary x"
    x_temp::Array{FT,1} = deepcopy(x_inis)
    "Initial step"
    Δ_inis::Array{FT,1} = FT[0.1, 0.1]
    "Operation step"
    Δ_oper::Array{FT,1} = deepcopy(Δ_inis)
    "Array of judges"
    Δjd::Array{Bool,1} = [false for i in 1:length(x_inis)]
end


text:
@testset "OneDAtX" begin

  # On vérifie que le constructeur par défaut fonctionne
  ls_at_t = OneDAtX(0.0)

  @test scoretype(ls_at_t) == Float64
  @test xtype(ls_at_t) == Float64

  @test ls_at_t.x == 0.0
  @test isnan(ls_at_t.fx)
  @test isnan(ls_at_t.gx)
  @test isnan(ls_at_t.f₀)
  @test isnan(ls_at_t.g₀)

  @test isnan(ls_at_t.current_time)
  @test isnan(ls_at_t.current_score)

  # On test la fonction update!(...)
  update!(ls_at_t, x = 1.0, fx = 1.0, gx = 1.0, f₀ = 1.0)
  update!(ls_at_t, g₀ = 1.0, current_time = 0.0, current_score = 0.0)

  @test ls_at_t.x == 1.0
  @test ls_at_t.fx == 1.0
  @test ls_at_t.gx == 1.0
  @test ls_at_t.f₀ == 1.0
  @test ls_at_t.g₀ == 1.0
  @test ls_at_t.current_time == 0.0
  @test ls_at_t.current_score == 0.0

  # on vérifie que la fonction copy fonctionne
  ls_at_t_2 = copy(ls_at_t)

  @test scoretype(ls_at_t_2) == Float64
  @test xtype(ls_at_t_2) == Float64

  @test ls_at_t_2.x == 1.0
  @test ls_at_t_2.fx == 1.0
  @test ls_at_t_2.gx == 1.0
  @test ls_at_t_2.f₀ == 1.0
  @test ls_at_t_2.g₀ == 1.0
  @test ls_at_t_2.current_time == 0.0
  @test ls_at_t_2.current_score == 0.0

  ls_64 = OneDAtX(0.0)

  @test scoretype(ls_64) == Float64
  @test xtype(ls_64) == Float64

  update!(ls_64, x = 1.0, fx = 1.0, gx = 1.0, f₀ = 1.0)

  reinit!(ls_64)
  @test ls_64.x == 1.0
  @test isnan(ls_64.fx)
  @test isnan(ls_64.current_time)
end


text:
"""
	CustomNetwork <: AbstractNetwork

A custom network created from an adjacency matrix ou vector.

# Fields
- N              :: Integer -> Number of nodes
"""
struct CustomNetwork <: AbstractNetwork
	N :: Integer

	seed :: Integer
	_adjMat :: AbstractMatrix
	_props :: NetworkProperties
end

"""
	CustomNetwork(adjMatrx::AbstractMatrix; directed=nothing)

Create a network with `N` nodes from the adjacency matrix.
"""
function CustomNetwork(mat::AbstractMatrix; directed::Union{Bool,Nothing}=nothing)
	(size(mat, 1) != size(mat, 2)) &&
		throw(ArgumentError("The adjacency matrix must be a square matrix!"))

	try
		BitArray(mat)
	catch
		throw(ArgumentError("Weighted matrices are not supported"))
	end

	if isnothing(directed)
		directed = !issymmetric(mat)
	else
		if !directed && !issymmetric(mat)
			throw(ArgumentError(
				"The adjacency matrix must be symmetric for a undirected network"
			))
		end
	end

	CustomNetwork(
		size(mat, 1),
		rand(1:99999999),
		SparseArrays.issparse(mat) ? copy(mat) : BitArray(mat),
		NetworkProperties(size(mat, 1), directed)
	)
end

"""
	CustomNetwork(adjVet::Vector{<:Integer}, N::Integer; directed=nothing)

Create a network with `N` nodes from the adjacency vector.

A keywork argument `directed` can be passed to ensure a (un)directed network.
If `nothing` this property will be derived from the network's symmetry.
"""
function CustomNetwork(vet::Vector{<:Integer}, N::Integer; directed::Union{Bool,Nothing}=nothing)
	(N <= 0) && throw(ArgumentError("Number of nodes must be a positive integer!"))
	CustomNetwork(adjVetToMat(vet, N, sparse=true))
end

function show(network::CustomNetwork)
	print(isdirected(network) ? "Directed " : "Undirected ")
	println("Custom Network")
	println("- N = $(network.N)")
end


text:
module Base_threads

end


text:
import Base: first, last, length,
    getindex, lastindex, eachindex, iterate,
    eltype, similar,
    show

# * AbstractKnotSet
"""
    AbstractKnotSet{k,ml,mr,T}

Abstract base for B-spline knot sets. `T` is the `eltype` of the knot
set, `k` is the order of the piecewise polynomials (`order = degree +
1`) and `ml` and `mr` are the knot multiplicities of the left and
right endpoints, respectively.  """
abstract type AbstractKnotSet{k,ml,mr,T} <: AbstractVector{T} end

distribution(::AbstractKnotSet) = NonUniform()

"""
    assert_multiplicities(k,ml,mr,t)

Assert that the multiplicities at the endpoints, `ml` and `mr`,
respectively, are consistent with the order `k`. Also check that the
amount of knots in the knot set `t` are enough to support the
requested order `k`.
"""
function assert_multiplicities(k,ml,mr,t)
    k > 0 || throw(ArgumentError("Polynomial order must be a positive integer, got k = $k"))
    ml ∈ 1..k || throw(ArgumentError("Left multiplicity has to be in range 1..$(k), got $ml"))
    mr ∈ 1..k || throw(ArgumentError("Right multiplicity has to be in range 1..$(k), got $mr"))
    n = length(t)
    nreq = max(2, k + 3 - ml - mr)
    n ≥ nreq || throw(ArgumentError("Polynomial order k = $(k) and left,right multiplicities $(ml),$(mr) require a knot sequence of at least $(nreq) points"))
end

"""
    order(t)

Returns the order `k` of the knot set `t`.
"""
order(t::AbstractKnotSet{k}) where k = k

"""
    numinterval(t)

Returns the number of intervals generated by the knot set `t`.

# Examples

```jldoctest
julia> numintervals(LinearKnotSet(3, 0, 1, 2))
2
```
"""
numintervals(t::AbstractKnotSet) = length(t.t)-1

"""
    numfunctions(t)

Returns the number of basis functions generated by knot set `t`.

# Examples

```jldoctest
julia> numfunctions(LinearKnotSet(3, 0, 1, 2))
4

julia> numfunctions(LinearKnotSet(5, 0, 1, 2))
6
```
"""
numfunctions(t::AbstractKnotSet) = length(t) - order(t)

"""
    leftmultiplicity(t)

Return the multiplicity of the left endpoint.
"""
leftmultiplicity(::AbstractKnotSet{k,ml,mr}) where {k,ml,mr} = ml

"""
    rightmultiplicity(t)

Return the multiplicity of the right endpoint.
"""
rightmultiplicity(::AbstractKnotSet{k,ml,mr}) where {k,ml,mr} = mr

"""
    first(t)

Return the first knot of `t`.
"""
first(t::AbstractKnotSet) = first(t.t)

"""
    last(t)

Return the last knot of `t`.
"""
last(t::AbstractKnotSet) = last(t.t)

"""
    length(t)

Return the number of knots of `t`.

# Examples

```jldoctest
julia> length(LinearKnotSet(3, 0, 1, 3))
8

julia> length(LinearKnotSet(3, 0, 1, 3, 1, 1))
4
```
"""
length(t::AbstractKnotSet) = length(t.t) + leftmultiplicity(t) + rightmultiplicity(t) - 2
# size(t) is not very useful for knot sets, but necessary to fulfill
# the interface expected from AbstractVectors.
size(t::AbstractKnotSet) = (length(t),)

"""
    getindex(t, i)

Return the `i`th knot of the knot set `t`, accounting for the
multiplicities of the endpoints.

# Examples

```jldoctest
julia> LinearKnotSet(3, 0, 1, 3)[2]
0.0

julia> LinearKnotSet(3, 0, 1, 3, 1, 1)[2]
0.3333333333333333
```

"""
function getindex(t::AbstractKnotSet{k,ml,mr}, i::Integer) where {k,ml,mr}
    i < 1 && throw(BoundsError(t, i))

    ni = numintervals(t)
    if i < ml
        first(t)
    elseif i < ml + ni
        t.t[i-ml+1]
    elseif i < ml + ni + mr
        last(t)
    else
        throw(BoundsError("Trying to access knot set of length $(length(t)) at index $i"))
    end
end
lastindex(t::AbstractKnotSet) = length(t)
eachindex(t::AbstractKnotSet) = 1:length(t)
# function iterate(t::AbstractKnotSet, (element,i)=(t[1],0))
#     i ≥ length(t) && return nothing
#     element, (t[i+1], i+1)
# end

eltype(t::AbstractKnotSet{k,ml,mr,T}) where {k,ml,mr,T} = T
similar(t::AbstractKnotSet{k,ml,mr,T}) where {k,ml,mr,T} = Vector{T}(undef, length(t))

function show_order_multiplicities(io::IO, t::AbstractKnotSet{k,ml,mr}) where {k,ml,mr}
    write(io, "order k = $(k)")
    k < 7 && write(io, " (",
                  ["constant", "linear", "parabolic", "cubic", "quartic", "quintic"][k],
                  ")")
    ml == k && mr == k && return
    write(io, " (left, right multiplicities: $(ml), $(mr))")
end

function show(io::IO, t::AbstractKnotSet)
    write(io, "$(typeof(t).name)($(eltype(t))) of ")
    show_order_multiplicities(io, t)
    write(io, " on $(first(t)..last(t)) ($(numintervals(t)) intervals)")
end

# TODO: Implement via the `iterate` interface?
"""
    nonempty_intervals(t)

Return the indices of all intervals of the knot set `t` that are
non-empty.

# Examples

```jldoctest
julia> nonempty_intervals(ArbitraryKnotSet(3, [0.0, 1, 1, 3, 4, 6], 1, 3))
4-element Array{Int64,1}:
 1
 3
 4
 5
```

"""
nonempty_intervals(t::AbstractKnotSet) =
    [j for j in 1:length(t)-1
     if t[j] ≠ t[j+1]]

"""
    find_interval(t, x[, i=ml])

Find the interval in the knot set `t` that includes `x`, starting from
interval `i` (which by default is the first non-zero interval of the
knot set). The search complexity is linear, but by storing the result
and using it as starting point for the next call to `find_interval`,
the knot set need only be traversed once.
"""
function find_interval(t::AbstractKnotSet{k,ml,mr,T}, x, i=ml) where {T,k,ml,mr}
    (x < first(t) || x > last(t) || i > length(t) || x < t[i]) && return nothing
    x == last(t) && return length(t) - mr
    for r ∈ i:length(t)
        t[r] > x && return r-1
    end
    @assert false
end

const RightContinuous{T} = Interval{:closed,:open,T}

"""
    within_support(x, t, j)

Return the indices of the elements of `x` that lie withing the compact
support of the `j`th basis function (enumerated `1..n`), given the
knot set `t`. For each index of `x` that is covered, the index `k` of
the interval within which `x[i]` falls is also returned.
"""
function within_support(x::AbstractVector, t::AbstractKnotSet, j::Integer)
    isempty(x) && return 1:0
    k = order(t)
    # The last interval includes the right endpoint as well, whereas
    # the other intervals are only right-continuous.
    IntervalKind(i) = i == numintervals(t) ? ClosedInterval : RightContinuous
    ml = leftmultiplicity(t)
    supports = [(findall(in(IntervalKind(i-ml+1)(t[i], t[i+1])), x), i)
                for i = j:j+k-1]
    filter(s -> !isempty(s[1]), supports)
end

# * Specific knot sets
# ** Arbitrary knot set
"""
    ArbitraryKnotSet{k,ml,mr}(t)

An arbitrary knot set of order `k` and left and right multiplicities
of `ml` and `mr`, respectively. The knot set is specified by the
non-decreasing vector `t`; the same knot may appear multiple times,
which influences the continuity of the B-splines at that location.
"""
struct ArbitraryKnotSet{k,ml,mr,T,V<:AbstractVector{T}} <: AbstractKnotSet{k,ml,mr,T}
    t::V
    ArbitraryKnotSet{k,ml,mr}(t::V) where {k,ml,mr,T,V<:AbstractVector{T}} =
        new{k,ml,mr,T,V}(sort(t))
end

"""
    ArbitraryKnotset(k, t[, ml=k, mr=k])

Construct an order-`k` arbitrary knot set, with the locations of the
knots given by non-decreasing vector `t`.
"""
function ArbitraryKnotSet(k::Integer, t::AbstractVector,
                          ml::Integer=k, mr::Integer=k)
    assert_multiplicities(k, ml, mr, t)
    ArbitraryKnotSet{k,ml,mr}(t)
end

# ** Linear knot set
"""
    LinearKnotSet{k,ml,mr}(t)

A knot set of order `k` and left and right multiplicities of `ml` and
`mr`, respectively, whose knots are uniformly distributed according to
the range `t`. The interior basis functions are thus Cᵏ⁻²-continuous.
"""
struct LinearKnotSet{k,ml,mr,T,R<:AbstractRange{T}} <: AbstractKnotSet{k,ml,mr,T}
    t::R
    LinearKnotSet{k,ml,mr}(t::R) where {k,ml,mr,T,R<:AbstractRange{T}} =
        new{k,ml,mr,T,R}(t)
end

"""
    LinearKnotSet(k, a, b, N[, ml=k, mr=k])

Construct an order-`k` linear knot set spanning from `a` to `b`, with
`N` intervals.

# Examples

```jldoctest
julia> LinearKnotSet(2, 0, 1, 3)
6-element LinearKnotSet{2,2,2,Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}:
 0.0
 0.0
 0.3333333333333333
 0.6666666666666666
 1.0
 1.0
```
"""
function LinearKnotSet(k::Integer, a, b, N::Integer,
                       ml::Integer=k, mr::Integer=k)
    t = range(a, stop=b, length=N+1)
    assert_multiplicities(k, ml, mr, t)
    LinearKnotSet{k,ml,mr}(t)
end

distribution(::LinearKnotSet) = Uniform()

# ** Exponential knot set
"""
    ExpKnotSet{k,ml,mr}(exponents, base, t, include0)

A knot set of order `k` and left and right multiplicities of `ml` and
`mr`, respectively, whose knots are exponentially distributed
according to `t = base .^ exponents`, optionally including 0 as the
left endpoint."""
struct ExpKnotSet{k,ml,mr,T,R<:AbstractRange{T},TV<:AbstractVector{T}} <: AbstractKnotSet{k,ml,mr,T}
    exponents::R
    base::T
    t::TV
    include0::Bool
    ExpKnotSet{k,ml,mr}(exponents::R, base::T, t::TV, include0) where {k,ml,mr,T,R,TV} =
        new{k,ml,mr,T,R,TV}(exponents, base, t, include0)
end

"""
    ExpKnotSet(k, a, b, N[, ml=k, mr=k, base=10, include0=true])

Construct an order-`k` knot spanning from `base^a` to `base^b` in `N`
intervals, optionally including 0 as the left endpoint.

# Examples

```jldoctest
julia> ExpKnotSet(2, -4, 2, 7)
10-element ExpKnotSet{2,2,2,Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Array{Float64,1}}:
   0.0
   0.0
   0.0001
   0.001
   0.01
   0.1
   1.0
  10.0
 100.0
 100.0
```
"""
function ExpKnotSet(k::Integer, a::T, b::T, N::Integer,
                    ml::Integer=k, mr::Integer=k;
                    base::T=T(10), include0::Bool=true) where T
    exponents = range(a, stop=b, length=include0 ? N : N+1)
    t = base .^ exponents
    assert_multiplicities(k, ml, mr, t)
    ExpKnotSet{k,ml,mr}(exponents, eltype(t)(base), include0 ? vcat(0,t) : t, include0)
end

function show(io::IO, t::ExpKnotSet)
    write(io, "$(typeof(t).name)($(eltype(t))) of ")
    show_order_multiplicities(io, t)
    write(io, " on ")
    t.include0 && write(io, "0,")
    write(io, "$(t.base^first(t.exponents))..$(t.base^last(t.exponents)) ($(numintervals(t)) intervals)")
end

# function arcsin_knot_set(k::Integer, a::Integer, b::Integer, N::Integer)
#     N2 = N/2
#     arcsin.(range(-N2, stop=N2, lengthN+1)/N2)*(b-a)/π+(a+b)/2
# end

# arcsin_half_knot_set(a, b, N) = arcsin(linspace(0,N,N+1)/N)*(b-a)*2.0/π+a

# function exp_linear_knot_set(a, b, N)
#     ap = a != 0 ? a : 1e-1
#     @assert t[2][1]+t[2][2] == N
#     [logspace(log10(ap),log10(t[1]), t[2][1]); linspace(t[1],b,t[2][2]+1)[2:end]]
# end

@recipe function plot(t::AbstractKnotSet)
    markershape --> :circle
    y = similar(t)
    y′ = 1
    t′ = -Inf
    for (i,e) in enumerate(t)
        y′ = e == t′ ? (y′ + 1) : 1
        t′ = e
        y[i] = y′
    end
    collect(t),y
end

export ArbitraryKnotSet, LinearKnotSet, ExpKnotSet, order, numintervals, numfunctions, nonempty_intervals


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "libsixel"
version = v"1.8.6"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/saitoha/libsixel.git", "5db717dfef6fa327cd4025e7352550f63d20699c"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libsixel/

# This patch fixes https://github.com/johnnychen94/Sixel.jl/issues/5
if [[ "${target}" == *-mingw* ]]; then
    atomic_patch -p1 ../patches/01-libsixel-win32.patch
fi

./configure --prefix=${prefix} \
    --build=${MACHTYPE} \
    --host=${target} \
    --includedir=$WORKSPACE/destdir/include \
    --libdir=$WORKSPACE/destdir/lib \
    --enable-python=no
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libsixel", :libsixel)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
    GMC_NS_Ensemble

Abstract supertype for model ensembles intended to be subjected to Galilean Monte Carlo Nested Sampling ('GMC_NS'). Subtyped ensemble structs must implement the following fields as their interface with 'GMC_NS':

    mutable struct Minimal_Ensemble <: GMC_NS_Ensemble
        path::String #ensemble records serialized here

        models::Vector{<:GMC_NS_Model_Record} #records of model paths & lhs
        model_initλ::Function #GMC_NS_Model constructor

        contour<:AbstractFloat #current log_Li[end], init minimum([record.log_Li for record in models])
        log_Li::Vector{<:AbstractFloat} #likelihood of lowest-ranked model at iterate i, init [-Inf]
        log_Xi::Vector{<:AbstractFloat} #amt of prior mass included in ensemble contour at Li, init [0]
        log_wi::Vector{<:AbstractFloat} #width of prior mass covered in this iterate, init [-Inf]
        log_Liwi::Vector{<:AbstractFloat} #evidentiary weight of the iterate, init [-Inf]
        log_Zi::Vector{<:AbstractFloat} #ensemble evidence, init [typemin(FloatType)]
        Hi::Vector{<:AbstractFloat} #ensemble information, init [0]

        obs::Any #observations formatted for GMC_NS_Model likelihood Function
        priors::Vector{<:Distribution} #prior distributions on GMC_NS_Model parameters, formatted for model_initλ
        constants::Vector #fixed parameters to be passed to GMC_NS_Model constructor

        sample_posterior::Bool #switch for posterior_samples collection
        posterior_samples::Vector{<:GMC_NS_Model_Record}

        #GMC tuning settings
        GMC_tune_ι<:Integer #maximum number of initial τ tuning steps
        GMC_tune_ρ<:AbstractFloat #use ensemble size*ρ GMC samples per initial tuning step
        GMC_tune_μ<:Integer #tune τ on the basis of this number of the most recent steps
        GMC_tune_α<:AbstractFloat #tune τ such that the acceptance rate is between this value and 1.
        GMC_tune_β<:AbstractFloat #multiply τ by 1+β (sampling acceptance rate=1.) or 1-β (rate <α)

        #GMC Sampler settings
        GMC_timestep_η<:AbstractFloat #>0. to apply an perturbation term randomly drawn from Normal(0, η*τ) to τ at each step (can assist equilibration)
        GMC_reflect_η<:AbstractFloat #>0. to apply perturbation to reflection vectors, keep value small (can assist equilibration)
        GMC_exhaust_σ<:AbstractFloat #apply this scale parameter to the ensemble size to obtain the maximum number of GMC iterates to perform before terminating nested sampling
        GMC_chain_κ<:Integer #maximum chain length

        t_counter<:Integer #counter for next instantiated trajectory number
    end

"""
abstract type GMC_NS_Ensemble end

"""
    show(e<:GMC_NS_Ensemble)
    show(io,e<:GMC_NS_Ensemble)

Display the ensemble path, likelihood histogram, contour, maximum log likelihood, and log evidence (marginal likelihood).
"""
function Base.show(io::IO, e::GMC_NS_Ensemble; progress=false)
	livec=[model.log_Li for model in e.models]
	maxLH=maximum(livec)
	printstyled(io, "$(typeof(e)) @ $(e.path)\n", bold=true)
	msg = @sprintf "Contour: %3.6e MaxLH:%3.3e log Evidence:%3.6e" e.contour maxLH e.log_Zi[end]
    println(io, msg)
    if length(unique(livec)) > 1
        try
            hist=UnicodePlots.histogram(livec, title="Ensemble Likelihood Distribution")
            show(io, hist)
            println()
            progress && return(nrows(hist.graphics)+6)
        catch
            printstyled("HISTOGRAM UNAVAILABLE", color=:green)
            println()
            progress && return 3
        end
    else
        printstyled("HISTOGRAM UNAVAILABLE", color=:green)
        println()
        progress && return 3
    end
end


text:
using DelimitedFiles

puzzleinput = readlines("inputs/day16.txt")

Ranges = NamedTuple{(:fld, :rng),Tuple{String,Array{UnitRange{Int}}}}

function readinput(inpt::Vector{String})

    ranges = []

    function torange(s1, s2)::UnitRange{Int}
        atoi(s) = parse(Int, s)
        UnitRange(atoi(s1), atoi(s2))
    end

    for l in inpt[1:20]
        m = match(r"([^:]+): (\d+)-(\d+) or (\d+)-(\d+)", l)
        rng = [torange(m[2], m[3]), torange(m[4], m[5])]
        push!(ranges, (fld = m[1], rng = rng))
    end

    yourtkt = parse.(Int, split(inpt[23], ","))

    nearbytix = []
    for l in inpt[26:end]
        push!(nearbytix, parse.(Int, split(l, ",")))
    end

    return ranges, yourtkt, nearbytix
end

ranges, yourtkt, neartix = readinput(puzzleinput)

# First Half

function tkt_scan_err_rate(tix, rngs)
    # Returns the ticket scanning error for tickets in tix and
    # rules of ticket fields in rngs

    isvalid(fval) = foldl(|, [fval ∈ rs.rng[1] || fval ∈ rs.rng[2] for rs in rngs])
    # The field in a ticket is invalid nearby if contains values
    # that are not valid for any field

    sumofnotvalid(tkt) = sum(push!([fval for fval in tkt if !isvalid(fval)], 0))
    # Sum of tha values in invalid fields for a ticket

    sum(sumofnotvalid.(tix))
end

tkt_scan_err_rate(neartix, ranges)

# Second Half

function fldspos(tix, rngs, fldrngs)
    # Returns an array with the positions in the ticket
    # for the fields in field ranges

    isvalid(fval) = foldl(|, [fval ∈ rs.rng[1] || fval ∈ rs.rng[2] for rs in rngs])
    isok(ticket) = foldl(&, [isvalid(fval) for fval in ticket])
    # True if the ticket does not contain invalid values

    tixok = tix[isok.(tix)]
    # Valid tickets

    function check(c, f, tickets)
        # For every value in column c of ticket matrix checks
        # if satisfy the ranges for field f

        col = hcat(tickets...)[c, :]
        # Select column c of the matrix

        foldl(&, [v ∈ f.rng[1] || v ∈ f.rng[2] for v in col])
    end

    colsgoodforfld = Dict{String,Set{Int}}()
    numofcols = length((hcat(tix...)[:, 1]))
    for f in rngs
        # For every field in the rules collects, in a set, the columns
        # in the ticket matrix that have all values valid

        colsgoodforfld[f.fld] = Set()
        for c = 1:numofcols
            if check(c, f, tixok)
                push!(colsgoodforfld[f.fld], c)
            end
        end
    end

    function sievecols(colsgoodforfld)
        # Filter the cols good for a field untils only one remains for field

        function sieve_(todo, done)
            if length(todo) == 0
                # No more fiels to check, we have done

                return done
            end

            s = findfirst(x -> length(x) == 1, todo)
            # Search first field with only a column good for it

            if (isnothing(s))
                # No more fields that have only one column, we have done

                return done
            else
                singleton = todo[s]
                u = unique(singleton)[1]
                done[s] = u
                delete!(todo, s)
                # This field is good for shure, put it on done pool and
                # remove it from the ones to do

                for (f, cols) in todo
                    delete!(cols, u)
                    # The column of this field is taken, remove it from

                end
                sieve_(todo, done)
            end
        end
        ∅ = Dict()
        sieve_(colsgoodforfld, ∅)
    end

    positions = sievecols(colsgoodforfld)

    [unique(positions[f.fld])[1] for f in fldrngs if haskey(positions, f.fld)]
    # Extract only the positions of the requested fields in fldrngs

end

@time prod(yourtkt[fldspos(neartix, ranges, ranges[1:6])])


text:
using GaussianProcesses
using Base.Test

# Add tests here...
println("Running test_utils.jl...")
include("test_utils.jl")

println("Running test_GP.jl...")
include("test_GP.jl")

println("Running test_kernels.jl...")
include("test_kernels.jl")


text:
julia> [i + j for i in 1:5, j in 1:5]
5×5 Array{Int64,2}:
 2  3  4  5   6
 3  4  5  6   7
 4  5  6  7   8
 5  6  7  8   9
 6  7  8  9  10

julia> [i + j for i in 1:5, j in 1:5, k in 1:2]
5×5×2 Array{Int64,3}:
[:, :, 1] =
 2  3  4  5   6
 3  4  5  6   7
 4  5  6  7   8
 5  6  7  8   9
 6  7  8  9  10

[:, :, 2] =
 2  3  4  5   6
 3  4  5  6   7
 4  5  6  7   8
 5  6  7  8   9
 6  7  8  9  10


text:
using ITensors
using ITensorGLMakie
using ReferenceTests
using Test

@testset "Basic test for ITensorGLMakie" begin
  extension = "png"

  N = 10
  s(n) = Index([QN("Sz", 0) => 1, QN("Sz", 1) => 1]; tags="S=1/2,Site,n=$n")
  l(n) = Index([QN("Sz", 0) => 10, QN("Sz", 1) => 10]; tags="Link,l=$n")
  h(n) = Index([QN("Sz", 0) => 5, QN("Sz", 1) => 5]; tags="ham,Link,l=$n")
  s⃗ = [s(n) for n in 1:N]
  l⃗ = [l(n) for n in 1:(N - 1)]
  h⃗ = [h(n) for n in 1:(N - 1)]

  # Add some more indices between two of the tensors
  x = Index([QN("Sz", 0) => 2]; tags="X")
  y = Index([QN("Sz", 0) => 2]; tags="Y")

  n = 2
  ψn1n2 = randomITensor(l⃗[n - 1], s⃗[n], s⃗[n + 1], l⃗[n + 1], dag(x), dag(y))
  hn1 = randomITensor(dag(h⃗[n - 1]), s⃗[n]', dag(s⃗[n]), h⃗[n], x, y)
  hn2 = randomITensor(dag(h⃗[n]), s⃗[n + 1]', dag(s⃗[n + 1]), h⃗[n + 1])
  ELn0 = randomITensor(l⃗[n - 1]', h⃗[n - 1], dag(l⃗[n - 1]))
  ERn2 = randomITensor(l⃗[n + 1]', dag(h⃗[n + 1]), dag(l⃗[n + 1]))

  tn = [ELn0, ψn1n2, hn1, hn2, ERn2]

  R = @visualize figR ELn0 * ψn1n2 * hn1 * hn2 * ERn2
  R1 = @visualize figR1 ELn0 * ψn1n2 * hn1
  R2 = @visualize figR2 R1 * hn2 * ERn2 vertex_labels = ["T1", "T2", "T3"]
  T = @visualize figT ELn0

  @test R ≈ ELn0 * ψn1n2 * hn1 * hn2 * ERn2
  @test R1 ≈ ELn0 * ψn1n2 * hn1
  @test R2 ≈ ELn0 * ψn1n2 * hn1 * hn2 * ERn2
  @test T == ELn0

  fig_tn = @visualize_noeval tn

  by = extension == "png" ? psnr_equality(0.5) : isequal

  @test_reference "references/R.$extension" figR by = by
  @test_reference "references/R1.$extension" figR1 by = by
  @test_reference "references/R2.$extension" figR2 by = by
  @test_reference "references/tn.$extension" fig_tn by = by
  @test_reference "references/T.$extension" figT by = by

  R = @visualize fig_grid ELn0 * ψn1n2 * hn1 * hn2 * ERn2
  R1 = @visualize! fig_grid[1, 2] ELn0 * ψn1n2 * hn1
  R2 = @visualize! fig_grid[2, 1] R1 * hn2 * ERn2 vertex_labels = ["T1", "T2", "T3"]
  @visualize_noeval! fig_grid[2, 2] tn

  # XXX: Broken, passes locally but fails on CI with:
  # Warning: test fails because PSNR -0.6602330207824707 < 1
  #@test_reference "references/grid.$extension" fig_grid by=by

  @test_throws DimensionMismatch @visualize fig R1 * hn2 * ERn2 vertex_labels = ["T1", "T2"]
end


text:
module TensorOperators

using Base.LinAlg.BLAS
using ArrayViews

import Base: rand, randn!, rand!, zeros, push!, length, size

# --- Devices and types

typealias Float AbstractFloat

abstract Device

denseRealArray(::TypeVar) = Union{DenseArray{Float64}, DenseArray{Float32}}
denseRealArray(dim::Int) = Union{DenseArray{Float32,dim}, DenseArray{Float64, dim}}
const DenseRealTensor = denseRealArray(TypeVar(:N))
const DenseRealMatrix = denseRealArray(2)
const DenseRealVector = denseRealArray(1)

denseTensor(::TypeVar, F::TypeVar, N::Union{TypeVar,Int}) = DenseArray{F, N}
denseTensor{D<:Device}(::Type{D}, F, ::Int) =  error("Cannot compute type of dense matrix of ($D, $F)")

denseMatrix(D, F) = denseTensor(D, F, 2)
denseVector(D, F) = denseTensor(D, F, 1)

denseRealMatrix(::TypeVar, ::TypeVar) = DenseRealMatrix
denseRealMatrix{F<:AbstractFloat}(D, ::Type{F}) = denseMatrix(D, F)
denseRealVector(::TypeVar, ::TypeVar) = DenseRealVector
denseRealVector{F<:AbstractFloat}(D, ::Type{F}) = denseVector(D, F)

const IntVector = Array{Int, 1}

arrayOf{D<:Device}(::Type{D}, F, ::Integer) =  error("Cannot compute matrix of ($D, $F)")

vectorOf{D<:Device}(::Type{D}, F) = arrayOf(D, F, 1)
matrixOf{D<:Device}(::Type{D}, F) = arrayOf(D, F, 2)

# --- CPU

immutable CPUDevice <: Device end
export CPUDevice, cpu

cpu = CPUDevice()

arrayOf{F<:Float}(::Type{CPUDevice}, ::Type{F}, dims::Int64) = DenseArray{F, dims}

denseTensor{F}(::Type{CPUDevice}, ::Type{F}, dims::Int64) = DenseArray{F, dims}

array{F<:Float}(d::CPUDevice, ::Type{F}, dims::Int64...) = Array(F, dims...)
rand{F<:Float}(d::CPUDevice, ::Type{F}, dims::Int64...) = rand(F, dims...)
zeros{F<:Float}(d::CPUDevice, ::Type{F}, dims::Int64...) = zeros(F, dims...)

# --- Layers ----

abstract Layer

"""Computes the gradient wrt the parameters and the input. Returns the gradient wrt to the inputs"""
function backward!(m::Layer, input, gradOutput, scale::Float64=1.)
    update_gradient!(m, input, gradOutput, scale)
    return compute_inputgradient!(m, input, gradOutput)
end

""" Generic function for initializing the gradient """
function init_gradient!(m::Layer)
    for v in parameters(m)
        if !isnull(v.gradient)
            init_gradient!(v)
        end
    end
end

function init!(m::Layer)
    for v in parameters(m)
        init!(v)
    end
end

export init_gradient!, forward!, backward!, init!, update_gradient!, compute_inputgradient!

# --- Includes

# Misc
include("utils.jl")
include("parameters.jl")

#
include("cuda.jl")

# Optimization
include("optimization.jl")

# Container modules
include("containers.jl")

# Layers
include("layers/linear.jl")
include("layers/convolution.jl")
include("layers/bhsm.jl")
include("layers/lookup.jl")

# Transfer functions
include("transfer.jl")

# Costs
include("cost.jl")

end


text:
using Test
using FWFTables
using FixedSizeStrings

bla = readbla("../testdata/testbla.bla")

@test length(bla) == 3
@test bla[1] == Varspec("var1", Union{Missing, Int64}, 1, 9, 0)
@test bla[2] == Varspec("var2", FixedSizeStrings.FixedSizeString, 10, 8, 0)
@test bla[3] == Varspec("var3", Float64, 18, 8, 2)

max_stars_repo_name: HenricoWitvliet/FWFTables


text:
#####
##### Game trace
#####

"""
    Trace{Game, State}

An object that collects all states visited during a game, along with the
rewards obtained at each step and the successive player policies to be used
as targets.

# Constructor

    Trace{Game}(initial_state)

"""
mutable struct Trace{Game, State}
  states :: Vector{State}
  policies :: Vector{Vector{Float64}}
  rewards :: Vector{Float32}
  function Trace{G}(init_state) where G
    return new{G, GI.State(G)}([init_state], [], [])
  end
end

GameType(::Trace{Game}) where Game = Game

function trace_invariant(t::Trace)
  return length(t.policies) == length(t.rewards) == length(t.states) - 1
end

"""
    Base.push!(t::Trace, π, r, s)

Add a (target policy, reward, new state) triple to a trace.
"""
function Base.push!(t::Trace, π, r, s)
  push!(t.states, s)
  push!(t.policies, π)
  push!(t.rewards, r)
end

function Base.length(t::Trace)
  return length(t.rewards)
end

function total_reward(t::Trace, gamma=1.)
  return sum(gamma^(i-1) * r for (i, r) in enumerate(t.rewards))
end

function debug_trace(t::Trace{Game}) where Game
  n = length(t)
  for i in 1:n
    println("Transition $i:")
    game = Game(t.states[i])
    GI.render(game)
    for (a, p) in zip(GI.available_actions(game),  t.policies[i])
      print("$(GI.action_string(Game, a)): $(fmt(".3f", p))  ")
    end
    println("")
    println("Obtained reward of: $(t.rewards[i]).")
    println("")
  end
  println("Showing final state:")
  GI.render(Game(t.states[n + 1]))
end


text:
"""

Convert XXX.mmd + XXX.jl into a XXX.md file.
This allows values in XXX.jl to be inserted into the mustache template XXX.mmd to
dynamically render an md file. Examples of use include randomizing questions, including fancy graphics, ...

If converting into many formats this should only be run once, as otherwise the different formats will possible differences when randomization is involved.

"""
function mmd_to_md(fname::AbstractString; kwargs...)

    bname = basename(fname)
    occursin(r"\.mmd$", bname) || error("this is for mmd template files")
    bname = replace(bname, r"\.mmd$" => "")

    jl = replace(fname,".mmd" => ".jl")
    md = replace(fname,".mmd" => ".md")
    mmd = fname

    ## do this only if md file older than either .mmd or .jl
    if !isfile(md) || (mtime(mmd) > mtime(md)) || (mtime(jl) > mtime(md))
        include("$bname.jl")

        tpl = Mustache.template_from_file(fname)

        io = open(md, "w")
        write(io, Mustache.render(tpl, getfield(Main,Symbol(bname))))
        close(io)
    end

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct AvailableProvidersListCity <: SwaggerModel
    cityName::Any # spec type: Union{ Nothing, String } # spec name: cityName
    providers::Any # spec type: Union{ Nothing, Vector{String} } # spec name: providers

    function AvailableProvidersListCity(;cityName=nothing, providers=nothing)
        o = new()
        validate_property(AvailableProvidersListCity, Symbol("cityName"), cityName)
        setfield!(o, Symbol("cityName"), cityName)
        validate_property(AvailableProvidersListCity, Symbol("providers"), providers)
        setfield!(o, Symbol("providers"), providers)
        o
    end
end # type AvailableProvidersListCity

const _property_map_AvailableProvidersListCity = Dict{Symbol,Symbol}(Symbol("cityName")=>Symbol("cityName"), Symbol("providers")=>Symbol("providers"))
const _property_types_AvailableProvidersListCity = Dict{Symbol,String}(Symbol("cityName")=>"String", Symbol("providers")=>"Vector{String}")
Base.propertynames(::Type{ AvailableProvidersListCity }) = collect(keys(_property_map_AvailableProvidersListCity))
Swagger.property_type(::Type{ AvailableProvidersListCity }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_AvailableProvidersListCity[name]))}
Swagger.field_name(::Type{ AvailableProvidersListCity }, property_name::Symbol) =  _property_map_AvailableProvidersListCity[property_name]

function check_required(o::AvailableProvidersListCity)
    true
end

function validate_property(::Type{ AvailableProvidersListCity }, name::Symbol, val)
end


text:
export heaviside, nf, nb

function heaviside(x)
   0.5 * (sign(x) + 1)
end

function nf(x,T=0,d=0)
    if d == 0
        return nfBrute(x,T)
    end
    if d == 1
        return nf1Brute(x,T)
    end
    if d == -1
        return nfiBrute(x,T)
    end
    throw(ArgumentError)
end

function nfBrute(x,T)
    if T == 0
        if x>0 return 0
        elseif x<0 return 1
        else return 0.5
        end
    end
    return 1/(1+exp(x/T))
end

function nf1Brute(x,T)
    if T == 0
        if x == 0
            return NaN
        else
            return 0
        end
    end
    return -1/( (exp(x/T)+2+exp(-x/T)) * T)
end

function nfiBrute(x,T)
    if T == 0
        if x == 0
            return NaN
        else
            return  1. /(x*pi)
        end
    end
    return imag(trigamma(0.5-im*x/(2*pi*T))/(2*pi^2*T))
end

function nb(x,T=0)
    if T == 0
        if x>0 return 0
        elseif x<0 return -1
        else return NaN
        end
    end
    return 1/(exp(x/T)-1)
end


text:
"""
```
compute_coeffs(coeff_func, current_coeffs)
```
calculates the coefficients for the stencil of UpwindDifference operators.
"""
function compute_coeffs!(coeff_func::Number, current_coeffs::AbstractVector{T}) where {T<:Number}
    return current_coeffs .+= coeff_func
end

function compute_coeffs!(coeff_func::AbstractVector{T}, current_coeffs::AbstractVector{T}) where {T<:Number}
     return current_coeffs[:] += coeff_func
end

# Coefficient functions when coeff_func is a Function and current_coeffs exists
function compute_coeffs!(coeff_func::Function, current_coeffs::AbstractVector{T}) where {T<:Number}
    if hasmethod(coeff_func, (Vector{T},))
        current_coeffs[:] = coeff_func(current_coeffs)
    else
        map!(coeff_func, current_coeffs, current_coeffs)
    end
    return current_coeffs
    # if hasmethod(coeff_func, (Int,)) # assume we want to provide the index of coefficients to coeff_func
    #     for i = 1:length(current_coeffs)
    #         current_coeffs[i] += coeff_func(i)
    #     end
    #     return current_coeffs
    # elseif hasmethod(coeff_func, (Vector{Int},))  # assume we want to provide the index of coefficients to coeff_func
    #     current_coeffs[:] += coeff_func(collect(1:length(current_coeffs)))
    #     return current_coeffs
    # elseif hasmethod(coeff_func, (UnitRange{Int},)) # assume we want to provide the index of coefficients to coeff_func
    #     current_coeffs[:] += coeff_func(1:length(current_coeffs))
    #     return current_coeffs
    # elseif hasmethod(coeff_func, (Float64,)) # assume we want coeff_func to operate on existing coefficients
    #     map!(coeff_func, current_coeffs, current_coeffs)
    #     return current_coeffs
    # elseif hasmethod(coeff_func, (Vector{Float64},)) # assume we want to coeff_func to operate on existing coefficients
    #     current_coeffs[:] = coeff_func(current_coeffs)
    #     return current_coeffs
    # else
    #     error("Coefficient functions with the arguments of $coeff_func have not been implemented.")
    # end
end

compute_coeffs(coeff_func::Number, current_coeffs::AbstractVector{T}) where {T<:Number} = current_coeffs .+ coeff_func
compute_coeffs(coeff_func::AbstractVector{T}, current_coeffs::AbstractVector{T}) where {T<:Number} = coeff_func + current_coeffs

function compute_coeffs(coeff_func::Function, current_coeffs::AbstractVector{T}) where {T<:Number}
    if hasmethod(coeff_func, (Vector{T},))
        return coeff_func(current_coeffs)
    else
        return map(coeff_func, current_coeffs)
    end
    # if hasmethod(coeff_func, (Int,))
    #     return current_coeffs + map(coeff_func, collect(1:length(current_coeffs)))
    # elseif hasmethod(coeff_func, (Vector{Int},))
    #     return current_coeffs + coeff_func(collect(1:length(current_coeffs)))
    # elseif hasmethod(coeff_func, (UnitRange{Int},))
    #     return current_coeffs + coeff_func(1:length(current_coeffs))
    # elseif hasmethod(coeff_func, (Float64,)) # assume we want coeff_func to operate on existing coefficients
    #     return map(coeff_func, current_coeffs)
    # elseif hasmethod(coeff_func, (Vector{Float64},)) # assume we want to coeff_func to operate on existing coefficients
    #     return coeff_func(current_coeffs)
    # else
    #     error("Coefficient functions with the arguments of coeff_func $coeff_func have not been implemented.")
    # end
end


text:
module Forecast

using CSV, Distributions, ColorSchemes, DataFrames, DataFramesMeta, Dates, CodecZlib,
    HypothesisTests, LinearAlgebra, Optim, Plots, PrettyTables, RecipesBase, Statistics, StatsBase

# types
export AR, CCF, STL

# methods
export acf, ar, arsim, boxcox, ccf, d, hma, iboxcox, loess, p, pacf, forecast,
       setnames!, sma, stl, transform, summarize

# datasets
export air, co2, london, quakes, seaborne

# types
include("AR.jl")
include("CCF.jl")
include("FORECAST.jl")
include("STL.jl")
include("SUMMARIZE.jl")

# source files
include("acf.jl")
include("ar.jl")
include("arsim.jl")
include("boxcox.jl")
include("ccf.jl")
include("d.jl")
include("datasets.jl")
include("hma.jl")
include("loess.jl")
include("p.jl")
include("pacf.jl")
include("forecast_ar.jl")
include("sma.jl")
include("stl.jl")
include("summarize.jl")
include("utils.jl")
include("utils_datetime.jl")

# plot recipes
##  type
include("plot_CCF.jl")
include("plot_DataFrame.jl")
include("plot_FORECAST.jl")
include("plot_STL.jl")

## user
include("splot.jl")

"""
Collection of methods for Time Series analysis

Featured Methods:

    acf:        Auto-correlation or auto-covariance of univariate series.
    ar:         Multivariate Autoregressive Model.
    arsim:      Simulated Multivariate Autoregressive Model.
    boxcox:     Box-Cox Transformations.
    ccf:        Cross-correlation or cros-covariance of two univariate series.
    d:          Lagged differences of a given order for Vector and Array.
    forecast:   Forecast values of fitted time series models.
    hma:        Henderson moving average filter.
    iboxcox:    Inverse Box-Cox Transformations.
    loess:      Locally estimated scatterplot smoothing.
    p:          Reverse lagged differences of a given order for types Vector and Array.
    pacf:       Partial Auto-correlation function.
    sma:        Simple moving average.
    splot:      Plot a seasonal plot for types Vector and TimeArray.
    stl:        Seasonal and Trend decomposition using loess.
    summarize:  Statistical summary.
"""
Forecast

end


text:
"""
    set_dual_model_sense!(dual_model::MOI.ModelLike, model::MOI.ModelLike)

Set the dual model objective sense.
"""
function set_dual_model_sense(dual_model::MOI.ModelLike, primal_model::MOI.ModelLike)::Nothing where T
    # Get model sense
    primal_sense = MOI.get(primal_model, MOI.ObjectiveSense())
    if primal_sense == MOI.FEASIBILITY_SENSE
        error(primal_sense, " is not supported") # Feasibility should be supported?
    end
    # Set dual model sense
    dual_sense = (primal_sense == MOI.MIN_SENSE) ? MOI.MAX_SENSE : MOI.MIN_SENSE
    MOI.set(dual_model, MOI.ObjectiveSense(), dual_sense)
    return
end

# Primals
"""
    PrimalObjective{T}

Primal objective is defined as a `MOI.ScalarAffineFunction`
"""
struct PrimalObjective{T}
    saf::SAF{T}

    function PrimalObjective{T}(obj::SAF{T}) where T
        canonical_obj = MOIU.canonical(obj)
        if isempty(canonical_obj.terms)
            error("Dualization does not support models with no variables in the objective function.")
        end
        return new(canonical_obj)
    end
end

# Duals
"""
    DualObjective{T}

Dual objective is defined as a `MOI.ScalarAffineFunction`.
"""
struct DualObjective{T}
    saf::SAF{T}
end

function get_saf(objective::Union{PrimalObjective{T}, DualObjective{T}})::SAF{T} where T
    return objective.saf
end

"""
    get_primal_obj_coeffs(model::MOI.ModelLike)

Get the coefficients from the primal objective function and return a `PrimalObjective{T}`.
"""
function get_primal_objective(primal_model::MOI.ModelLike)
    T = MOI.get(primal_model, MOI.ObjectiveFunctionType())
    return _get_primal_objective(MOI.get(primal_model, MOI.ObjectiveFunction{T}()))
end

function _get_primal_objective(obj_fun::SAF{T}) where T
    return PrimalObjective{T}(obj_fun)
end

# Float64 is default while I don't know how to take other types
_get_primal_objective(obj_fun::SVF) = _get_primal_objective(obj_fun, Float64)
function _get_primal_objective(obj_fun::SVF, T::Type)
    return PrimalObjective{T}(SAF{T}(obj_fun))
end

# You can add other generic _get_primal_obj_coeffs functions here

"""
    set_dual_objective(dual_model::MOI.ModelLike, dual_objective::DualObjective{T})::Nothing where T

Add the objective function to the dual model.
"""
function set_dual_objective(dual_model::MOI.ModelLike, dual_objective::DualObjective{T})::Nothing where T
    # Set dual model objective function
    MOI.set(dual_model, MOI.ObjectiveFunction{SAF{T}}(),
            get_saf(dual_objective))
    return
end

"""
    get_dual_objective(dual_model::MOI.ModelLike, dual_obj_affine_terms::Dict,
                       primal_objective::PrimalObjective{T})::DualObjective{T} where T

build the dual model objective function from the primal model.
"""
function get_dual_objective(dual_model::MOI.ModelLike, dual_obj_affine_terms::Dict,
                            primal_objective::PrimalObjective{T})::DualObjective{T} where T

    sense = MOI.get(dual_model, MOI.ObjectiveSense()) # Get dual model sense

    num_objective_terms = length(dual_obj_affine_terms)
    term_vec = Vector{T}(undef, num_objective_terms)
    vi_vec   = Vector{VI}(undef, num_objective_terms)
    for (i, var) in enumerate(keys(dual_obj_affine_terms)) # Number of constraints of the primal model
        term = dual_obj_affine_terms[var]
        # Add positive terms bi if dual model sense is max
        term_vec[i] = (sense == MOI.MAX_SENSE ? -1 : 1) * term
        # Variable index associated with term bi
        vi_vec[i] = var
    end
    saf_dual_objective = MOI.ScalarAffineFunction(
                         MOI.ScalarAffineTerm.(term_vec,
                                               vi_vec),
                                               MOI.constant(get_saf(primal_objective)))
    return DualObjective{T}(saf_dual_objective)
end