text:
"""
    add_quadratic!(Λ::PiecewiseQuadratic{T}, μ::QuadraticPolynomial{T}) where {T}
Inserts a quadratic polynomial μ into the linked list `Λ`, which represents a piecewise
quadratic polynomial, so that the new `Λ` satisfies
`Λ(y) := min{ Λ(y) ,μ(y) } ∀ y`
"""
function add_quadratic!(Λ::PiecewiseQuadratic{T}, μ::QuadraticPolynomial{T}) where T

    if Λ.next.left_endpoint == Inf # I.e. the piecewise quadratic object is empty, perhaps better to add dummy polynomial
        μ.has_been_used = true
        insert(Λ, μ, -Inf) # FIXME: Is specific insert function needed?
        return
    end

    λ_prev = Λ # Points to the list head which is just a node with NaN data
    λ_curr = Λ.next # The leftmost segment of the list

    while λ_curr.left_endpoint != Inf #???
        #global counter2 += 1
        DEBUG && println(Λ)

        left_endpoint = λ_curr.left_endpoint
        right_endpoint = get_right_endpoint(λ_curr)

        Δa = μ.a - λ_curr.p.a
        Δb = μ.b - λ_curr.p.b
        Δc = μ.c - λ_curr.p.c

        b2_minus_4ac =  Δb^2 - 4*Δa*Δc

        if Δa > 0 # μ has greater curvature, i.e., μ is smallest in the middle if intersect
            if b2_minus_4ac <= ACCURACY
                # Zero (or one) intersections, old quadratic is smallest, just step forward
                DEBUG && println("No intersections, old quadratic is smallest, Δa > 0, breaking.")
                break
            else

                # Compute the intersections
                term1 = -(Δb / 2 / Δa)
                term2 = sqrt(b2_minus_4ac) / 2 / Δa # Δa > 0
                root1, root2 = term1-term2, term1+term2

                DEBUG && println("Δa > 0   root1:", root1, "   root2:", root2)

                # Check where the intersections are and act accordingly
                if root1 >= right_endpoint
                    DEBUG && println("Two intersections to the right: old quadratic is smallest")
                    λ_prev, λ_curr = update_segment_do_nothing(λ_curr)
                elseif root2 <= left_endpoint
                    DEBUG && println("Two intersections to the left: old quadratic is smallest")
                    break # There will be no more intersections since Δa > 0
                elseif root1 <= left_endpoint && root2 >= right_endpoint
                    DEBUG && println("One intersections on either side, new quadratic is smallest")
                    λ_prev, λ_curr = update_segment_new(λ_prev, λ_curr, μ)
                elseif root1 > left_endpoint && root2 < right_endpoint
                    DEBUG && println("Two intersections: old-new-old")
                    λ_prev, λ_curr = update_segment_old_new_old(λ_curr, μ, root1, root2)
                    break # There will be no more intersections since Δa > 0
                elseif root1 > left_endpoint
                    DEBUG && println("Root 1 within the interval: old-new")
                    λ_prev, λ_curr = update_segment_old_new(λ_curr, μ, root1)
                elseif root2 < right_endpoint
                    DEBUG && println("Root 2 within the interval: new-old")
                    λ_prev, λ_curr = update_segment_new_old(λ_prev, λ_curr, μ, root2)
                    break # There will be no more intersections since Δa > 0
                else
                    error("Shouldn't end up here")
                end
            end

        elseif Δa < 0 # μ has lower curvature, i.e., μ is smallest on the sides
            if b2_minus_4ac <= ACCURACY
                # Zero (or one) roots
                λ_prev, λ_curr = update_segment_new(λ_prev, λ_curr, μ)
            else
                # Compute the intersections
                term1 = -(Δb / 2 / Δa)
                term2 = sqrt(b2_minus_4ac) / 2 / Δa # < 0
                root1, root2 = term1+term2, term1-term2
                DEBUG && println("Δa < 0   root1:", root1, "   root2:", root2)

                # Check where the intersections are and act accordingly
                if root1 >= right_endpoint || root2 <= left_endpoint
                    # No intersections, μ is smallest
                    λ_prev, λ_curr = update_segment_new(λ_prev, λ_curr, μ)
                elseif root1 <= left_endpoint && root2 >= right_endpoint
                    # No intersections, old quadratic is smallest, just step forward
                    λ_prev, λ_curr = update_segment_do_nothing(λ_curr)
                elseif root1 > left_endpoint && root2 < right_endpoint
                    # Two intersections within the interval
                    λ_prev, λ_curr = update_segment_new_old_new(λ_prev, λ_curr, μ, root1, root2)
                elseif root1 > left_endpoint
                    λ_prev, λ_curr = update_segment_new_old(λ_prev, λ_curr, μ, root1)
                elseif root2 < right_endpoint
                    λ_prev, λ_curr = update_segment_old_new(λ_curr, μ, root2)
                else
                    error("Shouldn't end up here")
                end
            end
        else # Δa == 0.0
            DEBUG && println("Δa == 0")
            DEBUG2 && println("Δa == 0 : $μ")
            if Δb == 0
                if Δc >= 0
                    λ_prev, λ_curr = update_segment_do_nothing(λ_curr)
                else
                    λ_prev, λ_curr = update_segment_new(λ_prev, λ_curr, μ)
                end
                continue
            end

            root = -Δc / Δb
            if Δb > 0
                if root < left_endpoint
                    λ_prev, λ_curr = update_segment_do_nothing(λ_curr)
                elseif root > right_endpoint
                    λ_prev, λ_curr = update_segment_new(λ_prev, λ_curr, μ)
                else
                    λ_prev, λ_curr = update_segment_new_old(λ_prev, λ_curr, μ, root)
                end
            else
                if root < left_endpoint
                    λ_prev, λ_curr = update_segment_new(λ_prev, λ_curr, μ)
                elseif root > right_endpoint
                    λ_prev, λ_curr = update_segment_do_nothing(λ_curr)
                else
                    DEBUG2 && println("Special case")
                    λ_prev, λ_curr = update_segment_old_new(λ_curr, μ, root)
                end
            end
        end

    end
    return
end

@inline function update_segment_do_nothing(λ_curr)
    return λ_curr, λ_curr.next
end

@inline function update_segment_new(λ_prev, λ_curr, μ)
    μ.has_been_used = true # TODO: Could be moved to the second clause
    if λ_prev.p === μ
        λ_prev.next = λ_curr.next
        v1, v2 = λ_prev, λ_curr.next
    else
        λ_curr.p = μ
        v1, v2 = λ_curr, λ_curr.next
    end
    return v1, v2
end

@inline function update_segment_old_new(λ_curr, μ, root)
    # ... λ_curr | λ_new ...
    μ.has_been_used = true
    λ_new =  PiecewiseQuadratic(μ, root, λ_curr.next)
    λ_curr.next = λ_new
    return λ_new, λ_new.next
end

@inline function update_segment_new_old(λ_prev, λ_curr, μ, root)
    # ... λ_prev | (new segment) | λ_curr ...
    if λ_prev.p === μ
        λ_curr.left_endpoint = root
    else
        μ.has_been_used = true
        λ_prev.next = PiecewiseQuadratic(μ, λ_curr.left_endpoint, λ_curr)
        λ_curr.left_endpoint = root
    end
    return λ_curr, λ_curr.next
end

@inline function update_segment_new_old_new(λ_prev, λ_curr, μ, root1, root2)
    # Insert new segments with μ before and after λ_curr
    # ... λ_prev : (new segment) | λ_curr | (new segment) ...
    # ( μ.has_been_used = true is set inside the called funcitons )
    update_segment_new_old(λ_prev, λ_curr, μ, root1)
    return update_segment_old_new(λ_curr, μ, root2)
end

@inline function update_segment_old_new_old(λ_curr, μ, root1, root2)
    # ... λ_curr | λ_1  | λ_2 ...
    # The second new piece contains a copy of the polynomial in λ_curr

    μ.has_been_used = true
    λ_2 =  PiecewiseQuadratic(λ_curr.p, root2, λ_curr.next)
    λ_1 =  PiecewiseQuadratic(μ, root1, λ_2)
    λ_curr.next = λ_1
    return λ_2, λ_2.next
end

"""
    μ = minimize_wrt_x2(qf::QuadraticForm{T}, p::QuadraticPolynomial{T}, μ=QuadraticPolynomial{T}()) where {T}
Takes a quadratic form in `[x₁; x₂]` and a polynomial in `x₂`
and returns the minimum of the sum wrt to `x₂`,
i.e. μ(x₁) = min_x₂{ qf(x₁,x₂) +  p(x₂) }`

The input `μ` can be pre-allocated on input and will then be changed.
"""
@inline function minimize_wrt_x2(qf::QuadraticForm{T}, p::QuadraticPolynomial{T}, μ=QuadraticPolynomial{T}()) where T
    P = qf.P
    q = qf.q
    r = qf.r

    P22_new = P[2,2] + p.a

    if P22_new > 0
        μ.a = P[1,1] - P[1,2]^2 / P22_new
        μ.b = q[1] - P[1,2]*(q[2] + p.b) / P22_new
        μ.c = (r + p.c) - (q[2] + p.b)^2 / P22_new / 4
    elseif P22_new == 0
        μ.a = P[1,1]
        μ.b = q[1]
        μ.c = r + p.c
    else
        error("Piecewise quadratic cost should be positive (semi-)definite. Please submit a bug report.")
    end
    return μ
end

global counter1
global counter2

"""
    Λ = pwq_dp_constrained(ℓ::AbstractTransitionCost{T}, V_N::QuadraticPolynomial{T}, M::Integer, upper_bound=Inf) where {T}

Given the transition costs `ℓ[i,j](y_i,y_j)` and the cost at the endpoint `V_N(y_N)` find all solutions `f` with up to `M` segments for the problem

V_i^m = minimize_f^M [ Σ_{k=1}^i { ℓ[k,k+1](f(k),f(k+1)) } + V_N(f(N)) ]
s.t          f(k) being continuous piecewise linear with `m` segements.

i.e. `Λ[m,i]` contains the best (in `ℓ` cost) continuous piecewise linear function `f` with up to `M` segments over the interval `i` to `N`
"""
function pwq_dp_constrained(ℓ::AbstractTransitionCost{T}, V_N::QuadraticPolynomial{T}, M::Integer, upper_bound=Inf) where T
    #global counter1
    #global counter2
    #counter1 = 0
    #counter2 = 0
    N = size(ℓ, 2)
    @assert M <= N-1 "Cannot have more segments than N-1."

    V_N = deepcopy(V_N)
    V_N.time_index = N

    Λ = Array{PiecewiseQuadratic{T}}(undef, M, N)

    for i=1:N-1
        p = minimize_wrt_x2(ℓ[i, N], V_N)
        p.time_index = i
        p.ancestor = V_N
        Λ[1, i] = create_new_pwq(p)
    end

    μ = QuadraticPolynomial{T}()
    upper_bound_inner = Inf
    global times
    for m=2:M
        #println("m: $m")
        for i=1:N-m
            Λ_new = create_new_pwq(T)
            if OPTIMIZE
                if min(upper_bound, upper_bound_inner) < Inf
                    add_quadratic!(Λ_new, QuadraticPolynomial{T}(0.0, 0.0, min(upper_bound,upper_bound_inner)))
                end
            end

            for ip=i+1:N-m+1
                DEBUG && println("(m:$m, i:$i, ip:$ip)")

                for λ in Λ[m-1, ip]

                    minimize_wrt_x2(ℓ[i,ip], λ.p, μ)

                    DEBUG && println("Obtained μ = $μ")

                    if OPTIMIZE
                        μmin = unsafe_minimum(μ)
                        if (μmin > upper_bound || μmin > upper_bound_inner)
                            DEBUG && println("Breaking due to that $μmin > max($upper_bound, $upper_bound_inner)")
                            continue
                        end
                    end

                    DEBUG && println("Inserting...")

                    add_quadratic!(Λ_new, μ)

                    if μ.has_been_used == true
                        μ.time_index = i
                        μ.ancestor = λ.p
                        μ = QuadraticPolynomial{T}()
                        μ.has_been_used = false
                    end
                end
            end
            #remove_over(Λ_new, min(upper_bound_inner,upper_bound))
            Λ[m, i] = Λ_new
            if OPTIMIZE
                if i == 1
                    upper_bound_inner = find_minimum_value(Λ[m,1])
                    upper_bound_inner += sqrt(eps())
                end
            end
        end
    end
    return Λ
end

"""
TODO Update docstring:

Finds the set I=(i_1, ..., i_M) that is the solution to the regularization problem

minimize ∑ ℓ_(i, i+1)(y, y+1)  +  V_N(y_M)  +  ζ⋅card(I)

where ℓ are positive-definite quadratic forms.
"""
function pwq_dp_regularized(ℓ::AbstractTransitionCost{T}, V_N::QuadraticPolynomial{T}, ζ::T) where T
    N = size(ℓ, 2)

    Λ = Vector{PiecewiseQuadratic{T}}(undef, N)

    V_N = deepcopy(V_N)
    V_N.time_index = N
    Λ[N] = create_new_pwq(V_N)

    μ = QuadraticPolynomial{T}()

    for i=N-1:-1:1
        Λ_new = create_new_pwq(T)
        for ip=i+1:N

            ζ_level_insertion = false
            for λ in Λ[ip]
                #counter1 += 1

                minimize_wrt_x2(ℓ[i,ip], λ.p, μ)
                μ.c += ζ # add cost for break point

                add_quadratic!(Λ_new, μ)

                if μ.has_been_used
                    μ.time_index = i
                    μ.ancestor = λ.p
                    μ = QuadraticPolynomial{T}()
                    μ.has_been_used = false

                    ζ_level_insertion = true
                else
                    if ζ_level_insertion == false
                        if !poly_minus_constant_is_greater(Λ_new, μ, ζ)
                            ζ_level_insertion = true
                        end
                    end
                end
            end

            if ζ_level_insertion == false
                break
            end
        end
        Λ[i] = Λ_new
    end

    return Λ
end

function recover_optimal_index_set(Λ::PiecewiseQuadratic{T}, first_index=1) where T

    p, y, f = find_minimum(Λ)

    I = Vector{Int}(undef, 0)

    while true
        push!(I, p.time_index)

        if !isdefined(p, :ancestor);
            break;
        end

        p = p.ancestor
    end

    return I, y, f
end

