text:
"""
    function addUpdateMedium_!(p::Program_,platform::String)

Generate the functions related with Medium Update.
"""
function addIntegratorMediumFTCS_!(p::Program_,platform::String)

    if "UpdateMedium" in keys(p.agent.declaredUpdates)

        #Construct functions

            #Make function to compute everything inside the volume

        f = p.agent.declaredUpdates["UpdateMedium"]
        for (i,j) in enumerate(p.agent.declaredSymbols["Medium"]) # Change symbol for update
            f = postwalk(x -> @capture(x,g_(s_)=v_) && s == j && g == DIFFMEDIUMSYMBOL ? :($j = $j + $v*dt) : x, f)
        end
        f = postwalk(x->@capture(x,s_(v_)) ? :($(adaptOperatorsMediumFTCS_(v,s,p))) : x, f) # Adapt
        f = postwalk(x->@capture(x,s_) ? adaptSymbolsMediumFTCS_(s,p) : x, f) # Adapt

        f = vectorizeMedium_(p.agent,f,p)

        f = simpleGridLoopWrapInFunction_(platform,:mediumInnerStep_!, f, p.agent.dims)

        #Remove count over boundaries
        if platform == "cpu"
            f = postwalk(x->@capture(x,1:Nx_) && Nx == :Nx_ ? :(2:Nx_-1) : x, f)
            f = postwalk(x->@capture(x,1:Ny_) && Ny == :Ny_ ? :(2:Ny_-1) : x, f)
            f = postwalk(x->@capture(x,1:Nz_) && Nz == :Nz_ ? :(2:Nz_-1) : x, f)
        elseif platform == "gpu"
            f = postwalk(x->@capture(x,indexX_:strideX_:Nx_) && Nx == :Nx_ ? :(indexX_+1:strideX_:Nx_-1) : x, f)
            f = postwalk(x->@capture(x,indexY_:strideY_:Ny_) && Ny == :Ny_ ? :(indexY_+1:strideY_:Ny_-1) : x, f)
            f = postwalk(x->@capture(x,indexZ_:strideZ_:Nz_) && Nz == :Nz_ ? :(indexZ_+1:strideZ_:Nz_-1) : x, f)
        end
        push!(p.declareF.args, ## Add it to code
                f
            )

        #Add boundary computation
        boundariesFunctionDefinition(p.agent.boundary, p, platform)

        if "UpdateMediumInteraction" in keys(p.agent.declaredUpdates)
            fWrap = wrapInFunction_(:mediumStep_!,
                                    :(begin
                                        @platformAdapt mediumInnerStep_!(ARGS_)
                                        @platformAdapt mediumBoundaryStep_!(ARGS_)
                                        @platformAdapt mediumInteractionStep_!(ARGS_)
                                    end)
                                    )
        else
            fWrap = wrapInFunction_(:mediumStep_!,
                                    :(begin
                                        @platformAdapt mediumInnerStep_!(ARGS_)
                                        @platformAdapt mediumBoundaryStep_!(ARGS_)
                                    end)
                                    )
        end
        push!(p.declareF.args, ## Add it to code
            fWrap
            )
        # println(prettify(code))

        push!(p.execInloop.args,
                :(mediumStep_!(ARGS_))
            )
    end

    return nothing
end

function adaptSymbolsMediumFTCS_(s,p)

    if s == :xₘ

        return :(ic1_*dxₘ_+box[1,1])

    elseif s == :yₘ

        return :(ic2_*dyₘ_+box[2,1])

    elseif s == :zₘ

        return :(ic1_*dzₘ_+box[3,1])

    end

    return s

end

function adaptOperatorsMediumFTCS_(f,op,p)

    f = vectorizeMedium_(p.agent,f,p)

    if op == :δx

        f = :(δMedium_(round(Int,($f-simulationBox[1,1])/dxₘ_)+1,ic1_))

    elseif op == :δy

        if p.agent.dims < 2
            error("Operator δy cannot be used in Agent with dimension 1.")
        end

        f = :(δMedium_(round(Int,($f-simulationBox[2,1])/dyₘ_)+1,ic2_))

    elseif op == :δz

        if p.agent.dims < 3
            error("Operator δz cannot be used in Agent with dimension 1 or 2.")
        end

        f = :(δMedium_(round(Int,($f-simulationBox[3,1])/dzₘ_)+1,ic3_))

    elseif op == :Δ
        if p.agent.dims == 1
            add = postwalk(x->@capture(x,s_[ic1_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,$j]) : x,f)

            f = :(($add+$subs-2*$f)/dxₘ_^2)
        elseif p.agent.dims == 2
            add = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,ic2_,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,ic2_,$j]) : x,f)
            addy = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_+1,$j]) : x,f)
            subsy = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_-1,$j]) : x,f)

            f = :(($add+$subs-2*$f)/dxₘ_^2+($addy+$subsy-2*$f)/dyₘ_^2)
        elseif p.agent.dims == 3
            add = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,ic2_,ic3_,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,ic2_,ic3_,$j]) : x,f)
            addy = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_+1,ic3_,$j]) : x,f)
            subsy = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_-1,ic3_,$j]) : x,f)
            addz = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_,ic3_+1,$j]) : x,f)
            subsz = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_,ic3_-1,$j]) : x,f)

            f = :(($add+$subs-2*$f)/dxₘ_^2+($addy+$subsy-2*$f)/dyₘ_^2+($addz+$subsz-2*$f)/dzₘ_^2)

        end
    elseif op == :Δx
        if p.agent.dims == 1
            add = postwalk(x->@capture(x,s_[ic1_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,$j]) : x,f)

            f = :(($add+$subs-2*$f)/dxₘ_^2)
        elseif p.agent.dims == 2
            add = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,ic2_,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,ic2_,$j]) : x,f)

            f = :(($add+$subs-2*$f)/dxₘ_^2)
        elseif p.agent.dims == 3
            add = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,ic2_,ic3_,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,ic2_,ic3_,$j]) : x,f)

            f = :(($add+$subs-2*$f)/dxₘ_^2)
        end
    elseif op == :Δy
        if p.agent.dims == 1
            error("Operator Δy cannot exist in Agent with dimension 1.")
        elseif p.agent.dims == 2
            addy = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_+1,$j]) : x,f)
            subsy = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_-1,$j]) : x,f)

            f = :(($addy+$subsy-2*$f)/dyₘ_^2)
        elseif p.agent.dims == 3
            addy = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_+1,ic3_,$j]) : x,f)
            subsy = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_-1,ic3_,$j]) : x,f)

            f = :(($addy+$subsy-2*$f)/dyₘ_^2)
        end
    elseif op == :Δz
        if p.agent.dims == 1
            error("Operator Δz cannot exist in Agent with dimension 1.")
        elseif p.agent.dims == 2
            error("Operator Δz cannot exist in Agent with dimension 2.")
        elseif p.agent.dims == 3
            addz = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_,ic3_+1,$j]) : x,f)
            subsz = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_,ic3_-1,$j]) : x,f)

            f = :(($addz+$subsz-2*$f)/dzₘ_^2)
        end
    elseif op == :∇
        if p.agent.dims == 1
            add = postwalk(x->@capture(x,s_[ic1_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,$j]) : x,f)

            f = :(($add-$subs)/(2*dxₘ_))
        elseif p.agent.dims == 2
            add = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,ic2_,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,ic2_,$j]) : x,f)
            addy = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_+1,$j]) : x,f)
            subsy = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_-1,$j]) : x,f)

            f = :(($add-$subs)/(2*dxₘ_)+($addy-$subsy)/(2*dyₘ_))
        elseif p.agent.dims == 3
            add = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,ic2_,ic3_,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,ic2_,ic3_,$j]) : x,f)
            addy = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_+1,ic3_,$j]) : x,f)
            subsy = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_-1,ic3_,$j]) : x,f)
            addz = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_,ic3_+1,$j]) : x,f)
            subsz = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_,ic3_-1,$j]) : x,f)

            f = :(($add-$subs)/(2*dxₘ_)+($addy-$subsy)/(2*dyₘ_)+($addz-$subsz)/(2*dzₘ_))
        end
    elseif op == :∇x
        if p.agent.dims == 1
            add = postwalk(x->@capture(x,s_[ic1_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,$j]) : x,f)

            f = :(($add-$subs)/(2*dxₘ_))
        elseif p.agent.dims == 2
            add = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,ic2_,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,ic2_,$j]) : x,f)

            f = :(($add-$subs)/(2*dxₘ_))
        elseif p.agent.dims == 3
            add = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,ic2_,ic3_,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,ic2_,ic3_,$j]) : x,f)

            f = :(($add-$subs)/(2*dxₘ_))
        end
    elseif op == :∇y
        if p.agent.dims == 1
            error("Operator Δy cannot exist in Agent with dimension 1.")
        elseif p.agent.dims == 2
            addy = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_+1,$j]) : x,f)
            subsy = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_-1,$j]) : x,f)

            f = :(($addy-$subsy)/(2*dyₘ_))
        elseif p.agent.dims == 3
            addy = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_+1,ic3_,$j]) : x,f)
            subsy = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_-1,ic3_,$j]) : x,f)

            f = :(($addy-$subsy)/(2*dyₘ_))
        end
    elseif op == :∇z
        if p.agent.dims == 1
            error("Operator Δz cannot exist in Agent with dimension 1.")
        elseif p.agent.dims == 2
            error("Operator Δy cannot exist in Agent with dimension 2.")
        elseif p.agent.dims == 3
            addz = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_,ic3_+1,$j]) : x,f)
            subsz = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_,ic3_-1,$j]) : x,f)

            f = :(($addz-$subsz)/(2*dzₘ_))
        end
    end

    return f

end


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

mutable struct FlowLogStatusParameters <: SwaggerModel
    targetResourceId::Any # spec type: Union{ Nothing, String } # spec name: targetResourceId

    function FlowLogStatusParameters(;targetResourceId=nothing)
        o = new()
        validate_property(FlowLogStatusParameters, Symbol("targetResourceId"), targetResourceId)
        setfield!(o, Symbol("targetResourceId"), targetResourceId)
        o
    end
end # type FlowLogStatusParameters

const _property_map_FlowLogStatusParameters = Dict{Symbol,Symbol}(Symbol("targetResourceId")=>Symbol("targetResourceId"))
const _property_types_FlowLogStatusParameters = Dict{Symbol,String}(Symbol("targetResourceId")=>"String")
Base.propertynames(::Type{ FlowLogStatusParameters }) = collect(keys(_property_map_FlowLogStatusParameters))
Swagger.property_type(::Type{ FlowLogStatusParameters }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_FlowLogStatusParameters[name]))}
Swagger.field_name(::Type{ FlowLogStatusParameters }, property_name::Symbol) =  _property_map_FlowLogStatusParameters[property_name]

function check_required(o::FlowLogStatusParameters)
    (getproperty(o, Symbol("targetResourceId")) === nothing) && (return false)
    true
end

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

ext: jl


text:
function rouwenhorst(μ,ρ,σ,n)
    n==1 && (return [μ],[1.0])

	mu_eps=0
	q = (ρ+1)/2
	nu = ((n-1)/(1-ρ^2))^(1/2) * σ
	P = [q 1-q; 1-q q]
	for i = 2:n-1
		P = q*[P zeros(i,1);zeros(1,i+1)] + (1-q)*[zeros(i,1) P; zeros(1,i+1)]+ (1-q)*[zeros(1,i+1); P zeros(i,1)] + q*[zeros(1,i+1); zeros(i,1) P]
		P[2:i,:] = P[2:i,:]/2
	end
	x = [linspace(mu_eps/(1-ρ).-nu,mu_eps/(1-ρ).+nu,n);].+μ
	return x,P
end

function cdf_normal(x) :inline
    c = 0.5 * erfc(-x/sqrt(2))
end

function tauchen(μ,ρ,σ,N,m=3)
	if N==1
		return [μ],[1.0]
	end
	Z     = zeros(N)
	Zprob = zeros(N,N)
	a     = (1-ρ)*μ

	Z[N]  = m * sqrt(σ^2 / (1 - ρ^2))
	Z[1]  = -Z[N]
	zstep = (Z[N] - Z[1]) / (N - 1)

	for i=2:(N-1)
	    Z[i] = Z[1] + zstep * (i - 1)
	end

	Z = Z .+ a / (1-ρ)

	for j = 1:N
	    for k = 1:N
	        if k == 1
	            Zprob[j,k] = cdf_normal((Z[1] - a - ρ * Z[j] + zstep / 2) / σ)
	        elseif k == N
	            Zprob[j,k] = 1 - cdf_normal((Z[N] - a - ρ * Z[j] - zstep / 2) / σ)
	        else
	            Zprob[j,k] = cdf_normal((Z[k] - a - ρ * Z[j] + zstep / 2) / σ) - cdf_normal((Z[k] - a - ρ * Z[j] - zstep / 2) / σ)
	        end
	    end
	end
	return Z,Zprob
