text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

include("../../src/Sparso.jl")
include("../../src/simple-show.jl")
include("./utils.jl")
using Sparso

set_options(SA_ENABLE, SA_VERBOSE, SA_USE_SPMP, SA_CONTEXT)

# The original pcg_symgs
#include("./pcg-symgs.jl")
function pcg_symgs(x, A, b, tol, maxiter)
    total_time = time()

    trsv_time = 0.
    spmv_time = 0.
    dot_time = 0.
    blas1_time = 0.

    L = tril(A)
    U  = SparseMatrixCSC{Cdouble, Cint}(spdiagm(1./diag(A))*triu(A))

    spmv_time -= time()
    #r = b - A * x
    r = b - Sparso.SpMV(A,x)
    spmv_time += time()

    blas1_time -= time()
    #normr0 = norm(r)
    normr0 = Sparso.norm(r)
    z = copy(r)
    blas1_time += time()

    rel_err = 1

    trsv_time -= time()
    Base.SparseMatrix.fwdTriSolve!(L, z)
    Base.SparseMatrix.bwdTriSolve!(U, z)
    trsv_time += time()

    blas1_time -= time()
    p = copy(z) #NOTE: do not write "p=z"! That would make p and z aliased (the same variable)
    #rz = dot(r, z)
    rz = Sparso.dot(r,z)
    blas1_time += time()

    k = 1
    while k <= maxiter

        set_matrix_property(Dict(
            :A => SA_CONST_VALUED|SA_MAXIMAL_STRUCTURED|SA_SYMM_VALUED|SA_SYMM_STRUCTURED,
            :L => SA_CONST_VALUED|SA_MAXIMAL_STRUCTURED,
            :U => SA_CONST_VALUED|SA_MAXIMAL_STRUCTURED
            )
        )

        old_rz = rz

        spmv_time -= time()
        #Ap = A*p
        Ap = Sparso.SpMV(A,p)
        spmv_time += time()

        blas1_time -= time()
        #alpha = old_rz / dot(p, Ap)
        alpha = old_rz / Sparso.dot(p, Ap)
        #x += alpha * p
        Sparso.WAXPBY!(x,1,x,alpha,p)
        #r -= alpha * Ap
        Sparso.WAXPBY!(r,1,r,-alpha,Ap)
        #rel_err = norm(r)/normr0
        rel_err = Sparso.norm(r)/normr0
        blas1_time += time()

        if rel_err < tol
            break
        end

        blas1_time -= time()
        #z = copy(r)
        Sparso.copy!(z, r)
        blas1_time += time()

        trsv_time -= time()
        Base.SparseMatrix.fwdTriSolve!(L, z)
        Base.SparseMatrix.bwdTriSolve!(U, z)
        trsv_time += time()

        blas1_time -= time()
        #rz = dot(r, z)
        rz = Sparso.dot(r,z)
        beta = rz/old_rz
        #p = z + beta * p
        Sparso.WAXPBY!(p,1,z,beta,p)
        blas1_time += time()

        k += 1
    end
    total_time = time() - total_time
    println("total = $(total_time)s trsv_time = $(trsv_time)s ($((12.*(nnz(L) + nnz(U)) + 2.*8*(size(L, 1) + size(L, 2)))*k/trsv_time/1e9) gbps) spmv_time = $(spmv_time)s ($((12.*nnz(A) + 8.*(size(A, 1) + size(A, 2)))*(k + 1)/spmv_time/1e9) gbps) blas1_time = $blas1_time")

    return x, k, rel_err
end

if length(ARGS) == 0
    m = 1000
    A = generate_symmetric_nonzerodiagonal_sparse_matrix(m)
#    A = matrix_market_read("tiny-diag.mtx", true, true)
else
    A = matrix_market_read(ARGS[1], true, true)
end
m       = size(A, 1)
x       = zeros(Float64, m)
b       = ones(Float64, m)
tol     = 1e-7
maxiter = 20000

println("Original: ")
x, k, rel_err = @acc pcg_symgs(x, A, b, tol, maxiter)
#println("\tsum of x=", sum(x))
println("\tk=", k)
println("\trel_err=", rel_err)

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
module PlotlyJS

using Reexport
@reexport using PlotlyBase
using JSON
using Base.Iterators
using REPL, Pkg

# need to import some functions because methods are meta-generated
import PlotlyBase:
    restyle!, relayout!, update!, addtraces!, deletetraces!, movetraces!,
    redraw!, extendtraces!, prependtraces!, purge!, to_image, download_image,
    restyle, relayout, update, addtraces, deletetraces, movetraces, redraw,
    extendtraces, prependtraces, prep_kwargs, sizes, savefig, _tovec

using WebIO
using JSExpr
using JSExpr: @var, @new
using Blink
using Requires

export plot

# globals for this package
const _pkg_root = dirname(dirname(@__FILE__))
const _js_path = joinpath(_pkg_root, "assets", "plotly-latest.min.js")
const _js_cdn_path = "https://cdn.plot.ly/plotly-latest.min.js"
const _mathjax_cdn_path =
    "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS-MML_SVG"

const _autoresize = [true]
autoresize(b::Bool) = (_autoresize[1] = b; b)
autoresize() = _autoresize[1]

abstract type AbstractPlotlyDisplay end

# include the rest of the core parts of the package
include("display.jl")
include("util.jl")
# include("savefig.jl")

function docs()
    schema_path = joinpath(dirname(dirname(@__FILE__)), "deps", "schema.html")
    if !isfile(schema_path)
        msg = "schema docs not built. Run `Pkg.build(\"PlotlyJS\")` to generate"
        error(msg)
    end
    w = Blink.Window()
    wait(w.content)
    for f in [
        "https://ajax.googleapis.com/ajax/libs/jquery/1.12.4/jquery.min.js",
        "http://maxcdn.bootstrapcdn.com/bootstrap/3.3.7/css/bootstrap.min.css",
        "http://maxcdn.bootstrapcdn.com/bootstrap/3.3.7/js/bootstrap.min.js"
        ]
        Blink.load!(w.content, f)
        wait(w.content)
    end
    Blink.content!(w, "html", open(f->read(f, String), schema_path), fade=false)
end

function PlotlyBase.savefig(p::SyncPlot, args...)
    has_orca = haskey(Pkg.installed(), "ORCA")
    if has_orca
        error("Please call `using ORCA` to save figures")
    end

    if Base.isinteractive()
        msg = "Saving figures requires the ORCA package."
        msg *= " Would you like to install it? (Y/n) "
        print(msg)
        answer = readline()
        if length(answer) == 0
            answer = "y"
        end
        if lowercase(answer)[1] == 'y'
            println("here!!")
            println("Ok. Installing ORCA now...")
            Pkg.add("ORCA")
            info("Please call `using ORCA` and try saving your plot again")
            return
        end
    end
    msg = "Please install ORCA separately, then call `using ORCA` and try again"
    error(msg)
end

function __init__()
    @require ORCA="200b8544-ab2f-11e8-2d2a-470a6868b879" include("savefig_orca.jl")
    if !isfile(_js_path)
        info("plotly.js javascript libary not found -- downloading now")
        include(joinpath(_pkg_root, "deps", "build.jl"))
    end
end

end # module

ext: jl


text:
module FSimPlots

using FlightSims
using Plots
import Plots: plot, plot!
using Transducers
using UnPack
using LinearAlgebra

export plot, plot!

include("multicopters.jl")

end


text:
using LinearAlgebra
using CompScienceMeshes
using SauterSchwab3D

const pI   = point(1,5,3)
const pII  = point(2,5,3)
const pIII = point(7,1,0)
const pIV  = point(3,4,0)

const qI  = point(10,11,12)
const qII  = point(10,11,13)
const qIII = point(11,11,12)

const P = simplex(pI,pII,pIII,pIV)
const Q = simplex(qI,qII,qIII)

const sing = SauterSchwab3D.singularity_detection(P,Q)

Accuracy2 = 16
pd_ref = PositiveDistance5D(sing,SauterSchwab3D._legendre(Accuracy2,0.0,1.0))

function integrand(x,y)
			return(exp(-im*1*norm(x-y))/(4pi*norm(x-y)))
end

function INTEGRAND(u,v)
   j1 = volume(P) * factorial(dimension(P))
   j2 = volume(Q) * factorial(dimension(Q))
   x = barytocart(P,u)
   y = barytocart(Q,v)
   output = integrand(x,y)*j1*j2
   return(output)
end

print("Ref: ")
ref = sauterschwab_parameterized(INTEGRAND, pd_ref)
println(ref)
println()

res_tp =[]
res_sp =[]
res_gm =[]
n1 = []
n2 = []
n3 = []

for i in 2:1:14
   Accuracy = i
   pd = PositiveDistance5D(sing,SauterSchwab3D._legendre(Accuracy,0.0,1.0))

   int_tp = sauterschwab_parameterized(INTEGRAND, pd)

   num_pts = length(pd.qps)^5
   push!(n1,num_pts)
   push!(res_tp,int_tp)
end

for i in 2:1:7
   Accuracy = i
   pd_s = PositiveDistance5D_S(sing,(SauterSchwab3D._shunnham3D(Accuracy),
                                SauterSchwab3D._shunnham2D(Accuracy)))

   int_sp = sauterschwab_parameterized(INTEGRAND, pd_s)

   num_pts = length(pd_s.qps[1])*length(pd_s.qps[2])
   push!(n2,num_pts)
   push!(res_sp,int_sp)
end

for i in 2:1:12
   Accuracy = i
   pd_gm = PositiveDistance5D_S(sing,(SauterSchwab3D._grundmannMoeller3D(2*Accuracy-1),
                                SauterSchwab3D._grundmannMoeller2D(2*Accuracy-1)))

   int_gm = sauterschwab_parameterized(INTEGRAND, pd_gm)

   num_pts = length(pd_gm.qps[1])*length(pd_gm.qps[2])
   push!(n3,num_pts)
   push!(res_gm,int_gm)
end

err_tp = norm.(res_tp.-ref)/norm(ref)
err_sp = norm.(res_sp.-ref)/norm(ref)
err_gm = norm.(res_gm.-ref)/norm(ref)

println(err_tp)
println(err_sp)
println(err_gm)

using Plots
plot( yaxis=:log, xaxis=:log, fontfamily="Times")
plot!(n1,err_tp, label="Gauss Tensor-Product",markershape=:circle)
plot!(n2,err_sp, label="Simplex Tensor-Product",markershape=:rect)
#plot!(n3,err_gm, label="Simplex-Product GM",markershape=:x)
plot!(xlims=(5,10^5),ylims=(1e-10,1e-1))
plot!(xlabel="Quad. points/Func. evals.", ylabel="Rel. Error.", title="Positive Distance 5D",legend=:bottomleft)

using BenchmarkTools

ref = 0.0024685313290508802 + 0.002663638878229659im

pd = PositiveDistance5D(sing,SauterSchwab3D._legendre(4,0.0,1.0))

int_tp = sauterschwab_parameterized(INTEGRAND, pd)

num_pts = length(pd.qps)^5
println(num_pts)
err_tp = norm.(int_tp.-ref)/norm(ref)
println(err_tp)
@time for i in 1:100 sauterschwab_parameterized(INTEGRAND, pd) end

pd_s = PositiveDistance5D_S(sing,(SauterSchwab3D._shunnham3D(4),
SauterSchwab3D._shunnham2D(4)))

int_sp = sauterschwab_parameterized(INTEGRAND, pd_s)

num_pts = length(pd_s.qps[1])*length(pd_s.qps[2])
println(num_pts)
err_sp = norm.(int_sp.-ref)/norm(ref)
println(err_sp)
@time for i in 1:100 sauterschwab_parameterized(INTEGRAND, pd_s) end

#=

pd_gm = PositiveDistance5D_S((SauterSchwab3D._grundmannMoeller3D(2*4-1),
SauterSchwab3D._grundmannMoeller2D(2*4-1)))

int_gm = sauterschwab_parameterized(INTEGRAND, pd_gm)

num_pts = length(pd_gm.qps[1])*length(pd_gm.qps[2])
println(num_pts)
err_gm = norm.(int_gm.-ref)/norm(ref)
println(err_gm)
=#


text:
function EspTornado(EspPump; ns=15, rfqs=30:5:60, stgs=false, fqs=false)
    stgs==false ? Stages=EspPump.RefStg : Stages=stgs
    fqs==false ? Freq=EspPump.RefFreq : Freq=fqs

    # ESP Head Polynomial Construction
    HeadCoef=convert(Vector,EspPump[[:H0,:H1,:H2,:H3,:H4,:H5,:H6]]).*1   # Extract Polynomial Head Coefficient
    HeadEq=Poly(HeadCoef)                                                  # Make the Polynomial equation

    # ESP Power Polynomial Construction
    PowCoef=convert(Vector,EspPump[[:P0,:P1,:P2,:P3,:P4,:P5,:P6]]).*1  #Extract Polynomial Head Coefficient
    PowEq=Poly(PowCoef)                                                  # Make the Polynomial equation

    CapRangeRef=range(0,stop=EspPump.AOF,length=ns)
    HeadRef=HeadEq(CapRangeRef)

    CapRange=zeros(ns,size(rfqs,1))
    Head=zeros(ns,size(rfqs,1))
    CapMin=zeros(size(rfqs,1))
    CapMax=zeros(size(rfqs,1))
    HeadMin=zeros(size(rfqs,1))
    HeadMax=zeros(size(rfqs,1))
    for i=1:size(rfqs,1)
        CapRange[:,i]=range(50,stop=EspPump.AOF,length=ns).*(rfqs[i]/EspPump.RefFreq)  # Capacity Range in bbl/d
        Head[:,i]=HeadEq(CapRangeRef).*Stages.*(rfqs[i]/EspPump.RefFreq).^2             # Pump Head  stg at Reference Frequency
        CapMin[i]=EspPump.Min.*(rfqs[i]/EspPump.RefFreq)
        CapMax[i]=EspPump.Max.*(rfqs[i]/EspPump.RefFreq)
        HeadMin[i]=HeadEq(EspPump.Min).*Stages.*(rfqs[i]/EspPump.RefFreq).^2
        HeadMax[i]=HeadEq(EspPump.Max).*Stages.*(rfqs[i]/EspPump.RefFreq).^2
    end

    #custom frecuency
    CapCustom=range(50,stop=EspPump.AOF,length=ns).*(fqs/EspPump.RefFreq)  # Capacity Range in bbl/d
    HeadCustom=HeadEq(CapRangeRef).*Stages.*(fqs/EspPump.RefFreq).^2             # Pump Head  stg at custom Frequency

    return CapRange,Head,CapMin,CapMax,HeadMin,HeadMax, CapCustom, HeadCustom

