text:
"""
    load_genome_data(chr; directory="./1000genomes/parsed")

Read parsed plain-text genome data and return it as a boolean `SparseMatrixCSC`, where the
`i, j`-th entry indicates if there is a mutation in the `j`-th position of the genom of the `i`-th
subject (i.e., rows correspond to subjects and columns to genome positions).
"""
function load_genome_data(chr; directory="./1000genomes/parsed")
    filename = joinpath(directory, "chr$(chr).txt")
    df = DataFrame(CSV.File(filename))
    sparse(df[:, 2], df[:, 1], true)
end

"""

Read parsed plain-text genome data for all chromosomes and write it to a single matrix in order.
"""
function genome_to_hdf5(filename="./1000genomes/parsed/1000genomes.h5")
    println("Parsing chromosome 1")
    X = load_genome_data(1)
    h5writecsc(filename, "X", X)
    for chr in 2:22
        GC.gc() # force garbage collection to avoid running out of memory
        println("Parsing chromosome $chr")
        h5appendcsc(filename, "X", load_genome_data(chr))
    end
end

"""

Write a randomly generated sparse matrix.
"""
function write_sprand_matrix(m, n, p, filename, name; nblocks::Integer=ceil(Int, n/10000), overwrite=true)
    firstcol = 1
    lastcol = floor(Int, 1/nblocks*n)
    println("Block 1 / $nblocks")
    h5open(filename, "cw") do fid
        h5writecsc(filename, name, sprand(m, lastcol-firstcol+1, p); overwrite)
        for i in 2:nblocks
            firstcol = floor(Int, (i-1)/nblocks*n+1)
            lastcol = floor(Int, i/nblocks*n)
            println("Block $i / $nblocks")
            h5appendcsc(filename, name, sprand(m, lastcol-firstcol+1, p))
            GC.gc()
        end
    end
    return
end


text:
using ATA
using JuMP
# add Cbc by running import Pkg; Pkg.add("Cbc")
using Cbc
# cd("where your input files are")
cd("examples")
# 1. Start ATA and add file with custom settings (Needed)
ata_model = start_ata(;
    settings_file = "settings_ata_maximin.jl",
    bank_file = "data/bank.csv",
    bank_delim = ";",
)

# 2. Add categorical constraints (Optional)
add_constraints!(ata_model; constraints_file = "constraints.csv", constraints_delim = ";")

# 3. Add overlap maxima (Optional)
add_overlap!(ata_model; overlap_file = "overlap_matrix.csv", overlap_delim = ";")

# 4. Add objective function (Optional)
add_obj_fun!(ata_model)
print_infos(ata_model)

# Assembly settings

# Set the solver, "siman" for simulated annealing, "jump" for MILP solver.
solver = "jump"

# MILP (Not suggested for large scale ATA)
# Select the solver, Cbc as open-source is a good option.
optimizer_constructor = "Cbc"
# #Optimizer attributes
optimizer_attributes = [("seconds", 100), ("logLevel", 1)]

# 5. assemble
assemble!(
    ata_model;
    solver = solver,
    optimizer_attributes = optimizer_attributes,
    optimizer_constructor = optimizer_constructor,
)

# All the settings and outputs from optimization are in ata_model object.
# See the struct in ATA.jl to understand how to retrieve all the information.
# A summary of the resulting tests is available in results_folder/results.txt
# If siman is chosen, the optimality and feasibility of the best neighbourhood
# is reported in "results/results_ata.jl"

print_results(ata_model; results_folder = "results")

#]add https://github.com/giadasp/ATAPlot.jl
using ATAPlot
plot_results(ata_model; plots_folder = "plots")

ext: jl


text:
using Robotlib, SymPy
Tworld2base = I4
# """
# `τ = gravity{P}(q::VecOrMat{P}, rm, m, dh::DH, Tworld2base = I4)`
# calculates the gravity torques given a joint angles `q`, (distances to center of mass × mass) `rm` and masses `m`. The DH-parameters for the robot must also be provided
# """
function gravity(q::VecOrMat{P}, rm, m, dh::DH, Tworld2base = I4) where P
    Tn  = dh2Tn(dh,q)
    gravity(q, rm, m, Tn, Tworld2base)
end

"""
    gravity{P}(q, rm, m, Tn, Tworld2base)

`q` is the joint coordinates, `rm` are distances to center of masses times masses, `m` are the masses and `Tn` are the transformations between consecutive joints.
"""
function gravity(q::VecOrMat{P}, rm, m, Tn::AbstractArray, Tworld2base = I4) where P
    n   = size(Tn,3)-1
    # Find elements that are ≈ 0 and make them zero, smiplifies a lot!
    if P == SymPy.Sym
        Tnz = Bool[round(Float64(evalf(subs(Tn[i,j,k],(q[1],1),(q[2],1),(q[3],1),(q[4],1),(q[5],1),(q[6],1),(q[7],1)))), digits=10) == 0 for i=1:4, j=1:4, k=1:n+1]
        Tn[Tnz] = 0
        # Replace the numerical value for π/2 by the SymPi constant PI, this allows for nice simplifications
        Tn = Sym[subs(Tn[i,j,k],("1.5707963267948966",SymPy.PI/2))  for i=1:4, j=1:4, k=1:n+1]
    end

    Ti = zeros(typeof(q[1]),4,4,n+1)
    Tn = cat(3,Tworld2base,Tn)
    Ti[:,:,1] = Tn[:,:,1]

    for i = 2:n
        Ti[:,:,i] = Ti[:,:,i-1]*Tn[:,:,i] # Cumulative forward kinematics transformations
    end

    gv         = [0, 0, -9.82]
    τ          = Array(P,3,n+1) # Extend by one to run everything in loop
    τ[:,end]   = 0
    τhat       = Array(P,n)
    force      = Array(P,3,n) # All forces are given in base frame

    for i = n:-1:1
        Ri     = Tn[1:3,1:3,i+1] # Rotation from i to i+1
        gi     = Ti[1:3,1:3,i]'gv # Gravity vector in i
        τ[:,i] = skew(rm[:,i]) * gi + Ri*τ[:,i+1] # Torque in i is gravity torque + rotated torque from i+1
        τi     = τ[:,i]
        for k = (i+1):n # I can not get it to work with accumulated forces which would be much faster
            rik = (trinv(Ti[:,:,i])*Ti[1:4,4,k])[1:3] # Vector from i to k, given in i
            τi += skew(rik)*Ti[1:3,1:3,i]'force[:,k] # Torque around i, given in i
        end

        τhat[i]     = τi[3] # The motor torque is the local torque around z
        force[:,i]  = gv*m[i]
    end
    return -τhat # Minus sign to match sensor torque
end

function create_gravmatrix(dh::DH)

    n_joints = size(dh.dhpar,1)

    rm       = Sym[symbols("rm$i$j",real=true) for i = 1:3, j=1:n_joints]
    m        = Sym[symbols("m$j",real=true) for j=1:n_joints]
    q        = Sym[symbols("q$j",real=true) for j=1:n_joints]
    tau      = Sym[symbols("tau$j",real=true) for j=1:n_joints]

    tauhat = gravity(q, rm, m, dh, Tbase);
    w      = [rm[:];m[:]];

    AA     = Array(Sym,n_joints, 4n_joints);
    AAz    = falses(size(AA))

    eq     = tauhat;

    println("Finding non-zero coefficients")
    for i = 1:n_joints
        for j = 1:size(AA,2)
            AA[i,j] = (coeff(expand(eq[i]),w[j])) # This is a nice place to put simplify
            AAz[i,j] = abs(evalf(subs(AA[i,j],(q[1],1),(q[2],1),(q[3],1),(q[4],1),(q[5],1),(q[6],1),(q[7],1)))) < 1e-5
        end
        print(i)
    end
    println("")

    # non_identifiable    = (sum(AA,1) .== 0)[:]
    # # non_identifiable = all(AAz,1)[:]
    # w                   = w[!non_identifiable];
    # AA                  = AA[:,!non_identifiable];
    n_joints,n_params   = size(AA)

    println("Trying to find expressions common in all terms")
    change,res  = cse(AA[:]);
    res         = reshape(res,size(AA));

    println("Printing results to file")
    fid = open("gravityFridaLS.jl","w");
    println(fid,"@fastmath function gravityFridaLS(q)")
    for i = 1:n_joints
        println(fid,"q$i = q[$i]")
    end
    for i = 0:length(change)-1
        sei = change[i+1][2]
        println(fid,"x$i = $sei")
    end

    println(fid,"A = zeros($n_joints,$n_params)")

    for j = 1:n_joints
        for p = 1:n_params
            if res[j,p] != 0 #&& !AAz[j,p]
                println(fid,"A[$j,$p] = $(res[j,p])")
            end
        end
    end
    println(fid,"return A")
    println(fid,"end")
    close(fid)
    println("Done")
end
# dh       = DHYuMi()
# baseAnglesLeft = [-0.63 , 0.95 , -0.18]
# Rbase          = rpy2R(baseAnglesLeft,"xyz")
# Tbase          = I4
# Tbase[1:3,1:3] = Rbase

"""
    create_gravmodel(filename, r,m, dh, Tbase=I4)
Given known masses and vectors to center of masses, this function creates a new function to calculate the gravity torque.
If only `r*m` is known, call the function with `m=1` and `r=rm`
`r ∈ ℜ(3 × n_joints)`
`m ∈ ℜ(n_joints)`

The created function is written to a file with and is called like `τ = gravity(q)`
"""
function create_gravmodel(filename, r::AbstractMatrix,m, dh::DH, Tbase=I4)
    n_joints       = size(dh.dhpar,1)
    rm             = (r.*m)'
    q              = Sym[symbols("q$j",real=true) for j=1:n_joints]
    tauhat         = gravity(q, rm, m, dh, Tbase);
    n_joints       = size(r,1)

    println("Trying to find expressions common in all terms")
    change,res  = cse(tauhat);

    println("Printing results to file")
    fid = open(filename*".jl","w");
    println(fid,"@fastmath function gravity(q)")
    for i = 1:n_joints
        println(fid,"q$i = q[$i]")
    end
    for i = 0:length(change)-1
        sei = change[i+1][2]
        println(fid,"x$i = $sei")
    end

    println(fid,"tau = zeros($n_joints)")

    for j = 1:n_joints
            if res[j] != 0 #&& !AAz[j,p]
                println(fid,"tau[$j] = $(res[j])")
            end
    end
    println(fid,"return tau")
    println(fid,"end")
    close(fid)
    println("Done")
end


text:
export ProjMPO,
       lproj,
       rproj,
       product

mutable struct ProjMPO
  lpos::Int
  rpos::Int
  nsite::Int
  H::MPO
  LR::Vector{ITensor}
  ProjMPO(H::MPO) = new(0,length(H)+1,2,H,fill(ITensor(),length(H)))
end

nsite(pm::ProjMPO) = pm.nsite
length(pm::ProjMPO) = length(pm.H)

function lproj(pm::ProjMPO)::ITensor
  (pm.lpos <= 0) && return ITensor()
  return pm.LR[pm.lpos]
end

function rproj(pm::ProjMPO)::ITensor
  (pm.rpos >= length(pm)+1) && return ITensor()
  return pm.LR[pm.rpos]
end

function product(pm::ProjMPO,
                 v::ITensor)::ITensor
  Hv = v
  if isnull(lproj(pm))
    if !isnull(rproj(pm))
      Hv *= rproj(pm)
    end
    for j=pm.rpos-1:-1:pm.lpos+1
      Hv *= pm.H[j]
    end
  else #if lproj is not null
    Hv *= lproj(pm)
    for j=pm.lpos+1:pm.rpos-1
      Hv *= pm.H[j]
    end
    if !isnull(rproj(pm))
      Hv *= rproj(pm)
    end
  end
  return noprime(Hv)
end

