text:
#
# This file is part of the YAActL.jl Julia package, MIT license
#
# Paul Bayer, 2020
#

function _send!(chn::Channel, m::Message)
    # reimplements Base.put_buffered with a modification
    lock(chn)
    try
        while length(chn.data) ≥ chn.sz_max  # modification: allow buffer overflow
            Base.check_channel_state(chn)
            wait(chn.cond_put)
        end
        push!(chn.data, m)
        # notify all, since some of the waiters may be on a "fetch" call.
        notify(chn.cond_take, nothing, true, false)
    finally
        unlock(chn)
    end
    return m
end
function _send!(rch::RemoteChannel, m::Message)
    if rch.where != myid()
        # change any local links in m to remote links
        m = typeof(m)((_rlink(getfield(m,i)) for i in fieldnames(typeof(m)))...)
    end
    put!(rch, m)
end

"""
```
send!(lk::Link, m::Message)
send!(name::Symbol, m::Message)
```
Send a message `m` to an actor `lk` or `name` (if registered).
"""
send!(lk::L, m::Message) where L<:Link = _send!(lk.chn, m)
send!(name::Symbol, m::Message) = _send!(whereis(name).chn, m)

# """
# ```
# send!(lks::Tuple{Link,Vararg{Link}}, m::M) where M<:Message
# send!(lks::Vector{Link}, m::M) where M<:Message
# ```
# Send a message `m` to a `Vector` or `Tuple` of [`Link`](@ref)s.
# """
# send!(lks::Tuple{Link,Vararg{Link}}, m::M) where M<:Message =
#     map(x->send!(x, m), lks)
# send!(lks::Vector{Link}, m::M) where M<:Message =
#     map(x->send!(x, m), lks)

_match(msg::Message, ::Nothing, ::Nothing) = true
_match(msg::Message, Msg::Type{<:Message}, ::Nothing) = msg isa Msg
_match(msg::Message, ::Nothing, from::Link) =
    :from in fieldnames(typeof(msg)) ? msg.from == from : false
function _match(msg::Message, Msg::Type{<:Message}, from::Link)
    if :from in fieldnames(typeof(msg))
        return msg isa Msg && msg.from == from
    else
        return false
    end
end

"""
```
receive!(lk; timeout=5.0)
receive!(lk, from; timeout=5.0)
receive!(lk, Msg; timeout=5.0)
receive!(lk, Msg, from; timeout=5.0)
```
Receive a message over a link `lk`.

If `Msg` or `from` are provided, `receive!` returns only a
matching message. Other messages in `lk` are restored to it in their
previous order.

# Parameters
- `lk::Link`: local or remote link over which the message is received,
- `Msg::Type{<:Message}`: [`Message`](@ref) type,
- `from::Link`: local or remote link of sender. If `from` is
    provided, only messages with a `from` field can be matched.
- `timeout::Real=5.0`: maximum waiting time in seconds.
    - If `timeout==0`, `lk` is scanned only for existing messages.
    - Set `timeout=Inf` if you don't want to timeout.

# Returns
- received message or `Timeout()`.
"""
receive!(lk::L; kwargs...) where L<:Link = receive!(lk, nothing, nothing; kwargs...)
receive!(lk::L, from::Link; kwargs...) where L<:Link = receive!(lk, nothing, from; kwargs...)
receive!(lk::L, Msg::Type{<:Message}; kwargs...) where L<:Link = receive!(lk, Msg, nothing; kwargs...)
function receive!(lk::L1, Msg::M, from::L2;
    timeout::Real=5.0) where {L1<:Link,M<:Union{Nothing,Type{<:Message}},L2<:Union{Nothing,Link}}

    done = [false]
    msg = Timeout()
    stash = Message[]
    ev = Base.Event()
    timeout > 0 && !isinf(timeout) && Timer(x->notify(ev), timeout)

    @async begin
        while !done[1]
            timeout == 0 && !isready(lk.chn) && break
            _match(fetch(lk.chn), Msg, from) && break
            done[1] || push!(stash, take!(lk.chn))
        end
        notify(ev)
    end

    wait(ev)
    done[1] = true
    isready(lk.chn) && (msg = take!(lk.chn))
    while !isempty(stash) && isready(lk.chn)
        push!(stash, take!(lk.chn))
    end
    foreach(x->put!(lk.chn, x), stash)
    return msg
end

"""
```
request!(lk::Link, msg::Message; full=false, timeout::Real=5.0)
request!(lk::Link, Msg::Type{<:Message}, args...; kwargs...)
request!(name::Symbol, args...; kwargs...)
```
Send a message to an actor, block, receive and return the result.

# Arguments
- `lk::Link`: actor link, or `name::Symbol` (if registered),
- `msg::Message`: a message,
- `Msg::Type{<:Message}`: a message type,
- `args...`: optional arguments to `Msg`,
- `full`: if `true` return the full [`Response`](@ref) message.
- `timeout::Real=5.0`: timeout in seconds after which a
    [`Timeout`](@ref) is returned,
- `kwargs...`: `full` or `timeout`.

"""
function request!(lk::L, msg::Message;
                full=false, timeout::Real=5.0) where L<:Link
    send!(lk, msg)
    resp = receive!(msg.from, timeout=timeout)
    return resp isa Timeout || full ? resp : resp.y
end
function request!(lk::L, Msg::Type{<:Message}, args...; kwargs...)  where L<:Link
    me = lk isa Link{Channel{Message}} ?
            Link(1) :
            Link(RemoteChannel(()->Channel{Message}(1)), myid(), :remote)
    return Msg in (Exec, Query) ?
        request!(lk, Msg(args..., me); kwargs...) :
        request!(lk, isempty(args) ? Msg(me) : Msg(args, me); kwargs...)
end
request!(name::Symbol, args...; kwargs...) = request!(whereis(name), args...; kwargs...)


text:
"""
Sparcity of functions
"""

abstract type AbstractSparsityPattern end

struct DensePattern <: AbstractSparsityPattern end

struct SparsePattern{TI} <: AbstractSparsityPattern
    rows::Vector{TI}
    cols::Vector{TI}
end

"""
    SparsePattern(A::SparseMatrixCSC)
construct sparse pattern from representative sparse matrix
# Arguments
- `A::SparseMatrixCSC`: sparse jacobian
"""
function SparsePattern(A::SparseMatrixCSC)
    rows, cols, _ = findnz(A)
    return SparsePattern(rows, cols)
end

"""
    SparsePattern(A::Matrix)
construct sparse pattern from representative matrix
# Arguments
- `A::Matrix`: sparse jacobian
"""
function SparsePattern(A::Matrix)
    return SparsePattern(sparse(A))
end

"""
    SparsePattern(::FD, func!, ng, x1, x2, x3)
detect sparsity pattern by computing derivatives (using finite differencing)
at three different locations. Entries that are zero at all three
spots are assumed to always be zero.
# Arguments
- `func!::Function`: function of form f = func!(g, x)
- `ng::Int`: number of constraints
- `x1,x2,x3::Vector{Float}`:: three input vectors.
"""
function SparsePattern(dtype::FD, func!, ng, x1, x2, x3)

    fdtype = finitediff_type(dtype)
    cache = FiniteDiff.JacobianCache(x1, zeros(ng), fdtype)

    nx = length(x1)
    J1 = zeros(ng, nx)
    J2 = zeros(ng, nx)
    J3 = zeros(ng, nx)
    FiniteDiff.finite_difference_jacobian!(J1, func!, x1, cache)
    FiniteDiff.finite_difference_jacobian!(J2, func!, x2, cache)
    FiniteDiff.finite_difference_jacobian!(J3, func!, x3, cache)

    @. J1 = abs(J1) + abs(J2) + abs(J3)
    Jsp = sparse(J1)

    return SparsePattern(Jsp)
end

#  used internally to get rows and cols for dense jacobian
function _get_sparsity(::DensePattern, nx, nf)
    len = nf*nx
    rows = [i for i = 1:nf, j = 1:nx][:]
    cols = [j for i = 1:nf, j = 1:nx][:]
    return rows, cols
end

#  used internally to get rows and cols for sparse jacobian
_get_sparsity(sp::SparsePattern, nx, nf) = sp.rows, sp.cols

ext: jl

max_stars_repo_name: Yuricst/joptimise

max_issues_repo_name: Yuricst/joptimise

max_forks_repo_name: Yuricst/joptimise


text:
#InspectDR: Plot generation templates
#-------------------------------------------------------------------------------

#Want more resolution on y-axis than default:
#TODO: is there a better way???
_yaxisscale(s::Symbol) = AxisScale(s, tgtmajor=8, tgtminor=2)

function Plot2D(xscale::Symbol, yscalelist::Vector{Symbol}; title::String="",
	xlabel::String="", ylabels::Vector{String}=String[], kwargs...)
	plot = InspectDR.Plot2D(;kwargs...)
	plot.xscale = AxisScale(xscale)

	#TODO: remove restriction once InspectDR deals better with 0 strips:
	if length(yscalelist) < 1
		yscalelist = [:lin]
	end

	plot.strips = [] #Reset
	for yscale in yscalelist
		strip = InspectDR.GraphStrip()
		push!(plot.strips, strip)
		strip.yscale = _yaxisscale(yscale)
	end

	a = plot.annotation
		a.title = title
		a.xlabel = xlabel
		a.ylabels = ylabels

	return plot
end
#Single y strip:
Plot2D(xscale::Symbol, yscale::Symbol; kwargs...) = Plot2D(xscale, [yscale]; kwargs...)

#Generate empty Bode plot:
function bodeplot(;title="Bode Plot", xlabel="Frequency (Hz)",
	ylabels=["Magnitude (dB)", "Phase (°)"], kwargs...)
	plot = InspectDR.Plot2D(:log10, [:dB20, :lin], title=title,
		xlabel=xlabel, ylabels=ylabels; kwargs...
	)

	#Control grid visibility:
	for strip in plot.strips
		strip.grid = InspectDR.GridRect(vmajor=true, vminor=true, hmajor=false)
	end

	return plot
end

#Generate empty transient plot:
function transientplot(yscalelist::Vector{Symbol}; title="Time-Domain Plot",
	xlabel="Time (s)", kwargs...)
	plot = Plot2D(:lin, yscalelist, title=title, xlabel=xlabel; kwargs...)
	plot.layout[:enable_legend] = true
	#plot.layout[:halloc_legend] = 150

	for graph in plot.strips
		graph.grid = InspectDR.GridRect(vmajor=true, vminor=true, hmajor=true)
	end

	#Use SI notation on both x & y axes:
	plot.layout[:format_xtick] = TickLabelStyle(NumericIO.UEXPONENT_SI)
	plot.layout[:format_ytick] = TickLabelStyle(NumericIO.UEXPONENT_SI)
	return plot
end
transientplot(yscale::Symbol=:lin; kwargs...) = transientplot([yscale]; kwargs...)

#Generate empty Smith chart:
function smithchart(t::Symbol=:Z; ref::Number=1.0, title="Smith Chart",
	xlabel = "Real(Γ)", ylabels = ["Imaginary(Γ)"], kwargs...)
	smithext = InspectDR.PExtents1D(min=-1.1,max=1.1) #WANTCONST A bit of padding
	plot = InspectDR.Plot2D(:lin, :lin, title=title, xlabel=xlabel, ylabels=ylabels; kwargs...)
	graph = plot.strips[1]
	graph.grid = InspectDR.GridSmith(t, ref=ref)
	plot.xext_full = smithext
	graph.yext_full = smithext
	return plot
end

#Last line


text:
export qromp
# Based on the algorithm developed by
# Baptista, R., Stolbunov, V., & Nair, P. B. (2019).
# Some greedy algorithms for sparse polynomial chaos expansions.
# Journal of Computational Physics, 387, 303-325.

# Solve ψc = u using the orthogonal matching pursuit with QR formulation