end


text:
module SmartTransitionSim

#packages
using OpenStreetMapX
using StatsBase
using SparseArrays
using Base.Iterators
using LightGraphs
using DataFrames
using DataStructures
using Random

export get_map_data

#types
export Rect, Agent

#functions
#generate_agents.jl
export generate_agents, pick_random_node, get_nodes_set
#rerouting.jl
export k_shortest_path_rerouting!, yen_a_star
#simulations.jl
export simulation_run, gather_statistics, run_parameter_analysis
#traffic_model.jl
export get_max_densities, traffic_constants, init_traffic_variables, next_edge
export update_weights!, update_event_agent!, update_agents_position!, update_weights_and_events!
#utilities
export print_map_statistics

#files
include("types.jl")
include("generate_agents.jl")
include("rerouting.jl")
include("simulations.jl")
include("traffic_model.jl")
include("utilities.jl")

end


text:
#############################
# Resampling uncertain tuples
#############################
include("uncertain_tuples/test_uncertain_tuples.jl")


text:
module PrecisionAquaculture

using DataFrames
using CSV
using Dates

include("utils.jl")
using .Utils
export readdata,groupedmean,extractgrmean,hourly,extractsensor
export FishDType, FVal

greet() = print("Hello World!")

end # module

max_stars_repo_name: IBM/precision-aquaculture

max_issues_repo_name: IBM/precision-aquaculture

max_forks_repo_name: IBM/precision-aquaculture


text:
#############################
#   Programming utilities   #
#############################

for N = 1:4
    # Equivalent to `integral_type(f, x) = typeof(p₀[1] * f(x...)` but inferrable
    @eval function integral_type(f, x::NTuple{$N,T}) where {T}
        Base.Cartesian.@nexprs $N d->(x_d = x[d])
        return typeof(p₀[1] * Base.Cartesian.@ncall($N, f, x))
    end

    # Equivalent to `evaluate(f, wp) = wp.w * f(wp.p...)`, but faster
    @eval function evaluate(f::F, v::SVector{$N,T}) where {F, T}
        Base.Cartesian.@nexprs $N d->(x_d = v[d])
        return Base.Cartesian.@ncall($N, f, x)
    end
end

idem(x, t) = (x, t)

### Getters and setters
weights(wp::WPoints) = wp.w
points(wp::WPoints) = wp.p

###############################
#   Array related utilities   #
###############################

Base.length(r::Regions) = length(r.v)
Base.getindex(r::Regions, i::Int) = r.v[i]

function Base.push!(wp::WPoints{N,T,R}, t::Tuple{R,SVector{N,T}}) where {N,T,R}
    push!(wp.w, t[1])
    push!(wp.p, t[2])
end

######################
#   Math utilities   #
######################

chop(x::T) where {T} = ifelse(x < eps(T), zero(T), x)

Base.isless(r₁::Region{N}, r₂::Region{N}) where {N} = isless(r₁.E.x, r₂.E.x)
Base.isequal(r₁::Region{N}, r₂::Region{N}) where {N} = isequal(r₁.E.x, r₂.E.x)

####################
#   IO utilities   #
####################

Base.show(io::IO, r::Region) = print(io, '(', r.x, ", ", r.h, ')')

Base.show(io::IO, m::MIME"text/plain", r::Regions) = show(io, m, r.v)
Base.show(io::IO, r::Regions) =
    (n = length(r.v); print(io, n, " subregion", n == 1 ? "" : "s"))


text:
# color vision deficiency (CVD) -  (color blind) - friendly
# https://zenodo.org/record/3381072
# Creative Commons Attribution 4.0 International
# 8, 12, and 15 color palettes for deuteranopia color blindness
# Martin Krzywinski martink@bcgsc.ca
# Methods and details: http://mkweb.bcgsc.ca/colorblind

#From Paul Tol: https://personal.sron.nl/~pault/
loadcolorscheme(:tol_bright, [
    colorant"#EE6677",
    colorant"#228833",
    colorant"#4477AA",
    colorant"#CCBB44",
    colorant"#66CCEE",
    colorant"#AA3377",
    colorant"#BBBBBB"
], "cvd", "cvd friendly, color blind friendly")

loadcolorscheme(:tol_muted, [
    colorant"#88CCEE",
    colorant"#44AA99",
    colorant"#117733",
    colorant"#332288",
    colorant"#DDCC77",
    colorant"#999933",
    colorant"#CC6677",
    colorant"#882255",
    colorant"#AA4499",
    colorant"#DDDDDD"
], "cvd", "cvd friendly, color blind friendly")

loadcolorscheme(:tol_light, [
    colorant"#BBCC33",
    colorant"#AAAA00",
    colorant"#77AADD",
    colorant"#EE8866",
    colorant"#EEDD88",
    colorant"#FFAABB",
    colorant"#99DDFF",
    colorant"#44BB99",
    colorant"#DDDDDD"
], "cvd", "cvd friendly, color blind friendly")

#From Color Universal Design (CUD): https://jfly.uni-koeln.de/color/
loadcolorscheme(:okabe_ito, [
    colorant"#E69F00",
    colorant"#56B4E9",
    colorant"#009E73",
    colorant"#F0E442",
    colorant"#0072B2",
    colorant"#D55E00",
    colorant"#CC79A7",
    colorant"#000000"
], "cvd", "cvd friendly, color blind friendly")

# 8, 12, and 15 color palettes for deuteranopia color blindness
# Martin Krzywinski martink@bcgsc.ca
# Methods and details: http://mkweb.bcgsc.ca/colorblind

loadcolorscheme(:mk_8, [
    colorant"#000000",
    colorant"#2271B2",
    colorant"#3DB7E9",
    colorant"#F748A5",
    colorant"#359B73",
    colorant"#d55e00",
    colorant"#e69f00",
    colorant"#f0e442",
], "cvd", "cvd friendly, color blind friendly")

loadcolorscheme(:mk_12, [
    colorant"#9F0162",
    colorant"#009F81",
    colorant"#FF5AAF",
    colorant"#00FCCF",
    colorant"#8400CD",
    colorant"#008DF9",
    colorant"#00C2F9",
    colorant"#FFB2FD",
    colorant"#A40122",
    colorant"#E20134",
    colorant"#FF6E3A",
    colorant"#FFC33B",
], "cvd", "cvd friendly, color blind friendly")

loadcolorscheme(:mk_15, [
    colorant"#68023F",
    colorant"#008169",
    colorant"#EF0096",
    colorant"#00DCB5",
    colorant"#FFCFE2",
    colorant"#003C86",
    colorant"#9400E6",
    colorant"#009FFA",
    colorant"#FF71FD",
    colorant"#7CFFFA",
    colorant"#6A0213",
    colorant"#008607",
    colorant"#F60239",
    colorant"#00E307",
    colorant"#FFDC3D",
], "cvd", "cvd friendly, color blind friendly")


text:
module JCModules

using StaticArrays
using CxxWrap
using JuliaChem.JERI

Base.include(@__MODULE__,"Globals.jl")

Base.include(@__MODULE__,"BasisStructs.jl")

Base.include(@__MODULE__,"MolStructs.jl")

end


text:
"""
```
DSGEVAR{T} <: AbstractDSGEVARModel{T}
```
implements a simple interface for combining a given DSGE model
with a VAR to create a DSGE-VAR. Confer with Del Negro and Schorfheide (2004),
Del Negro and Schorfheide (2006), Del Negro, Schorfheide, Smets, and Wouters (2007),
and/or Del Negro and Schorfheide (2009) for details about DSGE-VARs.
We recommend the first two papers as initial introductions to DSGE-VARs.

The recommended constructor requires the user to provide
(1) an `AbstractDSGEModel` object, (2) which structural shocks from
the DSGE to use, and (3) the subspec (optional, defaults to "ss0").
If the subspec "ss0" is used, then the result is a `DSGEVAR`
whose VAR component is "empty" in that the
observables, lags, and λ weight are not specified. The reason why
this constructor requires the user to specify which
structural shocks of DSGE to use is that this information is
DSGE-specific rather than information about the VAR.

However, we can also construct a `DSGEVAR` without having to
specify the structural shocks when calling the constructor,
although we still need to give an instance of an `AbstractDSGEModel`.

### Example
The code below instantiates an empty `DSGEVAR`
with `AnSchorfheide` as the underlying DSGE and then
calls `update!` on the empty DSGE-VAR
to add information about the desired DSGE-VAR spec.

```jldoctest; output = false
dsgevar = DSGEVAR(AnSchorfheide())
DSGE.update!(dsgevar, shocks = [:rm_sh, :z_sh, :g_sh],
    observables = [:obs_gdp, :obs_cpi, :obs_nominalrate],
    λ = 1., lags = 4)

# output

DSGE-VAR Model
observables:      3
data_vintage:     200310
DSGE model:       an_schorfheide
DSGE description: Julia implementation of model defined in 'Bayesian Estimation of DSGE Models' by Sungbae An and Frank Schorfheide: AnSchorfheide, ss0
```

### Fields

#### DSGE object
* `dsge::AbstractDSGEModel{T}`: underlying DSGE model object

#### DSGE-VAR Information
* `observables::OrderedDict{Symbol,Int}`: dictionary mapping observables
    of the VAR to their index in the matrices representing the DSGE-VAR
* `shocks::OrderedDict{Symbol,Int}`: dictionary mapping structural
    shocks in the DSGE to their index in the matrices representing the DSGE-VAR
* `lags::Int`: number of lags in the VAR
* `λ::T`: weight on the DSGE prior

#### Auxiliary Information
* `spec::String`: concatenates `dsgevar` with the spec of the DSGE, e.g.
    for `AnSchorfheide`, we have `dsgevar_an_schorfheide`.
* `subspec::String`: specifies the model subspecification.
    Cached here for filepath computation.
* `testing::Bool`: indicates whether the model is in testing mode.
    Currently, this setting has no uses in practice
"""
mutable struct DSGEVAR{T} <: AbstractDSGEVARModel{T}
    dsge::AbstractDSGEModel{T}
    observables::OrderedDict{Symbol,Int}
    shocks::OrderedDict{Symbol,Int}
    lags::Int
    λ::T
    spec::String
    subspec::String
    testing::Bool
end

function Base.show(io::IO, m::DSGEVAR)
    @printf io "DSGE-VAR Model\n"
    @printf io "observables:      %i\n" n_observables(m)
    @printf io "data_vintage:     %s\n" data_vintage(m.dsge)
    @printf io "DSGE model:       %s\n" spec(get_dsge(m))
    @printf io "DSGE description: %s\n" description(get_dsge(m))
end

function DSGEVAR(dsge::AbstractDSGEModel{T}, shocks::Vector{Symbol}, subspec::String = "ss0";
                 custom_settings::Dict{Symbol, Setting} = Dict{Symbol, Setting}(),
                 copy_dsge::Bool = false, testing = false) where {T<:Real}

    # Initialize specs
    spec     = "dsgevar_" * ModelConstructors.spec(dsge)
    subspec  = subspec

    # Initialize empty DSGEVAR
    if copy_dsge
        dsge = deepcopy(dsge)
    end
    m = DSGEVAR{T}(dsge, OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(), 0,
                   0., spec, subspec, testing)

    # Initialize shocks from DSGE
    update!(m; shocks = shocks, check_valid = false)

    # Initialize subspec
    init_subspec!(m)

    # Do checks of observables, shocks, and lags
    check_valid_dsgevar(m)

    return m
end

# Empty constructor
function DSGEVAR(dsge::AbstractDSGEModel{T}, subspec::String = "ss0";
                 custom_settings::Dict{Symbol, Setting} = Dict{Symbol, Setting}(),
                 copy_dsge::Bool = false, testing = false) where {T<:Real}

    # Initialize specs
    spec     = "dsgevar_" * ModelConstructors.spec(dsge)
    subspec  = subspec

    # Initialize empty DSGEVAR
    if copy_dsge
        dsge = deepcopy(dsge)
    end
    m = DSGEVAR{T}(dsge, OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(), 0,
                   0., spec, subspec, testing)

    # Initialize shocks from DSGE
    update!(m; check_valid = false)

    # Initialize subspec
    init_subspec!(m)

    # Do checks of observables, shocks, and lags
    check_valid_dsgevar(m)

    return m
end


text:
module TweedieDistributions

using Distributions
using Random

import Distributions: mean, var, succprob, failprob, insupport
import Base: rand, minimum, maximum

export Tweedie,
    CompoundPoissonGamma,
    rand,
    mean,
    var,
    succprob,
    failprob

abstract type AbstractTweedie{T} <: Distribution{Univariate,Continuous} end

"""
    Tweedie{T}
       μ::T
       p::T
       ϕ::T

Tweedie distribution. Typically used with 1 < ξ < 2 to give compound
gamma-Poisson. Has positive mass at zero, and is continuous for positive
Reals.
"""
struct Tweedie{T} <: AbstractTweedie{T}
    μ::T
    p::T
    ϕ::T
end