function Base.eltype(pm::ProjMPO)
  elT = eltype(pm.H[pm.lpos+1])
  for j = pm.lpos+2:pm.rpos-1
    elT = promote_type(elT,eltype(pm.H[j]))
  end
  if !isnull(lproj(pm))
    elT = promote_type(elT,eltype(lproj(pm)))
  end
  if !isnull(rproj(pm))
    elT = promote_type(elT,eltype(rproj(pm)))
  end
  return elT
end

(pm::ProjMPO)(v::ITensor) = product(pm,v)

function size(pm::ProjMPO)::Tuple{Int,Int}
  d = 1
  if !isnull(lproj(pm))
    for i in inds(lproj(pm))
      plev(i) > 0 && (d *= dim(i))
    end
  end
  for j=pm.lpos+1:pm.rpos-1
    for i in inds(pm.H[j])
      plev(i) > 0 && (d *= dim(i))
    end
  end
  if !isnull(rproj(pm))
    for i in inds(rproj(pm))
      plev(i) > 0 && (d *= dim(i))
    end
  end
  return (d,d)
end

function makeL!(pm::ProjMPO,
                psi::MPS,
                k::Int)
  while pm.lpos < k
    ll = pm.lpos
    if ll <= 0
      pm.LR[1] = psi[1]*pm.H[1]*dag(prime(psi[1]))
      pm.lpos = 1
    else
      pm.LR[ll+1] = pm.LR[ll]*psi[ll+1]*pm.H[ll+1]*dag(prime(psi[ll+1]))
      pm.lpos += 1
    end
  end
end

function makeR!(pm::ProjMPO,
                psi::MPS,
                k::Int)
  N = length(pm.H)
  while pm.rpos > k
    rl = pm.rpos
    if rl >= N+1
      pm.LR[N] = psi[N]*pm.H[N]*dag(prime(psi[N]))
      pm.rpos = N
    else
      pm.LR[rl-1] = pm.LR[rl]*psi[rl-1]*pm.H[rl-1]*dag(prime(psi[rl-1]))
      pm.rpos -= 1
    end
  end
end

function position!(pm::ProjMPO,
                   psi::MPS,
                   pos::Int)
  makeL!(pm,psi,pos-1)
  makeR!(pm,psi,pos+nsite(pm))

  #These next two lines are needed
  #when moving lproj and rproj backward
  pm.lpos = pos-1
  pm.rpos = pos+nsite(pm)
end

# Return a "noise tensor" a la PRB 72, 180403.
function noisetensor(M   :: MPS,
                     PH  :: ProjMPO,
                     b   :: Int,
                     noise :: Float64,
                     dir :: String)

  if (PH.lpos < 1 && dir=="fromleft") || (PH.rpos > length(M) &&  dir=="fromright")
    return ITensor(0)
  end

  if dir=="fromleft"
    nt = lproj(PH)
    ind = commonindex(nt, PH.H[PH.lpos + 1])
  elseif dir=="fromright"
    nt = rproj(PH)
    ind = commonindex(nt, PH.H[PH.rpos - 1])
  else throw(ArgumentError("In noisetensor(), no dir = $dir supported. Use fromleft or fromright."))
  end

  #Do we want complex to break TR symmetry?
  nt *= randomITensor(ind)
  nt *= noise
  return nt
end


text:
using HttpServer

function fibonacci(n)
  if n == 1 return 1 end
  if n == 2 return 1 end
  prev = BigInt(1)
  pprev = BigInt(1)
  for i=3:n
    curr = prev + pprev
    pprev = prev
    prev = curr
  end
  return prev
end

http = HttpHandler() do req::Request, res::Response
    m = match(r"^/fibo/(\d+)/?$",req.resource)
    if m == nothing return Response(404) end
    number = BigInt(m.captures[1])
    if number < 1 || number > 100_000 return Response(500) end
    return Response(string(fibonacci(number)))
end

http.events["error"]  = (client, err) -> println(err)
http.events["listen"] = (port)        -> println("Listening on $port...")

server = Server(http)
run(server, 8000)


text:
module TrafficWeavingPlanner

using RobotOS
using StaticArrays

@rosimport std_msgs.msg: Float32MultiArray, Float64MultiArray,
                         Int32MultiArray, Int64MultiArray,
                         MultiArrayLayout, MultiArrayDimension
@rosimport traffic_weaving_prediction.msg: prediction_input, prediction_output
@rosimport vtd_interface.msg: planner_result, rdb_object_state, scenario_start
rostypegen()

import std_msgs.msg: Float32MultiArray, Float64MultiArray,
                     Int32MultiArray, Int64MultiArray,
                     MultiArrayLayout, MultiArrayDimension
import traffic_weaving_prediction.msg: prediction_input, prediction_output
import vtd_interface.msg: planner_result, rdb_object_state, scenario_start

reversedims(A::Array) = permutedims(A, ndims(A):-1:1) # @generated probably the "right" way

for MAType in (Float32MultiArray, Float64MultiArray, Int32MultiArray, Int64MultiArray)
    @eval function Base.convert(::Type{$MAType}, A::Array{<:Real})
        Arev = reversedims(A)
        msg = $MAType()
        msg.layout.dim = [MultiArrayDimension("",d,s) for (d,s) in zip(size(A), reverse(cumprod([size(Arev)...])))]
        msg.data = Arev[:]
        msg
    end

    @eval function Base.convert(::Type{Array{T}}, msg::$MAType) where {T}
        if isempty(msg.data)
            T[]
        else
            reversedims(reshape(msg.data, (Int(d.size) for d in msg.layout.dim[end:-1:1])...))
        end
    end
end

mutable struct PredictionOutput
    stale::Bool
    y
    z
    r
end

function prediction_callback(msg::prediction_output, pred_container::PredictionOutput)
    pred_container.stale = false
    pred_container.y = Array{Float32}(msg.y)
    pred_container.z = Array{Float32}(msg.z)
    pred_container.r = Array{Float32}(msg.r)
end

init_node("julia_laneswap", anonymous=true)

const pred_container = PredictionOutput(true,nothing,nothing,nothing)
const publisher = Publisher{prediction_input}("prediction_input", queue_size=10)
const subscriber = Subscriber{prediction_output}("prediction_output", prediction_callback, (pred_container,), queue_size=10)

# rossleep(Duration(5))

function predict!(car1, car2, extras, traj_lengths, car1_future, sample_ct,
                  output::PredictionOutput=pred_container, timeout=1000)
    # println("predict! function running")
    output.stale = true
    publish(publisher, prediction_input(Float32MultiArray(car1),
                                        Float32MultiArray(car2),
                                        Float32MultiArray(extras),
                                        Int32MultiArray(traj_lengths),
                                        Float32MultiArray(car1_future),
                                        Float32MultiArray(Float32[]),
                                        Float32MultiArray(Float32[]),
                                        Int32(sample_ct)))
    # println("published data")
    for i in 1:timeout
        rossleep(Duration(0.001))
        output.stale || break
    end
    output
end

function predict!(car1, car2, extras, traj_lengths, car1_future_x, car1_future_y, sample_ct,
                  output::PredictionOutput=pred_container, timeout=1000)
    # println("predict! function running")
    output.stale = true
    publish(publisher, prediction_input(Float32MultiArray(car1),
                                        Float32MultiArray(car2),
                                        Float32MultiArray(extras),
                                        Int32MultiArray(traj_lengths),
                                        Float32MultiArray(Float32[]),
                                        Float32MultiArray(car1_future_x),
                                        Float32MultiArray(car1_future_y),
                                        Int32(sample_ct)))
    # println("published data")
    for i in 1:timeout
        rossleep(Duration(0.001))
        output.stale || break
    end
    output
end

####

mutable struct HumanCarStates
    t::Vector{Float32}
    X::Vector{SVector{6,Float32}}
end

function human_state_callback(msg::rdb_object_state, human_states::HumanCarStates)
    if isempty(human_states.t) || msg.simTime > human_states.t[end] + 1e-4
        push!(human_states.t, Float32(msg.simTime))
        push!(human_states.X, SVector{6,Float32}(msg.base.pos.x,
                                                 msg.base.pos.y,
                                                 msg.ext.speed.x,
                                                 msg.ext.speed.y,
                                                 msg.ext.accel.x,
                                                 msg.ext.accel.y))
    end
end

const human_states = HumanCarStates(Float32[], SVector{6,Float32}[])
const human_state_sub = Subscriber{rdb_object_state}("/object_state/1", human_state_callback, (human_states,), queue_size=10)
const planner_result_pub = Publisher{planner_result}("/planner_path", queue_size=10)
const scenario_setup_pub = Publisher{scenario_start}("/scenario_setup", queue_size=10)

function clear_human_states()
    human_states.t = Float32[]
    human_states.X = SVector{6,Float32}[]
    nothing
end

reset_simulator() = publish(planner_result_pub, planner_result(Float32[0], Float32[0], Float32[0]))
function setup_simulator(;xh=-136.0, yh=-1.83, vh=30.0, xr=-136.0, yr=-6.09, vr=30.0)
    publish(scenario_setup_pub, scenario_start(xh, yh, vh, xr, yr, vr))
end
publish_robot_path_plan(t, x, y) = publish(planner_result_pub, planner_result(t, x, y))

@spawn spin()

end # module

max_issues_repo_name: codingblazes/NHumanModeling

max_forks_repo_name: codingblazes/NHumanModeling


text:
export ImmutablePolynomial

"""
    ImmutablePolynomial{T, X, N}(coeffs::AbstractVector{T})

Construct an immutable (static) polynomial from its coefficients
`a₀, a₁, …, aₙ`,
lowest order first, optionally in terms of the given variable `x`
where `x` can be a character, symbol, or string.

If ``p = a_n x^n + \\ldots + a_2 x^2 + a_1 x + a_0``, we construct
this through `ImmutablePolynomial((a_0, a_1, ..., a_n))` (assuming
`a_n ≠ 0`). As well, a vector or number can be used for construction.

The usual arithmetic operators are overloaded to work with polynomials
as well as with combinations of polynomials and scalars. However,
operations involving two non-constant polynomials of different variables causes an
error. Unlike other polynomials, `setindex!` is not defined for `ImmutablePolynomials`.

As the degree of the polynomial (`+1`) is a compile-time constant,
several performance improvements are possible. For example, immutable
polynomials can take advantage of faster polynomial evaluation
provided by `evalpoly` from Julia 1.4; similar methods are also used
for addtion and multiplication.

However, as the degree is included in the type, promotion between
immutable polynomials can not promote to a common type. As such, they
are precluded from use in rational functions.

!!! note
    `ImmutablePolynomial` is not axis-aware, and it treats `coeffs` simply as a list of coefficients with the first
    index always corresponding to the constant term.

# Examples

```jldoctest
julia> using  Polynomials

julia> ImmutablePolynomial((1, 0, 3, 4))
ImmutablePolynomial(1 + 3*x^2 + 4*x^3)

julia> ImmutablePolynomial((1, 2, 3), :s)
ImmutablePolynomial(1 + 2*s + 3*s^2)

julia> one(ImmutablePolynomial)
ImmutablePolynomial(1.0)
```

!!! note
    This was modeled after [StaticUnivariatePolynomials](https://github.com/tkoolen/StaticUnivariatePolynomials.jl) by `@tkoolen`.

"""
struct ImmutablePolynomial{T, X, N} <: StandardBasisPolynomial{T, X}
    coeffs::NTuple{N, T}
    function ImmutablePolynomial{T,X,N}(coeffs::NTuple{N,T}) where {T, X, N}
        N == 0 && return new{T,X, 0}(coeffs)
        iszero(coeffs[end]) &&  throw(ArgumentError("Leading  term must  be  non-zero"))
        new{T,X,N}(coeffs)
    end
end

@register ImmutablePolynomial

## Various interfaces
## Abstract Vector coefficients
function ImmutablePolynomial{T,X, N}(coeffs::AbstractVector{T})  where {T,X, N}
    cs = NTuple{N,T}(coeffs[i] for i ∈ firstindex(coeffs):N)
    ImmutablePolynomial{T, X, N}(cs)

end