function qromp(ψ::AbstractMatrix{T}, u::AbstractVector{T}; invert::Bool=true, verbose::Bool = true, ϵrel::Float64 = 1e-1, maxterms::Int64=typemax(Int64)) where {T}
    m, n = size(ψ)

    residue = copy(u)
    idxset = Int64[]
    dict = collect(1:n)
    ϵu = norm(u)
    ϵhist = [ϵu]
    ϵrel *= ϵu
    ϕ = zeros(T, n)

    @inbounds for j=1:n
        for i=1:m
        ϕ[j] += ψ[i,j]^2
        end
    end

    # F = qrfactUnblocked(zeros(0,0))
    # F = elasticqrfact(zeros(m,0))
    cache = zeros(T, m)
    F = elasticqrfact(reshape(cache, length(cache),1))
    ek = zeros(T, n)
    # ektest = zeros(T, n)
    q = zeros(T, m)
    new_idx = 1

    @inbounds for k=1:n
        # Step 6 & 7 update denominators and choose best candidate

        entry = 0.0
        @fastmath @inbounds for (j, idx) in enumerate(dict)
            ψj = view(ψ,:,idx)
            ratio = dot(residue, ψj)
            ratio = ratio^2
            ratio /= ϕ[idx]

            if ratio >= entry
                entry = copy(ratio)
                new_idx = copy(idx)
            end
        end

        # if k>1
        #     for (j, idx) in enumerate(dict)
        #         ψj = view(ψ,:,idx)
        #         @show norm(ϕ[idx] - norm(ψj-Matrix(F.Q)[:,1:k-1]*(Matrix(F.Q)[:,1:k-1]'*ψj))^2)
        #     end
        # end

        # Step 8 Update set of selected basis
        push!(idxset, new_idx)

        # Step 9 Update candidate dictionary
        filter!(x-> x != new_idx, dict)

        # Step 10 Update Q(k-1) and R(k-1) with ψi(k)
        if k>1
            # F = updateqrfactUnblocked!(F, view(ψ,:,new_idx))
            updateelasticqrfact!(F, view(ψ,:,new_idx), cache)

            ek[k-1] = 0.0
        else
            # F = qrfactUnblocked(ψ[:,new_idx:new_idx])
            F = elasticqrfact(view(ψ,:,new_idx:new_idx))
        end

        # Step 11 Update residual
        ek[k] = 1.0
        mul!(q, F.Q, view(ek,1:k))

        @inbounds for (j, idx) in enumerate(dict)
            ψj = view(ψ,:,idx)
            ϕj = dot(q, ψj)
            ϕ[idx] -= ϕj^2
        end
        factor = dot(q, u)
        residue -= factor*q

        # Step 12 Calculate stopping critera
        ϵ = norm(residue)

        if verbose == true
            push!(ϵhist, copy(ϵ))
        end
        if ϵ < ϵrel || k == maxterms
            break
        end
    end

    # Solve the system with the set of indices
    if verbose == true
        if invert == true
            c = zeros(length(idxset))
            ldiv!(c, F, u)
            return idxset, c, ϵhist
        else
            return idxset, ϵhist
        end
    else
        if invert == true
            c = zeros(length(idxset))
            ldiv!(c, F, u)
            return idxset, c
        else
            return idxset
        end
    end
end

max_forks_repo_name: mleprovost/QROMP


text:
module TestPlotter

using Test
using Plots
using StatsBase: sample
using TSML

default(show=false, reuse=true)

function generatedf()
    Random.seed!(123)
    gdate = DateTime(2014,1,1):Dates.Minute(15):DateTime(2014,1,5)
    gval = Array{Union{Missing,Float64}}(rand(length(gdate)))
    gmissing = floor(0.30*length(gdate)) |> Integer
    gndxmissing = Random.shuffle(1:length(gdate))[1:gmissing]
    df = DataFrame(Date=gdate,Value=gval)
    df.Value[gndxmissing] .= missing
    return df
end;

function test_artificialdataplotter()
  Random.seed!(123)
  mdates = DateTime(2017,1,1):Dates.Hour(1):DateTime(2017,6,1)
  mvals = rand(1:1000,length(mdates))
  # create some outliers
  soutliers = rand([500:10000;-10000:500],div(length(mdates),10))
  soutndx = sample(1:length(mdates),length(soutliers))
  mvals[soutndx] = soutliers
  df = DataFrame(Date=mdates,Value=mvals)
  pltr = TSML.Plotter(Dict(:interactive => false))
  fit!(pltr,df)
  myplot=transform!(pltr,df);
  @test isa(myplot,Plots.Plot) == true;
  df = generatedf()
  fit!(pltr,df)
  myplot1=transform!(pltr,df);
  @test isa(myplot1,Plots.Plot) == true;
end
@testset "Plotter: using artificial data" begin
  test_artificialdataplotter()
end

function test_realdataplotter()
  Random.seed!(123)
  fname = joinpath(dirname(pathof(TSML)),"../data/testdata.csv")
  csvfilter = CSVDateValReader(Dict(:filename=>fname,:dateformat=>"dd/mm/yyyy HH:MM"))
  pltr = TSML.Plotter(Dict(:interactive => false))
  mpipeline1 = @pipeline csvfilter |> pltr
  myplot = fit_transform!(mpipeline1);
  @test isa(myplot,Plots.Plot) == true;
end
@testset "Plotter: readcsv |> plotter" begin
  test_realdataplotter()
end

end


text:
using BenchmarkTools

##

function benchmark_broadcasting()
    shape, wcs = fullsky_geometry(deg2rad(1); dims=(3,))
    A, B = rand(shape...), rand(shape...)
    ma = Enmap(A, wcs)
    mb = Enmap(B, wcs)

    @btime $mb .* $ma .* exp.($ma.^2)
    @btime $B .* $A .* exp.($A.^2)

    return (mb .* ma .* exp.(ma.^2)) == parent(mb .* ma .* exp.(ma.^2))

end

benchmark_broadcasting()

##
using Pixell, BenchmarkTools

function benchmark_sky2pix()
    shape, wcs = fullsky_geometry(deg2rad(1); dims=(3,))
    A, B = rand(shape...), rand(shape...)
    ma = Enmap(A, wcs)
    mb = Enmap(B, wcs)

    @btime sky2pix($ma, 5, 6)

end

benchmark_sky2pix()


text:
#############################################################################
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#
# Author:  Louis Luangkesorn <lugerpitt@gmail.com>
# Date:  February 26, 2015
#
# This model determines a set of workforce levels that will most economically
# meet demands and inventory requirements over time. The formulation is
# motivated by the experiences of a large producer in the United States.  The
# data are for three products and 13 numperiods
#
# Problem taken from the Appendix C of the expanded version of Fourer, Gay,
# and Kernighan, A Modeling Language for Mathematical Programming
#############################################################################

function PrintSolution(status, CREWS, HIRE, LAYOFF)
    println("RESULTS:")
    if status == :Optimal
      println("Crews")
      for t = 0:length(CREWS)-1
        print(" $(getvalue(CREWS[t])) ")
      end
      println()
      println("Hire")
      for t = 1:length(HIRE)
        print(" $(getvalue(HIRE[t])) ")
      end
      println()
      println("Layoff")
      for t = 1:length(LAYOFF)
        print(" $(getvalue(LAYOFF[t])) ")
      end
      println()
    else
        println("  No solution")
    end
    println("")
end

using JuMP

####  PRODUCTION SETS AND PARAMETERS  ###

prd = ["18REG"  "24REG" "24PRO"]
# Members of the product group
numprd = length(prd)
pt=	[1.194,	1.509,	1.509]
# Crew-hours to produce 1000 units
pc=	[2304,	2920,	2910]
# Nominal production cost per 1000, used
# to compute inventory and shortage costs

###  TIME PERIOD SETS AND PARAMETERS  ###

firstperiod = 1
# Index of first production period to be modeled
lastperiod  = 13
# Index of last production period to be modeled
numperiods = firstperiod:lastperiod
# 'planning horizon' := first..last;

###  EMPLOYMENT PARAMETERS  ###

cs = 18
# Workers per crew
sl =  8
# Regular-time hours per shift
rtr = 16.00
# Wage per hour for regular-time labor
otr = 43.85
# Wage per hour for overtime labor
iw =  8
# Crews employed at start of first period
dpp =	 [19.5,	19,	20,	19,	19.5,	19,	19,	20,	19,	20,	20,	18,	18]
# Regular working days in a production period
ol =	 [96,	96,	96,	96,	96,	96,	96,	96,	96,	96,	96,	96,	96]
# Maximum crew-hours of overtime in a period
cmin =	[0,	0,	0,	0,	0,	0,	0,	0,	0,	0,	0,	0,	0]
# Lower limit on average employment in a period
cmax =	[8,	8,	8,	8,	8,	8,	8,	8,	8,	8,	8,	8,	8]
# Upper limit on average employment in a period
hc =	 [7500,	7500,	7500,	7500,	15000,	15000,	15000,	15000,	15000,	15000,	7500,	7500,	7500]
# Penalty cost of hiring a crew
lc =	 [7500,	7500,	7500,	7500,	15000,	15000,	15000,	15000,	15000,	15000,	7500,	7500,	7500]
# Penalty cost of laying off a crew

###  DEMAND PARAMETERS  ###

d18REG=[63.8,	76,	88.4,	913.8,	115,	133.8,	79.6,	111,	121.6,	470,	78.4,	99.4,	140.4,	63.8]
d24REG = [1212,	306.2,	319,	208.4,	298,	328.2,	959.6,	257.6,	335.6,	118,	284.8,	970,	343.8,	1212]
d24PRO = [0,	0,	0,	0,	0,	0,	0,	0,	0,	1102,	0,	0,	0,	0]
dem = Array[d18REG, d24REG, d24PRO]
# Requirements (in 1000s)
# to be met from current production and inventory

pro = Array[
  [0,	0,	0,	1,	0,	0,	0,	0,	0,	1,	0,	0,	0,	0],
  [1,	0,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0,	1,	1],
  [0,	0,	0,	0,	0,	0,	0,	0,	0,	1,	0,	0,	0,	0]]
# true if product will be the subject
# of a special promotion in the period

###  INVENTORY AND SHORTAGE PARAMETERS  ###

rir =  0.75
# Proportion of non-promoted demand
# that must be in inventory the previous period
pir =  0.80
# Proportion of promoted demand
# that must be in inventory the previous period
life  =  2
# Upper limit on number of periods that
# any product may sit in inventory
cri	=[0.015,	0.015,	0.015]
# Inventory cost per 1000 units is
# cri times nominal production cost
crs	=[1.1,	1.1,	1.1]
# Shortage cost per 1000 units is
# crs times nominal production cost
iinv	=[82,	792.2,	0]
# Inventory at start of first period; age unknown
iil = [[max(0,iinv[p] - sum([dem[p][v] for v=firstperiod:t]))
        for t = numperiods]
        for p=1:numprd]
# Initial inventory still available for allocation
# at end of period t
function checkpro(product, timeperiod, production, promotionalrate, regularrate)
  if production[product][timeperiod + 1]==1
    value=promotionalrate
  else
    value =regularrate
  end
  value
end
minv = [[dem[p][t+1] * checkpro(p,t, pro, pir, rir) for t=numperiods] for p=1:numprd]
# Lower limit on inventory at end of period t

prod = Model()

###  VARIABLES  ###
@variable(prod, Crews[0:lastperiod] >= 0)
# Average number of crews employed in each period

@variable(prod, Hire[numperiods] >= 0)
# Crews hired from previous to current period

@variable(prod, Layoff[numperiods]>= 0)
# Crews laid off from previous to current period

@variable(prod, Rprd[1:numprd, numperiods] >=0)
# Production using regular-time labor, in 1000s

@variable(prod, Oprd[1:numprd, numperiods]>=0)
# Production using overtime labor, in 1000s

@variable(prod, Inv[1:numprd, numperiods, 1:life] >=0)
# a numperiods old -- produced in period (t+1)-a --
# and still in storage at the end of period t
@variable(prod, Short[1:numprd, numperiods]>=0)
			# Accumulated unsatisfied demand at the end of period t

###  CONSTRAINTS  ###

@constraint(prod, xyconstr[t=numperiods],
               sum{pt[p] * Rprd[p,t], p=1:numprd} <= sl * dpp[t] * Crews[t])
# Hours needed to accomplish all regular-time
# production in a period must not exceed
# hours available on all shifts

@constraint(prod, xyconstr[t=numperiods],
               sum{pt[p] * Oprd[p,t], p=1:numprd}  <= ol[t])
# Hours needed to accomplish all overtime
# production in a period must not exceed
# the specified overtime limit

@constraint(prod, Crews[firstperiod-1]==iw)
# Use given initial workforce

@constraint(prod,xyconstr[t=numperiods],
               Crews[t]== Crews[t-1] + Hire[t]-Layoff[t])
# Workforce changes by hiring or layoffs

@constraint(prod, xyconstr[t=numperiods], cmin[t] <= Crews[t])
@constraint(prod, xyconstr[t=numperiods], Crews[t] <= cmax[t])
# Workforce must remain within specified bounds

@constraint(prod, xyconstr[p=1:numprd],
               Rprd[p, firstperiod] + Oprd[p,firstperiod] + Short[p, firstperiod]-Inv[p, firstperiod, 1]
                 == max(0,dem[p][firstperiod]-iinv[p]))
# 'first demand requirement

# NOTE: JuMP xyconstr[] requires that indices be integer at compile time,
# so firstperiod +1 could not be an index within xycontr or triconstr
for t=(firstperiod+1:lastperiod)
  @constraint(prod, xyconstr[p=1:numprd],
                 Rprd[p,t] + Oprd[p,t] + Short[p,t] - Short[p,t-1] +
                   sum{Inv[p, t-1, a] - Inv[p,t,a], a=1:life} ==
                   max(0, dem[p][t]-iil[p][t-1]))
end
# Production plus increase in shortage plus
# decrease in inventory must equal demand

@constraint(prod, xyconstr[p=1:numprd, t=numperiods],
               sum{Inv[p,t,a] + iil[p][t], a=1:life} >= minv[p][t])
# Inventory in storage at end of period t
# must meet specified minimum

@constraint(prod, triconstr[p=1:numprd, v=1:(life-1), a=v+1:life], Inv[p, firstperiod+v-1, a] ==0)

# In the vth period (starting from first)
# no inventory may be more than v numperiods old
# (initial inventories are handled separately)