"""
    CompoundPoissonGamma{T} <: AbstractTweedie
        μ::T
        p::T
        ϕ::T

A Tweedie distribution with ``1 < p < 2`` is a compound Poisson-Gamma
distribution, making it more straightforward to draw random samples from.
"""
struct CompoundPoissonGamma{T} <: AbstractTweedie{T}
    μ::T
    p::T
    ϕ::T

    function CompoundPoissonGamma(μ::T, p::T, ϕ::T) where T
        μ ≥ 0 || throw("Mean must be nonnegative")
        ϕ > 0 || throw("Variance parameter must be positive")
        1 < p < 2 ||
            throw("CPG requires 1 < p < 2, otherwise use Tweedie.")
        new{T}(μ, p, ϕ)
    end
end

function CompoundPoissonGamma(Tw::Tweedie)
    CompoundPoissonGamma(Tw.μ, Tw.p, Tw.ϕ)
end

function CompoundPoissonGamma(μ, p, ϕ)
    μ_prom, p_prom, ϕ_prom = promote(μ, p, ϕ)
    CompoundPoissonGamma(μ_prom, p_prom, ϕ_prom)
end

@inline function _cpg_poissonmean(CPG::CompoundPoissonGamma)
    CPG.μ^(2 - CPG.p) / ((2 - CPG.p) * CPG.ϕ)
end

@inline function _cpg_gammashape(CPG::CompoundPoissonGamma)
    gam_shape = (2 - CPG.p) / (CPG.p - 1)
end

@inline function _cpg_gammarate(CPG::CompoundPoissonGamma)
    CPG.ϕ * (CPG.p - 1) * CPG.μ ^ (CPG.p - 1)
end

sampler(cpg::CompoundPoissonGamma) = CompoundPoissonGammaSampler(cpg)
rand(rng::AbstractRNG, cpg::CompoundPoissonGamma) = rand(rng, sampler(cpg))

# function rand(rng::Random.AbstractRNG, CPG::CompoundPoissonGamma{T}) where T
#     pois_mean = _cpg_poissonmean(CPG)
#     N = rand(rng, Poisson(pois_mean))
#     N == 0 && return zero(T)
#     gam_shape = N * _cpg_gammashape(CPG)
#     gam_rate = _cpg_gammarate(CPG)
# end

mean(Tw::AbstractTweedie) = Tw.μ
var(Tw::AbstractTweedie) = Tw.ϕ * Tw.μ ^ Tw.p

# rand(CPG::CompoundPoissonGamma) = rand(Random.GLOBAL_RNG, CPG)

function succprob(CPG::CompoundPoissonGamma{T}) where T
    N = _cpg_poissonmean(CPG)
    ccdf(Poisson(N), zero(T))
end

function failprob(CPG::CompoundPoissonGamma{T}) where T
    N = _cpg_poissonmean(CPG)
    cdf(Poisson(N), zero(T))
end

minimum(CPG::CompoundPoissonGamma{T}) where T = zero(T)
maximum(CPG::CompoundPoissonGamma) = Inf
# Need to restrict `x` to ::Real or dispatch is ambiguous
insupport(CPG::CompoundPoissonGamma, x::Real) = 0 ≤ x < Inf

struct CompoundPoissonGammaSampler{T} <: Sampleable{Univariate, Continuous}
    λ::T
    α::T
    θ::T
end

function CompoundPoissonGammaSampler(CPG::CompoundPoissonGamma{T}) where T
    λ = _cpg_poissonmean(CPG)
    α = _cpg_gammashape(CPG)
    θ = _cpg_gammarate(CPG)
    CompoundPoissonGammaSampler{T}(λ, α, θ)
end

function rand(rng::Random.AbstractRNG, cpgs::CompoundPoissonGammaSampler{T}) where T
    N = rand(rng, Poisson(cpgs.λ))
    N == 0 && return zero(T)
    rand(rng, Gamma(N * cpgs.α, cpgs.θ))
end
rand(cpgs::CompoundPoissonGammaSampler) = rand(Random.GLOBAL_RNG, cpgs)

end # module


text:
include("./logging.jl")

"""Annealing solver

Args:
    inst (Instance): instance

    temp_init (Float64): température courante
    temp_mini (Float64): température mini avant arrêt
    temp_coef (Float64): coefficiant de refroidissement
    temp (Float64): température courante

    nb_test (Int): Nombre total de voisins testés
    nb_move (Int): Nombre de voisins acceptés (améliorant ou non)
    nb_reject (Int): Nombre de voisins refusés
    nb_steps (Int): Nombre de paliers à température constante effectués
    step_size (Int): Nombre d'itérations à température constante

    nb_cons_reject (Int): Nombre de refus consécutifs
    nb_cons_reject_max (Int): Nombre maxi de refus consécutifs

    nb_cons_no_improv (Int): Nombre de tests non améliorants
    nb_cons_no_improv_max (Int): Nombre de maxi tests non améliorants

    min_acceptance_ratio (Float64): Taux d'acceptation minimal avant l'arrêt

    duration (Float64): Durée réelle (mesurée) de l'exécution
    durationmax (Float64): Durée max de l'exécution (--duration)
    starttime (Float64): Heure de début d'une résolution

    cursol (Solution): Solution courante
    bestsol (Solution): meilleure Solution rencontrée
    testsol (Solution): nouvelle solution courante potentielle
"""
mutable struct AnnealingSolver
    inst::Instance

    temp_init::Float64
    temp_mini::Float64
    temp_coef::Float64
    temp::Float64

    nb_test::Int
    nb_move::Int
    nb_reject::Int
    nb_steps::Int
    step_size::Int

    nb_cons_reject::Int
    nb_cons_reject_max::Float64

    nb_cons_no_improv::Int
    nb_cons_no_improv_max::Int

    min_acceptance_ratio::Float64

    duration::Float64
    durationmax::Float64
    starttime::Float64

    cursol::Solution
    bestsol::Solution
    testsol::Solution
end

"""Annealing solver outer constructor"""
function AnnealingSolver(inst::Instance;
                         temp_init=nothing, temp_init_rate=0.8, temp_mini=1e-6,
                         temp_coef=0.999_999_9, step_size=1,
                         nb_cons_reject_max=nothing,
                         nb_cons_no_improv_max=nothing,
                         min_acceptance_ratio=0.01,
                         startsol=nothing)

    if nb_cons_reject_max === nothing
        nb_cons_reject_max = inst.nb_planes^2
    end
    if nb_cons_no_improv_max === nothing
        nb_cons_no_improv_max = inst.nb_planes^2
    end

    cursol = startsol === nothing ? Solution(inst) : startsol

    # On calcule éventuellement la température initiale automatiquement
    if temp_init === nothing
        temp_init = guess_temp_init_cost(cursol, temp_init_rate)
    end

    bestsol = Solution(cursol)
    testsol = Solution(cursol)

    durationmax = 15*60
    duration = 0.0
    starttime = 0.0

    solver = AnnealingSolver(inst,
                             temp_init, temp_mini, temp_coef, temp_init,
                             0, 0, 0, 0, step_size,
                             0, nb_cons_reject_max,
                             0, nb_cons_no_improv_max,
                             min_acceptance_ratio,
                             durationmax, duration, starttime,
                             cursol, bestsol, testsol)
    return solver
end

function solve(sv::AnnealingSolver, neighbour_operator!;
               durationmax::Int = 0)
    ln2("BEGIN solve(AnnealingSolver)")
    lg2() && println(get_stats(sv))

    if durationmax != 0
        sv.durationmax = durationmax
    end

    sv.starttime = time_ns()/1_000_000_000
    while ! finished(sv)
        for iter_in_step in 1:sv.step_size
            copy!(sv.testsol, sv.cursol)
            neighbour_operator!(sv.testsol)
            sv.nb_test += 1
            if sv.testsol.cost <= sv.cursol.cost
                if sv.testsol.cost < sv.cursol.cost
                    copy!(sv.cursol, sv.testsol)
                    if sv.cursol.cost < sv.bestsol.cost
                        print_log(sv)
                        copy!(sv.bestsol, sv.cursol)
                    end
                    sv.nb_cons_reject = 0
                    sv.nb_cons_no_improv = 0
                    sv.nb_move += 1
                else
                    sv.nb_cons_reject += 1
                    sv.nb_cons_no_improv += 1
                end
            elseif rand() < exp(-(sv.testsol.cost - sv.cursol.cost) / sv.temp)
                ln4("Move ", sv.testsol.cost - sv.cursol.cost, " temp ", sv.temp,
                    " Proba ", exp(-(sv.testsol.cost - sv.cursol.cost) / sv.temp))
                copy!(sv.cursol, sv.testsol)
                sv.nb_cons_reject = 0
                sv.nb_move += 1
                sv.nb_cons_no_improv += 1
            else
                # ln5("No move ", sv.testsol.cost - sv.cursol.cost,
                    # "proba ", exp(-(sv.testsol.cost - sv.cursol.cost) / sv.temp))
                sv.nb_cons_reject += 1
                sv.nb_reject += 1
                sv.nb_cons_no_improv += 1
            end
        end
        sv.nb_steps += 1
        sv.temp = max(sv.temp_coef * sv.temp, sv.temp_mini)
    end
    ln2("END solve(AnnealingSolver)")
end

"""
Retourne true ssi l'état justifie l'arrêt de l'algorithme.
"""
function finished(sv::AnnealingSolver)
    # return sv.nb_cons_reject >= sv.nb_cons_reject_max
    sv.duration = time_ns()/1_000_000_000 - sv.starttime
    too_long = sv.duration >= sv.durationmax

    too_many_no_improv = sv.nb_cons_no_improv >= sv.nb_cons_no_improv_max
    too_many_cons_reject = sv.nb_cons_reject >= sv.nb_cons_reject_max

    # ratio = sv.nb_move > 10 ? sv.nb_move / sv.nb_test : 1.0
    # acceptance_ratio_too_low = ratio < sv.min_acceptance_ratio
    # stop = too_long || acceptance_ratio_too_low
    stop = too_long || too_many_no_improv || too_many_cons_reject
    if stop && lg1()
        println("\nSTOP car :")
        println("     sv.nb_cons_reject=$(sv.nb_cons_reject)")
        println("     sv.nb_cons_reject_max=$(sv.nb_cons_reject_max)")
        println("     sv.nb_cons_no_improv=$(sv.nb_cons_no_improv)")
        println("     sv.nb_cons_no_improv_max=$(sv.nb_cons_no_improv_max)")
        println("     sv.duration=$(sv.duration)")
        println("     sv.durationmax=$(sv.durationmax)")
        # println("     acceptance_ratio=$ratio")
        # println("     sv.min_acceptance_ratio=$(sv.min_acceptance_ratio)")
        println(get_stats(sv))
    end
    return stop
end

"""
Calcul d'une température initiale de manière à avoir un taux d'acceptation τ en démarrage

Args:
    taux_cible: pourcentage représentant le taux d'acceptation cible(e.g. 0.8)
    nb_degrad_max: nbre de degradation à accepter pour le calcul de la moyenne

On lance une suite de mutations (succession de mouvement systématiquement
acceptés). On relève le nombre et la moyenne des mouvements conduisant à une
dégradation du coût de la solution.

degrad : dégradation moyenne du coût pour deux mutations consécutives de coût croissant

La probabilité standard d'acceptation d'une mauvaise solution est :

..math::
    p = e^{ -degrad/T } = 0.8    =>    T = t_init = -degrad / ln(p)

avec :

- p = taux_cible = proba(t_init)
- degrad = moyenne des dégradations de l'énergie
- T = t_init = la température initiale à calculer

Example:
    On va lancer des mutations jusqu'à avoir 1000 dégradations.
    Si par exemple le coût des voisins forme une suite de la forme :
    990, 1010, 990, 1010, 990,...
    On devra faire 2000 mutations pour obtenir 1000 dégradations de valeur 20,
    d'où t_init = -degrad / ln(proba)
    proba = 0.8   =>  t_init = degrad * 4.5
    proba = 0.37  =>  t_init = degrad

Attention:
    Cette fonction n'est **pas** une méthode de AnnealingSolver.
    Elle a juste besoin d'une solution et du type de mouvement à effectuer.
    Ici, on suppose que le seul mouvement possible est swap!(sol::Solution)
    Mais il faudra pouvoir paramétrer cette méthode pour des voisinages différents.

"""
function guess_temp_init_degrad(sol::Solution, nb_degrad_max::Int;
                                taux_cible=0.8)
    # we apply the formula : exp(-delta_h/T_0) = taux_cible
    # from the formula we get -delta_h / ln(taux_cible) = T_0
    #t_init = 0    # stupide : pour faire une descente pure !
    nb_degrad = 0
    last_cost = sol.cost
    degrads = []
    cursol = Solution(sol)
    while nb_degrad < nb_degrad_max
        # We go back from initial solution
        swap!(cursol)
        degrad = last_cost - cursol.cost
        last_cost = cursol.cost
        if degrad > 0
            nb_degrad += 1
            push!(degrads, degrad)
        end
    end

    delta = mean(degrads)

    t_init = - delta / log(taux_cible)
    ln2("Temp init : ", t_init)
    return t_init
end

"""
Suppose que les variations de coût sont de l'ordre du coût moyen par avion.
"""
function guess_temp_init_cost(sol::Solution, taux_cible=0.8)
    delta = mean(sol.costs)/10
    t_init = - delta / log(taux_cible)
    ln1("Temp init : ", t_init)
    return t_init
end