function ImmutablePolynomial{T,X, N}(coeffs::AbstractVector{S})  where {T,X, N, S}
    cs = NTuple{N,T}(coeffs[i] for i ∈ firstindex(coeffs):N)
    ImmutablePolynomial{T, X, N}(cs)

end

function ImmutablePolynomial{T,X}(coeffs::AbstractVector{S})  where {T,X,S}
    R = promote_type(T,S)

    if Base.has_offset_axes(coeffs)
        @warn "ignoring the axis offset of the coefficient vector"
    end
    N = findlast(!iszero, coeffs)
    N == nothing && return ImmutablePolynomial{R,X,0}(())
    N′ = N + 1 - firstindex(coeffs)
    ImmutablePolynomial{T, X, N′}([coeffs[i] for i ∈ firstindex(coeffs):N])
end

## -- Tuple arguments
function ImmutablePolynomial{T,X}(coeffs::Tuple)  where {T,X}
    N = findlast(!iszero, coeffs)
    N == nothing && return zero(ImmutablePolynomial{T,X})
    ImmutablePolynomial{T,X,N}(NTuple{N,T}(coeffs[i] for i in 1:N))
end

ImmutablePolynomial{T}(coeffs::Tuple, var::SymbolLike=:x)  where {T} = ImmutablePolynomial{T,Symbol(var)}(coeffs)

function ImmutablePolynomial(coeffs::Tuple, var::SymbolLike=:x)
    cs = NTuple(promote(coeffs...))
    T = eltype(cs)
    ImmutablePolynomial{T, Symbol(var)}(cs)
end

##
## ----
##
# overrides from common.jl due to  coeffs being non mutable, N in type parameters

Base.copy(p::P) where {P <: ImmutablePolynomial} = P(coeffs(p))

# degree, isconstant
degree(p::ImmutablePolynomial{T,X, N}) where {T,X,N} = N - 1 # no trailing zeros
isconstant(p::ImmutablePolynomial{T,X,N}) where {T,X,N}  = N <= 1

Base.setindex!(p::ImmutablePolynomial, val,  idx::Int) = throw(ArgumentError("ImmutablePolynomials are immutable"))

for op in [:isequal, :(==)]
    @eval function Base.$op(p1::ImmutablePolynomial{T,N}, p2::ImmutablePolynomial{S,M}) where {T,N,S,M}
        check_same_variable(p1,p2) || return false
        p1s, p2s = coeffs(p1), coeffs(p2)
        (N == M  && $op(p1s,p2s)) &&  return  true
        n1 = findlast(!iszero, p1s) # now trim out zeros
        n2 = findlast(!iszero, p2s)
        (n1 == nothing && n2 == nothing) && return true
        (n1 == nothing || n2  == nothing) && return false
        $op(p1s[1:n1],p2s[1:n2]) &&  return true
        false
    end
end

# in common.jl these call chop! and truncate!
function Base.chop(p::ImmutablePolynomial{T,X};
              rtol::Real = Base.rtoldefault(real(T)),
                   atol::Real = 0)  where {T,X}
    ps = chop(p.coeffs; rtol=rtol, atol=atol)
    return ImmutablePolynomial{T,X}(ps)
end

function Base.truncate(p::ImmutablePolynomial{T,X};
                  rtol::Real = Base.rtoldefault(real(T)),
                  atol::Real = 0)  where {T,X}
    ps = truncate(p.coeffs; rtol=rtol, atol=atol)
    ImmutablePolynomial{T,X}(ps)
end

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

## Addition
# scalar ops
function Base.:+(p::P, c::S) where {T, X, N, P <: ImmutablePolynomial{T,X,N}, S<:Number}
    R = promote_type(T,S)

    iszero(c) && return ImmutablePolynomial{R,X,N}(convert(NTuple{N,R},p.coeffs))
    N == 0 && return ImmutablePolynomial{R,X,1}(NTuple{1,R}(c))
    N == 1 && return ImmutablePolynomial((p[0]+c,), X)

    cs = ⊕(P, convert(NTuple{N,R},p.coeffs), NTuple{1,R}(c))
    q = ImmutablePolynomial{R,X,N}(cs)

    return q

end

Base.:-(p::ImmutablePolynomial{T,X,N}) where {T,X,N} = ImmutablePolynomial{T,X,N}(.-p.coeffs)

function Base.:+(p1::P, p2::Q) where {T,X,N,P<:ImmutablePolynomial{T,X,N},
                                      S,  M,Q<:ImmutablePolynomial{S,X,M}}

    R = promote_type(T,S)
    P′ = ImmutablePolynomial{R,X}
    if  N == M
        cs = ⊕(P, p1.coeffs, p2.coeffs)
        return P′(R.(cs))
    elseif N < M
        cs = ⊕(P, p2.coeffs, p1.coeffs)
        return P′{M}(R.(cs))
    else
        cs = ⊕(P, p1.coeffs, p2.coeffs)
        return P′{N}(R.(cs))
    end

end

## multiplication

function scalar_mult(p::ImmutablePolynomial{T,X,N}, c::S) where {T, X,N, S <: Number}
    R = eltype(p[0] * c * 0)
    (N == 0  || iszero(c)) && return zero(ImmutablePolynomial{R,X})
    cs = p.coeffs .* c
    return ImmutablePolynomial(cs, X)
 end

function scalar_mult(c::S, p::ImmutablePolynomial{T,X,N}) where {T, X,N, S <: Number}
    R = eltype(p[0] * c * 0)
    (N == 0  || iszero(c)) && return zero(ImmutablePolynomial{R,X})
    cs = p.coeffs .* c
    return ImmutablePolynomial(cs, X)
end

function Base.:/(p::ImmutablePolynomial{T,X,N}, c::S) where {T,X,N,S}
    R = eltype(one(T)/one(S))
    P = ImmutablePolynomial{R,X}
    (N == 0  || isinf(c)) && return zero(P)
    cs = p.coeffs ./ c
    iszero(cs[end]) ? P(cs) : P{N}(cs) # more performant to specify when N is known
end

function Base.:*(p1::ImmutablePolynomial{T,X,N}, p2::ImmutablePolynomial{S,X,M}) where {T,S,X,N,M}
    R = promote_type(T,S)
    P = ImmutablePolynomial{R,X}

    (iszero(N) || iszero(M)) && return zero(P)

    cs = ⊗(ImmutablePolynomial, p1.coeffs, p2.coeffs) #(p1.coeffs) ⊗ (p2.coeffs)
    iszero(cs[end]) ? P(cs) : P{N+M-1}(cs)  # more performant to specify when N is known

end

Base.to_power_type(p::ImmutablePolynomial{T,X,N}) where {T,X,N} = p

## more performant versions borrowed from StaticArrays
## https://github.com/JuliaArrays/StaticArrays.jl/blob/master/src/linalg.jl
LinearAlgebra.norm(q::ImmutablePolynomial{T,X,0}) where {T,X} = zero(real(float(T)))
LinearAlgebra.norm(q::ImmutablePolynomial) = _norm(q.coeffs)
LinearAlgebra.norm(q::ImmutablePolynomial, p::Real) = _norm(q.coeffs, p)

@generated function _norm(a::NTuple{N,T}) where {T, N}

    expr = :(abs2(a[1]))
    for j = 2:N
        expr = :($expr + abs2(a[$j]))
    end

    return quote
        $(Expr(:meta, :inline))
        @inbounds return sqrt($expr)
    end

end

_norm_p0(x) = iszero(x) ? zero(x) : one(x)
@generated function _norm(a::NTuple{N,T}, p::Real) where {T, N}
    expr = :(abs(a[1])^p)
    for j = 2:N
        expr = :($expr + abs(a[$j])^p)
    end

    expr_p1 = :(abs(a[1]))
    for j = 2:N
        expr_p1 = :($expr_p1 + abs(a[$j]))
    end

    return quote
        $(Expr(:meta, :inline))
        if p == Inf
            return mapreduce(abs, max, a)
        elseif p == 1
            @inbounds return $expr_p1
        elseif p == 2
            return norm(a)
        elseif p == 0
            return mapreduce(_norm_p0, +, a)
        else
            @inbounds return ($expr)^(inv(p))
        end
    end

end


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using CSV, CategoricalArrays, StatsPlots
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing4.28"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	df = CSV.read(datadir("companyData.csv"))
	years = levels(df.Year)
	data  = reshape(df.MarketCap, 5, 3)
end

# ╔═╡ fabeade4-2522-11eb-21d7-f160a9d3f5b7
begin
	p1 = groupedbar(years, data, bar_position=:stack)
	p2 = groupedbar(years, data, bar_position=:dodge)
	plot(p1, p2, bar_width=0.7, fill=[:blue :red :green], label=["A" "B" "C"],
		ylims=(0,6), xlabel="Year", ylabel="Market Cap (MM)",
		legend=:topleft, size=(800,400))
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing4.28"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═fabeade4-2522-11eb-21d7-f160a9d3f5b7
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
VERSION >= v"1.0.0"

module AlphaVantage

const alphavantage_api = "https://www.alphavantage.co/"

using ArgCheck
using DelimitedFiles
using HTTP
using HttpCommon
using JSON

function __init__()
    if !haskey(ENV, "ALPHA_VANTAGE_API_KEY")
        @warn "No API key found"
    end
end

include("utils.jl")
include("stock_time_series.jl")
include("digital_currency.jl")
include("foreign_exchange_currency.jl")
include("stock_technical_indicators.jl")
include("sector_performance.jl")
include("fundamentals.jl")

end # module


text:
for T ∈ (:Float16,:Int16,:UInt16)
    @eval @inline sizeequivalentfloat(::Type{$T}) = Float16
    @eval @inline sizeequivalentint(::Type{$T}) = Int16
    @eval @inline sizeequivalentuint(::Type{$T}) = UInt16
    @eval @inline Base.$T(v::SVec{W}) where {W} = SVec(vconvert(Vec{W,$T}, extract_data(v)))
end
for T ∈ (:Float32,:Int32,:UInt32)
    @eval @inline sizeequivalentfloat(::Type{$T}) = Float32
    @eval @inline sizeequivalentint(::Type{$T}) = Int32
    @eval @inline sizeequivalentuint(::Type{$T}) = UInt32
    @eval @inline Base.$T(v::SVec{W}) where {W} = SVec(vconvert(Vec{W,$T}, extract_data(v)))
end
for T ∈ (:Float64,:Int64,:UInt64)
    @eval @inline sizeequivalentfloat(::Type{$T}) = Float64
    @eval @inline sizeequivalentint(::Type{$T}) = Int64
    @eval @inline sizeequivalentuint(::Type{$T}) = UInt64
    @eval @inline Base.$T(v::SVec{W}) where {W} = SVec(vconvert(Vec{W,$T}, extract_data(v)))
end
@inline sizeequivalentfloat(::Type{T}, x::T) where {T<:FloatingTypes} = x
@inline sizeequivalentint(::Type{T}, x::T) where {T<:Signed} = x
@inline sizeequivalentuint(::Type{T}, x::T) where {T<:Unsigned} = x
@inline function sizeequivalentfloat(::Type{T}, x) where {T}
    convert(sizeequivalentfloat(T), x)
end
@inline function sizeequivalentint(::Type{T}, x) where {T}
    convert(sizeequivalentint(T), x)
end
@inline function sizeequivalentint(::Type{T}, x::Integer) where {T}
    x % sizeequivalentint(T)
end
@inline function sizeequivalentuint(::Type{T}, x) where {T}
    convert(sizeequivalentint(T), x)
end
@inline function sizeequivalentuint(::Type{T}, x::Integer) where {T}
    x % sizeequivalentuint(T)
end

# Element-wise access

@generated function vsetindex(v::_Vec{_W,T}, x::Number, ::Val{I}) where {_W,T,I}
    W = _W + 1
    @assert isa(I, Integer)
    1 <= I <= W || throw(BoundsError())
    typ = llvmtype(T)
    ityp = llvmtype(Int)
    vtyp = "<$W x $typ>"
    # decls = String[]
    instrs = String[]
    push!(instrs, "%res = insertelement $vtyp %0, $typ %1, $ityp $(I-1)")
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        llvmcall($(join(instrs, "\n")), Vec{$W,$T}, Tuple{Vec{$W,$T}, $T}, extract_data(v), $T(x))
    end