@constraint(prod, xyconstr[p=1:numprd, t=numperiods],
               Inv[p,t,1] <= Rprd[p,t]+Oprd[p,t])
# New inventory cannot exceed
# production in the most recent period

secondperiod = firstperiod + 1
@constraint(prod, triconstr[p=1:numprd, t=2:lastperiod, a=2:life],
                     Inv[p,t,a] <= Inv[p,t-1,a-1])
# Inventory left from period (t+1)-p
# can only decrease as time goes on

###  OBJECTIVE  ###

@objective(prod, Min,
                sum{rtr * sl * dpp[t] * cs * Crews[t], t=numperiods} +
                sum{hc[t] * Hire[t], t=numperiods} +
                sum{lc[t] * Layoff[t], t=numperiods} +
                sum{sum{ otr * cs * pt[p] * Oprd[p,t], t=numperiods}, p=1:numprd} +
                sum{sum{sum{cri[p] * pc[p] * Inv[p,t,a], t=numperiods}, p=1:numprd}, a=1:life} +
                sum{sum{ crs[p] * pc[p] * Short[p,t], t=numperiods},p=1:numprd})

# Full regular wages for all crews employed, plus
# penalties for hiring and layoffs, plus
# wages for any overtime worked, plus
# inventory and shortage costs

# (All other production costs are assumed
# to depend on initial inventory and on demands,
# and so are not included explicitly.)

production = solve(prod)

PrintSolution(production, Crews, Hire, Layoff)


text:
module Space

using ..Basis
import ..Tagging: hastag, traithastag, tag
import ..Traits: traits

# Sample / Parameter Space
# include("abstractomega.jl")     # Abstract sample/Paramter Spaces
include("pdf.jl")               # Pdf
# include("scope.jl")
include("simpleomega.jl")       # Sample Space / Distributions
include("lazyomega.jl")         # Sample Space / Distributions

# Defaults
Basis.defΩ(args...) = LazyΩ{EmptyTags, Dict{Any, Any}}
Basis.defω(args...) = tagrng(defΩ()(), Random.GLOBAL_RNG)

end