"""
    Y, f = find_optimal_y_values(ℓ, V_N, I)
Given transition costs `ℓ`, cost of right endpoint `V_N`, and breakpoint indicies `I`
the optimal y-values `Y` and the optimal cost `f` are computed.
"""
function find_optimal_y_values(ℓ, V_N::QuadraticPolynomial, I)
    m = length(I) - 1

    P = zeros(m+1, m+1)
    q = zeros(m+1)

    # Add cost for the right endpoint
    P[m+1, m+1] = V_N.a
    q[m+1] = V_N.b
    r = V_N.c

    # Form quadratic cost function Y'*P*Y + q'*Y + r
    # corresponding to the y-values in the vector Y
    for j=1:m
        P[j:j+1,j:j+1] .+= ℓ[I[j], I[j+1]].P
        q[j:j+1] .+= ℓ[I[j], I[j+1]].q
        r += ℓ[I[j], I[j+1]].r
    end

    # find the optimal Y-vector, and compute the correspinding error
    Y = -(P \ q) / 2
    f = Y' * P * Y + q' * Y + r
    return Y, f
end

# TODO: Include mζ in the cost?!
"""
    I, Y, f = recover_solution(Λ::PiecewiseQuadratic{T}, ℓ, V_N::QuadraticPolynomial, first_index=1)

"""
function recover_solution(Λ::PiecewiseQuadratic, ℓ, V_N::QuadraticPolynomial, ζ=0.0)
    I, _, f_expected = recover_optimal_index_set(Λ, 1)
    Y, f = find_optimal_y_values(ℓ, V_N::QuadraticPolynomial, I)

    f_regularized = f + ζ*(length(I)-1) # Include regularization cost

    if !isapprox(f_regularized, f_expected, atol=1e-10)
        @warn "Recovered cost ($f_regularized) is not what was expected from value function ($f_expected). Solution might be incorrect."
    end

    if f_regularized < 0
        @warn "Computed cost ($f_regularized) < 0, if ≈ 0, this is probably due to numerical errors and nothing to worry about."
    end

    return I, Y, f
end

"""
Given a piecewise quadratic V_Λ represented by Λ, a polynomial μ, and a real number ζ
this function evaluates if
μ(y) > V_Λ(y) + ζ  ∀ y
""" # FIXME: Group polynomial and constant in tuple?
function poly_minus_constant_is_greater(Λ::PiecewiseQuadratic{T}, μ::QuadraticPolynomial{T}, ζ::Real) where T

    if Λ.next.left_endpoint == Inf
        return false
    end

    for λ_curr in Λ

        left_endpoint = λ_curr.left_endpoint
        right_endpoint = get_right_endpoint(λ_curr)

        Δa = μ.a - λ_curr.p.a
        Δb = μ.b - λ_curr.p.b
        Δc = (μ.c - ζ) - λ_curr.p.c

        b2_minus_4ac =  Δb^2 - 4*Δa*Δc

        if Δa > 0 # μ has greater curvature, i.e., μ is smallest in the middle if intersect
            #println("Δa > 0")
            if b2_minus_4ac <= 0
                # Zero (or one) intersections, old quadratic is smallest, just step forward
                #println("No intersections, old quadratic is smallest, Δa > 0, breaking.")
                return true
            else

                # Compute the intersections
                term1 = -(Δb / 2 / Δa)
                term2 = sqrt(b2_minus_4ac) / 2 / Δa # Δa > 0
                root1, root2 = term1-term2, term1+term2

                # Check where the intersections are and act accordingly
                if root1 >= right_endpoint
                    continue
                elseif root2 <= left_endpoint
                    return true # There will be no more intersections since Δa > 0
                else
                    return false
                end
            end

        elseif Δa < 0 # μ has lower curvature, i.e., μ is smallest on the sides
            #println("Δa < 0")
            if b2_minus_4ac <= 0
                # Zero (or one) roots
                return false
            else
                # Compute the intersections
                term1 = -(Δb / 2 / Δa)
                term2 = sqrt(b2_minus_4ac) / 2 / Δa # < 0
                root1, root2 = term1+term2, term1-term2

                # Check where the intersections are and act accordingly
                if root1 <= left_endpoint && root2 >= right_endpoint
                    # No intersection on either side of the interval, and μ is smallest
                    continue
                else
                    return false
                end
            end
        else # Δa == 0.0
            DEBUG2 && println("Δa == 0 : $μ")
            if Δb == 0
                if Δc >= 0

                else
                    return false
                end
                continue
            end

            root = -Δc / Δb
            if Δb > 0
                if root < left_endpoint

                elseif root > right_endpoint
                    return false
                else
                    return false
                end
            else
                if root < left_endpoint
                    return false
                elseif root > right_endpoint

                else
                    DEBUG2 && println("Special case")
                    return false
                end
            end
        end

    end
    return true
end

# Continuous case
""" `get_transition_costs(g, t, lazy; tol=1e-3)`
    Return `(ℓ, cost_last)`: the transition cost and end-cost.
    Slightly unsafe since `tol` is ignored (not needed) on discrete problems
"""
function get_transition_costs(g, t, lazy; tol=1e-3)
    T = Float64
    ℓ = lazy ? TransitionCostContinuous{T}(g, t, tol) :
               compute_transition_costs(g, t, tol)
    # Continouous case, no cost at endpoint
    cost_last = zero(QuadraticPolynomial{Float64})
    return ℓ, cost_last
end

# Discrete case, tol is not used here, but in signature to enable dispatch
function get_transition_costs(g::AbstractArray, t, lazy; tol=1e-3)
    T = promote_type(eltype(g), Float64)
    ℓ = lazy ? TransitionCostDiscrete{T}(g, t) :
               compute_discrete_transition_costs(g, t)
    if t[1] != 1 || t[end] != length(g)
        @warn "In fit_pwl_constrained: The supplied grid t only covers the range ($(t[1]),$(t[end])) while the range of indices for g is (1,$(length(g))). No costs will be considered outside the range of t."
    end
    # Discrete case, so cost at endpoint is quadratic
    cost_last = QuadraticPolynomial(1.0, -2*g[t[end]], g[t[end]]^2)
    return ℓ, cost_last
end

## Public interface
"""
    I, Y, v = fit_pwl_regularized(g::AbstractArray, ζ; t=1:length(g), lazy=true)
    I, Y, v = fit_pwl_regularized(g, t, ζ, tol=1e-3; lazy=true)
Approximate `g[k]` (or `g(t)`) with a continuous piecewise linear function `f` according to
    v = min_f `||f-g||₂^2 + ζ⋅(length(I)-2)`
where the norm is `sum((f[k]-g[k])²)` (or integral over `(f(t)-g(t))²` from `t[1]` to `t[end]`).
Returns:
    `I`: Vector of length `M`
    `Y`: Vector of length `M`
    `v`: Float64
such that the optimal function has breakpoints in `I` (or `t[I]`) and satisfies
`f(I) .== Y` (or `f(t[I]) .== Y`).
Kwargs:
`t` is optional parameter in the discrete case, restricting the set of possible gridpoints,
i.e. so that `f[t[I]] .== Y`.
`lazy` = true, means that the internal transition costs `ℓ[i,j]` will be calculated when needed.
`tol` specifies the relative tolerance sent to `quadg` kused when calculating the integrals (continuous case).
"""
function  fit_pwl_regularized(g::AbstractArray, ζ; t=1:length(g), lazy=true)
    ℓ, cost_last = get_transition_costs(g, t, lazy)
    fit_pwl_regularized_internal(ℓ, cost_last, ζ)
end

# Continuous function
function  fit_pwl_regularized(g, t, ζ, tol=1e-3; lazy=true)
    ℓ, cost_last = get_transition_costs(g, t, lazy, tol=tol)
    fit_pwl_regularized_internal(ℓ, cost_last, ζ)
end

function  fit_pwl_regularized_internal(ℓ, cost_last, ζ)
    Λ_reg = pwq_dp_regularized(ℓ, cost_last, ζ)
    #Get solution that starts at first index
    I, Y, f = recover_solution(Λ_reg[1], ℓ, cost_last, ζ)
    return I, Y, f
end

"""
    I, Y, v = fit_pwl_constrained(g::AbstractArray, M; t=1:length(g), lazy=true)
    I, Y, v = fit_pwl_constrained(g, t, M, tol=1e-3; lazy=true)
Approximate `g[k]` (or `g(t)`) with a continuous piecewise linear function `f` according to
    v = min_f `||f-g||₂^2`
    s.t. length(I)-2 = M
where the norm is `sum((f[k]-g[k])²)` (or integral over `(f(t)-g(t))²` from `t[1]` to `t[end]`).
Returns:
    `I`: Vector of length `M`
    `Y`: Vector of length `M`
    `v`: Float64
such that the optimal function has breakpoints in `I` (or `t[I]`) and satisfies
`f(I) .== Y` (or `f(t[I]) .== Y`).
Kwargs:
`t` is optional parameter in the discrete case, restricting the set of possible gridpoints,
i.e. so that `f[t[I]] .== Y`.
`lazy` = true, means that the internal transition costs `ℓ[i,j]` will be calculated when needed.
`tol` specifies the relative tolerance sent to `quadg` kused when calculating the integrals (continuous case).
"""
function  fit_pwl_constrained(g::AbstractArray, M; t=1:length(g), lazy=false)
    ℓ, cost_last = get_transition_costs(g, t, lazy)
    fit_pwl_constrained_internal(ℓ, cost_last, M)
end

# Continuous function
function  fit_pwl_constrained(g, t, M, tol=1e-3; lazy=false)
    ℓ, cost_last = get_transition_costs(g, t, lazy, tol=tol)
    fit_pwl_constrained_internal(ℓ, cost_last, M)
end

function  fit_pwl_constrained_internal(ℓ::AbstractTransitionCost{T}, cost_last, M) where T
    Λ = pwq_dp_constrained(ℓ, cost_last, M);

    Ivec = Vector{Vector{Int}}(undef, M)
    Yvec = Vector{Vector{T}}(undef, M)
    fvec = Vector{T}(undef, M)

    for m=1:M
        Ivec[m], Yvec[m], fvec[m] = recover_solution(Λ[m, 1], ℓ, cost_last)
    end

    return Ivec, Yvec, fvec
end


text:
"""
The Git module allows you to use command-line Git in your Julia
packages. This is implemented with the `git` function, which returns a
`Cmd` object giving access to command-line Git.
"""
module Git

import Git_jll

export git

include("git_function.jl")

end # module

ext: jl


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule VOLK_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("VOLK")
JLLWrappers.@generate_main_file("VOLK", UUID("13737630-5ab7-50a7-b423-42a4a7ec78b3"))
end  # module VOLK_jll


