text:
module Mongo

include("../deps/deps.jl")

using LibBSON

import Base.show,
    Base.convert,
    Base.start,
    Base.next,
    Base.count,
    Base.find,
    Base.done

ccall(
    (:mongoc_init, libmongoc),
    Void, ()
    )
atexit() do
    ccall((:mongoc_cleanup, libmongoc), Void, ())
end

const NakedDict = Union{Pair,Tuple}

include("MongoClient.jl")
include("MongoCollection.jl")
include("MongoCursor.jl")
include("query.jl")
include("bulk_insert.jl")
include("aggregate.jl")
include("conversions.jl")

end


text:
using CohortIdentification
using Test

@testset "CohortIdentification.jl" begin
    # Write your tests here.
end

ext: jl

max_stars_repo_name: dorbarker/cohort-identification

max_issues_repo_name: dorbarker/cohort-identification

max_forks_repo_name: dorbarker/cohortid


text:
n = 10

#for elty in (Float32, Float64, ComplexF32, ComplexF64)

#    A0 = randn(elty, n, n)
A0 = randn(n, n)

A = copy(A0)
ws = QrpWs(A)

geqp3!(A, ws)

@static if VERSION < v"1.7"
    target = qr(A0, Val(true))
else
    target = qr(A0, ColumnNorm())
end

A = randn(n, 15)
vA = zeros(n, 2)
vA .= view(A, :, [7, 9])
ws = QrWs(vA)
geqrf_core!(vA, ws)
ormqr_core!('L', vA', A, ws)


text:
# DEPRECATED, no longer needed
module Iterables
export IterateEmpty, IterateSingleton, Iterable
import ProxyInterfaces
using DataTypesBasic
using TypeClasses.Utils

# Iterable Wrapper
# ================

# TODO get rid of TypeTag ElemType
"""
wrapper to clearly indicate that something should be treated as an Iterable
"""
struct Iterable{IterType}
  iter::IterType  # we follow the naming of Base.Generator
end
Iterable() = Iterable(IterateEmpty())  # empty iter
Iterable(it::Iterable) = it  # don't nest Iterable wrappers

ProxyInterfaces.iterator(::Type{Iterable{IterT}}) where {IterT} = IterT
ProxyInterfaces.iterator(it::Iterable) = it.iter
ProxyInterfaces.@iterator Iterable  # includes map

# generic convert method
Base.convert(::Type{Iterable{T}}, x::Iterable{T}) where T = x
function Base.convert(::Type{<:Iterable}, x)
  @assert(isiterable(x), "Only iterables can be converted to Iterable, please overload `Base.isiterable` respectively")
  Iterable(x)
end

# Iterable Helpers
# ================

struct IterateEmpty{ElType} end
# Union{} makes most sense as default type, as it has no element at all
# also typeinference should work correctly as `promote_type(Union{}, Int) == Int`
IterateEmpty() = IterateEmpty{Union{}}()

Base.iterate(::IterateEmpty) = nothing

Base.IteratorSize(::Type{<:IterateEmpty}) = Base.HasLength()
Base.length(::IterateEmpty) = 0
Base.IteratorEltype(::Type{<:IterateEmpty}) = Base.HasEltype()
Base.eltype(::Type{IterateEmpty{T}}) where T = T

struct IterateSingleton{T}
  value::T
end
Base.iterate(iter::IterateSingleton) = iter.value, nothing
Base.iterate(iter::IterateSingleton, state) = nothing

Base.IteratorSize(::Type{<:IterateSingleton}) = Base.HasLength()
Base.length(::IterateSingleton) = 1
Base.IteratorEltype(::Type{<:IterateSingleton}) = Base.HasEltype()
Base.eltype(::Type{IterateSingleton{T}}) where T = T

end # module


text:
#=

  Simple PERT model in ConstraintSolver.jl

  From Pascal van Hentenryck
  "Scheduling and Packing In the Constraint Language cc(FD)", page 7f
  http://citeseer.ist.psu.edu/300151.html

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function pert(problem, print_solutions=true,all_solutions=true,timeout=6)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   # "all_solutions"=> all_solutions,
                                                            "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>timeout,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))
    times        = problem[:times]
    dependencies = problem[:dependencies]
    max_time     = problem[:max_time]
    n = length(times)

    @variable(model, 0 <= start[1:n] <= max_time, Int)
    @variable(model, 0 <= s_end <= max_time, Int)

    for (d1,d2) in eachrow(dependencies)
        @constraint(model, start[d1] >= start[d2]+times[d2])
    end

    @constraint(model, s_end == start[n]) # to minimize

    @objective(model,Min,s_end)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    num_sols = 0
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                start_val = convert.(Integer,JuMP.value.(start; result=sol))
                s_end_val = convert.(Integer,JuMP.value.(s_end; result=sol))
                println("start:$start_val")
                println("s_end:$s_end_val")
                println()

            end
        end
    else
        println("status:$status")
    end

    return status, num_sols
end

problem = Dict(
             # Note: There is no i
             # a  b  c  d  e  f  g  h  j  k  Send
    :times => [7, 3, 1, 8, 1, 1, 1, 3, 2, 1, 1],
    # Dependencies
    :dependencies => resize_matrix(
            [[2,1],  # Sb >= Sa + 7
             [4,1],  # Sd >= Sa + 7
             [3,2],  # Sc >= Sb + 3
             [5,3],  # Se >= Sc + 1
             [5,4],  # Se >= Sd + 8
             [7,3],  # Sg >= Sc + 1
             [7,4],  # Sg >= Sd + 8
             [6,4],  # Sf >= Sd + 8
             [6,3],  # Sf >= Sc + 1
             [8,6],  # Sh >= Sf + 1
             [9,8],  # Sj >= Sh + 3
             [10,7], # Sk >= Sg + 1
             [10,5], # Sk >= Se + 1
             [10,9], # Sk >= Sj + 2
             [11,10] # Send >= Sk + 1
            ]),
    :max_time => 30,
)

@time pert(problem)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
mutable struct HopfieldQAP <: HopfieldNet
    # s::Array{Float64, N} where {N}
    s::BitArray{N} where {N}
    W::Array{Float64, N} where {N}
    λ::Float64
    # D::Array{Float64, N} where {N}
    # F::Array{Float64, N} where {N}
    # A::Float64
    # B::Float64
    # q::Float64
    ϵ::Float64
end

function δ(ind...)
    all(y -> y==first(ind), ind)
end

function HopfieldQAP(D::Array{Float64, Dim}, F::Array{Float64, Dim}
                    # , A=0.9, B=0.9, q=140
                    , ϵ=0.02
                    , λ = 1.1
    ) where {Dim}
    # A = B = 2

    # D & F are square matrices always
    N = (size(D)[1],size(F)[1]) #dimensions of the Hopfield Network

    tnet = DiscreteHopfieldNet(N)

    iters = Base.Iterators.product([collect(1:i) for i in N]...)
    for i in iters
        for j in iters
            # tnet.W[i..., j...] = -2(A * (1-δ(i[1],j[1]) * δ(i[2],j[2])) + B * δ(i[1], j[1]) * (1 - δ(i[2], j[2])) + (D[i[1], j[1]] * F[i[2], j[2]])/q)
            tnet.W[i..., j...] = D[i[1], j[1]] * F[i[2], j[2]]
        end
    end

    return HopfieldQAP(tnet.s , tnet.W
                        # , D, F
                        # , A, B, q
                        , ϵ
                        , λ
                        )
end

# This is a test function to return variable wieghts at each call
function weights(net::HopfieldQAP, i::Tuple)
    N = size(net.s)
    w = Array{Float64}(undef, N...)
    A = net.A
    B = net.B
    q = net.q
    D = net.D
    F = net.F
    iters = Base.Iterators.product([collect(1:i) for i in N]...)
    for j in iters
        w[j...] = -2(A * (1-δ(i[1],j[1]) * δ(i[2],j[2])) + B * δ(i[1], j[1]) * (1 - δ(i[2], j[2])) + (D[i[1], j[1]] * F[i[2], j[2]])/q)
    end

    return w

end

function energy(net::HopfieldQAP)
    N = size(net.s)
    N_2 = size(net.W)
    iters = Base.Iterators.product([collect(1:i) for i in N]...)
    e = 0.0
    for ind in iters
        a = []
        for i in 1:length(N_2)
            if i <= length(N)
                push!(a, ind[i])
            else
                push!(a, :)
            end
        endMMo
        # e += (dot(net.W[a...], net.s) + net.A + net.B) *net.s[ind...]
        e += dot(net.W[a...], net.s) * net.s[ind...] + net.λ * ((sum(net.s[:,ind[2]]) - 1)^2 + (sum(net.s[ind[1],:]) - 1)^2)
        # e += (dot(weights(net, ind), net.s) + net.A + net.B) * net.s[ind...]
    end
    return e

end

function update!(net::HopfieldQAP, ind::Tuple
                 , θ = 50
                 )
    N = size(net.s)
    N_2 = size(net.W)
    a = []
    for i in 1:length(N_2)
        if i <= length(N)
            push!(a, ind[i])
        else
            push!(a, :)
        end
    end
    # dotp = dot(net.W[a...], net.s) * net.s[ind...] - net.A - net.B
    # dotp = dot(weights(net, ind), net.s) * net.s[ind...] - net.A - net.B
    # net.s[ind...] = (1/ (1 + exp(dotp/net.ϵ)) )#> 0.5
    # net.s[ind...] = dot(net.W[a...], net.s) > θ #? +1 : -1
    fin = dot(net.W[a...], net.s) + net.λ * 2 * (sum(net.s[ind[1], : ]) + sum(net.s[:, ind[2]])) - 4 * net.λ
    net.s[ind...] = round(1 / (1 + exp(fin/net.ϵ)))
    return nothing
end

# Perform an update cycle on the whole network
function cycle!(net::HopfieldQAP
                , θ::Number = 50
                )
    N = size(net.s)
    # N_2 = size(net.W)
    randInd = [shuffle(1:i) for i in N]
    for i in Base.Iterators.product(randInd...)
        update!(net, i
                , θ
                )
    end
    return nothing
end

function isFeasible(net::HopfieldQAP)
    N = size(net.s)
    for i in 1:N[1]
        if sum(net.s[i,:]) != 1
            return false
        end
    end
    for j in i:N[2]
        if sum(net.s[:,j]) != 1
            return false
        end
    end
    return true
end

function lowEFeasible(net::HopfieldQAP, cyc::Int)
    e = 0.0
    f = isFeasible(net)

    for i in 1:cyc
        e = energy(net)
        cycle!(hnqap1)
        f = isFeasible(hnqap1)
end


text:
# Copyright (c) 2021, Joseph M. Rutherford
# All rights reserved.
#
# Code provided under the license contained in the LICENSE file.
#
# This file depends upon the Sudoku module being defined first.
using Test
using Random

function multiple_random_permutations!(seed,results)
    rng = MersenneTwister(seed)
    num_tests = size(results)[1]
    puzzle = Sudoku.SolvablePuzzle(3)
    for i = 1:num_tests
        permutation = Sudoku.random_permutation!(puzzle,rng)
        results[i,1:length(permutation)] = permutation[:]
    end
end

function random_permutations(num_tests)
    all_results = zeros(Int,(num_tests,4))

    # Plan to split the work across all worker
    num_tasks = Threads.nthreads() # number of workers
    # Make a separate RNG seed for each worker
    seed_rng = MersenneTwister(12345)
    rng_seeds = zeros(Int,(num_tasks,1))
    for i = 1:num_tasks
        rng_seeds[i] = rand(seed_rng,1000:100000)
    end
    workload = Int(num_tests/num_tasks)
    Threads.@threads for i = 1:num_tasks
        if i == num_tasks
            multiple_random_permutations!(
                rng_seeds[i],
                @view all_results[((i-1)*workload+1):num_tests,:])
        else
            multiple_random_permutations!(
                rng_seeds[i],
                @view all_results[((i-1)*workload+1):(i*workload),:])
        end
    end
    MAX_RESIDUAL_ERROR = 2.0/sqrt(num_tests)
    # Verify that the number of each permutation type is consistent
    column_1_counts = [count(i->i==j,@view all_results[:,1]) for j in 1:7]
    @test abs(column_1_counts[1] - (num_tests/7))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_1_counts[2] - (num_tests/7))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_1_counts[3] - (num_tests/7))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_1_counts[4] - (num_tests/7))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_1_counts[5] - (num_tests/7))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_1_counts[6] - (num_tests/7))/num_tests < MAX_RESIDUAL_ERROR

    # Verify that the number of each block selected for permutation 1:4 is consistent
    column_2_counts = zeros(Int,(3,))
    for i in 1:size(all_results)[1]
        if all_results[i,1] in 1:4
            column_2_counts[all_results[i,2]] += 1
        end
    end
    # 4/7 of the time we selected a block, and 1/3 of the those were 1, 2, and 3
    @test abs(column_2_counts[1] - (4*num_tests/21))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_2_counts[2] - (4*num_tests/21))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_2_counts[3] - (4*num_tests/21))/num_tests < MAX_RESIDUAL_ERROR

    # Verify that the subrow selected for permutations 1, 3 is consistent
    column_3_subrow_counts = zeros(Int,(3,))
    for i in 1:size(all_results)[1]
        if all_results[i,1] == 1 || all_results[i,1] == 3
            column_3_subrow_counts[all_results[i,3]] += 1
        end
    end
    # 2/7 of the time we selected a subrow, and 1/3 of the those were 1, 2, and 3
    @test abs(column_3_subrow_counts[1] - (2*num_tests/21))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_3_subrow_counts[2] - (2*num_tests/21))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_3_subrow_counts[3] - (2*num_tests/21))/num_tests < MAX_RESIDUAL_ERROR

    # Verify that the second subrow selected for permutations 1, 3 is consistent
    column_4_subrow_counts = zeros(Int,(3,))
    for i in 1:size(all_results)[1]
        if all_results[i,1] == 1 || all_results[i,1] == 3
            @test all_results[i,4] != all_results[i,3]
            column_4_subrow_counts[all_results[i,4]] += 1
        end
    end
    # 2/7 of the time we selected a subrow, and 1/3 of the those were 1, 2, and 3
    @test abs(column_4_subrow_counts[1] - (2*num_tests/21))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_4_subrow_counts[2] - (2*num_tests/21))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_4_subrow_counts[3] - (2*num_tests/21))/num_tests < MAX_RESIDUAL_ERROR

    # Verify that the second block selected for permutations 2, 4 is consistent
    column_3_block_counts = zeros(Int,(3,))
    for i in 1:size(all_results)[1]
        if all_results[i,1] == 2 || all_results[i,1] == 4
            @test all_results[i,3] != all_results[i,2]
            column_3_block_counts[all_results[i,3]] += 1
        end
    end
    # 2/7 of the time we selected a second block, and 1/3 of the those were 1, 2, and 3
    @test abs(column_3_block_counts[1] - (2*num_tests/21))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_3_block_counts[2] - (2*num_tests/21))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_3_block_counts[3] - (2*num_tests/21))/num_tests < MAX_RESIDUAL_ERROR

    # Verify that the first, second symbol selected for permutation 7 is consistent
    column_2_symbol_counts = zeros(Int,(9,))
    column_3_symbol_counts = zeros(Int,(9,))
    for i in 1:size(all_results)[1]
        if all_results[i,1] == 7
            @test all_results[i,3] != all_results[i,2]
            column_2_symbol_counts[all_results[i,2]] += 1
            column_3_symbol_counts[all_results[i,3]] += 1
        end
    end
    # 1/7 of the time we selected a symbol, and 1/9 of the those were 1:9
    @test abs(column_2_symbol_counts[1] - (1*num_tests/63))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_2_symbol_counts[2] - (1*num_tests/63))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_2_symbol_counts[3] - (1*num_tests/63))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_2_symbol_counts[4] - (1*num_tests/63))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_2_symbol_counts[5] - (1*num_tests/63))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_2_symbol_counts[6] - (1*num_tests/63))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_2_symbol_counts[7] - (1*num_tests/63))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_2_symbol_counts[8] - (1*num_tests/63))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_2_symbol_counts[9] - (1*num_tests/63))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_3_symbol_counts[1] - (1*num_tests/63))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_3_symbol_counts[2] - (1*num_tests/63))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_3_symbol_counts[3] - (1*num_tests/63))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_3_symbol_counts[4] - (1*num_tests/63))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_3_symbol_counts[5] - (1*num_tests/63))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_3_symbol_counts[6] - (1*num_tests/63))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_3_symbol_counts[7] - (1*num_tests/63))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_3_symbol_counts[8] - (1*num_tests/63))/num_tests < MAX_RESIDUAL_ERROR
    @test abs(column_3_symbol_counts[9] - (1*num_tests/63))/num_tests < MAX_RESIDUAL_ERROR