text:
requiredTypes = [ "AcademicLevel", "Affiliation", "Assignment", "CareerStatus", "Job", "Rank", "SubJob", "CareerPath" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( dirname( Base.source_path() ), "functions", reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end

XlsxFile = joinpath( dirname( Base.source_path() ), "..", "SimulationInput.xlsx" )#____________________________________________ Excel File

include( "types/DataStructures.jl" )

include( "OptCareerPaths/ReadDataXlsx.jl" )

include("OptCareerPaths/PLresolution.jl")

include("OptCareerPaths/WrittingResults.jl")

if PlottingResults
    include("OptCareerPaths/PlottingResults.jl")
end

if AgeDistPlot
    include("OptCareerPaths/ComputeAgeDist.jl")
end

max_stars_repo_name: Omazaria/CareerPathBeDef

max_issues_repo_name: Omazaria/CareerPathBeDef

max_forks_repo_name: Omazaria/CareerPathBeDef


text:
@testset "SemialgebraicSets" begin
    solver = SemialgebraicSetsHCSolver(; compile = false)
    @polyvar x y
    V = SemialgebraicSets.@set x^2 == 1 && y^2 == 2 solver
    S = collect(V)
    S = sort!(map(s -> round.(s, digits = 2), S))
    @test S == [[-1.0, -1.41], [-1.0, 1.41], [1.0, -1.41], [1.0, 1.41]]
end


text:
# Optimisation parameters
mutable struct Options{F <: Function}
    f::F                    # Function to optimise
    N::Int64                # Dimension of x
    Ns::Int64               # Number of cycles before step length adjustment
    Nt::Int64               # Number of cycles before temperature reduction/termination check
    Ne::Int64               # Number of temperature reductions without increase > tol before terminating
    lb::Vector{Float64}     # Lower bound of x
    ub::Vector{Float64}     # Upper bound of x
    R::Float64              # Temperature reduction factor
    c::Vector{Float64}      # Step length adjustment factor
    max::Bool               # Maximise/minimise function
    tol::Float64            # Tolerance of improvement on cost
    max_eval::Int64         # Maximum number of evaluations before termination
    print_status::Bool      # true if print status at each temp reduction
end
# Default constructor
function Options(;func=f,
                  N=9,
                  Ns=20,
                  Nt=10,
                  Ne=4,
                  lb=repeat([-10.0], N),
                  ub=repeat([10.0], N),
                  R=0.75,
                  c=repeat([2.0], N),
                  max=false,
                  tol=1e-6,
                  max_eval=1_000_000,
                  print_status=false)

    return Options(func, N, Ns, Nt, Ne, lb, ub, R, c, max, tol, max_eval, print_status)
end

# Result of optimisation
mutable struct Result
    fopt::Float64
    xopt::Vector{Float64}
    terminate::Symbol           # Termination code
    n_moves::Int64              # Total number of moves
    n_accepted::Int64           # Total number of accepted points
    n_optimums::Int64           # Total number of global optimums found

    Result(fopt, xopt) = new(fopt, xopt, :none, 0, 0, 0)
end
# Default constructor
function Result(;fopt=3208,
                 xopt=repeat([2.0], 9))
    return Result(fopt, copy(xopt))
end

# Current state of optimisation
mutable struct State
    f::Float64                   # Current f
    x::Vector{Float64}           # Current x
    v::Vector{Float64}           # Current step length
    T::Float64                   # Current temperature
    n_accepted::Vector{Float64}  # Number of accepted points for current step length (reset after every step length adjustment)
    ftemp::Float64
    xtemp::Vector{Float64}

    State(f, x, v, T) = new(f, x, v, T, zeros(Int64, length(x)), 0.0, zeros(length(x)))
end
# Default constructor
function State(;f=3208,
                x=repeat([2.0], 9),
                v=repeat([2.0], 9),
                T=1.0)
    return State(f, copy(x), copy(v), T)
end

Base.copy(state::State) = State(state.f, state.x, state.v, state.T)
Base.copy(result::Result) = Result(result.fopt, result.xopt)
Base.copy(options::Options) = Options(options.f, options.N, options.Ns, options.Nt, options.Ne, options.lb, options.ub, options.R, options.c, options.max, options.tol, options.max_eval, options.print_status)

Base.show(io::IO, state::State) = print(io,
"State:
    f: $(state.f)
    x: $(state.x)
    v: $(state.v)
    T: $(state.T)"
)

# function Base.show(io::IO, state::State)
#     compact = get(io, :compact, false)

#     if compact
#         print(io,
#         "State:
#             f: $(state.f)
#             x: [$(state.x[1]), $(state.x[2]), ... , $(state.x[end - 1]), $(state.x[end])]
#             v: [$(state.v[1]), $(state.v[2]), ... , $(state.v[end - 1]), $(state.v[end])]
#             T: $(state.T)"
#         )
#     else
#         print(io,
#         "State:
#             f: $(state.f)
#             x: $(state.x)
#             v: $(state.v)
#             T: $(state.T)"
#         )
#     end
# end

# Base.show(io::IO, ::MIME"text/plain", state::State) = print(io, state)

Base.show(io::IO, result::Result) = print(io,
"Result:
    fopt: $(result.fopt)
    xopt: $(result.xopt)
    n moves: $(result.n_moves)
    n accepted: $(result.n_accepted)
    n optimums: $(result.n_optimums)")

Base.show(io::IO, options::Options) = print(io,
"Options:
    f: $(options.f)
    N: $(options.N)
    Ns: $(options.Ns)
    Nt: $(options.Nt)
    Ne: $(options.Ne)
    lb: $(options.lb)
    ub: $(options.ub)
    R: $(options.R)
    c: $(options.c)
    max: $(options.max)
    tol: $(options.tol)
    max_eval: $(options.max_eval)
    print_status: $(options.print_status)")


text:
## Load Packages
using Omega
using StatsBase
using Random
using Plots
using Distributions

Random.seed!(123)

## Define Functions
function gillespie_(rng, transitions, N, t0=0.0)

    """
    Implementation of the Gillespie algorithm for Lotka-Volterra model

    Args:
    transitions(dict): A dictionary with keys "prey" and "pred", and each element
        is a list containing the row of the transition matrix.
    N(int): Number of desired simulation iterations.
    t0(float): initial time point.
    """

    initial = Dict("prey" => prey_init(rng), "pred" => pred_init(rng))
    trajectory = Dict("prey" => [initial["prey"]], "pred" => [initial["pred"]])
    times = [t0]
    ecology = initial

    theta = Dict("spawn_prey" => spawn_prey(rng),
                 "prey2pred" => prey2pred(rng),
                 "pred_dies" => pred_dies(rng))
    t = times[1]

    for i = 1:N
        hazards = get_hazards(ecology, theta)
        # Can we replace this with a random variable?
        # It already samples one of three hazards according to some weights
        transition = transitions[sample(collect(keys(hazards)), Weights(collect(values(hazards))))]
        t = t + sum(values(hazards))

        ecology["prey"] += transition[1]
        ecology["pred"] += transition[2]
        # Enforce only positive integers
        ecology["prey"] = max(1, ecology["prey"])
        ecology["pred"] = max(1, ecology["pred"])

        append!(trajectory["prey"], ecology["prey"])
        append!(trajectory["pred"], ecology["pred"])
        append!(times, t)
    end

    return Dict(
        "prey" => trajectory["prey"],
        "pred" => trajectory["pred"],
        "times" => times
        )
end

function get_hazards(ecology, theta)
    """
    Compute the hazard function given the current states.  Note that in this
        model the function depends only on state, not time, though this is not
        the case in general. "spawn_prey" represents the event of a prey being born,
        "prey2pred" represents a predator consuming a new prey and consequently spawning
        a new predator, "pred_dies" represents the death of a predator.

    args:
        ecology(dict): a dictionary containing species counts.
        theta(dict): A dictionary where keys represent events (reactions) and
            values are the location hyperparameters for rate parameter constants
            corresponding to those events.
    """

    return Dict(
        "spawn_prey" => theta["spawn_prey"] * ecology["prey"],
        "prey2pred" => theta["prey2pred"] * ecology["prey"] * ecology["pred"],
        "pred_dies" => theta["pred_dies"] * ecology["pred"]
        )
end

# Define gillespie parameters
Pre = [[1, 0], [1, 1], [0, 1]]
Post = [[2, 0], [0, 2], [0, 0]]
transition_mat = Post - Pre
transitions = Dict("spawn_prey" => transition_mat[1,],
                    "prey2pred" => transition_mat[2,],
                    "pred_dies" => transition_mat[3,])
t0=0.0
N = 10000

# Initiate starting values randomely
prey_init = normal(100, 10)
pred_init = normal(50, 10)

# Should these be random variables? maybe not? they are in the paper
spawn_prey = normal(1, .1)
prey2pred = normal(.01, .001)
pred_dies = normal(.5, .05)

# Make gillespie into random variable
gillespie = ciid(gillespie_, transitions, 10000, 0.0)

# Sample from random variables
samples = rand((prey_init, pred_init, spawn_prey, prey2pred, pred_dies, gillespie),
                5, alg = RejectionSample)

plot_vals = hcat(values(samples[2][6]["prey"]), values(samples[2][6]["pred"]))
plot(1:10001, plot_vals[1:10001,1:2],
        title = "Basic Simulation",
        label = ["Prey" "Pred"],
        lw = 1.25)

## Sample with intervention
samples_int = rand((prey_init, pred_init, spawn_prey, prey2pred, pred_dies, gillespie),
                prey_init > 120, 5, alg = RejectionSample)

plot_vals = hcat(values(samples[2][6]["prey"]), values(samples_int[2][6]["prey"]))
plot(1:10001, plot_vals[1:10001,1:2],
        title = "Intervention Comparison",
        label = ["Prey No Intervention" "Prey > 120"],
        lw = 1.25)

plot_vals = hcat(values(samples_int[2][6]["prey"]), values(samples_int[2][6]["pred"]))
plot(1:10001, plot_vals[1:10001,1:2],
        title = "Intervention Simulation",
        label = ["Prey" "Pred"],
        lw = 1.25)

# Histogram
function intervention_hist(sim1, sim2, b)

    sim1_final_vals = zeros(0)
    sim2_final_vals = zeros(0)
    for i in 1:length(sim1)
        append!(sim1_final_vals, sim1[1][6]["prey"][i])
        append!(sim2_final_vals, sim2[1][6]["prey"][i])
    end

    final_val_diff = sim1_final_vals - sim2_final_vals
    histogram(final_val_diff, bins = b)
end

samples_100 = rand((prey_init, pred_init, spawn_prey, prey2pred, pred_dies, gillespie),
                100, alg = RejectionSample)
samples_int_100 = rand((prey_init, pred_init, spawn_prey, prey2pred, pred_dies, gillespie),
                prey_init > 120, 100, alg = RejectionSample)

intervention_hist(samples,samples_int, 20)

# Sample Rate
samples_rate_int = rand((prey_init, pred_init, spawn_prey, prey2pred, pred_dies, gillespie),
                    spawn_prey > 1.2, 5, alg = RejectionSample)

plot_vals = hcat(values(samples[2][6]["prey"]), values(samples_rate_int[2][6]["prey"]))
plot(1:5001, plot_vals[1:5001,1:2],
        title = "Intervention Comparison",
        label = ["Prey No Intervention" "Prey - spawn_prey_rate > 1.2"],
        lw = 1.25)

## Idea for intervention at t_now
function sample_weights(rng, hazards)
    result = sample(collect(keys(hazards)), Weights(collect(values(hazards))))
    return result
end

function initial_simulation(rng, t, transitions)

    ecology = Dict("prey" => prey_init(rng), "pred" => pred_init(rng))
    theta = Dict("spawn_prey" => spawn_prey(rng),
                 "prey2pred" => prey2pred(rng),
                 "pred_dies" => pred_dies(rng))

    hazards = get_hazards(ecology, theta)

    # Does this make sense?
    hazard_sample = ciid(sample_weights, hazards)
    transition = transitions[hazard_sample(rng)]
    t_new = t + sum(values(hazards))

    ecology_new_prey = ecology["prey"] + transition[1]
    ecology_new_pred = ecology["pred"] + transition[2]
    # Enforce only positive integers
    ecology_new_prey = max(1, ecology_new_prey)
    ecology_new_pred = max(1, ecology_new_pred)

    return Dict("prey" => ecology_new_prey,
                "pred" => ecology_new_pred,
                "t" => t_new)
end

temp = ciid(initial_simulation, t0, transitions)
samples = rand((prey_init, pred_init, spawn_prey, prey2pred, pred_dies, temp),
                5, alg = RejectionSample)

function get_hazards_fix(ecology, theta)
    """
    A function to fix annoying julia things
    """
    return Dict(
        "spawn_prey" => theta["spawn_prey"][1] * ecology["prey"],
        "prey2pred" => theta["prey2pred"][1] * ecology["prey"] * ecology["pred"],
        "pred_dies" => theta["pred_dies"][1] * ecology["pred"]
        )
end

function one_simulation(rng, n, theta, transitions)
    ecology = functions_list[n](rng)
    ecology_temp = copy(ecology)
    t = ecology["t"]
    delete!(ecology_temp, "t")
    hazards = get_hazards_fix(ecology_temp, theta)
    # Does this make sense?
    #hazard_sample = ciid(sample_weights, hazards)
    transition = transitions[sample(collect(keys(hazards)), Weights(collect(values(hazards))))]
    t_new = t + sum(values(hazards))

    ecology_new_prey = ecology_temp["prey"] + transition[1]
    ecology_new_pred = ecology_temp["pred"] + transition[2]
    # Enforce only positive integers
    ecology_new_prey = max(1, ecology_new_prey)
    ecology_new_pred = max(1, ecology_new_pred)

    return Dict("prey" => ecology_new_prey,
                "pred" => ecology_new_pred,
                "t" => t_new)
end

temp = ciid(initial_simulation, t0, transitions)
functions_list = Any[]
push!(functions_list, temp)
N = 5
for f in 2:N
    last = f - 1
    temp = ciid(one_simulation, last, theta, transitions)
    push!(functions_list, temp)
end

push!(functions_list,prey_init)
push!(functions_list, pred_init)
push!(functions_list, spawn_prey)
push!(functions_list, prey2pred)
push!(functions_list, pred_dies)

samples = rand((prey_init, pred_init, spawn_prey, prey2pred, pred_dies,
functions_list[1], functions_list[2]),
                5, alg = RejectionSample)

samples = rand(Tuple(x for x in functions_list),
                5, alg = RejectionSample)


text:
export Rules, Lexicon, Grammar,
        rewrites_for, is_category, cnf_rules, rewrite, generate_random_sentence,
        ProbabilityRules, ProbabilityLexicon, ProbabilityGrammar,
        Chart, parse_sentence, cyk_parse,
        Page, load_page_html, determine_inlinks, find_outlinks, init_pages, only_wikipedia_urls,
        expand_pages, relevant_pages, normalize_pages,
        ConvergenceDetector, detect_convergence, get_inlinks, get_outlinks, HITS;

"""
    Rules(rules_array::Array{T}) where {T <: Pair}

Return a Dict of mappings for symbols (lexical categories) to alternative sequences.
"""
function Rules(rules_array::Array{T, 1}) where {T <: Pair}
    local rules::Dict = Dict();
    for (lhs, rhs) in rules_array
        rules[lhs] = collect(map(String, split(strip(ss))) for ss in map(String, split(rhs, ['|'])));
    end
    return rules;
end

"""
    Lexicon(rules_array::Array{T}) where {T <: Pair}

Return a Dict of mappings for symbols (lexical categories) to alternative words.

The lexicon is the list of allowable words.
"""
function Lexicon(rules_array::Array{T, 1}) where {T <: Pair}
    local rules::Dict = Dict();
    for (lhs, rhs) in rules_array
        rules[lhs] = collect(strip(ss) for ss in map(String, split(rhs, "|")));
    end
    return rules;
end

#=

    Grammar consists of a set of rules and a lexicon.

=#
struct Grammar
    name::String
    rules::Dict
    lexicon::Dict
    categories::Dict

    function Grammar(name::String, rules::Dict, lexicon::Dict)
        local ng::Grammar = new(name, rules, lexicon, Dict());
        for (category, words) in ng.lexicon
            for word in words
                ng.categories[word] = push!(get!(ng.categories, word, Array{String, 1}()), category);
            end
        end
        return ng;
    end
end

"""
    rewrites_for(g::Grammar, cat::String)

Return an Array of possible rhs's that category 'cat' can be rewritten as.
"""
function rewrites_for(g::Grammar, cat::String)
    return get(g.rules, cat, Array{String, 1}());
end

"""
    is_category(g::Grammar, word::String, cat::String)

Return whether the given word 'word' is of category 'cat'.
"""
function is_category(g::Grammar, word::String, cat::String)
    return (cat in g.categories[word]);
end

"""
    cnf_rules(g::Grammar)

Return the rules of grammar 'g' as an Array of Tuples (X, Y, Z) such that X -> Y Z.
"""
function cnf_rules(g::Grammar)
    local cnf::AbstractVector = [];
    for (x, rules) in g.rules
        for (y, z) in rules
            push!(cnf, (x, y, z));
        end
    end
    return cnf;
end

"""
    rewrite(g::Grammar, tokens::AbstractVector, into::AbstractVector)

Return the resulting array of words by replacing each token in 'tokens' with a random word.
"""
function rewrite(g::Grammar, tokens::AbstractVector, into::AbstractVector)
    for token in tokens
        if (token in keys(g.rules))
            rewrite(g, rand(RandomDeviceInstance, g.rules[token]), into);
        elseif (token in keys(g.lexicon))
            push!(into, rand(RandomDeviceInstance, g.lexicon[token]));
        else
            push!(into, token);
        end
    end
    return into;
end

"""
    generate_random_sentence(g::Grammar, categories::String)
    generate_random_sentence(g::Grammar)

Return a randomly generated sentence as a String by using the given categories 'categories'.
"""
function generate_random_sentence(g::Grammar, categories::String)
    return join(rewrite(g, split(categories), []), " ");
end

function generate_random_sentence(g::Grammar)
    return generate_random_sentence(g, "S");
end

"""
    ProbabilityRules(rules_array::Array{T}) where {T <: Pair}

Return a Dict of mappings for symbols (lexical categories) to alternative sequences with probabilities.
"""
function ProbabilityRules(rules_array::Array{T, 1}) where {T <: Pair}
    local rules::Dict = Dict();
    for (lhs, rhs) in rules_array
        rules[lhs] = [];
        local rhs_split::AbstractVector = collect(map(String, split(strip(ss))) for ss in map(String, split(rhs, ['|'])));
        for rule in rhs_split
            local rule_probability::Float64 = parse(Float64, rule[end][2:(end - 1)]);
            local rule_tuple::Tuple = (rule[1:(end - 1)], rule_probability);
            push!(rules[lhs], rule_tuple);
        end
    end
    return rules;
end

"""
    ProbabilityLexicon(rules_array::Array{T}) where {T <: Pair}

Return a Dict of mappings for symbols (lexical categories) to alternative words with probabilities.

The lexicon is the list of allowable words.
"""
function ProbabilityLexicon(rules_array::Array{T, 1}) where {T <: Pair}
    local rules::Dict = Dict();
    for (lhs, rhs) in rules_array
        rules[lhs] = [];
        local rhs_split::AbstractVector = collect(map(String, split(strip(word))) for word in map(String, split(rhs, ['|'])));
        for rule in rhs_split
            local word_probability::Float64 = parse(Float64, rule[end][2:(end - 1)]);
            local rule_word::String = rule[1];
            local rule_tuple::Tuple = (rule_word, word_probability);
            push!(rules[lhs], rule_tuple);
        end
    end
    return rules;
end

#=

    ProbabilityGrammar consists of a set of rules and a lexicon.

=#
struct ProbabilityGrammar
    name::String
    rules::Dict
    lexicon::Dict
    categories::Dict

    function ProbabilityGrammar(name::String, rules::Dict, lexicon::Dict)
        local npg::ProbabilityGrammar = new(name, rules, lexicon, Dict());
        for (category, words) in npg.lexicon
            for (word, p) in words
                npg.categories[word] = push!(get!(npg.categories, word, []), (category, p));
            end
        end
        return npg;
    end
end

"""
    rewrites_for(pg::ProbabilityGrammar, cat::String)

Return an Array of possible rhs's that category 'cat' can be rewritten as.
"""
function rewrites_for(pg::ProbabilityGrammar, cat::String)
    return get(pg.rules, cat, []);
end

"""
    is_category(pg::ProbabilityGrammar, word::String, cat::String)

Return whether the given word 'word' is of category 'cat'.
"""
function is_category(pg::ProbabilityGrammar, word::String, cat::String)
    return (cat in map(first, pg.categories[word]));
end

"""
    cnf_rules(pg::ProbabilityGrammar)

Return the rules of grammar 'pg' as an Array of Tuples (X, Y, Z, p) such that X -> Y Z [p].
"""
function cnf_rules(pg::ProbabilityGrammar)
    local cnf::AbstractVector = [];
    for (x, rules) in pg.rules
        for ((y, z), p) in rules
            push!(cnf, (x, y, z, p));
        end
    end
    return cnf;
end

"""
    rewrite(pg::ProbabilityGrammar, tokens::AbstractVector, into::AbstractVector)

Return the resulting array of words by replacing each token in 'tokens' with a random word.
"""
function rewrite(pg::ProbabilityGrammar, tokens::AbstractVector, into::AbstractVector)
    local p::Float64;

    for token in tokens
        if (token in keys(pg.rules))
            local non_terminal_symbols::AbstractVector;
            non_terminal_symbols, p = weighted_choice(pg.rules[token]);
            into[2] = into[2] * p;
            rewrite(pg, non_terminal_symbols, into);
        elseif (token in keys(pg.lexicon))
            local terminal_symbol::String;
            terminal_symbol, p = weighted_choice(pg.lexicon[token]);
            push!(into[1], terminal_symbol);
            into[2] = into[2] * p;
        else
            push!(into[1], token);
        end
    end
    return into;
end

"""
    generate_random_sentence(pg::ProbabilityGrammar, categories::String)
    generate_random_sentence(pg::ProbabilityGrammar)

Return a randomly generated sentence as a String by using the given categories 'categories'.
"""
function generate_random_sentence(pg::ProbabilityGrammar, categories::String)
    local rewritten_sentence::AbstractVector;
    local p::Float64;
    rewritten_sentence, p = rewrite(pg, split(categories), [[], 1.0]);
    return (join(rewritten_sentence, " "), p);
end

function generate_random_sentence(pg::ProbabilityGrammar)
    return generate_random_sentence(pg, "S");
end

# Create a submodule to contain the following example grammars.
# These grammars are defined in the module path "aimajulia.nlp".

module nlp

using ..aimajulia;

epsilon_0 = Grammar("ε_0",
                    Rules(      # Grammar for ε_0 (Fig. 22.4 2nd edition)
                        ["S"=>"NP VP | S Conjunction S",
                        "NP"=>"Pronoun | Name | Noun | Article Noun | Digit Digit | NP PP | NP RelClause",
                        "VP"=>"Verb | VP NP | VP Adjective | VP PP | VP Adverb",
                        "PP"=>"Preposition NP",
                        "RelClause"=>"That VP"]),
                    Lexicon(    # Lexicon for ε_0 (Fig. 22.3 2nd edition)
                        ["Noun"=>"stench | breeze | glitter | nothing | wumpus | pit | pits | gold | east",
                        "Verb"=>"is | see | smell | shoot | fell | stinks | go | grab | carry | kill | turn | feel",
                        "Adjective"=>"right | left | east | south | back | smelly",
                        "Adverb"=>"here | there | nearby | ahead | right | left | east | south | back",
                        "Pronoun"=>"me | you | I | it",
                        "Name"=>"John | Mary | Boston | Aristotle",
                        "Article"=>"the | a | an",
                        "Preposition"=>"to | in | on | near",
                        "Conjunction"=>"and | or | but",
                        "Digit"=>"0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9"
                        ]));

# An example grammar for testing
epsilon_ = Grammar("ε_",
                    Rules(["S"=>"NP VP",
                            "NP"=>"Art N | Pronoun",
                            "VP"=>"V NP"]),
                    Lexicon(["Art"=>"the | a",
                            "N"=>"man | woman | table | shoelace | saw",
                            "Pronoun"=>"I | you | it",
                            "V"=>"saw | liked | feel"]));

# An example grammar for testing
epsilon_np_ = Grammar("ε_NP_",
                    Rules(["NP"=>"Adj NP | N"]),
                    Lexicon(["Adj"=>"happy | handsome | hairy", "N"=>"man"]));

# ε_probability is a probabilistic grammar found in the Python notebook.
# The rules (Fig 23.2 3rd edition) for the language 'ε_probability' use the probabilities from the Python notebook.
epsilon_probability = ProbabilityGrammar("ε_probability",
                                        ProbabilityRules(["S"=>"NP VP [0.6] | S Conjunction S [0.4]",
                                                        "NP"=>("Pronoun [0.2] | Name [0.05] | Noun [0.2] | "*
                                                                "Article Noun [0.15] | Article Adjs Noun [0.1] | "*
                                                                "Digit [0.05] | NP PP [0.15] | NP RelClause [0.1]"),
                                                        "VP"=>("Verb [0.3] | VP NP [0.2] | VP Adjective [0.25] | "*
                                                                "VP PP [0.15] | VP Adverb [0.1]"),
                                                        "Adjs"=>"Adjective [0.5] | Adjective Adjs [0.5]",
                                                        "PP"=>"Preposition NP [1]",
                                                        "RelClause"=>"RelPro VP [1]"]),
                                        ProbabilityLexicon(["Verb"=>"is [0.5] | say [0.3] | are [0.2]",
                                                            "Noun"=>"robot [0.4] | sheep [0.4] | fence [0.2]",
                                                            "Adjective"=>"good [0.5] | new [0.2] | sad [0.3]",
                                                            "Adverb"=>"here [0.6] | lightly [0.1] | now [0.3]",
                                                            "Pronoun"=>"me [0.3] | you [0.4] | he [0.3]",
                                                            "RelPro"=>"that [0.5] | who [0.3] | which [0.2]",
                                                            "Name"=>"john [0.4] | mary [0.4] | peter [0.2]",
                                                            "Article"=>"the [0.5] | a [0.25] | an [0.25]",
                                                            "Preposition"=>"to [0.4] | in [0.3] | at [0.3]",
                                                            "Conjuction"=>"and [0.5] | or [0.2] | but [0.3]",
                                                            "Digit"=>"0 [0.35] | 1 [0.35] | 2 [0.3]"]));

epsilon_chomsky = Grammar("ε_chomsky",
                        Rules(["S"=>"NP VP",
                                "NP"=>"Article Noun | Adjective Noun",
                                "VP"=>"Verb NP | Verb Adjective"]),
                        Lexicon(["Article"=>"the | a | an",
                                "Noun"=>"robot | sheep | fence",
                                "Adjective"=>"good | new | sad",
                                "Verb"=>"is | say | are"]));

epsilon_probability_chomsky = ProbabilityGrammar("ε_probability_chomsky",
                                                ProbabilityRules(["S"=>"NP VP [1]",
                                                                "NP"=>"Article Noun [0.6] | Adjective Noun [0.4]",
                                                                "VP"=>"Verb NP [0.5] | Verb Adjective [0.5]"]),
                                                ProbabilityLexicon(["Article"=>"the [0.5] | a [0.25] | an [0.25]",
                                                                    "Noun"=>"robot [0.4] | sheep [0.4] | fence [0.2]",
                                                                    "Adjective"=>"good [0.5] | new [0.2] | sad [0.3]",
                                                                    "Verb"=>"is [0.5] | say [0.3] | are [0.2]"]));

end

#=

    Chart is a data structure used in the process of analyzing a string

    of words to uncover the phrase structure.

=#
mutable struct Chart
    trace::Bool
    grammar::Grammar
    chart::AbstractVector

    function Chart(g::Grammar; trace::Bool= false)
        return new(trace, g);
    end
end

"""
    parse_sentence(chart::Chart, words::String, categories::String)
    parse_sentence(chart::Chart, words::String)
    parse_sentence(chart::Chart, words::Array{String, 1}, categories::String)
    parse_sentence(chart::Chart, words::Array{String, 1})

Return an array of parses given the sentence 'words' and categories 'categories'.
"""
function parse_sentence(chart::Chart, words::String, categories::String)
    local words_array::Array{String, 1} = map(String, split(words));
    return parse_sentence(chart, words_array, categories);
end

function parse_sentence(chart::Chart, words::String)
    local words_array::Array{String, 1} = map(String, split(words));
    return parse_sentence(chart, words_array, "S");
end

function parse_sentence(chart::Chart, words::Array{String, 1}, categories::String)
    parse_words(chart, words, categories);
    return collect([i, j, categories, found, []]
                    for (i, j, lhs, found, expects) in chart.chart[length(words)]);
end

function parse_sentence(chart::Chart, words::Array{String, 1})
    return parse_sentence(chart, words, "S");
end

"""
    parse_words(chart::Chart, words::Array{String, 1}, categories::String)
    parse_words(chart::Chart, words::Array{String, 1})

Return an array of words given the array of 'words' and categories 'categories'.
"""
function parse_words(chart::Chart, words::Array{String, 1}, categories::String)
    chart.chart = collect([] for i in 1:(length(words) + 1));
    add_edge(chart, [1, 1, "S_", [], [categories]]);
    for i in 1:length(words)
        scanner(chart, i, words[i]);
    end
    return chart.chart;
end

function parse_words(chart::Chart, words::Array{String, 1})
    return parse_words(chart, words, "S");
end

"""
    add_edge(chart::Chart, edge::AbstractVector)

Add the given edge 'edge' to the chart 'chart'.
"""
function add_edge(chart::Chart, edge::AbstractVector)
    local start_index::Int64;
    local end_index::Int64;
    local lhs::String;
    local found::AbstractVector;
    local expects::AbstractVector;

    start_index, end_index, lhs, found, expects = edge;
    if (!(edge in chart.chart[end_index]))
        push!(chart.chart[end_index], edge);
        if (chart.trace)
            println("Chart: added ", edge);
        end
        if (length(expects) == 0)
            extender(chart, edge);
        else
            predictor(chart, edge);
        end
    end
    return nothing;
end

"""
    scanner(chart::Chart, index::Int64, word::String)

Extend the edge for given index 'index' if the given word 'word' and its category are expected.
"""
function scanner(chart::Chart, index::Int64, word::String)
    for (i, j, A, alpha, Bb) in chart.chart[index]
        if ((length(Bb) != 0) && (is_category(chart.grammar, word, Bb[1])))
            add_edge(chart, [i, (j + 1), A, vcat(alpha, [(Bb[1], word)]), Bb[2:end]]);
        end
    end
    return nothing;
end

"""
    predictor(chart::Chart, edge::AbstractVector)

Add edges with rules for 'B' that may help in extending the given edge 'edge'.
"""
function predictor(chart::Chart, edge::AbstractVector)
    local i::Int64;
    local j::Int64;
    local A::String;
    local alpha::AbstractVector;
    local Bb::AbstractVector;

    i, j, A, alpha, Bb = edge;
    local B::String = Bb[1];
    if (haskey(chart.grammar.rules, B))
        for rhs in rewrites_for(chart.grammar, B)
            add_edge(chart, [j, j, B, [], rhs]);
        end
    end
    return nothing;
end

"""
    extender(chart::Chart, edge::AbstractVector)

Extend whatever edge can be extended from the given edge 'edge' to the chart 'chart'.
"""
function extender(chart::Chart, edge::AbstractVector)
    local m::Int64;
    local n::Int64;
    local B::String;
    local found::AbstractVector;
    local expects::AbstractVector;

    m, n, B, found, expects = edge;
    for (i, j, A, alpha, B1b) in chart.chart[m]
        if ((length(B1b) != 0) && (B == B1b[1]))
            add_edge(chart, [i, n, A, vcat(alpha, [edge]), B1b[2:end]]);
        end
    end
    return nothing;
end

"""
    cyk_parse(words::Array{String, 1}, grammar::ProbabilityGrammar)

Return the resulting table as a Dict by applying the CYK algorithm (Fig. 23.5) on the
given sequence of words 'words' and grammar 'grammar'.
"""
function cyk_parse(words::Array{String, 1}, grammar::ProbabilityGrammar)
    local N::Int64 = length(words);
    local P::Dict = Dict();

    for (i, word) in enumerate(words)
        for (X, p) in get!(grammar.categories, word, [])
            P[(X, i, 1)] = p;
        end
    end

    for l in 2:N
        for start_index in 1:(N - l + 1)
            for len1 in 1:(l - 1)
                local len2::Int64 = l - len1;
                for (X, Y, Z, p) in cnf_rules(grammar)
                    P[(X, start_index, l)] = max(get!(P, (X, start_index, l), 0.0),
                                                (get!(P, (Y, start_index, len1), 0.0)
                                                * get!(P, (Z, start_index + len1, len2), 0.0) * p));
                end
            end
        end
    end
    return P;
end

#=

    Page consists of an address, an array of inlinks, an array of outlinks, a hub score, and an authority score.

=#
mutable struct Page
    address::String
    inlinks::AbstractVector
    outlinks::AbstractVector
    hub::Float64
    authority::Float64

    function Page(address::String; inlinks::AbstractVector=[], outlinks::AbstractVector=[], hub::Int64=0, authority::Int64=0)
        return new(address, inlinks, outlinks, hub, authority);
    end
end

"""
    load_page_html(addresses::AbstractVector)

Return a Dict of page content for the given URLs 'addresses' by downloading the HTML pages.
"""
function load_page_html(addresses::AbstractVector)
    local content::Dict = Dict();
    for address in addresses
        local tmp_file::Tuple = mktemp(pwd());
        close(tmp_file[2]);
        rm(tmp_file[1]);
        download(address, tmp_file[1]);
        local raw_html::String = String(read(tmp_file[1]));
        rm(tmp_file[1]);
        local html::String = replace(raw_html, @r_str("<head>(.*)</head>", "s")=>"");
        content[address] = html;
    end
    return content;
end

"""
    determine_inlinks(page::Page, pages_index::Dict)

Return an Array of inlinks for the given page 'page'.
"""
function determine_inlinks(page::Page, pages_index::Dict)
    local inlinks::AbstractVector = [];
    for (address, index_page) in pages_index
        if (page.address == index_page.address)
            continue;
        elseif (page.address in index_page.outlinks)
            push!(inlinks, address);
        end
    end
    return inlinks;
end

"""
    find_outlinks(page::Page, pages_content::Dict)
    find_outlinks(page::Page, pages_content::Dict, handle_urls::Function)

Return an Array of outlinks to other pages for the given page 'page'.

If the argument 'handle_urls' is given, the resulting array is returned
after the function is applied to the array of outlinks.
"""
function find_outlinks(page::Page, pages_content::Dict)
    local urls::AbstractVector = Array{String, 1}();
    for regex_m in eachmatch(r"href=['\"]?([^'\" >]+)", pages_content[page.address])
        push!(urls, String(regex_m.captures[1]));
    end
    return urls;
end

function find_outlinks(page::Page, pages_content::Dict, handle_urls::Function)
    local urls::AbstractVector = Array{String, 1}();
    for regex_m in eachmatch(r"href=['\"]?([^'\" >]+)", pages_content[page.address])
        push!(urls, String(regex_m.captures[1]));
    end
    return handle_urls(urls);
end

"""
    init_pages(addresses::AbstractVector)

Return a Dict of pages from the given array of URLs 'addresses'.
"""
function init_pages(addresses::AbstractVector)
    local pages::Dict = Dict();
    for address in addresses
        pages[address] = Page(address);
    end
    return pages;
end

"""
    only_wikipedia_urls(urls::AbstractVector)

Return an Array of wikipedia URLs where relative wiki URLs are converted to their absolute URL.
"""
function only_wikipedia_urls(urls::AbstractVector)
    local wiki_urls::AbstractVector = collect(url for url in urls if (startswith(url, "/wiki/")));
    return collect(("https://en.wikipedia.org" * url) for url in wiki_urls);
end

"""
    expand_pages(pages::Dict, pages_index::Dict)

Return the Dict of expanded pages by adding in every page that links to
or is linked from one of the relevant pages.
"""
function expand_pages(pages::Dict, pages_index::Dict)
    local expanded::Dict = Dict();
    for (address, page) in pages
        if (!haskey(expanded, address))
            expanded[address] = page;
        end
        for inlink in page.inlinks
            if (!haskey(expanded, inlink))
                expanded[inlink] = pages_index[inlink];
            end
        end
        for outlink in page.outlinks
            if (!haskey(expanded, outlink))
                expanded[outlink] = pages_index[outlink];
            end
        end
    end
    return expanded;
end

"""
    relevant_pages(query::String, pages_index::Dict, pages_content::Dict)

Return a Dict of pages that contain all of the query words in 'query'.
These pages are found by intersecting the hit lists of the query words.
"""
function relevant_pages(query::String, pages_index::Dict, pages_content::Dict)
    local hit_intersection::Set = Set(collect(keys(pages_index)));
    local query_words::AbstractVector = map(String, split(query));
    for query_word in query_words
        local hit_list::Set = Set();
        for address in keys(pages_index)
            if (occursin(lowercase(query_word), lowercase(pages_content[address])))
                push!(hit_list, address);
            end
        end
        hit_intersection = intersect(hit_intersection, hit_list);
    end
    return Dict(collect((address, pages_index[address]) for address in hit_intersection));
end

"""
    normalize_pages(pages::Dict)

Divide the scores of each page in 'pages' by the sum of all the squares of all pages' scores
(separately for both the authority and hubs scores).
"""
function normalize_pages(pages::Dict)
    local summed_hub::Float64 = sum(collect(page.hub^2 for page in values(pages)));
    local summed_authority::Float64 = sum(collect(page.authority^2 for page in values(pages)));
    local sqrt_summed_hub::Float64 = sqrt(summed_hub);
    local sqrt_summed_authority::Float64 = sqrt(summed_authority);
    for address in keys(pages)
        pages[address].hub = pages[address].hub / sqrt_summed_hub;
        pages[address].authority = pages[address].authority / sqrt_summed_authority;
    end
    return nothing;
end

#=

    ConvergenceDetector contains the hub history and authorities history for a set of pages.

=#
struct ConvergenceDetector
    hub_history::AbstractVector
    authority_history::AbstractVector

    function ConvergenceDetector()
        return new([], []);
    end
end

"""
    detect_convergence(cd::ConvergenceDetector, pages_index::Dict)

Return a boolean indicating if both the 'hub' and 'authority' values of the pages
in 'pages_index' have converged.
"""
function detect_convergence(cd::ConvergenceDetector, pages_index::Dict)
    local current_hubs::AbstractVector = collect(page.hub for page in values(pages_index));
    local current_authorities::AbstractVector = collect(page.authority for page in values(pages_index));
    if (length(cd.hub_history) != 0)
        local diffs_hub::AbstractVector = collect(abs(x - y) for (x, y) in zip(current_hubs, cd.hub_history[end]));
        local diffs_authority::AbstractVector = collect(abs(x - y) for (x, y) in zip(current_authorities, cd.authority_history[end]));
        local average_delta_hub::Float64 = sum(diffs_hub)/length(pages_index);
        local average_delta_authority::Float64 = sum(diffs_authority)/length(pages_index);
        if ((average_delta_hub < 0.01) && (average_delta_authority < 0.01))
            return true;
        end
    end
    if (length(cd.hub_history) > 2)
        deleteat!(cd.hub_history, 1);
        deleteat!(cd.authority_history, 1);
    end
    push!(cd.hub_history, current_hubs);
    push!(cd.authority_history, current_authorities);
    return false;
end

"""
    get_inlinks(page::Page, pages_index::Dict)

Return an Array of addresses where each address is in both the page's inlinks
and index of pages 'pages_index'.
"""
function get_inlinks(page::Page, pages_index::Dict)
    if (length(page.inlinks) == 0)
        page.inlinks = determine_inlinks(page, pages_index);
    end
    return collect(address for (address, p) in pages_index if (address in page.inlinks));
end

"""
    get_outlinks(page::Page, pages_index::Dict, pages_content::Dict)

Return an Array of addresses where each address is in both the page's outlinks
and index of pages 'pages_index'.
"""
function get_outlinks(page::Page, pages_index::Dict, pages_content::Dict)
    if (length(page.outlinks) == 0)
        page.outlinks = find_outlinks(page, pages_content);
    end
    return collect(address for (address, p) in pages_index if (address in page.outlinks));
end

"""
    HITS(query::String, pages_index::Dict, pages_content::Dict)

Return the computed hubs and authorities with respect to the query as a Dict by using the
HITS algorithm (Fig. 22.1) to the given query 'query', index of pages 'pages_index',
and the content of the pages 'pages_content'.
"""
function HITS(query::String, pages_index::Dict, pages_content::Dict)
    local pages::Dict = expand_pages(relevant_pages(query, pages_index, pages_content), pages_index);
    for p in values(pages)
        p.authority = 1.0;
        p.hub = 1.0;
    end
    local convergence::ConvergenceDetector = ConvergenceDetector();
    while (detect_convergence(convergence, pages_index))
        local authority::Dict = Dict(collect(Pair(p, pages[p].authority) for p in pages));
        local hub::Dict = Dict(collect(Pair(p, pages[p].hub) for p in pages));
        for p in values(pages)
            p.authority = sum(hub[x] for x in get_inlinks(p, pages_index));
            p.hub = sum(authority[x] for x in get_outlinks(p, pages_index, pages_content));
        end
        normalize_pages(pages);
    end
    return pages;
end

max_stars_repo_name: mikhail-j/aima-julia-aimacode-fork

max_issues_repo_name: mikhail-j/aima-julia-aimacode-fork

max_forks_repo_name: mikhail-j/aima-julia-aimacode-fork


text:
  function sample_space!(sp,domain,nsvt,nframes)
    for i in 1:nsvt
      dt = 0
      for t in 1:nframes
        for j in 1:nframes
          if ((j + dt) <= nframes )
            if domain[j,i] == 1
              sp[t,i] = sp[t,i] + 1
            end
          end
        end
        dt = dt + 1
      end
    end
  end


text:
## waveform and times
ramp =  [  # -0.0400066666667    0.5
			-0.0400000000000    0.0
			-0.0399933333333   -0.5
			-0.0200066666667   -0.5
			-0.0200000000000    0.0
			-0.0199933333333    0.5
			-0.0000066666667    0.5
			 0.0000000000000    0.0
			 0.0000066666667   -0.5]

times = vec(10 .^mean(log10.([
            0.0000066667	0.0000200000
			0.0000333333	0.0000466667
			0.0000600000	0.0000733333
			0.0000866667	0.0001266667
			0.0001400000	0.0002066667
			0.0002200000	0.0003400000
			0.0003533333	0.0005533333
			0.0005666667	0.0008733333
			0.0008866667	0.0013533333
			0.0013666667	0.0021000000
			0.0021133333	0.0032733333
			0.0032866667	0.0051133333
			0.0051266667	0.0079933333
			0.0080066667	0.0123933333
			0.0124066667	0.0199933333]), dims=2))

## High alt noise
Hx_add_noise = [0.010619 0.009453 0.008506 0.006687 0.007244 0.005554 0.004701 0.004353 0.003539 0.003493 0.003035 0.002875 0.002343 0.001613 0.001304][:]
Hz_add_noise = [0.005554 0.005280 0.004101 0.003093 0.002969 0.002723 0.002696 0.002429 0.002377 0.002188 0.002018 0.001818 0.001557 0.001106 0.000906][:]

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
# =========================================================================== #
# Types
# =========================================================================== #

abstract type AbstractCurve{T}
end

struct StaticCurve{T} <: AbstractCurve{T}
    x::Vector{T}
    y::Vector{T}
end
function StaticCurve(x::Vector{T}, y::Vector{T}) where {T}
    return StaticCurve{T}(copy(x), copy(y))
end

mutable struct MutableCurve{T} <: AbstractCurve{T}
    x::Vector{T}
    y::Vector{T}
end
function MutableCurve(x::Vector{T}, y::Vector{T}) where {T}
    return MutableCurve{T}(x, y)
end

# =========================================================================== #
# Base functions
# =========================================================================== #

function Base.eltype(curve::T) where {T <: AbstractCurve}
    return eltype(curve.x)
end

function Base.length(curve::AbstractCurve)
    return length(curve.x)
end

function Base.size(curve::AbstractCurve)
    return (length(curve.x), 2)
end

function Base.show(io::IO, curve::T) where {T <: AbstractCurve}
    println("$(T) - curve with $(length(curve)) elements:")
    println("\tx\ty")
    for kk = 1:length(curve)
        println("\t$(curve.x[kk])\t$(curve.y[kk])")
    end
    return nothing
end

function Base.push!(curve::AbstractCurve, x, y)
    push!(curve.x, x)
    push!(curve.y, y)
    return curve
end

# =========================================================================== #
# Math functions
# =========================================================================== #

function add_x!(curve::AbstractCurve, factor)
    curve.x .+= factor
    return curve
end
function add_y!(curve::AbstractCurve, factor)
    curve.y .+= factor
    return curve
end

function scale_x!(curve::AbstractCurve, factor)
    curve.x .*= factor
    return curve
end
function scale_y!(curve::AbstractCurve, factor)
    curve.y .*= factor
    return curve
end

function map_x!(curve::AbstractCurve, fcn, args...)
    for kk = 1:length(curve)
        curve.x[kk] = fcn(curve.x[kk], args...)
    end
    return curve
end

function map_y!(curve::AbstractCurve, fcn, args...)
    for kk = 1:length(curve)
        curve.y[kk] = fcn(curve.y[kk], args...)
    end
    return curve
end

function map!(curve::AbstractCurve, fcn, args...)
    fcn(curve, args...)
    return curve
end

# =========================================================================== #
# Special creation functions
# =========================================================================== #

function ramp_curve(T, x_start, x_delta, x_end)
    x = ramp(x_start, x_delta, x_end)
    y = copy(x)
    return T(x, y)
end

function linspace_curve(T, x_start, x_end, x_length)
    x = linspace( x_start, x_end, x_length )
    y = copy(x)
    return T(x, y)
end

function rand_curve(T, x_length)
    x = collect(1:x_length)
    y = rand(x_length)
    return T(x, y)
end

function Base.copy(curve::T) where {T <: AbstractCurve}
    new_curve = T( copy(curve.x), copy(curve.y) )
    return new_curve
end

function inverse(curve::T) where {T <: AbstractCurve}
    new_curve = T(curve.y, curve.x)
    return new_curve
end

## some issue with the implementation
# function initialize(::Type{T}, n::Int) where {T <: AbstractCurve}
#     return T( Vector{Float64}(undef, n), Vector{Float64}(undef, n) )
# end

# function Base.similar(curve::T) where {T <: AbstractCurve}
#     new_curve = T( copy(curve.x), copy(curve.y) )
# end

# =========================================================================== #
# Interpolations
# =========================================================================== #

function linear_search_interpolate(curve::AbstractCurve{T}, xq::T; init_index::Int = 2) where {T}
    # linear search
    # it is assumed that x is strictly monotonically increasing
    # only the first match is found
    #
    # constant extrapolation
    # linear interpolation

    L = length(curve)

    # handle extrapolation
    if xq <= curve.x[1]
        return (curve.y[1], 1)
    elseif xq >= curve.x[L]
        return (curve.y[L], L)
    end

    # handle interpolation
    ii = init_index
    yq = NaN
    left_idx = 0
    while ii <= L
        if xq <= curve.x[ii]
            left_idx = ii-1
            break
        end
        ii += 1
    end

    # interpolation
    yq = (curve.y[ii] - curve.y[ii-1]) / (curve.x[ii] - curve.x[ii-1]) * (xq - curve.x[ii-1]) + curve.y[ii-1]

    return (yq, left_idx)

end

function binary_search_interpolate(curve::AbstractCurve, xq; left_init_index = 1, right_init_index = length(curve), linear_search_crit_index = 0 )
    # binary search
    # it is assumed that x is strictly monotonically increasing
    #
    # constant extrapolation
    # linear interpolation

    L = length(curve)

    # handle extrapolation
    if xq <= curve.x[1]
        return (curve.y[1], 1)
    elseif xq >= curve.x[L]
        return (curve.y[L], L)
    end

    # handle interpolation
    left_idx  = left_init_index
    right_idx = right_init_index
    mid_idx   = div(right_idx + left_idx, 2)

    ii = 0
    yq = NaN
    while true

        if xq <= curve.x[mid_idx]
            right_idx = mid_idx
        else # xq > curve.x[mid_idx]
            left_idx = mid_idx
        end
        d_idx = right_idx - left_idx

        if d_idx <= linear_search_crit_index
            (yq, right_idx) = linear_search_interpolate(curve, xq; init_index = left_idx)
            break
        end

        mid_idx = div(right_idx + left_idx, 2)
        if mid_idx == left_idx
            # right_idx = left_idx + 1
            ii = left_idx + 1
            break
        end
    end

    # interpolation
    yq = (curve.y[ii] - curve.y[ii-1]) / (curve.x[ii] - curve.x[ii-1]) * (xq - curve.x[ii-1]) + curve.y[ii-1]

    return (yq, left_idx)

end

function vector_linear_search_interpolate(crv::AbstractCurve{T}, x_vec::Vector{T}; init_idx::Int = 2) where {T}
    x_query = sort(x_vec)
    left_idx = init_idx
    y_vec = similar(x_vec)
    for kk = 1:length(x_query)
        xq = x_query[kk]
        (y_vec[kk], left_idx) = linear_search_interpolate(crv, xq, init_index = left_idx)
    end
    return y_vec
end

function vector_binary_search_interpolate(crv::AbstractCurve{T}, x_vec::Vector{T}; init_idx::Int = 1) where {T}
    x_query = sort(x_vec)
    left_idx = init_idx
    y_vec = similar(x_vec)
    for kk = 1:length(x_query)
        xq = x_query[kk]
        (y_vec[kk], left_idx) = binary_search_interpolate(crv, xq, left_init_index = left_idx)
    end
    return y_vec
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using BinaryBuilder

name = "VMEC"
version = v"1.0.0"

sources = [
    ArchiveSource("https://gitlab.com/wistell/VMEC2000/-/archive/v1.0.0.tar", "179761d5f02cb530fbf6da4121c394a1d585c49c24391ebd87fe5c014e0addcf"),
]

script = raw"""
mkdir -p ${libdir}
cd ${WORKSPACE}/srcdir/VMEC2000*
make USE_MKL=1
install -m755 lib/libvmec_mkl.so ${libdir}
make USE_MKL=1 vmec_clean ; make USE_MKL=0
install -m755 lib/libvmec_openblas.so ${libdir}
make USE_MKL=0 vmec_clean
"""

platforms = [
    Platform("x86_64", "linux"; libc = "glibc", libgfortran_version="4.0.0"),
    Platform("x86_64", "linux"; libc = "glibc", libgfortran_version="5.0.0"),
]
#platforms = expand_gfortran_versions(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libvmec_mkl", :libvmec_mkl, dont_dlopen = true),
    LibraryProduct("libvmec_openblas", :libvmec_openblas, dont_dlopen = true),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("MPICH_jll"),
    Dependency("OpenBLAS_jll", v"0.3.17"),
    Dependency("SCALAPACK_jll"),
    Dependency("MKL_jll"),
    Dependency("oneTBB_jll"),
    Dependency("NetCDF_jll"),
    Dependency("NetCDFF_jll"),
    Dependency("CompilerSupportLibraries_jll")
]

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

max_stars_repo_name: AbhimanyuAryan/Yggdrasil

max_issues_repo_name: AbhimanyuAryan/Yggdrasil

max_forks_repo_name: AbhimanyuAryan/Yggdrasil


text:
using Base.BinaryPlatforms
function augment_platform!(platform)
    llvm_version = Base.libllvm_version

    # does our LLVM build use assertions?
    llvm_assertions = try
        cglobal((:_ZN4llvm24DisableABIBreakingChecksE, Base.libllvm_path()), Cvoid)
        false
    catch
        true
    end

    platform["llvm"] = if llvm_assertions
        "$(llvm_version.major).asserts"
    else
        "$(llvm_version.major)"
    end

    return platform
end


text:
module GenerativeModels

using POMDPs

export generate_s,
       generate_o,
       generate_sr,
       generate_so,
       generate_or,
       generate_sor,
       initial_state

warn("GenerativeModels.jl has been merged into POMDPs.jl. Please import the functions from POMDPs rather than from GenerativeModels.")

end # module


text:
module DarkChaser

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Badeline Chaser" => Ahorn.EntityPlacement(
        Maple.DarkChaser
    ),

    "Badeline Chaser Barrier" => Ahorn.EntityPlacement(
        Maple.DarkChaserEnd,
        "rectangle"
    ),
)

Ahorn.minimumSize(entity::Maple.DarkChaserEnd) = 8, 8
Ahorn.resizable(entity::Maple.DarkChaserEnd) = true, true

# This sprite fits best, not perfect, thats why we have a offset here
chaserSprite = "characters/badeline/sleep00.png"

function Ahorn.selection(entity::Maple.DarkChaser)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(chaserSprite, x + 4, y, jx=0.5, jy=1.0)
end

function Ahorn.selection(entity::Maple.DarkChaserEnd)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return Ahorn.Rectangle(x, y, width, height)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.DarkChaserEnd)
    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    Ahorn.drawRectangle(ctx, 0, 0, width, height, (0.4, 0.0, 0.4, 0.4), (0.4, 0.0, 0.4, 1.0))
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.DarkChaser) = Ahorn.drawSprite(ctx, chaserSprite, 4, 0, jx=0.5, jy=1.0)