function get_stats(sv::AnnealingSolver)
    # temp_init_rate=    $(sv.temp_init_rate)
    txt = "
    Paramètres de l'objet AnnealingSolver :
    step_size=         $(sv.step_size)
    temp_init=         $(sv.temp_init)
    temp_mini=         $(sv.temp_mini)
    temp_coef=         $(sv.temp_coef)
    nb_cons_reject_max=$(sv.nb_cons_reject_max)
    Etat de l'objet AnnealingSolver :
    nb_steps=$(sv.nb_steps) step_size=$(sv.step_size)
    nb_cons_reject=$(sv.nb_cons_reject) nb_cons_reject_max=$(sv.nb_cons_reject_max)
    nb_cons_no_improv=$(sv.nb_cons_no_improv) nb_cons_no_improv_max=$(sv.nb_cons_no_improv_max)
    nb_test=$(sv.nb_test)
    nb_move=$(sv.nb_move)
    nb_reject=$(sv.nb_reject)
    temp=$(sv.temp) temp_init=$(sv.temp_init)
    testsol.cost=$(sv.testsol.cost)
    cursol.cost=$(sv.cursol.cost)
    bestsol.cost=$(sv.bestsol.cost)
    sv.testsol.solver.nb_infeasable=$(sv.testsol.solver.nb_infeasable)
    "
    # temp_init_rate=    $(sv.temp_init_rate)"
    return replace(txt, r"^ {4}" => "")
end


text:
function test_simulate()
    d = 2
    lims = repeat([-1.0 1.0], d, 1)
    cdist = UniformContext(lims)
    rdist = GaussianRewards([x -> x[1], x -> -x[1]], 1.0)
    policy = FixedPolicy(1)
    traj = Trajectory(d)
    simulate!(10, traj, cdist, rdist, policy; verbose=false)
    simulate!(10, traj, cdist, rdist, policy; verbose=false)
    @test traj.t == 20
    @test size(traj.X) == (20, d)
    @test traj.a == fill(1, 20)

    # Test matrix extension
    # traj = Trajectory(d)
    # ContextualBandits.add_memory!(traj, 10)
    # X = rand(10, d)
    # a = ones(Int, 10)
    # r = rand(10)
    # total_regret = zeros(10)
    # ContextualBandits.extend!(traj, X, a, r, total_regret)
    # @test traj.X == X
    # @test traj.r == r
    # @test traj.a == a
    # @test traj.total_regret == total_regret
    # @test traj.t == 10
end

test_simulate()


text:
# Mel-Log Spectrum Approximation filter (MLSADF)

# MLSABaseFilter represents a base filter of the MLSADF.
mutable struct MLSABaseFilter <: Filter
    order::Int
    α::Float64
    delay::Vector{Float64}

    function MLSABaseFilter(order::Int, α::Float64)
        new(order, α, zeros(order+1))
    end
end

delay(f::Filter) = f.delay
allpass_alpha(f::MLSABaseFilter) = f.α

function filt!(f::MLSABaseFilter, x::Float64, coef::Vector{Float64})
    d = delay(f)
    α = allpass_alpha(f)

    d[1] = x
    d[2] = (1.0-α*α)*d[1] + α*d[2]

    result = 0.0
    for i=3:length(coef)
        @inbounds d[i] = d[i] + α*(d[i+1] - d[i-1])
        @inbounds result += d[i] * coef[i]
    end

    # special case
    if length(coef) == 2
        result += d[2] * coef[2]
    end

    # t <- t+1 in time
    for i=length(d):-1:3
        @inbounds d[i] = d[i-1]
    end

    result
end

# MLSACascadeFilter represents a cascade filter which contains MLSA base
# filters.
mutable struct MLSACascadeFilter <: Filter
    filters::Vector{MLSABaseFilter}
    padecoef::Vector{Float64}
    delay::Vector{Float64}

    function MLSACascadeFilter(order::Int, α::Float64, pade::Int)
        padecoef = Vector{Float64}(pade+1)

        if pade == 4
            padecoef = padecoef4th
        elseif pade == 5
            padecoef = padecoef5th
        else
            error("MLSADF: Order of pade approximation 4 or 5 is only supported.")
        end

        filters = Vector{MLSABaseFilter}(pade+1)
        for i=1:length(filters)
            filters[i] = MLSABaseFilter(order, α)
        end

        new(filters, padecoef, zeros(pade+1))
    end
end

allpass_alpha(f::MLSACascadeFilter) = allpass_alpha(f.filters[1])
padecoef(f::MLSACascadeFilter) = f.padecoef

function filt!(f::MLSACascadeFilter, x::Float64, coef::Vector{Float64})
    d = delay(f)
    pade = padecoef(f)
    result, feedback = 0.0, 0.0

    for i=length(pade):-1:2
        @inbounds d[i] = filt!(f.filters[i], d[i-1], coef)
        @inbounds val = d[i] * pade[i]
        if iseven(i)
            feedback += val
        else
            feedback -= val
        end
        result += val
    end

    d[1] = feedback + x
    result += d[1]
    result
end

mutable struct MLSADF <: MelGeneralizedCepstrumSynthesisFilter
    filters::Vector{MLSACascadeFilter}

    function MLSADF(order::Int, α::Float64; pade::Int=5)
        filters = Vector{MLSACascadeFilter}(2)
        filters[1] = MLSACascadeFilter(2, α, pade)
        filters[2] = MLSACascadeFilter(order+1, α, pade)
        new(filters)
    end
end

first(f::MLSADF) = f.filters[1]
last(f::MLSADF) = f.filters[2]
allpass_alpha(f::MLSADF) = allpass_alpha(first(f))
glog_gamma(f::MLSADF) = 0.0

function filt!(f::MLSADF, x::Float64, coef::Vector{Float64})
    filt!(last(f), filt!(first(f), x, [0.0, coef[2]]), coef)
end

function to_filtcoef{T<:Union{MelCepstrum,LinearCepstrum}}(f::MLSADF,
                                                           mc::SpectralParamState{T})
    @assert allpass_alpha(paramdef(mc)) == allpass_alpha(f)
    rawdata(mgc2b(mc))
end


text:
###########################################################
## Generic Entity functions
###########################################################
abstract type Entity end

radius(entity::Entity) = entity.radius
id(entity::Entity) = entity.id

"""
    calculate_distance_between(e1, e2)

Calculates the eucledean distance between entities `e1` and `e2`.
"""
calculate_distance_between(e1::Entity, e2::Entity) = sqrt((e1.x - e2.x)^2 + (e1.y - e2.y)^2)

"""
    calculate_angle_between(e1, e2)

Calculates the angle between entities `e1` and `e2`.
"""
calculate_angle_between(e1::Entity, e2::Entity) = rad2deg(atan2(e2.y - e1.y, e2.x - e1.x)) % 360

"""
    closest_point_to(e1, e2, min_distance = 3.0)

Find the closest point to the given entity `e1` near the given target entity `e2`, outside its given radius, with an added fudge of min_distance.
"""
function closest_point_to(e1::Entity, e2::Entity, min_distance::Float64 = 3.0)
    angle = calculate_angle_between(e1, e2)
    r = radius(e2) + min_distance
    x = e2.x + r * cos(deg2rad(angle))
    y = e2.y + r * sin(deg2rad(angle))

    return Position(x, y)
end

###########################################################
## Planet
###########################################################

struct Planet <: Entity
    id::String
    x::Float64
    y::Float64
    health::Int
    radius::Float64
    num_docking_spots::Int
    current_production::Int
    remaining_resources::Int
    owned::Bool
    owner_id::String
    docked_ships_ids::Vector{String}

    function Planet(tokens::Vector{String})
        id = shift!(tokens)
        x = parse(Float64, shift!(tokens))
        y = parse(Float64, shift!(tokens))
        hp = parse(Int, shift!(tokens))
        r = parse(Float64, shift!(tokens))
        docking = parse(Int, shift!(tokens))
        current = parse(Int, shift!(tokens))
        remaining = parse(Int, shift!(tokens))
        owned = shift!(tokens) != "0"
        owner_id = shift!(tokens)
        num_docked_ships = parse(Int, shift!(tokens))
        docked_ships_ids = Vector{String}(num_docked_ships)
        for i in 1:num_docked_ships
            docked_ships_ids[i] = shift!(tokens)
        end

        new(id, x, y, hp, r, docking, current, remaining, owned, owner_id, docked_ships_ids)
    end
end

"""
    isowned(planet)

Determines if the planet has an owner.
"""
isowned(planet::Planet) = planet.owned

"""
    isfull(planet)

Determines if the planet has been fully occupied (all possible ships are docked)
"""
isfull(planet::Planet) = length(planet.docked_ships_ids) >= planet.num_docking_spots

###########################################################
## Ship
###########################################################

struct Ship <: Entity
    owner_id::String
    id::String
    x::Float64
    y::Float64
    hp::Int
    radius::Float64
    vel_x::Float64
    vel_y::Float64
    docked::DockedStatus
    docked_planet_id::String
    progress::Int
    cooldown::Int

    function Ship(owner_id::String, tokens::Vector{String})
        id = shift!(tokens)
        x = parse(Float64, shift!(tokens))
        y = parse(Float64, shift!(tokens))
        hp = parse(Int, shift!(tokens))
        vel_x = parse(Float64, shift!(tokens))
        vel_y = parse(Float64, shift!(tokens))
        docked = DockedStatus(parse(Int, shift!(tokens)))
        docked_planet = shift!(tokens)
        progress = parse(Int, shift!(tokens))
        cooldown = parse(Int, shift!(tokens))

        new(owner_id, id, x, y, hp, Constants.SHIP_RADIUS, vel_x, vel_y, docked, docked_planet, progress, cooldown)
    end

end

"""
    isdocked(ship)

Test if current ship is undocked. Returns `true` if ship is `UNDOCKED`, `false` otherwise.
"""
isdocked(ship::Ship) = ship.docked != UNDOCKED

"""
    can_dock(ship, planet)

Determine whether a ship is close enough to planet so it can dock.
"""
function can_dock(ship::Ship, planet::Planet)
    calculate_distance_between(ship, planet) <= radius(planet) + Constants.DOCK_RADIUS
end

"""
    thrust(ship, speed, angle)

Generate a command to accelerate this ship.
"""
function thrust(ship::Ship, speed, angle)
    @sprintf("t %s %s %s", ship.id, floor(Int, speed), round(Int, angle))
end

"""
    dock(ship, planet)

Generate a command to dock to a planet.
"""
function dock(ship::Ship, planet::Planet)
    @sprintf("d %s %s", ship.id, planet.id)
end

"""
    undock(ship)

Generate a command to undock from the current planet.
"""
function undock(ship::Ship)
    @sprintf("u %s", ship.id)
end

function nearest_unoccupied_planet(ship::Ship, planets)
    d = Inf
    p = nothing
    for planet in planets
        if isowned(planet)
            if planet.owner_id == ship.owner_id && isfull(planet)
                continue
            end
        end
        dist = calculate_distance_between(ship, planet)
        if dist < d
            p = planet
            d = dist
        end
    end

    return p
end

###########################################################
## Position
###########################################################

struct Position <: Entity
    x::Float64
    y::Float64
end
radius(p::Position) = 0
id(p::Position) = "p"

max_stars_repo_name: xinnosuke/Halite-II

max_issues_repo_name: xinnosuke/Halite-II

max_forks_repo_name: xinnosuke/Halite-II


text:
using BallroomSkatingSystem
using Test

# fixture for tests
function BallroomSkatingSystem.get_rule11_table(::Nothing)
    return nothing
end

@testset "Individual Dances" begin
    results1 = DataFrame(Number = [11, 12, 13, 14, 15],
        JudgeA = [1, 2, 3, 4, 5],
        JudgeB = [1, 3, 2, 4, 5],
        JudgeC = [1, 3, 2, 5, 4],
        JudgeD = [5, 1, 2, 3, 4],
        JudgeE = [2, 1, 5, 3, 4])

    results2 = DataFrame(Number = [11, 12, 13],
        JudgeA = [1, 2, 3],
        JudgeB = [1, 2, 3],
        JudgeC = [2, 1, 3],
        JudgeD = [2, 1, 3],
        JudgeE = [3, 2, 1])

    results3 = DataFrame(Number = [11, 12, 13],
        JudgeA = [1, 2, 3],
        JudgeB = [1, 2, 3],
        JudgeC = [2, 1, 3],
        JudgeD = [1, 3, 2],
        JudgeE = [2, 1, 3])

    results4 = DataFrame(Number = [11, 12, 13, 14],
        JudgeA = [1, 2, 3, 4],
        JudgeB = [1, 2, 3, 4],
        JudgeC = [3, 2, 4, 1],
        JudgeD = [3, 2, 1, 4],
        JudgeE = [4, 2, 1, 3])

    results5 = DataFrame(Number = [11, 12, 13, 14, 15],
        JudgeA = [1, 2, 5, 3, 4],
        JudgeB = [1, 4, 2, 3, 5],
        JudgeC = [1, 5, 2, 3, 4],
        JudgeD = [2, 1, 5, 3, 4],
        JudgeE = [5, 2, 1, 3, 4])

    # Multiple majorities
    skating1_text, skating1 = skating_single_dance(results1)
    @test skating1[!, :Number] == [11, 12, 13, 14, 15]
    @test skating1[!, :Place] == [1, 2, 3, 4, 5]

    # Clear majorities
    skating2_text, skating2 = skating_single_dance(results2)
    @test skating2[!, :Number] == [11, 12, 13]
    @test skating2[!, :Place] == [2, 1, 3]

    # No majority
    skating3_text, skating3 = skating_single_dance(results3)
    @test skating3[!, :Number] == [11, 12, 13]
    @test skating3[!, :Place] == [1, 2, 3]

    # Shared places and wins with 2nd place
    skating4_text, skating4 = skating_single_dance(results4)
    @test skating4[!, :Number] == [11, 12, 13, 14]
    @test skating4[!, :Place] == [2.5, 1, 2.5, 4]

    # Only evaluate a specific set of competitors
    skating5_text, skating5 = skating_single_dance(results5)
    @test skating5[!, :Number] == [11, 12, 13, 14, 15]
    @test skating5[!, :Place] == [1, 2, 3, 4, 5]
end

