text:
function compress_ids(cell_to_bgcell::AbstractVector{<:Integer})
  nccells = 0
  c = -1
  for bgcell in cell_to_bgcell
    if c != bgcell
      nccells += 1
      c = bgcell
    end
  end
  ccell_to_first_cell  = zeros(Int32,nccells+1)
  ccell = 0
  c = -1
  for (cell,bgcell) in enumerate(cell_to_bgcell)
    if c != bgcell
      ccell += 1
      c = bgcell
      ccell_to_first_cell[ccell] = cell
    end
  end
  ccell_to_first_cell[end] = length(cell_to_bgcell)+1
  ccell_to_first_cell
end

function compress_ids(cell_ids,cell_to_bgcell::AbstractVector{<:Integer})
  ccell_to_first_cell = compress_ids(cell_to_bgcell)
  nccells = length(ccell_to_first_cell)-1
  ccell_to_cell = view(ccell_to_first_cell,1:nccells)
  ccell_ids = view(cell_ids,ccell_to_cell)
  ccell_ids
end

function compress_contributions(
  cell_to_mat,
  cell_to_bgcell::AbstractVector{<:Integer},
  ccell_to_first_cell::AbstractVector{<:Integer}=compress_ids(cell_to_bgcell))

  nccells = length(ccell_to_first_cell)-1
  ccell_to_cell = view(ccell_to_first_cell,1:nccells)
  ccell_to_mat = CompressedCellArray(cell_to_mat,ccell_to_first_cell)
  ccell_to_mat
end

struct CompressedCellArray{T,A,B} <: AbstractVector{T}
  cell_to_array::A
  ccell_to_first_cell::B
  function CompressedCellArray(cell_to_array::A, ccell_to_first_cell::B) where {A,B}
    T = eltype(cell_to_array)
    new{T,A,B}(cell_to_array,ccell_to_first_cell)
  end
end

Base.size(a::CompressedCellArray) = (length(a.ccell_to_first_cell)-1,)
Base.IndexStyle(::Type{<:CompressedCellArray}) = IndexLinear()

function Base.getindex(a::CompressedCellArray,ccell::Integer)
  cell_first = a.ccell_to_first_cell[ccell]
  cell_last = a.ccell_to_first_cell[ccell+1]-1
  array = a.cell_to_array[cell_first]
  cell_range = (cell_first+1):cell_last
  _compress!(array,a.cell_to_array,cell_range)
end

@inline function _compress!(array,cell_to_array,cell_range)
  r = copy(array)
  for cell in cell_range
    carray = cell_to_array[cell]
    _addto!(r,carray)
  end
  r
end

@inline function _compress!(matvec::Tuple,cell_to_matvec,cell_range)
  mat, vec = matvec
  rm = copy(mat)
  rv = copy(vec)
  for cell in cell_range
    cmat, cvec = cell_to_matvec[cell]
    _addto!(rm,cmat)
    _addto!(rv,cvec)
  end
  (rm,rv)
end

function _addto!(a,b)
  @check size(a) == size(b)
  for k in eachindex(a)
    a[k] = a[k] + b[k]
  end
  a
end

function _addto!(a::ArrayBlock,b::ArrayBlock)
  @check size(a) == size(b)
  for k in eachindex(a.array)
    @check a.touched[k] == b.touched[k]
    if a.touched[k]
      _addto!(a.array[k],b.array[k])
    end
  end
  a
end

function array_cache(a::CompressedCellArray)
  array = testitem(a.cell_to_array)
  c = _compress_cache(array)
  ccache = array_cache(a.cell_to_array)
  c,ccache
end

function _compress_cache(array)
  c1 = CachedArray(copy(array))
  c2 = return_cache(Fields.unwrap_cached_array,c1)
  c1,c2
end

function _compress_cache(matvec::Tuple)
  mat, vec = matvec
  cmat = _compress_cache(mat)
  cvec = _compress_cache(vec)
  cmat, cvec
end

function getindex!(cache,a::CompressedCellArray,ccell::Integer)
  c, ccache = cache
  cell_first = a.ccell_to_first_cell[ccell]
  cell_last = a.ccell_to_first_cell[ccell+1]-1
  array = getindex!(ccache,a.cell_to_array,cell_first)
  cell_range = (cell_first+1):cell_last
  _compress!(c,array,ccache,a.cell_to_array,cell_range)
end

@inline function _compress!(c,array,ccache,cell_to_array,cell_range)
  c1,c2 = c
  _setsize_compress!(c1,array)
  r = evaluate!(c2,Fields.unwrap_cached_array,c1)
  copyto!(r,array)
  for cell in cell_range
    carray = getindex!(ccache,cell_to_array,cell)
    _addto!(r,carray)
  end
  r
end

@inline function _compress!(c,matvec::Tuple,ccache,cell_to_matvec,cell_range)
  cm, cv = c
  cm1,cm2 = cm
  cv1,cv2 = cv
  mat, vec = matvec
  _setsize_compress!(cm1,mat)
  _setsize_compress!(cv1,vec)
  rm = evaluate!(cm2,Fields.unwrap_cached_array,cm1)
  rv = evaluate!(cv2,Fields.unwrap_cached_array,cv1)
  copyto!(rm,mat)
  copyto!(rv,vec)
  for cell in cell_range
    cmat, cvec = getindex!(ccache,cell_to_matvec,cell)
    _addto!(rm,cmat)
    _addto!(rv,cvec)
  end
  (rm, rv)
end

function _setsize_compress!(a::CachedArray,b::AbstractArray)
  setsize!(a,size(b))
end

function _setsize_compress!(a::ArrayBlock,b::ArrayBlock)
  @check size(a) == size(b)
  for k in eachindex(a.array)
    @check a.touched[k] == b.touched[k]
    if a.touched[k]
    _setsize_compress!(a.array[k],b.array[k])
    end
  end
end


text:
using OASA,
      RDatasets,
      UnicodePlots,
      Base.Test

# Load data
# =========
d = Matrix(dataset("datasets", "iris"))

X = Float64.(d[:, 2:4])
Y = reshape(Float64.(d[:, 1]), (1, 150))

# Define model
# ============
# Variable "X" is an Mxn matrix of independent (endogenic) features
obj = L2_Obj()                          # L2 objective: ∑ (y - ŷ)^2
alg = Linear()                          # Linear model
sol = GradientDescent()                 # Full batch gradient descent
arc = SingleLayer(zeros(1,size(X,2)))   # 1xn matrix of weights

# OASA models also contain a key-value store for miscellaneous information
# generated during training and other events
inf = Dict()

# Instantiate model
reg = Model(obj, alg, sol, arc, inf)

# Test model
# ==========
# Ensure everything is set up properly
@testset "Model Initialized" begin
  # Dimensions of inference == dimensions of Y
  @test size(infer(reg, X)) == size(Y)
  # loss is producing a floating point number greater than 0
  @test assess(reg, Y, infer(reg, X)) > 0.0
  # lossgrad dimensions == 1xn matrix
  @test size(assessgrad(reg, Y, infer(reg, X), X)) == (1, size(X,2))
end

# Train model
# ===========
train!(reg, X, Y; learn_rate = 1e-6)

# Plot errors in terminal
lineplot(reg.inf["errors"], canvas=AsciiCanvas, border=:ascii, title="Objective")

# Also note that the components of `reg` are pointers to the variables
# defined earlier. We can access this information more directly:
# lineplot(inf["errors"])

ext: jl


text:
module Example_Household_Getter

using Model_Household
using DataFrames
using CSV

include( "hhld_from_frame.jl" )

export  initialise, get_household

EXAMPLE_HOUSEHOLDS = Dict{String,Household}()

KEYMAP = Vector{AbstractString}()

"""
return number of households available
"""
function initialise(
    ;
    household_name :: AbstractString = "example_households",
    people_name :: AbstractString = "example_people" ) :: Vector{AbstractString}

    global KEYMAP
    global EXAMPLE_HOUSEHOLDS

    hh_dataset = CSV.File("$(MODEL_DATA_DIR)/$(household_name).tab", delim='\t' ) |> DataFrame
    people_dataset = CSV.File("$(MODEL_DATA_DIR)/$(people_name).tab", delim='\t' ) |> DataFrame
    npeople = size( people_dataset)[1]
    nhhlds = size( hh_dataset )[1]
    for hseq in 1:nhhlds
        hhf = hh_dataset[hseq,:]
        push!( KEYMAP, hhf.name )
        EXAMPLE_HOUSEHOLDS[hhf.name] = load_hhld_from_frame( hseq, hhf, people_dataset )
    end
    KEYMAP
end

function example_names()
    global KEYMAP
    KEYMAP
end

function get_household( pos :: Integer ) :: Household
    global EXAMPLE_HOUSEHOLDS
    global KEYMAP
    key = KEYMAP[pos]
    EXAMPLE_HOUSEHOLDS[key]
end

function get_household( name :: AbstractString ) :: Household
    global EXAMPLE_HOUSEHOLDS
    EXAMPLE_HOUSEHOLDS[name]
end

end


text:
module Timeout

function timeout(f, timeout_sec)
    t = @async begin
        task = current_task()
        function timeout_cb(timer)
            Base.throwto(task, InterruptException())
        end
        timeout = Timer(timeout_cb, timeout_sec)
        try
            f()
        catch e
            if typeof(e) == InterruptException
                @error "Timed out after $(timeout_sec)s"
            end
            throw(e)
        end
        close(timeout)
    end
    try
        wait(t)
    catch ex
        throw(ex.task.exception)
    end
end

"""
    @timeout limit expr

Runs `expr` as and throws an `InterruptException` after `timeout_sec`.
## Example:
```julia
julia> @timeout 0.5 sleep(1)
┌ Error: Timed out after 0.5s
└ @ Main REPL[1]:12
ERROR: InterruptException:
Stacktrace:
 [1] timeout(::var"#8#9", ::Float64) at ./REPL[1]:21
 [2] top-level scope at REPL[3]:1
```
"""
macro timeout(limit, expr)
    :(timeout(()->$expr, $limit))
end

export @timeout

end # module


text:
"""
Trait that signals if a CellDatum type is implemented in the physical or the reference domain
"""
abstract type DomainStyle end
struct ReferenceDomain <: DomainStyle end
struct PhysicalDomain <: DomainStyle end

"""
Data associated with the cells of a Triangulation.
CellDatum objects behave as if they are defined in the physical space of the triangulation.
But in some cases they are implemented as reference quantities plus some transformation to the physical domain.
"""
abstract type CellDatum <: GridapType end

"""
Get the stored array of cell-wise data. It can be defined in the physical or the reference domain.
"""
get_data(a::CellDatum) = @abstractmethod

"""
Tell if the stored array is in the reference or physical domain
"""
DomainStyle(::Type{<:CellDatum}) = @abstractmethod

"""
Return the underlying Triangulation object
"""
get_triangulation(a::CellDatum) = @abstractmethod

"""
Change the underlying data to the target domain
"""
change_domain(a::CellDatum,target_domain::DomainStyle) = change_domain(a,DomainStyle(a),target_domain)
change_domain(a::CellDatum,input_domain::T,target_domain::T) where T<: DomainStyle = a
change_domain(a::CellDatum,input_domain::DomainStyle,target_domain::DomainStyle) = @abstractmethod

# Tester
"""
"""
function test_cell_datum(a::CellDatum)
  @test isa(get_data(a),AbstractArray)
  @test isa(get_triangulation(a),Triangulation)
  @test isa(DomainStyle(a),DomainStyle)
end

# Some API

DomainStyle(::T) where T<:CellDatum = DomainStyle(T)

"""
Get the raw array of cell data defined in the physical space.
"""
get_array(a::CellDatum) = get_data(change_domain(a,PhysicalDomain()))

"""
"""
#Base.length(a::CellDatum) = num_cells(get_triangulation(a))
num_cells(a::CellDatum) = num_cells(get_triangulation(a))


text:
##### deal with nutrients uptake
function gpu_calc_consume_kernel!(ctsdic, ctsdoc, ctsnh4, ctsno3, ctspo4, plank, ac, x, y, z, ΔT)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    for i = index:stride:size(ac,1)
        @inbounds CUDA.@atomic ctsdic[x[i], y[i], z[i]] += (plank.resp[i] - plank.PS[i]) * ΔT * ac[i]
        @inbounds CUDA.@atomic ctsdoc[x[i], y[i], z[i]] += (plank.exu[i] - plank.VDOC[i]) * ΔT * ac[i]
        @inbounds CUDA.@atomic ctsnh4[x[i], y[i], z[i]] += -plank.VNH4[i] * ΔT * ac[i]
        @inbounds CUDA.@atomic ctsno3[x[i], y[i], z[i]] += -plank.VNO3[i] * ΔT * ac[i]
        @inbounds CUDA.@atomic ctspo4[x[i], y[i], z[i]] += -plank.VPO4[i] * ΔT * ac[i]
    end
    return nothing
end
function calc_consume!(ctsdic, ctsdoc, ctsnh4, ctsno3, ctspo4, plank, ac, x, y, z, ΔT, ::GPU)
    @cuda threads=256 blocks=ceil(Int, size(ac,1)/256) gpu_calc_consume_kernel!(ctsdic, ctsdoc, ctsnh4,
                                                        ctsno3, ctspo4, plank, ac, x, y, z, ΔT)
    return nothing
end
function calc_consume!(ctsdic, ctsdoc, ctsnh4, ctsno3, ctspo4, plank, ac, x, y, z, ΔT, ::CPU)
    for i in 1:size(ac,1)
        @inbounds ctsdic[x[i], y[i], z[i]] += (plank.resp[i] - plank.PS[i]) * ΔT * ac[i]
        @inbounds ctsdoc[x[i], y[i], z[i]] += (plank.exu[i] - plank.VDOC[i]) * ΔT * ac[i]
        @inbounds ctsnh4[x[i], y[i], z[i]] += -plank.VNH4[i] * ΔT * ac[i]
        @inbounds ctsno3[x[i], y[i], z[i]] += -plank.VNO3[i] * ΔT * ac[i]
        @inbounds ctspo4[x[i], y[i], z[i]] += -plank.VPO4[i] * ΔT * ac[i]
    end
    return nothing