end


text:
#=
This model is taken from [1]. The theoretical derivations of the analytic solution can be found in [2].

[1] Malakiyeh, Mohammad Mahdi, Saeed Shojaee, and Klaus-Jürgen Bathe. "The Bathe time integration method revisited for prescribing desired numerical dissipation." Computers & Structures 212 (2019): 289-298.

[2] Mechanical Vibrations, Gerardin et al, page 250-251.
=#

using ReachabilityAnalysis, LinearAlgebra, LazySets, SparseArrays
using StructuralDynamicsODESolvers
using LazySets.Arrays

include("../krylov_concept.jl")
using FastExpm

LazySets.set_ztol(Float64, 1e-15)
LazySets.set_atol(Float64, 1e-15)
LazySets.set_rtol(Float64, 1e-15)

function clamped_bar(; N=1000, E=30e6, ρ=7.3e-4, A=1, L=200)
    ℓ = L / N
    K = (E * A / ℓ) * SymTridiagonal(fill(2.0, N), fill(-1.0, N))
    K[end, end] = E * A / ℓ

    M = (ρ * A * ℓ / 2) * Diagonal(vcat(fill(2.0, N-1), 1.0))
    M[end, end] = ρ * A * ℓ / 2

    return M, K
end

"""
    clamped_free(; N)

Free (no forcing term) instance of the clamped bar.
"""
function clamped_free(; N)
    M, K = clamped_bar(N=N)
    C = spzeros(N, N) # no damping
    sys = SecondOrderLinearContinuousSystem(M, C, K)