@testset "Multiple Dances w/o Rule 11" begin
    dances = ["Waltz", "Tango", "Viennese Waltz", "Slowfox", "Quickstep"]

    places1 = DataFrame(Number = [20, 30, 40, 50, 60, 70],
        Waltz = [2, 1, 6, 4, 5, 3],
        Tango = [1, 2, 6, 4, 3, 5],
        Viennese = [5, 2, 1, 4, 3, 6],
        Slowfox = [3, 6, 1, 4, 2, 5],
        Quickstep = [2, 6, 3, 1, 4, 5])

    places2 = DataFrame(Number = [19, 29, 39, 49, 59, 69],
        Waltz = [1, 3, 2, 6, 5, 4],
        Tango = [1, 4, 2, 3, 5, 6],
        Viennese = [2, 1, 3, 4, 5, 6],
        Slowfox = [3, 2, 1, 5, 4, 6],
        Quickstep = [2, 1, 3, 5, 4, 6])

    places3 = DataFrame(Number = [18, 28, 38, 48, 58, 68],
        Waltz = [3, 2, 1, 4, 5, 6],
        Tango = [1, 2, 4, 5, 3, 6],
        Viennese = [1, 3, 2, 6, 5, 4],
        Slowfox = [1, 2, 5, 4, 3, 6],
        Quickstep = [1, 6, 3, 2, 5, 4])

    places4 = DataFrame(Number = [17, 27, 37, 47, 57, 67],
        Waltz = [1, 3, 2, 4, 6, 5],
        Tango = [1, 3, 5, 2, 6, 4],
        Slowfox = [2, 1, 4, 3, 5, 6],
        Quickstep = [5, 2, 4, 6, 1, 3])

    places5 = DataFrame(Number = [16, 26, 36, 46, 56, 66],
        Waltz = [1, 3, 2, 4, 5, 6],
        Tango = [2, 1, 3, 4, 5, 6],
        Viennese = [1, 2, 3, 4, 5, 6],
        Slowfox = [2, 1, 3, 4, 5, 6],
        Quickstep = [2, 4, 3, 1, 6, 5])

    # Rule 10c
    places_text, rule10_text, rule11_text, pl1 = skating_combined(dances, nothing, places1, [0, 0, 0, 0, 0])
    @test pl1[!, :Number] == places1[!, :Number]
    @test pl1[!, :Place] == [1, 2, 3, 4, 5, 6]

    # Rule 10b (second part) and 10d
    places_text, rule10_text, rule11_text, pl2 = skating_combined(dances, nothing, places2, [0, 0, 0, 0, 0])
    @test pl2[!, :Number] == places2[!, :Number]
    @test pl2[!, :Place] == [1, 2, 3, 4, 5, 6]

    # Rule 10b (first part) and 10d
    places_text, rule10_text, rule11_text, pl3 = skating_combined(dances, nothing, places3, [0, 0, 0, 0, 0])
    @test pl3[!, :Number] == places3[!, :Number]
    @test pl3[!, :Place] == [1, 2, 3, 4, 5, 6]

    # Rule 10a and 10b
    places_text, rule10_text, rule11_text, pl4 = skating_combined(dances, nothing, places4, [0, 0, 0, 0, 0])
    @test pl4[!, :Number] == places4[!, :Number]
    @test pl4[!, :Place] == [1, 2, 4, 3, 6, 5]

    # Rule 9
    places_text, rule10_text, rule11_text, pl5 = skating_combined(dances, nothing, places5, [0, 0, 0, 0, 0])
    @test pl5[!, :Number] == places5[!, :Number]
    @test pl5[!, :Place] == [1, 2, 3, 4, 5, 6]
end

@testset "Rule 11a, 11b, 11c" begin
    dances = ["Waltz", "Tango", "Viennese Waltz", "Slowfox", "Quickstep"]

    resultsW = DataFrame(Number = [25, 35, 45, 55, 65, 75],
        JudgeA = [1, 2, 4, 3, 6, 5],
        JudgeB = [1, 2, 5, 3, 6, 4],
        JudgeC = [2, 5, 3, 6, 1, 4],
        JudgeD = [5, 2, 1, 3, 4, 6],
        JudgeE = [3, 5, 1, 6, 4, 2])

    resultsT = DataFrame(Number = [25, 35, 45, 55, 65, 75],
        JudgeA = [1, 5, 3, 4, 2, 6],
        JudgeB = [2, 3, 1, 4, 5, 6],
        JudgeC = [2, 1, 3, 4, 5, 6],
        JudgeD = [3, 2, 1, 4, 5, 6],
        JudgeE = [4, 2, 3, 1, 5, 6])

    resultsV = DataFrame(Number = [25, 35, 45, 55, 65, 75],
        JudgeA = [1, 3, 5, 2, 4, 6],
        JudgeB = [1, 2, 3, 6, 4, 5],
        JudgeC = [2, 1, 3, 6, 5, 4],
        JudgeD = [2, 1, 6, 3, 5, 4],
        JudgeE = [3, 2, 1, 4, 6, 5])

    resultsS = DataFrame(Number = [25, 35, 45, 55, 65, 75],
        JudgeA = [2, 1, 5, 3, 6, 4],
        JudgeB = [2, 1, 6, 3, 4, 5],
        JudgeC = [5, 2, 3, 1, 6, 4],
        JudgeD = [1, 5, 2, 3, 6, 4],
        JudgeE = [5, 2, 3, 6, 4, 1])

    resultsQ = DataFrame(Number = [25, 35, 45, 55, 65, 75],
        JudgeA = [2, 1, 3, 5, 4, 6],
        JudgeB = [6, 5, 2, 1, 4, 3],
        JudgeC = [2, 1, 4, 5, 6, 3],
        JudgeD = [6, 5, 2, 1, 3, 4],
        JudgeE = [1, 2, 4, 3, 6, 5])

    results = [resultsW, resultsT, resultsV, resultsS, resultsQ]

    places_text, rule10_text, rule11_text, pl = skating_combined(dances, results)
    @test pl[!, :Number] == resultsW[!, :Number]
    @test pl[!, :Place] == [2, 1, 3, 4, 6, 5]
end

@testset "Rule 11b" begin
    dances = ["Waltz", "Tango", "Viennese Waltz", "Slowfox", "Quickstep"]

    resultsW = DataFrame(Number = [26, 36, 46, 56, 66, 76],
        JudgeA = [1, 2, 3, 4, 6, 5],
        JudgeB = [1, 2, 3, 5, 6, 4],
        JudgeC = [2, 5, 6, 3, 1, 4],
        JudgeD = [5, 2, 3, 1, 4, 6],
        JudgeE = [3, 5, 6, 1, 4, 2])

    resultsT = DataFrame(Number = [26, 36, 46, 56, 66, 76],
        JudgeA = [1, 5, 4, 3, 2, 6],
        JudgeB = [2, 3, 4, 1, 5, 6],
        JudgeC = [2, 1, 4, 3, 5, 6],
        JudgeD = [3, 2, 4, 1, 5, 6],
        JudgeE = [4, 2, 1, 3, 5, 6])

    resultsV = DataFrame(Number = [26, 36, 46, 56, 66, 76],
        JudgeA = [1, 2, 3, 4, 5, 6],
        JudgeB = [1, 3, 4, 5, 2, 6],
        JudgeC = [1, 4, 5, 6, 2, 3],
        JudgeD = [1, 5, 6, 2, 3, 4],
        JudgeE = [1, 6, 2, 3, 4, 5])

    resultsS = DataFrame(Number = [26, 36, 46, 56, 66, 76],
        JudgeA = [1, 5, 2, 3, 4, 6],
        JudgeB = [1, 6, 2, 3, 5, 4],
        JudgeC = [2, 3, 5, 1, 4, 6],
        JudgeD = [5, 2, 1, 3, 4, 6],
        JudgeE = [2, 3, 5, 6, 1, 4])

    resultsQ = DataFrame(Number = [26, 36, 46, 56, 66, 76],
        JudgeA = [1, 3, 2, 5, 6, 4],
        JudgeB = [5, 2, 6, 1, 3, 4],
        JudgeC = [1, 4, 2, 5, 3, 6],
        JudgeD = [5, 2, 6, 1, 4, 3],
        JudgeE = [2, 4, 1, 3, 5, 6])

    results = [resultsW, resultsT, resultsV, resultsS, resultsQ]

    places_text, rule10_text, rule11_text, pl = skating_combined(dances, results)
    @test pl[!, :Number] == resultsW[!, :Number]
    @test pl[!, :Place] == [1, 2, 3, 4, 5, 6]
end

@testset "Rule 11d" begin
    dances = ["Waltz", "Tango", "Viennese Waltz", "Slowfox", "Quickstep"]

    resultsW = DataFrame(Number = [25, 35, 45, 55, 65, 75],
        JudgeA = [1, 2, 4, 3, 6, 5],
        JudgeB = [1, 2, 5, 3, 6, 4],
        JudgeC = [2, 5, 3, 6, 1, 4],
        JudgeD = [5, 2, 1, 3, 4, 6],
        JudgeE = [3, 5, 1, 6, 4, 2])

    resultsT = DataFrame(Number = [25, 35, 45, 55, 65, 75],
        JudgeA = [1, 5, 3, 4, 2, 6],
        JudgeB = [2, 3, 1, 4, 5, 6],
        JudgeC = [2, 1, 3, 4, 5, 6],
        JudgeD = [3, 2, 1, 4, 5, 6],
        JudgeE = [4, 2, 3, 1, 5, 6])

    resultsV = DataFrame(Number = [25, 35, 45, 55, 65, 75],
        JudgeA = [1, 3, 5, 2, 4, 6],
        JudgeB = [1, 2, 3, 6, 4, 5],
        JudgeC = [2, 1, 3, 6, 5, 4],
        JudgeD = [2, 1, 6, 3, 5, 4],
        JudgeE = [3, 2, 1, 4, 6, 5])

    resultsS = DataFrame(Number = [25, 35, 45, 55, 65, 75],
        JudgeA = [2, 1, 5, 3, 6, 4],
        JudgeB = [2, 1, 6, 3, 4, 5],
        JudgeC = [5, 2, 3, 1, 6, 4],
        JudgeD = [1, 5, 2, 3, 6, 4],
        JudgeE = [5, 2, 3, 6, 4, 1])

    resultsQ = DataFrame(Number = [25, 35, 45, 55, 65, 75],
        JudgeA = [2, 1, 3, 5, 4, 6],
        JudgeB = [6, 5, 2, 1, 4, 3],
        JudgeC = [2, 1, 4, 5, 6, 3],
        JudgeD = [6, 5, 2, 1, 3, 4],
        JudgeE = [1, 2, 4, 3, 6, 5])

    results = [resultsW, resultsT, resultsV, resultsS, resultsQ]

    places_text, rule10_text, rule11_text, pl = skating_combined(dances, results)
    @test pl[!, :Number] == resultsW[!, :Number]
    @test pl[!, :Place] == [2, 1, 3, 4, 6, 5]
end

@testset "Rule 11b" begin
    dances = ["Waltz", "Tango", "Viennese Waltz", "Slowfox", "Quickstep"]

    resultsW = DataFrame(Number = [27, 37, 47, 57, 67, 77],
        JudgeA = [2, 3, 4, 1, 6, 5],
        JudgeB = [2, 3, 5, 1, 6, 4],
        JudgeC = [5, 6, 3, 2, 1, 4],
        JudgeD = [2, 3, 1, 5, 4, 6],
        JudgeE = [5, 6, 1, 3, 4, 2])

    resultsT = DataFrame(Number = [27, 37, 47, 57, 67, 77],
        JudgeA = [1, 3, 2, 4, 6, 5],
        JudgeB = [2, 1, 5, 4, 6, 3],
        JudgeC = [2, 3, 5, 4, 6, 1],
        JudgeD = [3, 1, 5, 4, 6, 2],
        JudgeE = [4, 3, 5, 1, 6, 2])

    resultsV = DataFrame(Number = [27, 37, 47, 57, 67, 77],
        JudgeA = [1, 2, 3, 4, 5, 6],
        JudgeB = [1, 3, 4, 5, 2, 6],
        JudgeC = [1, 4, 5, 6, 2, 3],
        JudgeD = [1, 5, 6, 2, 3, 4],
        JudgeE = [1, 6, 2, 3, 4, 5])

    resultsS = DataFrame(Number = [27, 37, 47, 57, 67, 77],
        JudgeA = [3, 1, 4, 5, 2, 6],
        JudgeB = [3, 1, 5, 6, 2, 4],
        JudgeC = [1, 2, 4, 3, 5, 6],
        JudgeD = [3, 5, 4, 2, 1, 6],
        JudgeE = [6, 2, 1, 3, 5, 4])

    resultsQ = DataFrame(Number = [27, 37, 47, 57, 67, 77],
        JudgeA = [2, 4, 1, 6, 5, 3],
        JudgeB = [6, 4, 5, 3, 1, 2],
        JudgeC = [2, 6, 1, 3, 5, 4],
        JudgeD = [6, 3, 5, 4, 1, 2],
        JudgeE = [1, 6, 2, 5, 4, 3])

    results = [resultsW, resultsT, resultsV, resultsS, resultsQ]

    places_text, rule10_text, rule11_text, pl = skating_combined(dances, results)
    @test pl[!, :Number] == resultsW[!, :Number]
    @test pl[!, :Place] == [1, 2, 3, 4, 5, 6]
end

@testset "Tie with Rule 11" begin
    # Artificial example to hit code
    dances = ["Waltz", "Tango", "Viennese Waltz"]

    resultsW = DataFrame(Number = [10, 20, 30],
        JudgeA = [1, 1, 1],
        JudgeB = [1, 1, 1],
        JudgeC = [1, 1, 1],
        JudgeD = [1, 1, 1],
        JudgeE = [1, 1, 1])

    resultsT = DataFrame(Number = [10, 20, 30],
        JudgeA = [1, 1, 1],
        JudgeB = [1, 1, 1],
        JudgeC = [1, 1, 1],
        JudgeD = [1, 1, 1],
        JudgeE = [1, 1, 1])

    resultsV = DataFrame(Number = [10, 20, 30],
        JudgeA = [1, 1, 1],
        JudgeB = [1, 1, 1],
        JudgeC = [1, 1, 1],
        JudgeD = [1, 1, 1],
        JudgeE = [1, 1, 1])

    results = [resultsW, resultsT, resultsV]

    places_text, rule10_text, rule11_text, pl = skating_combined(dances, results)
    @test pl[!, :Number] == resultsW[!, :Number]
    @test pl[!, :Place] == [2, 2, 2]