end

##### deal with grazed or dead individuals
function gpu_calc_loss_kernel!(ctsdoc, ctspoc, ctsdon, ctspon, ctsdop, ctspop, plank, ac, x, y, z,
                                   loss, lossFracC, lossFracN, lossFracP, R_NC, R_PC)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    for i = index:stride:size(ac,1)
        @inbounds CUDA.@atomic ctsdoc[x[i], y[i], z[i]] += (plank.Bm[i] + plank.Cq[i]) * lossFracC * ac[i] * loss[i]
        @inbounds CUDA.@atomic ctsdon[x[i], y[i], z[i]] += (plank.Bm[i]*R_NC + plank.Nq[i]) * lossFracN * ac[i] * loss[i]
        @inbounds CUDA.@atomic ctsdop[x[i], y[i], z[i]] += (plank.Bm[i]*R_PC + plank.Pq[i]) * lossFracP * ac[i] * loss[i]
        @inbounds CUDA.@atomic ctspoc[x[i], y[i], z[i]] += (plank.Bm[i] + plank.Cq[i]) * (1.0-lossFracC) * ac[i] * loss[i]
        @inbounds CUDA.@atomic ctspon[x[i], y[i], z[i]] += (plank.Bm[i]*R_NC + plank.Nq[i]) * (1.0-lossFracN) * ac[i] * loss[i]
        @inbounds CUDA.@atomic ctspop[x[i], y[i], z[i]] += (plank.Bm[i]*R_PC + plank.Pq[i]) * (1.0-lossFracP) * ac[i] * loss[i]
    end
    return nothing
end
function calc_loss!(ctsdoc, ctspoc, ctsdon, ctspon, ctsdop, ctspop, plank, ac, x, y, z,
                    loss, lossFracC, lossFracN, lossFracP, R_NC, R_PC, ::GPU)
    @cuda threads=256 blocks=ceil(Int, size(ac,1)/256) gpu_calc_loss_kernel!(ctsdoc, ctspoc, ctsdon, ctspon,
                                                        ctsdop, ctspop, plank, ac, x, y, z, loss,
                                                        lossFracC, lossFracN, lossFracP, R_NC, R_PC)
    return nothing
end
function calc_loss!(ctsdoc, ctspoc, ctsdon, ctspon, ctsdop, ctspop, plank, ac, x, y, z,
                           loss, lossFracC, lossFracN, lossFracP, R_NC, R_PC, ::CPU)
    for i in 1:size(ac,1)
        @inbounds ctsdoc[x[i], y[i], z[i]] += (plank.Bm[i] + plank.Cq[i]) * lossFracC * ac[i] * loss[i]
        @inbounds ctsdon[x[i], y[i], z[i]] += (plank.Bm[i]*R_NC + plank.Nq[i]) * lossFracN * ac[i] * loss[i]
        @inbounds ctsdop[x[i], y[i], z[i]] += (plank.Bm[i]*R_PC + plank.Pq[i]) * lossFracP * ac[i] * loss[i]
        @inbounds ctspoc[x[i], y[i], z[i]] += (plank.Bm[i] + plank.Cq[i]) * (1.0-lossFracC) * ac[i] * loss[i]
        @inbounds ctspon[x[i], y[i], z[i]] += (plank.Bm[i]*R_NC + plank.Nq[i]) * (1.0-lossFracN) * ac[i] * loss[i]
        @inbounds ctspop[x[i], y[i], z[i]] += (plank.Bm[i]*R_PC + plank.Pq[i]) * (1.0-lossFracP) * ac[i] * loss[i]
    end
    return nothing
end


text:
using JuMP
using StructJuMP
using CPLEX
using ECOS
using FactCheck

include("../src/BendersBridge.jl")

misocp_solver = CplexSolver()
socp_solver = ECOS.ECOSSolver()

facts("[Benders] Empty scenario test") do

    m = StructuredModel(num_scenarios=0)
    @defVar(m, x, Int)
    @addConstraint(m, x <= 4)
    @setObjective(m, :Min, -5*x)

    output = BendersBridge(m, misocp_solver, socp_solver)

    @fact output[1] --> :Optimal
    @fact output[2] --> roughly(-20.0)

end

facts("[Benders] Infeasible problem test") do

    numScen = 1
    m = StructuredModel(num_scenarios=numScen)

    @defVar(m, x, Int)

    @addConstraint(m, x <= 1)
    @setObjective(m, :Min, -5*x)

    bl = StructuredModel(parent=m)
    @defVar(bl, y1 >= 2)
    @defVar(bl, y2 <= 2)
    @addConstraint(bl, x >= y1)
    @addConstraint(bl, norm(y1) <= y2)
    @setObjective(bl, Min, 2*y1 + y2)

    output = BendersBridge(m, misocp_solver, socp_solver)

    @fact output[1] --> :Infeasible

end

facts("[Benders] Infeasibility cut execution test #1") do

    numScen = 1
    m = StructuredModel(num_scenarios=numScen)

    @defVar(m, x, Int)

    @addConstraint(m, x <= 4)
    @setObjective(m, :Min, -5*x)

    bl = StructuredModel(parent=m)
    @defVar(bl, y1 >= 0)
    @defVar(bl, y2 <= 2)
    @addConstraint(bl, x <= y1)
    @addConstraint(bl, norm(y1) <= y2)
    @setObjective(bl, Min, 2*y1 + y2)

    output = BendersBridge(m, misocp_solver, socp_solver)

    @fact output[1] --> :Optimal
    @fact output[2] --> roughly(-4.0)

end

facts("[Benders] Optimality cut execution test #1") do

    numScen = 1
    m = StructuredModel(num_scenarios=numScen)

    @defVar(m, x, Int)

    @addConstraint(m, x <= 4)
    @setObjective(m, :Min, -5*x)

    bl = StructuredModel(parent=m)
    @defVar(bl, y1 >= 2)
    @defVar(bl, y2 <= 4)
    @addConstraint(bl, x <= y1)
    @addConstraint(bl, norm(y1) <= y2)
    @setObjective(bl, Min, 2*y1 + y2)

    output = BendersBridge(m, misocp_solver, socp_solver)

    @fact output[1] --> :Optimal
    @fact output[2] --> roughly(-8.0)

end


text:
module Assignment07

export normalizeDNA,
        composition,
        gc_content,
        complement,
        reverse_complement,
        parse_fasta

# # uncomment the following line if you intend to use BioSequences types
# using BioSequences

"""
    normalizeDNA(::AbstractString)

Ensures that a sequence only contains valid bases
(or `'N'` for unknown bases).
Returns a String.
"""
function normalizeDNA(seq)
    seq = uppercase(string(seq))
    for base in seq
        # note: `N` indicates an unknown base
        occursin(base, "AGCTN") || error("invalid base $base")
    end
    return seq # change to `return LongDNASeq(seq)` if you want to try to use BioSequences types
end

# Your code here.
# Don't forget to export your functions!

function composition(seq)
    seq = normalizeDNA(seq) # make uppercase string, check invalid bases

    base_dict = Dict('A'=>  0, 'C' => 0, 'G' => 0, 'T'=> 0, 'N' => 0) #create dictionary w/ keys A,C,T,G,N
    for base in seq
         if base == 'A'
            base_dict['A'] = base_dict['A'] +1
         elseif base == 'C'
            base_dict['C'] = base_dict['C'] +1
         elseif base == 'G'
            base_dict['G'] = base_dict['G']+1
        elseif base == 'T'
            base_dict['T'] = base_dict['T']+1
         elseif base == 'N'
            base_dict['N'] = base_dict['N']+1
         end
        end
    return base_dict

end

function gc_content(seq)
    base_dict=composition(seq)
    c=base_dict['C']
    g=base_dict['G']

    total=base_dict['N'] + base_dict['G'] + base_dict['T'] + base_dict['C'] + base_dict['A']   #same as total=content[1]+content[2]+content[3]+content[4]
    #c in second position of the tuple
    #g in third position
    #total is sum of all positions
    return (c+g)/total
end

 function complement(seq::String)
    seq= normalizeDNA(seq)
    complements= Dict('A' => 'T', 'T' => 'A', 'G' => 'C', 'C' => 'G', 'N' => 'N')
    my_complement=[]
    for base in seq
        if occursin(base, "AGTCN")
        push!(my_complement, complements[base])
        end
    end

return join(my_complement)
end

function reverse_complement(seq)
        sequence = reverse(uppercase(seq))
        arr=[]
        for i in sequence
        push!(arr, string(i))  #pushes each base into the array as a string (b/c complement only takes strings)
        end

        comp = map(complement, arr)
        join(comp) #turns array back into a string
end

function parse_fasta(path)

    header_vector = []
    combined_seq_vector = []
    currentvector = [] #for sequences
    for line in eachline(path)
        if startswith(line, ">")
            header = lstrip(line, ['>'])
            push!(header_vector, header) #pushes the header string into the vector
            joined_string = join(currentvector)
            push!(combined_seq_vector, joined_string)
            currentvector = []
        else
            line = strip(line)
            normalizeDNA(line) #only allow valid DNA sequences including N to enter seq vector
            push!(currentvector, line)

        end
    end
    joined_string = join(currentvector)
    push!(combined_seq_vector, joined_string)#need to deal with last sequence
    return header_vector, combined_seq_vector[2:end] #dont return the empty string from the first header
end

end # module Assignment07


text:
# Some radial basis functions
rbf_thin_plate{T <: FloatingPoint}(r::T, c::T = 1.0) = (r == 0.0) ? 0.0 : (r^2 * log(c*r))

rbf_gaussian{T <: FloatingPoint}(r::T, c::T = 1.0) = exp(-c*r^2)

rbf_cubic{T <: FloatingPoint}(r::T, c::T = 1.0) = r^3

rbf_multiquadric{T <: FloatingPoint}(r::T, c::T = 1.0) = sqrt(r^2 + c^2)

identity_function{T <: FloatingPoint}(r::T) = r

euclidean_distance{T <: FloatingPoint}(x1::Array{T, 1}, x2::Array{T, 1}) = norm(x1 - x2)

# Create a kernel function from a RBF, a distance function and a c constant that
# is passed on to the RBF.
function create_kernel_func(rbffunc, distancefunc = euclidean_distance, c = 1.0)
  (x1, x2) -> rbffunc(distancefunc(x1, x2), c)
end

# Calculate the pairwise kernel matrix for the points in x taken to be the
# individual columns. The kernel matrix can be padded with zeros by specifying
# num_padding_columns and num_padding_rows
function calc_kernel_matrix{T <: FloatingPoint}(x::Array{T, 2}, centers::Array{T, 2},
  kernelfunc = create_kernel_func(rbf_gaussian, euclidean_distance, 1/size(x, 1)))

  num_points = size(x, 2)
  num_centers = size(centers, 2)
  k_matrix = zeros(Float64, num_points, num_centers)

  for i in 1:num_points
    xi = x[:,i]
    for j in 1:num_centers
      k_matrix[i, j] = kernelfunc(xi, centers[:,j])
    end
  end

  k_matrix

end

# Scale and center with the mean and standard deviation (or with specific
# center and value vectors given as input (if given)).
function scale_and_center{T <: FloatingPoint}(x::Array{T, 2}, center = true, scale = true)

  if center == true
    center = mean(x, 2)
    xres = broadcast(-, x, center)
  elseif typeof(center) == Array{T, 2}
    if size(center, 1) != size(x, 1)
      throw("When giving a vector of center values they must have as many rows as the number of rows per point; here they differ!")
    end
    xres = broadcast(-, x, center)
  else
    xres = x
  end

  if scale == true
    scale = std(x, 2)
    xres = broadcast(/, xres, scale)
  elseif typeof(scale) == Array{T, 2}
    if size(scale, 1) != size(x, 1)
      throw("When giving a vector of scale values they must have as many rows as the number of rows per point; here they differ!")
    end
    xres = broadcast(/, xres, scale)
  end

  return xres, center, scale

end

type RadialBasisFunctionModel{T <: FloatingPoint}
  kernel::Function

  xcenters::Array{T, 2}     # X values for the centers. One center per column.
  x_row_means::Array{T, 2}  # Average value for the original, uncentered and unscaled xcenters.
  x_row_stds::Array{T, 2}   # Standard dev for the original, uncentered and unscaled xcenters.

  fcenters::Array{T, 2}     # Y values for the centers. One value per as many columns as in xcenters.
  f_mean::T                 # Mean value for the original, uncentered and unscaled function values.
  f_std::T                  # Standard dev for the original, uncentered and unscaled function values.

  lambdas::Array{T, 2}      # Coefficients in the RBF model.
end