end

"""
    clamped_forced(; N, F=10e3, E=30e6, A=1)

Forced (constant force `F`) instance of the clamped bar.
"""
function clamped_forced(; N, F=10e3, E=30e6, A=1)
    M, K = clamped_bar(N=N)
    C = spzeros(N, N) # no damping
    F = vcat(zeros(N-1), F) # the right-most node is excited
    sys = SecondOrderAffineContinuousSystem(M, C, K, F)
end

# "nominal" step size
const dtn = 9.88e-7

# problem instances
s5 = clamped_forced(N=5);
s10 = clamped_forced(N=10);
s20 = clamped_forced(N=20);
s30 = clamped_forced(N=30);
s50 = clamped_forced(N=50);
s100 = clamped_forced(N=100);
s200 = clamped_forced(N=200);
s500 = clamped_forced(N=500);
s1000 = clamped_forced(N=1000);

# box initial condition
function prob_box(s, e)
    two_n = 2statedim(s)
    return InitialValueProblem(s, X0(two_n) + BallInf(two_n, e))
end

# solve with GLGM06
function solve_glg(s; δ=dtn, NSTEPS=10, model=Forward())
    # singleton initial condition
    prob = InitialValueProblem(s, Singleton(zeros(2statedim(s))))
    alg = GLGM06(δ=δ, approx_model=model)
    sol = solve(ivp, NSTEPS=NSTEPS, alg=alg, homogeneize=true)