text:
@testset "Testing Lots of things" begin
    localcache = joinpath(@__DIR__, "cache", "fuzz.jls")

    backend = AutoTM.Backend("CPU")
    cache = AutoTM.Profiler.CPUKernelCache(localcache)

    # Queue up a collection of models and optimizers.
    batchsize = 16
    fns = [
        () -> AutoTM.Zoo.resnet_training(AutoTM.Zoo.Resnet50(), batchsize),
        () -> AutoTM.Zoo.vgg_training(AutoTM.Zoo.Vgg19(), batchsize),
        () -> AutoTM.Zoo.inception_v4_training(batchsize),
        # TODO: Finish profiling to add this test.
        #() -> AutoTM.Zoo.densenet_training(batchsize)
    ]

    # Go through the static and synchronous optimizers.
    #
    # Just do one where we know we'll get move nodes / a mix of DRAM and PMM
    # Keep the total number down to facilitate quicker testing.
    optimizers = [
        AutoTM.Optimizer.Static(1 // 1),
        AutoTM.Optimizer.Synchronous(1 // 1),
    ]

    for (i, f) in enumerate(fns)
        for (j, opt) in enumerate(optimizers)
            printstyled("Function $i -- Optimizer $j\n"; color = :green)
            @test AutoTM.Verifier.verify(backend, f, opt, cache)
        end
    end
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
include("../IO/PolymeshReader.jl")
include("MeshProcessor.jl")

function readopenfoammesh(casedir)
    polymeshdir = joinpath(casedir, "constant/polyMesh")
    polymeshfiles = getpolymeshfiles(polymeshdir)
    polymesh = readpolymesh(polymeshfiles)
    return processopenfoammesh(polymesh)
end


text:
#!/usr/bin/env julia
#
# bisection script to find at which Julia commit a package started failing
#
# start by setting-up a bisection:  git bisect start BAD_COMMIT GOOD_COMMIT
# then run this script:             git run ./path_to_script.jl PackageName

using DataFrames
using PkgEval
using LibGit2

repo = GitRepo(pwd())
head = LibGit2.head_oid(repo)
hash = string(head)
shorthash = chomp(read(`git rev-parse --short $hash`, String))

try
    julia_version = PkgEval.perform_julia_build(hash; precompile=false)

    # gather results. run each test three times to catch flaky errors
    isempty(ARGS) && error("You should specify at least a single package")
    results = PkgEval.run([julia_version, julia_version, julia_version], ARGS;
                             update_registry=false, retries=0)
    results[results[!, :status] .== :kill, :status] .= :fail    # merge killed and failed

    # analyze
    has_issues = false
    for group in groupby(results, :name)
        pkg = first(group).name
        if all(eachrow(group).status .== :ok)
            @info "$pkg passed tests"
        elseif any(eachrow(group).status .== :skip)
            error("$pkg was skipped")
        elseif any(eachrow(group).status .== :fail)
            has_issues = true
            log = "$pkg-$shorthash.log"
            test = first(filter(row->row.status == :fail, group))
            @warn "$pkg failed tests: $(PkgEval.reasons[test.reason]) (see $log)"
            open(log, "w") do io
                println(io, test.log)
            end
        end
    end

    # report to Git
    exit(has_issues ? 1 : 0)
catch ex
    @error "Could not check commit" exception=(ex, catch_backtrace())
    exit(125)
end


text:
module MidiTS
    using DelimitedFiles, LinearAlgebra

    include("constr_series.jl")

    export get_time_series, get_pitch_seq
    export pitches_gaps, notas_hertz, get_onon_notes, serie_notas, get_new_voice, get_onoff_notes, silence_gaps, min_voces, max_tempo, rounding, filter_undef
    export indice, find_more_voices, get_subs

end


text:
module FESpacesTests

include("runtests_1.jl")
include("runtests_2.jl")

end # module


text:
# using CuArrays
using WAV
using MFCC
using JLD2
using Base.Iterators
using IterTools

# make dictionary to map from phones to class numbers
const PHONEMES, PHN2IDX = let
   phonemes = split("h# q eh dx iy r ey ix tcl sh ow z s hh aw m t er l w aa hv ae dcl y axr d kcl k ux ng gcl g ao epi ih p ay v n f jh ax en oy dh pcl ah bcl el zh uw pau b uh th ax-h em ch nx eng")
   phonemes = convert(Vector{String}, phonemes)
   phn2idx = Dict(phoneme => i for (i, phoneme) ∈ enumerate(phonemes))
   phn2idx["sil"] = phn2idx["h#"]
   phonemes, phn2idx
end

# make dictionary to perform class folding
const FOLDINGS = Dict(
   "ao" => "aa",
   "ax" => "ah",
   "ax-h" => "ah",
   "axr" => "er",
   "hv" => "hh",
   "ix" => "ih",
   "el" => "l",
   "em" => "m",
   "en" => "n",
   "nx" => "n",
   "eng" => "ng",
   "zh" => "sh",
   "pcl" => "sil",
   "tcl" => "sil",
   "kcl" => "sil",
   "bcl" => "sil",
   "dcl" => "sil",
   "gcl" => "sil",
   "h#" => "sil",
   "pau" => "sil",
   "epi" => "sil",
   "ux" => "uw"
)

const WINDOW = 0.025 # s
const STEP   = 0.010 # s

"""
   build_features(wavfile::AbstractString, phnfile::AbstractString)

Extracts Mel filterbanks and associated labels from `wavfile` and `phnfile`.
"""
function build_features(wavfile::AbstractString, phnfile::AbstractString, Δorder::Integer=2)
   samples, sample_rate = wavread(wavfile)
   (size(samples, 2) == 1) || throw("input .wav file must be mono channeled")
   samples = vec(samples)
   mfccs, _, _ = mfcc(samples, sample_rate, :rasta; wintime=WINDOW, steptime=STEP)

   lines = readlines(phnfile)
   labels = similar(lines)
   boundaries = similar(lines, Int)
   for (i, line) ∈ enumerate(lines)
      _, boundary, labels[i] = split(line)
      boundaries[i] = parse(Int, boundary)
   end

   sampl_length = WINDOW * sample_rate
   sampl_interval = STEP * sample_rate
   half_frame_length = WINDOW / 2

   # begin generating sequence labels by looping through the MFCC frames
   nsegments = length(labels)
   label_idx = 1
   label_sequence = similar(labels, size(mfccs,1))
   label, boundary = labels[label_idx], boundaries[label_idx]
   for i ∈ axes(mfccs, 1)
      win_end = sampl_length + (i-1)sampl_interval
      # move on to next label if current frame of samples is more than half
      # way into next labeled section and there are still more labels to iterate through
      if (label_idx < nsegments) && (win_end - boundary > half_frame_length)
         label_idx += 1
         label, boundary = labels[label_idx], boundaries[label_idx]
      end
      label_sequence[i] = label
   end
   # get rid of the frames that were labeld as "q"
   idxs_to_keep = filter(i -> label_sequence[i] != "q", eachindex(label_sequence))
   label_sequence = label_sequence[idxs_to_keep]
   mfccs = mfccs[idxs_to_keep,:]
   # compute filterbank derivates, fitted over 2 consecutive frames
   Δmfccs = (Δmfcc for Δmfcc ∈ take(iterated(x -> deltas(x, 2), mfccs), Δorder))
   featuremat = mapreduce(x -> x', vcat, (mfccs, Δmfccs...))
   # convert from double precision to single precision for Flux
   featuremat = convert(Matrix{Float32}, featuremat)
   features = [featuremat[:,j] for j ∈ axes(featuremat,2)]
   # generate class numbers, there are 61 total classes, but only 39 are used after folding
   target = [PHN2IDX[label] for label ∈ label_sequence]
   return (features, target)
end

"""
   build_dataset(dir_data::AbstractString, path_out::AbstractString)

Extracts data from files in `dir_data`, builds featuremat, concatenates produced outputs and saves results as `path_out`.
"""
function build_dataset(dir_data::AbstractString, path_out::AbstractString; Δorder::Integer=2)
   Xys = mapreduce(vcat, walkdir(dir_data)) do (root, _, files)
      print("$(root)\r")
      phnfiles = (file for file ∈ files if !startswith(file, "SA") && endswith(file, ".PHN"))
      wavfiles = (file for file ∈ files if !startswith(file, "SA") && endswith(file, ".wav"))

      map(wavfiles, phnfiles) do wavfile, phnfile
         wavfilepath = joinpath(root, wavfile)
         phnfilepath = joinpath(root, phnfile)
         build_features(wavfilepath, phnfilepath, Δorder)
      end
   end
   Xs = first.(Xys)
   ys = last.(Xys)

   ispath(path_out) || mkpath(dirname(path_out))
   JLD2.@save path_out Xs ys PHONEMES
   return Xs, ys
end

const DIR_DATA_TST  = "data/TIMIT/TIMIT_wav/TEST"
const PATH_OUT_TST  = "data/TIMIT/TIMIT_MFCC/data_test.jld2"

const DIR_DATA_TRN = "data/TIMIT/TIMIT_wav/TRAIN"
const PATH_OUT_TRN = "data/TIMIT/TIMIT_MFCC/data_train.jld2"

@time build_dataset(DIR_DATA_TST, PATH_OUT_TST)
@time build_dataset(DIR_DATA_TRN, PATH_OUT_TRN)


text:
function PrepRD()

println("RD")

end


text:
# by Lazaro Alonso
#using GLMakie, Random # HIDE
using GLMakie, Proj4, NPZ, FileIO
using Colors, ColorSchemes
GLMakie.activate!()
let
    # download the data from here:
    # https://github.com/lazarusA/BeautifulMakie/tree/main/_assets/data

    x = npzread("./data/aggregationOSMx.npy")
    y = npzread("./data/aggregationOSMy.npy")
    aggOSM = npzread("./data/aggregationOSM.npy")'
    img = load("./data/osm_kbc.png") # osm_bmy.png

    # getting back to lon and lat coordinates
    source = Projection("+proj=webmerc +datum=WGS84")
    dest = Projection("+proj=longlat +datum=WGS84")
    lonlat = transform(source, dest, [x y])

    function toCartesian(lon, lat; r = 1, cxyz = (0,0,0) )
        lat, lon = lat*π/180, lon*π/180
        cxyz[1] + r * cos(lat) * cos(lon), cxyz[2] + r * cos(lat) * sin(lon), cxyz[3] + r *sin(lat)
    end
    function lonlat3D(lon, lat, data; cxyz = (0,0,0))
        xyzw = zeros(size(data)..., 3)
        for (i,lon) in enumerate(lon), (j,lat) in enumerate(lat)
            x, y, z = toCartesian(lon, lat; cxyz = cxyz)
            xyzw[i,j,1] = x
            xyzw[i,j,2] = y
            xyzw[i,j,3] = z
        end
        xyzw[:,:,1], xyzw[:,:,2], xyzw[:,:,3]
    end

    x, y, z = lonlat3D(lonlat[:,1], lonlat[:,2], aggOSM)

    fig = Figure(figure_padding = -650, resolution = (1400,1400),
        backgroundcolor = "#141414")
    ax = Axis3(fig, aspect = :data, viewmode = :fit, perspectiveness = 0.5,
        azimuth = 6.49150301535897, elevation = 0.8516713267948967)
    mesh!(Sphere(Point3(0),1), color = :black, shading = false)
    surface!(ax, x, y, z, color = rotr90(img), shading = false,
        transparency = false)
    hidedecorations!(ax)
    hidespines!(ax)
    fig[1,1] = ax
    fig
    save(joinpath(@__DIR__, "output", "proj4_to_Sphere.png"), fig) # HIDE
end
using Pkg # HIDE
Pkg.status(["GLMakie","FileIO", "ImageMagick", "QuartzImageIO", "Proj4", "Colors", "ColorSchemes"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
#function sync{AA,I}(FA::Follower, B::ArrayIndexingWrapper{AA,I,true})
function sync(FA::Follower, B)
    A = parent(FA)
    checksame_inds(A, B)
    iter, Bind, Bfunc = syncable(B)
    SyncedIterator(each(iter), (stripwrapper(A), Bind), (synciterfunc(A, iter), Bfunc))
end

# value-iterator
syncable(A) = index(A), stripwrapper(A), (A, i) -> (@inbounds ret = A[i]; ret)
# syncable{AA,I}(A::ArrayIndexingWrapper{AA,I,false}) = A, parent(A), (A, i) -> (println("2"); @inbounds ret = A[i]; ret)
# index-iterator
syncable(A::ArrayIndexingWrapper{AA,I,true}) where {AA,I} = A, parent(A), (A, i) -> i

synciterfunc(A, B) = _synciterfunc(A, extent_dims(B))
# value-iterator
_synciterfunc(A, d::Tuple{Int}) = (A, i) -> A[i[d[1]]]
# index-iterator
_synciterfunc(A::ArrayIndexingWrapper{AA,I,true}, d::Tuple{Int}) where {AA,I} = (A, i) -> i[d[1]]


text:
module test_bukdu_anonymous

using Test # @test
using Bukdu # Conn

get("/") do conn::Conn
    42
end

result = Router.call(get, "/")
@test result.resp.status == 200
@test result.got == 42
@test result.route.C === Bukdu.System.AnonymousController

post("/") do conn::Conn
    42
end

result = Router.call(post, "/")
@test result.resp.status == 200
@test result.got == 42
@test result.route.C === Bukdu.System.AnonymousController

post("/:year", :year=>Int) do conn::Conn
    conn.params.year
end

result = Router.call(post, "/42")
@test result.resp.status == 200
@test result.got == 42
@test result.route.C === Bukdu.System.AnonymousController

Routing.reset!()

end # module test_bukdu_anonymous


text:
using LinearAlgebra
using SimpleSolvers
using SimpleSolvers: initialize!, objective, solver_step!
using Test

include("optimizers_problems.jl")

# struct OptimizerTest{T} <: Optimizer end

# test_optim = OptimizerTest{Float64}()

# @test_throws ErrorException config(test_optim)
# @test_throws ErrorException status(test_optim)
# @test_throws ErrorException objective(test_optim)
# @test_throws ErrorException initialize!(test_optim)
# @test_throws ErrorException solver_step!(test_optim)

# function test_optimizer(opt) end

for method in (Newton(), BFGS(), DFP())
    for linesearch in (Static(0.8), Backtracking(), Quadratic(), Bisection())
        n = 1
        x = ones(n)
        opt = Optimizer(x, F; algorithm = method, linesearch = linesearch)

        @test config(opt) == opt.config
        @test status(opt) == opt.result.status

        solve!(x, opt)
        # println(opt)
        @test norm(minimizer(opt)) ≈ 0 atol=1E-7
        @test norm(minimum(opt)) ≈ F(0) atol=1E-7

        x = ones(n)
        opt = Optimizer(x, F; ∇F! = ∇F!, algorithm = method, linesearch = linesearch)
        solve!(x, opt)
        # println(opt)
        @test norm(minimizer(opt)) ≈ 0 atol=1E-7
        @test norm(minimum(opt)) ≈ F(0) atol=1E-7
    end
end


text:
using BenchmarkTools
using Ipopt
const TO = TrajectoryOptimization

T = Float64

# options
max_con_viol = 1.0e-8
verbose=false

opts_ilqr = iLQRSolverOptions{T}(verbose=verbose,live_plotting=:off)

opts_al = AugmentedLagrangianSolverOptions{T}(verbose=verbose,
    opts_uncon=opts_ilqr,
    cost_tolerance=1.0e-4,
    cost_tolerance_intermediate=1.0e-2,
    constraint_tolerance=max_con_viol,
    penalty_scaling=1000.,
    penalty_initial=1.)

opts_altro = ALTROSolverOptions{T}(verbose=verbose,
    opts_al=opts_al,
    projected_newton=false)

opts_ipopt = DIRCOLSolverOptions{T}(verbose=verbose,
    nlp=Ipopt.Optimizer(),
    feasibility_tolerance=max_con_viol)

# AL-iLQR
prob = copy(Problems.doubleintegrator_static)
U0 = deepcopy(controls(prob))
alilqr = AugmentedLagrangianSolver(prob, opts_al)
solve!(alilqr)
max_violation(alilqr)

@btime begin
    initial_controls!($alilqr, $U0)
    solve!($alilqr)
end

# ALTRO
altro = ALTROSolver(prob, opts_altro)
solve!(altro)
max_violation(altro)

@btime begin
    initial_controls!($altro, $U0)
    solve!($altro)
end

# Ipopt
prob_ipopt = copy(Problems.doubleintegrator_static)
prob_ipopt = TO.change_integration(prob_ipopt, HermiteSimpson)
rollout!(prob_ipopt)
ipopt = DIRCOLSolver(prob_ipopt)
ipopt.opts.verbose = false
solve!(ipopt)
max_violation(ipopt)

@btime solve!($ipopt)


text:
# define the type of data

mutable struct fixedData
    # static network data
    baseMVA :: Float64
    bType :: Dict{Int64,Any}

    IDList :: Array{Int64,1}
    genIDList :: Array{Int64,1}
    brList :: Array{Any,1}
    brRev :: Dict{Any,Any}

    Loc :: Dict{Int64,Any}
    LocRev :: Dict{Int64,Any}
    Vmax :: Dict{Int64,Any}
    Vmin :: Dict{Int64,Any}
    Pmax :: Dict{Int64,Any}
    Pmin :: Dict{Int64,Any}
    Qmax :: Dict{Int64,Any}
    Qmin :: Dict{Int64,Any}
    gs :: Dict{Int64,Any}
    bs :: Dict{Int64,Any}
    Vmag :: Dict{Int64,Any}
    Vang :: Dict{Int64,Any}
    Pd :: Dict{Int64,Any}
    Qd :: Dict{Int64,Any}
    Pg :: Dict{Int64,Any}
    Qg :: Dict{Int64,Any}
    RU :: Dict{Int64,Any}
    RD :: Dict{Int64,Any}

    g :: Dict{Tuple{Int64,Int64,Int64},Any}
    b :: Dict{Tuple{Int64,Int64,Int64},Any}
    bc :: Dict{Tuple{Int64,Int64,Int64},Any}
    θmax :: Dict{Tuple{Int64,Int64,Int64},Any}
    θmin :: Dict{Tuple{Int64,Int64,Int64},Any}
    rateA :: Dict{Tuple{Int64,Int64,Int64},Any}
    τ1 :: Dict{Tuple{Int64,Int64,Int64},Any}
    τ2 :: Dict{Tuple{Int64,Int64,Int64},Any}
    σ :: Dict{Tuple{Int64,Int64,Int64},Any}

    cp :: Dict{Any,Any}
    cq :: Dict{Any,Any}
    cz :: Any

    busInd :: Dict{Any,Any}
    branchDict1 :: Dict{Any,Any}
    branchDict2 :: Dict{Any,Any}
    connectPair :: Array{Any,1}
    connectDict :: Dict{Any,Any}
    kpDict :: Dict{Any,Any}

end

struct costDataType
    model :: Int64
    upCost :: Float64
    downCost :: Float64
    n :: Int64
    params :: Array{Float64,1}
end

struct probDistrn
    # probability distribution of disruption timing
    tDistrn :: Dict{Int64,Any}

    # probability distribution of disruption magnitude
    # ωDistrn: keys are the indices, values are probability
    # ωDict: keys are the indices, values are lists of broken components
    ωDistrn :: Dict{Any,Any}
    ωDict :: Dict{Any,Any}
    ωτ :: Dict{Any,Any}
end

mutable struct batteryData
    # battery information: charging/discharging factor, capacity
    IDList :: Array{Any,1}
    Loc :: Dict{Int64,Any}
    bInv :: Dict{Int64,Any}
    ηα :: Dict{Int64,Any}
    ηβ :: Dict{Int64,Any}
    cap :: Dict{Int64,Any}
    cost :: Dict{Int64,Any}
    uCap :: Dict{Int64,Any}
end

struct solData
    # solution information: generaion level, storage level, maximum generation output
    sp :: Dict{Any,Any}
    sq :: Dict{Any,Any}
    w :: Dict{Any,Any}
    u :: Dict{Any,Any}
    lp :: Dict{Any,Any}
    lq :: Dict{Any,Any}
    zp :: Dict{Any,Any}
end

struct demandData
    # demand data: active demand and reactive demand
    T :: Int64
    pd :: Dict{Any,Any}
    qd :: Dict{Any,Any}
end

struct cutData
    # cut data: slope for sp, w and u and the intercept
    solStatus :: Any
    λ :: Dict{Any,Any}
    γ :: Dict{Any,Any}
    μ :: Dict{Any,Any}
    vhat :: Float64
end


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

# ╔═╡ 2ea1ccbc-8e0d-11eb-29c3-05a91ea18a6f
using TreeKnit
using TreeKnit.SplitGraph
using Test
using TreeTools

println("##### Basic resolving #####")

t1 = node2tree(TreeTools.parse_newick("(((A,B),C),D)"))
t2 = node2tree(TreeTools.parse_newick("(B,C,(A,D))"))
@testset "1" begin
	newsplits = resolve!(t1,t2)
	@test length(newsplits[1]) == 0
	@test length(newsplits[2]) == 0
end

t1 = read_tree("$(dirname(pathof(TreeKnit)))/..//test/resolving/tree1.nwk")
t2 = read_tree("$(dirname(pathof(TreeKnit)))/..//test/resolving/tree2.nwk")
@testset "2" begin
	newsplits = resolve!(t1,t2)
	@test newsplits[1].splits == [Split([1,2])]
	@test newsplits[2] == [
		["A1", "A2", "B", "C"],
		["D1", "D2"],
		["E", "F", "G", "H"],
	]
end

t1 = node2tree(TreeTools.parse_newick("(A,B,C,D)"))
t2 = node2tree(TreeTools.parse_newick("(A,(B,C,D))"))
t3 = node2tree(TreeTools.parse_newick("(A,B,(C,D))"))

@testset "Node naming" begin
	ns12 = resolve!(t1,t2)
	ns13 = resolve!(t1,t3)
	@test isempty(ns12[2])
	@test isempty(ns12[2])
	@test haskey(t1.lnodes, "RESOLVED_1")
	@test haskey(t1.lnodes, "RESOLVED_2")
	@test !haskey(t1.lnodes, "RESOLVED_3")
end

println("##### Resolve using pre-inferred MCCs #####")
t3 = read_tree("$(dirname(pathof(TreeKnit)))/..//test/resolving/tree3.nwk")
t4 = read_tree("$(dirname(pathof(TreeKnit)))/..//test/resolving/tree4.nwk")
mccs = read_mccs("$(dirname(pathof(TreeKnit)))/..//test/resolving/mccs34.dat")
@testset "Resolve with pre-inferred MCCs" begin
	rS = resolve!(t3, t4, mccs)
	@test rS[1] == [["F", "G"]]
	@test rS[2] == [["E", "F", "G"]]
end

# println("##### Resolve using MCC inference #####")

# # ╔═╡ 93ccb338-8e0d-11eb-0f36-49bebaf5570b
# function infer_mccs(trees;kwargs...)
# 	TreeKnit.runopt(
# 		OptArgs(;seq_lengths = Dict(s=>1 for s in keys(trees)), kwargs...),
# 		trees
# 	)
# end

# # ╔═╡ 45e695a6-8e0d-11eb-178d-97bf0cb45379
# nwk1 = "(A:1.,B:1.,(C:1.,D:1.):1.)"

# # ╔═╡ 6165e6e0-8e0d-11eb-1f87-b94b13b3ec95
# nwk2 = "(((A:1.,B:1.):1.,C:1.):1.,D:100.)"

# # ╔═╡ 6f0f4310-8e0d-11eb-1110-9d6d12d0cbae
# t1 = node2tree(TreeTools.parse_newick(nwk1))

# # ╔═╡ 7bdceb1a-8e0d-11eb-0160-0d60f81fd810
# t2= node2tree(TreeTools.parse_newick(nwk2))

# # ╔═╡ 90c5c30c-8e0e-11eb-00a0-45ccc42a2500
# trees = Dict(1=>deepcopy(t1), 2=>deepcopy(t2))

# # ╔═╡ fa1be92c-8e0d-11eb-0f5a-21209d7ffc68
# @testset "Inner function _resolve_from_mccs!" begin
# 	ct = deepcopy(trees)
# 	nS = TreeKnit._resolve_from_mccs!(infer_mccs, ct)
# 	@test isempty(nS[2])
# 	@test nS[1] == [["A", "B"]]
# end

# # ╔═╡ ef1e0a34-8e10-11eb-198d-1d31d25d8167
# nwk3 = "(((A1:1,A2:1):1.,A3:1,A4:1):1.,B1:1.,B2:1,(C:1.,D:1.):1)"

# # ╔═╡ 117d3fe6-8e11-11eb-02bf-6b21357867d7
# t3 = node2tree(TreeTools.parse_newick(nwk3))

# # ╔═╡ 1be36e60-8e11-11eb-2ca0-e59c20fa35a4
# nwk4 = "(((((A1:1,A2:1,A3:1):1,A4:1):1.,(B1:1.,B2:1):1):1,C:1.):1.,D:100.)"

# # ╔═╡ 50b097f8-8e11-11eb-1a43-ed66b32e8ee4
# t4 = node2tree(TreeTools.parse_newick(nwk4))

# # ╔═╡ 5fb91e64-8e11-11eb-3ac7-4906f631bd80
# trees2 = Dict(1=>deepcopy(t3), 2=>deepcopy(t4))

# # ╔═╡ 8d4d48b4-8e11-11eb-2445-f5196c3666b3
# @testset "Outer function resolve_from_mccs!" begin
# 	ct2 = deepcopy(trees2)
# 	nS2 = TreeKnit.resolve_from_mccs!(infer_mccs, ct2)
# 	@test nS2[2] == [["A1", "A2"]]
# 	@test nS2[1] == [["A1", "A2", "A3"], ["B1", "B2"], ["A1", "A2", "A3", "A4", "B1", "B2"]]
# end

# # # Testing max-clique splits
# # # trees 2 and 3 try to resolve 1 in a different way.
# # nwk1 = "((A,B,C),(D,E,F,G))"
# # nwk2 = "(((A,B),C),(D,(E,(F,G))))"
# # nwk3 = "(((A,B),C),(D,(E,F),G))"
# # trees = Dict(i=>node2tree(parse_newick(nwk)) for (i,nwk) in enumerate([nwk1, nwk2, nwk3]))
# # MCCs = TreeKnit._computeMCCs(infer_mccs, trees)
# # resolvable_splits = TreeKnit.new_splits(trees, MCCs)
# # @testset "max-clique splits" begin
# # 	S1 = TreeKnit.max_clique_splits(resolvable_splits[1])
# # 	S2 = TreeKnit.max_clique_splits(resolvable_splits[2])
# # 	S3 = TreeKnit.max_clique_splits(resolvable_splits[3])
# # 	@test length(S1) == 3 # would be 0 with compat
# # 	@test length(S2) == 0 # is already resolved
# # 	@test length(S3) == 1 # one split from tree 2 (depending if E or G is an MCC)
# # end

# ╔═╡ Cell order:
# ╠═2ea1ccbc-8e0d-11eb-29c3-05a91ea18a6f
# ╠═93ccb338-8e0d-11eb-0f36-49bebaf5570b
# ╠═45e695a6-8e0d-11eb-178d-97bf0cb45379
# ╠═6165e6e0-8e0d-11eb-1f87-b94b13b3ec95
# ╠═6f0f4310-8e0d-11eb-1110-9d6d12d0cbae
# ╠═7bdceb1a-8e0d-11eb-0160-0d60f81fd810
# ╠═90c5c30c-8e0e-11eb-00a0-45ccc42a2500
# ╠═fa1be92c-8e0d-11eb-0f5a-21209d7ffc68
# ╠═ef1e0a34-8e10-11eb-198d-1d31d25d8167
# ╠═117d3fe6-8e11-11eb-02bf-6b21357867d7
# ╠═1be36e60-8e11-11eb-2ca0-e59c20fa35a4
# ╠═50b097f8-8e11-11eb-1a43-ed66b32e8ee4
# ╠═5fb91e64-8e11-11eb-3ac7-4906f631bd80
# ╠═8d4d48b4-8e11-11eb-2445-f5196c3666b3

max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/TreeKnit

max_forks_repo_name: PierreBarrat/TreeKnit


text:
load("debug.jl")

module TestInstrument
using Base, Debug
import Debug.trap

trap(line, file) = println("trap: line = $line, file = $(repr(file))")

code = quote
    function()
        i = 1
        while i < 4
            println(i)
            i = i+1
        end
#         for i=1:3
#             println(i)
#         end
    end
end

icode = instrument(code)
println(icode)

eval(icode)()

acode = analyze(code)

end


text:
function bullsandcows ()
    bulls = cows = turns = 0
    result = (s = [] ; while length(unique(s))<4 push!(s,rand('1':'9')) end; unique(s))
    println("A game of bulls and cows!")
    while (bulls != 4)
      print("Your guess? ")
      guess = collect(chomp(readline(STDIN)))
      if ! (length(unique(guess)) == length(guess) == 4 && all(isdigit,guess))
         println("please, enter four distincts digits")
         continue
      end
      bulls = sum(map(==, guess, result))
      cows = length(intersect(guess,result)) - bulls
      println("$bulls bulls and $cows cows!") ; turns += 1
    end
    println("You win! You succeeded in $turns guesses.")
  end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#=
# MNIST

We begin by importing the required packages.
We load MNIST via the MLDatasets.jl package.
=#

import Makie
import CairoMakie
using Random: bitrand
using Statistics: mean
import MLDatasets
import Flux
import RestrictedBoltzmannMachines as RBMs
using ValueHistories: MVHistory
import CenteredRBMs
using RestrictedBoltzmannMachines: BinaryRBM
nothing #hide

#=
Useful function to plot MNIST digits.
=#

"""
    imggrid(A)

Given a four dimensional tensor `A` of size `(width, height, ncols, nrows)`
containing `width x height` images in a grid of `nrows x ncols`, this returns
a matrix of size `(width * ncols, height * nrows)`, that can be plotted in a heatmap
to display all images.
"""
function imggrid(A::AbstractArray{<:Any,4})
    reshape(permutedims(A, (1,3,2,4)), size(A,1)*size(A,3), size(A,2)*size(A,4))
end

#=
Load the MNIST dataset. We only work with 0s and 1s for speed.
=#

Float = Float32
train_x, train_y = MLDatasets.MNIST.traindata()
train_x = Array{Float}(train_x[:, :, train_y .∈ Ref((0,1))] .> 0.5)
println(size(train_x, 3), " train samples, with ", count(train_y .== 0), " zeros and ", count(train_y .== 1), " ones.")
nothing #hide

#=
Initialize and train a centered RBM
=#

rbm_c = CenteredRBMs.center(BinaryRBM(Float, (28,28), 400))
RBMs.initialize!(rbm_c, train_x) # centers from data
batchsize = 256
optim = Flux.ADAM()
vm = bitrand(28, 28, batchsize) # fantasy chains
history_c = MVHistory()
push!(history_c, :lpl, mean(RBMs.log_pseudolikelihood(rbm_c, train_x)))
@time for epoch in 1:100 # track pseudolikelihood every 5 epochs
    RBMs.pcd!(rbm_c, train_x; epochs=5, vm, history=history_c, batchsize, optim)
    push!(history_c, :lpl, mean(RBMs.log_pseudolikelihood(rbm_c, train_x)))
end
rbm_c = CenteredRBMs.uncenter(rbm_c) # convert to equivalent RBM (without offsets)
nothing #hide

#=
For comparison, we also train a normal (uncentered) RBM.
=#

rbm_u = BinaryRBM(Float, (28,28), 400)
RBMs.initialize!(rbm_u, train_x)
vm = bitrand(28, 28, batchsize)
history_u = MVHistory()
push!(history_u, :lpl, mean(RBMs.log_pseudolikelihood(rbm_u, train_x)))
@time for epoch in 1:100 # track pseudolikelihood every 5 epochs
    RBMs.pcd!(rbm_u, train_x; epochs=5, vm, history=history_u, batchsize, optim)
    push!(history_u, :lpl, mean(RBMs.log_pseudolikelihood(rbm_u, train_x)))
end
nothing #hide

# Plot log-pseudolikelihood of train data during learning.

fig = Makie.Figure(resolution=(600, 300))
ax = Makie.Axis(fig[1,1], xlabel="epochs", ylabel="pseudolikelihood")
Makie.lines!(ax, get(history_u, :lpl)..., label="normal")
Makie.lines!(ax, get(history_c, :lpl)..., label="centered")
Makie.axislegend(ax, position=:rb)
fig

# Seconds per epoch.

fig = Makie.Figure(resolution=(600, 300))
ax = Makie.Axis(fig[1,1], xlabel="epoch", ylabel="seconds")
Makie.lines!(ax, get(history_u, :Δt)..., label="normal")
Makie.lines!(ax, get(history_c, :Δt)..., label="centered")
Makie.axislegend(ax, position=:rt)
fig

# Log-pseudolikelihood vs. computation time instead of epoch count.

fig = Makie.Figure(resolution=(600, 300))
ax = Makie.Axis(fig[1,1], xlabel="seconds", ylabel="pseudolikelihood")
Makie.lines!(ax, cumsum([0; get(history_u, :Δt)[2]])[1:5:end], get(history_u, :lpl)[2], label="normal")
Makie.lines!(ax, cumsum([0; get(history_c, :Δt)[2]])[1:5:end], get(history_c, :lpl)[2], label="centered")
Makie.axislegend(ax, position=:rb)
fig

# Now we do the Gibbs sampling to generate RBM digits.

nrows, ncols = 10, 15
@time fantasy_x_c = RBMs.sample_v_from_v(rbm_c, bitrand(28,28,nrows*ncols); steps=10000)
@time fantasy_x_u = RBMs.sample_v_from_v(rbm_u, bitrand(28,28,nrows*ncols); steps=10000)
nothing #hide

# Plot the resulting samples.

# Normal RBM.

fig = Makie.Figure(resolution=(40ncols, 40nrows))
ax = Makie.Axis(fig[1,1], yreversed=true)
Makie.image!(ax, imggrid(reshape(fantasy_x_u, 28, 28, ncols, nrows)), colorrange=(1,0))
Makie.hidedecorations!(ax)
Makie.hidespines!(ax)
fig

# Centered RBM.

fig = Makie.Figure(resolution=(40ncols, 40nrows))
ax = Makie.Axis(fig[1,1], yreversed=true)
Makie.image!(ax, imggrid(reshape(fantasy_x_c, 28, 28, ncols, nrows)), colorrange=(1,0))
Makie.hidedecorations!(ax)
Makie.hidespines!(ax)
fig


text:
# a simple implementation avoiding vectors should make it simpler to port to Python and Julia
function CrossProduct(v, i, j, result)
    @inbounds result[1] = v[i][2]*v[j][3]-v[i][3]*v[j][2]
	@inbounds result[2] = v[i][3]*v[j][1]-v[i][1]*v[j][3]
    @inbounds result[3] = v[i][1]*v[j][2]-v[i][2]*v[j][1]
	@inbounds result[4] = 0
    return 0
end

function Thrust(particleArray)
    entries = length(particleArray)
    PartVector = MVector{4, Float64}(0., 0., 0., 0.)
    RefVector = MVector{4, Float64}(0., 0., 0., 0.)
    SumVector = MVector{4, Float64}(0., 0., 0., 0.)
    maxVal = 0.0

    if (entries<2)
        println(stderr, "Thrust : Not enough entries")
        return -1.0
    end
    # start calcualting the thrust
    @inbounds for i in 1:entries
        SumVector += particleArray[i]
    end

    # now do the thrust
    @inbounds for i = 1:entries-1
        @inbounds for j = i+1:entries
            CrossProduct(particleArray, i, j, RefVector)
			RefVector ./= sqrt(sum(RefVector.^2))
			PartVector .*= 0
            # Add all momenta with sign; two choices for each reference particle.
            @inbounds for k = 1:entries
                if k != i && k != j
                    dotProduct = particleArray[k] ⋅ RefVector
					#print(sign(dotProduct), "\t", particleArray[k], "\t", PartVector)
					#println(sign(dotProduct) * particleArray[k])
                    PartVector .+= sign(dotProduct) * particleArray[k]
                end
            end
            c1 = 0
            c2 = 0
            for l = 1:4
                if l == 0
                    c1 = 1
                    c2 = 1
                elseif l == 1
                    c1 = 1
                    c2 = -1
                elseif l == 2
                    c1 = -1
                    c2 = 1
                else
                    c1 = -1
                    c2 = -1
                end

				tempVal = 0.0
                @inbounds for n=1:3 # only copy three of the four as we overwrite the fourth element again ...
                    tempVal += (PartVector[n] + c1*particleArray[i][n] + c2*particleArray[j][n])^2
                end
                tempVal = sqrt(tempVal)
                if (tempVal > maxVal)
                    maxVal = tempVal
                end
            end
        end
    end
    return maxVal / SumVector[4]
end


text:
"""
    FeMat{T,S}

Term with an explicit, constant matrix representation

# Fields
* `x`: matrix
* `wtx`: weighted matrix
* `piv`: pivot `Vector{Int}`` for pivoted Cholesky factorization of `wtx'wtx`
* `rank`: computational rank of `x`
* `cnames`: vector of column names
"""
mutable struct FeMat{T,S<:AbstractMatrix}
    x::S
    wtx::S
    piv::Vector{Int}
    rank::Int
    cnames::Vector{String}
end

function FeMat(X::AbstractMatrix, cnms)
    T = eltype(X)
    ch = statscholesky(Symmetric(X'X))
    pivot = ch.piv
    Xp = all(pivot .== 1:size(X, 2)) ? X : X[:, ch.piv]
    FeMat{T,typeof(X)}(Xp, Xp, pivot, ch.rank, cnms[pivot])
end

function reweight!(A::FeMat{T}, sqrtwts::Vector{T}) where {T}
    if !isempty(sqrtwts)
        if A.x === A.wtx
            A.wtx = similar(A.x)
        end
        mul!(A.wtx, Diagonal(sqrtwts), A.x)
    end
    A
end

Base.adjoint(A::FeMat) = Adjoint(A)

Base.eltype(A::FeMat{T}) where {T} = T

fullrankwtx(A::FeMat) = rank(A) == size(A, 2) ? A.wtx : A.wtx[:, 1:rank(A)]

Base.length(A::FeMat) = length(A.wtx)

LinearAlgebra.rank(A::FeMat) = A.rank

Base.size(A::FeMat) = size(A.wtx)

Base.size(A::Adjoint{T,<:FeMat{T}}) where {T} = reverse(size(A.parent))

Base.size(A::FeMat, i) = size(A.wtx, i)

Base.copyto!(A::FeMat{T}, src::AbstractVecOrMat{T}) where {T} = copyto!(A.x, src)

*(adjA::Adjoint{T,<:FeMat{T}}, B::FeMat{T}) where {T} =
    fullrankwtx(adjA.parent)'fullrankwtx(B)

LinearAlgebra.mul!(R::StridedVecOrMat{T}, A::FeMat{T}, B::StridedVecOrMat{T}) where {T} =
    mul!(R, A.x, B)

LinearAlgebra.mul!(C, adjA::Adjoint{T,<:FeMat{T}}, B::FeMat{T}) where {T} =
    mul!(C, fullrankwtx(adjA.parent)', fullrankwtx(B))


text:
module MaxCut

using LinearAlgebra
using Convex
using SCS

"The Goemans-Williamson algorithm for the MAXCUT problem."

function maxcut(W::Matrix{<:Real}; iter::Int=100, tol::Real=0)
	"Partition a graph into two disjoint sets such that the sum of the edge weights
	which cross the partition is as large as possible (known to be NP-hard)."

	"A cut of a graph can be produced by assigning either 1 or -1 to each vertex. The Goemans-Williamson
	algorithm relaxes this binary condition to allow for vector assignments drawn from the (n-1)-sphere
	(choosing an n-1 dimensional space will ensure seperability). This relaxation can then be written as
	a semidefinite program. Once the optimal vector assignments are found, origin centered hyperplanes are generated and
	their corresponding cuts evaluated. After 'iter' trials, or when the desired tolerance is reached,
	the hyperplane with the highest corresponding binary cut is used to partition the vertices."

	"W:		Adjacency matrix."
	"tol:	Maximum acceptable distance between a cut and the MAXCUT upper bound."
	"iter:	Maximum number of hyperplane iterations before a cut is chosen."

	LinearAlgebra.checksquare(W)
	issymmetric(W)					|| throw(ArgumentError("Adjacency matrix must be symmetric."))
	all(W .>= 0)					|| throw(ArgumentError("Adjacency matrix must be nonnegative."))
	all(iszero.(diag(W)))			|| throw(ArgumentError("Diagonal of adjacency matrix must be zero (no self loops)."))
	(tol >= 0)						|| throw(ArgumentError("The tolerance must be nonnegative."))
	(iter > 0)						|| throw(ArgumentError("The number of iterations must be a positive integer."))

	"This is the standard SDP Relaxation of the MAXCUT problem, a reference can be found at,
	http://www.sfu.ca/~mdevos/notes/semidef/GW.pdf"
	k = size(W, 1)
	S = Semidefinite(k)

	expr = dot(W, S)
	constr = [S[i,i] == 1.0 for i in 1:k]
	problem = minimize(expr, constr...)
	solve!(problem, SCS.Optimizer(verbose=false))

	## ensure symmetric positive-definite
	A = 0.5 * (S.value + S.value')
	A += (max(0, -eigmin(A)) + 1e-14) * Matrix(I, size(A, 1), size(A, 1))

	X = Matrix(cholesky(A))

	## a non-trivial upper bound on MAXCUT
	upperbound = (sum(W) - dot(W, S.value)) / 4

	## random origin-centered hyperplanes, generated to produce partitions of the graph
	max_cut = -1
	max_partition = nothing

	for i in 1:iter
		gweval = X' * randn(k)
		partition = (findall(gweval .>= 0), findall(gweval .< 0))
		cut = sum(W[partition...])

		if cut > max_cut
			max_partition = partition
			max_cut = cut
		end

		(upperbound - max_cut < tol) && break
		(i == iter) && println("Max iterations reached.")
	end
	return max_cut, max_partition
end

export maxcut

end


text:
"""
    StaticUnitRange

Supertype for `UnitSRange` and `UnitMRange`. It's subtypes should behave
identically to `UnitRange`.
"""
abstract type StaticUnitRange{T<:Real} <: AbstractUnitRange{T} end

Base.firstindex(::StaticUnitRange) = 1

Base.lastindex(r::StaticUnitRange) = length(r)

"""
    UnitSRange

A static range parameterized by a `start` and `stop` of type `T`, filled with
elements spaced by `1` from `start` until `stop` is exceeded. The syntax `a:b`
with `a` and `b` both `Integer`s creates a `UnitRange`.
"""
struct UnitSRange{T,F,L} <: StaticUnitRange{T}

    function UnitSRange{T}(start, stop) where {T<:Real}
        return new{T,start,Base.unitrange_last(start,stop)}()
    end
end

UnitSRange{T}(r::AbstractUnitRange) where {T<:Real} = UnitSRange{T}(first(r), last(r))

function Base.getproperty(r::StaticUnitRange, s::Symbol)
    if s === :start
        return first(r)
    elseif s === :stop
        return last(r)
    else
        error("type $(typeof(r)) has no property $s")
    end
end

"""
    UnitMRange

A mutable range parameterized by a `start` and `stop` of type `T`, filled with
elements spaced by `1` from `start` until `stop` is exceeded. The syntax `a:b`
with `a` and `b` both `Integer`s creates a `UnitRange`.
"""
mutable struct UnitMRange{T<:Real} <: StaticUnitRange{T}
    start::T
    stop::T

    function UnitMRange{T}(start, stop) where {T<:Real}
        return new(start, Base.unitrange_last(start,stop))
    end
end

UnitMRange{T}(r::AbstractUnitRange) where {T<:Real} = UnitMRange{T}(first(r), last(r))

function Base.setproperty!(r::UnitMRange, s::Symbol, val)
    if s === :start
        return set_first!(r, val)
    elseif s === :stop
        return set_last!(r, val)
    else
        error("type $(typeof(r)) has no property $s")
    end
end

for (F,f) in ((:M,:m), (:S,:s))
    UR = Symbol(:Unit, F, :Range)
    frange = Symbol(f, :range)
    @eval begin
        Base.AbstractUnitRange{T}(r::R) where {T,R<:$(UR)} = $(UR){T}(r)
        $(UR)(start::T, stop::T) where {T<:Real} = $(UR){T}(start, stop)
        $(UR)(r::AbstractUnitRange) = $(UR)(first(r), last(r))

        $(UR){T}(r::$(UR){T}) where {T<:Real} = r
        $(UR){T}(r::$(UR)) where {T<:Real} = $(UR){T}(first(r), last(r))
    end
end


text:
module Input

using ..FEStructure
using ..LoadCase

function parse_add_node(s::String)
    words=split(a,",")
    id=words[2]
    x=words[3]
    y=words[4]
    z=words[5]
    add_node!(structure, id, x, y, z)
end

function parse_add_uniaxial_metal(s::String)
    add_uniaxial_metal!(structure,id,E,ν,ρ,α)
end

function parse_add_general_section!(s::String)
end

function parse_add_beam_section!(s::String)
end

function parse_add_node!(s::String)
end

function parse_set_nodal_mass!(s::String)
end

function set_nodal_restraint!(s::String)
end

function set_nodal_spring!(s::String)
end

function add_beam!(s::String)
end

function set_beam_release!(s::String)
end

function set_beam_orient!(s::String)
end

function add_quad!(s::String)
end

function set_quad_rotation!(s::String)
end

function add_tria!(s::String)
end

function set_tria_rotation!(s::String)
end

function add_static_case!(s::String)
end

function add_modal_case!(s::String)
end

function add_buckling_case!(s::String)
end

function add_time_history_case!(s::String)
end

function add_nodal_force!(s::String)
end

function add_beam_strain!(s::String)
end

function add_beam_distributed!(s::String)
end

function set_modal_params!(s::String)
end

function set_buckling_params!(s::String)
end

function set_time_history_params!(s::String)
end


text:
# ---
# title: 173. Binary Search Tree Iterator
# id: problem173
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Stack, Tree, Design
# link: <https://leetcode.com/problems/binary-search-tree-iterator/description/>
# hidden: true
# ---
#
# Implement the `BSTIterator` class that represents an iterator over the **[in-
# order traversal](https://en.wikipedia.org/wiki/Tree_traversal#In-
# order_\(LNR\))**  of a binary search tree (BST):
#
#   * `BSTIterator(TreeNode root)` Initializes an object of the `BSTIterator` class. The `root` of the BST is given as part of the constructor. The pointer should be initialized to a non-existent number smaller than any element in the BST.
#   * `boolean hasNext()` Returns `true` if there exists a number in the traversal to the right of the pointer, otherwise returns `false`.
#   * `int next()` Moves the pointer to the right, then returns the number at the pointer.
#
# Notice that by initializing the pointer to a non-existent smallest number, the
# first call to `next()` will return the smallest element in the BST.
#
# You may assume that `next()` calls will always be valid. That is, there will
# be at least a next number in the in-order traversal when `next()` is called.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2018/12/25/bst-tree.png)
#
#
#
#     **Input**
#     ["BSTIterator", "next", "next", "hasNext", "next", "hasNext", "next", "hasNext", "next", "hasNext"]
#     [[[7, 3, 15, null, null, 9, 20]], [], [], [], [], [], [], [], [], []]
#     **Output**
#     [null, 3, 7, true, 9, true, 15, true, 20, false]
#
#     **Explanation**
#     BSTIterator bSTIterator = new BSTIterator([7, 3, 15, null, null, 9, 20]);
#     bSTIterator.next();    // return 3
#     bSTIterator.next();    // return 7
#     bSTIterator.hasNext(); // return True
#     bSTIterator.next();    // return 9
#     bSTIterator.hasNext(); // return True
#     bSTIterator.next();    // return 15
#     bSTIterator.hasNext(); // return True
#     bSTIterator.next();    // return 20
#     bSTIterator.hasNext(); // return False
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[1, 105]`.
#   * `0 <= Node.val <= 106`
#   * At most `105` calls will be made to `hasNext`, and `next`.
#
#
#
# **Follow up:**
#
#   * Could you implement `next()` and `hasNext()` to run in average `O(1)` time and use `O(h)` memory, where `h` is the height of the tree?
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
Determines if of Image Type
"""
function isa_image{T<:Matrix}(x::Type{T})
    eltype(T) <: Union{Colorant, AbstractFloat}
end
isa_image(x::Matrix) = isa_image(typeof(x))
isa_image(x::Images.Image) = true
isa_image(x) = false

# Splits a dictionary in two dicts, via a condition
function Base.split(condition::Function, associative::Associative)
    A = similar(associative)
    B = similar(associative)
    for (key, value) in associative
        if condition(key, value)
            A[key] = value
        else
            B[key] = value
        end
    end
    A, B
end

function assemble_robj(data, program, bb, primitive, pre_fun, post_fun)
    pre = if pre_fun != nothing
        () -> (GLAbstraction.StandardPrerender(); pre_fun())
    else
        GLAbstraction.StandardPrerender()
    end
    robj = RenderObject(data, program, pre, nothing, bb, nothing)
    post = if haskey(data, :instances)
        GLAbstraction.StandardPostrenderInstanced(data[:instances], robj.vertexarray, primitive)
    else
        GLAbstraction.StandardPostrender(robj.vertexarray, primitive)
    end
    robj.postrenderfunction = if post_fun != nothing
        () -> begin
            post()
            post_fun()
        end
    else
        post
    end
    robj
end

function assemble_shader(data)
    shader = data[:shader]
    delete!(data, :shader)
    default_bb = Signal(centered(AABB))
    bb  = get(data, :boundingbox, default_bb)
    if bb == nothing || isa(bb, Signal{Void})
        bb = default_bb
    end
    glp = get(data, :gl_primitive, GL_TRIANGLES)
    robj = assemble_robj(
        data, shader, bb, glp,
        get(data, :prerender, nothing),
        get(data, :postrender, nothing)
    )
    Context(robj)
end

function y_partition(area, percent)
    amount = percent / 100.0
    p = const_lift(area) do r
        (SimpleRectangle{Int}(r.x, r.y, r.w, round(Int, r.h*amount)),
            SimpleRectangle{Int}(r.x, round(Int, r.h*amount), r.w, round(Int, r.h*(1-amount))))
    end
    return const_lift(first, p), const_lift(last, p)
end
function x_partition(area, percent)
    amount = percent / 100.0
    p = const_lift(area) do r
        (SimpleRectangle{Int}(r.x, r.y, round(Int, r.w*amount), r.h ),
            SimpleRectangle{Int}(round(Int, r.w*amount), r.y, round(Int, r.w*(1-amount)), r.h))
    end
    return const_lift(first, p), const_lift(last, p)
end

glboundingbox(mini, maxi) = AABB{Float32}(Vec3f0(mini), Vec3f0(maxi)-Vec3f0(mini))
function default_boundingbox(main, model)
    main == nothing && return Signal(AABB{Float32}(Vec3f0(0), Vec3f0(1)))
    const_lift(*, model, AABB{Float32}(main))
end
@compat (::Type{AABB})(a::GPUArray) = AABB{Float32}(gpu_data(a))
@compat (::Type{AABB{T}}){T}(a::GPUArray) = AABB{T}(gpu_data(a))

"""
Returns two signals, one boolean signal if clicked over `robj` and another
one that consists of the object clicked on and another argument indicating that it's the first click
"""
function clicked(robj::RenderObject, button::MouseButton, window::Screen)
    @materialize mouse_hover, mouse_buttons_pressed = window.inputs
    leftclicked = const_lift(mouse_hover, mouse_buttons_pressed) do mh, mbp
        mh[1] == robj.id && mbp == Int[button]
    end
    clicked_on_obj = keepwhen(leftclicked, false, leftclicked)
    clicked_on_obj = const_lift((mh, x)->(x,robj,mh), mouse_hover, leftclicked)
    leftclicked, clicked_on_obj
end

"""
Returns a boolean signal indicating if the mouse hovers over `robj`
"""
is_hovering(robj::RenderObject, window::Screen) =
    droprepeats(const_lift(is_same_id, mouse2id(window), robj))

function dragon_tmp(past, mh, mbp, mpos, robj, button, start_value)
    diff, dragstart_index, was_clicked, dragstart_pos = past
    over_obj = mh[1] == robj.id
    is_clicked = mbp == Int[button]
    if is_clicked && was_clicked # is draggin'
        return (dragstart_pos-mpos, dragstart_index, true, dragstart_pos)
    elseif over_obj && is_clicked && !was_clicked # drag started
        return (Vec2f0(0), mh[2], true, mpos)
    end
    return start_value
end

"""
Returns a signal with the difference from dragstart and current mouse position,
and the index from the current ROBJ id.
"""
function dragged_on(robj::RenderObject, button::MouseButton, window::Screen)
    @materialize mouse_buttons_pressed, mouseposition = window.inputs
    mousehover = mouse2id(window)
    mousedown = const_lift(GLAbstraction.singlepressed, mouse_buttons_pressed, GLFW.MOUSE_BUTTON_LEFT)
    condition = const_lift(is_same_id, mousehover, robj)
    dragg = GLAbstraction.dragged(mouseposition, mousedown, condition)
    filterwhen(mousedown, (value(dragg), 0), map(dragg) do d
        d, value(mousehover).index
    end)
end

points2f0{T}(positions::Vector{T}, range::Range) = Point2f0[Point2f0(range[i], positions[i]) for i=1:length(range)]

extrema2f0{T<:Intensity,N}(x::Array{T,N}) = Vec2f0(extrema(reinterpret(Float32,x)))
extrema2f0{T,N}(x::Array{T,N}) = Vec2f0(extrema(x))
extrema2f0(x::GPUArray) = extrema2f0(gpu_data(x))
function extrema2f0{T<:Vec,N}(x::Array{T,N})
    _norm = map(norm, x)
    Vec2f0(minimum(_norm), maximum(_norm))
end

"""
Converts index arrays to the OpenGL equivalent.
"""
to_indices(x::GLBuffer) = x
to_indices(x::TOrSignal{Int}) = x
to_indices(x::VecOrSignal{UnitRange{Int}}) = x
"""
For integers, we transform it to 0 based indices
"""
to_indices{I<:Integer}(x::Vector{I}) = indexbuffer(map(i-> Cuint(i-1), x))
function to_indices{I<:Integer}(x::Signal{Vector{I}})
    x = map(x-> Cuint[i-1 for i=x], x)
    gpu_mem = GLBuffer(value(x), buffertype = GL_ELEMENT_ARRAY_BUFFER)
    preserve(const_lift(update!, gpu_mem, x))
    gpu_mem
end
"""
If already GLuint, we assume its 0 based (bad heuristic, should better be solved with some Index type)
"""
to_indices{I<:GLuint}(x::Vector{I}) = indexbuffer(x)
function to_indices{I<:GLuint}(x::Signal{Vector{I}})
    gpu_mem = GLBuffer(value(x), buffertype = GL_ELEMENT_ARRAY_BUFFER)
    preserve(const_lift(update!, gpu_mem, x))
end
to_indices(x) = error(
    "Not a valid index type: $x.
    Please choose from Int, Vector{UnitRange{Int}}, Vector{Int} or a signal of either of them"
)

function mix_linearly{C<:Colorant}(a::C, b::C, s)
    RGBA{Float32}((1-s)*comp1(a)+s*comp1(b), (1-s)*comp2(a)+s*comp2(b), (1-s)*comp3(a)+s*comp3(b), (1-s)*alpha(a)+s*alpha(b))
end

color_lookup(cmap, value, mi, ma) = color_lookup(cmap, value, (mi, ma))
function color_lookup(cmap, value, color_norm)
    mi,ma = color_norm
    scaled = clamp((value-mi)/(ma-mi), 0, 1)
    index = scaled * (length(cmap)-1)
    i_a, i_b = floor(Int, index)+1, ceil(Int, index)+1
    mix_linearly(cmap[i_a], cmap[i_b], scaled)
end


text:
@testset "PCA reduction" begin
    wv = fake_wordvectors()
	for d in [1, 10]
		for outdim in [1, 10]
            for do_pca in [false, true]
                wv_pca = pca_reduction(wv, d, outdim, do_pca=do_pca)
                @test wv isa typeof(wv)
            end
		end
	end
end

max_stars_repo_name: zgornel/EmbeddingsAnalysis

max_issues_repo_name: zgornel/EmbeddingsAnalysis

max_forks_repo_name: zgornel/EmbeddingsAnalysis


text:
module gss

function goldss(f,a,b;ϵ=1e-6,τ=0.618)
    fa,fb=f(a),f(b)
    x1,x2=b-(τ)*(b-a),a+(τ)*(b-a)
    err=Inf
    while err>ϵ
        if f(x2)>f(x1)
            b=x2
            x2=x1
            x1=a+(1-τ)*(b-a)
        elseif f(x2)<f(x1)
            a=x1
            x1=x2
            x2=b-(1-τ)*(b-a)
        else
            a=(x1+x2)/2
            b=a
        end
        err=2*abs((b-a)/(b+a))
    end
    xsr=(x1+x2)/2
    wart=f(xsr)
    return(round(xsr,digits=2),round(wart,digits=3))
end

export goldss

end # module


text:
module Forward

import ..Zygote
import ..Zygote: __new__, __splatnew__

export pushforward

include("compiler.jl")
include("interface.jl")
include("lib.jl")
include("number.jl")
include("array.jl")
include("broadcast.jl")

end


text:
## ----------------  Bounded learning

"""
    $(SIGNATURES)

For a single college (not a ModelObject).
Each college is endowed with `maxLearn`. Once a student has learned this much, learning productivity falls to 0 (or a constant).

`dh = exp(aScale * a) * studyTime ^ timeExp * A`

The functional form for `A` is governed by the `tfpSpec`. It depends on how much has been learned.

The way learning is defined is governed by `learnRelativeToH0`.
- false: `h learned = h - h0`
- true: `h learned = (h / h0 - 1)`. Then a college limits the percentage increase in the h endowment.

Options should be type parameters +++

A potential alternative with better scaling would be
`A = hExp .* ( log(maxLearn) .- log.(max.(1.0, h learned)) )`
`hExp` governs the slope. `maxLearn` governs the intercept. But the shape is fixed.
"""
mutable struct HcProdFctBounded <: AbstractHcProdFct
    minTfp :: Double  # Additive minimum tfp
    tfp :: Double
    maxLearn  ::  Double
    timeExp  ::  Double
    # Curvature: how strongly does learning decline as (h-h0) → maxLearn
    hExp :: Double
    # Depreciation rate
    deltaH :: Double
    # Ability scale
    aScale  ::  Double
    # Fixed time cost per course
    timePerCourse :: Double
    # Study time per course minimum (this is assigned when study time very low)
    minTimePerCourse :: Double
    # Learning as percentage of endowment? Or as (h - h0).
    learnRelativeToH0 :: Bool
    # TFP can be computed in several ways. See `base_tfp`.
    tfpSpec :: AbstractTfpSpec
end

## -------------  All colleges

Base.@kwdef mutable struct HcProdBoundedSwitches <: AbstractHcProdSwitches
    # Same exponents on time and h. If yes, ignore `hExp`.
    minTfp :: Double = 1.0
    calMinTfp :: Bool = true
    tfp :: Double = 1.0
    calTfpBase :: Bool = true
    sameExponents :: Bool = true
    timeExp :: Double = 0.6
    calTimeExp :: Bool = true
    hExp :: Double = 0.9
    # hExpLb :: Double = 0.5
    calHExp :: Bool = true
    deltaH :: Double = 0.0
    calDeltaH :: Bool = false
    aScale :: Double = 0.2
    calAScale :: Bool = true
    # Learning as percentage of endowment?
    learnRelativeToH0 :: Bool = false
    # TFP from (max learning - learning)
    tfpSpec :: AbstractTfpSpec = TfpMaxLearnMinusLearn()
end

"""
    $(SIGNATURES)

Since all colleges share some parameters, we need a model object that keeps
track of parameters that are common or differ by college.
"""
mutable struct HcProdBoundedSet <: AbstractHcProdSet
    objId :: ObjectId
    switches :: HcProdBoundedSwitches
    nc :: CollInt

    # Calibrated parameters
    minTfp :: Double
    tfp :: Double
    maxLearnV  ::  BoundedVector
    timeExp  ::  Double
    hExp  ::  Double
    deltaH :: Double
    # Ability scale
    aScale  ::  Double

    # Fixed time cost per course
    timePerCourse :: Double
    # Study time per course minimum (this is assigned when study time very low)
    minTimePerCourse :: Double

    pvec :: ParamVector
end

## H production: Bounded learning

max_learn(hs :: HcProdBoundedSet, ic) = ModelParams.values(hs.maxLearnV, ic);
max_learn(h :: HcProdFctBounded) = h.maxLearn;

learning_relative_to_h0(h :: HcProdFctBounded) = h.learnRelativeToH0;
learning_relative_to_h0(h :: HcProdBoundedSwitches) = h.learnRelativeToH0;
learning_relative_to_h0(h :: HcProdBoundedSet) = learning_relative_to_h0(h.switches);

tfp_spec(h :: HcProdFctBounded) = h.tfpSpec;
tfp_spec(h :: HcProdBoundedSwitches) = h.tfpSpec;
tfp_spec(h :: HcProdBoundedSet) = tfp_spec(h.switches);

## ----------  Construction

# Initialize with defaults
function make_hc_prod_set(objId :: ObjectId, nc :: Integer,
    switches :: HcProdBoundedSwitches)

    st = symbol_table(); # eventually use preconstructed +++
    @assert validate_hprod(switches);

    pTimePerCourse = init_time_per_course();
    pMaxLearn = init_max_learn(objId, switches, nc);

    minTfp = switches.minTfp;
    pMinTfp = Param(:minTfp, "Min tfp", "A_{min}",
        minTfp, minTfp, 0.0, 2.0, switches.calMinTfp);

    tfpBase = switches.tfp;
    pTfpBase = Param(:tfp, ldescription(:hTfpNeutral), lsymbol(:hTfpNeutral),
        tfpBase, tfpBase, 0.1, 2.0, switches.calTfpBase);

    timeExp = switches.timeExp;
    pTimeExp = Param(:timeExp, ldescription(:hTimeExp), lsymbol(:hTimeExp),
        timeExp, timeExp, 0.2, 0.9, switches.calTimeExp);

    deltaH = delta_h(switches);
    pDeltaH = Param(:deltaH, ldescription(:ddh), lsymbol(:ddh),
        deltaH, deltaH, 0.0, 0.5, cal_delta_h(switches));

    # Governs slope inside of TFP (should be inside of TFP spec +++)
    hExp = switches.hExp;
    pHExp = Param(:hExp, "TFP slope coefficient", lsymbol(:hHExp),
        hExp, hExp, gma_range(tfp_spec(switches))..., switches.calHExp);

    aScale = switches.aScale;
    tfpSpec = tfp_spec(switches);
    pAScale = Param(:aScale, ldescription(:hAScale), lsymbol(:hAScale),
        aScale, aScale, gma_range(tfpSpec)..., switches.calAScale);

    pvec = ParamVector(objId,  [pMinTfp, pTfpBase, pTimeExp, pHExp, pDeltaH, pAScale, pTimePerCourse]);
    # Min study time required per course. Should never bind.
    minTimePerCourse =
        hours_per_week_to_mtu(0.1 / data_to_model_courses(1));

    h = HcProdBoundedSet(objId, switches, nc,
        minTfp, tfpBase, pMaxLearn,
        timeExp, hExp, deltaH, aScale,
        pTimePerCourse.value, minTimePerCourse, pvec);
    @assert validate_hprod_set(h)
    return h
end

# Upper bound should depend on whether learning is relative to h0.
function init_max_learn(objId :: ObjectId, switches, nc :: Integer)
    ownId = make_child_id(objId, :tfpV);

    dMaxLearnV = fill(0.2, nc);
    if learning_relative_to_h0(switches)
        ub = 3.0;
    else
        ub = 5.0;
    end
    b = BoundedVector(ownId, ParamVector(ownId), :increasing, 0.2, ub, dMaxLearnV);
    set_pvector!(b; description = ldescription(:maxLearn),
        symbol = lsymbol(:maxLearn));
    return b
end

make_test_hc_bounded_set(; learnRelativeToH0 = true,
    tfpSpec = TfpMaxLearnMinusLearn()) =
    make_hc_prod_set(ObjectId(:HProd), 4,
        HcProdBoundedSwitches(
            deltaH = 0.05,
            learnRelativeToH0 = learnRelativeToH0,
            tfpSpec = tfpSpec
            ));

# Make h production function for one college
function make_h_prod(hs :: HcProdBoundedSet, iCollege :: Integer)
    return HcProdFctBounded(hs.minTfp, hs.tfp,
        max_learn(hs, iCollege),
        time_exp(hs), h_exp(hs),
        delta_h(hs), hs.aScale,
        hs.timePerCourse,  hs.minTimePerCourse,
        learning_relative_to_h0(hs), tfp_spec(hs));
end

function make_test_hprod_bounded(;
    learnRelativeToH0 = true, tfpSpec = TfpMaxLearnMinusLearn())

    minTfp = 0.7;
    gma = sum(gma_range(tfpSpec)) / 2;
    hS = HcProdFctBounded(minTfp, 0.6, 3.1, 0.7, gma, 0.1, 0.3, 0.01, 0.005,
        learnRelativeToH0, tfpSpec);
    @assert validate_hprod(hS);
    return hS
end

## ----------  One college

function validate_hprod(hS :: HcProdFctBounded)
    isValid = (max_learn(hS) > 0.05)  &&  (0.0 < time_exp(hS) ≤ 1.0);
    gmaMin, gmaMax = gma_range(tfp_spec(hS));
    gma = h_exp(hS);
    isValid = isValid  &&  (gmaMin <= gma <= gmaMax);
    return isValid
end

"""
    $(SIGNATURES)

H produced (before shock is realized). Nonnegative.

# Arguments
- nTriedV
    number of courses attempted this period.
- h0V
    h endowments, so that `hV - h0V` is learning.
"""
function dh(hS :: HcProdFctBounded, abilV,  hV, h0V, timeV, nTriedV)
    sTimeV = study_time_per_course(hS, timeV, nTriedV);
    # deltaHV = (max_learn(hS) ^ h_exp(hS) .- learned_h(hS, hV, h0V) .^ h_exp(hS));
    # tfpV = hS.tfp .* max.(0.0, deltaHV) .^ (1.0 / h_exp(hS));
    return nTriedV .* base_tfp(hS, hV, h0V) .* (sTimeV .^ hS.timeExp) .*
        exp.(hS.aScale .* abilV);
end

## ----------  TFP specs

# Base TFP: the term in front of (sTime ^ beta * exp(ability))
function base_tfp(hS :: HcProdFctBounded, hV, h0V)
    tfpSpec = tfp_spec(hS);
    learnV = learned_h(hS, hV, h0V);
    tfpV = hS.minTfp .+ hS.tfp .* tfp(tfpSpec, learnV, max_learn(hS), h_exp(hS));
    return tfpV
end

# Expected range of TFP
function tfp_range(hS :: HcProdFctBounded)
    return hS.minTfp .+ tfp(hS) .*
        tfp_range(tfp_spec(hS), h_exp(hS), max_learn(hS));
end

# Learned h, scaled for the production function
function learned_h(hS :: HcProdFctBounded, hV, h0V)
    if learning_relative_to_h0(hS)
        dh = max.(0.0, hV .- h0V) ./ h0V;
    else
        dh = max.(0.0, hV .- h0V);
    end
    return dh
end

# function show_string(hS :: HcProdFctBounded)
#     fs = Formatting.FormatExpr("dh = {1:.2f}  h ^ {2:.2f}  t ^ {3:.2f}  exp({3:.2f} a)");
#     return format(fs,   hS.tfp, h_exp(hS), hS.timeExp, hS.aScale);
# end

function Base.show(io :: IO, hS :: HcProdFctBounded)
    maxLearn = round(max_learn(hS), digits = 2);
    print(io,  "H prod fct:  Bounded learning < $maxLearn");
end

## ---------------------  For all colleges

function Base.show(io :: IO, switches :: HcProdBoundedSwitches)
    print(io, "H production: bounded learning.");
end

function settings_table(h :: HcProdBoundedSwitches)
    ddh = delta_h(h);
    cal_delta_h(h)  ?  deprecStr = "calibrated"  :  deprecStr = "fixed at $ddh";
    h.learnRelativeToH0  ?  learnStr = "h/h0 - 1"  :  learnStr = "h - h0";
    ownSettings = [
        "H production function" "Bounded learning";
        "Depreciation"  deprecStr
        "Learning of the form"  learnStr
        ];
    return vcat(ownSettings, settings_table(h.tfpSpec))
end

function settings_list(h :: HcProdBoundedSwitches, st)
    eqnHChange = ["H production function", "eqnHChange",  eqn_hchange(h)];
    return [eqnHChange]
end

function validate_hprod(s :: HcProdBoundedSwitches)
    isValid = true;
    return isValid
end

function validate_hprod_set(h :: HcProdBoundedSet)
    isValid = (h.nc > 1)  &&  (h.timeExp > 0.0)  &&
        (h.aScale > 0.0)  &&  (h.timePerCourse > 0.0);
    isValid = isValid  &&  (1.0 > delta_h(h) >= 0.0);
    return isValid
end

function eqn_hchange(h :: HcProdBoundedSwitches)
    "\\hTfp \\sTimePerCourse^{\\hTimeExp} e^{\\hAScale \\abil}"
end

# --------------


text:
using Soss, MeasureTheory

hmm = @model begin
    ε ~ Exponential() #  transition
    σ ~ Exponential() # Observation noise
    x ~ Chain(Normal()) do xj
        Normal(xj, ε)
    end

    y ~ For(x) do xj
        Normal(xj, σ)
    end
end

using Random

_rng = Random.GLOBAL_RNG

# sourceRand(hmm)

σ = rand(_rng, Exponential())
ε = rand(_rng, Exponential())
x = rand(_rng, Chain(Normal()) do xj
            #= REPL[2]:5 =#
            Normal(xj, ε)
        end)
y = rand(_rng, For(x) do xj
            #= REPL[2]:8 =#
            Normal(xj, σ)
        end)

rand(hmm())


text:
module TestEquations

using Modia
using ModiaMath: plot

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
    using LinearAlgebra
    eye(n) = Matrix{Float64}(I, n, n)
    diagm(v) = Matrix(Diagonal(v))
    q(qr) = qr.Q
    r(qr) = qr.R
end

# @testset "Equations" begin

@model Test begin
    a = -1; b = 1; c = 1
    n = 10;
    M = Float()
    u = Float(); x = Float(start=0.0)
    on = Variable(); s = Variable()
    QR = Float(); Q = Float(); R = Float()
    dummy = Variable(size=())
    @equations begin
        u = if s == "on"; 10 else M end # Varying size of u
        der(x) - a * x = b * (s == "on" ? u : 10 * u[1,1])
        on = time > 1 && time < 3
        s := on ? "on" : "off"
        dummy = if time < 0.00001 || abs(time - 2.0) < 0.01; println(" Time=", time, ":  size of u: ", size(u)) else nothing end
        ((9 + 1) * diagm(1:n) + ones(n, n))^3 * M * diagm(1:n) = eye(n)
#        QR = qr(M)
#       Q = QR.Q; R = QR.R # ERROR: LoadError: No substitution provided for: this.QR.Q
#       Q = q(QR); R = r(QR) # ERROR: LoadError: MethodError: no method matching length(::LinearAlgebra.QRCompactWY{Float64,Array{Float64,2}})
        QR = Tuple{Array,Array}(qr(M)) # Tuple{Array,Array}() needed for Julia 0.7
        Q = QR[1]; R = QR[2]
    end
end;

result = simulate(Test, 5, storeEliminated=false, logTranslation=true, removeSingularities=false)
plot(result, ("x"), heading="TestEquations", figure=2)
@test result["x"][end] == 1.1776785083961023

end


text:
using PyPlot
using IterativeSolvers

include("../src/FastConvolution.jl")
include("../src/Preconditioner.jl")

# setting the number of threads for the FFT and BLAS
# libraries (please set them to match the number of
# physical cores in your system)
FFTW.set_num_threads(4);
BLAS.set_num_threads(4);

#Defining Omega
hInit = 1/8;
nInit = 8;
k = 18.0
uRho = []

sigma = 0.05
renorm_factor = 1/(sigma^2*(2*pi));
rho(x,y) = renorm_factor*exp((-1/(2*sigma^2))*(x.^2 + y.^2 ));

nu(x,y) = 0.3*exp(-40*(x.^2 + y.^2)).*(abs(x).<0.48).*(abs(y).<0.48);

nSamples = 8

for ii = 1:nSamples

    #h even
    h = hInit/2.0^((ii-1))
    # h odd
    #h = 1/(nInit*2.0^((ii-1))+1)

    # size of box
    a = 1
    # x = collect(-a/2:h:a/2)
    # y = collect(-a/2:h:a/2)
    x = collect(-a/2:h:a/2-h)
    y = collect(-a/2:h:a/2-h)
    (n,m) = length(x), length(y)
    println(n)
    N = n*m
    X = repmat(x, 1, m)[:]
    Y = repmat(y', n,1)[:]
    # we solve \triangle u + k^2(1 + nu(x))u = 0

    Rho = rho(X,Y) + 0*1im;
    RHO = reshape(Rho, n,m);

    figure(5); clf();
    imshow(real(RHO))

    # Greengard Vico Quadrature (this is not optimized and is 2-3 times slower)
    fastconv = buildFastConvolution(x,y,h,k,nu, quadRule = "Greengard_Vico");

    u_rho = FFTconvolution(fastconv, Rho);
    u_rho = reshape(u_rho,n,m)

    u_inc = exp(k*im*X);
    u_conv = FFTconvolution(fastconv, u_inc );
    u_conv = reshape(u_conv,n,m)

    push!(uRho, u_rho[1:round(Integer,2.0^((ii-1))):end, 1:round(Integer,2.0^((ii-1))):end])

end

error = zeros(nSamples-1)
for ii = 1:nSamples-1
    error[ii] = norm( uRho[end][:] - uRho[ii][:])/norm(uRho[end][:])
end

println(error)

max_stars_repo_name: Forgotten/Lippmann-Schingwer

max_issues_repo_name: Forgotten/Lippmann-Schingwer

max_forks_repo_name: Forgotten/Lippmann-Schingwer


text:
include("SimulationsUtils.jl")
include("MethodUtils.jl")
include("WassersteinVariance.jl")
using Serialization

function time_power_wvm_icp(nsim, numS, ne=500, noisetype="normal", seed=123)
    Random.seed!(seed)

    ######################## SIMULATION ##############################
    # These will store a four tuple for each simulation: nfp, nfn, length(S), length(Ŝ)
    total_icp_errs, total_time_icp = [], []
    # total_wvm_errs, total_time_wvm = [], []

    simtype = "numS-$(numS)/ne-$(ne)/$(noisetype)"
    for i in 1:nsim
        println("--------- Sim $(i) ----------")
        X_full, S, nₑs = deserialize("../data/Simulations/$(simtype)/sim-$(i).jls")
        p = size(X_full,2)
        X = X_full[:, 1:end-1]
        y = X_full[:, end]
        n = size(X, 1)

        icp_errs, time_icp = [], []
        # wvm_errs, time_wvm = [], []
        for numvar in 2:2:30
            print("------ Number of Selected Vars: $(numvar) -------")
            ################# lasso #################
            cf = select_lasso_by_num_var(X, y, numvar)
            select = filter(x -> cf[x] != 0, 1:p-1)
            if length(select) <= 1; continue end
            X_select = X[:,select]

            ################# ICP #################
            α = 0.1
            if numvar <= 18
                t_icp = @elapsed begin
                    Ŝ_icp = select[icp(X_select, y, nₑs, α, false, "all")]
                end
                push!(icp_errs, error_count(S, Ŝ_icp))
                push!(time_icp, t_icp)
            end

            ################# WVM #################
            t_wvm = @elapsed begin
                wv = DirectWassersteinVariance(X_select, y, nₑs, nₑs / n)
                res = WV_screening(wv, α=α, solver=LBFGS(m=100), B=50)
                Ŝ_wvm = select[findall_nonzero(res[2])]
            end
            push!(wvm_errs, error_count(S, Ŝ_wvm))
            push!(time_wvm, t_wvm)

            println("###### SIMULATION $i IS DONE ! ######")
        end
        push!(total_icp_errs, icp_errs)
        push!(total_time_icp, time_icp)

        push!(total_wvm_errs, wvm_errs)
        push!(total_time_wvm, time_wvm)
    end
    serialize("../data/results/$(simtype)/icp-time-power.jls", (total_icp_errs, total_time_icp))
    serialize("../data/results/$(simtype)/wvm-time-power.jls", (total_wvm_errs, total_time_wvm))
end
time_power_wvm_icp(100, 6)


text:
#=
Aalens additive model
=#

function aalen(X, t, status; entry=[], id=[], beta=[], opt...)

end
(


text:
export upgrade_efficiency, vegetation, nutrients

# Pre-defined helpers for data collection
nutrients(model) = model.lake.p.nutrients

# A couple of different vegetation extractors, depending on the shape of your data
# Good for `model.lake`
function vegetation(
    lake::OrdinaryDiffEq.ODEIntegrator{A,B,C,D,E,SchefferParameters},
) where {A,B,C,D,E}
    bream = lake.sol[1, :]
    @. lake.p.K * (lake.p.H₃^2 / (lake.p.H₃^2 + bream^2))
end

# Good for a singular bream vector
function vegetation(B::Vector{Float64}, p::SchefferParameters)
    @. p.K * (p.H₃^2 / (p.H₃^2 + B^2))
end

# TODO: This should now be at the municipality level, not WC.
function upgrade_efficiency(m::ABM)
    if m.year > m.outcomes_year_when_informed > 0
        if 0 < m.outcomes_year_of_full_upgrade < m.year
            100.0
        else
            t1 = m.outcomes_year_when_informed
            t2 = m.outcomes_year_of_full_upgrade > 0 ? m.outcomes_year_of_full_upgrade :
                m.year
            # Calculate current efficiency (with optimal efficiency in denominator)
            (m.outcomes_upgraded_households_sum * 100.0) / ((t2 - t1) * nagents(m))
        end
    else
        0.0
    end
end


text:
function counting_sort_letters(A, pos) # A:Array, p:position
    n = length(A)
    k = 26 # number of letters in alphabet
    start_ascii = Int('a'[1]) # Ascii decimal of a
    C = fill(0,k) # initialize counter array to 0
    # Storing count for each letter
    for i = 1:n
        if (length(A[i]) < pos) return A end
        counterIdx = Int(A[i][pos])-start_ascii+1 # -a_ascii+1 to start from 1
        C[counterIdx] = C[counterIdx] + 1
    end
    # Changing C[j] so it contains actual position of all each element
    for j = 2:(k)
        C[j] = C[j] + C[j-1]
    end
    B = fill("z",n) # initialize return array
    for i = n:-1:1
        counterIdx = Int(A[i][pos])-start_ascii+1
        B[C[counterIdx]] = A[i] # Always adds next string to first element of B
        C[counterIdx] = C[counterIdx] - 1 # Reduce the count of current element by 1
    end
    return B
end

function counting_sort_length(A)
    n = length(A)
    k = 0 # Max string length in array A
    for i = 1:n # find the longest string
        if (length(A[i]) > k) k = length(A[i]) end
    end
    C = fill(0,k) # initialize counter array to 0
    for i = 1:n # counter
        counterIdx = length(A[i])
        C[counterIdx] = C[counterIdx] + 1
    end
    # Changing C[j] so it contains actual position of all each element
    for j = 2:(k)
        C[j] = C[j] + C[j-1]
    end
    B = fill("z",n) # initialize return array
    for i = n:-1:1
        counterIdx = length(A[i])
        B[C[counterIdx]] = A[i]  # Always adds next string to first element of B
        C[counterIdx] = C[counterIdx] - 1 # Reduce the count of current element by 1
    end
    return B
end

function flexradix(A, max_length)
    A = counting_sort_length(A)

    # SPLIT UP ARRAY AND RETURN ONLY THE SIDE WHERE EVERY ELEMENT HAS
    # STRING SIZE BIGGER THAN P??

    for p = 1:max_length # p: position
        A = counting_sort_letters(A, p)
    end
    return A

    #=
    n = length(A)
    C = fill(0,max_length) # initialize counter array to 0
    for i = 1:n # counter
        counterIdx = length(A[i])
        C[counterIdx] = C[counterIdx] + 1
    end
    # Now we have a sorted list by length, and the number of lengths. We can then
    # iterate through the lengths and ommit the previous smaller arrays as
    # the problem developes
    for i = 1:max_length
        A =
        C = fill(0,n) # Counter for strings with different sizes
        A = counting_sort_letters(A, i)
        C[counterIdx] = C[counterIdx] - 1 # Reduce the count of current element by 1
    end
    =#
end

max_stars_repo_name: evenlwanvik-student/algdat

max_issues_repo_name: evenlwanvik-student/algdat

max_forks_repo_name: evenlwanvik-student/algdat


text:
for func in (:daily, :weekly, :monthly)
    x = "digital_currency_$(func)"
    fname = Symbol(x)
    @eval begin
        function ($fname)(symbol::String, market::String="USD"; outputsize::String="compact", datatype::String="json")
            @argcheck in(outputsize, ["compact", "full"])
            @argcheck in(datatype, ["json", "csv"])
            uri = "$(alphavantage_api)query?function="* uppercase($x) * "&symbol=$symbol&market=$market&outputsize=$outputsize&datatype=$datatype&apikey=" * ENV["ALPHA_VANTAGE_API_KEY"]
            data = _get_request(uri)
            return _parse_response(data, datatype)
        end
        export $fname
    end
end

function crypto_rating(symbol::String)
    uri = "$(alphavantage_api)query?function=CRYPTO_RATING&symbol=$symbol&apikey=" * ENV["ALPHA_VANTAGE_API_KEY"]
    data = _get_request(uri)
    return _parse_response(data, "json")
end
export crypto_rating


text:
abstract type AbstractIntegratorVPRK{DT,TT,D,S} <: AbstractIntegratorIRK{DT,TT} end
abstract type AbstractIntegratorVPRKwProjection{DT,TT,D,S} <: AbstractIntegratorVPRK{DT,TT,D,S} end

@inline parameters(integrator::AbstractIntegratorVPRK) = integrator.params

@inline nstages(integrator::AbstractIntegratorVPRK{DT,TT,D,S}) where {DT,TT,D,S} = S
@inline Base.ndims(integrator::AbstractIntegratorVPRK{DT,TT,D,S}) where {DT,TT,D,S} = D

@inline eachstage(integrator::AbstractIntegratorVPRK) = 1:nstages(integrator)
@inline eachdim(integrator::AbstractIntegratorVPRK) = 1:ndims(integrator)


text:
# test new sampling API

using Test
using IncrementalInference
using DistributedFactorGraphs

# import IncrementalInference: getSample

struct SpecialPrior{T} <: AbstractPrior where T <: Distribution
  z::T
  specialSampler::Function
end
specialSample(s::SpecialPrior, N::Int, fmd::FactorMetadata, Xi) = (reshape(rand(s.z,N),1,:), )
SpecialPrior(z::T) where {T <: SamplableBelief}= SpecialPrior{T}(z, specialSample)

struct SpecialLinearOffset{T <: SamplableBelief} <: AbstractRelativeRoots
  z::T
  specialSampler::Function
end

specialSample(s::SpecialLinearOffset, N::Int, fmd::FactorMetadata, Xi, Xj) = (reshape(rand(s.z,N),1,:), )

SpecialLinearOffset(z::T) where {T <: SamplableBelief} = SpecialLinearOffset{T}(z, specialSample)

function (s::SpecialLinearOffset)(res::Array{Float64},
                                  userdata::FactorMetadata,
                                  idx::Int,
                                  meas::Tuple{Array{Float64, 2}},
                                  X1::Array{Float64,2},
                                  X2::Array{Float64,2}  )
  #
  res[1] = meas[1][idx] - (X2[1,idx] - X1[1,idx])
end

@testset "test specialSampler functionality..." begin

##

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addFactor!(fg, [:x0], SpecialPrior(Normal()))

addVariable!(fg, :x1, ContinuousScalar)
addFactor!(fg, [:x0;:x1], SpecialLinearOffset(Normal(10,1)))

tree, smt, hist = solveTree!(fg)

@test getPPE(fg, :x0).suggested[1] |> abs < 1.0

@test getPPE(fg, :x1).suggested[1] - 10 |> abs < 3.0

## special test for IIF #568

## Singleton (Prior)

fcm, = IIF._getCCW(fg, :x0f1).measurement |> deepcopy
pts = approxConv(fg, :x0f1, :x1)
fcm2, = IIF._getCCW(fg, :x0f1).measurement
fcm3, = IIF._getCCW(fg, :x0f1).measurement |> deepcopy

@test 0.1 < norm(fcm - fcm2)
@test norm(fcm2 - fcm3) < 1e-5

## Pairwise

# forward direction
fcm, = IIF._getCCW(fg, :x0x1f1).measurement |> deepcopy
pts = approxConv(fg, :x0x1f1, :x1)
fcm2, = IIF._getCCW(fg, :x0x1f1).measurement

@test 0.1 < norm(fcm - fcm2)

# reverse direction
fcm, = IIF._getCCW(fg, :x0x1f1).measurement |> deepcopy
pts = approxConv(fg, :x0x1f1, :x0)
fcm2, = IIF._getCCW(fg, :x0x1f1).measurement

@test 0.1 < norm(fcm - fcm2)

0
end

#


text:
function t_operators()

N = 10	# variables
M = 25  # samples

vectordata = rand(M)
matrixdata = rand(N,M)
labels = round.(2*rand(M))
multilabels = round.(2*rand(2,M))

A01 = j4pr.datacell(vectordata)
A02 = j4pr.datacell(vectordata, labels)
A03 = j4pr.datacell(vectordata, multilabels)
A04 = j4pr.datacell(matrixdata)
A05 = j4pr.datacell(matrixdata, labels)
A06 = j4pr.datacell(matrixdata, multilabels)
AV = [A01, A02, A03, A04, A05, A06]

# Data cells
#println("Checking that data cell operators work... ")
#print("\t|> operator...")
for (i,a) in enumerate(AV)
	for (j,b) in enumerate(AV)
		try
			a|>b
			Test.@test true
		catch
			Test.@test false
		end
	end
end
#println("PASSED")

#print("\tvcat operator...")
for (i,a) in enumerate(AV)
	for (j,b) in enumerate(AV)
		try
			[a;b]
			Test.@test true
		catch
			Test.@test false
		end
	end
end
#println("PASSED")

#print("\thcat operator...")
for (i,a) in enumerate(AV)
	for (j,b) in enumerate(AV)
		try
			[a;b]
			Test.@test true
		catch
			Test.@test false
		end
	end
end
#println("PASSED")

# Function cells
W = j4pr.functioncell(x->x)
PT = [W;W]
PPV = j4pr.pipeparallel((W,W), DataStructures.SortedDict(Dict(1=>1, 2=>1)))
PPM = [W W]
PS = W+W

#println("Checking that function cell operators work... ")
#print("\t|> for stacked pipes...")
for (i,a) in enumerate(AV)
	try
		a |> PT
		Test.@test true
	catch
		Test.@test false
	end
end
#println("PASSED")

#print("\t|> for parallel pipes...")
for (i,a) in enumerate(AV)
	try
		if i < 4
			a |> PPV # Send data cells with vector content to their corresponding parallel pipe
		else
			a |> PPM # Send data cells with matrix content to their corresponding parallel pipe
		end
		Test.@test true
	catch
		Test.@test false
	end
end
#println("PASSED")

#print("\t|> for serial pipes...")
for (i,a) in enumerate(AV)
	try
		a |> PS
		Test.@test true
	catch
		Test.@test false
	end
end
#println("PASSED")

end


text:
"""
    heun(f,a,b,y0,M)

berisi program untuk mencari solusi persamaan differensial `y'=f(t,y)`
dengan masalah nilai awal `y(a) = y0` pada interval `[a, b]`. Program ini secara default
berisi 5 masukan, yaitu fungsi `f(t,y)`, titik ujung interval penyelesaian `[a,b]`, nilai awal
`y0`, dan jumlah sub-interval `M`.

# Examples
```jldoctest
julia> heun((t,y)->(t-y)/2,0,3,1,6)
7×2 Array{Float64,2}:
 0.0  1.0
 0.5  0.84375
 1.0  0.831055
 1.5  0.930511
 2.0  1.11759
 2.5  1.37311
 3.0  1.68212
```
return solusi masalah nilai awal `sol`.
"""
function heun(f,a,b,y0,M)
    M = Int(M)
    h = (b-a)/M;
    T = a:h:b;
    Y = Array{Float64}(undef,length(T),1)
    P = Array{Float64}(undef,length(T),1)
    Y[1] = y0;
    # Mulai langkah iterasi Heun
    for k = 1:M
        #% Rumus iterasi Heun
        P[k+1]=Y[k]+h*f(T[k],Y[k]);
        Y[k+1]=Y[k]+h/2*(f(T[k],Y[k])+f(T[k+1],P[k+1]));
    end
    sol = [T Y];
end


text:
"""
Common exceptions
https://github.com/frictionlessdata/datapackage-jl#exceptions
"""

struct PackageError <: Exception
    message::String
end