function RBFModel{T <: FloatingPoint}(x::Array{T,2}, y::Array{T,2};
  rbf = rbf_gaussian,
  c = nothing,
  distance_func = euclidean_distance)

  # Set c to ok values unless it was specified
  if c == nothing && rbf == rbf_gaussian
    num_dims = size(x, 1)
    c = 1/(num_dims * num_dims)
  elseif c == nothing
    c = 1.0
  end

  kernel = create_kernel_func(rbf, distance_func, c)

  # Scale and center x and y (per row, i.e. covariate)
  xcenters, x_row_means, x_row_stds = scale_and_center(x)
  fcenters, f_mean, f_std = scale_and_center(y)

  # Create kernel matrix.
  A = calc_kernel_matrix(xcenters, xcenters, kernel)

  # Solve for the lambda coefficients.
  lambdas = A\(fcenters')

  RadialBasisFunctionModel(kernel, xcenters, x_row_means, x_row_stds, fcenters, f_mean[1], f_std[1], lambdas)

end

function predict{T <: FloatingPoint}(rbfm::RadialBasisFunctionModel, xx::Array{T, 2})
  # Scale xx
  xx_scaled, c, s = scale_and_center(xx, rbfm.x_row_means, rbfm.x_row_stds)

  # Predict for the test points in xx
  B = calc_kernel_matrix(xx_scaled, rbfm.xcenters, rbfm.kernel)
  ypred_scaled = (B * rbfm.lambdas)'

  # Transform back to original range
  ypred = rbfm.f_mean + rbfm.f_std * ypred_scaled

  return ypred
end

# Convenience function that builds the model from x and y and then uses the model
# to predict the y values for the points in xx.
function rbf_predict{T <: FloatingPoint}(x::Array{T, 2}, y::Array{T, 2},
  xx::Array{T, 2}; options...)

  rbfm = RBFModel(x, y; options...)
  ypred = predict(rbfm, xx)
  return ypred, rbfm

end


text:
using SumOfSquares
using DynamicPolynomials
using Test

include("sparse_polynomials.jl")

function test(test_function, optimizer_constructor)
    sparse_model = sos_lower_bound(test_function, optimizer_constructor, VariableSparsity())
    dense_model  = sos_lower_bound(test_function, optimizer_constructor, NoSparsity())
    return abs(objective_value(sparse_model) - objective_value(dense_model))
end

function chordal_sos(optimizer_constructor, tol=1e-5)
    val = test(chained_singular(4), optimizer_constructor)
    @test val < tol

    val = test(broyden_banded(2), optimizer_constructor)
    @test val < tol

    val = test(broyden_tridiagonal(4), optimizer_constructor)
    @test val < tol

    val = test(chained_wood(4), optimizer_constructor)
    @test val < tol

    val = test(generalized_rosenbrock(4), optimizer_constructor)
    @test val < tol
end

import CSDP
chordal_sos(CSDP.Optimizer)


text:
# This file was generated, do not modify it. # hide
scitype_union(X.age)

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
using Logging, Test, HSM

struct ActiveStateTestState <: HSM.AbstractHsmState
    state_info::HSM.HsmStateInfo

    ActiveStateTestState(parent_state) = new(HSM.HsmStateInfo(parent_state))
end

"""
    basic_state_machine()

Builds and returns the following state machine:

```plantuml
@startuml
state A {
    state B
    state C {
        state D
        state E
    }
}
@enduml
```
"""
function active_state_test_state_machine()
    A = ActiveStateTestState(nothing)
    B = ActiveStateTestState(A)
    C = ActiveStateTestState(A)
    D = ActiveStateTestState(C)
    E = ActiveStateTestState(C)

    return A, B, C, D, E
end

@testset "active_state()" begin
    A, B, C, D, E = active_state_test_state_machine()
    HSM.transition_to_state!(A, B)

    @test HSM.active_state(A) == B
    @test HSM.active_state(B) == B
    @test HSM.active_state(C) == B
    @test HSM.active_state(D) == B
    @test HSM.active_state(E) == B
end


text:
# ---
# title: 141. Linked List Cycle
# id: problem141
# author: zhwang
# date: 2022-02-15
# difficulty: Easy
# categories: Linked List, Two Pointers
# link: <https://leetcode.com/problems/linked-list-cycle/description/>
# hidden: true
# ---
#
# Given `head`, the head of a linked list, determine if the linked list has a
# cycle in it.
#
# There is a cycle in a linked list if there is some node in the list that can
# be reached again by continuously following the `next` pointer. Internally,
# `pos` is used to denote the index of the node that tail's `next` pointer is
# connected to.  **Note that  `pos` is not passed as a parameter**.
#
# Return `true` _if there is a cycle in the linked list_. Otherwise, return
# `false`.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2018/12/07/circularlinkedlist.png)
#
#
#
#     Input: head = [3,2,0,-4], pos = 1
#     Output: true
#     Explanation: There is a cycle in the linked list, where the tail connects to the 1st node (0-indexed).
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2018/12/07/circularlinkedlist_test2.png)
#
#
#
#     Input: head = [1,2], pos = 0
#     Output: true
#     Explanation: There is a cycle in the linked list, where the tail connects to the 0th node.
#
#
# **Example 3:**
#
# ![](https://assets.leetcode.com/uploads/2018/12/07/circularlinkedlist_test3.png)
#
#
#
#     Input: head = [1], pos = -1
#     Output: false
#     Explanation: There is no cycle in the linked list.
#
#
#
#
# **Constraints:**
#
#   * The number of the nodes in the list is in the range `[0, 104]`.
#   * `-105 <= Node.val <= 105`
#   * `pos` is `-1` or a **valid index** in the linked-list.
#
#
#
# **Follow up:** Can you solve it using `O(1)` (i.e. constant) memory?
#
#
## @lc code=start
using LeetCode
function has_cycle(head::Union{ListNode,Nothing})::Bool
    isnothing(head) && return false
    swift = slow = head
    while !isnothing(swift.next) && !isnothing(swift.next.next)
        slow, swift = slow.next, swift.next.next
        slow === swift && return true
    end
    return false
end
## @lc code=end


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

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Distributions, Random, StatsBase, DataFrames, Plots
	using QuadGK, SpecialFunctions
	Random.seed!(1)
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 3.25"

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
g(x) = x^(0.5-1) * MathConstants.e^-x

# ╔═╡ f8ef8fb4-21bf-11eb-2f3e-73cc8d30faaf
quadgk(g,0,Inf)[1], sqrt(pi), gamma(1/2),  factorial(1/2-1)

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 3.25"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╠═f8ef8fb4-21bf-11eb-2f3e-73cc8d30faaf
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3


text:
#========================================================================================#
#	Laboratory 300
#
# Welcome to course 300: Constructing Dynamical Models!
#
# Author: Niall Palfreyman, 24/04/2022
#========================================================================================#
[
	Activity(
		"""
		???
		""",
		"???",
		x -> true
	),
]


text:
################################################################################
#
#	CONCRETE TYPE
#
#   ReciprocalPoint{D} <: AbstractReciprocalPoint{D}
#   --> D is the dimension of embedding space
#
#   FILE CONTAINS
#       - concrete struct definition
#       - interface implementation
#
################################################################################

################################################################################
#
#   CONCRETE STRUCT DEFINITION
#
################################################################################
mutable struct ReciprocalPoint{D} <: AbstractReciprocalPoint{D}

	# point
	point       :: Vector{Float64}

	# label
	label	    :: String
	# LaTeX label
	label_LaTeX :: String

end

# export the concrete type
export ReciprocalPoint

################################################################################
#
#	IMPLEMENTATION OF INTERFACE FOR CONCRETE RECIPROCAL POINT TYPE
#	(functions that had to be overwritten by concrete type)
#
################################################################################

# default constructor interface
# used for creation of new reciprocal points
function newReciprocalPoint(
            :: Type{ReciprocalPoint{D}},
            point   :: Vector{<:Real},
            label   :: AbstractString
        ) :: ReciprocalPoint{D} where {D}

    # check if correct dimension is given
    @assert length(point) == D

    # return a new object
    return ReciprocalPoint{D}(point, label, label)
end
function newReciprocalPoint(
            :: Type{ReciprocalPoint{D}},
            point       :: Vector{<:Real},
            label       :: AbstractString,
            label_LaTeX :: AbstractString
        ) :: ReciprocalPoint{D} where {D}

    # check if correct dimension is given
    @assert length(point) == D

    # return a new object
    return ReciprocalPoint{D}(point, label, label_LaTeX)
end

# get label
function label(
            p :: ReciprocalPoint{D}
        ) :: String where {D}

    # return the label
    return p.label
end
# set label
function label!(
            p :: ReciprocalPoint{D},
            l :: AbstractString
        ) where {D}

    # set the label
    p.label = l
end

# get label (LaTeX)
function labelLaTeX(
            p :: ReciprocalPoint{D}
        ) :: String where {D}

    # return the label
    return p.label_LaTeX
end
# set label (LaTeX)
function labelLaTeX!(
            p :: ReciprocalPoint{D},
            l :: AbstractString
        ) where {D}

    # set the label
    p.label_LaTeX = l
end

# get point
function point(
            p :: ReciprocalPoint{D}
        ) :: Vector{Float64} where {D}

    # return the point
    return p.point
end
# set point
function point!(
            p :: ReciprocalPoint{D},
            r :: Vector{<:Real}
        ) where {D}

    # set the point
    p.point = r
end


text:
@testset "Integer Conic" begin
    mock = MOIU.MockOptimizer(Model{Float64}())
    config = MOIT.TestConfig()

    @testset "SOC" begin
        mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 1.0, 0.0])
        MOIT.intsoc1test(mock, config)
    end
end


text:
# Nos ubicamos en la direccion del repositorio
cd(@__DIR__)
# Haciendo uso de Pkg.jl
using Pkg
# Inicializamos el entorno de este repositorio
Pkg.activate(".")
# Inicializamos los tests
Pkg.test()

max_stars_repo_name: vntr-CC/Forward-Football

max_issues_repo_name: vntr-CC/Forward-Football

max_forks_repo_name: vntr-CC/Forward-Football


text:
## define features traits

# most general feature(roadway, scenes, actions, vehid)
abstract type AbstractFeature end

# feature(roadway, veh)
struct EntityFeature <: AbstractFeature end
# feature(roadway, scene, veh)
struct SceneFeature <: AbstractFeature end
# feature(roadway, scenes, veh)
struct TemporalFeature <: AbstractFeature end
# feature(actions, veh)
# struct TemporalActionFeature <: AbstractFeature end
# feature(action, veh)
# struct ActionFeature <: AbstractFeature end

"""
    featuretype(f::Function)
function used for dispatching feature types depending on the method of the feature function
"""
function featuretype(f::Function)
    if static_hasmethod(f, Tuple{Roadway, Entity})
        return EntityFeature()
    elseif static_hasmethod(f, Tuple{Roadway, Scene, Entity})
        return SceneFeature()
    elseif static_hasmethod(f, Tuple{Roadway, Vector{<:Scene}, Entity})
        return TemporalFeature()
    else
        error(""""
              unsupported feature type
              try adding a new feature trait or implement one of the existing feature type supported
              """)
    end
    return
end

## Feature extraction functions

# Top level extract several features
"""
    extract_features(features, roadway::Roadway, scenes::Vector{<:Scene}, ids::Vector{I}) where I

Extract information from a list of scenes. It returns a dictionary of DataFrame objects.

# Inputs

- `features`: a tuple of feature functions. The feature types supported are `EntityFeature`, `SceneFeature` and `TemporalFeature`. Check the documentation for the list of available feature functions, and how to implement you own feature function.
- `roadway::Roadway` : the roadway associated to the scenes.
- `scenes::Vector{<:Scene}` : the simulation data from which we wish to extract information from. Each scene in the vector corresponds to one time step.
- `ids::Vector{I}`: a list of entity IDs for which we want to extract the information.

# Output

A dictionary mapping IDs to `DataFrame` objects. For a given ID, the DataFrame's columns correspond to the name of the input feature functions.
The row correspond to the feature value for each scene (time history).

# Example:

```julia
roadway = gen_straight_roadway(4, 100.0)

scene_0 = Scene([
            Entity(VehicleState(VecSE2( 0.0,0.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 1),
            Entity(VehicleState(VecSE2(10.0,0.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 2),
        ])

scene_1 = Scene([
            Entity(VehicleState(VecSE2( 10.0,0.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 1),
            Entity(VehicleState(VecSE2(20.0,0.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 2),
        ])

dfs = extract_features((posgx, iscolliding), roadway, [scene_0, scene_1], [1,2])

dfs[1].posgx # history of global x position for vehicle of ID 1
```
"""
function extract_features(features, roadway::Roadway, scenes::Vector{<:Scene}, ids::Vector{I}) where I
    dfs_list = broadcast(f -> extract_feature(featuretype(f), f, roadway, scenes, ids), features)
    # join all the dataframes
    dfs = Dict{I, DataFrame}()
    for id in ids
        dfs[id] = DataFrame([Symbol(features[i]) => dfs_list[i][id] for i=1:length(features)])
    end
    return dfs
end

# feature is a general (temporal) feature, broadcast on id only
function extract_feature(ft::AbstractFeature, feature, roadway::Roadway, scenes::Vector{<:Scene}, ids)
    values = feature.(Ref(roadway), Ref(scenes), ids)  # vector of vectors
    return Dict(zip(ids, values))
end

# feature is a scene feature, broadcast on ids and scenes
function extract_feature(ft::SceneFeature, feature, roadway::Roadway, scenes::Vector{<:Scene}, ids)
    values = broadcast(ids) do id
        veh_hist = get_by_id.(scenes, Ref(id))
        feature.(Ref(roadway), scenes, veh_hist)
    end
    return Dict(zip(ids, values))
end

# feature is an entity feature, broadcast on ids and scenes
function extract_feature(ft::EntityFeature, feature, roadway::Roadway, scenes::Vector{<:Scene}, ids)
    values = broadcast(ids) do id
        veh_hist= get_by_id.(scenes, Ref(id))
        feature.(Ref(roadway), veh_hist)
    end
    return Dict(zip(ids, values))
end

## feature functions

# vehicle only features
"""
    posgx(roadway::Roadway, veh::Entity)
returns x position in the global frame
"""
posgx(roadway::Roadway, veh::Entity) = posg(veh).x

"""
    posgy(roadway::Roadway, veh::Entity)
returns y position in the global frame
"""
posgy(roadway::Roadway, veh::Entity) = posg(veh).y

"""
    posgθ(roadway::Roadway, veh::Entity)
returns heading in the global frame
"""
posgθ(roadway::Roadway, veh::Entity) = posg(veh).θ

"""
    posfs(roadway::Roadway, veh::Entity)
returns longitudinal position in Frenet frame
"""
posfs(roadway::Roadway, veh::Entity) = posf(veh).s

"""
    posft(roadway::Roadway, veh::Entity)
returns lateral position in Frenet frame
"""
posft(roadway::Roadway, veh::Entity) = posf(veh).t

"""
    posfϕ(roadway::Roadway, veh::Entity)
returns heading relative to centerline of the lane
"""
posfϕ(roadway::Roadway, veh::Entity) = posf(veh).ϕ

"""
    vel(roadway::Roadway, veh::Entity)
returns the velocity of `veh`
"""
vel(roadway::Roadway, veh::Entity) = vel(veh)