end


text:
"""
    ContinuousDomain{T <: Real} <: AbstractDomain
An abstract supertype for all continuous domains.
"""
abstract type ContinuousDomain{T<:Real} <: AbstractDomain end

"""
    Intervals{T <: Real} <: ContinuousDomain{T}
An encapsuler to store a vector of `PatternFolds.Interval`. Dynamic changes to `Intervals` are not handled yet.
"""
struct Intervals{T<:Real,L,R} <: ContinuousDomain{T}
    domain::Vector{Interval{T,L,R}}
end

"""
    domain(a::Tuple{T, Bool}, b::Tuple{T, Bool}) where {T <: Real}
    domain(intervals::Vector{Tuple{Tuple{T, Bool},Tuple{T, Bool}}}) where {T <: Real}
Construct a domain of continuous interval(s).
```julia
d1 = domain((0., true), (1., false)) # d1 = [0, 1)
d2 = domain([ # d2 = [0, 1) ∪ (3.5, 42]
    (0., true), (1., false),
    (3.5, false), (42., true),
])
"""
function domain(intervals::Vector{Interval{T,L,R}}) where {T<:Real,L,R}
    return Intervals(map(i -> Interval(i), intervals))
end
domain(interval::Interval) = domain([interval])
domain(intervals::Vector{I}) where {I<:Interval} = Intervals(intervals)

"""
    Base.length(itv::Intervals)
Return the sum of the length of each interval in `itv`.
"""
Base.length(itv::Intervals) = sum(size, get_domain(itv); init=0)

"""
    Base.rand(itv::Intervals)
    Base.rand(itv::Intervals, i)
Return a random value from `itv`, specifically from the `i`th interval if `i` is specified.
"""
Base.rand(itv::Intervals, i::Int) = rand(get_domain(itv)[i])
function Base.rand(itv::Intervals{T}) where {T<:Real}
    r = length(itv) * rand()
    weight = 0.0
    result = zero(T)
    for i in get_domain(itv)
        weight += size(i)
        if weight > r
            result = rand(i)
        end
    end
    return result
end

"""
    Base.in(x, itv::Intervals)
Return `true` if `x ∈ I` for any 'I ∈ itv`, false otherwise. `x ∈ I` is equivalent to
- `a < x < b` if `I = (a, b)`
- `a < x ≤ b` if `I = (a, b]`
- `a ≤ x < b` if `I = [a, b)`
- `a ≤ x ≤ b` if `I = [a, b]`
"""
Base.in(x, itv::Intervals) = any(i -> x ∈ i, get_domain(itv))

"""
    domain_size(itv::Intervals)
Return the difference between the highest and lowest values in `itv`.
"""
function domain_size(itv::Intervals)
    return maximum(last, get_domain(itv)) - minimum(first, get_domain(itv))
end

function merge_domains(d1::D, d2::D) where {D<:ContinuousDomain}
end

function intersect_domains(i₁::I1, i₂::I2) where {I1<:Interval,I2<:Interval}
    a = a_ismore(i₁, i₂) ? i₁.a : i₂.a
    b = b_isless(i₁, i₂) ? i₁.b : i₂.b
    return Interval(a, b)
end

function intersect_domains!(is::IS, i::I, new_itvls) where {IS<:Intervals,I<:Interval}
    for interval in get_domain(is)
        intersection = intersect_domains(interval, i)
        !isempty(intersection) && push!(new_itvls, intersection)
    end
end

function intersect_domains(d₁::D, d₂::D) where {T<:Real,D<:ContinuousDomain{T}}
    new_itvls = Vector{Interval{T}}()
    for i in get_domain(d₂)
        intersect_domains!(d₁, i, new_itvls)
    end
    return Intervals(new_itvls)
end

size(i::I) where {I <: Interval} = span(i)


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020-2021, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

"""
Methods described in:

    Xavier, Álinson S., Feng Qiu, and Shabbir Ahmed. "Learning to solve
    large-scale security-constrained unit commitment problems." INFORMS
    Journal on Computing 33.2 (2021): 739-756. DOI: 10.1287/ijoc.2020.0976
"""
module XavQiuAhm2021

using Distributions
import ..UnitCommitmentInstance

"""
    struct Randomization
        cost = Uniform(0.95, 1.05)
        load_profile_mu = [...]
        load_profile_sigma = [...]
        load_share = Uniform(0.90, 1.10)
        peak_load = Uniform(0.6 * 0.925, 0.6 * 1.075)
        randomize_costs = true
        randomize_load_profile = true
        randomize_load_share = true
    end

Randomization method that changes: (1) production and startup costs, (2)
share of load coming from each bus, (3) peak system load, and (4) temporal
load profile, as follows:

1. **Production and startup costs:**
    For each unit `u`, the vectors `u.min_power_cost` and `u.cost_segments`
    are multiplied by a constant `α[u]` sampled from the provided `cost`
    distribution. If `randomize_costs` is false, skips this step.

2. **Load share:**
   For each bus `b` and time `t`, the value `b.load[t]` is multiplied by
   `(β[b] * b.load[t]) / sum(β[b2] * b2.load[t] for b2 in buses)`, where
   `β[b]` is sampled from the provided `load_share` distribution. If
   `randomize_load_share` is false, skips this step.

3. **Peak system load and temporal load profile:**
    Sets the peak load to `ρ * C`, where `ρ` is sampled from `peak_load` and `C`
    is the maximum system capacity, at any time. Also scales the loads of all
    buses, so that `system_load[t+1]` becomes equal to `system_load[t] * γ[t]`,
    where `γ[t]` is sampled from `Normal(load_profile_mu[t], load_profile_sigma[t])`.

    The system load for the first time period is set so that the peak load
    matches `ρ * C`. If `load_profile_sigma` and `load_profile_mu` have fewer
    elements than `instance.time`, wraps around. If `randomize_load_profile`
    is false, skips this step.

The default parameters were obtained based on an analysis of publicly available
bid and hourly data from PJM, corresponding to the month of January, 2017. For
more details, see Section 4.2 of the paper.
"""
Base.@kwdef struct Randomization
    cost = Uniform(0.95, 1.05)
    load_profile_mu::Vector{Float64} = [
        1.0,
        0.978,
        0.98,
        1.004,
        1.02,
        1.078,
        1.132,
        1.018,
        0.999,
        1.006,
        0.999,
        0.987,
        0.975,
        0.984,
        0.995,
        1.005,
        1.045,
        1.106,
        0.981,
        0.981,
        0.978,
        0.948,
        0.928,
        0.953,
    ]
    load_profile_sigma::Vector{Float64} = [
        0.0,
        0.011,
        0.015,
        0.01,
        0.012,
        0.029,
        0.055,
        0.027,
        0.026,
        0.023,
        0.013,
        0.012,
        0.014,
        0.011,
        0.008,
        0.008,
        0.02,
        0.02,
        0.016,
        0.012,
        0.014,
        0.015,
        0.017,
        0.024,
    ]
    load_share = Uniform(0.90, 1.10)
    peak_load = Uniform(0.6 * 0.925, 0.6 * 1.075)
    randomize_load_profile::Bool = true
    randomize_costs::Bool = true
    randomize_load_share::Bool = true
end

function _randomize_costs(
    instance::UnitCommitmentInstance,
    distribution,
)::Nothing
    for unit in instance.units
        α = rand(distribution)
        unit.min_power_cost *= α
        for k in unit.cost_segments
            k.cost *= α
        end
        for s in unit.startup_categories
            s.cost *= α
        end
    end
    return
end

function _randomize_load_share(
    instance::UnitCommitmentInstance,
    distribution,
)::Nothing
    α = rand(distribution, length(instance.buses))
    for t in 1:instance.time
        total = sum(bus.load[t] for bus in instance.buses)
        den = sum(
            bus.load[t] / total * α[i] for
            (i, bus) in enumerate(instance.buses)
        )
        for (i, bus) in enumerate(instance.buses)
            bus.load[t] *= α[i] / den
        end
    end
    return
end

function _randomize_load_profile(
    instance::UnitCommitmentInstance,
    params::Randomization,
)::Nothing
    # Generate new system load
    system_load = [1.0]
    for t in 2:instance.time
        idx = (t - 1) % length(params.load_profile_mu) + 1
        gamma = rand(
            Normal(params.load_profile_mu[idx], params.load_profile_sigma[idx]),
        )
        push!(system_load, system_load[t-1] * gamma)
    end
    capacity = sum(maximum(u.max_power) for u in instance.units)
    peak_load = rand(params.peak_load) * capacity
    system_load = system_load ./ maximum(system_load) .* peak_load

    # Scale bus loads to match the new system load
    prev_system_load = sum(b.load for b in instance.buses)
    for b in instance.buses
        for t in 1:instance.time
            b.load[t] *= system_load[t] / prev_system_load[t]
        end
    end

    return
end

end

"""
    function randomize!(
        instance::UnitCommitment.UnitCommitmentInstance,
        method::XavQiuAhm2021.Randomization,
    )::Nothing

Randomize costs and loads based on the method described in XavQiuAhm2021.
"""
function randomize!(
    instance::UnitCommitment.UnitCommitmentInstance,
    method::XavQiuAhm2021.Randomization,
)::Nothing
    if method.randomize_costs
        XavQiuAhm2021._randomize_costs(instance, method.cost)
    end
    if method.randomize_load_share
        XavQiuAhm2021._randomize_load_share(instance, method.load_share)
    end
    if method.randomize_load_profile
        XavQiuAhm2021._randomize_load_profile(instance, method)
    end
    return
end

export randomize!


text:
config = MOI.Test.Config(atol=1e-5, rtol=1e-5)
α = 0.7291971148804172
γ = √2 + 2e-1
β = (γ^2 - 2) / γ^2 * (1 + α)
# We have the solution
# P_0 = [α 0
#        0 α]
# P   = [1+α 0
#        0   1+α]
# and the slack of the LMI's
# P - A1'P*A*1 = [β   0
#                 0   1+α]
# P - A2'P*A*2 = [1+α 0
#                 0   β]
optimize!(mock) = MOIU.mock_optimize!(mock, [α, α, 0.0, β, 1 + α, 0.0, 1 + α, β, 0.0])
for mock in mocks(optimize!)
    Tests.quadratic_feasible_lyapunov_switched_system_test(mock, config)
    Tests.quadratic_feasible_scaled_lyapunov_switched_system_test(mock, config)
end
# TODO quadratic infeasible and quartic

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
include("partial_sum.jl")

img = draw_sums(50000, width=1800, height=1800)
imwrite(colorim(img), "partial_sums.png")

max_stars_repo_name: kalvotom/julia-examples

max_issues_repo_name: kalvotom/julia-examples

max_forks_repo_name: kalvotom/julia-examples


text: Dict{String, Any}("plain_table" => Dict{String, Any}("with.dot" => Dict{String, Any}("value" => "4", "type" => "integer"), "plain" => Dict{String, Any}("value" => "3", "type" => "integer")), "with.dot" => Dict{String, Any}("value" => "2", "type" => "integer"), "table" => Dict{String, Any}("withdot" => Dict{String, Any}("key.with.dots" => Dict{String, Any}("value" => "6", "type" => "integer"), "plain" => Dict{String, Any}("value" => "5", "type" => "integer"))), "plain" => Dict{String, Any}("value" => "1", "type" => "integer"))


text:
module sptests

using Base.Test
using TypedPolynomials: StaticPolynomials

macro wrappedallocs(expr)
    @assert expr.head == :call
    f = expr.args[1]
    args = expr.args[2:end]
    argnames = [gensym() for a in args]
    quote
        function g($(argnames...))
            @allocated $(Expr(:call, esc(f), argnames...))
        end
        $(Expr(:call, :g, [esc(a) for a in args]...))
    end
end

@testset "variables" begin
    @polyvar x y

    @test x * x == x^2
    @test x^1 == x
    @test x^0 == 1
    @test (@wrappedallocs x * x) == 0
    @test (@wrappedallocs x^2) == 0
    @test (@wrappedallocs x^1) == 0
#     @test (@wrappedallocs x^0) == 0
end

@testset "orderings" begin
    @polyvar x y z

#     @test Monomial(y) < Monomial(x)
#     @test Monomial(x) > Monomial(y)
#     @test Monomial(y) < x
#     @test y < Monomial(x)
    @test y < x
    @test x > y
#     @test x > Monomial(y)
#     @test Monomial(x) > y
end

@testset "monomials" begin
    @polyvar x y z

    @test typeof(@inferred x * x) == Monomial{(x,), (2,)}
    @test typeof(@inferred x * y) == Monomial{(x, y), (1, 1)}
    @test typeof(@inferred y * x) == Monomial{(x, y), (1, 1)}
    @test typeof(@inferred x * z * y) == Monomial{(x, y, z), (1, 1, 1)}
end

# include("substitution.jl")

end