end

function symbol_swap()
    puzzle_2_input = Array{Int16}(undef,(4,4))
    puzzle_2_input[:] = [1,2,0,3, 0,3,1,2, 2,0,3,1, 3,1,2,0][:]
    p2 = Sudoku.SolvablePuzzle(2)
    Sudoku.assign_values!(p2,puzzle_2_input)
    for i = 1:length(p2.grid)
        if Sudoku.get_value(p2.grid[i]) == 0
            # force possibilities to be either index 4 or 1
            p = BitVector(undef,4)
            p.chunks[1] = 0x9
            p2.grid[i] = Sudoku.PuzzleEntry(p)
        end
    end
    # Interchange symbols 4 and 2
    puzzle_2_reference = Array{Int16}(undef,(4,4))
    puzzle_2_reference[:] = [1,4,0,3, 0,3,1,4, 4,0,3,1, 3,1,4,0][:]
    Sudoku.symbol_swap!(p2,4,2)
    for i = length(p2.grid)
        @test Sudoku.get_value(p2.grid[i]) == puzzle_2_reference[i]
        if puzzle_2_reference[i] == 0
            @test p2.grid[i].possibilities == 0x3
        else
            @test p2.grid[i].possibilities == (0x1 << (puzzle_2_reference[i]-1))
        end
    end
end

function random_puzzle()
    rng = MersenneTwister(123456)
    for i = 1:1000
        solution, puzzle = Sudoku.random_puzzle(3,rng,1000,50)
        @test Sudoku.valid_puzzle(Sudoku.as_values(solution))
        @test Sudoku.satisfies(Sudoku.as_values(puzzle),Sudoku.as_values(solution))
    end
end

max_stars_repo_name: joseph-m-rutherford/sudoku-julia

max_issues_repo_name: joseph-m-rutherford/sudoku-julia

max_forks_repo_name: joseph-m-rutherford/sudoku-julia


text:
function controlflow_1(a, b)
    if a == 10
        x = if b == 22
            7
        else
            8
        end
        for i=1:100
            x += i
            x -= 77
            if i == 77
                continue
            elseif i == 99
                break
            end
        end
        return x
    else
        return 77
    end
end

function staticparams_1{T1 <: Integer, T2}(a::T1, b::T2)
    x = if T2 <: Number
        T1(b)
    else
        one(T1)
    end
    f = rand(T2)
    f * x
end

"""
Will docs disturb macro_form's file/string based method to get source
"""
staticparams_short{T, T2}(a::T, b::T2) = a + b

#comment to test if this annoys macro_form file/string based method to get source
function static_annon{T}(a::T, b)
    c = a+b::T
    if c == 10
        for i=1:7
            c += rand(Bool) ? 1 : 2
        end
    end
    f = (x,y) -> (x::T)^y
    c, f
end


text:
using LinearAlgebra: norm, inv, eigvals
import ForwardDiff
using StaticArrays: setindex
using DynamicalSystems

# This code is written for `SVector`s, or out-of-place dynamical systems.
# Notation z is used for the vector (x..., p), i.e. mixed space
"""
    bifurcation_rule_form(ds) → f, J
Produce dynamic rule `f` and Jacobian function `J` in a two-argument form.
"""
function bifurcation_rule_form(ds::ContinuousDynamicalSystem)
    f = (x, p) -> ds.f(x, p, 0.0)
    J = (x, p) -> ds.jacobian(x, p, 0.0)
    return f, J