"""
    velfs(roadway::Roadway, veh::Entity)
returns the longitudinal velocity in the Frenet frame
"""
velfs(roadway::Roadway, veh::Entity) = velf(veh).s

"""
    velft(roadway::Roadway, veh::Entity)
returns the lateral velocity in the Frenet frame
"""
velft(roadway::Roadway, veh::Entity) = velf(veh).t

"""
    velgx(roadway::Roadway, veh::Entity)
returns the velocity in the x direction of the global frame
"""
velgx(roadway::Roadway, veh::Entity) = velg(veh).x

"""
    velgy(roadway::Roadway, veh::Entity)
returns the velocity in the y direction of the global frame
"""
velgy(roadway::Roadway, veh::Entity) = velg(veh).y

"""
    time_to_crossing_right(roadway::Roadway, veh::Entity)
time before crossing the right lane marker assuming constant lateral speed.
"""
function time_to_crossing_right(roadway::Roadway, veh::Entity)
    d_mr = markerdist_right(roadway, veh)
    velf_t = velf(veh).t
    if d_mr > 0.0 && velf_t < 0.0
        return -d_mr / velf_t
    else
        return Inf
    end
end

"""
    time_to_crossing_left(roadway::Roadway, veh::Entity)
time before crossing the left lane marker assuming constant lateral speed.
"""
function time_to_crossing_left(roadway::Roadway, veh::Entity)
    d_mr = markerdist_left(roadway, veh)
    velf_t = velf(veh).t
    if d_mr > 0.0 && velf_t > 0.0
        return -d_mr / velf_t
    else
        return Inf
    end
end

"""
    estimated_time_to_lane_crossing(roadway::Roadway, veh::Entity)
minimum of the time to crossing left and time to crossing right.
"""
function estimated_time_to_lane_crossing(roadway::Roadway, veh::Entity)
    ttcr_left = time_to_crossing_left(roadway, veh)
    ttcr_right = time_to_crossing_right(roadway, veh)
    return min(ttcr_left, ttcr_right)
end

"""
    iswaiting(roadway::Roadway, veh::Entity)
returns true if the vehicle is waiting (if the velocity is 0.0)
"""
iswaiting(roadway::Roadway, veh::Entity) = vel(veh) ≈ 0.0

## scene features

"""
    iscolliding(roadway::Roadway, scene::Scene, veh::Entity)
returns true if the vehicle is colliding with another vehicle in the scene.
"""
function iscolliding(roadway::Roadway, scene::Scene, veh::Entity)
    return collision_checker(scene, veh.id)
end

"""
    distance_to(egoid)
generate a feature function distance_to_\$egoid.

`distance_to_\$egoid(roadway, scene, veh)` returns the distance between veh and the vehicle of id egoid in the scene.
"""
function distance_to(egoid)
    @eval begin
        function ($(Symbol("distance_to_$(egoid)")))(roadway::Roadway, scene::Scene, veh::Entity)
            ego = get_by_id(scene, $egoid)
            return norm(VecE2(posg(ego) - posg(veh)))
        end
        export $(Symbol("distance_to_$(egoid)"))
        $(Symbol("distance_to_$(egoid)"))
    end
end

# temporal features

"""
    acc(roadway::Roadway, scenes::Vector{<:Scene}, vehid)
returns the history of acceleration of the entity of id vehid using finite differences on the velocity.
The first element is `missing`.
"""
function acc(roadway::Roadway, scenes::Vector{<:Scene}, vehid)
    result = zeros(Union{Missing,Float64},length(scenes))
    angles = broadcast(x->vel(get_by_id(x, vehid)), scenes)
    result[1] = missing
    result[2:end] .= diff(angles)
    return result
end

"""
    accfs(roadway::Roadway, scenes::Vector{<:Scene}, vehid)
returns the history of longitudinal acceleration in the Frenet frame of the entity of id vehid using finite differences on the velocity.
The first element is `missing`.
"""
function accfs(roadway::Roadway, scenes::Vector{<:Scene}, vehid)
    result = zeros(Union{Missing,Float64},length(scenes))
    angles = broadcast(x->velf(get_by_id(x, vehid)).s, scenes)
    result[1] = missing
    result[2:end] .= diff(angles)
    return result
end

"""
    accft(roadway::Roadway, scenes::Vector{<:Scene}, vehid)
returns the history of lateral acceleration in the Frenet frame of the entity of id vehid using finite differences on the velocity.
The first element is `missing`.
"""
function accft(roadway::Roadway, scenes::Vector{<:Scene}, vehid)
    result = zeros(Union{Missing,Float64},length(scenes))
    angles = broadcast(x->velf(get_by_id(x, vehid)).t, scenes)
    result[1] = missing
    result[2:end] .= diff(angles)
    return result
end

"""
    jerk(roadway::Roadway, scenes::Vector{<:Scene}, vehid)
returns the jerk history of the entity of id vehid using finite differences on acceleration (which uses finite differences on velocity).
The first two elements are missing.
"""
function jerk(roadway::Roadway, scenes::Vector{<:Scene}, vehid)
    return vcat(missing, diff(acc(roadway, scenes, vehid)))
end

"""
    jerkfs(roadway::Roadway, scenes::Vector{<:Scene}, vehid)
returns the longitudinal jerk history in the Frenet frame of the entity of id vehid using finite differences on acceleration (which uses finite differences on velocity).
The first two elements are missing.
"""
function jerkfs(roadway::Roadway, scenes::Vector{<:Scene}, vehid)
    return vcat(missing, diff(accfs(roadway, scenes, vehid)))
end

"""
    jerkft(roadway::Roadway, scenes::Vector{<:Scene}, vehid)
returns the lateral jerk history in the Frenet frame of the entity of id vehid using finite differences on acceleration (which uses finite differences on velocity).
The first two elements are missing.
"""
function jerkft(roadway::Roadway, scenes::Vector{<:Scene}, vehid)
    return vcat(missing, diff(accft(roadway, scenes, vehid)))
end

"""
    turn_rate_g(roadway::Roadway, scenes::Vector{<:Scene}, vehid)
returns the turn rate history in the Frenet frame of the entity of id vehid using finite differences on global heading.
The first element is missing.
"""
function turn_rate_g(roadway::Roadway, scenes::Vector{<:Scene}, vehid)
    result = zeros(Union{Missing,Float64},length(scenes))
    angles = broadcast(x->posg(get_by_id(x, vehid)).θ, scenes)
    result[1] = missing
    result[2:end] .= diff(angles)
    return result
end

"""
    turn_rate_g(roadway::Roadway, scenes::Vector{<:Scene}, vehid)
returns the turn rate history in the Frenet frame of the entity of id vehid using finite differences on heading in the Frenet frame.
The first element is missing.
"""
function turn_rate_f(roadway::Roadway, scenes::Vector{<:Scene}, vehid)
    result = zeros(Union{Missing,Float64},length(scenes))
    angles = broadcast(x->posf(get_by_id(x, vehid)).ϕ, scenes)
    result[1] = missing
    result[2:end] .= diff(angles)
    return result
end

"""
    isbraking(roadway::Roadway, scenes::Vector{<:Scene}, vehid)
history of braking events: true when the acceleration is negative.
The first element is missing.
"""
function isbraking(roadway::Roadway, scenes::Vector{<:Scene}, vehid)
    return acc(roadway, scenes, vehid) .< 0.0
end

"""
    isaccelerating(roadway::Roadway, scenes::Vector{<:Scene}, vehid)
history of acceleration events: true when the acceleration is positive
The first element is missing.
"""
function isaccelerating(roadway::Roadway, scenes::Vector{<:Scene}, vehid)
    return acc(roadway, scenes, vehid) .> 0.0
end


text:
function inverse(side::HorizonSide)
    return HorizonSide(mod(Int(side) + 2, 4))
end

function f()
    steps = []
    while !isborder(robot, Nord)
        move!(robot, Nord)
        push!(steps, inverse(Nord))
    end
    while !isborder(robot, Ost)
        move!(robot, Ost)
        push!(steps, inverse(Ost))
    end

    # Тут решение.
    a = 0
    if length(steps) % 2 == 0
        a -= 1
    end
    side = Ost
    while !isborder(robot, Sud)
        side = inverse(side)
        # putmarker!(robot)
        while !isborder(robot, side)
            move!(robot, side)
            a += 1
            if a % 2 == 1
                putmarker!(robot)
            end
        end
        move!(robot, Sud)
        a += 1
    end
    side = inverse(side)
    while !isborder(robot, side)
        move!(robot, side)
        a += 1
        if a % 2 == 1
            putmarker!(robot)
        end
    end

    while !isborder(robot, Nord)
        move!(robot, Nord)
    end
    while !isborder(robot, Ost)
        move!(robot, Ost)
    end
    while length(steps) > 0
        move!(robot, pop!(steps))
    end
end

f()

max_stars_repo_name: TheXlebushek/mirea-progs

max_issues_repo_name: TheXlebushek/mirea-progs

max_forks_repo_name: TheXlebushek/mirea-progs


text:
using LinearAlgebra, FFTW
import BSON, Statistics, Random
import PyPlot

import NMRSpectraSimulator

include("../src/NMRCalibrate.jl")
import .NMRCalibrate

# for loading something with Interpolations.jl
import OffsetArrays
import Interpolations

import PlotlyJS
import Plots
Plots.plotly()

#import Clustering

PyPlot.close("all")
fig_num = 1

Random.seed!(25)
PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])

include("./helpers/final_helpers.jl")

#projects_dir = save_folder_path
projects_dir = "/home/roy/MEGAsync/outputs/NMR/calibrate/final/D-(+)-Glucose-NRC-600"
#projects_dir = "/home/roy/MEGAsync/outputs/NMR/calibrate/final/D-(+)-Glucose-700"
#projects_dir = "/home/roy/MEGAsync/outputs/NMR/calibrate/final/Nam2022_Serine"
#projects_dir = "/home/roy/MEGAsync/outputs/NMR/calibrate/final/L-Serine-700"

### user inputs.
# projects_dir = "/home/roy/MEGAsync/outputs/NMR/calibrate/final"
#
# project_name = "D-(+)-Glucose-700"
# molecule_names = ["D-(+)-Glucose"; "DSS"]
#
# project_name = "D-(+)-Glucose-NRC-600"
# molecule_names = ["D-(+)-Glucose";]
#
# project_name = "L-Phenylalanine-700"
# molecule_names = ["L-Phenylalanine"; "DSS"]
#
# project_name = "L-Glutamine-700"
# molecule_names = ["L-Glutamine"; "DSS"]
#
# project_name = "D-(+)-Glucose-NRC-600"
# molecule_names = ["D-(+)-Glucose";]

### end user inputs.
plots_save_folder = joinpath(projects_dir, "plots")
isdir(plots_save_folder) || mkdir(plots_save_folder)
project_title = "test" # TODO change later.

### load block.
#load_path = joinpath(joinpath(projects_dir, project_name), "results_full.bson")
load_path = joinpath(projects_dir, "results_full.bson")
dict = BSON.load(load_path)
As = collect( dict[:As][i] for i = 1:length(dict[:As]) )
Δsys_cs = convert(Vector{Vector{Float64}}, dict[:Δsys_cs])
y = convert(Vector{Complex{Float64}}, dict[:y])
U_y = convert(Vector{Float64}, dict[:U_y])
SW = dict[:SW]
fs = dict[:fs]
ν_0ppm = dict[:ν_0ppm]

hz2ppmfunc = uu->(uu - ν_0ppm)*SW/fs
ppm2hzfunc = pp->(ν_0ppm + pp*fs/SW)
P_y = hz2ppmfunc.(U_y)

ΩS_ppm = NMRSpectraSimulator.getPsnospininfo(As, hz2ppmfunc)
ΩS_ppm_sorted = sort(NMRSpectraSimulator.combinevectors(ΩS_ppm))
u_offset = 0.5
u_min = ppm2hzfunc(ΩS_ppm_sorted[1] - u_offset)
u_max = ppm2hzfunc(ΩS_ppm_sorted[end] + u_offset)

P = LinRange(hz2ppmfunc(u_min), hz2ppmfunc(u_max), 50000)
U = ppm2hzfunc.(P)