text:
# Mostly copied over from Cassette in `src/overdub.jl`
# Return `Reflection` for signature `sigtypes` and `world`, if possible. Otherwise, return `nothing`.
function reflect(@nospecialize(sigtypes::Tuple), world::UInt = typemax(UInt))
    if length(sigtypes) > 2 && sigtypes[1] === typeof(invoke)
        @assert sigtypes[3] <: Type{<:Tuple}
        sigtypes = (sigtypes[2], sigtypes[3].parameters[1].parameters...)
    end
    # This works around a subtyping bug. Basically, callers can deconstruct upstream
    # `UnionAll` types in such a way that results in a type with free type variables, in
    # which case subtyping can just break.
    #
    # God help you if you try to use a type parameter here (e.g. `::Type{S} where S<:Tuple`)
    # instead of this nutty workaround, because the compiler can just rewrite `S` into
    # whatever it thinks is "type equal" to the actual provided value. In other words, if
    # `S` is defined as e.g. `f(::Type{S}) where S`, and you call `f(T)`, you should NOT
    # assume that `S === T`. If you did, SHAME ON YOU. It doesn't matter that such an
    # assumption holds true for essentially all other kinds of values. I haven't counted in
    # a while, but I'm pretty sure I have ~40+ hellish years of Julia experience, and this
    # still catches me every time. Who even uses this crazy language?
    S = Tuple{map(s -> Core.Compiler.has_free_typevars(s) ? typeof(s.parameters[1]) : s, sigtypes)...}
    (S.parameters[1]::DataType).name.module === Core.Compiler && return nothing
    _methods = Base._methods_by_ftype(S, -1, world)
    method_index = 0
    for i in 1:length(_methods)
        if _methods[i][1] === S
            method_index = i
            break
        end
    end
    method_index === 0 && return nothing
    type_signature, raw_static_params, method = _methods[method_index]
    method_instance = Core.Compiler.specialize_method(method, type_signature, raw_static_params, false)
    method_signature = method.sig
    static_params = Any[raw_static_params...]
    return method_instance, method_signature, static_params
end

# ugly hack to make differentiating `getproperty` infer a lot better
@generated function _pullback(cx::AContext, ::typeof(literal_getproperty), x, ::Val{f}) where f
    sig(x) = Tuple{x, typeof(f)}
    rrule_sig(x) = Tuple{typeof(getproperty), x, typeof(f)}
    pb_sig(x) = Tuple{cx, typeof(getproperty), x, typeof(f)}

    # either `getproperty` has a custom implementation or `_pullback(cx, getproperty, x, f)`
    # / `rrule(getproperty, x, f) is overloaded directly
    is_getfield_fallback = which(getproperty, sig(x)) == which(getproperty, sig(Any)) &&
        which(_pullback, pb_sig(x)) == which(_pullback, pb_sig(Any)) &&
        which(rrule, rrule_sig(x)) == which(rrule, rrule_sig(Any))

    #ccall(:jl_safe_printf, Cvoid, (Cstring,), "$is_getfield_fallback: $x\n")

    if is_getfield_fallback
        # just copy pullback of `literal_getfield`
        mi, _sig, sparams = reflect((typeof(_pullback), cx, typeof(literal_getfield), x, Val{f}))
        ci = copy(Core.Compiler.retrieve_code_info(mi))

        # we need to change the second arg to `_pullback` from `literal_getproperty` to
        # `literal_getfield`
        Meta.partially_inline!(
            ci.code, Any[_pullback, Core.SlotNumber(2), literal_getfield],
            _sig, sparams, 0, 0, :propagate,
        )
        ci.inlineable = true

        # backedge for `_pullback`, see https://docs.julialang.org/en/v1/devdocs/ast/#MethodInstance
        # this will cause a backedge to this particular MethodInstance to be attached to
        # `_pullback(cx, getproperty, x, f)`
        mi_pb_getproperty, _, _ = reflect((typeof(_pullback), pb_sig(x).parameters...))
        mi_getproperty, _, _ = reflect((typeof(getproperty), sig(x).parameters...))
        mi_rrule, _, _ = reflect((typeof(rrule), rrule_sig(x).parameters...))
        ci.edges = Core.MethodInstance[mi, mi_pb_getproperty, mi_getproperty, mi_rrule]

        return ci
    else
        # nothing to optimize here, need to recurse into `getproperty`
        return quote
            Base.@_inline_meta
            _pullback(cx, getproperty, x, $(QuoteNode(f)))
        end
    end
end


text:
using WGLMakie

x = collect(-4π:0.1:4π)
y = Base.sin.(x)

function plot_d18o(dist, d18o) # Create menu
    fig = Figure()

    colormenu =  Menu(fig, options=[:red, :blue, :black])
    legendmenu = Toggle(fig, active=true)
    fig[1, 1] = vgrid!(
        Label(fig, "Color", width=nothing),
        colormenu,
        Label(fig, "Legend", width=nothing),
        legendmenu,
        width=200, tellheight=false
    )

    ax = Axis(fig[1, 2], size=(800, 400),
            xlabel="Testlabel", ylabel="some quantity (kgm s⁻²)")
    plot = lines!(ax, dist, d18o, c=:black, label="testlegend")

    on(colormenu.selection) do s
        plot.color = s
    end
    # connect!(plot.legend, legendmenu.active)

    fig
end

plot_d18o(range(1,100), cumsum(randn(100)))

include("sinreg.jl")


text:
module Permutations
using LinearAlgebra

import Base: length, show, inv, reverse, ==, getindex, *, ^, sign, hash, getindex,
                Matrix, Array, AbstractMatrix, AbstractArray, Array, adjoint

#                SparseMatrixCSC, AbstractSparseMatrix, AbstractSparseArray, sparse

import Combinatorics: nthperm
import Random: randperm

export Permutation, RandomPermutation
export length, getindex, array, two_row
export inv, cycles, cycle_string
export order, matrix, fixed_points
export longest_increasing, longest_decreasing, reverse, sign
export hash, dict
export CoxeterDecomposition  # no reason to expose CoxeterGenerator

# Defines the Permutation class. Permutations are bijections of 1:n.

abstract type AbstractPermutation end

"""
* `Permutation(list)` creates a new `Permutation`. Here `list` must be a rearrangement of `1:n`.
* `Permutation(n)` creates the identity `Permutation` of `1:n`.
* `Permutation(n,k)` creates the `k`'th `Permutation` of `1:n`.
"""
struct Permutation <: AbstractPermutation
    data::Vector{Int}
    function Permutation(dat::Vector{Int})
        n = length(dat)
        if sort(dat) != collect(1:n)
            error("Improper array: must be a permutation of 1:n")
        end
        new(dat)
    end
end

Permutation(v::AbstractVector) = Permutation(Vector{Int}(v))

# create the k'th permutation of 1:n
function Permutation(n,k)
    return Permutation(nthperm(collect(1:n),k))
end

# create the identity permutation
Permutation(n::Int) = Permutation(collect(1:n))

"""
`RandomPermutation(n)` creates a random permutation of `1:n`,
each with probability `1/factorial(n)`.
"""
RandomPermutation(n::Int) = Permutation(randperm(n))

# Returns the number of elements in the Permtuation

"""
`length(p)` is the number of elements in the `Permutation` `p`.
"""
function length(p::Permutation)
    return length(p.data)
end

# Check for equality of permutations
==(p::Permutation, q::Permutation) = p.data==q.data

# Apply the Permutation to an element: p[k] or p(k)
getindex(p::Permutation, k::Int) = p.data[k]
(p::Permutation)(k::Int) = p.data[k]

# Create a two-row representation of this permutation
"""
`two_row(p)` creates a two-row representation of the `Permutation` `p`
in which the first row is `1:n` and the second row are the
values `p(1), p(2), ..., p(n)`.
"""
function two_row(p::Permutation)
    n = length(p)
    return [ (1:n)'; p.data']
end

# Composition of two permutations
*(p::Permutation) = p
function *(p::Permutation, q::Permutation)
    n = length(p)
    if n != length(q)
        error("Cannot compose Permutations of different lengths")
    end
    @inbounds dat = [ p.data[q.data[k]] for k in 1:n ]
    return Permutation(dat)
end

# Inverse of a permutation
"""
`inv(p)` gives the inverse of `Permutation` `p`. This may
also be computed with `p'`.
"""
function inv(p::Permutation)
    n = length(p)
    data = zeros(Int,n)
    for k=1:n
        @inbounds j = p.data[k]
        @inbounds data[j] = k
    end
    return Permutation(data)
end

adjoint(p::Permutation) = inv(p)

# Find the cycles in a permutation
"""
`cycles(p)` returns a list of the cycles in `Permutation` `p`.
"""
function cycles(p::Permutation)
    n = length(p)
    result = Array{Int,1}[]
    todo = trues(n)
    while any(todo)
        k = findall(todo)[1]
        todo[k] = false
        cycle = [k]
        j = p[k]
        while j != k
            append!(cycle,[j])
            todo[j] = false
            j = p[j]
        end
        push!(result, cycle)
    end
    return result
end

"""
`cycle_string(p)` creates a nice, prinatble string representation
from the cycle structure of the permutation `p`.
"""
function cycle_string(p::Permutation)
    if length(p)==0
        return "()"
    end
    str = ""
    cc = cycles(p)
    for c in cc
        str *= array2string(c)
    end
    return str
end

# helper function for cycle_string (not exposed). Converts an integer
# array such as [1,3,5,2] into the string "(1,3,5,2)".
function array2string(a::Array{Int,1})
    n = length(a)
    res = "("
    for k=1:n
        res *= string(a[k])
        if k<n
            res *= ","
        else
            res *= ")"
        end
    end
    return res
end

"""
`sign(p)` is `+1` is `p` is an even `Permtuation` and `-1` if p is odd.
"""
function sign(p::Permutation)
    cc = cycles(p)
    szs = [ length(c)+1 for c in cc ]
    par = sum(szs)%2
    return (par>0) ? -1 : 1
end

# When we print a permutation, it appears in disjoint cycle format.
function show(io::IO, p::Permutation)
    print(io, cycle_string(p))
end

# Find the smallest positive n such that p^n is the identity
"""
`order(p)` is the smallest positive integer `n`
such that `p^n` is the identity `Permutation`.
"""
function order(p::Permutation)
    result = 1
    clist = cycles(p)
    for c in clist
        result = lcm(result, length(c))
    end
    return result
end

# Extend p^n so negative exponents work too
function ^(p::Permutation, n::Int)
    if n==0
        return Permutation(length(p))
    end
    if n<0
        return inv(p)^(-n)
    end
    if n==1
        return p
    end

    m::Int = round(Int,floor(n/2)) # int(floor(n/2))  # m = floor(n/2)
    q::Permutation = p^m

    if n%2 == 0   # if even
        return q*q
    end

    return p*q*q
end

# Represent as a permtuation matrix.
"""
`Matrix(p)` returns the permutation matrix for the `Permutation` `p`.
"""
function Matrix{T}(p::Permutation) where T
    n = length(p)
    A = Matrix{T}(I,n,n)     # A = eye(T,n)   #  int(eye(n))
    return A[:,p.data]
end

Matrix(p::Permutation) = Matrix{Int}(p)
Array(p::Permutation) = Matrix(p)
AbstractMatrix(p::Permutation) = Matrix(p)
AbstractArray(p::Permutation) = Matrix(p)

Array{T}(p::Permutation) where T = Matrix{T}(p)
AbstractMatrix{T}(p::Permutation) where T = Matrix{T}(p)
AbstractArray{T}(p::Permutation) where T = Matrix{T}(p)

#
# """
# `SparseMatrixCSC(p)` returns the permutation matrix for the `Permutation` `p`.
# """
# function SparseMatrixCSC{T}(p::Permutation) where T
#     n = length(p)
#     A = speye(T,n)   #  int(eye(n))
#     return A[p.data,:]
# end
#
# SparseMatrixCSC(p::Permutation) = SparseMatrixCSC{Int}(p)
# AbstractSparseMatrix(p::Permutation) = SparseMatrixCSC{Int}(p)
# AbstractSparseArray(p::Permutation) = SparseMatrixCSC{Int}(p)
# AbstractSparseMatrix{T}(p::Permutation) where T = SparseMatrixCSC{T}(p)
# AbstractSparseArray{T}(p::Permutation) where T = SparseMatrixCSC{T}(p)
#
# sparse(p::Permutation) = SparseMatrixCSC(p)

# find the fixed points of a Permutation
"""
`fixed_points(p)` returns the list of values `k` for which `p(k)==k`.
"""
fixed_points(p::Permutation) = findall([ p[k]==k for k in 1:length(p)])