end

@generated function vsetindex(v::_Vec{_W,T}, x::Number, i::Int) where {_W,T}
    W = _W + 1
    typ = llvmtype(T)
    ityp = llvmtype(Int)
    vtyp = "<$W x $typ>"
    # decls = String[]
    instrs = String[]
    push!(instrs, "%res = insertelement $vtyp %0, $typ %2, $ityp %1")
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        @boundscheck 1 <= i <= $W || throw(BoundsError())
        llvmcall($(join(instrs, "\n")), Vec{$W,$T}, Tuple{Vec{$W,$T}, Int, $T}, extract_data(v), vsub(i, 1), $T(x))
    end
end

Base.@propagate_inbounds Base.setindex(v::AbstractStructVec{W,T}, x::Number, i) where {W,T} = @inbounds SVec{W}(vsetindex(extract_data(v), x, i % Int))
@inline Base.setindex(v::AbstractStructVec{W,T}, x::Number, ::Val{I}) where {W,T,I} = SVec{W}(vsetindex(estract_data(v), x, I))

@inline getvalindex(v::AbstractSIMDVector, ::Val{I}) where {I} = extract_data(v)[I].value
@inline Base.convert(::Type{SVec{W,T}}, v::Vec{W,T}) where {W,T} = SVec(v)
@inline Base.convert(::Type{SVec{W,T}}, x::SVec{W,T}) where {W,T} = x

@inline VectorizationBase.vbroadcast(::Val{W}, v::AbstractSIMDVector{W}) where {W} = v
@inline vconvert(::Type{_Vec{_W,T}}, s::T) where {_W, T <: Unsigned} = vbroadcast(_Vec{_W,T}, s)
@generated function vconvert(::Type{_Vec{_W,T1}}, v::_Vec{_W,T2}) where {_W,T1 <: FloatingTypes, T2 <: Signed}
    W = _W + 1
    typ1 = llvmtype(T1)
    typ2 = llvmtype(T2)
    vtyp1 = "<$W x $typ1>"
    vtyp2 = "<$W x $typ2>"
    instrs = """
    %res = sitofp $vtyp2 %0 to $vtyp1
    ret $vtyp1 %res
    """
    quote
        $(Expr(:meta, :inline))
        llvmcall($instrs, Vec{$W,$T1}, Tuple{Vec{$W,$T2}}, v)
    end
end
@generated function vconvert(::Type{_Vec{_W,T1}}, v::_Vec{_W,T2}) where {_W,T1 <: FloatingTypes, T2 <: Unsigned}
    W = _W + 1
    typ1 = llvmtype(T1)
    typ2 = llvmtype(T2)
    vtyp1 = "<$W x $typ1>"
    vtyp2 = "<$W x $typ2>"
    instrs = """
    %res = uitofp $vtyp2 %0 to $vtyp1
    ret $vtyp1 %res
    """
    quote
        $(Expr(:meta, :inline))
        llvmcall($instrs, Vec{$W,$T1}, Tuple{Vec{$W,$T2}}, v)
    end
end
@generated function vconvert(::Type{_Vec{_W,T1}}, v::_Vec{_W,T2}) where {_W,T1 <: Signed, T2 <: FloatingTypes}
    W = _W + 1
    typ1 = llvmtype(T1)
    typ2 = llvmtype(T2)
    vtyp1 = "<$W x $typ1>"
    vtyp2 = "<$W x $typ2>"
    instrs = """
    %res = fptosi $vtyp2 %0 to $vtyp1
    ret $vtyp1 %res
    """
    quote
        $(Expr(:meta, :inline))
        llvmcall($instrs, Vec{$W,$T1}, Tuple{Vec{$W,$T2}}, v)
    end
end
@generated function vconvert(::Type{_Vec{_W,T1}}, v::_Vec{_W,T2}) where {_W,T1 <: Unsigned, T2 <: FloatingTypes}
    W = _W + 1
    typ1 = llvmtype(T1)
    typ2 = llvmtype(T2)
    vtyp1 = "<$W x $typ1>"
    vtyp2 = "<$W x $typ2>"
    instrs = """
    %res = fptoui $vtyp2 %0 to $vtyp1
    ret $vtyp1 %res
    """
    quote
        $(Expr(:meta, :inline))
        llvmcall($instrs, Vec{$W,$T1}, Tuple{Vec{$W,$T2}}, v)
    end
end
@generated function vconvert(::Type{_Vec{_W,I1}}, v::_Vec{_W,I2}) where {_W,I1<:Integer,I2<:Integer}
    W = _W + 1
    typ1 = llvmtype(I1)
    typ2 = llvmtype(I2)
    b1 = 8sizeof(I1)
    b2 = 8sizeof(I2)
    op = if b1 < b2
        "trunc"
    elseif b1 > b2
        "zext"
    else
        return I1 === I2 ? :v : Expr(:block, Expr(:meta,:inline), Expr(:call, :vreinterpret, Vec{W,I1}, :v))
    end
    instrs = """
    %res = $op <$W x $typ2> %0 to <$W x $typ1>
    ret <$W x $typ1> %res
    """
    quote
        $(Expr(:meta, :inline))
        llvmcall($instrs, Vec{$W,$I1}, Tuple{Vec{$W,$I2}}, v)
    end
end
@generated function vconvert(::Type{_Vec{_W,T1}}, v::_Vec{_W,T2}) where {_W, T1 <: FloatingTypes, T2 <: FloatingTypes}
    W = _W + 1
    typ1 = llvmtype(T1)
    typ2 = llvmtype(T2)
    b1 = 8sizeof(T1)
    b2 = 8sizeof(T2)
    op = if b1 < b2
        "fptrunc"
    elseif b1 > b2
        "fpext"
    end
    instrs = """
    %res = $op <$W x $typ2> %0 to <$W x $typ1>
    ret <$W x $typ1> %res
    """
    quote
        $(Expr(:meta, :inline))
        llvmcall($instrs, Vec{$W,$T1}, Tuple{Vec{$W,$T2}}, v)
    end
end
@inline vconvert(::Type{_Vec{_W,T}}, v::_Vec{_W,T}) where {_W,T<:Integer} = v # specific definition
@inline vconvert(::Type{_Vec{_W,T}}, v::SVec) where {_W,T} = vconvert(_Vec{_W,T}, extract_data(v))
@inline vconvert(::Type{SVec{W,T}}, v) where {W,T} = SVec(vconvert(Vec{W,T}, extract_data(v)))
@inline vconvert(::Type{T}, v::T) where {T} = v
@inline vconvert(::Type{_Vec{_W,T}}, v::_Vec{_W,T}) where {_W,T<:FloatingTypes} = v
@inline vconvert(::Type{_Vec{_W,T}}, v::_Vec{_W,T}) where {_W,T} = v
@inline vconvert(::Type{SVec{W,T}}, v::SVec{W,T}) where {W,T} = v
@inline Base.convert(::Type{SVec{W,T}}, v) where {W,T} = SVec(vconvert(Vec{W,T}, extract_data(v)))
@inline VectorizationBase.SVec{W,T1}(v::SVec{W,T2}) where {W,T1<:FloatingTypes,T2<:IntegerTypes} = vconvert(SVec{W,T1}, v)
@inline VectorizationBase.SVec{W,T1}(v::SVec{W,T2}) where {W,T1<:IntegerTypes,T2<:FloatingTypes} = vconvert(SVec{W,T1}, v)
@inline vconvert(::Type{_Vec{_W,T1}}, s::T2) where {_W, T1, T2 <: FloatingTypes} = vbroadcast(_Vec{_W,T1}, convert(T1, s))
@inline vconvert(::Type{_Vec{_W,T1}}, s::T2) where {_W, T1 <: FloatingTypes, T2 <: Integer} = vbroadcast(_Vec{_W,T1}, convert(T1, s))
@inline vconvert(::Type{_Vec{_W,T1}}, s::T2) where {_W, T1 <: Integer, T2 <: Integer} = vbroadcast(_Vec{_W,T1}, s % T1)
@inline vconvert(::Type{_Vec{_W,T1}}, s::T1) where {_W, T1 <: Integer} = vbroadcast(_Vec{_W,T1}, s)
@inline vconvert(::Type{T}, s::Number) where {T <: Number} = convert(T, s)
@inline vconvert(::Type{T}, s::Integer) where {T <: Integer} = s % T
@inline vconvert(::Type{<:Vec{W}}, m::Mask{W}) where {W} = m
@inline vconvert(::Type{<:Vec{W}}, u::Unsigned) where {W} = Mask{W}(u)
@inline vconvert(::Type{SVec{W,T}}, m::Mask{W}) where {W,T} = m
@inline vconvert(::Type{SVec{W,T1}}, s::T1) where {W, T1 <: Unsigned} = vbroadcast(SVec{W,T1}, s)
@inline vconvert(::Type{SVec{W,T}}, u::Unsigned) where {W,T} = Mask{W}(u)

@inline promote_vtype(::Type{T}, ::Type{T}) where {T} = T
@inline promote_vtype(::Type{Mask{W,U}}, ::Type{V}) where {W, U, V <: AbstractSIMDVector{W}} = V
@inline promote_vtype(::Type{V}, ::Type{Mask{W,U}}) where {W, U, V <: AbstractSIMDVector{W}} = V
@inline promote_vtype(::Type{Mask{W,U}}, ::Type{SVec{W,T}}) where {W, U, T} = SVec{W,T}
@inline promote_vtype(::Type{SVec{W,T}}, ::Type{Mask{W,U}}) where {W, U, T} = SVec{W,T}
@inline promote_vtype(::Type{Mask{W,U}}, ::Type{Mask{W,U}}) where {W, U} = Mask{W,U}
@inline promote_vtype(::Type{Mask{W,U}}, ::Type{T}) where {W, U, T <: Number} = SVec{W,T}
@inline promote_vtype(::Type{T}, ::Type{Mask{W,U}}) where {W, U, T <: Number} = SVec{W,T}
@inline promote_vtype(::Type{T1}, ::Type{T2}) where {T1 <: Number, T2 <: Number} = promote_type(T1, T2)
@inline promote_vtype(::Type{Vec{W,T}}, ::Type{SVec{W,T}}) where {W,T} = SVec{W,T}

@inline Base.unsafe_trunc(::Type{I}, v::SVec{W,T}) where {W, I <: Signed, T <: FloatingTypes} = SVec(vconvert(Vec{W,I}, extract_data(v)))

@inline function promote_vtype(::Type{V1}, ::Type{V2}) where {W,T1,T2,V1<:AbstractStructVec{W,T1},V2<:AbstractStructVec{W,T2}}
    T = promote_type(T1, T2)
    SVec{W,T}
end
@inline function promote_vtype(::Type{_Vec{W,T1}}, ::Type{_Vec{W,T2}}) where {W,T1,T2}
    T = promote_type(T1, T2)
    _Vec{W,T}
end
@inline promote_type_vs(::Type{V}, ::Type{S}) where {V,S} = V # fallback
@inline promote_type_vs(::Type{V}, ::Type{S}) where {V<:Signed,S<:Signed} = V
@inline promote_type_vs(::Type{V}, ::Type{S}) where {V<:Unsigned,S<:Signed} = V
@inline promote_type_vs(::Type{V}, ::Type{S}) where {V<:Signed,S<:Unsigned} = unsigned(V)
@inline promote_type_vs(::Type{V}, ::Type{S}) where {V<:Unsigned,S<:Unsigned} = V
@inline promote_type_vs(::Type{V}, ::Type{S}) where {V<:FloatingTypes,S<:FloatingTypes} = V
@inline promote_type_vs(::Type{V}, ::Type{S}) where {V<:Unsigned,S<:FloatingTypes} = Base.uinttype(S)
@inline promote_type_vs(::Type{V}, ::Type{S}) where {V<:Signed,S<:FloatingTypes} = signed(Base.uinttype(S))
@inline promote_type_vs(::Type{<:Integer}, ::Type{Ptr{T}}) where {T} = Ptr{T}