function graphall(dict, Δ_shifts, Es, y, U_y, P_y, fs, SW::T,
    project_title::String,
    plots_save_folder::String;
    fig_num::Int = 1) where T <: Real

    p_star_set = dict[:p_star_set]
    κ_lb_default = dict[:κ_lb_default]
    κ_ub_default = dict[:κ_ub_default]
    κ_star_set = dict[:κ_star_set]
    d_star_set = dict[:d_star_set]
    β_star_set = dict[:β_star_set]
    λ_star_set = dict[:λ_star_set]
    cost_inds_set = dict[:cost_inds_set]
    w = dict[:w]

    As = collect( Es[n].core for n = 1:length(Es))

    for r = 1:length(cost_inds_set)

        y_cost = y[cost_inds_set[r]]
        U_cost = U_y[cost_inds_set[r]]
        P_cost = P_y[cost_inds_set[r]]

        LS_inds = 1:length(U_cost)

        q, updatedfunc, updateβfunc, updateλfunc, updateκfunc,
        κ_BLS, getshiftfunc, getβfunc, getλfunc,
        N_vars_set = NMRCalibrate.setupcostcLshiftLS(Es, As, fs, SW,
        LS_inds, U_cost, y_cost, Δ_shifts;
        w = w, κ_lb_default = κ_lb_default, κ_ub_default = κ_ub_default)

        obj_func = pp->NMRCalibrate.costcLshift(U_cost, y_cost,
        updatedfunc, updateβfunc, updateλfunc, updateκfunc, pp, Es, κ_BLS, q)

        #####

        ### end new.

        ### reference.

        # reference, zero shift, phase.
        N_d = sum( NMRCalibrate.getNd(As[n]) for n = 1:length(As) )
        #N_β = sum( NMRCalibrate.getNβ(As[n]) for n = 1:length(As) )
        N_λ = sum( NMRCalibrate.getNλ(As[n]) for n = 1:length(As) )
        shift_initial = zeros(T, N_d)
        β_initial = zeros(T, N_β)
        λ_initial = ones(T, N_λ)
        p_initial = [shift_initial; β_initial; λ_initial]

        initial_cost = obj_func(p_initial)
        NMRCalibrate.parseκ!(Es, ones(T, length(κ_BLS)))
        fill!(w, 1.0)
        q_initial_U = q.(U)
        #println("norm(q_initial_U) = ", norm(q_initial_U))

        p_test = copy(p_star_set[r])
        final_cost = obj_func(p_test)
        q_final_U = q.(U)

        PyPlot.figure(fig_num)
        fig_num += 1

        PyPlot.plot(P_y, real.(y), label = "y")
        PyPlot.plot(P_cost, real.(y_cost), "x")

        PyPlot.legend()
        PyPlot.xlabel("ppm")
        PyPlot.ylabel("real")
        PyPlot.title("r = $(r). cost vs. data (y)")

        PyPlot.figure(fig_num)
        fig_num += 1

        PyPlot.plot(P, real.(q_initial_U), label = "q initial")
        PyPlot.plot(P_y, real.(y), label = "y")
        PyPlot.plot(P, real.(q_final_U), "--", label = "q final")
        #PyPlot.plot(P_cost, real.(y_cost), "x")

        PyPlot.legend()
        PyPlot.xlabel("ppm")
        PyPlot.ylabel("real")
        PyPlot.title("r = $(r). data (y) vs. fit")

        plots_save_path = joinpath(plots_save_folder, "region_$(r)_real.html")
        title_string = "$(project_title): region $(r), real"
        savefigfitresult(plots_save_path, title_string,
         real.(q_final_U), P, P_cost, real.(y_cost);
        initial_fit = real.(q_initial_U))

        plots_save_path = joinpath(plots_save_folder, "region_$(r)_imaginary.html")
        title_string = "$(project_title): region $(r), imaginary"
        savefigfitresult(plots_save_path, title_string,
        imag.(q_final_U), P, P_cost, imag.(y_cost);
        initial_fit = imag.(q_initial_U))

        plots_save_path = joinpath(plots_save_folder, "region_$(r)_modulus.html")
        title_string = "$(project_title): region $(r), modulus"
        savefigfitresult(plots_save_path, title_string,
        abs.(q_final_U), P, P_cost, abs.(y_cost);
        initial_fit = abs.(q_initial_U))
    end

end

Δ_shifts = NMRSpectraSimulator.combinevectors(Δsys_cs)

Es = collect( NMRSpectraSimulator.κCompoundFIDType(As[i]) for i = 1:length(As) )
graphall(dict, Δ_shifts, Es, y, U_y, P_y, fs, SW,
project_title, plots_save_folder; fig_num = 1)


text:
using ALifeBenchmark
using Plots
using Measures

if length(ARGS) == 1
    trial_id = ARGS
else
    trial_id = "22578359168939802"
end

most_frequent_genotypes = []

snapshot_ids = get_snapshot_ids(trial_id)

for snapshot_id in snapshot_ids
    genotype_distribution = get_snapshot(snapshot_id) |> get_genotype_distribution

    current_frequent_genotypes = [g for (g, n) in genotype_distribution if n >= 0.01]
    append!(most_frequent_genotypes, current_frequent_genotypes)
end

most_frequent_genotypes = unique(most_frequent_genotypes)

println(length(most_frequent_genotypes))

"""
graph_data = load_graph_data()

nn_dict = Dict()
nns = []
for genotype in most_frequent_genotypes
    global cnt

    nn = 0
    for (i, neutral_network) in enumerate(graph_data.neutral_networks)
        if haskey(graph_data.genotype_vertex_mapping, genotype) && graph_data.genotype_vertex_mapping[genotype] in neutral_network
            nn = i
        end
    end

    if haskey(nn_dict, nn) == false
        nn_dict[nn] = length(nn_dict) + 1
    end

    push!(nns, nn_dict[nn])
end

@assert length(nns) == length(most_frequent_genotypes)

occurances = zeros(length(snapshot_ids), length(nn_dict))
times = zeros(length(snapshot_ids))

for (i, snapshot_id) in enumerate(snapshot_ids)
    genotype_distribution = get_snapshot(snapshot_id) |> get_genotype_distribution
    times[i] = get_snapshot(snapshot_id).time

    for (j, genotype) in enumerate(most_frequent_genotypes)
        if haskey(genotype_distribution, genotype)
            occurances[i, nns[j]] += genotype_distribution[genotype]
        end
    end
end

"""
occurances = zeros(length(snapshot_ids), length(most_frequent_genotypes))

for (i, snapshot_id) in enumerate(snapshot_ids)
    genotype_distribution = get_snapshot(snapshot_id) |> get_genotype_distribution

    for (j, genotype) in enumerate(most_frequent_genotypes)
        if haskey(genotype_distribution, genotype)
            occurances[i, j] += genotype_distribution[genotype]
        end
    end
end

plot(times, occurances,
        title = "",
        label = "",
        xguide = "Time",
        yguide = "Genotype Occurances",
        size = (900, 600),
        margin = 10mm,
        dpi = 1000)
savefig("genotype_distribution_$trial_id.png")

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
module RatPack
__precompile__(false)

import Base: copy, push!, append!, insert!, getindex, setindex!, length, size, lastindex, ==

using Random
using StatsBase
using Distributions
using CategoricalArrays
using DataFrames
using CSV
using LinearAlgebra

export update_ratings, predict_outcome, update_info, increment!, reset!, player_indexes, linesearch
export update_rule_list, update_rule_names, simulate_rule_list, simulate_rule_names, scoring_rule_list, scoring_rule_names
export UpdateRule, SimulateRule, GenerateRule, ScoringRule # parent abstract types -- instantiations are exported in their own files
export RatingsList, RatingsTable # main datastructures
export simulate, generate # simulation tools
export score_ratings, scoring_function, score_direction
export copy, push!, append!, insert!, getindex, setindex!, length, size, lastindex, ==
export cross_validate, optimise # higher-level wrapper functions

### Main data structures used here
include("DataStructures.jl")

### General purpose utilities
include("Utilities.jl")

### I/O routines
include("IO.jl")

###########################################################################3
### Models for update calculation rules and predictions
abstract type UpdateRule end

"""
    update_ratings()

 Compute one round of updates to ratings.

## Input arguments
* `rule::UpdateRule`: the type of update rule to use
* `input_ratings::RatingsList`: a list of current ratings to be used in recursive updates
* `input_competitions::DataFrame`: a DataFrame containing the outcomes of a series of competitions

```
"""
function update_ratings( rule::UpdateRule,
                         input_ratings::RatingsList,
                         input_competitions::DataFrame )
    error("undefined update rule") # this is the default that is called with an update rule that isn't properly defined
end

"""
    update_info()

 Return information about an update rule.

## Input arguments
* `rule::UpdateRule`: the type of update rule to describe

## Output
* `info::Dict(Symbol, Any)`: its fields describe the update rule
   + `:name => ` self reference for checking
   + `:reference =>` a reference to the source of the rule
   + `:computation =>` the mode in which computations are performed, options are
        - "simultaneous" = all input competitions are considered in one calculation
        - "sequential" = computations are made one batch of competitions at a time
   + `:state_model =>` how past ratings are used
        - "recursive" = new ratings are built on top of old ratings
        - "none" = new ratings ignore old ratings
   + `:input =>` the type of input competition data used
        - "none" = all input competitions are ignored
        - "outcome" = win/loss/tie listing for single competitions
        - "score" = the scores of the two players
        - "margin" = the margin between the two players
        - "result" = the totalled results of a set of competitions between two players
   + `:output =>`
        - "deterministic" model produces deterministic predictions
        - "probabilistic" model produces probabilistic predictions
        - "either" depends on sub-rule being used
   + `:model` => the type of model being used, at the moment only "single" is returned
   + `:ties` => true/false, does the model use information about ties appropriately
   + `:factors` => true/false, does the model all factors
   + `:parameters` => an array of strings describing the parameters of the update model (if any)

## Surprises
* `Elo` has `:computation => "simultaneous"` because this represents a single update round
   + sequential-Elo is performed by using the `Iterate` rule
* `Iterate` and `IterateSample `
   + take many qualities from a sub-rule that they iterate
   + only makes sense to use these for recursive sub-rules.
* Remember dictionaries are not sorted, so the order may change in output.

```
"""
function update_info( rule::UpdateRule )
    error("undefined update rule") # this is the default that is called with an update rule that isn't properly defined
end

"""
    predict_outcome()

 Predict outcome of a single match

## Input arguments
* `rule::UpdateRule`: the type of update rule to use
* `ratingA::Real`: rating of player A
* `ratingB::Real`: rating of player B
* `factorA::Union{Missing,Real}`: factors affecting player A where relevant (small int usually)
* `factorB::Union{Missing,Real}`: factors affecting player B where relevant (small int usually)

```
"""
function predict_outcome(rule::UpdateRule,
                         ratingA::Real, ratingB::Real,
                         factorA::Union{Missing,Real}, factorB::Union{Missing,Real} )
    if update_info(rule)[:input] != "outcome"
        error("the $rule update doesn't predict outcomes")
    else
        error("the $rule update doesn't have a prediction function")
    end
end

# actual instantiations of updates
update_rule_list = ["Colley", "Massey", "MasseyColley", "KeenerScores", "Revert", "EloF", "Elo", "Iterate", "SampleIterate"]
update_rule_names = Array{String,1}(undef,length(update_rule_list))
for (i,u) in enumerate(update_rule_list)
    u_file = "UpdateRules/$(u).jl"
    include(u_file)
    update_rule_names[i] = "Update$(u)"
end

###########################################################################3
### Simulation tools
abstract type SimulateRule end
abstract type GenerateRule end

"""
    simulate()

 Simulate a set of competition for rated players

## Input arguments
* `r::RatingsList`: the list of players and their ratings
* `model::SimulateRule`: the type of competition to simulate, and it's parameters
* `perf_model::ContinuousUnivariateDistribution`: the performance model (strength mapped to outcomes)

```
"""
function simulate(r::RatingsList, model::SimulateRule, perf_model::ContinuousUnivariateDistribution )
    error("undefined simulation rule") # this is the default that is called with rules that aren't properly defined
end

# actual instantiations of simulation rules
simulate_rule_list = ["RoundRobin","RoundRobinF",  "Elimination"]
simulate_rule_names = Array{String,1}(undef,length(simulate_rule_list))
for (i,u) in enumerate(simulate_rule_list)
    u_file = "Simulation/$(u).jl"
    include(u_file)
    simulate_rule_names[i] = "Sim$(u)"
end

"""
    generate()

 generate a set of rated players

## Input arguments
* `m::Int`: the number of players to simulate
* `model::???`: the model to use in generation

```
"""
function generate()
    error("undefined generation rule") # this is the default that is called with an update rule that isn't properly defined
end

###########################################################################3
### Scoring Rules
abstract type ScoringRule end

"""
    score()

 Compure "Scores" for a set of ratings and a model by comparing predictions from the ratings/model to contest data

## Input arguments
* `srule::ScoringRule`: the scoring rule to use
* `outcomes::DataFrame`: the outcomes of the contests (see I/O)
* `rule::UpdateRule`: the update/prediction rule
* `r::RatingsList`: a list of ratings (should be complete) of all players involved (typically derived using same rule)

```
"""
function score_ratings(srule::ScoringRule, outcomes::DataFrame, rule::UpdateRule, r::RatingsList)
    total_score = 0.0
    for i=1:size(outcomes,1)
        # form the predictions from the update rule's prediction function
        out_prob = predict_outcome(rule,
                                   r.ratings[outcomes[i,PlayerA]], r.ratings[outcomes[i,PlayerB]],
                                   outcomes[i,FactorA], outcomes[i,FactorB])

        # add scores over all contests
        o = outcomes[i,Outcome]
        if o==-1
            o = 2
        elseif o==0
            o = 3
        end
        total_score += scoring_function(srule, out_prob, o)
    end
    av_score = total_score/size(outcomes,1)

    return av_score
end

"""
    score_direction()

 Indicate whether score is positive (1) or negative (-1), i.e.,
* positive: means larger scores are better
* negative: means scores closer to zero are better

## Input arguments
* `srule::ScoringRule`: the scoring rule to use

```
"""
function score_direction(srule::ScoringRule)
    error("undefined generation rule") # this is the default that is called with scoring function that isn't properly defined
end

"""
    scoring_function()

 Compute "Scores" for predicted probabilities of outcome

## Input arguments
* `srule::ScoringRule`: the scoring rule to use
* `predicted_probabilities::Array{Float64,1}`: vector length C for C classes of outcome (must add to 1)
* `outcome::Int`: the outcome as an Integer from 1,...,C

```
"""
function scoring_function(srule::ScoringRule, predicted_probabilities::Array{Float64,1}, outcome::Int)
    error("undefined generation rule") # this is the default that is called with scoring function that isn't properly defined
end

function scoring_function(srule::ScoringRule, predicted_probabilities::Tuple, outcome::Int)
    return scoring_function(srule, [x for x in predicted_probabilities], outcome)
end

# actual instantiations of simulation rules
scoring_rule_list = ["Brier", "Logarithmic", "Quadratic", "Spherical"]
scoring_rule_names = Array{String,1}(undef,length(scoring_rule_list))
for (i,u) in enumerate(scoring_rule_list)
    u_file = "ScoringRules/$(u).jl"
    include(u_file)
    scoring_rule_names[i] = "Score$(u)"
end

###########################################################################3
### Wrapper functions
include("Wrappers.jl")

end # module


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

struct VirtualMachinesApi <: SwaggerApi
    client::Swagger.Client
end