end


text:
abstract type Temperature <: Number end
struct Celsius <: Temperature
    value
end
struct Fahrenheit <: Temperature
    value
end
struct Kelvin <: Temperature
    value
end
Base.show(io::IO, t::Temperature)   = print(io, "$(t.value)")
Base.show(io::IO, t::Kelvin)        = print(io, "$(t.value) K")
Base.show(io::IO, t::Fahrenheit)    = print(io, "$(t.value) F")
Base.show(io::IO, t::Celsius)       = print(io, "$(t.value) C")

promote_rule(::Type{Kelvin}, ::Type{Celsius})     = Kelvin
promote_rule(::Type{Fahrenheit}, ::Type{Celsius}) = Celsius
promote_rule(::Type{Fahrenheit}, ::Type{Kelvin})  = Kelvin

Base.convert(::Type{Kelvin},  t::Celsius)     = Kelvin(t.value + 273.15)
Base.convert(::Type{Kelvin},  t::Fahrenheit)  = Kelvin(Celsius(t))
Base.convert(::Type{Celsius}, t::Kelvin)      = Celsius(t.value - 273.15)
Base.convert(::Type{Celsius}, t::Fahrenheit)  = Celsius((t.value - 32)*5/9)
Base.convert(::Type{Fahrenheit}, t::Celsius)  = Fahrenheit(t.value*9/5 + 32)
Base.convert(::Type{Fahrenheit}, t::Kelvin)   = Fahrenheit(Celius(t))

Base. +(a::Temperature, b::Temperature) = +(promote(a,b)...)
Base. -(a::Temperature, b::Temperature) = -(promote(a,b)...)

Base. +(a::T, b::T) where {T <: Temperature} = T(a.value + b.value)
Base. -(a::T, b::T) where {T <: Temperature} = T(a.value - b.value)

Base. *(a::Number, b::T) where {T <: Temperature} = T(a * b.value)
Base. *(a::T, b::Number) where {T <: Temperature} = T(a.value * b)
Base. /(a::Number, b::T) where {T <: Temperature} = T(a / b.value)


text:
"""
    ArcHybrid()

Arc-Hybrid system for transition dependency parsing.

Described in [Kuhlmann et al, 2011](https://www.aclweb.org/anthology/P11-1068.pdf),
[Goldberg & Nivre, 2013](https://aclweb.org/anthology/Q13-1033.pdf).
"""
struct ArcHybrid <: AbstractTransitionSystem end

initconfig(::ArcHybrid, graph::DependencyTree) = ArcHybridConfig(graph)
initconfig(::ArcHybrid, words) = ArcHybridConfig(words)

projective_only(::ArcHybrid) = true

transition_space(::ArcHybrid, labels=[]) =
    isempty(labels) ? [LeftArc(), RightArc(), Shift()] :
    [LeftArc.(labels)..., RightArc.(labels)..., Shift()]

struct ArcHybridConfig <: AbstractParserConfiguration
    stack::Vector{Int}
    buffer::Vector{Int}
    A::Vector{Token}
end

ArcHybridConfig(sentence) = stack_buffer_config(ArcHybridConfig, sentence)

buffer(cfg::ArcHybridConfig) = cfg.buffer
stack(cfg::ArcHybridConfig)  = cfg.stack
tokens(cfg::ArcHybridConfig) = cfg.A

# transition operations: leftarc, rightarc, shift

function apply_transition(f, cfg::ArcHybridConfig, a...; k...)
    σ, β, A = f(cfg.stack, cfg.buffer, cfg.A, a...; k...)
    return ArcHybridConfig(σ, β, A)
end

leftarc(cfg::ArcHybridConfig, args...; kwargs...) =
    apply_transition(leftarc_reduce, cfg, args...; kwargs...)

rightarc(cfg::ArcHybridConfig, args...; kwargs...) =
    apply_transition(rightarc_reduce, cfg, args...; kwargs...)

shift(cfg::ArcHybridConfig) = apply_transition(shift, cfg)

isfinal(cfg::ArcHybridConfig) = all(has_head, tokens(cfg))

"""
    static_oracle(cfg::ArcHybridConfig, tree, arc=untyped)

Static oracle for arc-hybrid dependency parsing.

Return a gold transition (one of LeftArc, RightArc, or Shift)
for parser configuration `cfg`.

See [Kuhlmann et al, 2011](https://www.aclweb.org/anthology/P11-1068.pdf)?
"""
function static_oracle(cfg::ArcHybridConfig, tree, arc=untyped)
    l = i -> arc(token(tree, i))
    if stacklength(cfg) > 0
        σ, s = popstack(cfg)
        if bufferlength(cfg) > 0
            b, β = shiftbuffer(cfg)
            if has_arc(tree, b, s)
                return LeftArc(l(s)...)
            end
        end
        if length(σ) > 0
            s2 = σ[end]
            if has_arc(tree, s2, s) && !any(k -> has_arc(tree, s, k), buffer(cfg))
                return RightArc(l(s)...)
            end
        end
    end
    if bufferlength(cfg) > 0
        return Shift()
    end
end

"""
    dynamic_oracle(cfg::ArgHybridConfig, tree, arc)

Dynamic oracle function for arc-hybrid parsing.

For details, see [Goldberg & Nivre, 2013](https://aclweb.org/anthology/Q13-1033.pdf).
"""
dynamic_oracle(cfg::ArcHybridConfig, tree, arc) =
    filter(t -> cost(t, cfg, tree) == 0, possible_transitions(cfg, tree, arc))

function is_possible(::LeftArc, cfg::ArcHybridConfig)
    s = last(stack(cfg))
    return s != 0 && !has_head(token(cfg, s)) && bufferlength(cfg) > 0
end

is_possible(::RightArc, cfg::ArcHybridConfig) =
    stacklength(cfg) > 1 && !has_head(token(cfg, last(stack(cfg))))

is_possible(::Shift, cfg::ArcHybridConfig) = bufferlength(cfg) >= 1

function cost(t::LeftArc, cfg::ArcHybridConfig, gold)
    # number of arcs (s0, d) and (h, s0) for h ϵ H and d ϵ D
    (σ, s0), (b, β) = popstack(cfg), shiftbuffer(cfg)
    H = length(σ) > 1 ? [σ[end] ; β] : β
    D = buffer(cfg)
    count(d -> has_arc(gold, s0, d), D) + count(h -> has_arc(gold, h, s0), H)
end

function cost(t::RightArc, cfg::ArcHybridConfig, gold)
    # number of arcs (s0,d) and (h,s0) for h, d ϵ B
    s0, n = last(stack(cfg)), 0
    for k in buffer(cfg)
        n += has_arc(gold, s0, k)
        n += has_arc(gold, k, s0)
    end
    return n
end

function cost(t::Shift, cfg::ArcHybridConfig, gold)
    # num of arcs (b, d), (h, b) s.t. h ϵ H, d ϵ D
    b = first(buffer(cfg))
    s, s0 = popstack(cfg)
    D = stack(cfg)
    H = length(D) > 1 ? s : Int[]
    return count(h -> has_arc(gold, h, b), H) + count(d -> has_arc(gold, b, d), D)
end

function possible_transitions(cfg::ArcHybridConfig, arc=untyped)
    s = last(stack(cfg))
    l = arc(token(cfg, s))
    transitions = [LeftArc(l...), RightArc(l...), Shift()]
    return filter(t -> is_possible(t, cfg), transitions)
end

function possible_transitions(cfg::ArcHybridConfig, tree::DependencyTree, arc=untyped)
    s = last(stack(cfg))
    l = arc(token(tree, s))
    transitions = [LeftArc(l...), RightArc(l...), Shift()]
    return filter(t -> is_possible(t, cfg), transitions)
end

==(cfg1::ArcHybridConfig, cfg2::ArcHybridConfig) =
    cfg1.stack == cfg2.stack && cfg1.buffer == cfg2.buffer && cfg1.A == cfg2.A


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 02329bbc-15a4-11ec-2522-df175bfb53e0
begin
    import Pkg
    # careful: this is _not_ a reproducible environment
    # activate the global environment
    Pkg.activate()

    using PlutoUI
	using FlightGNC
	using ComponentArrays, Transducers
	using Plots
end

# ╔═╡ 50becda1-a0e6-4e4c-98ca-a72992a5fc42
include("main_missiles.jl")

# ╔═╡ b9c304e7-9034-44d4-a4a6-a5d34d0baf58
include("sim_missiles.jl")

# ╔═╡ b4aa1cc3-07de-4717-a39e-8a9c4b1c8512
md"""
# Missile Guidance Simulation with Reactive Interactivity
2D/3D Biased Proportional Navigation Guidance Law is implemented.

In 2D, lateral acceleration command is given as follows:

$\displaystyle a_{cmd} = NV\dot{\lambda} - \frac{K_{r}\left(r\right)K_{\eta}\left(\eta\right)}{r}e_{\gamma_{f}}$
"""

# ╔═╡ 9fcf0c9b-d559-402f-883c-c53f61815c94
begin
	md"""
	 $\alpha$ = $(@bind α Slider(0.01:0.01:1; default=1, show_value=true))

	 $n$ = $(@bind n Slider(0 : 0.1 : 10; default=1, show_value=true))

	 $\chi_{M_{0}}$ [deg] = $(@bind χ_M_0_deg Slider(0:1:180; default=30, show_value=true))

	 $\chi_{f_{d}}$ [deg] = $(@bind χ_f_d_deg Slider(0:1:180; default=180, show_value=true))

	 $\sigma_{M_{\lim}}$ [deg] = $(@bind σ_M_lim_deg Slider(0:1:90; default=60, show_value=true))
	"""
end

# ╔═╡ b75636f4-b7f2-4356-9e2a-d90942c0071e
begin
	s_Bias  = ComponentArray(δ = 0.01, n = n, r_ref = 10E3, k = 9, m = 10)
	sim_plot(3, 2, α, deg2rad(χ_M_0_deg), deg2rad(χ_f_d_deg), deg2rad(σ_M_lim), s_Bias)
end

# ╔═╡ bce435b0-c6d4-4ac7-93aa-727ed08c73e5
# begin
# 	# Initial Conditions / Design parameters
#     dim     = 2
#     N       = 3
#     s_Bias  = ComponentArray(δ = 0.01, n = n, r_ref = 10E3, k = 9, m = 10)

#     V_M_0 	= 300
# 	γ_M_0 	= deg2rad(45)
# 	χ_M_0 	= deg2rad(χ_M_0_deg)
# 	χ_f_d 	= deg2rad(χ_f_d_deg)
# 	σ_M_lim = deg2rad(σ_M_lim_deg)
#     if dim == 2
#         (p_M_0, v_M_0)  = (zeros(dim), V_M_0*[sin(χ_M_0); cos(χ_M_0)])
#         (p_T_0, v_T_0)  = ([5E3; 1E3], zeros(dim))
#         γ_f_d = 0
#     elseif dim == 3
#         (p_M_0, v_M_0)  = (zeros(dim), V_M_0*[cos(γ_M_0)*sin(χ_M_0); cos(γ_M_0)*cos(χ_M_0); sin(γ_M_0)])
#         (p_T_0, v_T_0)  = ([10E3; 5E3; 5E3], 100*[-1; 0; 0])
# 		γ_f_d = -0.01
#     end
#     v̂_f_d   = [cos(γ_f_d)*sin(χ_f_d); cos(γ_f_d)*cos(χ_f_d); sin(γ_f_d)]

# 	s_BPNG = BPNG(N, dim, α, σ_M_lim, v̂_f_d, Bias_zero, s_Bias)
# 	# Bias = Bias_zero, Bias_IACG_StationaryTarget, Bias_IACG_StationaryTarget_2D

# 	# Execute simulation
# 	df = main(p_M_0, v_M_0, p_T_0, v_T_0, s_BPNG; A_M_max = 100)

# 	ts = df.time
#     p_Ms = df.sol |> Map(datum -> datum.pursuer.p) |> collect
#     p_Ts = df.sol |> Map(datum -> datum.evador.p)  |> collect

#     p_Ms = hcat(p_Ms...)'
#     p_Ts = hcat(p_Ts...)'

# 	# Plotting
# 	legend_string = ["Missile" "Target"]

# 	if dim == 2
# 		f = fig_print([p_Ms[:,1] p_Ts[:,1]], [p_Ms[:,2] p_Ts[:,2]], "", legend_string, "x [m]", "y [m]"; ar_val = :equal, save_file = 0, N_markers = 10)
# 		plot!(f, xlims=(0, 6E3), ylims=(-3E3, 3E3))
# 		f = fig_print([], [], "filename"; fig_handle = f) 	# To just save the result
# 	elseif dim == 3
# 		f = plot([p_Ms[:,1] p_Ts[:,1]], [p_Ms[:,2] p_Ts[:,2]], [p_Ms[:,3] p_Ts[:,3]], label = legend_string)
# 		f = fig_print([], [], "filename"; fig_handle = f, ar_val = :equal)
# 	end
# end

# ╔═╡ Cell order:
# ╟─b4aa1cc3-07de-4717-a39e-8a9c4b1c8512
# ╟─02329bbc-15a4-11ec-2522-df175bfb53e0
# ╠═50becda1-a0e6-4e4c-98ca-a72992a5fc42
# ╠═b9c304e7-9034-44d4-a4a6-a5d34d0baf58
# ╟─9fcf0c9b-d559-402f-883c-c53f61815c94
# ╟─b75636f4-b7f2-4356-9e2a-d90942c0071e
# ╟─bce435b0-c6d4-4ac7-93aa-727ed08c73e5