end

# solve with LGG09 for nodeidx
function solve_lgg(s; δ=dtn, NSTEPS=100, nodeidx, model=Forward())
    # singleton initial condition
    prob = InitialValueProblem(s, Singleton(zeros(2statedim(s))))

    n = statedim(s)
    eplus = SingleEntryVector(nodeidx, 2n+1, 1.0)
    eminus = SingleEntryVector(nodeidx, 2n+1, -1.0)

    dirs = [eminus, eplus]

    alg = LGG09(δ=δ, template=dirs, approx_model=model)
    sol = solve(ivp, NSTEPS=NSTEPS, alg=alg, homogeneize=true)
end

# solve with BOX
# usage: sol = solve_box(s10, NSTEPS=100)
function solve_box(s; δ=dtn, NSTEPS=100, model=Forward())

    # ph = homogeneize(normalize(p))
    # @time pd = discretize(ph, dt, model) # ~35 sec for 1000 nodes

    alg = BOX(δ=δ, approx_model=model)

    # singleton initial condition
    prob = InitialValueProblem(s, Singleton(zeros(2statedim(s))))

    sbox = solve(p, alg=alg, NSTEPS=NSTEPS, homogeneize=true)
    return sbox
end

# Ec. (4.181) libro Geradin & Rixen (2015)
function _sol_analitica_clamped(x, t, Ncorte)

    F = 10E3
    E = 30E6
    dens = 7.3E-4
    L = 200
    A = 1

    # reescalar
    nMasses = 1000
    x = x*L/nMasses

    m = A*dens
    c = sqrt(E*A/m)

    α = 8*F*L/π^2/E/A

    β = π/2/L
    u = sin.( β * x ) * ( 1.0 .- cos.( β * c * t ) )
    v = sin.( β * x ) * β * c * sin.( β * c * t )

    for s in 2:Ncorte
        u .+= (-1)^(s-1)/(2*s-1)^2 * sin.( (2*s-1)* β * x ) * ( 1.0 .- cos.( (2*s-1)* β * c * t ) )
        v .+= (-1)^(s-1)/(2*s-1)^2 * sin.( (2*s-1)* β * x ) * ( (2*s-1)* β * c * sin.( (2*s-1)* β * c * t ) )
    end

    return u*α, v*α

end

function _sol_reachability_clamped(p; δt=dtn, NSTEPS=NSTEPS,
                                      nodeidx=[PosTwoThirds, VelTwoThirds],
                                      mk=150, tol=1e-12)

    # ==============================
    # Homogeneize
    # ==============================
    m = size(p.K, 1)
    Shom = _homogeneize(p)
    Aext = Shom.A
    @assert size(Aext, 1) == 2m+1

    # la condicion inicial es nula y se agrega la fuerza al final
    X0 = Vector{Float64}(undef, 2m+1)
    X0 = zeros(2m+1)
    X0[end] = 1.0

    # ==============================
    # Propagacion a tiempo continuo
    # ==============================
    #Aᵀδ = copy(transpose(Aext .* δt))
    num_nodes = isa(nodeidx, Int) ? 1 : length(nodeidx)

    #println("Calculando factor de bloating...")

    Eplus = _build_Eplus_optim(Aext, X0, δt; m=mk, tol=tol, cutoff=eps(Float64))

    #println("Propagando a tiempo continuo...")
    _resolver_clamped_expm(Aext, X0, m, δt, mk, tol, Eplus, num_nodes, NSTEPS, nodeidx)
end

function _resolver_clamped_expm(Aext, X0, m, δt, mk, tol, Eplus, num_nodes, NSTEPS, nodeidx)
    fastexpm = true

    if fastexpm
        Φᵀ = real.( fastExpm(transpose(Aext) .* δt, threshold=eps(Float64), nonzero_tol=eps(Float64)) )
    else
        Φᵀ = exp(Matrix(transpose(Aext) .* δt))
        Φᵀ = sparse(Φᵀ)
        droptol!(Φᵀ, eps(Float64))
    end

    q0 = transpose(transpose(X0) * Φᵀ) # para no calcular Φ y Φ^T

    # build Ω0 plus
    q0_msum_Eplus = Hyperrectangle(q0, Eplus.radius) # q0 ⊕ Eplus
    Ω0 = CH(Singleton(X0), q0_msum_Eplus)

    step_intersect = true
    if step_intersect
        Eplusstep = _build_Eplus_optim(Aext, q0, δt; m=mk, tol=tol, cutoff=eps(Float64))
        Ω0minus = CH(Singleton(q0), Singleton(X0) ⊕ Eplusstep)
        Ω0 = Ω0 ∩ Ω0minus
    end

    # Ω0box = box_approximation(Ω0, Hyperrectangle)

    Mout = Matrix{Float64}(undef, 2*num_nodes, NSTEPS)

    d = zeros(2m+1)
    @inbounds for j in 1:num_nodes
        node = nodeidx[j]
        d[node] = 1.0
    #    println("Propagando indice $node ...")
        idx = (2(j-1) + 1):2j
        reach_homog_dir_LGG09_modif!(view(Mout, idx, :), Ω0, Φᵀ, d, NSTEPS)
        d[node] = 0.0
    end
    δt0 = interval(0, δt)
    return Mout
end

function reach_homog_dir_LGG09_modif!(Mout, Ω₀, Φᵀ, d, NSTEPS)
    rᵢ = copy(d)
    rᵢ₊₁ = similar(rᵢ)

    @inbounds for i in 1:NSTEPS
        Mout[1, i] = ρ(rᵢ, Ω₀)
        Mout[2, i] = ρ(-rᵢ, Ω₀)

        # update cache for the next iteration
        mul!(rᵢ₊₁, Φᵀ, rᵢ)
        copy!(rᵢ, rᵢ₊₁)
    end
    return nothing
end

max_stars_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_issues_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_forks_repo_name: JuliaReach/SetPropagation-FEM-Examples


text:
function validate!(cur_plot::SimplePlot)
  if cur_plot.layout["xaxis"]["type"] == "log"
    for cur_data in cur_plot.data
      @assert all( cur_data["x"] .> 0 )
    end
  end

  if cur_plot.layout["yaxis"]["type"] == "log"
    for cur_data in cur_plot.data
      @assert all( cur_data["y"] .> 0 )
    end
  end

  cur_plot
end


text:
using QuantumSpin
using LinearAlgebra, SparseArrays
using Test