# Find a longest monotone subsequence of p in a given direction. This
# is not exposed, but it is used by longest_increasing and
# longest_increasing.
function longest_monotone(p::Permutation, order=<)
    n = length(p)
    if n==0
        return Int[]
    end
    scores = ones(Int,n)
    scores[n] = 1

    for k=n-1:-1:1
        for i=k+1:n
            if order(p[k],p[i]) && scores[k] <= scores[i]
                scores[k] = scores[i]+1
            end
        end
    end

    seq = Int[]
    mx = maximum(scores)
    a = findfirst([s == mx for s in scores])
    push!(seq,a)

    while scores[a] > 1
        val = scores[a]-1
        idx = findfirst([scores[k]==val && order(p[a],p[k]) for k=a+1:n ])
        a += idx
        push!(seq,a)
    end

    return [p[i] for i in seq]

end

"""
Let `p` be a `Permutation`. Thinking of this as a list of values,
`longest_increasing(p)` returns a longest subsequence of values
that are in ascending order
"""
longest_increasing(p::Permutation) = longest_monotone(p,<)

"""
Let `p` be a `Permutation`. Thinking of this as a list of values,
`longest_decreasing(p)` returns a longest subsequence of values
that are in descending order
"""
longest_decreasing(p::Permutation) = longest_monotone(p,>)

"""
Let `p` be a `Permutation`. Thinking of `p` as a list of values,
`reverse(p)` creates a `Permutation` with those values in reverse
sequence.
"""
function reverse(p::Permutation)
    d = reverse(p.data)
    return Permutation(d)
end

# hash function so Permutations can be keys in dictionaries, etc.
hash(p::Permutation, h::UInt64) = hash(p.data,h)

"""
`dict(p)` converts a permutation into a dictionary. If `d` is the
result then `d[k]` equals `p(k)` for all `k` in the domain of the
permutation.
"""
function dict(p::Permutation)
    n = length(p)
    d = Dict{Int,Int}()
    for i=1:n
        d[i] = p(i)
    end
    return d
end

#####
# Decomposing into Coxeter generators
#####

struct CoxeterGenerator <: AbstractPermutation
    n::Int
    i::Int
    function CoxeterGenerator(n::Int, i::Int)
        1 ≤ i ≤ n-1 || throw(ArgumentError("$i must be between 1 and $n-1"))
        new(n, i)
    end
end

length(sᵢ::CoxeterGenerator) = sᵢ.n
Permutation(P::CoxeterGenerator) = Permutation([1:P.i-1; P.i+1; P.i; P.i+2:P.n])

function show(io::IO, p::CoxeterGenerator)
    print(io, "length $(p.n) permutation: s_$(p.i)")
end

# _coxeter_reduce reduces a product of simple transpositions using the relationships
# https://en.wikipedia.org/wiki/Symmetric_group#Generators_and_relations
# combined with a sorting for uniqueness
function _coxeter_reduce!(terms::Vector{Int})
    for i = 1:length(terms)-1
        # sᵢ^2 = I
        if terms[i] == terms[i+1]
            return _coxeter_reduce!(deleteat!(terms, i:i+1))
        end
        # sort using s_is_j = s_js_i
        if terms[i+1] ≠ terms[i]-1 && terms[i+1] ≠ terms[i]+1 && terms[i] > terms[i+1]
            terms[i], terms[i+1] = terms[i+1], terms[i]
            return _coxeter_reduce!(terms)
        end
    end
    # (s_is_{i+1})^3 = I
    for i = 1:length(terms)-5
        if terms[i]+1 == terms[i+1] == terms[i+2]+1 == terms[i+3] == terms[i+4]+1 == terms[i+5]
            return _coxeter_reduce!(deleteat!(terms, i:i+5))
        end
    end

    terms
end

"""
`CoxeterDecomposition(p)` expresses the `Permutation` `p` as a
composition of transpositions of the form `(k,k+1)`.
"""
struct CoxeterDecomposition <: AbstractPermutation
    n::Int
    terms::Vector{Int}
    function CoxeterDecomposition(n::Int, terms::Vector{Int})
        for t in terms
            1 ≤ t ≤ n-1 || throw(ArgumentError("$t must be between 1 and $n-1"))
        end
        new(n, _coxeter_reduce!(terms))
    end
end

CoxeterDecomposition(n::Int, terms::AbstractVector) = CoxeterDecomposition(n, Vector{Int}(terms))
CoxeterDecomposition(sᵢ::CoxeterGenerator) = CoxeterDecomposition(sᵢ.n, [sᵢ.i])

length(P::CoxeterDecomposition) = P.n

function Permutation(P::CoxeterDecomposition)
    if isempty(P.terms)
        Permutation(1:P.n)
    else
        *(Permutation.(CoxeterGenerator.(P.n,P.terms))...)
    end
end
CoxeterDecomposition(P::Permutation) = CoxeterDecomposition!(Permutation(copy(P.data)))
CoxeterDecomposition!(P::Permutation) = CoxeterDecomposition(length(P), _coxeterdecomposition!(P))
function _coxeterdecomposition!(P::Permutation)
    n = length(P)
    data = P.data
    ret = Int[]
    while !issorted(data)
        for k=1:n-1
            if data[k] > data[k+1]
                data[k], data[k+1] =  data[k+1], data[k]
                push!(ret, k)
            end
        end
    end
    reverse!(ret)
end

==(A::CoxeterDecomposition, B::CoxeterDecomposition) = A.terms == B.terms

function *(A::CoxeterGenerator, B::CoxeterGenerator)
    length(A) == length(B) || throw(ArgumentError("Permutations must have same length to multiply"))
    CoxeterDecomposition(length(A), [A.i,B.i])
end
function *(A::CoxeterGenerator, B::CoxeterDecomposition)
    length(A) == length(B) || throw(ArgumentError("Permutations must have same length to multiply"))
    CoxeterDecomposition(length(A), [A.i; B.terms])
end
function *(A::CoxeterDecomposition, B::CoxeterGenerator)
    length(A) == length(B) || throw(ArgumentError("Permutations must have same length to multiply"))
    CoxeterDecomposition(length(A), [A.terms; B.i])
end
function *(A::CoxeterDecomposition, B::CoxeterDecomposition)
    length(A) == length(B) || throw(ArgumentError("Permutations must have same length to multiply"))
    CoxeterDecomposition(length(A), [A.terms; B.terms])
end

inv(A::CoxeterDecomposition) = CoxeterDecomposition(A.n, reverse(A.terms))
adjoint(A::CoxeterDecomposition) = inv(A)

function show(io::IO, p::CoxeterDecomposition)
    print(io, "Permutation of 1:$(length(p)): ")
    if length(p.terms)==0
        print(io,"()")
    end
    for i in p.terms
        # print(io, "s_$(i)")
        print(io,"($i,$(i+1))")
    end
end

# deprecations
@deprecate array(p::Permutation) p.data
@deprecate matrix(p::Permutation, sparse::Bool = false) sparse ? sparse(p) : Matrix(p)

include("sqrt.jl")

end # end of module Permutations


text:
function WeightingOp(in,adj;w=1)

	return in.*w

end

function WeightingOp(m::AbstractString,d::AbstractString,adj;w="NULL")

	if (adj==true)
		d1,h1,e1 = SeisRead(d)
		d2,h2,e2 = SeisRead(w)
		SeisWrite(m,d1[:,:].*d2[:,:],h1,e1)
	else
		d1,h1,e1 = SeisRead(m)
		d2,h2,e2 = SeisRead(w)
		SeisWrite(d,d1[:,:].*d2[:,:],h1,e1)
	end

end

function WeightingOp(m::Array{AbstractString,1},d::Array{AbstractString,1},adj;w="NULL")

	if !(isdefined(w,2))
		for j = 1 : length(m)
			WeightingOp(m[j],d[j],adj,w=w)
		end
	else
		for j = 1 : length(m)
			WeightingOp(m[j],d[j],adj,w=w[j])
		end
	end
end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
using Symbolics
using Test
using Latexify

@variables x y z

@test latexify(x^-1) == Latexify.L"$x^{-1}$"

@test latexify((z + x*y^-1) / sin(z)) ==
    Latexify.L"$\left( z + x \cdot y^{-1} \right) \cdot \sin^{-1}\left( z \right)$"

@test latexify((3x - 7y*z^23) * (z - z^2) / x) ==
    Latexify.L"x^{-1} \cdot \left( z -1 \cdot z^{2} \right) \cdot \left( 3 \cdot x -7 \cdot y \cdot z^{23} \right)"


text:
@testset "environments" begin

    function general_environment_test(env, steps=10^2)
        reset!(env)
        @test observe(env).isdone == false
        while steps > 0
            if observe(env).isdone
                reset!(env)
                @test observe(env).isdone == false
            end
            env(sample(actionspace(env)))
            @test observe(env).observation ∈ observationspace(env)
            steps -= 1
        end
    end

    @testset "SimpleRandomWalkEnv" begin
        env = SimpleRandomWalkEnv(N=7, actions=[-1, 1], start=3)
        general_environment_test(env)
    end

    @testset "CartPoleEnv" begin
        env = CartPoleEnv()
        general_environment_test(env)
    end
end


text:
@testset "540.single-element-in-a-sorted-array.jl" begin
    @test single_non_duplicate([1, 1, 2, 3, 3, 4, 4, 8, 8]) == 2
    @test single_non_duplicate([3, 3, 7, 7, 10, 11, 11]) == 10
    @test single_non_duplicate([3, 3, 7, 7, 10, 11, 11, 15, 15, 18, 18]) == 10
    @test single_non_duplicate([2, 3, 3, 4, 4, 8, 8]) == 2
    @test single_non_duplicate([1, 1, 2]) == 2
    @test single_non_duplicate([2]) == 2
end


text:
using Test, OTSIN, LinearAlgebra

include("optimaltransport.jl")
include("fitting.jl")

max_stars_repo_name: PoisotLab/OTSIN

max_issues_repo_name: PoisotLab/OTSIN

max_forks_repo_name: PoisotLab/OTSIN


text:
#=
function Base.show(io::IO, op_model::OperationModel)
    println(io, "Operation Model")
end
=#

function _organize_device_model(val::Dict{Symbol,DeviceModel}, field::Symbol, io::IO)

    println(io, "  $(field): ")
    for (i, ix) in val

        println(io, "      $(i):")
        for inner_field in fieldnames(DeviceModel)

            value = getfield(val[i], Symbol(inner_field))

            if !isnothing(value)
                println(io, "        $(inner_field) = $value")
            end
        end
    end

end

"""
    Base.show(io::IO, ::MIME"text/plain", op_model::OperationModel)

This function goes through the fields in OperationModel and then in ModelReference,
if the field contains a Device model dictionary, it calls organize_device_model() &
prints the data by field, key, value. If the field is not a Device model dictionary,
and a value exists for that field it prints the value.

"""
function Base.show(io::IO, ::MIME"text/plain", op_model::OperationModel)

    println(io, "\nOperation Model")
    println(io, "===============\n")

    for field in fieldnames(ModelReference)

        val = getfield(op_model.model_ref, Symbol(field))

        if typeof(val) == Dict{Symbol,DeviceModel}

            _organize_device_model(val, field, io)

        else
            if !isnothing(val)
                println(io, "  $(field):  $(val)")
            else
                println(io, "no data")
            end
        end
    end
end

function Base.show(io::IO, op_model::CanonicalModel)
    println(io, "Canonical()")
end

function Base.show(io::IO, op_model::Simulation)
    println(io, "Simulation()")
end
#=
function Base.show(io::IO, res_model::OperationModelResults)
    println(io, "Results Model")
 end

=#

function Base.show(io::IO, res_model::OperationModelResults)
    println(io, "\nResults Model")
    println(io, "===============\n")

    for (k, v) in res_model.variables
        time = DataFrames.DataFrame(Time = res_model.time_stamp[!, :Range])
        var = hcat(time, v)
        println(io, "$(k)")
        println(io, "==================")
        println(io, "$(var)\n")
    end

    println(io, "Optimizer Log")
    println(io, "-------------")
    for (k, v) in res_model.optimizer_log
        println(io, "        $(k) = $(v)")
    end
    println(io, "\n")
    for (k, v) in res_model.total_cost
        println(io, "Total Cost: $(k) = $(v)")
    end
 end

 function Base.show(io::IO, stage::Stage)
    println(io, "Stage()")
 end


text:
using HetaSimulator, Plots

p = load_platform("$HetaSimulatorDir/cases/story_3")

# p1 = load_jlplatform("./cases/story_3/dist/julia_platform/model.jl")

### load scenarios

scn_df = read_scenarios("$HetaSimulatorDir/cases/story_3/scenarios.csv")
# scn_df = read_scenarios("$HetaSimulatorDir/cases/story_3/scenarios.xlsx")
add_scenarios!(p, scn_df)

scenario1 = scenarios(p)[:dose_1]

### create scenarios

model = models(p)[:nameless]
new_scenario = Scenario(
    model,
    parameters = [:dose=>100.],
    events_active = [:sw1=>false, :sw1=>true],
    tspan = (0.,1000.),
    observables = [:A0, :C1, :C2, :v_abs, :v_el, :v_distr]
    )
push!(scenarios(p), :multiple_100=>new_scenario)

### create measurements