"""

Run command on the VM.
Param: resourceGroupName::String (required)
Param: vmName::String (required)
Param: parameters::RunCommandInput (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: RunCommandResult
"""
function virtualMachinesRunCommand(_api::VirtualMachinesApi, resourceGroupName::String, vmName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", RunCommandResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/virtualMachines/{vmName}/runCommand", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "vmName", vmName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export virtualMachinesRunCommand


text:
module Architectures

import GPUifyLoops
using CUDAapi: has_cuda

export
    @hascuda,
    CPU, GPU,
    device, architecture, array_type

"""
    AbstractArchitecture

Abstract supertype for architectures supported by Oceananigans.
"""
abstract type AbstractArchitecture end

"""
    CPU <: AbstractArchitecture

Run Oceananigans on a single-core of a CPU.
"""
struct CPU <: AbstractArchitecture end

"""
    GPU <: AbstractArchitecture

Run Oceananigans on a single NVIDIA CUDA GPU.
"""
struct GPU <: AbstractArchitecture end

"""
    @hascuda expr

A macro to compile and execute `expr` only if CUDA is installed and available. Generally used to
wrap expressions that can only be compiled if `CuArrays` and `CUDAnative` can be loaded.
"""
macro hascuda(expr)
    return has_cuda() ? :($(esc(expr))) : :(nothing)
end

device(::CPU) = GPUifyLoops.CPU()
device(::GPU) = GPUifyLoops.CUDA()

@hascuda using CuArrays

         architecture(::Array)   = CPU()
@hascuda architecture(::CuArray) = GPU()

         array_type(::CPU) = Array
@hascuda array_type(::GPU) = CuArray

end


text:
using AlgebraicDynamics
using AlgebraicDynamics.UWDDynam
using AlgebraicDynamics.DWDDynam

using Catlab.WiringDiagrams
using Catlab.CategoricalAlgebra

using OrdinaryDiffEq

using Test
const UWD = UndirectedWiringDiagram

approx_equal(u, v) = abs(maximum(u - v)) < .1

tspan = (0.0, 100.0)
tspan_short = (0.0, 10.0)

@testset "Resource sharers" begin
  dx(x) = [1 - x[1]^2, 2*x[1]-x[2]]
  dy(y) = [1 - y[1]^2]
  r = ContinuousResourceSharer{Real}(2, (u,p,t) -> dx(u))

  u0 = [-1.0, -2.0]
  prob = ODEProblem(r, u0, tspan)
  @test prob isa ODEProblem
  sol = solve(prob, Tsit5())
  @test all(s -> s == u0, sol.u)

  u0 = [.1, 10.0]
  prob = ODEProblem(r, u0, tspan)
  sol = solve(prob, Tsit5())
  @test approx_equal(last(sol), [1.0, 2.0])

  h = 0.2
  dr = euler_approx(r, h)
  sol = trajectory(dr, u0, nothing, tspan) # nothing because no parameters
  @test approx_equal(last(sol), [1.0, 2.0])

  dr = DiscreteResourceSharer{Real}(1, (u,p,t) -> dy(u))
  u0 = [1.0]
  sol = trajectory(dr, u0, nothing, tspan)
  @test all(((i,s),) -> s[1] == i%2, enumerate(sol.u))
end

@testset "Machines" begin
  @testset "Euler approximation" begin
    uf(u, x, p, t) = [x[1] - u[1]*x[2]]
    rf(u,p,t) = u
    u0 = [1.0]
    m = ContinuousMachine{Float64}(2,1, 1, uf, rf)

    xs = [t -> 1, t -> 1]
    prob = ODEProblem(m, u0, xs, tspan)
    @test prob isa ODEProblem
    sol = solve(prob, Tsit5())
    @test all(s -> s == u0, sol.u)

    xs = t->t
    prob = ODEProblem(m, u0, xs, tspan)
    sol = solve(prob, Tsit5())
    @test all(s -> s == u0, sol.u)

    xs = [2, 1]
    prob = ODEProblem(m, u0, xs, tspan)
    sol = solve(prob, Tsit5())
    @test approx_equal(last(sol), [2.0])

    h = 0.2
    dm = euler_approx(m,h)
    sol = trajectory(dm, u0, xs, nothing, tspan)
    @test approx_equal(last(sol), [2.0])
  end

  @testset "oapply" begin
    uf(u, x, p, t) = [x[1] - u[1], 0.0]
    rf(u,p,t) = u
    mf = ContinuousMachine{Float64}(2,2,2, uf, rf)
    d_id = singleton_diagram(Box(:f, [:A, :A], [:A, :A]))
    m_id = oapply(d_id, [mf])

    xs = [0.0, 0.0]
    u0 = [10.0, 2.0]
    @testset "" for m in [mf, m_id]
      prob = ODEProblem(m, u0, xs, tspan)
      sol = solve(prob, Tsit5(); dtmax = 1)
      @test approx_equal(last(sol), [0.0, 2.0])
    end
  end
end

@testset "Lotka-Volterra model" begin
  params = [0.3, 0.015, 0.015, 0.7]
  u0 = [10.0, 100.0]

  bounds(sol::ODESolution) = all(s -> all([0, 0] .< s .< [200.0, 150.0]), sol.u)

  @testset "Resource sharer implementation" begin
    dotr(u,p,t) = p[1]*u
    dotrf(u,p,t) = [-p[2]*u[1]*u[2], p[3]*u[1]*u[2]]
    dotf(u,p,t) = -p[4]*u

    r = ContinuousResourceSharer{Real}(1, dotr)
    rf_pred = ContinuousResourceSharer{Real}(2, dotrf)
    f = ContinuousResourceSharer{Real}(1, dotf)

    rf_pattern = UWD(0)
    add_box!(rf_pattern, 1); add_box!(rf_pattern, 2); add_box!(rf_pattern, 1)
    add_junctions!(rf_pattern, 2)
    set_junction!(rf_pattern, [1,1,2,2])

    lv = oapply(rf_pattern, [r, rf_pred, f])

    prob = ODEProblem(lv, u0, tspan, params)
    sol = solve(prob, Tsit5())
    @test bounds(sol)

    h = 0.002
    lv_discrete = oapply(rf_pattern, euler_approx([r, rf_pred, f], h))
    sol2 = trajectory(lv_discrete, u0, params, tspan; dt = h)
    @test bounds(sol2)
  end

  @testset "Machine implementation" begin
    dotr(u, x, p, t) = [p[1]*u[1] - p[2]*u[1]*x[1]]
    dotf(u, x, p, t) = [p[3]*u[1]*x[1] - p[4]*u[1]]

    rmachine = ContinuousMachine{Real}(1,1,1, dotr, (r,p,t) -> r)
    fmachine = ContinuousMachine{Real}(1,1,1, dotf, (f,p,t) -> f)

    rf_pattern = WiringDiagram([],[])
    boxr = add_box!(rf_pattern, Box(nothing, [nothing], [nothing]))
    boxf = add_box!(rf_pattern, Box(nothing, [nothing], [nothing]))
    add_wires!(rf_pattern, Pair[
      (boxr, 1) => (boxf, 1),
      (boxf, 1) => (boxr, 1)
    ])

    rf_machine = oapply(rf_pattern, [rmachine, fmachine])
    prob = ODEProblem(rf_machine, u0, tspan, params)
    sol = solve(prob, Tsit5(); dtmax = .1)
    @test bounds(sol)

    h = 0.005
    lv_discrete = oapply(rf_pattern, euler_approx([rmachine, fmachine], h))
    sol = trajectory(lv_discrete, u0, params, tspan_short; dt = h)
    @test bounds(sol)
  end
end

@testset "Ocean model" begin
  α, β, γ, δ, β′, γ′, δ′ = 0.3, 0.015, 0.015, 0.7, 0.017, 0.017, 0.35
  params = [α, β, γ, δ, β′, γ′, δ′]
  u0 = [100.0, 10.0, 2.0]

  bounds(sol::ODESolution) = all(s -> all([0, 0, 0] .< s .< [200.0, 75, 20]), sol.u)

  dotfish(f, x, p, t) = [p[1]*f[1] - p[2]*x[1]*f[1]]
  dotFISH(F, x, p, t) = [p[3]*x[1]*F[1] - p[4]*F[1] - p[5]*x[2]*F[1]]
  dotsharks(s, x, p, t) = [-p[7]*s[1] + p[6]*s[1]*x[1]]

  fish   = ContinuousMachine{Real}(1,1,1, dotfish,   (f,p,t)->f)
  FISH   = ContinuousMachine{Real}(2,1,2, dotFISH,   (F,p,t)->[F[1], F[1]])
  sharks = ContinuousMachine{Real}(1,1,1, dotsharks, (s,p,t)->s)

  ocean_pat = WiringDiagram([], [])
  boxf = add_box!(ocean_pat, Box(nothing, [nothing], [nothing]))
  boxF = add_box!(ocean_pat, Box(nothing, [nothing, nothing], [nothing, nothing]))
  boxs = add_box!(ocean_pat, Box([nothing], [nothing]))
  add_wires!(ocean_pat, Pair[
    (boxf, 1) => (boxF, 1),
    (boxs, 1) => (boxF, 2),
    (boxF, 1) => (boxf, 1),
    (boxF, 2) => (boxs, 1)
  ])

  ocean = oapply(ocean_pat, [fish, FISH, sharks])
  prob = ODEProblem(ocean, u0, tspan, params)
  sol = solve(prob, Tsit5(); dtmax = 0.1)
  @test bounds(sol)

  h = 0.01
  sol = trajectory(euler_approx(ocean, h), u0, params, tspan_short; dt = h)
  @test bounds(sol)
end


text:
import Dates

let sundayCount = 0
    date = Dates.Date(1901, 1, 1)
    endDate = Dates.Date(2000, 12, 31)
    while date <= endDate
        sundayCount += Dates.dayofweek(date) == Dates.Sunday
        date += Dates.Month(1)
    end
    println(sundayCount)
end

max_stars_repo_name: seanbreckenridge/poly-project-euler

max_issues_repo_name: seanbreckenridge/poly-project-euler

max_forks_repo_name: seanbreckenridge/poly-project-euler


text:
using Pkg; Pkg.develop(path=dirname(@__DIR__))
# headless GK to fix ci
ENV["GKSwstype"] = "100"

using Documenter
using Literate

# precompile stuff now so the output won't show up in the docs
using BlockSystems
using Plots
using OrdinaryDiffEq
using ModelingToolkit
using Graphs
using NetworkDynamics

# generate examples
examples = [
    joinpath(@__DIR__, "..", "examples", "spacecraft.jl"),
    joinpath(@__DIR__, "..", "examples", "kuramoto_network.jl"),
    joinpath(@__DIR__, "..", "examples", "kuramoto_without_nd.jl"),
    # joinpath(@__DIR__, "..", "examples", "pd_node.jl"),
]
OUTPUT = joinpath(@__DIR__, "src/generated")
isdir(OUTPUT) && rm(OUTPUT, recursive=true)
mkpath(OUTPUT)

for ex in examples
    Literate.markdown(ex, OUTPUT)
    Literate.script(ex, OUTPUT)
end

makedocs(;
    modules=[BlockSystems],
    authors="Hans Würfel <git@wuerfel.io> and contributors",
    repo="https://github.com/hexaeder/BlockSystems.jl/blob/{commit}{path}#L{line}",
    sitename="BlockSystems.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://hexaeder.github.io/BlockSystems.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Examples" => ["Spacecraft" => "generated/spacecraft.md",
                       # "PowerDynamics.jl Node" => "generated/pd_node.md",
                       "Kuramoto Network" => "generated/kuramoto_network.md",
                       "Kuramoto without ND.jl" => "generated/kuramoto_without_nd.md"]
    ],
)

deploydocs(;
    repo="github.com/hexaeder/BlockSystems.jl",
    devbranch="main",
    push_preview=true,
)


text:
@testset "NumField/NfAbs" begin
  include("NfAbs/ConjugatesNS.jl")
  include("NfAbs/MPolyGcd.jl")
  include("NfAbs/MPolyAbsFact.jl")
  include("NfAbs/NfAbs.jl")
end


text:
@testset "Constructing empty sequences" begin
    @test DNASequence() == BioSequence(DNA)
    @test RNASequence() == BioSequence(RNA)
    @test AminoAcidSequence() == BioSequence(AminoAcid)
    @test CharSequence() == BioSequence(Char)
end

@testset "Constructing uninitialized sequences" begin
    @test isa(BioSequence{DNAAlphabet{2}}(0), BioSequence)
    @test isa(BioSequence{DNAAlphabet{4}}(10), BioSequence)
    @test isa(BioSequence{RNAAlphabet{2}}(0), BioSequence)
    @test isa(BioSequence{RNAAlphabet{4}}(10), BioSequence)
    @test isa(BioSequence{AminoAcidAlphabet}(10), BioSequence)
end

@testset "Conversion from/to strings" begin
    # Check that sequences in strings survive round trip conversion:
    #   String → BioSequence → String
    function test_string_construction(A::Type, seq::AbstractString)
        @test convert(String, BioSequence{A}(seq)) == uppercase(seq)
    end

    function test_string_parse(A::Type, seq::AbstractString)
        @test parse(BioSequence{A}, seq) == BioSequence{A}(seq)
    end

    for len in [0, 1, 2, 3, 10, 32, 1000, 10000]
        test_string_construction(DNAAlphabet{4}, random_dna(len))
        test_string_construction(DNAAlphabet{4}, lowercase(random_dna(len)))
        test_string_construction(RNAAlphabet{4}, lowercase(random_rna(len)))
        test_string_construction(RNAAlphabet{4}, random_rna(len))
        test_string_construction(AminoAcidAlphabet, random_aa(len))
        test_string_construction(AminoAcidAlphabet, lowercase(random_aa(len)))

        test_string_parse(DNAAlphabet{4}, random_dna(len))
        test_string_parse(DNAAlphabet{4}, lowercase(random_dna(len)))
        test_string_parse(RNAAlphabet{4}, lowercase(random_rna(len)))
        test_string_parse(RNAAlphabet{4}, random_rna(len))
        test_string_parse(AminoAcidAlphabet, random_aa(len))
        test_string_parse(AminoAcidAlphabet, lowercase(random_aa(len)))

        probs = [0.25, 0.25, 0.25, 0.25, 0.00]
        test_string_construction(DNAAlphabet{2}, random_dna(len, probs))
        test_string_construction(DNAAlphabet{2}, lowercase(random_dna(len, probs)))
        test_string_construction(RNAAlphabet{2}, random_rna(len, probs))
        test_string_construction(RNAAlphabet{2}, lowercase(random_rna(len, probs)))

        test_string_parse(DNAAlphabet{2}, random_dna(len, probs))
        test_string_parse(DNAAlphabet{2}, lowercase(random_dna(len, probs)))
        test_string_parse(RNAAlphabet{2}, random_rna(len, probs))
        test_string_parse(RNAAlphabet{2}, lowercase(random_rna(len, probs)))
    end

    # non-standard string literal
    @test isa(dna"ACGTMRWSYKVHDBN-", DNASequence)
    @test isa(rna"ACGUMRWSYKVHDBN-", RNASequence)
    @test isa(aa"ARNDCQEGHILKMFPSTWYVBJZXOU*-", AminoAcidSequence)
    @test isa(char"いろは αβγ 甲乙丙", CharSequence)

    # Non-nucleotide characters should throw
    @test_throws Exception DNASequence("ACCNNCATTTTTTAGATXATAG")
    @test_throws Exception RNASequence("ACCNNCATTTTTTAGATXATAG")
    @test_throws Exception AminoAcidSequence("ATGHLMY@ZACAGNM")