include("test_pauli_matrices.jl")
include("test_kronecker_product.jl")
include("test_pbc.jl")
include("test_lattice.jl")
include("test_set_spins.jl")
include("test_operators.jl")


text:
import Base: getindex, length, iterate
import Distributions.logpdf

"""
    Priors

Struct
```
struct Priors
    priors
    indicesForUpdt::Array{Array{Int64,1},1}
end
```
holds relevant information about priors: distributions that need to be used
and the pattern in which MCMC makes use of those distributions. `priors`
is a list of distributions with all priors that are to be used. `indicesForUpdt`
is a list of lists---the reason is as follows. Markov chains can have multiple
transition kernels, each updating (possibly) multiple number of parameters
(think of Gibbs sampler, which updates parameters in blocks)---for each such
update, priors for all updated parameters are needed. The inner list of
`indicesForUpdt` gives indices of all priors that together (in a factories form)
make up a prior for a respective transition kernel of the Markov chain

    Priors(priors, indicesForUpdt)
Base constructor. `priors` is a list of priors, `indicesForUpdt` is a pattern
of indices indicating in which way parameters are being updated

    Priors(priors)
Convenience constructor. Most of the time each kernel of the Markov chain will
be updating only a single parameter and thus only a single prior will be needed
for each transition. In that case providing a list of priors in `priors` is
sufficient. This constructor takes care of the internal objects in such setting.
"""
struct Priors{S,T}
    priors::S
    coord_idx::T
    function Priors(priors, coord_idx)
        priors, coord_idx = (priors,), (coord_idx,)
        S, T = typeof(priors), typeof(coord_idx)
        new{S,T}(priors, coord_idx)
    end
    function Priors(priors::T1, coord_idx::T2) where {T1<:Tuple, T2<:Tuple}
        S, T = typeof(priors), typeof(coord_idx)
        new{S,T}(priors, coord_idx)
    end
end

"""
    getindex(p::Priors, i::Int)

Fetch the `i`-th list of priors (corresponding to the `i`-th transition kernel)
from the `p` object.
"""
getindex(p::Priors, i::Int) = p.priors[i]

"""
    logpdf(p::Priors, θ)
Compute the logarithm of a product of all priors in object `p`, evaluated at `θ`
"""
function logpdf(p::Priors, θ)
    total = 0.0
    for (prior,coords) in p
        total += logpdf(prior, coords, θ)
    end
    total
end

"""
    length(p::Priors)
Return the total number of priors held by the object `p`
"""
length(p::Priors) = length(p.priors)

"""
    iterate(iter::Priors, state=(1, 0))

Iterates over sets of priors defined for separate parameter update steps
"""
function iterate(iter::Priors, i=1)
    if i > length(iter)
        return nothing
    end
    return ((iter.priors[i], iter.coord_idx[i]), i + 1)
end

"""
    ExpUnif

Distribution of an exponential of a uniform random variable U,
where U ~ Unif([low, up])
"""
struct ExpUnif
    low::Float64 # lower bound on a log scale
    up::Float64 # upper bound on a log scale

    ExpUnif(low_bd, up_bd) = new(low_bd, up_bd)
end

function logpdf(p::ExpUnif, coord_idx, θ)
    θi = θ[[indices(coord_idx)...]]
    @assert length(θi) == 1
    lθi = log(θi[1])
    p.low <= lθi <= p.up ? -log(p.up-p.low) - lθi : -Inf
end


text:
@testset "centroid" begin
    p1 = SVector(0.0, 0.0, 0.0)
    p2 = SVector(1.0, 0.0, 0.0)
    p3 = SVector(1.0, 1.0, 0.0)
    p4 = SVector(0.0, 1.0, 0.0)
	n̂ = triangleNormal(p1, p2, p3)

    cent = centroid(poly_eight(4, (p1, p2, p3, p4, p1, p1, p1, p1)), n̂)
    @test cent[1] == 1.0
    @test all(SVector(0.5, 0.5, 0.0) .== cent[2])

    cent = centroid(poly_eight(8, (p1, p2, p3, p4, p1, p1, p1, p1)), n̂)
    @test cent[1] == 1.0
    @test all(SVector(0.5, 0.5, 0.0) .== cent[2])

    cent = centroid(poly_eight(5, (p1, p2, p2, p3, p4, p1, p1, p1)), n̂)
    @test cent[1] == 1.0
    @test all(SVector(0.5, 0.5, 0.0) .== cent[2])

    cent = centroid(poly_eight(3, (p1, p2, p4, p1, p1, p1, p1, p1)), n̂)
    @test cent[1] == 0.5
    @test all(SVector(1/3, 1/3, 0.0) .== cent[2])

    cent = centroid(poly_eight(3, (p1, p2, p2, p1, p1, p1, p1, p1)), n̂)
    @test cent[1] == 0.0
    @test !any(isnan.(cent[2]))
end

@testset "poly_eight" begin
    for i_size = 1:8
        # zero_small_coordinates
        for i_vert = 1:i_size  # each filled vertex
            for i_ind = 1:4
                A = rand(8, 4) .+ 0.5  # all positive numbers
                A[i_vert, i_ind] = (rand() - 0.5) * 3.0e-15  # small number
                tuple_vertex = Tuple(SVector{4,Float64}(A[k, :]) for k = 1:8)
                p8 = poly_eight(i_size, tuple_vertex)
                p8 = zero_small_coordinates(p8)
                for i_vert_check = 1:i_size
                    for i_ind_check = 1:4
                        val_now = p8[i_vert_check][i_ind_check]
                        if (i_ind_check == i_ind) && (i_vert_check == i_vert)
                            @test val_now == 0.0
                        else
                            @test val_now == A[i_vert_check, i_ind_check]
                        end
                    end
                end
            end
        end

        # one_pad_then_mul
        p8 = poly_eight(i_size, Tuple(rand(SVector{3,Float64}) for k = 1:8))
        the_scale = SVector(2.0, 3.0, 4.0)
        m = SArray{Tuple{4,4},Float64,2,16}(diagm(0=>[the_scale..., 0.0]))
        p8_post = one_pad_then_mul(m, p8)

        for i_vert = 1:i_size
            @test all(the_scale .* p8[i_vert] .== p8_post[i_vert][1:3])
        end
    end
end


text:
##cd(joinpath(dirname(pathof(CombinedParsers)),"../docs/"))
push!(LOAD_PATH,"../src")
using Pkg
Pkg.activate(".")

testpath = joinpath("..","test")

using CombinedParsers
using CombinedParsers.Regexp
using BenchmarkTools

@info "defining pcre testset parser"
include(joinpath(testpath,"pcretest-parser.jl"))

testfile = joinpath(testpath,"testoutput1")
@info "loading pcre testset $testfile"
tests_string=read(testfile,String);
tests = pcre_tests(tests_string)[1];

ignore_idx = Dict{Int,String}(
    ## 14 => "unicode escape in test parser needs to support \"\\u81\".",
    69 => "because of very long compile time. The complex pattern parses email adresses",
    70 => "because of very long compile time. The complex pattern parses email adresses",
    89 => "CombinedParsers brackets use unicode character ranges. Avoid such patterns for now, contributions very welcome.",
    90 => "CombinedParsers brackets use unicode character ranges. Avoid such patterns for now, contributions very welcome.",
    135 => "escaped characters needs investigation.",
    664 => "because of very long matching time for `Repeat(Repeat1('a'))` when there is no match. Avoid such patterns, contributions optimizing these cases are also very welcome.", ## slooow
    887 => """do I misunderstand [recursive backreferences](https://www.pcre.org/original/doc/html/pcrepattern.html#SEC19)?

```@repl session
pattern = "(?P<abn>(?P=abn)xxx|)+"
match(Regex(pattern),"xxx")
match(Regcomb(pattern),"xxx")
```

IMO the pcre behaviour is confusing after considering the logic of expanding the repeat
```@repl session
pattern = "(?P<abn>(?P=abn)xxx|)((?P=abn)xxx|)"
match(Regex(pattern),"xxx")
match(Regcomb(pattern),"xxx")
```

""",
    104 => "case-ignoring backreferences still need to be done, contributions are very welcome,",
    260 => "case-ignoring backreferences still need to be done, contributions are very welcome,",
    261 => "case-ignoring backreferences still need to be done, contributions are very welcome,",
    562 => "case-ignoring backreferences still need to be done, contributions are very welcome,",
    1146 => "case-ignoring backreferences still need to be done, contributions are very welcome,",
    1173 => "because of very long compile time. The complex pattern parses PCRE regex syntax."
)

suite = BenchmarkGroup()
i,tt=1,tests[59]
p = tt[2].pattern
s = tt[2].test[1].sequence
for (i,tt) in Iterators.take(enumerate(tests),100)
    if !haskey(ignore_idx,i)
        ts = tt[2]
        name = string(ts.pattern)
        @info "setting up $i: $name"
        nam = string(ts.pattern)[1:min(end,20)]
        bs = suite["$i"] = BenchmarkGroup([name])
        p = ts.pattern
        pc = try
            pc_ = optimize(Regcomb(p...))
            bs[["Regcomb","create"]] = @benchmarkable Regcomb($p...)
            pc_
        catch e
            @warn "ignore Regcomb $i" exception=e
            nothing
        end
        pcre = try
            pc_ = Regex(p...)
            bs[["Regex","create"]] = @benchmarkable Regex($p...)
            pc_
        catch
            nothing
        end
        ##(j,test_seq) = first( enumerate(ts.test))
        for (j,test_seq) in enumerate(ts.test)
            s = test_seq.sequence
            pcre != nothing && try
                bs[["Regex","match", "$j"]] = @benchmarkable match($pcre,$s)
            catch
                @warn "ignore Regex $i, $j" s
            end
            pc != nothing && try
                s_ = CombinedParsers.Regexp.SequenceWithCaptures(s,pc)
                _iterate(pc,s_)
                bs[["Regcomb","match", "$j"]] = @benchmarkable _iterate($pc,$s_)
            catch
                @warn "ignore Regcomb $i, $j" s
            end
        end
        for (j,s) in enumerate(ts.tests_nomatch)
            pcre != nothing && try
                bs[["Regex", "nomatch", "$j"]] = @benchmarkable match($pcre,$s)
            catch
                @warn "ignore Regex $i, $j" s
            end
            pc != nothing && try
                s_ = CombinedParsers.Regexp.SequenceWithCaptures(s,pc)
                _iterate(pc,s_)
                bs[["Regcomb", "nomatch", "$j"]] = @benchmarkable _iterate($pc,$s_)
            catch
                @warn "ignore Regcomb $i, $j" s
            end
        end
    end
end

@info "tune!"
paramsfile = "benchmark-params.json"
if isfile(paramsfile)
    loadparams!(suite, BenchmarkTools.load(paramsfile)[1],:evals,:samples)
else
    tune!(suite, verbose=true)
    BenchmarkTools.save(paramsfile,params(suite))
end

results = run(suite, verbose=true,seconds=1)

import Dates
datetimenow = Dates.format(Dates.now(),"Y-mm-dd_HHhMM")
sha1 = parse((!re"[^\n]+"*'\n')[1], read(`git rev-parse HEAD`,String))
resultfile = "benchmark-$datetimenow-$sha1.json"
BenchmarkTools.save(joinpath(dirname(pathof(CombinedParsers)),"..","benchmark",resultfile),results)

using StructArrays, DataFrames
StructArray([ (pattern=t[1][1], string = t[1][4], code=t[1][2], trial = t[2])
              for t in leaves(median(results))
              if t[1][3]=="match" ]) |> DataFrame

results= nothing
GC.gc()


text:
module CurrencyFormatting

using ..Baskets
using ..CurrenciesBase
using ..CurrencyData
using ..DeclarativeFormatting
import ..DeclarativeFormatting: conflict, reconcile

export format

include("rules.jl")
include("decimals.jl")
include("templates.jl")
include("render.jl")

# Single-line output in two-argument `show`
Base.show(io::IO, m::Monetary) =
        show(IOContext(io, :compact => true), "text/plain", m)

function Base.show(io::IO, b::Basket)
    write(io, "$(typeof(b))([")
    write(io, join([sprint(show, c; context=:compact => true) for c in b], ","))
    print(io, "])")
end

# Multi-line output with MIME types
function Base.show(io::IO, ::MIME"text/plain", m::Monetary)
    if get(io, :compact, false)
        print(io, m / majorunit(m), currency(m))
    else
        print(io, format(m; styles=[:plain]))
    end
end

function Base.show(io::IO, ::MIME"text/latex", m::Monetary)
    print(io, string('$', format(m; styles=[:latex]), '$'))
end

function Base.show(io::IO, ::MIME"text/plain", b::Basket)
    len = length(b)
    write(io, "$len-currency $(typeof(b)):")
    for val in b
        write(io, "\n ")
        show(io, "text/plain", val)
    end
end

function Base.show(io::IO, ::MIME"text/markdown", b::Basket)
    len = length(b)
    println(io, "\$$len\$-currency `$(typeof(b))`:")
    for val in b
        write(io, "\n - ")
        show(io, "text/latex", val)
    end
end

end  # module CurrencyFormatting