measurements_df = read_measurements("$HetaSimulatorDir/cases/story_3/measurements.csv")
add_measurements!(p, measurements_df)

### Simulation

res = sim(p)
# plotd = plot(res)
# savefig(plotd, "sim1.png")
# plotd = plot(res[1])
# savefig(plotd, "sim2.png")

res |> plot
res_df = DataFrame(res) # XXX: error here

plot(res[1])
res_df1 = DataFrame(res[1])

### Monte-Carlo

mc_res = mc(p, [:kabs=>Normal(10.,5e-1), :kel=>Normal(0.2,5e-3)], 1000)
mc_res |> DataFrame # error here
mc_res |> plot
# plotd = plot(mc_res)
# savefig(plotd, "sim3.png")

# DataFrame(mc_res[:dose_1])
# plot(mc_res[:dose_1])

### Fitting

res0 = sim(p)
plot(res0, yscale=:log, ylims=(1e-3,1e2))
# plotd = plot(res0, yscale=:log, ylims=(1e-3,1e2))
# savefig(plotd, "sim4.png")

to_fit = [
    :kabs => 8.0,
    :Q => 4.0,
    :kel => 2.2,
    :sigma1 => 0.1,
    :sigma2 => 0.1,
    :sigma3 => 0.1,
]
fit_res = fit(p, to_fit)

optim(fit_res)

res_optim = sim(p, parameters_upd = optim(fit_res))
plot(res_optim, yscale=:log, ylims=(1e-3,1e2))
# plotd = plot(res_optim, yscale=:log, ylims=(1e-3,1e2))
# savefig(plotd, "sim5.png")

### to check
res_x = sim(p, parameters_upd = [:kabs => 20.0, :Q => 1.0, :kel => 5.0e-01])
plot(res_x, yscale=:log, ylims=(1e-3,1e2))

### fitting with parameters.csv table

# all parameters are estimated (the same as fit_res)
params_df0 = read_parameters("$HetaSimulatorDir/cases/story_3/parameters0.csv")
fit_res0 = fit(p, params_df0)

# overwrite dose and estimate other parameters
params_df1 = read_parameters("$HetaSimulatorDir/cases/story_3/parameters.csv")
fit_res1 = fit(p, params_df1)


text:
include("Vtkio.jl")
@reexport using Gridap.Vtkio


text:
@testset "Dijkstra" begin
    g4 = path_digraph(5)
    d1 = float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0])
    d2 = sparse(float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0]))

    for g in testdigraphs(g4)
        y = @inferred(dijkstra_shortest_paths(g, 2, d1))
        z = @inferred(dijkstra_shortest_paths(g, 2, d2))

        @test y.parents == z.parents == [0, 0, 2, 3, 4]
        @test y.dists == z.dists == [Inf, 0, 6, 17, 33]

        y = @inferred(dijkstra_shortest_paths(g, 2, d1; allpaths=true))
        z = @inferred(dijkstra_shortest_paths(g, 2, d2; allpaths=true))
        @test z.predecessors[3] == y.predecessors[3] == [2]

        @test @inferred(enumerate_paths(z)) == enumerate_paths(y)
        @test @inferred(enumerate_paths(z))[4] ==
          enumerate_paths(z, 4) ==
          enumerate_paths(y, 4) == [2, 3, 4]
    end

    gx = path_graph(5)
    add_edge!(gx, 2, 4)
    d = ones(Int, 5, 5)
    d[2, 3] = 100
    for g in testgraphs(gx)
        z = @inferred(dijkstra_shortest_paths(g, 1, d))
        @test z.dists == [0, 1, 3, 2, 3]
        @test z.parents == [0, 1, 4, 2, 4]
    end

    # small function to reconstruct the shortest path; I copied it from somewhere, can't find the original source to give the credits
    # @Beatzekatze on github
    spath(target, dijkstraStruct, source) = target == source ? target : [spath(dijkstraStruct.parents[target], dijkstraStruct, source) target]
    spaths(ds, targets, source) = [spath(i, ds, source) for i in targets]

    G = LightGraphs.Graph()
    add_vertices!(G, 4)
    add_edge!(G, 2, 1)
    add_edge!(G, 2, 3)
    add_edge!(G, 1, 4)
    add_edge!(G, 3, 4)
    add_edge!(G, 2, 2)
    w = [0. 3. 0. 1.;
        3. 0. 2. 0.;
        0. 2. 0. 3.;
        1. 0. 3. 0.]

    for g in testgraphs(G)
        ds = @inferred(dijkstra_shortest_paths(g, 2, w))
      # this loop reconstructs the shortest path for vertices 1, 3 and 4
        @test spaths(ds, [1, 3, 4], 2) == Array[[2 1],
                                            [2 3],
                                            [2 1 4]]

    # here a selflink at source is introduced; it should not change the shortest paths
        w[2, 2] = 10.0
        ds = @inferred(dijkstra_shortest_paths(g, 2, w))
        shortest_paths = []
      # this loop reconstructs the shortest path for vertices 1, 3 and 4
        @test spaths(ds, [1, 3, 4], 2) == Array[[2 1],
                                            [2 3],
                                            [2 1 4]]
    end

    #615
    m = [0 2 2 0 0; 2 0 0 0 3; 2 0 0 1 2;0 0 1 0 1;0 3 2 1 0]
    G = SimpleGraph(5)
    add_edge!(G, 1, 2)
    add_edge!(G, 1, 3)
    add_edge!(G, 2, 5)
    add_edge!(G, 3, 5)
    add_edge!(G, 3, 4)
    add_edge!(G, 4, 5)
    for g in testgraphs(G)
        ds = @inferred(dijkstra_shortest_paths(g, 1, m; allpaths=true))
        @test ds.pathcounts   == [1, 1, 1, 1, 2]
        @test ds.predecessors == [[], [1], [1], [3], [3, 4]]
        @test ds.predecessors == [[], [1], [1], [3], [3, 4]]

        dm = @inferred(dijkstra_shortest_paths(g, 1; allpaths=true, trackvertices=true))
        @test dm.pathcounts       == [1, 1, 1, 1, 2]
        @test dm.predecessors     == [[], [1], [1], [3], [2, 3]]
        @test dm.closest_vertices == [1, 2, 3, 5, 4]
    end

    G = SimpleGraph(5)
    add_edge!(G, 1, 2)
    add_edge!(G, 1, 3)
    add_edge!(G, 4, 5)
    for g in testgraphs(G)
        dm = @inferred(dijkstra_shortest_paths(g, 1; allpaths=true, trackvertices=true))
        @test dm.closest_vertices == [1, 2, 3, 4, 5]
    end
end


text:
# fbp/z-test.jl

using Test: @testset

list = [
"image_geom.jl"
"sino_geom.jl"
#
"cuboid_im.jl"
"disk-phantom.jl"
"ellipse_im.jl"
"ellipsoid_im.jl"
"ellipse_sino.jl"
"rect_im.jl"
"rect_sino.jl"
"rotate2d.jl"
]

for file in list
	@testset "$file" begin
		include(file)
	end
end


text:
# printing for static taylor
function coeffstring(t::STaylor1, i, variable=:t)
    if i == 1  # order 0
        return string(t.coeffs[i])
    end
    if i == 2  # order 1
        return string(t.coeffs[i], variable)
    end
    return string(t.coeffs[i], variable, "^", i-1)
end

function print_taylor(io::IO, t::STaylor1, variable=:t)
    print(io, "(" * join([coeffstring(t, i, variable) for i in 1:length(t.coeffs)], " + ") * ")")
end

Base.show(io::IO, t::STaylor1) = print_taylor(io, t)
function Base.show(io::IO, t::STaylor1{N,T}) where {N, T<:STaylor1}
    print_taylor(io, t, :t)
end


text:
# using these is a little slower than the way we do it in the code, but ok
# can say:
# for i in 1:deg(mat, v)
#   nbr = nbri(mat, v, i)
#   wt = weighti(mat, v, i)
#
# but, it is faster to write:
#  for ind in colptr[curNode]:(colptr[curNode+1]-1)
#    nbr = rowval[ind]

# TODO: remove late
#=

deg(mat::SparseMatrixCSC{Tv,Ti}, v::Ti) where {Tv,Ti} = mat.colptr[v+1]-mat.colptr[v]
nbri(mat::SparseMatrixCSC{Tv,Ti}, v::Ti, i::Ti) where {Tv,Ti} = mat.rowval[mat.colptr[v]+i-1]
weighti(mat::SparseMatrixCSC{Tv,Ti}, v::Ti, i::Ti) where {Tv,Ti} = mat.nzval[mat.colptr[v]+i-1]
nbrs(mat::SparseMatrixCSC{Tv,Ti}, v::Ti) where {Tv,Ti} = mat.rowval[mat.colptr[v]:(mat.colptr[v+1]-1)]

" Finds the weighted degree of a vertex in the graph "
function wdeg(mat::SparseMatrixCSC{Tv,Ti}, v::Ti) where {Tv,Ti}
  sum = 0
  for i in 1:deg(mat,v)
    sum = sum + weighti(mat,v,i)
  end
  return sum
end

" Sets the value of a certain edge in a sparse graph; value can be 0 without the edges dissapearing "
function setValue(mat::SparseMatrixCSC{Tv,Ti}, v::Ti, i::Ti, a::Tv) where {Tv,Ti}
  mat.nzval[mat.colptr[v]+i-1] = a
end # setValue