end

@testset "Construction from vectors" begin
    function test_vector_construction(A, seq::AbstractString)
        T = eltype(A)
        xs = T[convert(T, c) for c in seq]
        @test BioSequence{A}(xs) == BioSequence{A}(seq)
    end

    for len in [0, 1, 10, 32, 1000, 10000]
        test_vector_construction(DNAAlphabet{4}, random_dna(len))
        test_vector_construction(RNAAlphabet{4}, random_rna(len))
        test_vector_construction(AminoAcidAlphabet, random_aa(len))

        probs = [0.25, 0.25, 0.25, 0.25, 0.00]
        test_vector_construction(DNAAlphabet{2}, random_dna(len, probs))
        test_vector_construction(RNAAlphabet{2}, random_rna(len, probs))
    end
end

@testset "Conversion between 2-bit and 4-bit encodings" begin
    function test_conversion(A1, A2, seq)
        @test convert(BioSequence{A1}, BioSequence{A2}(seq)) == convert(BioSequence{A1}, seq)
    end

    test_conversion(DNAAlphabet{2}, DNAAlphabet{4}, "")
    test_conversion(DNAAlphabet{4}, DNAAlphabet{2}, "")
    test_conversion(RNAAlphabet{4}, RNAAlphabet{2}, "")
    test_conversion(RNAAlphabet{2}, RNAAlphabet{4}, "")

    test_conversion(DNAAlphabet{2}, DNAAlphabet{4}, "ACGT")
    test_conversion(DNAAlphabet{4}, DNAAlphabet{2}, "ACGT")
    test_conversion(RNAAlphabet{4}, RNAAlphabet{2}, "ACGU")
    test_conversion(RNAAlphabet{2}, RNAAlphabet{4}, "ACGU")

    test_conversion(DNAAlphabet{2}, DNAAlphabet{4}, "ACGT"^100)
    test_conversion(DNAAlphabet{4}, DNAAlphabet{2}, "ACGT"^100)
    test_conversion(RNAAlphabet{4}, RNAAlphabet{2}, "ACGU"^100)
    test_conversion(RNAAlphabet{2}, RNAAlphabet{4}, "ACGU"^100)

    # ambiguous nucleotides cannot be stored in 2-bit encoding
    EncodeError = BioSequences.EncodeError
    @test_throws EncodeError convert(BioSequence{DNAAlphabet{2}}, dna"AN")
    @test_throws EncodeError convert(BioSequence{RNAAlphabet{2}}, rna"AN")

    # test promotion
    a = BioSequence{DNAAlphabet{2}}("ATCG")
    b = BioSequence{DNAAlphabet{4}}("ATCG")
    c = BioSequence{RNAAlphabet{2}}("AUCG")
    d = BioSequence{RNAAlphabet{4}}("AUCG")

    @test typeof(promote(a, b)) == Tuple{BioSequence{DNAAlphabet{4}},BioSequence{DNAAlphabet{4}}}
    @test typeof(promote(c, d)) == Tuple{BioSequence{RNAAlphabet{4}},BioSequence{RNAAlphabet{4}}}
    @test_throws ErrorException typeof(promote(a, d))
    @test_throws ErrorException typeof(promote(a, b, d))
end

@testset "Conversion between RNA and DNA" begin
    @test convert(RNASequence, DNASequence("ACGTN")) == rna"ACGUN"
    @test convert(DNASequence, RNASequence("ACGUN")) == dna"ACGTN"
end

@testset "Conversion to Matrices" begin
    dna = [dna"AAA", dna"TTT", dna"CCC", dna"GGG"]
    dnathrow = [dna"AAA", dna"TTTAAA", dna"CCC", dna"GGG"]

    rna = [rna"AAA", rna"UUU", rna"CCC", rna"GGG"]
    rnathrow = [rna"AAA", rna"UUU", rna"CCCUUU", rna"GGG"]

    prot = [aa"AMG", aa"AMG", aa"AMG", aa"AMG"]

    sitemajdna = [
        DNA_A  DNA_A  DNA_A
        DNA_T  DNA_T  DNA_T
        DNA_C  DNA_C  DNA_C
        DNA_G  DNA_G  DNA_G
    ]
    seqmajdna = [
        DNA_A  DNA_T  DNA_C  DNA_G
        DNA_A  DNA_T  DNA_C  DNA_G
        DNA_A  DNA_T  DNA_C  DNA_G
    ]
    sitemajrna = [
        RNA_A  RNA_A  RNA_A
        RNA_U  RNA_U  RNA_U
        RNA_C  RNA_C  RNA_C
        RNA_G  RNA_G  RNA_G
    ]
    seqmajrna = [
        RNA_A  RNA_U  RNA_C  RNA_G
        RNA_A  RNA_U  RNA_C  RNA_G
        RNA_A  RNA_U  RNA_C  RNA_G
    ]
    sitemajnucint = [
        0x01  0x01  0x01
        0x08  0x08  0x08
        0x02  0x02  0x02
        0x04  0x04  0x04
    ]
    seqmajnucint = [
        0x01  0x08  0x02  0x04
        0x01  0x08  0x02  0x04
        0x01  0x08  0x02  0x04
    ]
    sitemajaa = [
        AA_A AA_M AA_G
        AA_A AA_M AA_G
        AA_A AA_M AA_G
        AA_A AA_M AA_G
    ]
    seqmajaa = [
        AA_A AA_A AA_A AA_A
        AA_M AA_M AA_M AA_M
        AA_G AA_G AA_G AA_G
    ]

    @test seqmatrix(dna, :site) == sitemajdna
    @test seqmatrix(rna, :site) == sitemajrna
    @test seqmatrix(prot, :site) == sitemajaa
    @test seqmatrix(UInt8, dna, :site) == sitemajnucint
    @test seqmatrix(UInt8, rna, :site) == sitemajnucint

    @test seqmatrix(dna, :seq) == seqmajdna
    @test seqmatrix(rna, :seq) == seqmajrna
    @test seqmatrix(prot, :seq) == seqmajaa
    @test seqmatrix(UInt8, dna, :seq) == seqmajnucint
    @test seqmatrix(UInt8, rna, :seq) == seqmajnucint

    @test seqmatrix([dna"", dna"", dna""], :site) == Matrix{DNA}(undef, (3, 0))
    @test seqmatrix([dna"", dna"", dna""], :seq) == Matrix{DNA}(undef, (0, 3))
    @test seqmatrix([rna"", rna"", rna""], :site) == Matrix{RNA}(undef, (3, 0))
    @test seqmatrix([rna"", rna"", rna""], :seq) == Matrix{RNA}(undef, (0, 3))
    @test seqmatrix(UInt8, [dna"", dna"", dna""], :site) == Matrix{UInt8}(undef, (3, 0))
    @test seqmatrix(UInt8, [dna"", dna"", dna""], :seq) == Matrix{UInt8}(undef, (0, 3))
    @test seqmatrix(UInt8, [rna"", rna"", rna""], :site) == Matrix{UInt8}(undef, (3, 0))
    @test seqmatrix(UInt8, [rna"", rna"", rna""], :seq) == Matrix{UInt8}(undef, (0, 3))

    @test_throws ArgumentError seqmatrix(dnathrow, :site)
    @test_throws ArgumentError seqmatrix(rnathrow, :seq)
    @test_throws ArgumentError seqmatrix(dna, :lol)
    @test_throws MethodError seqmatrix(AminoAcid, dna, :site)
    @test_throws ArgumentError seqmatrix(DNASequence[], :site)
    @test_throws ArgumentError seqmatrix(DNASequence[], :seq)

end


text:
using DelimitedFiles

input = readdlm(ARGS[1], ' ', Int, '\n');

# Part 1
fuel_required(mass) = (mass ÷ 3) - 2;
part1 = reduce(+, map(fuel_required, input))
println("Part 1: ", part1)

# Part 2
function fuel_required_2(mass)
    fuel = fuel_required(mass)
    return (fuel > 0) ? fuel + fuel_required_2(fuel) : 0;
end
part2 = reduce(+, map(fuel_required_2, input))
println("Part 2: ", part2)

# Nothing too crazy here for those familiar with recursion, although reading the instructions carefully always helps.
# Notice something interesting: the part 2 function does not distribute. That is, you cannot compute
#   fuel_required_2(part1) + part1
# to get part2.
# https://twitter.com/wjholdentech/status/1201023263633219585


text:
function install_envs()
  cd("gym")
  if haskey(ENV, "GYM_ENVS")
    envs = split(ENV["GYM_ENVS"], ":")
    for env in envs
      run(`pip install --user -e .[$env]`)
      @info("$env environment has been installed")
    end
  else
    run(`pip install --user -e .`)
    @info("Minimal installation has been performed")
    @info("To install everything, set ENV[\"GYM_ENVS\"]=all")
    @info("Then build the package again, Pkg.build(\"Gym\")")
  end
  cd("..")
end

function which_gym()
  whichgym = ""
  filename = "whichgym"
  if isfile(filename)
    f = open(whichgym)
    whichgym = strip(readstring(f))
    close(f)
  else
    whichgym = "SYSTEM"
  end
  return whichgym
end

try
  using PyCall
  PyCall.pyimport("gym")
  if which_gym() == "SYSTEM"
    @warn("You are using the gym from the systemwise installation.")
    @warn("If you want to install new environments, you should install them using \"pip install -e .[env]\"")
  else
    install_envs()
  end
catch
  if !isdir("gym")
    @info("Downloading OpenAi gym")
    run(`git clone https://github.com/openai/gym.git`)
  end
  @info("Installing OpenAi gym")
  install_envs()
  filename = "whichgym"
  f = open(filename, "w")
  write(f, "OWN")
  close(f)
end


text:
using PlotlyJS
using DataFrames, RDatasets

function linescatter()
    trace1 = scatter(x=1:4, y=[10, 15, 13, 18])
    plot(trace1)
end
linescatter()

function multiple_scatter_traces()
    trace1 = scatter(;x=1:4, y=[10, 15, 13, 17], mode="markers", name="marker only")
    trace2 = scatter(;x=1:4, y=[16, 5, 11, 9], mode="lines", name="line")
    trace3 = scatter(;x=1:4, y=[12, 9, 15, 12], mode="lines+markers", name="line+marker")
    trace4 = scatter(;x=1:4, y=[5, 10, 8, 12], mode="lines", line_dash="dash", name="dash")
    plot([trace1, trace2, trace3, trace4])
end
multiple_scatter_traces()

function data_labels()
    trace1 = scatter(;x=1:5.5, y=[1, 6, 3, 6, 1],
                      mode="markers+text", name="Team A",
                      textposition="top center",
                      text=["A-1", "A-2", "A-3", "A-4", "A-5"],
                      marker_size=12, textfont_family="Raleway, sans-serif")

    trace2 = scatter(;x=1:5, y=[4, 1, 7, 1, 4],
                      mode="markers+text", name= "Team B",
                      textposition="bottom center",
                      text= ["B-a", "B-b", "B-c", "B-d", "B-e"],
                      marker_size=12, textfont_family="Times New Roman")

    data = [trace1, trace2]

    layout = Layout(;title="Data Labels on the Plot", xaxis_range=[0.75, 5.25],
                     yaxis_range=[0, 8], legend_y=0.5, legend_yref="paper",
                     legend=attr(family="Arial, sans-serif", size=20,
                                 color="grey"))
    plot(data, layout)
end
data_labels()

function area1()
    trace1 = scatter(;x=1:4, y=[0, 2, 3, 5], fill="tozeroy")
    trace2 = scatter(;x=1:4, y=[3, 5, 1, 7], fill="tonexty")
    plot([trace1, trace2])
end
area1()

function matrix_subplots()
    z = rand(10,10)
    trace0 = scatter(; y=z)

    trace1 = heatmap(; z=z, showscale=false)
    trace2 = contour(; z=z)
    trace3 = surface(; z=z)
    p = [plot(trace1) plot(trace2)]
end

matrix_subplots()

function advanced_layouts()
    trace1 = scatter(;y=rand(10), mode="markers")
    trace2 = bar(;y=rand(10), xaxis="x2", yaxis="y2")
    trace3 = scatter(;y=rand(10), xaxis="x3", yaxis="y3")

    data = [trace1, trace2, trace3]
    xdomains = [[0,0.3], [0.33, 0.53] , [0.56, 0.78], [0.8, 1] ]
    ydomains = [[0,1],   [0, 0.27] , [0.33, 0.63], [0.66, 1] ]
    layout = Layout(; xaxis_domain=xdomains[1],
                    xaxis2_domain=xdomains[2], yaxis2_domain=ydomains[2],yaxis2_anchor="x2", xaxis2_anchor="y2",
                    xaxis3_domain=xdomains[2], yaxis3_domain=ydomains[3],yaxis3_anchor="x3", xaxis3_anchor="y3"
                    )
    plot(data, layout)