text:
function add_drag(w, range, point_id, slider_length, slideridx_s)
    m2id = mouse2id(w)
    # interaction
    @materialize mouse_buttons_pressed, mouseposition = w.inputs
    isoverpoint = const_lift(is_same_id, m2id, point_id)
    # single left mousekey pressed (while no other mouse key is pressed)
    key_pressed = const_lift(GLAbstraction.singlepressed, mouse_buttons_pressed, GLFW.MOUSE_BUTTON_LEFT)
    # dragg while key_pressed. Drag only starts if isoverpoint is true
    mousedragg = GLAbstraction.dragged(mouseposition, key_pressed, isoverpoint)
    preserve(foldp(value(slideridx_s), mousedragg) do v0, dragg
        if dragg != Vec2f0(0)
            idx_steps = round(Int, (dragg[1]/slider_length)*length(range))
            new_idx = clamp(v0 + idx_steps, 1, length(range))
            push!(slideridx_s, new_idx)
            return v0
        else # dragging started
            return value(slideridx_s)
        end
    end)
    slideridx_s
end
function add_play(slideridx_s, play_signal, range, rate=30.0)
    play_s = fpswhen(play_signal, rate)
    preserve(map(play_s, init=nothing) do t
        push!(slideridx_s, mod1(value(slideridx_s)+1, length(range)))
        nothing
    end)
end
function slider(
        range, window;
        startidx::Int=1,
        play_signal=Signal(false), slider_length=50mm
    )
    startpos = 2.1mm
    height = value(icon_size)
    point_id = Signal((0,0))
    slideridx_s = Signal(startidx)
    slider_s = map(slideridx_s) do idx
        range[clamp(idx, 1, length(range))]
    end
    add_drag(window, range, point_id, slider_length, slideridx_s)
    add_play(slideridx_s, play_signal, range)
    bb = map(icon_size) do is
        AABB(Vec3f0(0), Vec3f0(slider_length, is, 1))
    end
    line_pos = map(icon_size) do is
        Point2f0[(startpos, is/2), (slider_length, is/2)]
    end
    line = visualize(
        line_pos, :linesegment,
        boundingbox=bb, thickness=1.5mm
    ).children[]
    i = Signal(0)
    pos = Point2f0[(0, 0)]
    position = map(slideridx_s) do idx
        x = ((idx-1)/length(range-1))*slider_length
        pos[1] = (x, 0)
        pos
    end
    knob_scale = map(is->Vec2f0(is/3), icon_size)
    offset = map(line_pos, icon_size) do lp, is
        p = first(lp)
        Vec2f0(p - (is/6)) # - minus half knob scale
    end
    point_robj = visualize(
        (Circle, position),
        scale_primitive=true,
        offset=offset, scale=knob_scale,
        boundingbox=bb
    ).children[]
    push!(point_id, (point_robj.id, line.id))

    slider_s, Context(point_robj, line)
end

function maxdigits(range)
    if eltype(range) <: AbstractFloat
        return 7
    else
        return max(ndigits(first(range)), ndigits(last(range)))
    end
end
function play_widget(
        range, window=GLPlot.widget_screen!();
        startidx::Int=1
    )
    glyph_scale = GLVisualize.glyph_scale!('X')
    numberbox = map(icon_size) do is
        AABB(Vec3f0(0,-1mm, 0), Vec3f0(2is, is, 1))
    end
    target_s = value(icon_size) * 0.4
    scale = (target_s ./ glyph_scale)
    sliderlen = 70mm-(3*value(icon_size))-2mm
    play_button, play_stop_signal = GLVisualize.toggle_button(
        rot180(GLPlot.imload("play.png")), GLPlot.imload("break.png"), window
    )
    play_s = map(!, play_stop_signal)
    slider_s, slider_w = slider(range, window,
        startidx=startidx, play_signal=play_s,
        slider_length=sliderlen
    )
    number = visualize(
        map(GLVisualize.printforslider, slider_s),
        color=RGBA{Float32}(0.6, 0.6, 0.6,1),
        boundingbox=numberbox,
        relative_scale=scale
    )
    GLPlot.add_widget!(play_button, number, slider_w, window=window)
    slider_s
end
export play_widget
export add_widget!
function add_widget!(widgets...;
        delete=Signal(false),
        window=widget_screen!(delete=delete),
        height=value(icon_size)
    )
    scalings = map(widgets) do widget
        bb = value(boundingbox(widget))
        w = widths(bb)
        if w[2] > height # only scale when to big
            s = height/w[2]
            w.*s, minimum(bb)
        else
            w, minimum(bb)
        end
    end
    last_x = 0f0
    foreach(zip(scalings, widgets)) do scaling_widget
        scaling, widget = scaling_widget
        scale, offset = scaling
        xwidth = scale[1]-1mm
        place = SimpleRectangle{Float32}(last_x-offset[1]+0.5mm, -offset[2]+0.5mm, xwidth, scale[2]-1mm)
        layout!(place, widget)
        _view(widget, window, camera=:fixed_pixel)
        last_x += xwidth
    end
end

function item_area(la, deleted, item_height)
    y = la.y-item_height-2
    deleted && return SimpleRectangle(la.x, la.y, la.w, 0)
    return SimpleRectangle(la.x, y, la.w, item_height)
end

function widget_screen!(parentscreen=edit_screen; left_gap=1.5mm, delete=Signal(false))
    scroll = parentscreen.inputs[:menu_scroll]
    if isempty(parentscreen.children)
        last_area = map(parentscreen.area, icon_size, scroll) do a, ih, s
            return SimpleRectangle{Int}(left_gap, a.h-ih+s, a.w-left_gap, ih)
        end
    else
        last_area = last(parentscreen.children).area
    end
    itemarea = map(GLPlot.item_area, last_area, delete, icon_size)
    Screen(parentscreen, area=itemarea)
end

function choices(x::Vector)
    signal, vis = GLVisualize.choice_widget(x, edit_screen, area=(48mm, 8mm))
    w = widget_screen!()
    _view(vis, w, camera=:fixed_pixel)
    signal
end

function edit_toggle()
    edit_button, no_edit_signal = toggle_button(
        imload("play.png"), rotr90(imload("play.png")), edit_screen
    )
end

function visible_toggle()
    toggle_button(
        imload("showing.png"), imload("notshowing.png"), edit_screen
    )
end
function delete_toggle()
    delete_button, del_signal = button(
        imload("delete.png"), edit_screen
    )
end


text:
#!/usr/bin/env julia
#=
# Searching all vertices in graph/tree in a depth-first motion,
# marking each vertex as visited/not-visited while avoiding cycles,
# uses queue to remember next vertex to start search, backtracks on dead-end
#
# We've implemented in 1 ways here
# * xdfs() using recursion and 2 global queues
=#

function mark_neighbours(vertices, linkings)
  result = Dict(v => Set{eltype(vertices)}() for v in vertices)
  for (vertex,neighbour) in linkings
    push!(result[vertex], neighbour)
    push!(result[neighbour], vertex)
  end
  result
end

xroute, xmarked = [], []
function xdfs(adjlist, source, dest)
    global xroute
    global xmarked
    append!(xroute, [source])
    if source == dest
      return xroute
    end
    append!(xmarked, [source])

    for vertex in adjlist[source]
      if in(vertex, xmarked)
        continue
      end
      if vertex == dest
        append!(xroute, [vertex])
        return xroute
      end
    end

    for vertex in adjlist[source]
      if in(vertex, xmarked)
        continue
      end
      tmp_xroute = xdfs(adjlist, vertex, dest)
      if tmp_xroute != []
        return tmp_xroute
      end
    end

    []
end

const people = Set([:jean, :javert, :cosette, :gavroche, :éponine, :marius, :alice])
const friendships = [
    (:jean, :cosette),
    (:jean, :marius),
    (:cosette, :éponine),
    (:cosette, :marius),
    (:gavroche, :éponine),
    (:gavroche, :alice),
]

println("-------------------------------------------")

xroute, xmarked = [], []
@time println(xdfs(mark_neighbours(people, friendships), :jean, :jean))
xroute, xmarked = [], []
@time println(xdfs(mark_neighbours(people, friendships), :jean, :cosette))
xroute, xmarked = [], []
@time println(xdfs(mark_neighbours(people, friendships), :jean, :gavroche))
xroute, xmarked = [], []
@time println(xdfs(mark_neighbours(people, friendships), :jean, :javert))
xroute, xmarked = [], []
@time println(xdfs(mark_neighbours(people, friendships), :jean, :alice))

println("-------------------------------------------")


text:
"""
Various macro helpers.
"""
module Macros

using Genie, Genie.Configuration

export @devtools, @ifdevtools
export @run_with_time, @unless, @psst, @in_repl, @location_in_file

"""
    devtools()

Injects modules to be used in development mode, such as the `Gallium` debugger.
To be used as a concise one liner to include all the dev tools.

# Examples
```julia
julia> @devtools()
```
"""
macro devtools()
  if Configuration.is_dev()
    quote
      using Gallium
    end
  end
end

"""
    ifdevtools(expr::Expr)

Conditionally injects `expr::Expr` if the app runs in development mode.
The `Gallium` debugger is automatically also included.

# Examples
```julia
julia> @ifdevtools :(println("dev mode")) |> eval
dev mode
```
"""
macro ifdevtools(expr::Expr)
  if Configuration.is_dev()
    quote
      using Gallium
      $(esc(expr))
    end
  end
end

"""
    runwithtime(expr::Expr)

Prepends `@time` to `expr` if the app runs in development mode and evals - otherwise it simply evals `expr`.
To be used to automatically disable timed execution outside `dev` mode.

# Examples
```julia
julia> @run_with_time rand(2,2)
  0.059437 seconds (39.18 k allocations: 1.586 MB)
2×2 Array{Float64,2}:
 0.294566  0.653612
 0.264837  0.337146
```
"""
macro runwithtime(expr::Expr)
  if Configuration.is_dev()
    quote
      @time $(esc(expr))
    end
  else
    expr
  end
end

"""
    macro location_in_file()

Returns the location (file name and line number) in the file where it's called.
"""
macro location_in_file()
  :(String(esc(@__FILE__) * ":" * esc(@__LINE__)))
end

end


text:
using BinaryBuilder, Pkg

name = "MUMPS_seq"
version = v"4.10.0"