@inline function promote_vtype(::Type{_Vec{W,T1}}, ::Type{T2}) where {W,T1,T2<:Number}
    T = promote_type_vs(T1,T2)
    _Vec{W,T}
end
@inline function promote_vtype(::Type{SVec{W,T1}}, ::Type{T2}) where {W,T1,T2<:Number}
    T = promote_type_vs(T1,T2)
    SVec{W,T}
end
@inline function promote_vtype(::Type{T2}, ::Type{_Vec{W,T1}}) where {W,T1,T2<:Number}
    T = promote_type_vs(T1,T2)
    _Vec{W,T}
end
@inline function promote_vtype(::Type{T2}, ::Type{SVec{W,T1}}) where {W,T1,T2<:Number}
    T = promote_type_vs(T1,T2)
    SVec{W,T}
end
@inline promote_vtype(::Type{<:_MM{W}}, ::Type{T}) where {W,T<:Number} = SVec{W,T}
@inline promote_vtype(::Type{T}, ::Type{<:_MM{W}}) where {W,T<:Number} = SVec{W,T}
@inline promote_vtype(::Type{<:_MM{W}}, ::Type{V}) where {W,V<:Vec{W}} = V
@inline promote_vtype(::Type{V}, ::Type{<:_MM{W}}) where {W,V<:Vec{W}} = V
@inline promote_vtype(::Type{<:_MM{W}}, ::Type{SVec{W,T}}) where {W,T<:Number} = SVec{W,T}
@inline promote_vtype(::Type{SVec{W,T}}, ::Type{<:_MM{W}}) where {W,T<:Number} = SVec{W,T}
@inline promote_vtype(::Type{T1}, ::Type{T2}, ::Type{T3}) where {T1,T2,T3} = promote_vtype(promote_vtype(T1, T2), T3)
@inline promote_vtype(::Type{T}, ::Type{T}, ::Type{T}) where {T} = T
@inline promote_vtype(::Type{_Vec{W,T}}, ::Type{_Vec{W,T}}) where {W,T} = _Vec{W,T}
@inline promote_vtype(::Type{SVec{W,T}}, ::Type{SVec{W,T}}) where {W,T} = SVec{W,T}

@inline function vpromote(args::Vararg{<:Any,N}) where {N}
    V = promote_vtype(typeof.(args)...)
    vconvert.(V, args)
end

@inline Base.float(v::AbstractStructVec{W,<:FloatingTypes}) where {W} = v
@inline Base.float(v::AbstractStructVec{W,I}) where {W, I <: Union{Int64,UInt64}} = vconvert(SVec{W,Float64}, v)
@inline Base.float(v::AbstractStructVec{W,I}) where {W, I <: Union{Int32,UInt32}} = vconvert(SVec{W,Float32}, v)

@inline Base.unsigned(v::AbstractStructVec{W,Int8}) where {W} = vconvert(SVec{W,UInt8}, v)
@inline Base.unsigned(v::AbstractStructVec{W,Int16}) where {W} = vconvert(SVec{W,UInt16}, v)
@inline Base.unsigned(v::AbstractStructVec{W,Int32}) where {W} = vconvert(SVec{W,UInt32}, v)
@inline Base.unsigned(v::AbstractStructVec{W,Int64}) where {W} = vconvert(SVec{W,UInt64}, v)
@inline Base.unsigned(v::AbstractStructVec{W,Int128}) where {W} = vconvert(SVec{W,UInt128}, v)

@generated function zeropad(v::_Vec{_W,T}) where {_W,T}
    W = _W + 1
    typ = llvmtype(T)
    W2 = W << 1
    sv = join((i for i ∈ 0:W2-1), ", i32 ")
    instr = """
    %res = shufflevector <$W x $typ> %0, <$W x $typ> zeroinitializer, <$W2 x i32> <i32 $(sv)>
    ret <$W2 x $typ> %res
    """
    :(llvmcall($instr, Vec{$W2,$T}, Tuple{Vec{$W,$T}}, v))
end

@inline Base.trunc(::Type{I}, v::SVec{W,T}) where {W,I<:Integer,T<:FloatingTypes} = vconvert(SVec{W,I}, v)

# Floating point formats

int_type(::Type{Float16}) = Int16
int_type(::Type{Float32}) = Int32
int_type(::Type{Float64}) = Int64
# int_type(::Type{Float128}) = Int128
# int_type(::Type{Float256}) = Int256

uint_type(::Type{Float16}) = UInt16
uint_type(::Type{Float32}) = UInt32
uint_type(::Type{Float64}) = UInt64
# uint_type(::Type{Float128}) = UInt128
# uint_type(::Type{Float256}) = UInt256

significand_bits(::Type{Float16}) = 10
significand_bits(::Type{Float32}) = 23
significand_bits(::Type{Float64}) = 52
# significand_bits(::Type{Float128}) = 112
# significand_bits(::Type{Float256}) = 136

exponent_bits(::Type{T}) where {T<:FloatingTypes} =
    8*sizeof(T) - 1 - significand_bits(T)
sign_bits(::Type{T}) where {T<:FloatingTypes} = 1

significand_mask(::Type{T}) where {T<:FloatingTypes} =
    uint_type(T)(uint_type(T)(1) << significand_bits(T) - 1)
exponent_mask(::Type{T}) where {T<:FloatingTypes} =
    uint_type(T)(uint_type(T)(1) << exponent_bits(T) - 1) << significand_bits(T)
sign_mask(::Type{T}) where {T<:FloatingTypes} =
    uint_type(T)(1) << (significand_bits(T) + exponent_bits(T))

# @generated function vecbool_to_unsigned(vb::Vec{N,Bool}) where {N}
#     Nout = VectorizationBase.nextpow2(max(N, 8))
#     Utype = VectorizationBase.mask_type(Nout)
#     btype = llvmtype(Bool)
#     vbtyp = "<$N x $btype>"
#     decls = String[]
#     instrs = String["%maskb = trunc $vbtyp %0 to <$N x i1>"]
#     u1name = Nout == N ? "%uout" : "%uabridged"
#     push!(instrs, "$u1name = bitcast <$N x i1> %maskb to i$N")
#     if Nout > N
#         push!(instrs, "%uout = zext i$N %uabridged to i$Nout")
#     end
#     push!(instrs, "ret i$Nout %uout")
#     quote
#         $(Expr(:meta,:inline))
#         llvmcall($(join(instrs, "\n")), $Utype, Tuple{Vec{$N,Bool}}, vb)
#     end
# end


text:
using SolidState
using LinearAlgebra
using StaticArrays
using Test

# types.jl

@test vec3d(1, 0, 0) ⋅ vec3d(0, 1, 0) ≈ 0
@test vec3d(2, 1, 3) ⋅ vec3d(2, 1, 3) ≈ 14
@test vec3d(1, 0, 0) × vec3d(0, 1, 0) ≈ vec3d(0, 0, 1)

@test norm2(vec3d(1, 0, 0)) ≈ 1
@test norm2(vec3d(1, 2, 3)) ≈ 14

# transformations.jl

@test rotate(π/6, vec3d(1, 0, 0)) ≈ rotatex(π/6)
@test rotate(π/6, vec3d(0, 1, 0)) ≈ rotatey(π/6)
@test rotate(π/6, vec3d(0, 0, 1)) ≈ rotatez(π/6)

@test translate(0.1, 0.2, 0.3) ≈ translate(vec3d(0.1, 0.2, 0.3))
@test translate(1, 0, 0) * translate(0, 2, 0) * translate(0, 0, 3) ≈ translate(1, 2, 3)

@test scalex(2) * scaley(3) * scalez(4) ≈ scale(2, 3, 4)

# Vectors should not vary when translated
@test translate(4, 5, 6) * vec3d(1, 2, 3) ≈ vec3d(1, 2, 3)
@test rotatez(π/2) * vec3d(1, 0, 0) ≈ vec3d(0, 1, 0)

# Remember that a ray transforms under the *inverse* transformation, so
# vector [1, 0, 0] rotates by -π/2 around the y axis
ray = rotatey(π/2) * Ray(vec3d(0, 0, 0), vec3d(1, 0, 0))
@test ray.origin ≈ vec3d(0, 0, 0)
@test ray_at_dist(ray, 10) ≈ vec3d(0, 0, 10)
@test ray(10) ≈ vec3d(0, 0, 10)

# cameras.jl

(u, v) = xy_to_uv(127, 348, 1024, 768)
(x, y) = uv_to_xy(u, v, 1024, 768)
@test [x, y] ≈ [127.0, 348.0]

camera = PerspectiveCamera(vec3d(0, 0, -2), # Position
                           vec3d(0, 0, 0),  # Look at
                           vec3d(0, 1, 0),  # Up
                           vec3d(1, 0, 0))  # Right

x, y, visible = project_point(camera, vec3d(0, 0, 0), 1024, 768)
@test x ≈ 1024 / 2
@test y ≈ 768 / 2
@test visible

x, y, visible = project_point(camera, vec3d(0, 0, -2), 1024, 768)
@test !visible

camera = OrthoCamera(vec3d(0, 0, -2), # Position
                     vec3d(0, 1, 0),  # Up
                     vec3d(1, 0, 0))  # Right

x, y, visible = project_point(camera, vec3d(0, 0, 0), 1024, 768)
@test x ≈ 1024 / 2
@test y ≈ 768 / 2
@test visible


text:
# Gun: variant P (polynomial case; only repeated nodes)

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test

include(joinpath("..", "rk_helper", "gun_test_utils.jl"))

@bench @testset "NLEIGS: Gun variant P" begin
    verbose = displaylevel

    nep, Σ, _, v, nodes, funres = gun_init()

    # solve nlep
    lambda, X, res, solution_info = nleigs(nep, Σ, logger=verbose > 0 ? 1 : 0, maxit=100, v=v, leja=0, nodes=nodes, reusefact=2, errmeasure=funres, return_details=verbose > 1)

    verify_lambdas(17, nep, lambda, X)

    if verbose > 1
        include("nleigs_residual_plot.jl")
        nleigs_residual_plot("Gun: variant P", solution_info, Σ; ylims=[1e-17, 1e-1])
    end
end


text:
module ProfileViewGtk

using InteractiveUtils
using Gtk.ShortNames, GtkReactive, Colors, FileIO, IntervalSets
import Cairo
using Graphics

using Gtk.GConstants.GdkModifierType: SHIFT, CONTROL, MOD1

mutable struct ZoomCanvas
    bb::BoundingBox  # in user-coordinates
    c::Canvas
end

function __init__()
    @eval import ProfileView
end

function closeall()
    for (w, _) in window_wrefs
        destroy(w)
    end
    nothing
end

const window_wrefs = WeakKeyDict{Gtk.GtkWindowLeaf,Nothing}()

function view(data = Profile.fetch(); lidict=nothing, kwargs...)
    bt, uip, counts, lidict, lkup = ProfileView.prepare_data(data, lidict)
    # Display in a window
    c = canvas(UserUnit)
    set_gtk_property!(widget(c), :expand, true)
    f = Frame(c)
    tb = Toolbar()
    bx = Box(:v)
    push!(bx, tb)
    push!(bx, f)
    tb_open = ToolButton("gtk-open")
    tb_save_as = ToolButton("gtk-save-as")
    push!(tb, tb_open)
    push!(tb, tb_save_as)
    signal_connect(open_cb, tb_open, "clicked", Nothing, (), false, (widget(c),kwargs))
    signal_connect(save_as_cb, tb_save_as, "clicked", Nothing, (), false, (widget(c),data,lidict,kwargs))
    win = Window(bx, "Profile")
    GtkReactive.gc_preserve(win, c)
    # Register the window with closeall
    window_wrefs[win] = nothing
    signal_connect(win, :destroy) do w
        delete!(window_wrefs, win)
    end

    if data != nothing && !isempty(data)
        viewprof(c, bt, uip, counts, lidict, lkup; kwargs...)
    end

    # Ctrl-w and Ctrl-q destroy the window
    signal_connect(win, "key-press-event") do w, evt
        if evt.state == CONTROL && (evt.keyval == UInt('q') || evt.keyval == UInt('w'))
            @async destroy(w)
            nothing
        end
    end

    Gtk.showall(win)