end
p1 = advanced_layouts()

function grouped_bar_example()
    trace1 = bar(x=["giraffes", "orangutans", "monkeys"],
                  y=[20, 14, 23],
                  name="SF Zoo")
    trace2 = bar(x=["giraffes", "orangutans", "monkeys"],
                  y=[12, 18, 29],
                  name="LA Zoo")
    data = [trace1, trace2]
    layout = Layout(;barmode="group")
    plot(data, layout)
end

grouped_bar_example()

function stacked_bar_example()
    trace1 = bar(;x=["giraffes", "orangutans", "monkeys"],
                  y=[20, 14, 23],
                  name="SF Zoo")
    trace2 = bar(x=["giraffes", "orangutans", "monkeys"],
                 y=[12, 18, 29],
                 name="LA Zoo")
    data = [trace1, trace2]
    layout = Layout(;barmode="stack")
    plot(data, layout)
end
stacked_bar_example()

function two_hists()
    x0 = randn(500)
    x1 = x0+1

    trace1 = histogram(x=x0, opacity=0.75)
    trace2 = histogram(x=x1, opacity=0.5)
    data = [trace1, trace2]
    layout = Layout(barmode="overlay")
    plot(data, layout)
end

two_hists()

function box_plot()
    x0 = ["day 1", "day 1", "day 1", "day 1", "day 1", "day 1",
          "day 2", "day 2", "day 2", "day 2", "day 2", "day 2"]
    trace1 = box(;y=[0.2, 0.2, 0.6, 1.0, 0.5, 0.4, 0.2, 0.7, 0.9, 0.1, 0.5, 0.3],
                  x=x0,
                  name="kale",
                  marker_color="#3D9970")
    trace2 = box(;y=[0.6, 0.7, 0.3, 0.6, 0.0, 0.5, 0.7, 0.9, 0.5, 0.8, 0.7, 0.2],
                  x=x0,
                  name="radishes",
                  marker_color="#FF4136")
    trace3 = box(;y=[0.1, 0.3, 0.1, 0.9, 0.6, 0.6, 0.9, 1.0, 0.3, 0.6, 0.8, 0.5],
                  x=x0,
                  name="carrots",
                  marker_color="#FF851B")
    data = [trace1, trace2, trace3]
    layout = Layout(;yaxis=attr(title="normalized moisture", zeroline=false),
                    boxmode="group")
    plot(data, layout)
end
box_plot()

iris = dataset("datasets", "iris");
display(head(iris))

function data_frame_scatter()
    iris = dataset("datasets", "iris");
#     display(head(iris))
#     my_trace = scatter(iris, x=:SepalLength, y=:SepalWidth, mode="markers", group=:Species)
#     plot(my_trace)
    p = Plot(iris, x=:SepalLength, y=:SepalWidth, mode="markers", marker_size=8, group=:Species)
    _p = JupyterPlot(p) #In Atom use _p = ElectronPlot(p)
    display(_p)
end

data_frame_scatter()


text:
# # Sphere sound-soft scattering

# In this example we solve the scattering by a sphere. First, we need to load
# the modules
using LinearAlgebra
using WaveProp
using WaveProp.Nystrom
using WaveProp.Geometry
using WriteVTK

# We now specify which *pde* we want to solve, as well as the geometry:
dim = 3; k=2π
pde = Helmholtz(dim=dim,k=k)
Ω,mesh   = WaveProp.IO.gmsh_sphere(dim=2,h=0.1)
WaveProp.Mesh.compute_quadrature!(mesh,dim=2,order=1,need_normal=true)

# We are interested in solving a sound-soft scattering problem. We will take an
# incident field `uᵢ` as a plane wave, and solve for the scattered field using a
# combined field integral equation:
uᵢ       = (x)   -> exp(im*k*dot((1,0,0),x))
rhs      = γ₀((x)-> -uᵢ(x),mesh)
S,D      = single_double_layer(pde,mesh)
L        = I/2 + D - im*k*S
ϕ        = L\rhs

# Now that the density ϕ has been found, we can represent the solution
𝒮        = SingleLayerPotential(pde,mesh)
𝒟        = DoubleLayerPotential(pde,mesh)
u(x)     = uᵢ(x) + 𝒟[ϕ](x) - im*k*𝒮[ϕ](x)

# and output it on a plane for visualization

## output the mesh
fname = joinpath(@__DIR__,"mesh")
vtkfile = WaveProp.IO.vtk_mesh_file(mesh,Ω,fname) |> vtk_save

## compute the solution on a plane
rect,out_mesh  =  WaveProp.IO.gmsh_rectangle(origin=(-5,-5,0),dx=10,dy=10,h=0.1,dim=2)
vals           = [u(x) for x in WaveProp.Mesh.nodes(out_mesh)]
fname = joinpath(@__DIR__,"output")
vtkfile = WaveProp.IO.vtk_mesh_file(out_mesh,rect,fname)
vtkfile["total field (real part)"] = real(vals)
vtkfile["total field (imag part)"] = imag(vals)
vtk_save(vtkfile)

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
const svg_rgb = Dict(
    "aliceblue" => "#F0F8FF",
    "antiquewhite" => "#FAEBD7",
    "aqua" => "#00FFFF",
    "aquamarine" => "#7FFFD4",
    "azure" => "#F0FFFF",
    "beige" => "#F5F5DC",
    "bisque" => "#FFEBCD",
    "black" => "#000000",
    "blanchedalmond" => "#FFEBCD",
    "blue" => "#0000FF",
    "blueviolet" => "#8A2BE2",
    "brown" => "#A52A2A",
    "burlywood" => "#DEB887",
    "cadetblue" => "#59EA0",
    "chartreuse" => "#7FFF00",
    "chocolate" => "#D2691E",
    "coral" => "#FF7F50",
    "cornflowerblue" => "#6495ED",
    "cornsilk" => "#FFF8DC",
    "crimson" => "#DC143C",
    "cyan" => "#00FFFF",
    "darkblue" => "#00008B",
    "darkcyan" => "#008B8B",
    "darkgoldenrod" => "#B8860B",
    "darkgray" => "#A9A9A9",
    "darkgreen" => "#006400",
    #"darkgrey" => "#A9A9A9",
    "darkkhaki" => "#BDB76B",
    "darkmagenta" => "#8B008B",
    "darkolivegreen" => "#556B2F",
    "darkorange" => "#FF8C00",
    "darkorchid" => "#9932CC",
    "darkred" => "#8B0000",
    "darksalmon" => "#E9967A",
    "darkseagreen" => "#8FBC8F",
    "darkslateblue" => "#483D8B",
    "darkslategray" => "#2F4F4F",
    #"darkslategrey" => "#2F4F4F",
    "darkturquoise" => "#00CED1",
    "darkviolet" => "#9400D3",
    "deeppink" => "#FF1493",
    "deepskyblue" => "#00BFFF",
    "dimgray" => "#696969",
    #"dimgrey" => "#696969",
    "dodgerblue" => "#1E90FF",
    "firebrick" => "#B22222",
    "floralwhite" => "#FFFAF0",
    "forestgreen" => "#22EB22",
    "fuchsia" => "#FF00FF",
    "gainsboro" => "#DCDCDC",
    "ghostwhite" => "#F8F8FF",
    "gold" => "#FFD700",
    "goldenrod" => "#DAA520",
    "gray" => "#808080",
    #"grey" => "#808080",
    "green" => "#008000",
    "greenyellow" => "#ADFF2F",
    "honeydew" => "#F0FFF0",
    "hotpink" => "#FF69B4",
    "indianred" => "#CD5C5C",
    "indigo" => "#4B0082",
    "ivory" => "#FFFFF0",
    "khaki" => "#F0E68C",
    "lavender" => "#E6E6FA",
    "lavenderblush" => "#FFF0F5",
    "lawngreen" => "#7CFC00",
    "lemonchiffon" => "#FFFACD",
    "lightblue" => "#ADD8E6",
    "lightcoral" => "#F08080",
    "lightcyan" => "#E0FFFF",
    "lightgoldenrodyellow" => "#FAFAD2",
    "lightgray" => "#D3D3D3",
    #"lightgrey" => "#D3D3D3",
    "lightgreen" => "#90EE90",
    "lightpink" => "#FFB6C1",
    "lightsalmon" => "#FFA07A",
    "lightseagreen" => "#20B2AA",
    "lightskyblue" => "#87CEFA",
    "lightslategray" => "#778899",
    #"lightslategrey" => "#778899",
    "lightsteelblue" => "#B0C4DE",
    "lightyellow" => "#FFFFE0",
    "lime" => "#00FF00",
    "limegreen" => "#32CD32",
    "linen" => "#FAF0E6",
    "magenta" => "#FF00FF",
    "maroon" => "#800000",
    "mediumaquamarine" => "#66CDAA",
    "mediumblue" => "#0000CD",
    "mediumorchid" => "#BA55D3",
    "mediumpurple" => "#9370D8",
    "mediumseagreen" => "#3CB371",
    "mediumslateblue" => "#7B68EE",
    "mediumspringgreen" => "#00FA9A",
    "mediumturquoise" => "#48D1CC",
    "mediumvioletred" => "#C71585",
    "midnightblue" => "#191970",
    "mintcream" => "#F5FFFA",
    "mistyrose" => "#FFE4E1",
    "moccasin" => "#FFE4B5",
    "navajowhite" => "#FFDEAD",
    "navy" => "#000080",
    "oldlace" => "#FDF5E6",
    "olive" => "#808000",
    "olivedrab" => "#6B8E23",
    "orange" => "#FFA500",
    "orangered" => "#FF4500",
    "orchid" => "#DA70D6",
    "palegoldenrod" => "#EEE8AA",
    "palegreen" => "#98FB98",
    "paleturquoise" => "#AFEEEE",
    "palevioletred" => "#D87093",
    "papayawhip" => "#FFEFD5",
    "peachpuff" => "#FFDAB9",
    "peru" => "#CD853F",
    "pink" => "#FFC0CB",
    "plum" => "#DDA0DD",
    "powderblue" => "#B0E0E6",
    "purple" => "#800080",
    "red" => "#FF0000",
    "rosybrown" => "#BC8F8F",
    "royalblue" => "#4169E1",
    "saddlebrown" => "#8B4513",
    "salmon" => "#FA8072",
    "sandybrown" => "#F4A460",
    "seagrean" => "#2E8B57",
    "seashell" => "#FFF5EE",
    "sienna" => "#A0522D",
    "silver" => "#C0C0C0",
    "skyblue" => "#87CEEB",
    "slateblue" => "#6A5ACD",
    "slategray" => "#708090",
    #"slategrey" => "#708090",
    "snow" => "#FFFAFA",
    "springgreen" => "#00FF7F",
    "steelblue" => "#4682B4",
    "tan" => "#D2B48C",
    "teal" => "#008080",
    "thistle" => "#D8BFD8",
    "tomato" => "#FF6347",
    "turquoise" => "#40E0D0",
    "violet" => "#EE82EE",
    "wheat" => "#F5DEB3",
    "white" => "#FFFFFF",
    "whitesmoke" => "#F5F5F5",
    "yellow" => "#FFFF00",
    "yellowgreen" => "#9ACD32"
)

max_stars_repo_name: gabrielgellner/OdinSon

max_issues_repo_name: gabrielgellner/OdinSon

max_forks_repo_name: gabrielgellner/OdinSon


text:
using Distributions
using Random
using StaticArrays

struct Kronecker_sampler{T, K} <: Random.Sampler{NTuple{K, T}} #TODO make edge type abstract
    power::Int
    dist::AliasTable{NTuple{K, T}, Vector{NTuple{K, T}}}
    initializer_size::T
    remainder::T
end

"""
    Kronecker_sampler(initializer, power; T=Int, space=1)

Create an efficient edge sampler whose probability distribution follows the specified
kroneker model.

Giving the sampler an adequite (but not excessive) amount of space may speed up generation
by up to 2x. `space = edge_count ÷ 100` may be a good heuristic. The sampler will take up
O(`sizeof(t)*ndims(initailizer)*space`) space with a constant factor around 4 bytes.

Draw a hypergraph from the sampler with `rand(sampler, edges)`.
"""
function Kronecker_sampler(initializer::AbstractArray, power::Integer; Type::Type{T}=Int, space=1,
    inner_power = max(1, min(power, Int(log(max(1,space)) ÷ log(max(2,length(initializer)))))),
    outer_power = cld(power, inner_power)) where T <: Integer

    Base.require_one_based_indexing(initializer)
    initializer_size = first(size(initializer))
    all(size(initializer) .== initializer_size) || ArgumentError("initializer must be square")

    inner_power = cld(power, outer_power)
    remainder = initializer_size ^ (power - inner_power*(outer_power-1))

    initializer = kronecker_power(initializer, inner_power)
    is = [Tuple(i).-one(T) for i in vec(eachindex(IndexCartesian(), initializer))]
    dist = AliasTable(vec(initializer), is)

    K = ndims(initializer)

    Kronecker_sampler{T, K}(outer_power, dist, initializer_size, remainder)
end

"""
    rand(s::Kronecker_sampler)
    rand(s::Kronecker_sampler, e::Integer)

Draw an edge or a hypergraph with `e` edges from the sampler `s`.
"""
function Base.rand(rng::AbstractRNG, s::Kronecker_sampler{T, K}) where {T, K}
    mult = s.remainder
    edge = MVector{K}(rand(rng, s.dist)) .% mult
    for _ in 1:(s.power-1)
        edge .+= rand(rng, s.dist) .* mult
        mult *= s.initializer_size
    end
    NTuple{K}(edge)
end

max_stars_repo_name: LilithHafner/FBD

max_issues_repo_name: LilithHafner/FBD

max_forks_repo_name: LilithHafner/FBD


text:
julia> sum(diff(Q(1:10^5)) .< 0)
49798

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


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

JLLWrappers.@generate_main_file_header("Chemfiles")
JLLWrappers.@generate_main_file("Chemfiles", UUID("78a364fa-1a3c-552a-b4bb-8fa0f9c1fcca"))
end  # module Chemfiles_jll