sources = [
  ArchiveSource("http://mumps.enseeiht.fr/MUMPS_4.10.0.tar.gz",
                "d0f86f91a74c51a17a2ff1be9c9cee2338976f13a6d00896ba5b43a5ca05d933"),
  DirectorySource("bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/MUMPS_4.10.0

# Patch from Coin-OR ThirdPartyMUMPS
(cd src && atomic_patch -p2 $WORKSPACE/srcdir/patches/mumps.patch)
# Patch from JuliaOpt CoinMumpsBuilder
(cd src && atomic_patch -p3 $WORKSPACE/srcdir/patches/quiet.diff)

cp Make.inc/Makefile.gfortran.SEQ Makefile.inc

make_args+=(OPTF=-O3
            CDEFS=-DAdd_
            LMETISDIR=${prefix}/lib
            IMETIS=-I${prefix}/include
            LMETIS='-L$(LMETISDIR) -lmetis'
            ORDERINGSF="-Dpord -Dmetis"
            CC="$CC -fPIC ${CFLAGS[@]}"
            FC="gfortran -fPIC ${FFLAGS[@]}"
            FL="gfortran -fPIC"
            LIBBLAS="-L${libdir} -lopenblas")

if [[ "${target}" == *-apple* ]]; then
  make_args+=(RANLIB=echo)
fi

# NB: parallel build fails
make alllib "${make_args[@]}"

mkdir -p ${prefix}/lib

cd libseq
mv libmpiseq.a ${prefix}/lib

cd ../lib
mv *.a ${prefix}/lib

cd ..
mkdir -p ${prefix}/include/mumps_seq
cp include/* ${prefix}/include/mumps_seq
cp libseq/*.h ${prefix}/include/mumps_seq
"""

platforms = expand_gfortran_versions(supported_platforms())

# The products that we will ensure are always built
products = [
    FileProduct("lib/libsmumps.a", :libsmumps_a),
    FileProduct("lib/libdmumps.a", :libdmumps_a),
    FileProduct("lib/libcmumps.a", :libcmumps_a),
    FileProduct("lib/libzmumps.a", :libzmumps_a),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency(PackageSpec(; name = "METIS_jll",
                                uuid = "d00139f3-1899-568f-a2f0-47f597d42d70",
                                version = v"4.0.3")),
    BuildDependency("OpenBLAS32_jll"),
    BuildDependency("CompilerSupportLibraries_jll"),
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=v"6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
#  function survivor_prob(evals,sp,nsvt)
#
#    n = size(evals)[1]
#    l = size(evals)[2]
#    pdt = zeros(Float64, n, l)
#
#    for i in 1:l
#      for j in 1:n
#        if evals[j,i] == 0
#          pdt[j,i] = 0
#        else
#          pdt[j,i] = evals[j,i] / sp[j,i]
#        end
#      end
#    end
#    final = sum(pdt,dims=2) / nsvt
#    return final
#
#  end
#

  function survivor_prob(evals,sp,nsvt)
    ev = sum(evals,dims=2) / nsvt
    ea = sum(sp, dims=2) / nsvt

    l     = length(ev)
    final = zeros(Float64, l)

    for i in 1:l
      if ev[i]==0
        final[i] = 0.
      else
        final[i] = ev[i] / ea[i]
      end
    end

   # final = ev ./ ea
    return final
  end


text:
module FinalBisc195

export normalizeDNA
export composition
export gc_content
export complement
export reverse_complement
export fasta_header
export 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(sequence)
    sequence = normalizeDNA(sequence) # make uppercase string, check invalid bases
    a = c = g = t = n = 0 # sets all 4 variables to `0`

    for base in sequence
        # add 1 to each base as it occurs
        #println(typeof(base))
        if base == 'A'
            a = a+1
        elseif base == 'G'
            g = g+1
        elseif base == 'T'
            t = t+1
        elseif base == 'C'
            c = c+1
        elseif base == 'N'
            n = n+1
        end
    end

    bases = Dict('A' => a,
                'G' => g,
                'T' => t,
                'N' => n,
                'C' => c)

    return bases
end

function gc_content(sequence)
    # Be sure to use `basecomposition()` in your answer.
    # Note: Since `basecomposition()` already calls `normalizeDNA`,
    # there's no need to call it here.
        ## throw an error if the string contains anything other than ACGT
    basecounts = composition(sequence)

    seqlength = length(sequence)
    gs = basecounts['G']
    cs = basecounts['C']

    return (gs + cs) / seqlength

end

function complement_char(base)
    complements =  Dict( 'A' => 'T',
                  'T' => 'A',
                  'G' => 'C',
                  'C' => 'G',
                  'N' => 'N')
    return complements[base]
end

function complement(base::String)
    base = normalizeDNA(base)
    outcome = []
    for letter in base
        letter = complement_char(letter)
        push!(outcome, letter)
    end
    outcome = join(outcome)
return outcome
end

function reverse_complement(sequence)
    sequence = normalizeDNA(sequence)
    backwards = reverse(sequence)
    backarray = []
    for i in 1:length(backwards)
       comp = complement_char(backwards[i])
       push!(backarray,comp)
       #@info backarray
    end
    return join(string.(backarray))
end

function fasta_header(header)
    if startswith(header, '>')
        header = chop(header, head = 1, tail = 0)
        header = string(header)
    return header
    else
        error("Invalid header (headers must start with '>')")
    end
end

function parse_fasta(path)
    ## Think through the components you need
    ## Does it make sense to define any containers at the beginning?
    ## How will you loop through the file?
    ## What do you need to get from each line?
    head_array = []
    body_array = []
    another_array = []

    for line in eachline(path)

        if startswith(line, '>')
            header_vector = fasta_header(line)
            #header_vector = Tuple([header_vector])
            push!(head_array, header_vector)
            another_array_string = join(another_array)
            if !isempty(another_array)
                push!(body_array, another_array_string)
            end
            another_array = []
        else
            body_line = line
            body_line = normalizeDNA(line)
            push!(another_array, body_line)
        end

    end
    another_array_string = join(another_array)
    push!(body_array, another_array_string)
    return head_array, body_array
end
end # module FinalBisc195


text:
module Utils
  files = ["ErrorAnalysis"]
  for f in files
    include("$f.jl")
  end

end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
"""
    normalized_laplacian(g)

Return the random-walk normalized laplacian of `g`, i.e. the matrix ``D^{-1}L``, where ``D`` is the diagonal degree matrix, and ``L`` is the standard laplacian matrix.
"""
function normalized_laplacian(g::SimpleGraph)
    L = float(laplacian_matrix(g))
    for i in vertices(g)
        if L[i,i] != 0
            L[i,:] /= L[i,i]
        end
    end
    return L
end

function normalized_laplacian(g::SimpleDiGraph)
    L = float(laplacian_matrix(g))
    for i in vertices(g)
        if L[i,i] != 0
            L[i,:] /= L[i,i]
        end
    end
    return L
end


text:
using .LatticePhysics.LatPhysBase: from, to, numSites, bonds, sites, latticeVectors, unitcell,
                   point, AbstractSite
import .LatticePhysics.LatPhysBase, .LatticePhysics

struct LatPhysLattice{LT <: LatPhysBase.AbstractLattice} <: AbstractLattice
    lattice::LT
    neighs::Matrix{Int}
    bonds::Vector{Vector{Int}}
end
export LatPhysLattice

function LatPhysLattice(lattice::LatPhysBase.AbstractLattice)
    # Build lookup table for neighbors
    # neighs[:, site] = list of neighoring site indices
    nested_bonds = [Int[] for _ in 1:numSites(lattice)]
    for b in LatPhysBase.bonds(lattice)
        push!(nested_bonds[LatPhysBase.from(b)], LatPhysBase.to(b))
    end
    max_bonds = maximum(length(x) for x in nested_bonds)

    neighs = fill(-1, max_bonds, numSites(lattice))
    for (src, targets) in enumerate(nested_bonds)
        for (idx, trg) in enumerate(targets)
            neighs[idx, src] = trg
        end
    end
    any(x -> x == -1, neighs) && @debug(
        "neighs is padded with -1 to indicated the lack of a bond. This is " *
        "due to the lattice having an irregular number of bonds per site."
    )
    LatPhysLattice(lattice, neighs, nested_bonds)
end

@inline Base.length(l::LatPhysLattice) = numSites(l.lattice)
@inline Base.ndims(l::LatPhysLattice) = ndims(l.lattice)
function Base.size(l::LatPhysLattice)
    lattice_extend = latticeVectors(l.lattice)
    uc_extend = latticeVectors(unitcell(l.lattice))
    Ls = map(lattice_extend, uc_extend) do lv, ucv
        scaling = lv ./ ucv
        if all(x -> x ≈ scaling[1], scaling)
            return scaling[1]
        else
            @warn "Irregular scaling of unitcell vectors $scaling"
            NaN
        end
    end

    return tuple(Int64.(Ls)...)
end

@inline function neighbors(l::LatPhysLattice, directed::Val{true})
    ((LatPhysBase.from(b), LatPhysBase.to(b)) for b in LatPhysBase.bonds(l.lattice))
end
@inline function neighbors(l::LatPhysLattice, directed::Val{false})
    (
        (LatPhysBase.from(b), LatPhysBase.to(b))
        for b in LatPhysBase.bonds(l.lattice)
        if LatPhysBase.from(b) < LatPhysBase.to(b)
    )
end
@inline function neighbors(l::LatPhysLattice, site_index::Integer)
    l.bonds[site_index]
end

@inline neighbors_lookup_table(l::LatPhysLattice) = copy(l.neighs)

positions(lattice::LatPhysLattice) = point.(sites(lattice.lattice))
lattice_vectors(l::LatPhysLattice) = latticeVectors(l.lattice)

function reciprocal_vectors(lattice::AbstractLattice, L)
    LatPhysBase.latticeVectors(
        LatticePhysics.getReciprocalUnitcell(
            LatPhysBase.unitcell(lattice.lattice)
        )
    )
end

################################################################################
### Saving & Loading
################################################################################

function save_lattice(file::JLDFile, lattice::LatPhysLattice, entryname::String)
    write(file, entryname * "/VERSION", 0)
    write(file, entryname * "/tag", "LatPhysLattice")
    _save_lattice(file, lattice.lattice, entryname * "/lattice")
    write(file, entryname * "/neighs", lattice.neighs)
    write(file, entryname * "/bonds", lattice.bonds)
    nothing
end
function _save_lattice(file::JLDFile, lattice::LatPhysBase.AbstractLattice, entryname::String)
    write(file, entryname * "/type", typeof(lattice))
    write(file, entryname * "/lv/N", length(lattice.lattice_vectors))
    for (i, v) in enumerate(lattice.lattice_vectors)
        write(file, entryname * "/lv/$i", v)
    end

    write(file, entryname * "/sites/type", typeof(lattice.sites[1])) # eltype of lattice
    write(file, entryname * "/sites/points", [s.point for s in lattice.sites])
    write(file, entryname * "/sites/labels", [s.label for s in lattice.sites])

    write(file, entryname * "/bonds/type", typeof(lattice.bonds[1])) # eltype
    write(file, entryname * "/bonds/froms",  [b.from for b in lattice.bonds])
    write(file, entryname * "/bonds/tos",    [b.to for b in lattice.bonds])
    write(file, entryname * "/bonds/labels", [b.label for b in lattice.bonds])
    write(file, entryname * "/bonds/wraps",  [b.wrap for b in lattice.bonds])

    save_unitcell(file, lattice.unitcell, entryname * "/unitcell")
    nothing
end
function save_unitcell(file::JLDFile, uc::LatPhysBase.AbstractUnitcell, entryname::String)
    write(file, entryname * "/type", typeof(uc))
    write(file, entryname * "/lv", uc.lattice_vectors)

    write(file, entryname * "/sites/type", typeof(uc.sites[1])) # eltype of lattice
    write(file, entryname * "/sites/points", [s.point for s in uc.sites])
    write(file, entryname * "/sites/labels", [s.label for s in uc.sites])

    write(file, entryname * "/bonds/type", typeof(uc.bonds[1])) # eltype
    write(file, entryname * "/bonds/froms",  [b.from  for b in uc.bonds])
    write(file, entryname * "/bonds/tos",    [b.to    for b in uc.bonds])
    write(file, entryname * "/bonds/labels", [b.label for b in uc.bonds])
    write(file, entryname * "/bonds/wraps",  [b.wrap  for b in uc.bonds])
    nothing
end

function _load(data, ::Val{:LatPhysLattice})
    @assert data["VERSION"] == 0
    LatPhysLattice(load_lattice(data["lattice"]), data["neighs"], data["bonds"])
end
function load_lattice(data)
    sites = load_sites(data["sites"])
    bonds = load_bonds(data["bonds"])
    lattice_vectors = [data["lv"]["$i"] for i in 1:data["lv"]["N"]]
    unitcell = load_unitcell(data["unitcell"])

    data["type"](lattice_vectors, sites, bonds, unitcell)
end
function load_unitcell(data)
    sites = load_sites(data["sites"])
    bonds = load_bonds(data["bonds"])
    lattice_vectors = data["lv"]
    data["type"](lattice_vectors, sites, bonds)
end
load_sites(data) = data["type"].(data["points"], data["labels"])
load_bonds(data) = data["type"].(data["froms"], data["tos"], data["labels"], data["wraps"])

to_tag(::Type{<: LatPhysLattice}) = Val(:LatPhysLattice)

################################################################################
### Extensions
################################################################################

function to_angle(v)
    angle = acos(v[1] / norm(v))
    (v[2] ≥ 0) ? angle : -angle
end

function get_sorted_corners(lattice::LatPhysLattice)
    bz = LatticePhysics.getBrillouinZone(
        LatticePhysics.getReciprocalUnitcell(
            LatticePhysics.unitcell(lattice.lattice)
        )
    )
    cs = sort(LatticePhysics.corners(bz), by = to_angle)
end

################################################################################
### Deprecated
################################################################################

function DistanceMask(lattice::LatPhysLattice)
    wrap = generate_combinations(latticeVectors(lattice.lattice))
    VerboseDistanceMask(lattice, wrap)
end

function directions(mask::VerboseDistanceMask, lattice::LatPhysLattice)
    pos = positions(lattice)
    dirs = [pos[src] - pos[trg] for (src, trg) in first.(mask.targets)]
    wrap = generate_combinations(latticeVectors(lattice.lattice))
    map(dirs) do _d
        d = round.(_d .+ wrap[1], digits=8)
        for v in wrap[2:end]
            new_d = round.(_d .+ v, digits=8)
            if directed_norm(new_d, 1e-6) < directed_norm(d, 1e-6)
                d .= new_d
            end
        end
        d
    end
end


text:
export tap_on_complete

import Base: show

"""
    tap_on_complete(tapFn::F) where { F <: Function }

Creates a tap operator, which performs a side effect
for only complete emission on the source Observable, but return an Observable that is identical to the source.

# Arguments
- `tapFn::Function`: side-effect tap function with `() -> Nothing` signature

# Producing

Stream of type `<: Subscribable{L}` where `L` refers to type of source stream

# Examples
```jldoctest
using Rocket

source = from([ 1, 2, 3 ])
subscribe!(source |> tap_on_complete(() -> println("Complete event received")), logger())
;

# output

[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Data: 3
[LogActor] Completed
Complete event received

```

See also: [`tap_on_subscribe`](@ref), [`tap`](@ref), [`logger`](@ref)
"""
tap_on_complete(tapFn::F) where { F <: Function } = TapOnCompleteOperator{F}(tapFn)

struct TapOnCompleteOperator{F} <: InferableOperator
    tapFn :: F
end

function on_call!(::Type{L}, ::Type{L}, operator::TapOnCompleteOperator{F}, source) where { L, F }
    return proxy(L, source, TapOnCompleteProxy{F}(operator.tapFn))
end

operator_right(operator::TapOnCompleteOperator, ::Type{L}) where L = L

struct TapOnCompleteProxy{F} <: ActorProxy
    tapFn :: F
end

actor_proxy!(::Type{L}, proxy::TapOnCompleteProxy{F}, actor::A) where { L, A, F } = TapOnCompleteActor{L, A, F}(proxy.tapFn, actor)

struct TapOnCompleteActor{L, A, F} <: Actor{L}
    tapFn :: F
    actor :: A
end

on_next!(actor::TapOnCompleteActor{L}, data::L) where L = next!(actor.actor, data)
on_error!(actor::TapOnCompleteActor, err)               = error!(actor.actor, err)
on_complete!(actor::TapOnCompleteActor)                 = begin complete!(actor.actor); actor.tapFn(); end

Base.show(io::IO, ::TapOnCompleteOperator)         = print(io, "TapOnCompleteOperator()")
Base.show(io::IO, ::TapOnCompleteProxy)            = print(io, "TapOnCompleteProxy()")
Base.show(io::IO, ::TapOnCompleteActor{L}) where L = print(io, "TapOnCompleteActor($L)")


text:
using Documenter, NewtonMethods

makedocs(modules=NewtonMethods,doctest=true)

deploydocs(deps = Deps.pip("mkdocs","python-markdown-math"),
    repo ="github.com/sonVishal/NewtonMethods.jl.git",
    julia = "0.4.7",
    osname = "linux")


text:
# shamelessly copied from Distributions.jl
macro check_args(D, cond)
    quote
        if !($(esc(cond)))
            throw(ArgumentError(string(
                $(string(D)), ": the condition ", $(string(cond)), " is not satisfied.")))
        end
    end
end