end

function viewprof(c, bt, uip, counts, lidict, lkup; C = false, colorgc = true, fontsize = 12, combine = true, pruned=[])
    img, lidict, imgtags = ProfileView.prepare_image(bt, uip, counts, lidict, lkup, C, colorgc, combine, pruned)
    img24 = Matrix(UInt32[reinterpret(UInt32, convert(RGB24, img[i,j])) for i = 1:size(img,1), j = size(img,2):-1:1]')
    fv = XY(0.0..size(img24,2), 0.0..size(img24,1))
    zr = Signal(ZoomRegion(fv, fv))
    sigrb = init_zoom_rubberband(c, zr)
    sigpd = init_pan_drag(c, zr)
    sigzs = init_zoom_scroll(c, zr)
    sigps = init_pan_scroll(c, zr)
    surf = Cairo.CairoRGBSurface(img24)
    sigredraw = draw(c, zr) do widget, r
        ctx = getgc(widget)
        set_coordinates(ctx, r)
        rectangle(ctx, BoundingBox(r.currentview))
        set_source(ctx, surf)
        p = Cairo.get_source(ctx)
        Cairo.pattern_set_filter(p, Cairo.FILTER_NEAREST)
        fill(ctx)
    end
    lasttextbb = BoundingBox(1,0,1,0)
    sigmotion = map(c.mouse.motion) do btn
        # Repair image from ovewritten text
        if c.widget.is_realized && c.widget.is_sized
            ctx = getgc(c)
            if Graphics.width(lasttextbb) > 0
                r = value(zr)
                set_coordinates(ctx, r)
                rectangle(ctx, lasttextbb)
                set_source(ctx, surf)
                p = Cairo.get_source(ctx)
                Cairo.pattern_set_filter(p, Cairo.FILTER_NEAREST)
                fill(ctx)
            end
            # Write the info
            xu, yu = btn.position.x, btn.position.y
            tag = gettag(xu, yu)
            if tag != ProfileView.TAGNONE
                li = lidict[tag.ip]
                str = ""
                for l in li
                    str = string(str, string(basename(string(l.file)), ", ", l.func, ": line ", l.line), "; ")
                end
                set_source(ctx, ProfileView.fontcolor)
                Cairo.set_font_face(ctx, "sans-serif $(fontsize)px")
                xi = value(zr).currentview.x
                xmin, xmax = minimum(xi), maximum(xi)
                lasttextbb = deform(Cairo.text(ctx, xu, yu, str, halign = xu < (2xmin+xmax)/3 ? "left" : xu < (xmin+2xmax)/3 ? "center" : "right"), -2, 2, -2, 2)
            end
            reveal(c)
        end
    end
    # From a given position, find the underlying tag
    function gettag(xu, yu)
        x = ceil(Int, Float64(xu))
        y = ceil(Int, Float64(yu))
        Y = size(imgtags, 2)
        x = max(1, min(x, size(imgtags, 1)))
        y = max(1, min(y, Y))
        imgtags[x,Y-y+1]
    end
    # Hover over a block and see the source line
    # Left-click prints the full path, function, and line to the console
    # Right-click calls the edit() function
    sigshow = map(c.mouse.buttonpress) do btn
        if btn.button == 1 || btn.button == 3
            ctx = getgc(c)
            xu, yu = btn.position.x, btn.position.y
            tag = gettag(xu, yu)
            if tag != ProfileView.TAGNONE
                li = lidict[tag.ip]
                if btn.button == 1
                    firstline = true
                    for l in li
                        if !firstline
                            print("  ")
                        else
                            firstline = false
                        end
                        println(l.file, ", ", l.func, ": line ", l.line)
                    end
                elseif btn.button == 3
                    if !isempty(li)
                        l = first(li)
                        edit(string(l.file),l.line)
                    end
                end
            end
        end
    end
    append!(c.preserved, [sigrb, sigpd, sigzs, sigps, sigredraw, sigmotion, sigshow])
    nothing
end

function open_cb(::Ptr, settings::Tuple)
    c, kwargs = settings
    selection = open_dialog("Load profile data", toplevel(c), ("*.jlprof","*"))
    isempty(selection) && return nothing
    data, lidict = load(File(format"JLD", selection), "li", "lidict")
    bt, uip, counts, lidict, lkup = ProfileView.prepare_data(data, lidict)
    viewprof(c, bt, uip, counts, lidict, lkup; kwargs...)
    nothing
end

function save_as_cb(::Ptr, profdata::Tuple)
    c, data, lidict = profdata
    selection = save_dialog("Save profile data as JLD file", toplevel(c), ("*.jlprof",))
    isempty(selection) && return nothing
    FileIO.save(File(format"JLD", selection), "li", data, "lidict", lidict)
    nothing
end

end


text:
@testset "Running WBS tests..." begin

    n = 1000;        # Number of samples
    λ = 100;         # Frequency of changes

    ########################
    # Normal mean segments #
    ########################
    μ, σ = Normal(0.0, 10.0), 1.0
    data, cps = @changepoint_sampler n λ Normal(μ, σ)
    seg_cost_sSIC = Changepoints.sSIC(data);

    @test_nowarn WBS(data);
    @test_nowarn get_WBS_changepoints(seg_cost_sSIC, WBS(data) );
    @test_nowarn WBS(data; do_seeded = true);

end


text:
module FFTViews

using Base: tail, unsafe_length, @propagate_inbounds
using FFTW

# A custom rangetype that will be used for indices and never throws a
# boundserror because the domain is actually periodic.
using CustomUnitRanges
include(CustomUnitRanges.filename_for_urange)

@static if isdefined(Base, :IdentityUnitRange)
    const indextypes = (URange, Base.IdentityUnitRange{<:URange})
    const FFTVRange{T} = Union{URange{T}, Base.Slice{URange{T}}, Base.IdentityUnitRange{URange{T}}}
    indrange(i) = Base.IdentityUnitRange(URange(first(i)-1, last(i)-1))
else
    const indextypes = (URange, Base.Slice{<:URange})
    const FFTVRange{T} = Union{URange{T}, Base.Slice{URange{T}}}
    indrange(i) = Base.Slice(URange(first(i)-1, last(i)-1))
end

for T in indextypes
    @eval begin
        Base.checkindex(::Type{Bool}, ::$T, ::Base.Slice) = true
        Base.checkindex(::Type{Bool}, ::$T, ::Base.LogicalIndex) = true
        Base.checkindex(::Type{Bool}, ::$T, ::Real) = true
        Base.checkindex(::Type{Bool}, ::$T, ::AbstractRange) = true
        Base.checkindex(::Type{Bool}, ::$T, ::AbstractVector{Bool}) = true
        Base.checkindex(::Type{Bool}, ::$T, ::AbstractArray{Bool}) = true
        Base.checkindex(::Type{Bool}, ::$T, ::AbstractArray) = true
    end
    if isdefined(Base, :IdentityUnitRange)
        @eval Base.checkindex(::Type{Bool}, ::$T, ::Base.IdentityUnitRange) = true
    end
end

export FFTView

abstract type AbstractFFTView{T,N} <: AbstractArray{T,N} end

struct FFTView{T,N,A<:AbstractArray} <: AbstractFFTView{T,N}
    parent::A

    function FFTView{T,N,A}(parent::A) where {T,N,A}
        new{T,N,A}(parent)
    end
end
FFTView(parent::AbstractArray{T,N}) where {T,N} = FFTView{T,N,typeof(parent)}(parent)
FFTView{T,N}(dims::Dims{N}) where {T,N} = FFTView(Array{T,N}(undef, dims))
FFTView{T}(dims::Dims{N}) where {T,N} = FFTView(Array{T,N}(undef, dims))

# Note: there are no bounds checks because it's all periodic
@inline @propagate_inbounds function Base.getindex(F::FFTView{T,N}, I::Vararg{Int,N}) where {T,N}
    P = parent(F)
    @inbounds ret = P[reindex(FFTView, axes(P), I)...]
    ret
end

@inline @propagate_inbounds function Base.setindex!(F::FFTView{T,N}, val, I::Vararg{Int,N}) where {T,N}
    P = parent(F)
    @inbounds P[reindex(FFTView, axes(P), I)...] = val
end

Base.parent(F::AbstractFFTView) = F.parent
Base.axes(F::AbstractFFTView) = map(indrange, axes(parent(F)))
Base.size(F::AbstractFFTView) = size(parent(F))

function Base.similar(A::AbstractArray, T::Type, shape::Tuple{FFTVRange,Vararg{FFTVRange}})
    all(x->first(x)==0, shape) || throw(BoundsError("cannot allocate FFTView with the first element of the range non-zero"))
    FFTView(similar(A, T, map(length, shape)))
end

function Base.similar(f::Union{Function,Type}, shape::Tuple{FFTVRange,Vararg{FFTVRange}})
    all(x->first(x)==0, shape) || throw(BoundsError("cannot allocate FFTView with the first element of the range non-zero"))
    FFTView(similar(f, map(length, shape)))
end

Base.reshape(F::FFTView{_,N}, ::Type{Val{N}}) where {_,N}   = F
Base.reshape(F::FFTView{_,M}, ::Type{Val{N}}) where {_,M,N} = FFTView(reshape(parent(F), Val(N)))

FFTW.fft(F::FFTView; kwargs...) = fft(parent(F); kwargs...)
FFTW.rfft(F::FFTView; kwargs...) = rfft(parent(F); kwargs...)
FFTW.fft(F::FFTView, dims; kwargs...) = fft(parent(F), dims; kwargs...)
FFTW.rfft(F::FFTView, dims; kwargs...) = rfft(parent(F), dims; kwargs...)

@inline reindex(::Type{V}, inds, I) where {V} = (_reindex(V, inds[1], I[1]), reindex(V, tail(inds), tail(I))...)
reindex(::Type{V}, ::Tuple{}, ::Tuple{}) where {V} = ()
_reindex(::Type{FFTView}, ind, i) = modrange(i+1, ind)

if VERSION >= v"1.7.0-beta4"
    # https://github.com/JuliaLang/julia/pull/40382
    modrange(i, rng::AbstractUnitRange) = mod(i-first(rng), length(rng))+first(rng)
else
    modrange(i, rng::AbstractUnitRange) = mod(i-first(rng), unsafe_length(rng))+first(rng)
end

end # module


text:
# test_with_actors2.jl

# Basically the same as with test_local_channels.jl

using Pkg; Pkg.activate(".")

@info("remove old stuff")
cmd = `rm test_invoicing.txt test_invoicing_paid.txt test_ledger.txt test_journal.txt invoicenbr.txt`
run(cmd)

using Rocket, Query, DataFrames, Query

# enable distrbuted computing
using Distributed
@info("Enable distributed computing")

# this should be the next step
np = addprocs(4; exeflags=`--project=$(Base.active_project())`)
@info("number of processes is $(length(np))")

# activate the packages (after the processes are created)
@everywhere begin
    using AppliSales
    using AppliGeneralLedger
    using AppliAR
end;

@info("Distributed computing enabled")

struct StmActor <: Actor{String} end
Rocket.on_next!(actor::StmActor, data::String) = begin
    if data == "READ_STMS"
        stms = AppliAR.read_bank_statements("./bank.csv")
        unpaid_inv = retrieve_unpaid_invoices()
        entries = @fetch AppliAR.process(unpaid_inv, stms)
        #println(entries)
        subscribe!(from(entries), gl_actor)
    end
end
Rocket.on_complete!(actor::StmActor) = println("StmActor completed!")
Rocket.on_error!(actor::StmActor, err) = error(err)

struct SalesActor <: Actor{String} end
Rocket.on_next!(actor::SalesActor, data::String) = begin
    if data == "START"
        #ar_actor = ARActor()
        orders = @fetch AppliSales.process()
        subscribe!(from(orders), ar_actor)
    end
end
Rocket.on_complete!(actor::SalesActor) = println("SalesActor completed!")
Rocket.on_error!(actor::SalesActor, err) = error(err)

struct ARActor <: Actor{AppliSales.Order}
    values::Vector{AppliGeneralLedger.JournalEntry}
    ARActor() = new(Vector{AppliGeneralLedger.JournalEntry}())
end
Rocket.on_next!(actor::ARActor, data::AppliSales.Order) = begin
        d = @fetch AppliAR.process([data])
        push!(actor.values, d[1])
end
Rocket.on_complete!(actor::ARActor) = begin
    #println(actor.values)
    println(typeof(actor.values))
    subscribe!(from(actor.values), gl_actor)
    println("ARActor Completed!")
end
Rocket.on_error!(actor::ARActor, err) = error(err)

struct GLActor <: Actor{Any} end
Rocket.on_next!(actor::GLActor, data::Any) = begin
    if data isa AppliGeneralLedger.JournalEntry
        result = @fetch AppliGeneralLedger.process([data])
    end
end
Rocket.on_complete!(actor::GLActor) = println("GLActor completed!")
Rocket.on_error!(actor::GLActor, err) = error(err)

sales_actor = SalesActor()
ar_actor = ARActor()
gl_actor = GLActor()
stm_actor = StmActor()

subscribe!(from(["START"]), sales_actor)

subscribe!(from(["READ_STMS"]), stm_actor)

r2 = AppliGeneralLedger.read_from_file("./test_ledger.txt")
df = DataFrame(r2)
println("\nGeneral Ledger mutations\n========================")
show(df)

r = DataFrame(report())
println("\nUnpaid invoices\n$("="^15)")
println(r)

df2 = r2 |> @filter(_.accountid == 1300) |> DataFrame
balance_1300 = sum(df2.debit - df2.credit)

df2 = df |> @filter(_.accountid == 8000) |> DataFrame
balance_8000 = sum(df2.credit - df2.debit)

df2 = df |> @filter(_.accountid == 1150) |> DataFrame
balance_1150 = sum(df2.debit - df2.credit)

df2 = df |> @filter(_.accountid == 4000) |> DataFrame
balance_4000 = sum(df2.credit - df2.debit)

println("")
println("Balance Accounts Receivable is $balance_1300. Should be 1210..")
println("Sales is $balance_8000. Should be 4000.0.")
println("Balance bank is $balance_1150. Should be 3630.0.")
println("Balance VAT is $balance_4000. Shouldbe 840.0.")

# end


text:
"""
Validation Eigenvalues with Andes
This case study defines an 11-bus that uses GENROU, SEXS and TGOV1.
We compare the obtained eigenvalues with Python based tool ANDES.
Results are slighty different due to the difference on shaft and
turbine governor power-torque approximations. Removing the division
for ω for output torque in TGOV1 model produces exact eigenvalue results.
"""

##################################################
############### SOLVE PROBLEM ####################
##################################################

raw_file = joinpath(TEST_FILES_DIR, "benchmarks/andes/test36/11BUS_KUNDUR.raw")
dyr_file = joinpath(TEST_FILES_DIR, "benchmarks/andes/test36/11BUS_KUNDUR_TGOV.dyr")
eigs_andes_csv = joinpath(TEST_FILES_DIR, "benchmarks/andes/test36/eigs_tgov_andes.csv")

@testset "Test 36 Eigenvalues" begin
    path = (joinpath(pwd(), "test-eigs"))
    !isdir(path) && mkdir(path)
    try
        sys = System(raw_file, dyr_file)

        # Define Simulation Problem
        sim = Simulation!(
            ResidualModel,
            sys, #system
            path,
            (0.0, 20.0), #time span
        )
        ss = small_signal_analysis(sim)
        eigs = ss.eigenvalues

        A = strip.(readdlm(eigs_andes_csv, ','))
        eigs_andes = parse.(Complex{Float64}, A)
        @test LinearAlgebra.norm(eigs - eigs_andes, 2) / length(eigs) < 1.0
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end


text:
using DelimitedFiles
function id(str)
	str = replace(str, "F" => "0")
	str = replace(str, "B" => "1")
	str = replace(str, "L" => "0")
	str = replace(str, "R" => "1")
	parse(Int, str, base=2)
end

# prints solution of both parts
function x()
	arr = sort(id.(vec(readdlm("5.txt"))))
	println(maximum(arr))
	for i = 1:length(arr)-1
		if arr[i] + 1 != arr[i+1]
			println(arr[i]+1)
		end
	end
end


text:
module StarDust

using ..Ahorn, Maple

placements = Maple.StarDust

function Ahorn.canFgBg(effect::Maple.StarDust)
    return true, true
end

end


text:
using Plots
pyplot()

A = [1.0 2.0
     2.0 1.0]

xmin, xmax = -3, 3
ymin, ymax = -3, 3

evals, evecs = eig(A)

# Extract x and y values of each eigenvector, arrange for plotting
x_vals = [zeros(1, size(A, 1)); evecs[1,:]]
y_vals = [zeros(1, size(A, 1)); evecs[2,:]]

# Image of each eigenvector
evecs =(evecs[:,1], evecs[:,2])
images = [[] []]'
for v in evecs
  v = A * v
  images = hcat(images, v)
end

# Extract x and y values of each image, arrange for plotting
x_ims = [zeros(1, size(A, 1)); images[1, :]]
y_ims = [zeros(1, size(A, 1)); images[2, :]]

x = linspace(xmin, xmax, 3)
as = []
for v in evecs
    a = v[2] / v[1]
    push!(as, a .* x)
end

# Plot each vector and its image
plot(x_ims, y_ims, arrow=true, color=:red,
     legend=:none, linewidth=1.8)
plot!(x_vals, y_vals, arrow=:true, color=:blue,
      xlims=(xmin, xmax), ylims=(ymin, ymax),
      legend=:none, linewidth=1.8)
plot!(x, as, color=:black, linewidth=0.3)
vline!([0], color=:black)
hline!([0], color=:black)
plot!(foreground_color_axis=:white, foreground_color_text=:white,
      foreground_color_border=:white)

max_stars_repo_name: mwhchen/quantecon

max_issues_repo_name: mwhchen/quantecon

max_forks_repo_name: mwhchen/quantecon


text:
using Plots

@fastmath function SIR_fs(N,bet,gamm)
    final_size = zeros(N+1,1);
    final_size[2] = 1;
    for Z2 in 0:N
        @inbounds for Z1 in Z2+1:N-1
            p1 = 1 / ( 1 + gamm/(bet*(N-Z1)));
            final_size[Z1+2] = final_size[Z1+2] + final_size[Z1+1]*p1;
            final_size[Z1+1] = final_size[Z1+1]*(1-p1);
        end
    end
    return final_size;
end

N = 20;
bet = 2/(N-1);
gamm = 1.0;

final_size = SIR_fs(N,bet,gamm);

@time final_size = SIR_fs(N,bet,gamm);

bar(0:N,final_size)

max_stars_repo_name: epimodels/epicookbook

max_issues_repo_name: epimodels/epicookbook

max_forks_repo_name: epimodels/epicookbook


text:
# Utilities for re-initializing initial conditions on the grid
"""
    CompleteAndReinit(complete_state, idxs, D)
Helper struct that completes a state and reinitializes the integrator once called
as a function with arguments `f(integ, y)` with `integ` the initialized dynamical
system integrator and `y` the projected initial condition on the grid.
"""
struct CompleteAndReinit{C, Y, R}
    complete_state::C
    u::Vector{Float64} # dummy variable for a state in full state space
    idxs::SVector{Y, Int}
    remidxs::R
end
function CompleteAndReinit(complete_state, idxs, D::Int)
    remidxs = setdiff(1:D, idxs)
    remidxs = isempty(remidxs) ? nothing : SVector(remidxs...)
    u = zeros(D)
    if complete_state isa AbstractVector
        @assert eltype(complete_state) <: Number
    end
    return CompleteAndReinit(complete_state, u, idxs, remidxs)
end
function (c::CompleteAndReinit{<: AbstractVector})(integ, y)
    c.u[c.idxs] .= y
    if !isnothing(c.remidxs)
        c.u[c.remidxs] .= c.complete_state
    end
    reinit!(integ, c.u)
end
function (c::CompleteAndReinit)(integ, y) # case where `complete_state` is a function
    c.u[c.idxs] .= y
    if !isnothing(c.remidxs)
        c.u[c.remidxs] .= c.complete_state(y)
    end
    reinit!(integ, c.u)
end

@generated function generate_ic_on_grid(grid::NTuple{B, T}, ind) where {B, T}
    gens = [:(grid[$k][ind[$k]]) for k=1:B]
    quote
        Base.@_inline_meta
        @inbounds return SVector{$B, Float64}($(gens...))
    end
end


text:
using Vinyl: @reset, shift
using Test

struct FakeArray end
Base.length(::FakeArray) = shift(k -> k)

@testset "Continuations" begin

k = @reset shift(k -> k)+5
@test k(1) == 6
@test k(2) == 7

@test @reset(2*(shift(k -> k(5))+1)) == 12

@test @reset(2*shift(k -> k(k(4)))) == 16

# TODO fix #1
@reset begin
  for i = 1:2
    x = shift(k -> (i,k()))
  end
  ()
end

f(x) = length(x) > 10 ? :big : :small
k = @reset f(FakeArray())

# see #2
@test k(5) == :small
@test k(15) == :big

end


text:
include("plot_recipes.jl")
include("interpolate.jl")
include("convert.jl")


text:
@testset "cpnn.jl" begin

    @testset "CPNN w/o global features" begin
        modelNN = SeaPearl.CPNN(
            graphChain = Flux.Chain(
                GeometricFlux.GraphConv(3 => 3),
                GeometricFlux.GraphConv(3 => 3),
            ),
            nodeChain = Flux.Chain(
                Flux.Dense(3, 3),
                Flux.Dense(3, 3),
            ),
            outputChain = Flux.Chain(
                Flux.Dense(3, 2)
            )
        )

        graphs = Matrix.(adjacency_matrix.([random_regular_graph(10, 4) for i = 1:4]))
        nodeFeatures = [rand(3, 10) for i = 1:4]
        featuredGraphs = [FeaturedGraph(g; nf=nf) for (g, nf) in zip(graphs, nodeFeatures)]

        trajectoryVector = SeaPearl.DefaultTrajectoryState.(featuredGraphs, rand(1:10, 4))
        nnInput = trajectoryVector |> cpu

        output = modelNN(nnInput)

        @test size(output) == (2, 4)

        grad = gradient(nnInput) do inp
            modelNN(inp)[2,2]
        end
        grad = grad[1]

        @test isnothing(grad.allValuesIdx)
        @test isnothing(grad.variableIdx)
        @test isnothing(grad.fg.ef)
        @test isnothing(grad.fg.gf)
        @test all(grad.fg.nf[:, :, [1, 3, 4]] .== 0)
        @test !all(grad.fg.nf[:, :, 2] .== 0)
        @test all(grad.fg.graph[:, :, [1, 3, 4]] .== 0)
        @test !all(grad.fg.graph[:, :, 2] .== 0)
    end

    @testset "CPNN w/ global features" begin
        modelNN = SeaPearl.CPNN(
            graphChain = Flux.Chain(
                GeometricFlux.GraphConv(3 => 3),
                GeometricFlux.GraphConv(3 => 3),
            ),
            nodeChain = Flux.Chain(
                Flux.Dense(3, 3),
                Flux.Dense(3, 3),
            ),
            globalChain = Flux.Chain(
                Flux.Dense(2, 3)
            ),
            outputChain = Flux.Chain(
                Flux.Dense(6, 2)
            ),
        )

        graphs = Matrix.(adjacency_matrix.([random_regular_graph(10, 4) for i = 1:4]))
        nodeFeatures = [rand(3, 10) for i = 1:4]
        globalFeatures = [rand(2) for i = 1:4]
        featuredGraphs = [FeaturedGraph(g; nf=nf, gf=gf) for (g, nf, gf) in zip(graphs, nodeFeatures, globalFeatures)]

        trajectoryVector = SeaPearl.DefaultTrajectoryState.(featuredGraphs, rand(1:10, 4))
        nnInput = trajectoryVector |> cpu

        output = modelNN(nnInput)

        @test size(output) == (2, 4)

        grad = gradient(nnInput) do inp
            modelNN(inp)[2,2]
        end
        grad = grad[1]

        @test isnothing(grad.allValuesIdx)
        @test isnothing(grad.variableIdx)
        @test isnothing(grad.fg.ef)
        @test all(grad.fg.gf[:, [1, 3, 4]] .== 0)
        @test !all(grad.fg.gf[:, 2] .== 0)
        @test all(grad.fg.nf[:, :, [1, 3, 4]] .== 0)
        @test !all(grad.fg.nf[:, :, 2] .== 0)
        @test all(grad.fg.graph[:, :, [1, 3, 4]] .== 0)
        @test !all(grad.fg.graph[:, :, 2] .== 0)
    end
end


text:
using Documenter, ArbFloats

makedocs(
    sitename = "ArbFloats.jl",
    modules  = [ArbFloats],
    format   = Documenter.Formats.HTML,
    clean    = false,
    pages    = Any["Home" => "index.md"],
    doctest  = true,
)

deploydocs(
    target = "build",
    deps   = Deps.pip("pygments", "mkdocs", "python-markdown-math"),
    make   = nothing,
    repo   = "github.com/JuliaArbTypes/ArbFloats.jl.git",
    julia  = "0.5.0-rc3",
    osname = "linux",
)


text:
module CornerSpaceRenorm

using QuantumOptics
include("types.jl");
export AbstractLattice, Lattice, SquareLattice, NdLattice, pbc_from_obc, obc_from_pbc
export eltype, vertices, edges, nv, ne, has_edge, has_vertex, inneighbors, outneighbors, is_directed
export AbstractSystem, SquareSystem, NdSystem
export gplot, plot_system
export CornerBasis
include("operators.jl");
export hamiltonian, dissipators
include("corner.jl")
export corner_subspace, cornerize, vmerge, hmerge, hermitianize, hermitianize!
include("steadystate.jl")
export steadystate_bicg, steadystate_bicg!
end # module


text:
####
#### Useful kernels
####

function velocity_div!(grid::RegularCartesianGrid, u, v, w, div)
    @loop for k in (1:grid.Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                @inbounds div[i, j, k] = div_f2c(grid, u, v, w, i, j, k)
            end
        end
    end
end

####
#### Horizontal averages
####

"""
    HorizontalAverage{F, R, P, I, Ω} <: AbstractDiagnostic

A diagnostic for computing horizontal average of a field.
"""
mutable struct HorizontalAverage{F, R, P, I, Ω, G} <: AbstractDiagnostic
          field :: F
         result :: P
      frequency :: Ω
       interval :: I
       previous :: Float64
    return_type :: R
           grid :: G
end

"""
    HorizontalAverage(model, field; frequency=nothing, interval=nothing, return_type=Array)

Construct a `HorizontalAverage` of `field`.

After the horizontal average is computed it will be stored in the `result` property.

The `HorizontalAverage` can be used as a callable object that computes and returns the
horizontal average.

A `frequency` or `interval` (or both) can be passed to indicate how often to run this
diagnostic if it is part of `model.diagnostics`. `frequency` is a number of iterations
while `interval` is a time interval in units of `model.clock.time`.

A `return_type` can be used to specify the type returned when the `HorizontalAverage` is
used as a callable object. The default `return_type=Array` is useful when running a GPU
model and you want to save the output to disk by passing it to an output writer.
"""
function HorizontalAverage(field; frequency=nothing, interval=nothing, return_type=Array)
    arch = architecture(field)
    result = zeros(arch, field.grid, 1, 1, field.grid.Tz)
    return HorizontalAverage(field, result, frequency, interval, 0.0, return_type, field.grid)
end

# Normalize a horizontal sum to get the horizontal average.
normalize_horizontal_sum!(havg, grid) = havg.result /= (grid.Nx * grid.Ny)

"""
    run_diagnostic(model, havg::HorizontalAverage{NTuple{1}})

Compute the horizontal average of `havg.field` and store the result in `havg.result`.
"""
function run_diagnostic(model, havg::HorizontalAverage)
    zero_halo_regions!(parent(havg.field), model.grid)
    sum!(havg.result, parent(havg.field))
    normalize_horizontal_sum!(havg, model.grid)
    return nothing
end

function (havg::HorizontalAverage{F, Nothing})(model) where F
    run_diagnostic(model, havg)
    return havg.result
end

function (havg::HorizontalAverage)(model)
    run_diagnostic(model, havg)
    return havg.return_type(havg.result)
end

####
#### NaN checker
####

"""
    NaNChecker{F} <: AbstractDiagnostic

A diagnostic that checks for `NaN` values and aborts the simulation if any are found.
"""
struct NaNChecker{F} <: AbstractDiagnostic
    frequency :: Int
       fields :: F
end

"""
    NaNChecker(model; frequency, fields)

Construct a `NaNChecker` for `model`. `fields` should be a `Dict{Symbol,Field}`. A
`frequency` should be passed to indicate how often to check for NaNs (in number of
iterations).
"""
function NaNChecker(model; frequency, fields)
    return NaNChecker(frequency, fields)
end

function run_diagnostic(model::Model, nc::NaNChecker)
    for (name, field) in nc.fields
        if any(isnan, field.data.parent)
            t, i = model.clock.time, model.clock.iteration
            error("time = $t, iteration = $i: NaN found in $name. Aborting simulation.")
        end
    end
end

####
#### Timeseries diagnostics
####

"""
    Timeseries{D, Ω, I, T, TT} <: AbstractDiagnostic

A diagnostic for collecting and storing timeseries.
"""
struct Timeseries{D, Ω, I, T, TT} <: AbstractDiagnostic
    diagnostic :: D
     frequency :: Ω
      interval :: I
          data :: T
          time :: Vector{TT}
end

"""
    Timeseries(diagnostic, model; frequency=nothing, interval=nothing)

A `Timeseries` `Diagnostic` that records a time series of `diagnostic(model)`.

Example
=======
```julia
julia> model = Model(grid=RegularCartesianGrid(size=(16, 16, 16), length=(1, 1, 1)));

julia> max_u = Timeseries(FieldMaximum(abs, model.velocities.u), model; frequency=1)

julia> model.diagnostics[:max_u] = max_u; data(model.velocities.u) .= π; time_step!(model, Nt=3, Δt=1e-16)

julia> max_u.data
3-element Array{Float64,1}:
 3.141592653589793
 3.1415926025389127
 3.1415925323439517
```
"""
function Timeseries(diagnostic, model; frequency=nothing, interval=nothing)
    TD = typeof(diagnostic(model))
    TT = typeof(model.clock.time)
    return Timeseries(diagnostic, frequency, interval, TD[], TT[])
end

@inline (timeseries::Timeseries)(model) = timeseries.diagnostic(model)

function run_diagnostic(model, diag::Timeseries)
    push!(diag.data, diag(model))
    push!(diag.time, model.clock.time)
    return nothing
end

"""
    Timeseries(diagnostics::NamedTuple, model; frequency=nothing, interval=nothing)

A `Timeseries` `Diagnostic` that records a `NamedTuple` of time series of
`diag(model)` for each `diag` in `diagnostics`.

Example
=======
```julia
julia> model = Model(grid=RegularCartesianGrid(size=(16, 16, 16), length=(1, 1, 1))); Δt = 1.0;

julia> cfl = Timeseries((adv=AdvectiveCFL(Δt), diff=DiffusiveCFL(Δt)), model; frequency=1);

julia> model.diagnostics[:cfl] = cfl; time_step!(model, Nt=3, Δt=Δt)

julia> cfl.data
(adv = [0.0, 0.0, 0.0, 0.0], diff = [0.0002688, 0.0002688, 0.0002688, 0.0002688])

julia> cfl.adv
4-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
```
"""
function Timeseries(diagnostics::NamedTuple, model; frequency=nothing, interval=nothing)
    TT = typeof(model.clock.time)
    TDs = Tuple(typeof(diag(model)) for diag in diagnostics)
    data = NamedTuple{propertynames(diagnostics)}(Tuple(T[] for T in TDs))
    return Timeseries(diagnostics, frequency, interval, data, TT[])
end

function (timeseries::Timeseries{<:NamedTuple})(model)
    names = propertynames(timeseries.diagnostic)
    return NamedTuple{names}(Tuple(diag(model) for diag in timeseries.diagnostics))
end

function run_diagnostic(model, diag::Timeseries{<:NamedTuple})
    ntuple(Val(length(diag.diagnostic))) do i
        Base.@_inline_meta
        push!(diag.data[i], diag.diagnostic[i](model))
    end
    push!(diag.time, model.clock.time)
    return nothing
end

Base.getproperty(ts::Timeseries{<:NamedTuple}, name::Symbol) = get_timeseries_field(ts, name)

"Returns a field of timeseries, or a field of `timeseries.data` when possible."
function get_timeseries_field(timeseries, name)
    if name ∈ propertynames(timeseries)
        return getfield(timeseries, name)
    else
        return getproperty(timeseries.data, name)
    end
end

"""
    FieldMaximum(mapping, field)

An object for calculating the maximum of a `mapping` function applied
element-wise to `field`.

Examples
=======
```julia
julia> model = Model(grid=RegularCartesianGrid(size=(16, 16, 16), length=(1, 1, 1)));

julia> max_abs_u = FieldMaximum(abs, model.velocities.u);

julia> max_w² = FieldMaximum(x->x^2, model.velocities.w);
```
"""
struct FieldMaximum{F, M}
    mapping :: M
      field :: F
end

(m::FieldMaximum)(args...) = maximum(m.mapping, m.field.data.parent)

(m::FieldMaximum{<:NamedTuple})(args...) =
    NamedTuple{propertynames(m.field)}(maximum(m.mapping, f.data.parent) for f in m.field)

"""
    CFL{D, S}

An object for computing the Courant-Freidrichs-Lewy (CFL) number.
"""
struct CFL{D, S}
           Δt :: D
    timescale :: S
end

"""
    CFL(Δt [, timescale=Oceananigans.cell_advection_timescale])

Returns an object for computing the Courant-Freidrichs-Lewy (CFL) number
associated with time step or `TimeStepWizard` `Δt` and `timescale`.

See also `AdvectiveCFL` and `DiffusiveCFL`
"""
CFL(Δt) = CFL(Δt, cell_advection_timescale)

(c::CFL{<:Number})(model) = c.Δt / c.timescale(model)
(c::CFL{<:TimeStepWizard})(model) = c.Δt.Δt / c.timescale(model)

"""
    AdvectiveCFL(Δt)

Returns an object for computing the Courant-Freidrichs-Lewy (CFL) number
associated with time step or `TimeStepWizard` `Δt` and the time scale
for advection across a cell.

Example
=======
```julia
julia> model = Model(grid=RegularCartesianGrid(size=(16, 16, 16), length=(8, 8, 8)));

julia> cfl = AdvectiveCFL(1.0);

julia> data(model.velocities.u) .= π;

julia> cfl(model)
6.283185307179586
```
"""
AdvectiveCFL(Δt) = CFL(Δt, cell_advection_timescale)

"""
    DiffusiveCFL(Δt)

Returns an object for computing the Courant-Freidrichs-Lewy (CFL) number
associated with time step or `TimeStepWizard` `Δt` and the time scale
for diffusion across a cell associated with `model.closure`.

Example
=======
```julia
julia> model = Model(grid=RegularCartesianGrid(size=(16, 16, 16), length=(1, 1, 1)));

julia> cfl = DiffusiveCFL(0.1);

julia> cfl(model)
2.688e-5
```
"""
DiffusiveCFL(Δt) = CFL(Δt, cell_diffusion_timescale)