end

function bifurcation_rule_form(ds::DiscreteDynamicalSystem)
    error("TODO: Convert `f` to `f - x`.")
    f = (x, p) -> ds.f(x, p, 0.0)
    J = (x, p) -> ds.jacobian(x, p, 0.0)
    return f, J
end

# Corrector: Newton's method in mixed-space with extra fixed condition
"""
    corrector(zpred, zprev, f, J; max_steps = 200, δ = 0.9, ε = 1e-3) → (z*, success)
Implement the Newton algorithm to converge to the fixed point for given predicted
mixed state `zpred`, previously found state `zprev`, function `f`, Jacobian `J`
and stepping factor `δ`.

If iteration occurs for more than `max_steps` without convergence better than `ε`,
then iteration stops and `success` will be `false`.
"""
function corrector(zpred, zprev, f, J; δ = 0.9, max_steps = 200, ε = 1e-6)
    Δz = zpred .- zprev
    # index of variable that changed the most
    # The `iszero` clause is just to keep the first iteration at given `p0`
    i = all(iszero, Δz) ? length(Δz) : argmax(Δz .^ 2)
    c = 0
    zⱼ = zpred
    zⱼ₊₁ = newton_step!(zⱼ, zpred, i, f, J, δ)
    while norm(zⱼ₊₁ - zⱼ) > ε
        zⱼ = zⱼ₊₁
        zⱼ₊₁ = newton_step!(zⱼ, zpred, i, f, J, δ)
        c += 1
        if c > max_steps
            @warn("Newton did not converge.")
            return (zⱼ₊₁, false)
        end
    end
    return zⱼ₊₁, true
end

function newton_step!(zⱼ, zpred, i, f, J, δ)
    Jfinal = mixed_jacobian(zⱼ, i, f, J)
    xⱼ = zⱼ[1:end-1]; pⱼ = zⱼ[end]
    g = f(xⱼ, pⱼ)
    gz = vcat(g, zⱼ[i] - zpred[i])
    zⱼ₊₁ = zⱼ - δ*(inv(Jfinal))*gz
    return zⱼ₊₁
end