=#
""" For a symmetric matrix, this gives the correspondance between pairs of entries in an ijv.
So, ai[ind] = aj[flip[ind]].  For example,

~~~
(ai,aj,av) = findnz(a);
fl = flipIndex(a)
ind = 10
@show backind = fl[10]
@show [ai[ind], aj[ind], av[ind]]
@show [ai[backind], aj[backind], av[backind]];

backind = fl[10] = 4
[ai[ind],aj[ind],av[ind]] = [2.0,4.0,0.7]
[ai[backind],aj[backind],av[backind]] = [4.0,2.0,0.7]
~~~
"""
function flipIndex(a::SparseMatrixCSC{Tval,Tind}) where {Tval,Tind}

  b = SparseMatrixCSC(a.m, a.n, copy(a.colptr), copy(a.rowval), collect(UnitRange{Tind}(1,nnz(a))) );
  bakMat = copy(b');
  return bakMat.nzval

end

# TODO: remove later
#=
" Computes the back indices in a graph in O(M+N). works if for every edge (u,v), (v,u) is also in the graph "
function backIndices(G::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
  n = max(G.n, G.m)

  # initialize values
  u,v,w = findEntries(G)

  backInfo = sparse(u, v, [(0,0) for i in 1:length(u)])
  backIndices = sparse(u, v, [-1 for i in 1:length(u)])
  index = zeros(Ti, n)

  # compute for each vertex the incoming edges
  for u in 1:n
    for i in 1:deg(G, u)
      v = nbri(G, u, i)
      index[v] = index[v] + 1
      setValue(backInfo, v, index[v], (u, i))
    end
  end

  # create the back edges
  incoming = zeros(Ti, n)
  for u in 1:n
    for i in 1:deg(G, u)
      incoming[weighti(backInfo, u, i)[1]] = weighti(backInfo, u, i)[2]
    end

    for i in 1:deg(G, u)
      v = nbri(G, u, i)
      setValue(backIndices, v, incoming[v], i)
    end
  end

  return backIndices

end # backIndices

" Same as the above, but now the graph is in adjacency list form "
function backIndices(G::Array{Array{Tuple{Tv1,Tv2},1},1}) where {Tv1,Tv2}
  n = length(G)

  # initialize values
  backInfo = [Tuple{Int64,Int64}[] for i in 1:n]
  backIndices = [(Int64)[] for i in 1:n]
  for u in 1:length(G)
    listU = G[u]
    backInfo[u] = [(0,0) for i in 1:length(listU)]
    backIndices[u] = [-1 for i in 1:length(listU)]
  end
  index = zeros(Int64, n)

  # compute for each vertex the incoming edges
  for u in 1:n
    for i in 1:length(G[u])
      v,w = G[u][i]
      index[v] = index[v] + 1
      backInfo[v][index[v]] = (u, i)
    end
  end

  # create the back edges
  incoming = zeros(Int64, n)
  for u in 1:n
    for i in 1:length(G[u])
      incoming[backInfo[u][i][1]] = backInfo[u][i][2]
    end

    for i in 1:length(G[u])
      v,w = G[u][i]
      backIndices[v][incoming[v]] = i
    end
  end

  return backIndices

end # backIndices

" Similar to findnz, but also returns 0 entries that have an edge in the sparse matrix "
function findEntries(G::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}

  n = max(G.n, G.m)

  u = copy(G.rowval)
  v = zeros(Ti, length(u))
  for i in 1:n
    for j in (G.colptr[i]):(G.colptr[i + 1] - 1)
      v[j] = i
    end
  end
  w = copy(G.nzval)

  return u,v,w

end # findEntries

"""
  Returns the quality of the cut for a given graph and a given cut set s.
  the result will be |outgoing edges| / min(|vertices in set|, |N - vertices in set|)
"""
function compConductance(G::SparseMatrixCSC{Tv,Ti}, s::Array{Int64,1}) where {Tv, Ti}

  n = max(G.n, G.m)

  vols = getVolume(G, s)
  obound = getObound(G, s)

  return obound / vols

end # compConductance

" Computes the volume of subset s in an unweighted graph G "
function getVolume(G::SparseMatrixCSC{Tv,Ti}, s::Array{Int64,1}) where {Tv,Ti}

  vol = 0.0
  for i in 1:length(s)
    deggsi = 0
    for j in 1:deg(G, s[i])
      deggsi = deggsi + weighti(G,s[i],j)
    end
    # vol = vol + deg(G, s[i])
    vol = vol + deggsi
  end

  return vol

end # getVolume

" Computes the number of edges leaving s "
function getObound(G::SparseMatrixCSC{Tv,Ti}, s::Array{Int64,1}) where {Tv,Ti}

  sets = BitSet(s)

  obound = 0
  for u in s
    for i in 1:deg(G, u)
      v = nbri(G, u, i)
      if !(v in sets)
        obound = obound + weighti(G, u, i)
      end
    end
  end

  return obound

end #getObound

=#


text:
using DataDeps

register(DataDep(
  "forest-fires",
  "https://archive.ics.uci.edu/ml/datasets/Forest+Fires",
  "https://archive.ics.uci.edu/ml/machine-learning-databases/forest-fires/forestfires.csv",
  "0d6586a1fa52f55bef48578aef14eb97273f1e9330e1a53423df497a77065253",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=13,
        feature_indices=1:12,
        categoric_indices=3:4,
        skipto=2,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
"Represents an infinite vector, accessible via getindex, setindex"
abstract type AbstractInfiniteVector{T} end


text:
module Parse
import DataStructures: OrderedDict
using ..Commands

export well_formed, matched_rule, parse_command

REGS = (
    text_object = r"^.?([ai][wWsp])$",
    # complete =
)
"""
Grammar for Vim's little language:
    1. motion
        w
        E
    2. operator [motion|textobject]
        dw
        cW
        yE
        ciw
    3. operator operator
        yy
        dd
        3dd
"""
repeat = "\\d*"
motion="[\$%^\\(\\)wWeE{}hjklGHLbB]"
# r_motion=r"([$%^\(\)wWeE{}hjklGHLbB]"
textobject="$repeat[ai][wWsp]"
operator="[ydc]"
rules = (
    # insert commands
    r"^(?<c>[aAiIoO])$",
    "^(?<n1>$repeat)(?<motion>$motion)\$" |> Regex,
    "^(?<n1>$repeat)(?<op>$operator)(?<n2>$repeat)(?:(?<motion>$motion)|(?<to>$textobject))\$" |> Regex,
    "^(?<n1>$repeat)(?<op>$operator)(\\k<op>)\$" |> Regex
)

"""
Determines whether the given string is accepted and successfully terminates
"""
function well_formed(cmd :: String) :: Bool
    for rule in rules
        if occursin(rule, cmd)
            return true
        end
    end
    return false
end

function matched_rule(cmd :: String)
    for rule in rules
        if occursin(rule, cmd)
            return rule
        end
    end
    return nothing
end

"""
    Get the typed value of `item`
"""
function parse_value(item :: Union{Nothing, AbstractString}) :: Union{Integer, Char, String, Nothing}
    if item === nothing
        return nothing
    end
    val = Meta.parse(item)
    if typeof(val) <: Integer
        return val
    elseif typeof(val) <: Symbol
        s = String(val)
        if length(s) == 1
            return s[1]
        else
            return s
        end
    end
end

"""
    Attempt to parse a command, return nothing if `s` could not be parsed into a command
"""
function parse_command(s :: AbstractString) :: Union{Command, Nothing}
    if !well_formed(s)
        @info "command not well formed", s
        return nothing
    end
    r = matched_rule(s)
    return command(match(r, s))
end

"""
    Return a struct corresponding to a regex match's Vim command
"""
function command(m :: RegexMatch) :: Command
    args = [ parse_value(capture) for capture in m.captures ]
    return command(args...)
end

function command(c :: Char)
    return InsertCommand(c)
end

function command(n1 :: Union{Integer, Nothing},
                 motion :: Char) :: Command
    r1 = if n1 === nothing 1 else n1 end
    return MotionCommand(r1, motion)
end

function command(n1 :: Union{Integer, Nothing},
                 operator1 :: Char,
                 operator2 :: Char) :: Command
    r1 = if n1 === nothing 1 else n1 end
    if operator1 != operator2
        error("operator1 is not equal to operator2: $operator1 != $operator2")
    end
    return LineOperatorCommand(r1, operator1)
end

function command(n1 :: Union{Integer, Nothing},
                 operator :: Char,
                 n2 :: Union{Integer, Nothing},
                 motion :: Union{Char, Nothing},
                 textobject :: Union{String, Nothing}) :: Command
    r1 = if n1 === nothing 1 else n1 end
    r2 = if n2 === nothing 1 else n2 end
    if motion === nothing && textobject === nothing
        error("Both `motion` and `textobject` are empty.")
    end
    action = if motion !== nothing motion else textobject end
    return OperatorCommand(r1,
                           operator,
                           r2,
                           action)
end

function parse(cmd :: String)
    rule = matched_rule(cmd)
    @assert rule !== nothing
    m = match(rule, cmd)
    return m
end

function Base.Dict(m :: RegexMatch)
    d = OrderedDict{Symbol, Any}()
    idx_to_capture_name = Base.PCRE.capture_names(m.regex.regex)
    if !isempty(m.captures)
        for i = 1:Base.length(m.captures)
            capture_name = get(idx_to_capture_name, i, i) |> Symbol
            d[capture_name] = m.captures[i]
        end
    end
    # Dict(Symbol(n)=>m[Symbol(n)] for n in values(Base.PCRE.capture_names(m.regex.regex)))
    return d
end

function text_object_part(cmd :: AbstractString) :: Union{String, Nothing}
    m = match(REGS.text_object, cmd)
    if m === nothing
        return nothing
    end
    return m.captures[1]
end

function verb_part(cmd :: AbstractString) :: Union{Char, Nothing}
    reg = Regex("\\d*($operator).*")
    m = match(reg, cmd)
    if m === nothing
        return nothing
    end
    return m.captures[1][1]
end
end


text:
function lazy_append(a::Triangulation,b::Triangulation)
  AppendedTriangulation(a,b)
end

struct AppendedTriangulation{Dc,Dp} <: Triangulation{Dc,Dp}
  a::Triangulation{Dc,Dp}
  b::Triangulation{Dc,Dp}
end

function get_cell_coordinates(trian::AppendedTriangulation)
  a = get_cell_coordinates(trian.a)
  b = get_cell_coordinates(trian.b)
  lazy_append(a,b)
end

function get_reffes(trian::AppendedTriangulation)
  vcat(get_reffes(trian.a),get_reffes(trian.b))
end

function get_cell_type(trian::AppendedTriangulation)
  a = get_cell_type(trian.a)
  b = get_cell_type(trian.b) .+ Int8(length(get_reffes(trian.a)))
  lazy_append(a,b)
end

function reindex(f::AbstractArray,trian::AppendedTriangulation)
  a = reindex(f,trian.a)
  b = reindex(f,trian.b)
  lazy_append(a,b)
end

function get_cell_id(trian::AppendedTriangulation)
  a = get_cell_id(trian.a)
  b = get_cell_id(trian.b)
  lazy_append(a,b)
end

function restrict(f::AbstractArray,trian::AppendedTriangulation)
  a = restrict(f,trian.a)
  b = restrict(f,trian.b)
  lazy_append(a,b)
end

function get_cell_reffes(trian::AppendedTriangulation)
  a = get_cell_reffes(trian.a)
  b = get_cell_reffes(trian.b)
  lazy_append(a,b)
end

function get_cell_shapefuns(trian::AppendedTriangulation)
  a = get_cell_shapefuns(trian.a)
  b = get_cell_shapefuns(trian.b)
  lazy_append(a,b)
end

function get_cell_map(trian::AppendedTriangulation)
  a = get_cell_map(trian.a)
  b = get_cell_map(trian.b)
  lazy_append(a,b)
end

function get_normal_vector(trian::AppendedTriangulation)
  cm = get_cell_map(trian)
  a = get_array(get_normal_vector(trian.a))
  b = get_array(get_normal_vector(trian.b))
  GenericCellField(lazy_append(a,b),cm)
end

function CellQuadrature(trian::AppendedTriangulation,degree1::Integer,degree2::Integer)
  quad1 = CellQuadrature(trian.a,degree1)
  quad2 = CellQuadrature(trian.b,degree2)
  lazy_append(quad1,quad2)
end

function CellQuadrature(trian::AppendedTriangulation,degree::Integer)
  CellQuadrature(trian,degree,degree)
end


text:
using Checkers

@testset "Tests that Pass" begin
    @testset "Tests using @test_forall" begin
        @test_forall x in -1:1, x*(x-1)*(x+1) == 0
        @test_forall i in 1:2, x in -1:3, y in x+1:4, (x+i)^2<(y+i)^2
        @test_forall x in [0,1,2], y in x:4, 2*x < y+4
        @test_forall x in ["a","b"], y in ["z","w"], x*y in Set(["az","aw","bz","bw"])
    end

    @testset "Tests using @test_formany" begin
        #For each x in (10,100), for each y in (x+10,1000) sampling near x+10,
        #for each z in (y-5,∞), x+5 < z
        @test_formany 10<x<100, x+10<y<<1000, y-5<z<Inf, x+5 < z
        #Test 1000 times that log is increasing on (0,∞).
        @test_formany ntests = 1000 Inf>x>0,Inf>y>0, x<y-->log(x)<log(y)
        #Try to test 100 times, but stop after 100_000 tests, that
        #f(x) = x^2 is strictly increasing. Strictly speaking this test could
        #error, but it probably won't.
        @test_formany ntests = 100 maxtests = 100_000 0<x<10, 0<y<100, x^2>y^2-->x>y
        #Test that f(x) = x^3 is convex on (0,100)
        @test_formany 0<a<1, 0<x<100, 0<y<100, (a*x+(1-a)*y)^3<=a*x^3+(1-a)*y^3
    end

    @testset "Tests using @test_exists" begin
        @test_exists ntests = 1000 -10<x<10, x^2>99  # test passed
        @test_exists ntests = 100000 0<x<10000, x>9999 # test passed
        @test_exists 0<x<2*pi,x<y<2*pi,y<z<2*pi, sin(x)<sin(y)<sin(z)
        @test_exists ntests = 10000 0<x<2*pi,0<y<2*pi,0<z<2*pi, x<y<z --> sin(x)<sin(y)<sin(z)
        @test_exists 0<x1<2,x1<x2<2, x1<y<x2, (x1-1)^2*(x1+1)^2<(y-1)^2*(y+1)^2<(x2-1)^2*(x2+1)^2
    end

end

#TODO: Make this a doctest too.
#Test Summary:   | Pass  Total
#  Tests that Pass |   13     13

@testset "Tests that Fail" begin
    #Fails because 1 is not a root
	@test_forall x in -1:1, x*(x-2)*(x+1) == 0
    #Fails because of a situation akin to 11=x, 21.1=y, 16.5=z
	@test_formany 10<x<100, x+10<y<<1000, y-5<z<Inf, x+6 < z
    #Fails because f(x) = x^3 is not convex over (-100,100)
	@test_formany 0<a<1,-100<x<100,-100<y<100,(a*x+(1-a)*y)^3<a*x^3+(1-a)*y^3
    #Fails because it does not do enough tests to find the unlikely witness
	@test_exists ntests = 100 0<x<10000, x>9999
end

#TODO: The witness to failure should be displayed.

@testset "Tests that Error During Test (without making it to Fail)" begin
    #(Usually) errors because the antecedent will not be met 100 times.
    @test_formany ntests = 100 maxtests = 100 0<x<10, 0<y<100, x^2>y^2-->x>y #will usually fail
    #Errors due to bad expression design:
        #throws UndefVarError: x not defined,
        #because iteration over y comes earlier yet refers to value of x
	@test_forall y in x:4, x in 0:2, 2*x < y+4
    #Errors with a DomainError because it will test log(negative values)
	@test_formany -Inf<x<Inf, -Inf<y<Inf, x<y-->log(x)<log(y)
end


text:
using ADCME

function neural_net(args)
    #=
    x:      flattened mesh coordinates x-direction (PyObject: tf.Tensor)
    t:      flattened mesh coordinates time (PyObject: tf.Tensor)
    μ_w:    viscosity of wetting fluid (PyObject: tf.Tensor)
    μ_o:    viscosity of production fluid (PyObject: tf.Tensor)
    config: configuration of neural network i.e. [20, 20, 20, 3] (PyObject: tf.Tensor)
    θ:      already initialised weights and biases
    =#

    sw_nn = squeeze(fc(hcat(args.x, args.t), args.config, args.Theta)) + 1  # return neural network result

    sw_t_nn = tf.gradients(sw_nn, args.t)[1]             # first time derivative of saturation
    sw_x_nn = tf.gradients(sw_nn, args.x)[1]             # first spatial derivative of saturation
    sw_xx_nn = tf.gradients(sw_x_nn, args.x)[1]          # second spatial derivative of saturation

    λ_w = sw_nn^2 / args.μ_w                             # mobility wetting fluid
    λ_o = (1 - sw_nn)^2 / args.μ_o                       # mobility production fluid

    f_w_nn = λ_w / (λ_w + λ_o)                      # flux-function
    f_w_x_nn = tf.gradients(f_w_nn, args.x)[1]           # spatial derivative of flux-function

    return sw_nn, sw_t_nn, sw_xx_nn, f_w_x_nn
end

function get_loss(args, sw_pred, sw_t_pred, sw_xx_pred, f_w_x_pred, prod)
    #=
    sw_train:       training saturation (PyObject: tf.Tensor)
    sw_pred:        predicted saturation (PyObject: tf.Tensor)
    sw_t_pred:      predicted saturation derivative with respect to time (PyObject: tf.Tensor)
    sw_xx_pred:     predicted saturation second derivative with respect to space (PyObject: tf.Tensor)
    f_w_x_pred:     predicted flux-fluncion result derivative with respect to space (PyObject: tf.Tensor)
    ϵ:              1/Pe (PyObject: tf.Tensor)
    calib1:         calibrate saturation term (Bool)
    calib2:         calibrate flux-funcion term (Bool)
    calib3:         calibrate diffusion term (Bool)
    =#

    # define first calibration parameter
    λ_1 = Variable(ones(args.Nx))
    λ_1_multi = repeat(λ_1, args.Nt)

    # define second calibration parameter
    λ_2 = Variable(ones(args.Nx))
    λ_2_multi = repeat(λ_2, args.Nt)

    # compute loss

    # initial loss
    loss_init = sum((sw_pred[1 : args.Nx] - args.sw[1 : args.Nx])^2) / args.Nx
    # injection loss
    loss_inj =  sum((sw_pred[1 : args.Nx : (args.Nx * args.Nt)] -
                                args.sw[1 : args.Nx : (args.Nx * args.Nt)])^2) / args.Nt
    # production loss
    loss_prod = sum((sw_pred[args.Nx : args.Nx : (args.Nx * (args.Nt))] -
                                args.sw[args.Nx : args.Nx : (args.Nx * (args.Nt))])^2)/ args.Nt
    # PDE regularization
    loss_PDE = sum((sw_t_pred * λ_1_multi + f_w_x_pred * λ_2_multi - ϵ * sw_xx_pred)^2) / (args.Nx * args.Nt)

    # full loss
    if prod
        loss = loss_init + loss_inj + loss_prod + loss_PDE
    else
        loss = loss_init + loss_inj + loss_PDE
    end

    return loss, λ_1, λ_2

end

mutable struct CtxNN
    x; t; sw; μ_w; μ_o; ϵ; Theta; config; Nx; Nt
  end

function tfCtxGenNN(x, t, sw, μ_w, μ_o, ϵ, config, Nx, Nt)
    x_tf = constant(x)                          # tf.Tensor of spatial mesh
    t_tf = constant(t)                          # tf.Tensor of time mesh
    μ_w_tf = constant(μ_w)                      # tf.Tensor of wetting viscosity
    μ_o_tf = constant(μ_o)                      # tf.Tensor of production viscosity
    Θ_tf = Variable(fc_init([2; config]))       # xavier initialise weights and biases
    sw_tf = constant(sw)                        # tf.Tensor of test data
    ϵ_tf = constant(ϵ)                          # tf.Tensor of dimensionless diffusion coefficient

    return CtxNN(x_tf, t_tf, sw_tf, μ_w_tf, μ_o_tf, ϵ_tf, Θ_tf, config, Nx, Nt)
end


text:
function meets_compat_for_all_deps(working_directory::AbstractString, pkg, version)
    deps = Pkg.TOML.parsefile(joinpath(working_directory, uppercase(pkg[1:1]), pkg, "Deps.toml"))
    compat = Pkg.TOML.parsefile(joinpath(working_directory, uppercase(pkg[1:1]), pkg, "Compat.toml"))
    # First, we construct a Dict in which the keys are the package's
    # dependencies, and the value is always false.
    dep_has_compat_with_upper_bound = Dict{String, Bool}()
    @debug("We always have julia as a dependency")
    dep_has_compat_with_upper_bound["julia"] = false
    for version_range in keys(deps)
        if version in Pkg.Types.VersionRange(version_range)
            for name in keys(deps[version_range])
                if (!is_jll_name(name)) & (!is_julia_stdlib(name))
                    @debug("Found a new (non-stdlib non-JLL) dependency: $(name)")
                    dep_has_compat_with_upper_bound[name] = false
                end
            end
        end
    end
    # Now, we go through all the compat entries. If a dependency has a compat
    # entry with an upper bound, we change the corresponding value in the Dict
    # to true.
    for version_range in keys(compat)
        if version in Pkg.Types.VersionRange(version_range)
            for compat_entry in compat[version_range]
                name = compat_entry[1]
                value = compat_entry[2]
                if value isa Vector
                    if !isempty(value)
                        value_ranges = Pkg.Types.VersionRange.(value)
                        each_range_has_upper_bound = _has_upper_bound.(value_ranges)
                        if all(each_range_has_upper_bound)
                            @debug("Dependency \"$(name)\" has compat entries that all have upper bounds")
                            dep_has_compat_with_upper_bound[name] = true
                        end
                    end
                else
                    value_range = Pkg.Types.VersionRange(value)
                    if _has_upper_bound(value_range)
                        @debug("Dependency \"$(name)\" has a compat entry with an upper bound")
                        dep_has_compat_with_upper_bound[name] = true
                    end
                end
            end
        end
    end
    meets_this_guideline = all(values(dep_has_compat_with_upper_bound))
    if meets_this_guideline
        return true, ""
    else
        bad_dependencies = Vector{String}()
        for name in keys(dep_has_compat_with_upper_bound)
            if !(dep_has_compat_with_upper_bound[name])
                @error("Dependency \"$(name)\" does not have a compat entry that has an upper bound")
                push!(bad_dependencies, name)
            end
        end
        sort!(bad_dependencies)
        message = string("The following dependencies ",
                         "do not have a ",
                         "compat entry that has ",
                         "an upper bound: ",
                         join(bad_dependencies,
                              ", "),
                         ". You may find ",
                         "[CompatHelper]",
                         "(https://github.com/bcbi/CompatHelper.jl) ",
                         "helpful for keeping ",
                         "your compat entries ",
                         "up-to-date.",
                         "Note: If your package works for the current version `x.y.z` of a dependency `foo`, ",
                         "then a compat entry `foo = x.y.z` implies a compatibility upper bound ",
                         "for packages following semver. You can additionally include earlier versions ",
                         "your package is compatible with. ",
                         "See https://julialang.github.io/Pkg.jl/v1/compatibility/ for details.")
        return false, message
    end
end

function meets_patch_release_does_not_narrow_julia_compat(pkg::String,
                                                          new_version::VersionNumber;
                                                          registry_head::String,
                                                          registry_master::String)
    old_version = latest_version(pkg, registry_master)
    julia_compats_for_old_version = julia_compat(pkg, old_version, registry_master)
    julia_compats_for_new_version = julia_compat(pkg, new_version, registry_head)
    if Set(julia_compats_for_old_version) == Set(julia_compats_for_new_version)
        return true, ""
    end
    meets_this_guideline = range_did_not_narrow(julia_compats_for_old_version, julia_compats_for_new_version)
    if meets_this_guideline
        return true, ""
    else
        if (old_version >= v"1") || (new_version >= v"1")
            msg = string("A patch release is not allowed to narrow the ",
                         "supported ranges of Julia versions. ",
                         "The ranges have changed from ",
                         "$(julia_compats_for_old_version) ",
                         "(in $(old_version)) ",
                         "to $(julia_compats_for_new_version) ",
                         "(in $(new_version)).")
            return false, msg
        else
            @info("Narrows Julia compat, but it's OK since package is pre-1.0")
            return true, ""
        end
    end
end

function meets_name_length(pkg)
    meets_this_guideline = length(pkg) >= 5
    if meets_this_guideline
        return true, ""
    else
        return false, "Name is not at least five characters long"
    end
end

function meets_normal_capitalization(pkg)
    meets_this_guideline = occursin(r"^[A-Z]\w*[a-z]\w*[0-9]?$", pkg)
    if meets_this_guideline
        return true, ""
    else
        return false, "Name does not meet all of the following: starts with an uppercase letter, ASCII alphanumerics only, not all letters are uppercase.**"
    end
end

function meets_repo_url_requirement(pkg::String; registry_head::String)
    url = Pkg.TOML.parsefile(joinpath(registry_head, uppercase(pkg[1:1]), pkg, "Package.toml"))["repo"]
    meets_this_guideline = url_has_correct_ending(url, pkg)
    if meets_this_guideline
        return true, ""
    else
        return false, "Repo URL does not end with /name.jl.git, where name is the package name"
    end
end

function _invalid_sequential_version(reason::AbstractString)
    return false, "Does not meet sequential version number guideline: $reason", :invalid
end

function _valid_change(old_version::VersionNumber, new_version::VersionNumber)
    diff = difference(old_version, new_version)
    @debug("Difference between versions: ", old_version, new_version, diff)
    if diff == v"0.0.1"
        return true, "", :patch
    elseif diff == v"0.1.0"
        return true, "", :minor
    elseif diff == v"1.0.0"
        return true, "", :major
    else
        return _invalid_sequential_version("increment is not one of: 0.0.1, 0.1.0, 1.0.0")
    end
end

function meets_sequential_version_number(existing::Vector{VersionNumber}, ver::VersionNumber)
    always_assert(!isempty(existing))
    if ver in existing
        return _invalid_sequential_version("version $ver already exists")
    end
    issorted(existing) || (existing = sort(existing))
    idx = searchsortedlast(existing, ver)
    idx > 0 || return _invalid_sequential_version("version $ver less than least existing version $(existing[1])")
    prv = existing[idx]
    always_assert(ver != prv)
    nxt = thismajor(ver) != thismajor(prv) ? nextmajor(prv) :
          thisminor(ver) != thisminor(prv) ? nextminor(prv) : nextpatch(prv)
    ver <= nxt || return _invalid_sequential_version("version $ver skips over $nxt")
    return _valid_change(prv, ver)
end

function _has_no_prerelease_data(version)
    result = version.prerelease == ()
    return result
end
function _has_no_build_data(version)
    result = version.build == ()
    return result
end
_has_prerelease_data(version) = !( _has_no_prerelease_data(version) )
_has_build_data(version) = !( _has_no_build_data(version) )
_has_prerelease_andor_build_data(version) = _has_prerelease_data(version) || _has_build_data(version)

function meets_sequential_version_number(pkg::String,
                                         new_version::VersionNumber;
                                         registry_head::String,
                                         registry_master::String)
    if _has_prerelease_andor_build_data(new_version)
        return false, "Version number is not allowed to contain prerelease or build data", :invalid
    end
    _all_versions = all_versions(pkg, registry_master)
    return meets_sequential_version_number(_all_versions, new_version)
end

function meets_standard_initial_version_number(version)
    if _has_prerelease_andor_build_data(version)
        return false, "Version number is not allowed to contain prerelease or build data"
    end
    meets_this_guideline = version == v"0.0.1" || version == v"0.1.0" || version == v"1.0.0" || _is_x_0_0(version)
    if meets_this_guideline
        return true, ""
    else
        return false, "Version number is not 0.0.1, 0.1.0, 1.0.0, or X.0.0"
    end
end

function _is_x_0_0(version::VersionNumber)
    if _has_prerelease_andor_build_data(version)
        return false
    end
    result = (version.major >= 1) && (version.minor == 0) && (version.patch == 0)
    return result
end

function _generate_pkg_add_command(pkg::String,
                                   version::VersionNumber)::String
    return "Pkg.add(Pkg.PackageSpec(name=\"$(pkg)\", version=v\"$(string(version))\"));"
end

function meets_version_can_be_pkg_added(working_directory::String,
                                        pkg::String,
                                        version::VersionNumber;
                                        registry_deps::Vector{<:AbstractString} = String[])
    pkg_add_command = _generate_pkg_add_command(pkg,
                                                version)
    _registry_deps = convert(Vector{String}, registry_deps)
    code = """
        import Pkg;
        Pkg.Registry.add(Pkg.RegistrySpec(path=\"$(working_directory)\"));
        _registry_deps = $(_registry_deps);
        for regdep in _registry_deps
            Pkg.Registry.add(Pkg.RegistrySpec(url = regdep))
        end
        @info("Attempting to `Pkg.add` package...");
        $(pkg_add_command)
        @info("Successfully `Pkg.add`ed package");
        """
    before_message = "Attempting to `Pkg.add` the package"
    success_message = "Successfully `Pkg.add`ed the package"
    success_return_1 = true
    success_return_2 = ""
    failure_message = "Was not able to successfully `Pkg.add` the package"
    failure_return_1 = false
    failure_return_2 = string("I was not able to install the package ",
                              "(i.e. `Pkg.add(\"$(pkg)\")` failed). ",
                              "See the CI logs for details.")
    return _run_pkg_commands(working_directory,
                             pkg,
                             version;
                             code = code,
                             before_message = before_message,
                             success_message = success_message,
                             success_return_1 = success_return_1,
                             success_return_2 = success_return_2,
                             failure_message = failure_message,
                             failure_return_1 = failure_return_1,
                             failure_return_2 = failure_return_2)
end

function meets_version_can_be_imported(working_directory::String,
                                       pkg::String,
                                       version::VersionNumber;
                                       registry_deps::Vector{<:AbstractString} = String[])
    pkg_add_command = _generate_pkg_add_command(pkg,
                                                version)
    _registry_deps = convert(Vector{String}, registry_deps)
    code = """
        import Pkg;
        Pkg.Registry.add(Pkg.RegistrySpec(path=\"$(working_directory)\"));
        _registry_deps = $(_registry_deps);
        for regdep in _registry_deps
            Pkg.Registry.add(Pkg.RegistrySpec(url = regdep))
        end
        @info("Attempting to `Pkg.add` package...");
        $(pkg_add_command)
        @info("Successfully `Pkg.add`ed package");
        @info("Attempting to `import` package");
        import $(pkg);
        @info("Successfully `import`ed package");
        """
    before_message = "Attempting to `import` the package"
    success_message = "Successfully `import`ed the package"
    success_return_1 = true
    success_return_2 = ""
    failure_message = "Was not able to successfully `import` the package"
    failure_return_1 = false
    failure_return_2 = string("I was not able to load the package ",
                              "(i.e. `import $(pkg)` failed). ",
                              "See the CI logs for details.")
    return _run_pkg_commands(working_directory,
                             pkg,
                             version;
                             code = code,
                             before_message = before_message,
                             success_message = success_message,
                             success_return_1 = success_return_1,
                             success_return_2 = success_return_2,
                             failure_message = failure_message,
                             failure_return_1 = failure_return_1,
                             failure_return_2 = failure_return_2)
end

function _run_pkg_commands(working_directory::String,
                           pkg::String,
                           version::VersionNumber;
                           code,
                           before_message,
                           success_message,
                           success_return_1,
                           success_return_2,
                           failure_message,
                           failure_return_1,
                           failure_return_2)
    original_directory = pwd()
    tmp_dir_1 = mktempdir()
    tmp_dir_2 = mktempdir()
    atexit(() -> rm(tmp_dir_1; force = true, recursive = true))
    atexit(() -> rm(tmp_dir_2; force = true, recursive = true))
    cd(tmp_dir_1)
    # We need to be careful with what environment variables we pass to the child
    # process. For example, we don't want to pass an environment variable containing
    # our GitHub token to the child process. Because if the Julia package that we are
    # testing has malicious code in its __init__() function, it could try to steal
    # our token. So we only pass four environment variables:
    # 1. PATH. If we don't pass PATH, things break. And PATH should not contain any
    #    sensitive information.
    # 2. PYTHON. We set PYTHON to the empty string. This forces any packages that use
    #    PyCall to install their own version of Python instead of using the system
    #    Python.
    # 3. JULIA_DEPOT_PATH. We set JULIA_DEPOT_PATH to the temporary directory that
    #    we created. This is because we don't want the child process using our
    #    real Julia depot. So we set up a fake depot for the child process to use.
    # 4. R_HOME. We set R_HOME to "*".
    cmd = Cmd(`$(Base.julia_cmd()) -e $(code)`;
              env = Dict("PATH" => ENV["PATH"],
                         "PYTHON" => "",
                         "JULIA_DEPOT_PATH" => tmp_dir_2,
                         "R_HOME" => "*"))
    # GUI toolkits may need a display just to load the package
    xvfb = Sys.which("xvfb-run")
    @info("xvfb: ", xvfb)
    if xvfb !== nothing
        pushfirst!(cmd.exec, "-a")
        pushfirst!(cmd.exec, xvfb)
    end
    @info(before_message)
    cmd_ran_successfully = success(pipeline(cmd, stdout=stdout, stderr=stderr))
    cd(original_directory)
    rm(tmp_dir_1; force = true, recursive = true)
    rm(tmp_dir_2; force = true, recursive = true)
    if cmd_ran_successfully
        @info(success_message)
        return success_return_1, success_return_2
    else
        @error(failure_message)
        return failure_return_1, failure_return_2
    end
end

url_has_correct_ending(url, pkg) = endswith(url, "/$(pkg).jl.git")


text:
using MittagLeffler
using Plots

function mittleffplot(α,β,t0,t1)
    plot( x ->  mittleff(α,β,-x), t0,t1)
end

function mittleffplot2(α,β,t0,t1,arg; kws...)
    st = (t1-t0)/20
    plot( r ->  abs(mittlefferr(α,β,r*exp(im*α*arg), 10*eps())), t0:st:t1; kws...)
end

function mittleffplot3(α,β,t0,t1)
    st = (t1-t0)/1000
    plot( r ->  abs(mittleff(α,β,r*exp(im*α*pi/2))), t0:st:t1, ylims = (1.3,1.38))
end


text: abstract type BaseScore end


text:
module GraphsOpt

using StatsBase

export GraphVertex, NodeType, Person, Address,
		generate_random_graph, get_random_person,
		get_random_address, generate_random_nodes,
		convert_to_graph,
		bfs, check_euler, partition,
		graph_to_str, node_to_str,
		test_graph

# Główne optymalizacje:
# 1. Generowanie grafu bezpośrednio do postaci użytkowej,
#    bez użycia macierzy.
# 2. Użycie IOBuffer dla zamiany na string.
# 3. Przekazywanie N i K jako parametry funkcji.
# 4. Podanie konkretnych typów.
# 5. Ustalenie struktur jako niemodyfikowalne.

# Czas wykonania przed optymalizacją: ~14.0s
# Czas wykonania po optymalizacji:    ~ 0.3s

abstract type NodeType end

struct Person <: NodeType
	name ::String
end

struct Address <: NodeType
	streetNumber ::Int8
end

function str_node(n::Person)
	string("Person: ", n.name)
end

function str_node(n::Address)
	string("Street nr: ", n.streetNumber)
end

mutable struct GraphVertex
	value ::NodeType
	neighbors ::Vector{GraphVertex}
end

function get_random_person()
	Person(randstring())
end

function get_random_address()
	Address(rand(1:100))
end

function random_graph(N::Int64, K::Int64)
	graph = Vector{GraphVertex}(N)

	for i = 1:length(graph)
		node = rand() > 0.5 ? get_random_person() : get_random_address()
		graph[i] = GraphVertex(node, GraphVertex[])
	end

	for s in sample(1:N*N, K, replace=false)
		i, j = ind2sub((N, N), s)
		push!(graph[i].neighbors, graph[j])
		push!(graph[j].neighbors, graph[i])
	end

	graph
end

function is_connected(graph::Vector{GraphVertex})
	remaining = Set(graph)
	bfs(Set{GraphVertex}(), remaining)
	return isempty(remaining)
end

function bfs{T}(visited::Set{T}, remaining::Set{T})
	first = next(remaining, start(remaining))[1]

	q::Set{T} = Set()
	push!(q, first)

	push!(visited, first)
	delete!(remaining, first)

	while !isempty(q)
		v = pop!(q)

		for n in v.neighbors
			if !(n in visited)
				push!(q, n)
				push!(visited, n)
				delete!(remaining, n)
			end
		end
	end
end

function check_euler(graph::Vector{GraphVertex})
	if is_connected(graph)
		return all(map(v -> iseven(length(v.neighbors)), graph))
	end
		"Graph is not connected"
end

function graph_to_str(buf::IOBuffer, graph::Vector{GraphVertex})
	graph_str::String = ""
	for v in graph::Vector{GraphVertex}
		println(buf, "****")
		println(buf, str_node(v.value))
		print(buf, "Neighbors: ")
		println(buf, string(length(v.neighbors)))
	end
end

function test_graph()
	buf::IOBuffer = IOBuffer()

	N = 800
	K = 10000

	for i = 1:100
		graph = random_graph(N, K)

		graph_to_str(buf, graph)
		println(check_euler(graph))
	end
end

test_graph()

Profile.clear()
@profile test_graph()
@time test_graph()
Profile.print(format=:flat)

end


text:
"""
    causality(source::AbstractUncertainIndexValueDataset,
        target::AbstractUncertainIndexValueDataset,
        test::ConstrainedTest{CT, CR}) where {CT, CR}

Apply a causality test of type `CT` to `test.n` independent
realisations of `source` and `target`, after first constraining
the supports of the uncertain values furnishing the datasets.

See also [`ConstrainedTest`](@ref).
"""
function causality(source::AbstractUncertainIndexValueDataset,
    target::AbstractUncertainIndexValueDataset,
    test::ConstrainedTest{CT, CR}) where {CT, CR}

     # Take the (index_constraints, values_constraints)-tuple for the source and target
    # and represent as separate `IndexValueResampling` instances.
    cs_src = ConstrainedIndexValueResampling(test.constrained_resampling[1])
    cs_targ = ConstrainedIndexValueResampling(test.constrained_resampling[2])

    constrained_source = constrain(source, cs_src)
    constrained_target = constrain(target, cs_targ)

    [causality(
            resample(constrained_source)[2],
            resample(constrained_target)[2],
            test.test) for i = 1:test.n]
end

export causality


text:
# Import APM package
include("../apm.jl")

# Solve optimization problem
sol = apm_solve("hs71",3)

println("--- Results of the Optimization Problem ---")
println("Solution Results")
println("  x[1] = " * string(sol["x[1]"]))
println("  x[2] = " * string(sol["x[2]"]))
println("  x[3] = " * string(sol["x[3]"]))
println("  x[4] = " * string(sol["x[4]"]))


text:
function build_url_query_string(base::String, options::Dict{String,String})::String

    # init -
    url_string = base

    parameters = ""
    for (key,value) in options
        parameters *= "$(key)=$(value)&"
    end

    # cut off trailing &
    query_parameters = parameters[1:end-1]

    # return -
    return url_string*query_parameters
end

function check(result::Some)::(Union{Nothing,T} where {T<:Any})

    # ok, so check, do we have an error object?
    # Yes: log the error if we have a logger, then throw the error.
    # No: return the result.value

    # Error case -
    if (isa(something(result), Exception) == true)

        # get the error object -
        error_object = result.value

        # get the error message as a String -
        error_message = sprint(showerror, error_object, catch_backtrace())
        @error(error_message)

        # throw -
        throw(result.value)
    end

    # default -
    return result.value
end

function process_csv_api_data(api_call_raw_data::String)::DataFrame

    # create a data table from the CSV data -
    tmp_data_table = CSV.read(IOBuffer(api_call_raw_data), DataFrame)

    # sort the table according to the timestamps -
    idx_sort = sortperm(tmp_data_table[:,1])

    # create a sorted data table -
    sorted_data_table = tmp_data_table[idx_sort,:]

    # return -
    return sorted_data_table
end

max_stars_repo_name: Pooksoft/BinomialEquityAssetPriceModel

max_issues_repo_name: Pooksoft/BinomialEquityAssetPriceModel

max_forks_repo_name: Pooksoft/BinomialEquityAssetPriceModel


text:
using FdeSolver
using Test
using SpecialFunctions

@testset "FdeSolver.jl" begin

    tSpan = [0, 5]
    y0 = [1, 0.5, 0.3]
    β = [0.5, 0.2, 0.6]

    function F(t, y)

        F1 = 1 / sqrt(pi) * (((y[ 2] - 0.5) * (y[3] - 0.3))^(1 / 6) + t^(1 / 2))
        F2 = gamma(2.2) * (y[1] - 1)
        F3 = gamma(2.8) / gamma(2.2) * (y[2] - 0.5)

        return [F1, F2, F3]

    end

    t, Yapp = FDEsolver(F, tSpan, y0, β)

    @test @isdefined(t)
    @test @isdefined(Yapp)

end


text:
module CLSolvers

using LinearAlgebra
using Random
using Distributions
#using CUDA
using LinearAlgebra
using SparseArrays
using ForwardDiff
using ShiftedArrays, LabelledArrays

using DocStringExtensions

include("model.jl")
include("contour.jl")
include("contourDiscretization.jl")
include("analysis.jl")
include("plot_observables.jl")
include("sol_helpers.jl")

##################
#   EXPORT
##################
#export Plots

end