function mixed_jacobian(z, i, f, J)
    x = z[1:end-1]; p = z[end]
    # start creating the mixed space jacobian
    j = J(x, p)
    # to the state space jacobian add one more column, derivative towards p
    pder = ForwardDiff.derivative(p -> f(x, p), p)
    Jmixed = hcat(j, pder)
    # add the last row, which is 1 for the `i` entry, 0 everywhere else
    last_row = setindex((@SVector zeros(length(z))), 1.0, i)
    Jfinal = vcat(Jmixed, last_row')
    return Jfinal
end

# Predictor: Secant (linear extrapolation of already found bifurcation curve)
function predictor(zs, dz0)
    if length(zs) == 1
        return zs[end]
    elseif length(zs) == 2 # 1 entry is z0, 2nd entry is 1st found fixed point
        return zs[end] .+ dz0
    else
        return 2zs[end] .- zs[end-1]
    end
end

# Continuation function: perform a step of predictor-corrector and save values
function continuation!(zs, f, J; dz0, pmin, pmax)
    zpred = predictor(zs, dz0)
    (pmin ≤ zpred[end] ≤ pmax) || return false
    zˣ, success = corrector(zpred, zs[end], f, J)
    push!(zs, zˣ)
    return success
end

# Continuation loop: do continuation for a given amount of steps
function continuation(f, J, x0, p0;
        pmin, pmax, dp0, dx0, N = 1000
    )

    z0 = vcat(x0, p0); zs = [z0]; dz0 = vcat(dx0, dp0)

    ps = [p0]
    xs = Dataset([x0])
    stability = Bool[]
    for i in 1:N
        success = continuation!(zs, f, J; dz0, pmin, pmax)
        # Stop iteration if we exceed given parameter margins
        success || break
        # Detect stability of found fixed point (needs `Array` coz of StaticArrays.jl)
        eigenvalues = eigvals(Array(J(zs[end][1:end-1], zs[end][end])))
        μ = maximum(real(v) for v ∈ eigenvalues)
        isstable = μ < 0
        push!(stability, isstable)
    end
    xs = Dataset([z[1:end-1] for z in zs])
    ps = [z[end] for z in zs]
    popfirst!(xs.data); popfirst!(ps) # remove initial guess
    return xs, ps, stability
end

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
export
    SolverOptions

abstract type AbstractSolverOptions{T<:Real} end
abstract type DirectSolverOptions{T} <: AbstractSolverOptions{T} end

"""
    set_options!(opts::AbstractSolverOptions; kwargs...)
    set_options!(solver::AbstractSolver; opts...)

Set solver options via keyword arguments, supporting either solvers or solver option
types directly. Will set any and all options that match the provided arguments. For
example, `set_options!(solver::ALTROSolver, constraint_tolerance=1e-4)` will set the
constraint tolerance option in the ALTRO, Augmented Lagrangian, and Project Newton
solvers.

The only exeption is the `verbose` setting, which always accepts a boolean, while
ALTRO and Augmented Lagrangian solvers accept integers 0-2, with 1 providing output
for the outer AL iterations but not the iLQR iterations.
"""
function set_options!(opts::OPT; d...) where OPT <: AbstractSolverOptions
    for (key, val) in pairs(d)
        if hasfield(OPT, key)
            val = convert(fieldtype(OPT, key), val)
            setfield!(opts, key, val)
        end
    end
end

@with_kw mutable struct SolverOptions{T} <: AbstractSolverOptions{T}
    # iLQR
    ilqr_max_iterations::Int = 300
    dJ_counter_limit::Int = 10
    line_search_lower_bound::T = 1e-8
    line_search_upper_bound::T = 10.0
    iterations_linesearch::Int = 10
    max_cost_value::T = 1.0e8
    max_state_value::T = 1.0e8
    max_control_value::T = 1.0e8
    ilqr_ctol::T = 1e-1
    ilqr_gtol::T = 1e-4
    # backward pass regularization
    bp_reg::Bool = false
    bp_reg_initial::T = 0.0
    bp_reg_increase_factor::T = 1.6
    bp_reg_max::T = 1.0e8
    bp_reg_min::T = 1.0e-8
    bp_reg_type::Symbol = :control
    bp_reg_fp::T = 10.0

    # AL
    al_max_iterations::Int = 30
    kickout_max_penalty::Bool = false
    al_vtol::T = 1e-4
    penalty_max::T = 1e8

    # PN
    projected_newton::Bool = true
    verbose_pn::Bool = false
    n_steps::Int = 2
    solve_type::Symbol = :feasible
    pn_vtol::T = 1e-8
    multiplier_projection::Bool = true
    ρ_chol::T = 1e-2     # cholesky factorization regularization
    ρ_primal::T = 1.0e-8 # primal regularization
    ρ_dual::T = 1.0e-8   # regularization for multiplier projection
    r_threshold::T = 1.1

    # general options
    iterations::Int = 10000   # max number of iterations
    show_summary::Bool = false
    verbose::Int = 0
end

@with_kw mutable struct SolverStats{T}
    # Iteration counts
    iterations::Int = 0
    iterations_outer::Int = 0
    iterations_pn::Int = 0
    # Iteration stats
    iteration::Vector{Int} = Int[]
    iteration_outer::Vector{Int} = Int[]
    iteration_pn::Vector{Bool} = Bool[]   # is the iteration a pn iteration
    cost::Vector{T} = Float64[]
    dJ::Vector{T} = Float64[]
    c_max::Vector{T} = Float64[]
    gradient::Vector{T} = Float64[]
    penalty_max::Vector{T} = Float64[]
    # iLQR
    dJ_zero_counter::Int = 0
    # Other
    tstart::Float64 = time()
    tsolve::Float64 = Inf
    status::TerminationStatus = UNSOLVED
    is_reset::Bool = false
    "Which solver is the top-level solver and responsible for resetting and trimming."
    parent::Symbol
end

function reset!(stats::SolverStats, N::Int, parent::Symbol)
    if parent == stats.parent
        stats.is_reset = false  # force a reset
        reset!(stats, N)
    end
    return nothing
end

function reset!(stats::SolverStats, N::Int=0)
    stats.is_reset && return nothing
    stats.iterations = 0
    stats.iterations_outer = 0
    stats.iterations_pn = 0
    function reset!(v::AbstractVector{T}, N::Int) where T
        resize!(v, N)
        v .= zero(T)
    end
    reset!(stats.iteration, N)
    reset!(stats.iteration_outer, N)
    reset!(stats.iteration_pn, N)
    reset!(stats.cost, N)
    reset!(stats.dJ, N)
    reset!(stats.c_max, N)
    reset!(stats.gradient, N)
    reset!(stats.penalty_max, N)
    stats.tstart = time_ns()
    stats.tsolve = Inf
    stats.dJ_zero_counter = 0
    stats.is_reset = true
    stats.status = UNSOLVED
    return nothing
end

function record_iteration!(stats::SolverStats;
        cost=Inf, dJ=Inf, c_max=Inf, gradient=Inf, penalty_max=Inf,
        is_pn::Bool=false, is_outer::Bool=false
    )
    # don't increment the iteration for an outer loop
    if is_outer
        stats.iterations_outer += 1
    else
        stats.iterations += 1
    end
    i = stats.iterations
    function record!(vec, val, i)
        if val == Inf
            val = i > 1 ? vec[i-1] : val
        end
        vec[i] = val
    end
    stats.iteration[i] = i
    stats.iteration_outer[i] = stats.iterations_outer
    stats.iteration_pn[i] = is_pn
    record!(stats.cost, cost, i)
    record!(stats.dJ, dJ, i)
    record!(stats.c_max, c_max, i)
    record!(stats.gradient, gradient, i)
    record!(stats.penalty_max, penalty_max, i)
    is_pn && (stats.iterations_pn += 1)
    return stats.iterations
end

function trim!(stats::SolverStats, parent::Symbol)
    if parent == stats.parent
        trim!(stats)
    end
end
function trim!(stats::SolverStats)
    N = stats.iterations
    resize!(stats.iteration, N)
    resize!(stats.iteration_outer, N)
    resize!(stats.iteration_pn, N)
    resize!(stats.cost, N)
    resize!(stats.dJ, N)
    resize!(stats.c_max, N)
    resize!(stats.gradient, N)
    resize!(stats.penalty_max, N)
    stats.is_reset = false
    return N
end

function Dict(stats::SolverStats)
    fields = (:iteration, :iteration_outer, :iteration_pn,
        :cost, :dJ, :c_max, :gradient, :penalty_max)
    pairs = [field=>getfield(stats,field) for field in fields]
    Dict(pairs...)
end


text:
function f(x,
y) end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
struct PoincareInvariant1st{ET,DT,TT,ΘT}
    equ::ET
    Θ::ΘT
    Δt::TT
    nloop::Int
    ntime::Int
    nsave::Int
    nt::Int
    I::OffsetArray{DT,1,Vector{DT}}
    J::OffsetArray{DT,1,Vector{DT}}
    L::OffsetArray{DT,1,Vector{DT}}
end

function PoincareInvariant1st(f_equ::Function, f_loop::Function, Θ::ΘT, Δt::TT, d::Int, nloop::Int, ntime::Int, nsave::Int=1, DT=Float64) where {TT,ΘT}

    if get_config(:verbosity) > 1
        println()
        println("First Euler-Poincaré Integral Invariant")
        println("=======================================")
        println()
        println(" nloop = ", nloop)
        println(" ntime = ", ntime)
        println(" nsave = ", nsave)
        println(" Δt    = ", Δt)
        println()
    end

    # compute initial conditions
    q₀ = [f_loop(i/nloop) for i in 1:nloop]

    # initialise euation
    equ = f_equ(q₀)

    # create arrays for results
    nt = div(ntime, nsave)

    I = OffsetArray(zeros(DT, nt+1), 0:nt)
    J = OffsetArray(zeros(DT, nt+1), 0:nt)
    L = OffsetArray(zeros(DT, nt+1), 0:nt)

    PoincareInvariant1st{typeof(equ),DT,TT,ΘT}(equ, Θ, Δt, nloop, ntime, nsave, nt, I, J, L)
end

function evaluate_poincare_invariant(pinv::PoincareInvariant1st, sol::Solution)
    p = zeros(size(sol.q[begin],1), size(sol.q,1), size(sol.q,2))
    g = zeros(size(sol.q[begin],1), size(sol.q,1), size(sol.q,2))
    v = zeros(size(sol.q[begin],1), size(sol.q,2))
    γ = zeros(size(sol.q[begin],1), size(sol.q,2))

    compute_one_form(sol.t, sol.q, p, pinv.Θ)

    if isdefined(sol, :λ)
        compute_correction(sol.t, sol.q, sol.λ, g, pinv.equ.g)
    end

    for i in axes(sol.q,1)
        compute_velocity(hcat(sol.q[i,:]...), v)
        pinv.I[i] = compute_loop_integral(hcat(sol.p[i,:]...), v)

        if isdefined(sol, :p)
            pinv.J[i] = compute_loop_integral(hcat(sol.d[i,:]...), v)
        end

        if isdefined(sol, :λ)
            compute_velocity(hcat(sol.λ[i,:]...), γ)
            pinv.L[i] = compute_loop_integral(p[:,i,:] .- pinv.Δt .* g[:,i,:], v .- pinv.Δt .* γ)
        end
    end

    (pinv.I, pinv.J, pinv.L)
end

function write_to_hdf5(pinv::PoincareInvariant1st, sol::Solution, output_file::String)
    # h5open(output_file, isfile(output_file) ? "r+" : "w") do h5
    h5open(output_file, "w") do h5

        write(h5, "t", sol.t)
        write(h5, "I", pinv.I)

        isdefined(sol, :p) ? write(h5, "J", pinv.J) : nothing
        isdefined(sol, :λ) ? write(h5, "L", pinv.L) : nothing

    end
end


text:
using Compat
using BinDeps

libdir_opt = ""
@static if is_windows()
  libdir_opt = Sys.WORD_SIZE==32 ? "32" : ""
end

@static if is_windows()
  # prefer building if requested
  if haskey(ENV, "BUILD_ON_WINDOWS") && ENV["BUILD_ON_WINDOWS"] == "1"
    saved_defaults = deepcopy(BinDeps.defaults)
    empty!(BinDeps.defaults)
    append!(BinDeps.defaults, [BuildProcess])
  end
end

@BinDeps.setup

function find_julia_lib(lib_suffix::AbstractString, julia_base_dir::AbstractString)
  julia_lib = joinpath(julia_base_dir, "lib", "julia", "libjulia.$lib_suffix")
  if !isfile(julia_lib)
    julia_lib = joinpath(julia_base_dir, "lib", "libjulia.$lib_suffix")
  end
  if !isfile(julia_lib)
    julia_lib = joinpath(julia_base_dir, "lib64", "julia", "libjulia.$lib_suffix")
  end
  if !isfile(julia_lib)
    julia_lib = joinpath(julia_base_dir, "lib", "x86_64-linux-gnu", "julia", "libjulia.$lib_suffix")
  end
  return julia_lib
end

# The base library, needed to wrap functions
cxx_wrap = library_dependency("cxx_wrap", aliases=["libcxx_wrap"])

prefix=joinpath(BinDeps.depsdir(cxx_wrap),"usr")
cxx_wrap_srcdir = joinpath(BinDeps.depsdir(cxx_wrap),"src","cxx_wrap")
cxx_wrap_builddir = joinpath(BinDeps.depsdir(cxx_wrap),"builds","cxx_wrap")
lib_prefix = @static is_windows() ? "" : "lib"
lib_suffix = @static is_windows() ? "dll" : (@static is_apple() ? "dylib" : "so")
julia_base_dir = splitdir(JULIA_HOME)[1]
julia_lib = ""
for suff in ["dll", "dll.a", "dylib", "so"]
  julia_lib = find_julia_lib(suff, julia_base_dir)
  if isfile(julia_lib)
    println("Found Julia library at ", julia_lib)
    break
  end
end

if !isfile(julia_lib)
  throw(ErrorException("Could not locate Julia library at $julia_lib"))
end

julia_include_dir = joinpath(julia_base_dir, "include", "julia")
if !isdir(julia_include_dir)  # then we're running directly from build
  julia_base_dir_aux = splitdir(splitdir(JULIA_HOME)[1])[1]  # useful for running-from-build
  julia_include_dir = joinpath(julia_base_dir_aux, "usr", "include" )
  julia_include_dir *= ";" * joinpath(julia_base_dir_aux, "src", "support" )
  julia_include_dir *= ";" * joinpath(julia_base_dir_aux, "src" )
end

makeopts = ["--", "-j", "$(Sys.CPU_CORES+2)"]

# Set generator if on windows
genopt = "Unix Makefiles"
@static if is_windows()
  makeopts = "--"
  if Sys.WORD_SIZE == 64
    genopt = "Visual Studio 14 2015 Win64"
  elseif !(haskey(ENV, "MINGW_CHOST") && ENV["MINGW_CHOST"] == "i686-w64-mingw32")
    genopt = "Visual Studio 14 2015"
  else
    lib_prefix = "lib" #Makefiles on windows do keep the lib prefix
  end
end

# Functions library for testing
example_labels = [:containers, :except, :extended, :functions, :hello, :inheritance, :parametric, :types]
examples = BinDeps.LibraryDependency[]
for l in example_labels
  @eval $l = $(library_dependency(string(l), aliases=["lib"*string(l)]))
  push!(examples, eval(:($l)))
end
examples_srcdir = joinpath(BinDeps.depsdir(functions),"src","examples")
examples_builddir = joinpath(BinDeps.depsdir(functions),"builds","examples")
deps = [cxx_wrap; examples]

cxx_steps = @build_steps begin
  `cmake -G "$genopt" -DCMAKE_INSTALL_PREFIX="$prefix" -DCMAKE_BUILD_TYPE="Release"  -DJULIA_INCLUDE_DIRECTORY="$julia_include_dir" -DJULIA_LIBRARY="$julia_lib" -DLIBDIR_SUFFIX=$libdir_opt $cxx_wrap_srcdir`
  `cmake --build . --config Release --target install $makeopts`
end

example_paths = AbstractString[]
for l in example_labels
  push!(example_paths, joinpath(prefix,"lib$libdir_opt", "$(lib_prefix)$(string(l)).$lib_suffix"))
end

examples_steps = @build_steps begin
  `cmake -G "$genopt" -DCMAKE_INSTALL_PREFIX="$prefix" -DCMAKE_BUILD_TYPE="Release" -DLIBDIR_SUFFIX=$libdir_opt $examples_srcdir`
  `cmake --build . --config Release --target install $makeopts`
end

# If built, always run cmake, in case the code changed
if isdir(cxx_wrap_builddir)
  BinDeps.run(@build_steps begin
    ChangeDirectory(cxx_wrap_builddir)
    cxx_steps
  end)
end
if isdir(examples_builddir)
  BinDeps.run(@build_steps begin
    ChangeDirectory(examples_builddir)
    examples_steps
  end)
end

provides(BuildProcess,
  (@build_steps begin
    CreateDirectory(cxx_wrap_builddir)
    @build_steps begin
      ChangeDirectory(cxx_wrap_builddir)
      FileRule(joinpath(prefix,"lib$libdir_opt", "$(lib_prefix)cxx_wrap.$lib_suffix"),cxx_steps)
    end
  end), cxx_wrap)

provides(BuildProcess,
  (@build_steps begin
    CreateDirectory(examples_builddir)
    @build_steps begin
      ChangeDirectory(examples_builddir)
      FileRule(example_paths, examples_steps)
    end
  end), examples)

provides(Binaries, Dict(URI("https://github.com/barche/CxxWrap.jl/releases/download/v0.1.5/CxxWrap-julia-$(VERSION.major).$(VERSION.minor)-win$(Sys.WORD_SIZE).zip") => deps), os = :Windows)

@BinDeps.install Dict([(:cxx_wrap, :_l_cxx_wrap),
                       (:containers, :_l_containers),
                       (:except, :_l_except),
                       (:extended, :_l_extended),
                       (:functions, :_l_functions),
                       (:hello, :_l_hello),
                       (:inheritance, :_l_inheritance),
                       (:parametric, :_l_parametric),
                       (:types, :_l_types)])

@static if is_windows()
  if haskey(ENV, "BUILD_ON_WINDOWS") && ENV["BUILD_ON_WINDOWS"] == "1"
    empty!(BinDeps.defaults)
    append!(BinDeps.defaults, saved_defaults)
  end
end


text:
# exp_mult.jl

using MIRT: exp_mult
using MIRT: exp_xform
using Test: @test

L = 10
N = 30
M = 20
A = randn(ComplexF32, N, L)
ur = randn(N)
ui = randn(N)
vr = randn(M)
vi = randn(M)
u = ur + im * ui
v = vr + im * vi

d1 = A' * exp.(-ur * transpose(vr))
@test d1 ≈ exp_mult(A, ur, vr)

d1 = A' * exp.(-ur * transpose(v))
@test d1 ≈ exp_mult(A, ur, v)

d1 = A' * exp.(-u * transpose(vr))
@test d1 ≈ exp_mult(A, u, vr)

d1 = A' * exp.(-u * transpose(v))
@test d1 ≈ exp_mult(A, u, v ; warnboth=false)

d1 = A' * exp.(-ur * transpose(v))
ur = reshape(ur, 1, :)
v = reshape(v, 1, :)
d2 = exp_xform(transpose(A'), ur, v, mode = :matrix)
d2 = transpose(d2)
@test d1 ≈ d2


text:
"""
An instance of the budgeted m-set problem. The budget is a minimum total
weight that must be taken for a solution to be feasible.

All weights are supposed to be positive integers, the budget too (this
hypothesis is crucial to get efficient algorithms). Values are not restricted.

It can be formalised as follows:

``\\max \\sum_i \\mathrm{values}_i x_i``
``\\mathrm{s.t.} \\sum_i x_i \\leq m, \\quad \\sum_i \\mathrm{weights}_i x_i \\geq \\mathrm{budget}, \\qquad x \\in \\{0, 1\\}^d``

`max_weight` is a parameter that may be used by the algorithm and corresponds
to the largest value a single weight may take. If not present, it is computed
as `max(weights)`.

`budget` can be ignored; in that case, it will be computed as the maximum
budget that may have a feasible solution, i.e. `d * max_weight`, `d` being
the dimension of the problem (the number of objects). If a `budget` is given
and the algorithm computes all solutions for lower budgets, then the algorithm
is allowed to stop at the given budget (i.e. it will not be possible to
retrieve a solution for a larger budget).
"""
struct BudgetedMSetInstance
  values::Vector{Float64}
  weights::Vector{Int}
  m::Int
  budget::Int # weights * solution >= budget
  max_weight::Int

  function BudgetedMSetInstance(values::Vector{Float64}, weights::Vector{Int}, m::Int;
                                budget::Union{Nothing, Int}=nothing,
                                max_weight::Union{Nothing, Int}=nothing)
    # Error checking.
    if m < 0
      error("m is less than zero: there is no solution.")
    end

    if m == 0
      error("m is zero: the only solution is to take no items, if the budget constraint is satisfied.")
    end

    if ! isnothing(budget) && budget < 0
      error("Budget is present and negative.")
    end

    if any(weights .< 0)
      error("At least a weight is negative.")
    end

    if length(values) != length(weights)
      error("Not the same number of values and weights; these two vectors must have the same size.")
    end

    # Complete the optional parameters.
    if isnothing(max_weight)
      max_weight = maximum(weights)
    end

    if isnothing(budget)
      d = length(values) # Dimension of the problem.
      budget = d * max_weight
    end

    # Return a new instance.
    new(values, weights, m, budget, max_weight)
  end
end

values(i::BudgetedMSetInstance) = i.values
weights(i::BudgetedMSetInstance) = i.weights
m(i::BudgetedMSetInstance) = i.m
budget(i::BudgetedMSetInstance) = i.budget
max_weight(i::BudgetedMSetInstance) = i.max_weight

value(i::BudgetedMSetInstance, o::Int) = values(i)[o]
values(i::BudgetedMSetInstance, o) = values(i)[o]
weight(i::BudgetedMSetInstance, o::Int) = weights(i)[o]
weights(i::BudgetedMSetInstance, o) = weights(i)[o]
dimension(i::BudgetedMSetInstance) = length(values(i))

struct BudgetedMSetSolution
  instance::BudgetedMSetInstance
  items::Vector{Int} # Indices to the chosen items.
  state::Array{Float64, 3} # Data structure built by the dynamic-programming recursion.
  solutions::Dict{Tuple{Int, Int, Int}, Vector{Int}} # From the indices of state to the corresponding solution.
end

function value(s::BudgetedMSetSolution)
  return sum(s.instance.values[i] for i in s.items)
end

function items(s::BudgetedMSetSolution, budget::Int)
  return s.solutions[s.instance.m, 0, budget]
end

function items_all_budgets(s::BudgetedMSetSolution, max_budget::Int)
  sol = Dict{Int, Vector{Int}}()
  m = s.instance.m
  for budget in 0:max_budget
    sol[budget] = s.solutions[m, 0, budget]
  end
  return sol
end

function value(s::BudgetedMSetSolution, budget::Int)
  its = items(s, budget)
  if -1 in its
    return -Inf
  end
  return sum(s.instance.values[i] for i in its)
end

function budgeted_msets_dp(i::BudgetedMSetInstance)
  # Consider items above δ.
  # Recurrence relation:
  #  V[m, t, δ] = max { V[m, t, δ + 1],  v_δ + V[m - 1, t - w_δ, δ + 1] }
  #                     \_ no item δ _/ \________ take item δ ________/
  V = Array{Float64, 3}(undef, m(i), dimension(i) + 1, budget(i) + 1)
  S = Dict{Tuple{Int, Int, Int}, Vector{Int}}()

  # Initialise: µ == 1, just take the best element among [δ+1, d], as long as
  # the budget constraint is satisfied; δ == d, take nothing, whatever the
  # budget; β == 0, a simple m-set problem.
  for µ in 1:m(i) # β == 0
    for δ in (dimension(i) - 1):-1:0 # δ < d
      # Due to δ, there may be too few items (with respect to µ).
      # i.e., there are d - δ items to consider.
      m_ = min(dimension(i) - δ, µ)

      items = collect(partialsortperm(values(i, (δ + 1):dimension(i)), 1:m_, rev=true)) .+ δ
      V[µ, δ + 1, 0 + 1] = sum(value(i, o) for o in items)
      S[µ, δ, 0] = items
    end
  end

  for β in 0:budget(i)
    # δ == d
    for µ in 1:m(i)
      if β == 0
        # Optimum solution: take no object.
        V[µ, dimension(i) + 1, β + 1] = 0.0
        S[µ, dimension(i), β] = Int[]
      else
        # No solution.
        V[µ, dimension(i) + 1, β + 1] = -Inf
        S[µ, dimension(i), β] = [-1] # Nonsensical solution.
      end
    end

    # δ < d
    for δ in (dimension(i) - 1):-1:0
      # Filter the available objects (i.e. above δ+1 and satisfying the
      # budget constraint) and retrieve just their indices.
      all_objects = collect(enumerate(weights(i)))
      all_objects = collect(filter(o -> o[1] >= δ + 1, all_objects))
      all_objects = collect(filter(o -> o[2] >= β, all_objects))
      all_objects = collect(t[1] for t in all_objects)

      # No feasible solution?
      if length(all_objects) == 0
        V[1, δ + 1, β + 1] = -Inf
        S[1, δ, β] = [-1] # Nonsensical solution.
        continue
      end

      # Take the best one available object.
      v, x = findmax(collect(value(i, t) for t in all_objects))
      x = all_objects[x]

      V[1, δ + 1, β + 1] = v
      S[1, δ, β] = [x]
    end
  end

  # Dynamic part.
  for β in 1:budget(i)
    for µ in 2:m(i)
      for δ in (dimension(i) - 1):-1:0
        remaining_budget_with_δ = max(0, β - weights(i)[δ + 1])
        take_δ = value(i, δ + 1) + V[µ - 1, δ + 1 + 1, remaining_budget_with_δ + 1]
        dont_take_δ = V[µ, δ + 1 + 1, β + 1]

        both_subproblems_infeasible = -Inf == take_δ && -Inf == dont_take_δ

        if both_subproblems_infeasible
          V[µ, δ + 1, β + 1] = -Inf
          S[µ, δ, β] = [-1]
        elseif take_δ >= dont_take_δ
          V[µ, δ + 1, β + 1] = take_δ
          S[µ, δ, β] = vcat(δ + 1, S[µ - 1, δ + 1, remaining_budget_with_δ])
        else
          V[µ, δ + 1, β + 1] = dont_take_δ
          S[µ, δ, β] = S[µ, δ + 1, β]
        end
      end
    end
  end

  return BudgetedMSetSolution(i, S[m(i), 0, budget(i)], V, S)
end

function _budgeted_msets_lp_sub(i::BudgetedMSetInstance, solver)
  model = Model(solver)
  @variable(model, x[1:length(values(i))], Bin)
  @objective(model, Max, dot(x, values(i)))
  @constraint(model, sum(x) <= m(i))
  @constraint(model, c, dot(x, weights(i)) >= 0)

  set_silent(model)

  return model, x, c
end

function budgeted_msets_lp(i::BudgetedMSetInstance; solver=nothing, β::Int=budget(i))
  # Solve for all budgets at once, even though it is not really more efficient than looping outside this function.
  model, x, c = _budgeted_msets_lp_sub(i, solver)
  set_normalized_rhs(c, β)
  optimize!(model)

  V = Array{Float64, 3}(undef, m(i), dimension(i) + 1, budget(i) + 1)
  S = Dict{Tuple{Int, Int, Int}, Vector{Int}}()

  if termination_status(model) == MOI.OPTIMAL
    sol = findall(JuMP.value.(x) .>= 0.5)
    V[m(i), 0 + 1, β + 1] = objective_value(model)
    S[m(i), 0, β] = sol
  else
    V[m(i), 0 + 1, β + 1] = -Inf
    S[m(i), 0, β] = Int[-1]
  end

  return BudgetedMSetSolution(i, S[m(i), 0, β], V, S)
end

function budgeted_msets_lp_select(i::BudgetedMSetInstance, budgets; solver=nothing)
  model, x, c = _budgeted_msets_lp_sub(i, solver)

  V = Array{Float64, 3}(undef, m(i), dimension(i) + 1, budget(i) + 1)
  S = Dict{Tuple{Int, Int, Int}, Vector{Int}}()
  for budget in budgets
    set_normalized_rhs(c, budget)
    optimize!(model)

    if termination_status(model) == MOI.OPTIMAL
      sol = findall(JuMP.value.(x) .>= 0.5)
      V[m(i), 0 + 1, budget + 1] = objective_value(model)
      S[m(i), 0, budget] = sol
    else
      V[m(i), 0 + 1, budget + 1] = -Inf
      S[m(i), 0, budget] = Int[-1]
    end
  end

  return BudgetedMSetSolution(i, S[m(i), 0, maximum(budgets)], V, S)
end

function budgeted_msets_lp_all(i::BudgetedMSetInstance; solver=nothing, max_budget::Int=budget(i))
  # Solve for all budgets at once, even though it is not really more efficient than looping outside this function.
  return budgeted_msets_lp_select(i, 0:max_budget, solver=solver)
end


text:
export AbstractSpatialEmbedding
export SpatioTemporalEmbedding
export outdim
export AbstractBoundaryCondition, PeriodicBoundary, ConstantBoundary
export indices_within_sphere
export light_cone_embedding, cubic_shell_embedding

#####################################################################################
#                          Spatio Temporal Delay Embedding                          #
#####################################################################################
"""
    AbstractSpatialEmbedding <: AbstractEmbedding
Super-type of spatiotemporal embedding methods. Valid subtypes:
* `SpatioTemporalEmbedding`
* `PCAEmbedding`
"""
abstract type AbstractSpatialEmbedding{Φ,BC,X} <: AbstractEmbedding end
const ASE = AbstractSpatialEmbedding

"""
    AbstractBoundaryCondition
Super-type of boundary conditions for [`SpatioTemporalEmbedding`](@ref).
Use `subtypes(AbstractBoundaryCondition)` for available methods.
"""
abstract type AbstractBoundaryCondition end

"""
	ConstantBoundary(b) <: AbstractBoundaryCondition
Constant boundary condition type. Enforces constant boundary conditions
when passed to [`SpatioTemporalEmbedding`](@ref)
by filling missing out-of-bounds values in the reconstruction with
parameter `b`.
"""
struct ConstantBoundary{T} <: AbstractBoundaryCondition
    b::T
end

"""
	PeriodicBoundary <: AbstractBoundaryCondition
Periodic boundary condition struct. Enforces periodic boundary conditions
when passed to [`SpatioTemporalEmbedding`](@ref) in the reconstruction.
"""
struct PeriodicBoundary <: AbstractBoundaryCondition end

"""
	Region{Φ}
Internal struct for efficiently keeping track of region far from boundaries of field.
Used to speed up reconstruction process.
"""
struct Region{Φ}
	mini::NTuple{Φ,Int}
	maxi::NTuple{Φ,Int}
end

Base.length(r::Region{Φ}) where Φ = prod(r.maxi .- r.mini .+1)
function Base.in(idx, r::Region{Φ}) where Φ
	for φ=1:Φ
		r.mini[φ] <= idx[φ] <= r.maxi[φ] || return false
 	end
 	return true
end
Base.CartesianIndices(r::Region{Φ}) where Φ =
	 CartesianIndices{Φ,NTuple{Φ,UnitRange{Int}}}(([r.mini[φ]:r.maxi[φ] for φ=1:Φ]...,))

function inner_region(βs::Vector{CartesianIndex{Φ}}, fsize) where Φ
	mini = Int[]
	maxi = Int[]
	for φ = 1:Φ
		js = map(β -> β[φ], βs) # jth entries
		mi,ma = extrema(js)
		push!(mini, 1 - min(mi, 0))
		push!(maxi,fsize[φ] - max(ma, 0))
	end
	return Region{Φ}((mini...,), (maxi...,))
end

function project_inside(α::CartesianIndex{Φ}, r::Region{Φ}) where Φ
	CartesianIndex(mod.(α.I .-1, r.maxi).+1)
end

"""
	SpatioTemporalEmbedding{Φ,BC,X} → embedding
A spatio temporal delay coordinates structure to be used as a functor. Applies
to data of `Φ` spatial dimensions and gives an embedding of dimensionality `X`.

	embedding(rvec, s, t, α)
Operates inplace on `rvec` (of length `X`) and reconstructs vector from spatial
timeseries `s` at timestep `t` and cartesian index `α`.
Note that there are no bounds checks for `t`.

It is assumed that `s` is a `Vector{<:AbstractArray{T,Φ}}`.

## Constructors
There are some convenience constructors that return intuitive embeddings here:
* [`cubic_shell_embedding`](@ref)
* [`light_cone_embedding`](@ref)

The "main" constructor is

	SpatioTemporalEmbedding{X}(τ, β, bc, fsize)

which allows full control over the spatio-temporal embedding.
* `Χ == length(τ) == length(β)` : dimensionality of resulting reconstructed space.
* `τ::Vector{Int}` = Vector of temporal delays **for each entry** of the reconstructed space
  (sorted in ascending order).
* `β::Vector{CartesianIndex{Φ}}` = vector of *relative* indices of spatial delays
  *for each entry* of the reconstructed space.
* `bc::BC` : boundary condition.
* `fsize::NTuple{Φ, Int}` : Size of each state in the timeseries.

For example, if you want to have spatial neighbors [0, ±1] for time delays
τ = 2, 4, then
```
τ = [2,2,2,4,4,4]
β = CartesianIndex.([1, 0, -1, 1, 0, -1])
```
"""
struct SpatioTemporalEmbedding{Φ,BC,X} <: AbstractSpatialEmbedding{Φ,BC,X}
  	τ::Vector{Int}
	β::Vector{CartesianIndex{Φ}}
	inner::Region{Φ}  #inner field far from boundary
	whole::Region{Φ}
    boundary::BC

	function SpatioTemporalEmbedding{X}(
			τ::Vector{Int}, β::Vector{CartesianIndex{Φ}}, bc::BC, fsize::NTuple{Φ, Int}
			) where {Φ,BC,X}
        @assert issorted(τ) "Delays need to be sorted in ascending order"
		#"ConstantBoundary condition value C needs to be the same type as values in s"
		inner = inner_region(β, fsize)
		whole = Region((ones(Int,Φ)...,), fsize)
		return new{Φ,BC,X}(τ,β,inner,whole, bc)
	end
end
const STE = SpatioTemporalEmbedding

#This function is not safe. If you call it directly with bad params - can fail
function (r::SpatioTemporalEmbedding{Φ,ConstantBoundary{T},X})(rvec,s,t,α) where {T,Φ,X}
	if α in r.inner
		@inbounds for n=1:X
			rvec[n] = s[ t + r.τ[n] ][ α + r.β[n] ]
		end
	else
		@inbounds for n=1:X
			rvec[n] = α + r.β[n] in r.whole ? s[ t+r.τ[n] ][ α+r.β[n] ] : r.boundary.b
		end
	end
	return nothing
end

function (r::SpatioTemporalEmbedding{Φ,PeriodicBoundary,X})(rvec,s,t,α) where {Φ,X}
	if α in r.inner
		@inbounds for n=1:X
			rvec[n] = s[ t + r.τ[n] ][ α + r.β[n] ]
		end
	else
		@inbounds for n=1:X
			rvec[n] = s[ t + r.τ[n] ][ project_inside(α + r.β[n], r.whole) ]
		end
	end
	return nothing
end

get_num_pt(em::AbstractSpatialEmbedding) = prod(em.whole.maxi)
get_τmax(em::AbstractSpatialEmbedding{Φ,BC,X}) where {Φ,BC,X} = em.τ[X]
outdim(em::AbstractSpatialEmbedding{Φ,BC,X}) where {Φ,BC,X} = X

get_usable_idxs(em::AbstractSpatialEmbedding{Φ,PeriodicBoundary,X}) where {Φ,X} =
			CartesianIndices(em.whole)
get_usable_idxs(em::AbstractSpatialEmbedding{Φ,<:ConstantBoundary,X}) where {Φ,X} =
			CartesianIndices(em.inner)

function Base.show(io::IO, em::SpatioTemporalEmbedding{Φ,BC, X}) where {Φ,BC,X}
	print(io, "$(Φ)D Spatio-Temporal Delay Embedding with $X Entries")
    if BC == PeriodicBoundary
        println(io, " and PeriodicBoundary condition.")
    else
        println(io, " and ConstantBoundary condition with c = $(em.boundary.c).")
    end
    println(io, "The included neighboring points are (forward embedding):")
    for (τ,β) in zip(em.τ,em.β)
        println(io, "τ = $τ , β = $(β.I)")
    end
end

#####################################################################################
#                                   CONSTRUCTORS                                    #
#####################################################################################
"""
    cubic_shell_embedding(s, γ, τ, B, k, bc=PeriodicBoundary()) → embedding
Create a [`SpatioTemporalEmbedding`](@ref) instance that
includes spatial neighbors in hypercubic *shells*.
The embedding is to be used with data from `s`.

## Description
Points are participating in the embedding by forming hypercubic shells around
the current point. The total shells formed are `B`. The points on the shells have
spatial distance `k ≥ 1` (distance in indices, like a cityblock metric).
`k = 1` means that all points of the shell participate.
The points of the hypercubic grid can be separated
by `k ≥ 1` points apart (i.e. dropping `k-1` in-between points).
In short, in each spatial dimension of the system the cartesian
offset indices are `-B*k : k : k*B`.

`γ` is the number of temporal steps in the past to be included in the embedding,
where each step in the past has additional delay time `τ::Int`.
`γ=0` corresponds to using only the present. Notice that **all** embedded
time frames have the same spatial structure, in contrast to [`light_cone_embedding`](@ref).

As an example, consider one of the `γ` embedded frames (all are the same) of a system
with 2 spatial dimensions (`□` = current point, (included *by definition* in the
embedding), `n` = included points in the embedding coming from `n`-th shell,
`.` = points not included in the embedding)

```
      B = 2,  k = 1        |        B = 1,  k = 2        |        B = 2,  k = 2
                           |                             |
.  .  .  .  .  .  .  .  .  |  .  .  .  .  .  .  .  .  .  |  2  .  2  .  2  .  2  .  2
.  .  .  .  .  .  .  .  .  |  .  .  .  .  .  .  .  .  .  |  .  .  .  .  .  .  .  .  .
.  .  2  2  2  2  2  .  .  |  .  .  1  .  1  .  1  .  .  |  2  .  1  .  1  .  1  .  2
.  .  2  1  1  1  2  .  .  |  .  .  .  .  .  .  .  .  .  |  .  .  .  .  .  .  .  .  .
.  .  2  1  □  1  2  .  .  |  .  .  1  .  □  .  1  .  .  |  2  .  1  .  □  .  1  .  2
.  .  2  1  1  1  2  .  .  |  .  .  .  .  .  .  .  .  .  |  .  .  .  .  .  .  .  .  .
.  .  2  2  2  2  2  .  .  |  .  .  1  .  1  .  1  .  .  |  2  .  1  .  1  .  1  .  2
.  .  .  .  .  .  .  .  .  |  .  .  .  .  .  .  .  .  .  |  .  .  .  .  .  .  .  .  .
.  .  .  .  .  .  .  .  .  |  .  .  .  .  .  .  .  .  .  |  2  .  2  .  2  .  2  .  2
```
"""
function cubic_shell_embedding(
		s::AbstractArray{<:AbstractArray{T,Φ}},
		γ, τ, B, k, boundary::BC = PeriodicBoundary()
		) where {T,Φ, BC<:AbstractBoundaryCondition}
    if (BC <: ConstantBoundary) && typeof(boundary.b) != T
	     throw(ArgumentError(
		"Boundary value must be same element type as the timeseries data."))
	end
    @assert k ≥ 1
	X = (γ+1)*(2B+1)^Φ
	τs = Vector{Int}(undef,X)
	βs = Vector{CartesianIndex{Φ}}(undef,X)
	n = 1
	for d=0:γ, α = Iterators.product([-B*k:k:B*k for φ=1:Φ]...)
		τs[n] = d*τ
		βs[n] = CartesianIndex(α)
		n +=1
	end
	return SpatioTemporalEmbedding{X}(τs, βs, boundary, size(s[1]))
end

"""
    indices_within_sphere(r, Φ) → β
Return all cartesian indices within a hypersphere or radius `r` and dimension `Φ`.
"""
function indices_within_sphere(radius, dimension)
    #Floor radius to nearest integer. Does not lose points!
    r = floor(Int,radius)
    #Hypercube of indices
    hypercube = CartesianIndices((repeat([-r:r], dimension)...,))
    #Select subset of hc which is in Hypersphere
    βs = [ β for β ∈ hypercube if norm(β.I) <= radius ]
end

"""
    light_cone_embedding(s, γ, τ, r0, c, bc=PeriodicBoundary()) → embedding
Create a [`SpatioTemporalEmbedding`](@ref) instance that
includes spatial and temporal neighbors of a point based on the notion of
a *light cone*.

The embedding is to be used with data from `s`.

## Description
Information does not travel instantly but with some finite speed `c ≥ 0.0`.
This constructor creates a cone-like embedding including all points in
space and time, whose value can influence a prediction based on the
information speed `c`. `γ` is the number of temporal steps in the past to be
included in the embedding, where each step in the past has additional delay time `τ::Int`.
`γ=0` corresponds to using only the present. `r0` is the initial radius at the top of
the cone, i.e. the radius of influence at the present. `bc` is the boundary condition.

The radius of the light cone evolves as: `r_i = i*τ*c + r0` for each step `i ∈ 0:γ`.

As an example, in a one-dimensional system with `γ = 1, τ = 2, r0 = 1`,
the embedding looks like (`□` = current point (included *by definition* in the embedding),
`o` point to be predicted using
[`temporalprediction`](@ref), `x` = points included in the embedding,
`.` = points not included in the embedding)

```
time  | c = 1.0               | c = 2.0               | c = 0.0

n + 1 | ..........o.......... | ..........o.......... | ..........o..........
n     | .........x□x......... | .........x□x......... | .........x□x.........
n - 1 | ..................... | ..................... | .....................
n - τ | .......xxxxxxx....... | .....xxxxxxxxxx...... | .........xxx.........
```
Besides this example, in the official documentation we show a function `explain_light_cone`
which produces a plot of the light cone for 2 spatial dimensions
(great for understanding!).
"""
function light_cone_embedding(
    s::AbstractArray{<:AbstractArray{T,Φ}},
    γ,
    τ,
    r0,
    c,
    bc::BC = PeriodicBoundary()
    ) where {T,Φ, BC<:AbstractBoundaryCondition}
    if (BC <: ConstantBoundary) && typeof(bc.b) != T
        throw(ArgumentError(
        "Boundary value must be same element type as the timeseries data."))
    end
    τs = Int[]
    βs = CartesianIndex{Φ}[]
    maxτ = τ*γ
    for delay = τ*(γ:-1:0) # Backwards for forward embedding
        radius = max(0, c*delay + r0)
        β = indices_within_sphere(radius, Φ)
        append!(βs, β)
        append!(τs, fill(maxτ - delay, length(β)))
    end
    X = length(τs)
    return SpatioTemporalEmbedding{X}(τs, βs, bc, size(s[1]))
end

SpatioTemporalEmbedding(s, p::NamedTuple{(:γ, :τ, :B, :k, :bc)}) =
    cubic_shell_embedding(s, p.γ, p.τ, p.B, p.k, p.bc)

SpatioTemporalEmbedding(s, p::NamedTuple{(:γ, :τ, :r, :c, :bc)}) =
    light_cone_embedding(s, p.γ, p.τ, p.r, p.c, p.bc)

SpatioTemporalEmbedding(s, p::NamedTuple{(:γ, :τ, :B, :k)}) =
    cubic_shell_embedding(s, p.γ, p.τ, p.B, p.k, PeriodicBoundary())

SpatioTemporalEmbedding(s, p::NamedTuple{(:γ, :τ, :r, :c)}) =
    light_cone_embedding(s, p.γ, p.τ, p.r, p.c, PeriodicBoundary())

SpatioTemporalEmbedding(s, p::NamedTuple{(:γ, :τ, :r0, :c)}) =
    light_cone_embedding(s, p.γ, p.τ, p.r0, p.c, PeriodicBoundary())

Base.:(==)(em1::T, em2::T) where {T <: AbstractSpatialEmbedding} =
    all(( eval(:($em1.$name == $em2.$name)) for name ∈ fieldnames(T)))


text:
module Utils

using Parameters, DifferentialEquations

export data, convert_functor_model, data, screen_domain, rangeStepsize

@with_kw mutable struct data
    x
    y
end

@with_kw struct MOD
    mod
end
function (set::MOD)(x, p)
    return @. set.mod(p...)(x)
end

function convert_functor_model(model_functor)
    return MOD(model_functor)
end

"Range with fixed step size"
rangeStepsize(start , stop ,resolution) = LinRange(start, stop, Int(abs(stop-start) ÷ resolution))

function screen_domainSS(f, domain; resolution=0.01, reset_step=100.0)
    x_lower, x_upper = domain[1], domain[2]

    if resolution> abs(x_upper-x_lower)
        res_old = resolution
        resolution = abs(x_upper-x_lower) / reset_step
        @warn "The resolution ($res_old) larger than the domain range($domain). Reset resolution to $(resolution)"
    end

    x_data = rangeStepsize(x_lower, x_upper, resolution)
    y_data = f.(x_data)

    return data(x=x_data, y=y_data)
end

"""Screening Domain with given input parameter"""
function screen_domainSS(mutable_functor, input_par::Symbol, domain, u0; resolution=0.01, method=DynamicSS(Rosenbrock23()))
    xdata = Utils.rangeStepsize(domain[1], domain[2], resolution)
    ydata = []
    for x in xdata
        sol = get_ss(mutable_functor, u0, input_par, method)(x)
        append!(ydata, sol)
    end
    return data(x=xdata, y=ydata)
end

@with_kw struct get_ss
    functor_model
    u0
    input_par
    method= DynamicSS(AutoTsit5(Rosenbrock23()), abstol=1e-8) # Recommended solver for positive_feedback model
end

function (set::get_ss)(x)
    @unpack functor_model, u0, input_par, method = set
    setfield!(functor_model, input_par, x)
    prob = SteadyStateProblem(functor_model, u0)
    sol = solve(prob, method)
    return sol.u
end

end


text:
"""
    withprogress(foldable) -> foldable′

Wrap a foldable so that progress is shown in logging-based progress meter
(e.g., Juno) during [`foldl`](@ref), [`foldxt`](@ref), [`foldxd`](@ref), etc.

For parallel reduction such as `foldxt` and `foldxd`, reasonably small
`basesize` and `threads_basesize` (for `foldxd`) must be used to ensure that
progress information is updated frequently.  However, it may slow down the
computation if `basesize` is too small.

# Keyword Arguments
- `interval::Real`: Minimum interval (in seconds) for how often progress is
  logged.

# Examples
```jldoctest
julia> using Transducers

julia> xf = Map() do x
           sleep(0.01)
           x
       end;

julia> foldl(+, xf, withprogress(1:100; interval=1e-3))  # see progress meter
5050
```

In `foldl` and `foldxt`, `withprogress` can be nested.  This is not supported
in `foldxd`.

```jldoctest; setup = :(using Transducers)
julia> xf = opcompose(
           MapCat() do x
               withprogress(1:x; interval=1e-3)  # nested progress
           end,
           Map() do x
               sleep(0.5)
               x
           end,
       );

julia> if VERSION >= v"1.3-alpha"
           # Calling `sleep` in thread is safe in Julia 1.3:
           foldxt(+, xf, withprogress(1:10; interval=1e-3); basesize=1)
       else
           foldl(+, xf, withprogress(1:10; interval=1e-3))
       end
220
```
"""
withprogress(
    foldable;
    interval::Real = 0.1,
) = ProgressLoggingFoldable(
    foldable,
    convert(Float64, interval),
)

struct ProgressLoggingFoldable{T} <: Foldable
    foldable::T
    interval::Float64
end

Base.IteratorSize(::Type{ProgressLoggingFoldable{T}}) where {T} = Base.IteratorSize(T)
Base.IteratorEltype(::Type{ProgressLoggingFoldable{T}}) where {T} = Base.IteratorEltype(T)
Base.length(foldable::ProgressLoggingFoldable) = length(foldable.foldable)
Base.eltype(::Type{ProgressLoggingFoldable{T}}) where {T} = eltype(T)

# Use Juno/Atom-compatible log-level.  See:
# https://github.com/JunoLab/Atom.jl/blob/v0.11.1/src/progress.jl#L75-L76
const PROGRESSLEVEL = LogLevel(-1)

# Juno.progress
function __progress(f; name = "")
    _id = gensym()
    @logmsg PROGRESSLEVEL name progress=0.0 _id=_id
    try
        f(_id)
    finally
        @logmsg PROGRESSLEVEL name progress="done" _id=_id
    end
end

function __foldl__(rf0, init, coll::ProgressLoggingFoldable)
    __progress() do id
        n = length(coll.foldable)
        progress_interval = coll.interval
        scaninit = (0, time())
        xf = ScanEmit(scaninit) do (i, t0), x
            t1 = time()
            if t1 - t0 > progress_interval
                @logmsg PROGRESSLEVEL "foldl" _id=id progress=i/n
                t0 = t1
            end
            x, (i + 1, t0)
        end
        rf = Reduction(xf, rf0)
        val = wrap(rf, scaninit, init)  # manually `start`ing `ScanEmit`
        return Transducers.__foldl__(rf, val, coll.foldable)
    end
end

Base.eltype(::ProgressLoggingFoldable{T}) where {T} = eltype(T)

struct LogProgressOnCombine{C} <: Transducer
    chan::C
    progress_interval::Float64
end

# Recrod started time:
start(rf::R_{LogProgressOnCombine}, result) =
    wrap(rf, (time(), 0), start(inner(rf), result))

# Count processed number of items:
@inline next(rf::R_{LogProgressOnCombine}, result, input) =
    wrapping(rf, result) do (t0, n), iresult
        (t0, n + 1), next(inner(rf), iresult, input)
    end

@inline complete(rf::R_{LogProgressOnCombine}, result) =
    complete(inner(rf), unwrap(rf, result)[2])
# TODO: `put!` on complete as well

# Send number of processed item to progress channel:
@inline function combine(rf::R_{LogProgressOnCombine}, a, b)
    (ta, na), ira = unwrap(rf, a)
    (tb, nb), irb = unwrap(rf, b)
    irc = combine(inner(rf), ira, irb)

    tc = time()
    t0 = min(ta, tb)
    nc = na + nb
    if tc - t0 > xform(rf).progress_interval
        put!(xform(rf).chan, nc)
        nc = 0
    else
        tc = t0
    end

    return wrap(rf, (tc, nc), irc)
end

function setup_logprogressoncombine(rf0, interval, chan)
    if rf0 isa R_{UseSIMD}
        xf0 = xform(rf0)
        rfinner = inner(rf0)
    else
        xf0 = IdentityTransducer()
        rfinner = rf0
    end

    xf = LogProgressOnCombine(chan, interval) ∘ xf0
    rf = Reduction(xf, rfinner)
    return rf
end

function _reduce_progress(reduce_impl, rf0, init, coll)
    chan = Channel{Int}(0)
    rf = setup_logprogressoncombine(rf0, coll.reducible.interval, chan)
    reducible = @set coll.reducible = coll.reducible.foldable
    progress_task = @async let n = length(coll.reducible.foldable)
        __progress() do id
            foreach(Scan(+), chan) do i
                @logmsg PROGRESSLEVEL "foldxt" _id=id progress=i/n
            end
        end
    end
    try
        result = reduce_impl(rf, init, reducible)
        result isa Reduced && return result
        # Manually unwrap LogProgressOnCombine's private state:
        _, iresult = unwrap(rf isa R_{UseSIMD} ? inner(rf) : rf, result)
        return iresult
    finally
        close(chan)
        wait(progress_task)
    end
end

_reduce(ctx, rf, init, coll::SizedReducible{<:ProgressLoggingFoldable}) =
    _reduce_progress(rf, init, coll) do rf, init, coll
        _reduce(ctx, rf, init, coll)
    end

if VERSION >= v"1.2"
    _reduce_threads_for(rf, init, coll::SizedReducible{<:ProgressLoggingFoldable}) =
        _reduce_progress(_reduce_threads_for, rf, init, coll)
else
    # In earlier versions, Channel was not thread-safe (?)
    # https://github.com/JuliaLang/julia/pull/30186
    _reduce_threads_for(rf, init, coll::SizedReducible{<:ProgressLoggingFoldable}) =
        _reduce_threads_for(
            rf,
            init,
            (@set coll.reducible = coll.reducible.foldable),
        )
end

if VERSION < v"1.3-alpha"
    maybe_collect(coll::ProgressLoggingFoldable) =
        @set coll.foldable = maybe_collect(coll.foldable)
end

struct RemoteReduceWithLogging{C} <: Function
    chan::C
    progress_interval::Float64
end
# Manually create a closure to make it work nicely with Revise.jl:
# https://github.com/timholy/Revise.jl/pull/157

function (foldl::RemoteReduceWithLogging)(rf0, init, coll, basesize)
    chan = _Channel(eltype(foldl.chan), 0) do chan
        while true
            put!(foldl.chan, take!(chan))
        end
    end
    local rf
    result = try
        rf = setup_logprogressoncombine(rf0, foldl.progress_interval, chan)
        _transduce_assoc_nocomplete(rf, init, coll, basesize)
    finally
        close(chan)
    end
    result isa Reduced && return result
    # Manually unwrap LogProgressOnCombine's private state:
    (_t0, n), iresult = unwrap(rf isa R_{UseSIMD} ? inner(rf) : rf, result)
    if n > 0
        put!(foldl.chan, n)
    end
    return iresult
end

function dtransduce(
    xform::Transducer, step, init, coll::ProgressLoggingFoldable;
    kwargs...,
)
    chan = Distributed.RemoteChannel()
    remote_reduce_with_logging = RemoteReduceWithLogging(chan, coll.interval)
    progress_task = @async let n = length(coll.foldable)
        __progress() do id
            i = 0
            while true
                i += try
                    take!(chan)
                catch
                    return
                end
                @logmsg PROGRESSLEVEL "foldxd" _id=id progress=i/n
            end
            #=
            foreach(Scan(+), chan) do i
                @logmsg PROGRESSLEVEL "foldxd" _id=id progress=i/n
            end
            =#
        end
    end

    try
        return dtransduce(
            xform, step, init, coll.foldable;
            _remote_reduce = remote_reduce_with_logging,
            kwargs...,
        )
    finally
        close(chan)
        wait(progress_task)
    end
end


text:
using BHAPtfem

ProjDir = dirname(@__FILE__)
ProjName = split(ProjDir, "/")[end]

segments = [
# Element type,  Material,    Length,     OD,         ID,        fc,     noofelements
  :bit,                    :steel ,        0.00,      2.75,      12.25,     0.0,         0,
  :collar ,              :monel ,    25.00,      2.75,       7.75,      0.0,       25,
  :stabilizer,          :steel,         0.00,       2.75,     12.125,    0.0,         0,
  :collar,                :steel,         30.00,      2.75,       7.75,      0.0,        30,
  :stabilizer,          :steel,          0.00,       2.75,     12.125,    0.0,          0,
  :collar,                :steel,        30.00,       2.75,      7.75,      0.0,         30,
  :stabilizer,          :steel,          0.00,       2.75,     12.125,    0.0,          0,
  :collar,               :steel,        90.00,       2.75,      7.75,      0.0,         90,
  :stabilizer,          :steel,          0.00,       2.75,     12.125,    0.0,          0,
  :pipe,                  :steel,     100.00,       2.75,      7.75,       0.0,       100
];

traj = [
#   Heading,      Bit diameter
        60.0,             12.25
]

wobrange = 45:5:55
inclinationrange = 35:5:55           # Or e.g. incls = [5, 10]

println()

prob = problem(segments, traj, wobrange, inclinationrange,
  medium=:lightmud, computefunction=PtFEM.p44, pdir=ProjDir)

@time results = solve!(prob)

sleep(1) # Wait for all processes to complete

println("\nSize of result  array of tuples = $(size(results))")
println("Size of results[1] in results tuple = $(size(results[3][2]))\n")

println(results[1][2])

#=
Needs to generate something like:

data = Dict(
  # Frame(nels, nn, ndim, nst, nip, finite_element(nod, nodof))
  :struc_el => Frame(10, 11, 2, 1, 1, Line(2, 3)),
  :properties => [
    1.07207e9 4.53116e9;],
  :x_coords =>0.0:1.0:10.0,
  :y_coords => [0.0 for i in 1:11],
  :g_num => [
    1 2 3 4 5 6 7 8 9 10;
    2 3 4 5 6 7 8 9 10 11],
  :support => [
    (1, [1 0 1]),
    (5, [1 0 1]),
    (11, [0, 1, 1])
    ],
  :loaded_nodes => [
    (1, [10000.0 -600.0 -600.0]),
    (2, [0.0 -600.0 0.0]),
    (3, [0.0 -600.0 0.0]),
...
    (10, [0.0 -600.0 0.0]),
    (11, [0.0 -600.0 -600.0])
    ],
  :penalty => 1e19,
  :eq_nodal_forces_and_moments => [
    (1, [0.0 -60.0 -60.0 0.0 -60.0 60.0]),
    (2, [0.0 -120.0 -140.0 0.0 -120.0 140.0]),
    (3, [0.0 -20.0 -6.67 0.0 -20.0 6.67])
  ]
)
  =#