text:
module YaoPastaQ

using YaoBase, YaoBlocks, PastaQ
export genlist, apply!, PastaQReg
flblock(blk::AbstractBlock) = YaoBlocks.Optimise.simplify(blk, rules=[YaoBlocks.Optimise.to_basictypes])
sublocs(subs, locs) = [locs[i] for i in subs]

function genlist(x::AbstractBlock{N}) where N
    plist = []
    genlist!(plist, flblock(x), [1:N...], Int[])
    return plist
end

function genlist!(plist, qc_simpl::ChainBlock, locs, controls)
    for block in subblocks(qc_simpl)
        genlist!(plist, block, locs, controls)
    end
end

function genlist!(plist, blk::PutBlock{N,M}, locs, controls) where {N,M}
    genlist!(plist, blk.content, sublocs(blk.locs, locs), controls)
end

genlist!(plist, ::XGate, locs, controls) = push!(plist, ("X", locs[1]))
genlist!(plist, ::YGate, locs, controls) = push!(plist, ("Y", locs[1]))
genlist!(plist, ::HGate, locs, controls) = push!(plist, ("H", locs[1]))
genlist!(plist, ::ZGate, locs, controls) = push!(plist, ("Z", locs[1]))
genlist!(plist, ::Scale{Val{im}, 1, YGate}, locs, controls) = push!(plist, ("iY", locs[1]))
genlist!(plist, ::TGate, locs, controls) = push!(plist, ("π/8", locs[1]))
genlist!(plist, ::ConstGate.SGate, locs, controls) = push!(plist, ("Phase", locs[1]))

function genlist!(plist, blk::ShiftGate{Float64}, locs, controls)
    if blk.theta == π/2
        push!(plist, ("Phase", locs[1]))
    elseif blk.theta == π/4
        push!(plist, ("π/8", locs[1]))
    end
end

genlist!(plist, blk::I2Gate, locs, controls) = push!(plist, ("Id", locs[1]))
genlist!(plist, blk::RotationGate{1, Float64, XGate}, locs, controls) = push!(plist, ("Rx", locs[1], (θ = blk.theta,)))
genlist!(plist, blk::RotationGate{1, Float64, YGate}, locs, controls) = push!(plist, ("Ry", locs[1], (θ = blk.theta,)))
genlist!(plist, blk::RotationGate{1, Float64, ZGate}, locs, controls) = push!(plist, ("Rz", locs[1], (ϕ = blk.theta,)))

function genlist!(plist, blk::ControlBlock{N, XGate, 1, 1}, locs, controls) where N
    push!(plist, ("CX", (blk.ctrl_locs[1], blk.locs[1])))
end

function genlist!(plist, blk::ControlBlock{N, YGate, 1, 1}, locs, controls) where N
    push!(plist, ("CY", (blk.ctrl_locs[1], blk.locs[1])))
end

function genlist!(plist, blk::ControlBlock{N, ZGate, 1, 1}, locs, controls) where N
    push!(plist, ("CZ", (blk.ctrl_locs[1], blk.locs[1])))
end

function genlist!(plist, blk::ControlBlock{N, RotationGate{1, Float64, ZGate}, 1, 1}, locs, controls) where N
    push!(plist, ("CRz", (blk.ctrl_locs[1], blk.locs[1]), (ϕ = blk.content.theta,)))
end

genlist!(plist, blk::SWAPGate, locs, controls) = push!(plist, ("SWAP", (locs[1], locs[2])))
genlist!(plist, blk::ControlBlock{3, XGate, 2, 1}, locs, controls) = push!(plist, ("Toffoli", (blk.ctrl_locs[1], blk.ctrl_locs[2], blk.locs[1])))
genlist!(plist, blk::ControlBlock{3, SWAPGate, 1, 2}, locs, controls) = push!(plist, ("Fredkin", (blk.ctrl_locs[1], blk.locs[1], blk.locs[2])))
genlist!(plist, blk::ControlBlock{4, XGate, 3, 1}, locs, controls) = push!(plist, ("CCCNOT", (blk.ctrl_locs[1], blk.ctrl_locs[2], blk.ctrl_locs[3], blk.locs[1])))

mutable struct PastaQReg{State <: Union{PastaQ.MPS, PastaQ.MPO}} <: AbstractRegister{1}
    state::State
end

PastaQReg(x::Int64) = PastaQReg(productstate(x))

function YaoBase.apply!(r::PastaQReg, x::AbstractBlock)
    r.state = runcircuit(r.state, genlist(x))
    return r
end

YaoBase.nqubits(r::PastaQReg) = length(r.state)
YaoBase.nactive(r::PastaQReg) = YaoBase.nqubits(r)
PastaQReg(x::YaoBlocks.BitStr) = PastaQReg(productstate(length(x), reverse(["$i" for i in x])))
PastaQReg(x::Union{PastaQ.MPS, PastaQ.MPO}) = PastaQReg(productstate(x))
Base.copy(r::PastaQReg) = PastaQReg(r)
PastaQReg(r::PastaQReg) = PastaQReg(copy(r.state))
YaoBase.fidelity(x::PastaQReg, y::PastaQReg) = PastaQ.fidelity(x.state, y.state)

function YaoBase.measure(x::PastaQReg, nshots::Int=1024)
    return getsamples(x.state, nshots)
end

end


text:
# Univariate stochastic volatility model
# The model is presented for instance in §2.3.4 of the following book and the
# proposal for the guided PF is introduced in Example 10.5
# N Chopin, O Papaspiliopoulos - An Introduction to Sequential Monte Carlo, 2020

module StochVolatility
using SequentialMonteCarlo
using RNGPool
import SMCExamples.Particles.Float64Particle

struct SVTheta
    σ::Float64
    μ::Float64
    ρ::Float64
end

mutable struct GuidedSVParticle
  x::Float64
  xprev::Float64
  GuidedSVParticle() = new()
end

function makeSVModel(theta::SVTheta, ys::Vector{Float64})
    n::Int64 = length(ys)
    log2pi = -0.5 * log(2 * π)
    @inline function lG(p::Int64, particle::Float64Particle, ::Nothing)
        @inbounds v::Float64 = particle.x + ys[p]^2 / exp(particle.x)
        return log2pi - 0.5 * v
    end
    @inline function M!(newParticle::Float64Particle, rng::RNG, p::Int64,
        particle::Float64Particle, ::Nothing)
        if p == 1
            newParticle.x = theta.μ + theta.σ*randn(rng)
        else
            newParticle.x = theta.μ + theta.ρ*(particle.x - theta.μ) + theta.σ*randn(rng)
        end
    end
    return SMCModel(M!, lG, n, Float64Particle, Nothing)
end

function makeguidedSVModel(theta::SVTheta, ys::Vector{Float64}; guide=false)
    n::Int64 = length(ys)
    log2pi = -0.5 * log(2 * π)
    @inline function computeξ(p::Int64, x::Float64)
        μstar::Float64 = if p == 1
            theta.μ
        else
            theta.μ + theta.ρ*(x - theta.μ)
        end
        @inbounds ξ::Float64 = μstar + 0.25*theta.σ^2*(ys[p]^2*exp(-μstar) - 2)
        return ξ
    end
    @inline function lG(p::Int64, particle::GuidedSVParticle, ::Nothing)
        @inbounds y::Float64 = ys[p]
        loggp::Float64 = log2pi -0.5*(particle.x + y^2/exp(particle.x))
        v::Float64 = if p == 1
            theta.μ
        else
            theta.μ + theta.ρ*(particle.xprev - theta.μ)
        end
        logfp::Float64 = log2pi -0.5*(particle.x - v)^2/theta.σ^2
        ξ::Float64 = computeξ(p, particle.xprev)
        logqp::Float64 = log2pi -0.5*(particle.x - ξ)^2/theta.σ^2
        return logfp + loggp - logqp
    end
    @inline function M!(newParticle::GuidedSVParticle, rng::RNG, p::Int64,
        particle::GuidedSVParticle, ::Nothing)
        newParticle.xprev = particle.x
        ξ::Float64 = computeξ(p, particle.x)
        newParticle.x = ξ + theta.σ*randn(rng)
    end
    return SMCModel(M!, lG, n, GuidedSVParticle, Nothing)
end

function simulateSVModel(theta::SVTheta, n::Int64)
  model = makeSVModel(theta, Vector{Float64}(undef, 0))
  ys = Vector{Float64}(undef, n)
  xParticle = Float64Particle()
  rng = getRNG()
  for p in 1:n
    model.M!(xParticle, rng, p, xParticle, nothing)
    ys[p] = exp(xParticle.x)*randn(rng)
  end
  return ys
end

function defaultSVModel(n::Int64)
    theta = SVTheta(0.2, -1, 0.97)
    ys = simulateSVModel(theta, n)

    svModel = makeSVModel(theta, ys)
    guidedSvModel = makeguidedSVModel(theta, ys)

    return svModel, guidedSvModel, theta, ys
end

end

ext: jl


text:
## ------------------------------------------------------------
# # Modify these as needed
# portname = "/dev/cu.usbserial-1420"
# baudrate = 9600

function read_port(portname, baudrate; wt = 5)

    ## ------------------------------------------------------------
    # Snippet from examples/mwe.jl
    while true

        @info("Serial reading",
            portname,
            baudrate
        )
        println()

        try
            LibSerialPort.open(portname, baudrate; mode = SP_MODE_READ) do sp
                while true
                    line = readline(sp)
                    println(line)
                end
            end
        catch err

            (err isa InterruptException) && return
            @error err
            println()

        end

        sleep(wt)
    end
end

function read_port(argv::Vector = ARGS)
    port, rate = _parse_port_args(argv)
    read_port(port, rate)
end

function _parse_port_args(argv = ARGS)
    argset = ArgParse.ArgParseSettings()
    ArgParse.@add_arg_table! argset begin
        "--port", "-p"
            help="Serial port name"
            arg_type = String
            required=true
        "--rate", "-r"
            help="Root directory"
            arg_type = Int
            required=true
    end

    parsed_args = ArgParse.parse_args(argv, argset)
    port = parsed_args["port"]
    rate = parsed_args["rate"]

    return port, rate
end


text:
using Printf
using Logging
#
#io = open("log.txt", "w+")
#logger = SimpleLogger(io)
#global_logger(logger)

struct WorkPkg
	num::Int64
end

struct Ctrl
	q::Channel{Bool}
	t::Float64
end

# Given Channels source and sink,
source = Channel{WorkPkg}(32)
sink = Channel{WorkPkg}(32)

# and a function `worker` which reads items from from source, processes the item read
# and writes a result to sink,
function worker(n::Int64, done::Ctrl)
	@info(@sprintf("worker[%d] starting -> %f ms\n", n, done.t))
    while true
	if isready(done.q)
		_ = take!(done.q)
		return 42
	end

	if isready(source)
        	data = take!(source)
		#@info(@sprintf("sending %d to sink\n", data))
		put!(sink, WorkPkg(data.num*2))    # write out result
    	end

	sleep(done.t)
	yield()

    end
    @info(@sprintf("worker [%d] ending \n", n))
end

# we can schedule `n` instances of `worker` to be active concurrently.
num_workers = 10
num_jobs = 10

qc = Array{Ctrl}(undef, num_workers)
for i = 1:size(qc)[1]
	qc[i]=Ctrl(Channel{Bool}(1), convert(Int64, floor(rand()*5)))
end

@info(@sprintf("creating %d tasks\n", num_workers))
ta = Array{Task}(undef,num_workers)
for i in 1:num_workers
	ta[i] = @task(worker(i, qc[i]))

	if !istaskstarted(ta[i])
		schedule(ta[i])
	end
end

@info(@sprintf("filling pipe...\n"))
for i = 1:num_jobs

	put!(source, WorkPkg(i))
end

@info(@sprintf("draining  pipe...\n"))
for i = 1:num_jobs
	info = take!(sink)
	@info(@sprintf("received: %d\n", info.num))
end

@info(@sprintf("quiting..\n"))
# quit all workers
for i=1:num_workers
	@info(@sprintf("signal worker %d to quit\n", i))
	put!(qc[i].q, true)
end

@info(@sprintf("waiting on done!\n"))
while true
        all_done = true

	yield()
        for i = 1:(size(ta)[1])
                if !istaskdone(ta[i])
			#@info(@sprintf("task %d isn't done\n", i))
                        all_done = false
                end
        end

        if (all_done)
		@info(@sprintf("done!\n"))
                break
        end
end

for i in 1:num_workers
	@info(@sprintf("worker %d quit, returned result: %d\n", i, ta[i].result))
end

#close(io)

max_stars_repo_name: rabarar/juliaSnippets

max_issues_repo_name: rabarar/juliaSnippets

max_forks_repo_name: rabarar/juliaSnippets


text:
# const strips = 2
# const brightness = 1

struct Star
    cardinality::String
    elevation::Int
    intensity::Int
    radius::Int
end
Star(a::String, b::String, c::String, d::String) = Star(a, parse(Int, b), parse(Int, c), parse(Int, d))

struct LED
    ind1::UInt8
    ind2::UInt8
    intensity::UInt8
    buff::Vector{UInt8}
    LED(ind1, ind2, intensity) = new(ind1, ind2, intensity, [ind1, ind2, 0x00, intensity, 0x00])
end

function index2led(c::Int, e)
    deadleds = 9
    ledsperstrip = 150
    liveleds = ledsperstrip - deadleds
    a = isodd(c) ? e : liveleds - e + 1
    b = c < 3 ? 0 : ledsperstrip
    pos = a + b - 1
    pos == (liveleds - 1)/2 ? ledsperstrip + pos : pos
end

function LED(s::Star)
    cardinals = ["NE", "SW", "SE", "NW"]
    c = findfirst(==(s.cardinality), cardinals)
    intensity = s.intensity
    map(-s.radius:s.radius) do i
        pos = index2led(c, s.elevation)
        ind = pos + i
        ind2, ind1 = reinterpret(UInt8, [UInt16(ind)])
        LED(ind1, ind2, intensity)
    end
end

parse2arduino(leds::Vector{LED}) = vcat((led.buff for led in leds)...)
parse2arduino(stars::Vector{Star}) = isempty(stars) ? zeros(UInt8, 5) : vcat((parse2arduino(LED(star)) for star in stars)...)

mutable struct LEDArduino <: AbstractArduino
    port::String
    sp::SerialPort
    pwm::Observable{Vector{UInt8}}
    function LEDArduino()
        ledport = nicolas[] ? "/dev/serial/by-id/usb-Arduino__www.arduino.cc__0043_757353036313519070B1-if00" : "/dev/serial/by-id/usb-1a86_USB2.0-Serial-if00-port0"
        sp = LibSerialPort.open(ledport, baudrate[])
        pwm = Observable(UInt8[0, 0])
        on(pwm) do x
            encode(sp, x)
        end
        new(ledport, sp, pwm)
    end
end


text:
# FWI on the 3D Overthrust model using spectral projected gradient descent
# Author: pwitte.slim@gmail.com
# Date: December 2018
#

using Statistics, Random, LinearAlgebra
using JUDI, SlimOptim, HDF5, SegyIO

# Load overthrust model
if ~isfile("$(JUDI.JUDI_DATA)/overthrust_3D_initial_model.h5")
    ftp_data("ftp://slim.gatech.edu/data/SoftwareRelease/WaveformInversion.jl/3DFWI/overthrust_3D_initial_model.h5")
end
n, d, o, m0 = read(h5open("$(JUDI.JUDI_DATA)/overthrust_3D_initial_model.h5", "r"), "n", "d", "o", "m0")

# Set up model structure
model0 = Model((n[1], n[2], n[3]), (d[1], d[2], d[3]), (o[1], o[2], o[3]), m0)

# Bound constraints
vmin = ones(Float32, model0.n) .+ 0.4f0
vmax = ones(Float32, model0.n) .+ 5.5f0
mmin = vec((1f0 ./ vmax).^2); vmax = []
mmax = vec((1f0 ./ vmin).^2); vmin = []

# Scan directory for segy files and create out-of-core data container
container = segy_scan("/path/to/shot/records/", "overthrust_3D_shot",
                     ["GroupX", "GroupY", "RecGroupElevation", "SourceSurfaceElevation", "dt"])
d_obs = judiVector(container)

# Set up source
src_geometry = Geometry(container; key = "source")
wavelet = ricker_wavelet(src_geometry.t[1], src_geometry.dt[1], 0.008)  # 8 Hz peak frequency
q = judiVector(src_geometry, wavelet)

########################################### FWI ####################################################

# Optimization parameters
fevals = 15
batchsize = 1080
count = 0
fvals = []
opt = Options(limit_m = true,   # model in only in area with receivers
              buffer_size = 500f0,  # w/ 500 m buffer zone
              optimal_checkpointing = true
              )

# Objective function for minConf library
function objective_function(x)

    # update model and save snapshot
    model0.m = reshape(x, model0.n);
    write(h5open(join(["snapshot_3D_FWI_iteration_", string(count),".h5"]), "w"), "v", sqrt.(1f0 ./ model0.m))

    # select batch
    idx = randperm(d_obs.nsrc)[1:batchsize]

    # fwi function value and gradient
    f, g = fwi_objective(model0, q[idx], d_obs[idx]; options = opt)
    g = reshape(g, model0.n); g[:, :, 1:21] .= 0f0    # reset gradient in water column to 0.
    g = .125f0*g/maximum(abs.(g))   # scale gradient to help line search

    global count; count += 1;
    global fvals; fvals = [fvals; f]
    return f, vec(g)
end

# Bound projection
ProjBound(x) = median([mmin x mmax], dims=2)

# FWI with SPG
spg_opt = spg_options(verbose = 3, maxIter = fevals, memory = 1)
x, fsave, funEvals = minConf_SPG(objective_function, vec(model0.m), ProjBound, spg_opt)

# Save results
h5open("result_3D_overthrust_spg.h5", "w") do file
    write(file, "x", sqrt.(1 ./reshape(x, model0.n)), "fsave", fsave, "fhistory", convert(Array{Float32, 1}, fvals))
end


text:
using MCTS
using POMDPModels

n_iter = 50
depth = 15
ec = 3.0

solver = MCTSSolver(n_iterations=n_iter, depth=depth, exploration_constant=ec)
mdp = GridWorld(10,10)

policy = MCTSPolicy(solver, mdp)
state = GridWorldState(1,1)

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
#
# @description Implements the Feistel block cipher algorithm.
#
# @author Alejandro Claro (alejandro.claro@gmail.com)
#
# Copyright 2017 All rights reserved.
# Use of this source code is governed by a MIT-style license that can be found in the LICENSE file.
#
export FeistelCipher
export block_size, key_size, encrypt, decrypt

# @description Defines the DES cipher data struture.
#
# f          The round function f(R, key).
# keys       The set of keys.The number of keys defines the number of iterations.
# key_size   The size of the key from input keys where derived.
# block_size The round function input block size.
struct FeistelCipher <: BlockCipher
  f::Function
  keys::Vector{Vector{UInt8}}
  key_size::UInt64
  block_size::UInt64
end

# @description Gets the cipher block size.
#
# @param {FeistelCipher} self The cipher data struture.
#
# @return {Integer} The block size in bytes.
function block_size(self::FeistelCipher)
  return self.block_size
end

# @description Gets the cipher key size.
#
# @param {DesCipher} self The cipher data struture.
#
# @return {Integer} The key size in bytes.
function key_size(self::FeistelCipher)
  return self.key_size
end

# @description Encrypts the given plaintext block with the key set at initialization.
#
# @param {FeistelCipher}     self      The cipher data struture.
# @param {Vector{UInt8}} plaintext The message to encrypt.
#
# @return {Vector{UInt8}} The resulting ciphertext block.
function encrypt(self::FeistelCipher, plaintext::Vector{UInt8})
  if length(plaintext) != block_size(self)
    throw(ArgumentError("Invalid input plaintext length."))
  end

  return compute_feistel_network(plaintext, self.keys, self.f);
end

# @description Decrypts the given ciphertext block with the key set at initialization.
#
# @param {FeistelCipher}     self       The cipher data struture.
# @param {Vector{UInt8}} ciphertext The message to decrypt.
#
# @return {Vector{UInt8}} The resulting plaintext block.
function decrypt(self::FeistelCipher, ciphertext::Vector{UInt8})
  if length(ciphertext) != block_size(self)
    throw(ArgumentError("Invalid input ciphertext length."))
  end

  return compute_feistel_network(ciphertext, reverse(self.keys), self.f);
end

# PRIVATE IMPLEMENTATION #######################################################

# @description Computes the Feistel network.
#
# @param {Vector{UInt8}}         block The block to transform.
# @param {Vector{Vector{UInt8}}} keys  The set of keys.
# @param {Function}             step   The round function f(R, key).
#
# @return {Vector{UInt8}} The Feistel network result.
function compute_feistel_network(block::Vector{UInt8}, keys::Vector{Vector{UInt8}}, round_function::Function)
  block = unpack_bits(block)
  half  = div(length(block), 2)
  ln    = block[1 : half]
  rn    = block[(half + 1) : end]

  for i in 1 : length(keys)
    r0 = copy(rn)
    f  = round_function(rn, keys[i])
    rn = map(xor, ln, f)
    ln = r0
  end

  return pack_bits(vcat(rn, ln))
end


text:
module CompressedArraysTests

using Test
using InplaceArrays.Arrays

values = [10,20,31]
ptrs = [1,2,3,3,2,2]
a = CompressedArray(values,ptrs)
r = values[ptrs]
test_array(a,r)

b = apply(-,a)
test_array(b,-r)
@test isa(b,CompressedArray)

c = apply(*,a,b)
test_array(c,r.*(-r))
@test isa(c,CompressedArray)

end # module


text:
using Cxx

include("/home/tom/.julia/v0.5/Qt5/src/Qt5.jl")
q.QApplication()
# const scene = q.Scene()

# # e1 = q.Ellipse(q.Pct(0.5), q.Pct(0.5),
# #                q.Pct(1), q.Pct(0.5),
# #                q.Pen("black"), q.Brush("red"))
# # push!(scene, e1)
# # display(scene)
# # e1.h = 0.8
# # display(scene)

# const vw = q.View(q.ViewBox(300,0,300,300), q.SceneItem[])
# push!(scene, vw)
# display(scene)

# # e2 = q.Ellipse(q.Pct(0.5), q.Pct(0.5),
# #                q.Pct(1), q.Pct(0.5),
# #                q.Pen("black"), q.Brush("blue"))
# # push!(vw, e2)
# # display(scene)

function addcircles(scene)
  pen = q.Pen("black")
  brush = q.Brush("orange")
  vw = scene.items[1]
  for i in 1:100000
     push!(vw, q.Ellipse(rand(),rand(),0.03,0.03,pen,brush))
  end
end

function dotest()
  scene = q.Scene()

  vw = q.View(q.ViewBox(300,0,300,300), q.SceneItem[])
  push!(scene, vw)
  display(scene)

  @time addcircles(scene)
  @time display(scene)
end


text:
function solve(digits, actual, easy)
    if easy
        return sum((x -> x in [2, 3, 4, 7]).((length).(actual)))
    else
        # Idea:
        #   find '1', '4', '7' and '8' through number of segments
        #
        #   find '3' by seeing which number with 5 segments has 2 in common with '1'
        #   find '2' and '5' by taking remaining 5 segment numbers, '5' has 3 segments in common with '4'
        #
        #   find '6' by seeing which number with 6 segments has 1 in common with '1'
        #   find '0' and '9' by taking remaining 6 segment numbers, '9' has 4 segments in common with '4'
        #
        sort_internal(x) = join(sort(collect(x)))
        ds = (sort_internal).(digits)
        as = (sort_internal).(actual)

        sort!(ds, by=length)
        zero, one, two, three, four, five, six, seven, eight, nine = "abcdefghij"
        one, seven, four, eight = ds[1], ds[2], ds[3], ds[end]
        for i in 4:6
            if length(intersect(ds[i], one)) == 2
                three = ds[i]
                two = ds[(i - 3) % 3 + 4]
                five = ds[(i - 2) % 3 + 4]
            end
        end
        if length(intersect(two, four)) == 3
            two, five = five, two
        end
        for i in 7:9
            if length(intersect(ds[i], one)) == 1
                six = ds[i]
                zero = ds[(i - 6) % 3 + 7]
                nine = ds[(i - 5) % 3 + 7]
            end
        end
        if length(intersect(zero, four)) == 4
            zero, nine = nine, zero
        end
        lookup = [zero, one, two, three, four, five, six, seven, eight, nine]

        ans = 0
        for a in as
            ans = 10ans + (findfirst(isequal(a), lookup) - 1)
        end
        return ans
    end
end

# 2 -> 1
# 3 -> 7
# 4 -> 4
# 5 -> 2, 3, 5
# 6 -> 6, 9, 0
# 7 -> 8

function main()
    s_easy = 0
    s_hard = 0
    for line in readlines()
        digits, actual = split(strip(line), " | ")
        digits = split(digits)
        actual = split(actual)
        s_easy += solve(digits, actual, true)
        s_hard += solve(digits, actual, false)
    end
    println(s_easy)
    println(s_hard)
end

main()

max_stars_repo_name: LiquidFun/adventofcode

max_issues_repo_name: LiquidFun/adventofcode

max_forks_repo_name: LiquidFun/adventofcode


text:
# returns vector of function value and gradient, using diff results
function value_gradient(f, x::AbstractVector, diffresults::Val{true})#  = Val(true))
    r = DiffResults.GradientResult(x) # this could take pre-allocated temporary storage
    ForwardDiff.gradient!(r, f, x)
    return vcat(r.value, r.derivs[1]) # d+1 # ... to avoid this
end
function value_gradient(f, x::AbstractVector, diffresults::Val{false} = Val(false))
    ∇ = ForwardDiff.gradient(f, x)
    return vcat(f(x), ∇)
end

function value_gradient_hessian(f, x::AbstractVector)
    d = length(x)
    dd = d^2 + d + 1
    fx = similar(x, dd)
    fx[1] = f(x)
    grad = @view fx[2:d+1]
    hess = @view fx[d+2:end]
    results = DiffResults.DiffResult(grad, hess)
    g(x) = ForwardDiff.gradient(f, x)
    ForwardDiff.jacobian!(results, g, x)
    # IDEA: try to minimize evaluations to a single pass
    # hessian = reshape(hessian, d, d) # needs to be d + 1  by d
    # result = DiffResults.DiffResult(valgrad, hessian)
    # println(result)
    # ForwardDiff.jacobian!(result, x->value_gradient(f, x), x)
    return fx

end

function value_derivative(f, x::Real, diffresults::Val{true} = Val(true))
    r = DiffResults.DiffResult(zero(x), zero(x)) # this could take pre-allocated temporary storage
    r = ForwardDiff.derivative!(r, f, x)
    return vcat(r.value, r.derivs[1])
end
function value_derivative(f, x::Real, diffresults::Val{false})
    ∇ = ForwardDiff.derivative(f, x)
    return vcat(f(x), ∇)
end

################################################################################
issquare(A::AbstractMatrix) = size(A, 1) == size(A, 2)
issquare(A::Real) = true
# converts y in various forms to a vector of numbers
_vec(Y::AbstractMatrix) = vec(Y)
_vec(y::AbstractVector{<:Number}) = y
# this needs some thinking if we want to make static arrays efficient
function _vec(y::VecOfVec{<:Number}) # different behavior to vec!
	d = length(y[1])
	u = zeros(length(y)*d)
	for i in eachindex(y)
		u[d*(i-1)+1:d*i] .= y[i]
	end
	return u
end

########################### Laplace Approximation ##############################
# IDEA: expand this into variational approximations (mean field etc.)
# IDEA: connect with GP classification
# ForwardDiff.hessian(f, x::Real) = FD.derivative(y->FD.derivative(f, y), x)
# approximates a pdf f with a normal distribution centered at x using laplace's method
# TODO: height of Gaussian?
function laplace(f, x::AbstractVector)
    H = ForwardDiff.hessian(y->log(f(y)), x)
    Gaussian(x, inverse(H))
end


text:
# using Dates

@testset "Gsl.jl" begin
    filename = "data\\test.gsl"
    dc = DDR2import.Gsl.read(filename)
    @test dc["BIRD"].ID == "BIRD"
    @test dc["BIRD"].Name == "_"
    @test dc["BIRD"].Category == "_"
    @test dc["BIRD"].Type == "FIR"
    @test dc["BIRD"].Airblocks[2].Name == "002BI"
    @test dc["BIRD"].Airblocks[2].LowerFL == 0
    @test dc["BIRD"].Airblocks[2].UpperFL == 999
    #LAST
    @test dc["BIRDBARDAR"].Type == "ES"
end


text:
# ------------------------------------------------------------------- #
# Analysis of stuff
# ------------------------------------------------------------------- #
using PyPlot

function plot_allocations(m::BCFL22C, allocations::Array{Float64, 2})
    # Create Plot Object
    fig, ax = subplots(3, 1)
    capT = size(allocations, 2)
    tTt = collect(range(1, capT))

    # Plot (a1, b2)
    ax[1][:plot](tTt, allocations[1, :][:], label="a1")
    ax[1][:plot](tTt, allocations[4, :][:], label="b2")
    ax[1][:set_title]("Home Goods")

    # Plot (a2, b1)
    ax[2][:plot](tTt, allocations[2, :][:], label="a2")
    ax[2][:plot](tTt, allocations[3, :][:], label="b1")
    ax[2][:set_title]("Foreign Goods")

    # Plot (c1, c2)
    ax[3][:plot](tTt, allocations[5, :][:], label="c1")
    ax[3][:plot](tTt, allocations[6, :][:], label="c2")
    ax[3][:set_title]("Consumption")

    fig[:tight_layout]()

    return fig, ax
end

"""
Simulate the economy and compute the allocations

NOTE: The `fsts` argument should have been initialized so that `l♠` field has
its first element set to the deisred initial condition. Also the `lzbar`,
`lzbarp`, and `lgp` fields should be set to include the entire simulation path.

The elements `fsts.l♠[2:end]` will be overwritten by this function
"""
function simulate_allocations!(m::BCFL22C, pf::PolicyFunction, fsts::FullState)

    # Allocate Space
    capT = length(fsts) - 1
    allocations = Array(Float64, 6, capT)

    for t=1:capT
        # extract time t state
        fst = fsts[t]
        # update ♠_{t+1}
        l♠p = evaluate(pf, fst)
        fsts.l♠[t+1] = l♠p

        # compute allocation
        alloc = get_allocation(m, fst)
        for (i, x) in enumerate(alloc)
            allocations[i, t] = x
        end

    end

    return allocations, fsts
end

function impulse_response(m::BCFL22C, pf::PolicyFunction, capT::Int)

    # Create a one standard deviation shock for each epsilon
    shock_at_1 = squeeze(eye(1, capT+1), 1)
    zero_shock = zeros(shock_at_1)

    # Simulate Exogenous processes
    lzbar_1, lg_1 = simulate_exog(m, shock_at_1, zero_shock)

    fsts_1 = FullState(zeros(capT+1), lzbar_1[1:capT+1],
                     [lzbar_1[2:capT+1]; NaN], [lg_1[2:capT+1]; NaN])

    lzbar_2, lg_2 = simulate_exog(m, zero_shock, shock_at_1)

    fsts_2 = FullState(zeros(capT+1), lzbar_2[1:capT+1],
                       [lzbar_2[2:capT+1]; NaN], [lg_2[2:capT+1]; NaN])

    # Get allocations
    allocations_1, fsts_1 = simulate_allocations!(m, pf, fsts_1)
    allocations_2, fsts_2 = simulate_allocations!(m, pf, fsts_2)

    # Plot allocations
    fig_1, ax_1 = plot_allocations(m, allocations_1)
    fig_1[:suptitle]("Impulse Response to e1")
    show()

    fig_2, ax_2 = plot_allocations(m, allocations_2)
    fig_2[:suptitle]("Impulse Response to e2")
    show()

    return allocations_1, allocations_2
end

max_stars_repo_name: NYUEcon/GrowthModels

max_issues_repo_name: NYUEcon/GrowthModels

max_forks_repo_name: NYUEcon/GrowthModels


text:
#!/home/matt/bin/julia

push!(LOAD_PATH, "/home/matt/Documents/julia/jsvg/")
push!(LOAD_PATH, "/home/matt/Documents/julia/3dprinting/")

using SVG
using Mesh

const THRESH = 0.0001

struct QFace
	i
	ab
	bc
	ca
	area
	function QFace(n, i)
		ai,bi,ci = abci(n, i)
		a,b,c = abc(n, i)
	        ab, bc, ca = b-a, c-b, a-c
		area = 0.5 * abs(magnitude(ab) * sin(angleXY(ab) - angleXY(ca))) * magnitude(ca)
#==
if area < 0.001
println("i ", i, " area ", area)
println(n.faces[i])
println("ai ", ai, " bi ", bi, " ci ", ci)
println("a ", a, " b ", b, " c ", c)
end
==#
		new(i, (angleXY(ab), magnitude(ab)), (angleXY(bc), magnitude(bc)), (angleXY(ca), magnitude(ca)), area)
	end
end

struct Quad
	qfaces::Tuple{QFace, QFace}
end

quad_area(q) = q.qfaces[1].area + q.qfaces[2].area

function qsvg(fid, n, q::Quad, styles, style)
	for qf in q.qfaces
		a,b,c = abc(n, qf.i)

		SVG.polyline(fid, [a, b], styles[style])
		SVG.polyline(fid, [b, c], styles[style])
		SVG.polyline(fid, [c, a], styles[style])
	end
end

function quad(n, a, b, c, d)
	s = SVG.blackline(1)
	f1 = face!(n, a, b, c, "visible", "visible", "hidden")
	f2 = face!(n, a, c, d, "hidden", "visible", "visible")
	Quad((QFace(n, f1), QFace(n, f2)))
end

function biggest_face(q)
	biga = -Inf
	bigi = 0
	for i in 1:length(q.qfaces)
		if q.qfaces[i].area > biga
			biga = q.qfaces[i].area
			bigi = i
		end
	end
	bigi
end

function random_p_on_qface(n, qface)
	a,_,c = abc(n, qface.i)
	rand() = 0.5
	p1 = a + rand() * qface.ab[2] * Vertex(cos(qface.ab[1]), sin(qface.ab[1]))
	p2 = c - p1
	vertex!(n, p1 + rand() * magnitude(p2) * Vertex(cos(angleXY(p2)), sin(angleXY(p2))))
end

function orth(n, rpi, fi)
	rp = n.vertices[rpi]
	face = n.faces[fi]
	new_faces = Vector{Face}()

	ai,bi,ci = abci(n, face)
	a,b,c = n.vertices[[ai,bi,ci]]

	function twoface!(vai, vbi, vpi)
		push!(new_faces, Face(Edge(vai, vpi, "visible"), Edge(vpi, rpi, "visible"), Edge(rpi, vpi, "hidden")))
		push!(new_faces, Face(Edge(vpi, vbi, "visible"), Edge(vbi, rpi, "hidden"), Edge(rpi, vpi, "visible")))
	end

	if face.ab.style != "hidden"
		p = intersect(a, b, rp)
		pi = vertex!(n, p)
		twoface!(ai, bi, pi)
	end

	if face.bc.style != "hidden"
		p = intersect(b, c, rp)
		pi = vertex!(n, p)
		twoface!(bi, ci, pi)
	end
	if face.ca.style != "hidden"
		p = intersect(c, a, rp)
		pi = vertex!(n, p)
		twoface!(ci, ai, pi)
	end

	faces = [fi]

	if length(new_faces) > 0
		n.faces[fi] = new_faces[1]
		for i in 2:length(new_faces)
			push!(faces, face!(n, new_faces[i]))
		end
	end
	faces

end

function intersect(a, b, c, d)

	ab = b - a
	cd = d - c

	a_ab = angleXY(ab)

	if abs(a_ab - angleXY(cd)) < THRESH
		println("parallel", a_ab - angleXY(cd) )
		return
	end

	function mag(t)
		p = a + t * ab
		magnitude(p - c) + magnitude(p - d)
	end

	intersect(a, b, mag)
end

intersect(a, b, v::Vertex) = intersect(a, b, t->magnitude((a + t * (b-a)) - v))

function intersect(a, b, mag::Function)
	ab = b - a

	s = 0
	sd = mag(s)
	f = 1.0
	fd = mag(f)
	dt = 0.5
	t = (s+f)/2
	dd = Inf
	k = 0
	while dd > THRESH && k < 1000
		td = mag(t)
		tdp = mag(t+dt)
		tdm = mag(t-dt)

		if td - tdp > 0
			s = t
			dd = abs(td - tdp)
		elseif td - tdm > 0
			f = t
			dd = abs(td - tdm)
		end
		dt /= 2
		t = (s+f)/2
		k += 1
	end
	a + t * ab
end

function divide_q(n, q)
	faces = Vector{Int}()
	rpi = random_p_on_qface(n, q.qfaces[biggest_face(q)])
	for i in 1:length(q.qfaces)
		faces = vcat(faces, orth(n, rpi, q.qfaces[i].i))
	end
	q1 = Quad((QFace(n, faces[1]), QFace(n, faces[8])))
	q2 = Quad((QFace(n, faces[2]), QFace(n, faces[3])))
	q3 = Quad((QFace(n, faces[4]), QFace(n, faces[5])))
	q4 = Quad((QFace(n, faces[6]), QFace(n, faces[7])))

	[q1, q2, q3, q4]
end

function divide(n, q)
	if quad_area(q) < 1000
		return
	end

	for qq in divide_q(n, q)
		divide(n, qq)
	end
end

fid = open("test.svg", "w+")
SVG.open(fid, 1350, 1350)

styles = Dict{String, String}()
styles["visible"] = SVG.blackline(1)
styles["hidden"] = SVG.style("none", "red", 0.1)
styles["blue"] = SVG.style("none", "blue", 3)
styles["green"] = SVG.style("none", "green", 3)
styles["magenta"] = SVG.style("none", "magenta", 3)
styles["red"] = SVG.style("none", "red", 3)
styles[""] = ""

n = Net()
vertex!(n, 50,50)
vertex!(n, 50,300)
vertex!(n, 300,300)
vertex!(n, 300,50)

vertex!(n, 250,790)
vertex!(n, 394,400)
vertex!(n, 936,657)
vertex!(n, 678,931)

q1 = quad(n, 1, 2, 3, 4)
q2 = quad(n, 5, 6, 7, 8)

divide(n, q1)

divide(n, q2)

SVG.meshXY(fid, n, styles)

SVG.close(fid)
close(fid)

max_stars_repo_name: lawless-m/jsvg


text:
abstract type AbstractIdData end

const AnyInput = Union{AbstractArray,AbstractIdData}

struct InputOutputData{Y,U,T} <: AbstractIdData
    y::Y
    u::U
    Ts::T
end

struct OutputData{Y,T} <: AbstractIdData
    y::Y
    Ts::T
end

struct InputOutputStateData{Y,U,X,T} <: AbstractIdData
    y::Y
    u::U
    x::X
    Ts::T
end

autodim(x::Vector{<:AbstractVector}) = x
autodim(x::AbstractVector) = x'
function autodim(x)
    r = size(x, 1)
    c = size(x, 2)
    if (c < 5 && c < r) || (r > 4c)
        @info "Transposing input. The convention used in ControlSystemIdentification is that input-output data is made out of either of 1) Vectors with scalars, 2) vectors of vectors or 3) matrices with time along the second dimension. The supplied input appears to be multidimensional and have time in the first dimension." maxlog =
            3
        return copy(x')
    end
    x
end

function Base.show(io::IO, d::OutputData)
    write(io, "Output data of length $(length(d)) with $(noutputs(d)) outputs")
end
function Base.show(io::IO, d::InputOutputData)
    write(
        io,
        "InputOutput data of length $(length(d)) with $(noutputs(d)) outputs and $(ninputs(d)) inputs",
    )
end

iddata(y::AbstractArray, Ts::Union{Real,Nothing} = nothing) = OutputData(autodim(y), Ts)
iddata(y::AbstractArray, u::AbstractArray, Ts::Union{Real,Nothing} = nothing) =
    InputOutputData(autodim(y), autodim(u), Ts)

"""
    iddata(y, u, x, Ts = nothing)

Returns the appropriate IdData object, depending on the input.

# Arguments
- `y::AbstractArray`: output data
- `u::AbstractArray`: input data
- `x::AbstractArray`: state data
- `Ts::Union{Real,Nothing} = nothing`: optional sample time

# Examples
```jldoctest
julia> iddata(randn(10))
Output data of length 10 with 1 outputs

julia> iddata(randn(10), randn(10), 1)
InputOutput data of length 10 with 1 outputs and 1 inputs
```
"""
iddata(
    y::AbstractArray,
    u::AbstractArray,
    x::AbstractArray,
    Ts::Union{Real,Nothing} = nothing,
) = InputOutputStateData(autodim(y), autodim(u), x, Ts)

output(d::AbstractIdData)                        = d.y
input(d::AbstractIdData)                         = d.u
LowLevelParticleFilters.state(d::AbstractIdData) = d.x
output(d::AbstractArray)                         = d
input(d::AbstractArray)                          = d
LowLevelParticleFilters.state(d::AbstractArray)  = d
hasinput(::OutputData)                           = false
hasinput(::AbstractIdData)                       = true
hasinput(::AbstractArray)                        = true
hasinput(::ControlSystems.LTISystem)             = true
ControlSystems.noutputs(d::AbstractIdData)       = obslength(d.y)
ControlSystems.ninputs(d::AbstractIdData)        = hasinput(d) ? obslength(d.u) : 0
ControlSystems.nstates(d::AbstractIdData)        = 0
ControlSystems.nstates(d::InputOutputStateData)  = obslength(d.x)
obslength(d::AbstractIdData)                     = ControlSystems.noutputs(d)
sampletime(d::AbstractIdData)                    = d.Ts === nothing ? 1.0 : d.Ts
function Base.length(d::AbstractIdData)
    y = output(d)
    y isa AbstractMatrix && return size(y, 2)
    return length(y)
end

Base.lastindex(d::AbstractIdData) = length(d)

function Base.getproperty(d::AbstractIdData, s::Symbol)
    if s === :fs || s === :Fs
        return 1 / d.Ts
    end
    return getfield(d, s)
end

function Base.:(==)(d1::T, d2::T) where {T<:AbstractIdData}
    all(fieldnames(T)) do field
        getfield(d1, field) == getfield(d2, field)
    end
end

timevec(d::AbstractIdData) = range(0, step = sampletime(d), length = length(d))
timevec(d::AbstractVector, h::Real) = range(0, step = h, length = length(d))
timevec(d::AbstractMatrix, h::Real) = range(0, step = h, length = maximum(size(d)))

function apply_fun(fun, d::OutputData, Ts = d.Ts)
    iddata(fun(d.y), Ts)
end

"""
	apply_fun(fun, d::InputOutputData)

Apply `fun(y)` to all time series `y[,u,[x]] ∈ d` and return a new `iddata` with the transformed series.
"""
function apply_fun(fun, d::InputOutputData, Ts = d.Ts)
    iddata(fun(d.y), fun(d.u), Ts)
end

function apply_fun(fun, d::InputOutputStateData, Ts = d.Ts)
    iddata(fun(d.y), fun(d.u), fun(d.x), Ts)
end

torange(x::Number) = x:x
torange(x) = x

function Base.getindex(d::Union{InputOutputData,InputOutputStateData}, i, j)
    iddata(d.y[torange(i), :], d.u[torange(j), :], d.Ts)
end

function Base.getindex(d::AbstractIdData, i)
    apply_fun(d) do y
        y[:, i]
    end
end

"""
dr = resample(d::InputOutputData, f)

Resample iddata `d` with fraction `f`, e.g., `f = fs_new / fs_original`.
"""
function DSP.resample(d::AbstractIdData, f)
    apply_fun(d, d.Ts / f) do y
        yr = mapslices(y, dims = 2) do y
            resample(y, f)
        end
        yr
    end
end

function DSP.resample(M::AbstractMatrix, f)
    mapslices(M, dims = 1) do y
        resample(y, f)
    end
end

function Base.hcat(d1::InputOutputData, d2::InputOutputData)
    @assert d1.Ts == d2.Ts
    iddata([d1.y d2.y], [d1.u d2.u], d1.Ts)
end

struct StateSpaceNoise{T,MT<:AbstractMatrix{T}} <:
       ControlSystems.AbstractStateSpace{Discrete{Float64}}
    A::MT
    B::MT
    K::MT
    Ts::Discrete{Float64}
    nx::Int
    nu::Int
    ny::Int
    function StateSpaceNoise(
        A::MT,
        B::MT,
        K::MT,
        Ts::Union{Real,Discrete{Float64}},
    ) where {MT}
        Ts = Ts isa Real ? Discrete(Float64(Ts)) : Float64(Ts)
        nx, nu, ny = ControlSystems.state_space_validation(
            A,
            B,
            K',
            zeros(size(K', 1), size(B, 2)),
        )
        new{eltype(A),typeof(A)}(A, B, K, Ts, nx, nu, ny)
    end
end

ControlSystems.isstable(s::StateSpaceNoise) =
    all(abs(e) <= 1 for e in eigvals(s.A - s.K * s.C))

# Funtions for number of intputs, outputs and states
# ControlSystems.ninputs(sys::StateSpaceNoise) = sys.nu
# ControlSystems.noutputs(sys::StateSpaceNoise) = sys.ny
# ControlSystems.nstates(sys::StateSpaceNoise) = sys.nx
#
# Base.ndims(::StateSpaceNoise) = 2 # NOTE: Also for SISO systems?
# Base.size(sys::StateSpaceNoise) = (noutputs(sys), ninputs(sys)) # NOTE: or just size(get_D(sys))
# Base.size(sys::StateSpaceNoise, d) = d <= 2 ? size(sys)[d] : 1
Base.eltype(::Type{S}) where {S<:StateSpaceNoise} = S
ControlSystems.numeric_type(::Type{<:StateSpaceNoise{T}}) where {T} = T
Base.convert(::Type{StateSpace}, sys::StateSpaceNoise) = ss(sys.A, sys.B, sys.C, 0, sys.Ts)
ControlSystems.ss(sys::StateSpaceNoise) = convert(StateSpace, sys)
ControlSystems.tf(sys::StateSpaceNoise) = tf(ss(sys))

function Base.getproperty(sys::StateSpaceNoise, p::Symbol)
    if p === :C
        return [I zeros(sys.ny, sys.nx - sys.ny)]
    elseif p === :D
        return zeros(sys.ny, sys.nu)
    elseif p === :timeevol
        return sys.Ts
    end
    return getfield(sys, p)
end

function Base.getindex(sys::StateSpaceNoise, inds...)
    if size(inds, 1) != 2
        error("Must specify 2 indices to index statespace model")
    end
    rows, cols = ControlSystems.index2range(inds...) # FIXME: ControlSystems.index2range(inds...)
    return StateSpaceNoise(copy(sys.A), sys.B[:, cols], sys.K[:, rows], sys.Ts)
end

struct SysFilter{T<:AbstractStateSpace{<:Discrete},FT}
    sys::T
    state::Vector{FT}
    yh::Vector{FT}
end
SysFilter(sys::LTISystem, x0 = zeros(sys.nx)) =
    SysFilter(sys, x0, zeros(eltype(x0), sys.ny))

(s::SysFilter)(y, u) = sysfilter!(s.state, s.sys, y, u)
(s::SysFilter)(u) = sysfilter!(s.state, s.sys, u)
sysfilter!(s::SysFilter, y, u) = sysfilter!(s.state, s.sys, y, u)
sysfilter!(s::SysFilter, u) = sysfilter!(s.state, s.sys, u)

function sysfilter!(state::AbstractVector, sys::StateSpaceNoise, y, u)
    @unpack A, B, K, ny = sys
    yh = state[1:ny] #vec(sys.C*state)
    e = y .- yh
    state .= vec(A * state + B * u + K * e)
    yh
end

function sysfilter!(state::AbstractVector, sys::StateSpaceNoise, u)
    @unpack A, B, K, ny = sys
    yh = state[1:ny] #vec(C*state)
    state .= vec(A * state + B * u)
    yh
end

function sysfilter!(state::AbstractVector, sys::StateSpace, y, u)
    @unpack A, B, C, D = sys
    yh = vec(C * state + D * u)
    state .= vec(A * state + B * u)
    yh
end

function sysfilter!(state::AbstractVector, sys::StateSpace, u)
    @unpack A, B, C, D = sys
    yh = vec(C * state + D * u)
    state .= vec(A * state + B * u)
    yh
end

struct OberservationIterator{T}
    y::T
    u::T
end

observations(y, u) = OberservationIterator(y, u)

function Base.iterate(i::OberservationIterator{<:AbstractMatrix}, state = 1)
    state > length(i) && return nothing
    ((i.y[:, state], i.u[:, state]), state + 1)
end
Base.length(i::OberservationIterator{<:AbstractMatrix}) = size(i.y, 2)

function Base.iterate(
    i::OberservationIterator{<:AbstractVector{<:Union{AbstractVector,Number}}},
    state = 1,
)
    state > length(i) && return nothing
    ((i.y[state], i.u[state]), state + 1)
end
Base.length(i::OberservationIterator{<:AbstractVector{<:Union{AbstractVector,Number}}}) =
    length(i.y)

struct PredictionErrorIterator{T,MT}
    model::MT
    oi::OberservationIterator{T}
end

prediction_errors(model, y, u) = PredictionErrorIterator(model, observations(y, u))

function Base.iterate(i::PredictionErrorIterator, state = 1)
    state >= length(i) && return nothing
    (y, u), state1 = iterate(i.oi, state)
    yh = i.model(y, u)
    (y .- yh, state1)
end
Base.length(i::PredictionErrorIterator) = length(i.oi)

mutable struct SimulationErrorIterator{T,MT,YT}
    model::MT
    oi::OberservationIterator{T}
    yh::YT # Stores the last prediction
end

simulation_errors(model, y, u) = SimulationErrorIterator(
    model,
    observations(y, u),
    zeros(eltype(model.sys.A), obslength(y)),
)

function Base.iterate(i::SimulationErrorIterator, state = 1)
    state >= length(i) && return nothing
    (y, u), state1 = iterate(i.oi, state)
    i.yh = i.model(u)
    (y - i.yh, state1)
end
Base.length(i::SimulationErrorIterator) = length(i.oi)

@recipe function plot(d::AbstractIdData)
    y = time1(output(d))
    n = noutputs(d)
    if hasinput(d)
        u = time1(input(d))
        n += ninputs(d)
    end
    layout --> (n, 1)
    legend --> false
    xguide --> "Time"
    link --> :x
    xvec = range(0, step = sampletime(d), length = length(d))

    for i = 1:size(y, 2)
        @series begin
            title --> "Output $i"
            label --> "Output $i"
            xvec, y[:, i]
        end
    end
    if hasinput(d)
        for i = 1:size(u, 2)
            @series begin
                title --> "Input $i"
                label --> "Input $i"
                xvec, u[:, i]
            end
        end
    end
end


text:
module Splines
    using EllipsisNotation

    include("Utils.jl")
    include("Polynomials.jl")
    include("Visuals.jl")

    import .Polynomials: bernstein, N_basis
    import .Visuals: plot_curve, plot_nodes, plot_surface
    import .Utils: get_curve_parameters, ∑, dimensions, nnodes

    export plot_curve, plot_nodes, plot_surface
    export Bezier, BSpline, AbstractSpline, AbstractBezier
    export ⊗, ∑

    # ---------------------------------------------------------------------------- #
    #                                Abstract types                                #
    # ---------------------------------------------------------------------------- #
    abstract type AbstractSpline end

    abstract type AbstractBezier <: AbstractSpline end
    abstract type AbstractBSpline <: AbstractSpline end

    # ---------------------------------- Methods --------------------------------- #

    function Base.show(io::IO, s::AbstractSpline)
        print(io, "Spline (nodes: $(s.η), N: $(s.N))")
    end

    """
        Spline(t, v, ...)

    Evaluates an N dimensional spline (with N > 1) at a set of N paramters t, v, ...
    values with each paramter `t ∈ [0, 1]`
    """
    function (bez::AbstractSpline)(params...)
        if length(params) != bez.N
            error("Number of parameters doesnt match the dimensionality of the spline surface")
        end

        sum((n)->bez.inner_fn[n](params[n]), 1:length(params))
    end

    # ---------------------------------------------------------------------------- #
    #                                   BSPLINES                                   #
    # ---------------------------------------------------------------------------- #

    struct BSpline <: AbstractBSpline
        nodes::AbstractArray            # coordinates of nodes (d x N (x M ...))
        coordinates::AbstractArray      # coordinates of points along the spline
        weights::AbstractArray          # weights, shape N x M ...
        knots::AbstractArray
        degree::Int
        η                               # shape of the nodes -1  | degree (linera, quadratic, cubic..)
        d::Int                          # dimensionality of the Euclidean space
        N::Int                          # spline dimensionality: 1 - curve, 2 - surface...
        inner_fn                        # a collection of lower dimensionality bezier curves
    end

    # --------------------------- knots initialization --------------------------- #
    """
        uniform(n, d)

    Defines d+n+1 uniform knots ∈[0,1] for b-spline interpolation.
    For i ∈ 0:n+d+1:
    if i ≤ d: k=0
    if d+1 ≤ i ≤ n: k=(i-d)/(n+1-d)
    else: k = 1
    See: https://www.geometrictools.com/Documentation/BSplineCurveLeastSquaresFit.pdf
    """
    function uniform(n::Int, d::Int)::AbstractArray
        if d > n
            @warn "When using B-splines, the degree should be < the number of control nodes"
        end

        knots = zeros(n+d+1+1)  # second 1 is because start from 0
        for i in 0:n+d+1
            if i ≤ d
                knots[i+1] = 0
            elseif (d+1) ≤ i ≤ n
                knots[i+1] = (i-d)/(n+1-d)
            else
                knots[i+1] = 1
            end
        end
        return knots
    end

    """
        periodic(n, d)

    Defines d+n+1 periodic knots ∈[0,1] for b-spline interpolation.
    See: https://www.geometrictools.com/Documentation/BSplineCurveLeastSquaresFit.pdf
    """
    function periodic(n::Int, d::Int)::AbstractArray
        if d > n
            @warn "When using B-splines, the degree should be < the number of control nodes"
        end
        map((i)->(i-d)/(n+1-d), 0:n+d+1)
    end

    # ------------------------------- Constructors ------------------------------- #
    """
        BSpline(nodes::AbstractArray; knots_type::Symbol=:uniform, δt=0.025)::BSpline

    Constructor for a (non rational) B-Spline curve (N=1).
    """
    function BSpline(nodes::AbstractArray, degree::Int; knots_type::Symbol=:uniform, δt=0.025)::BSpline
        weights = ones(size(nodes)...)
        return BSpline(nodes, weights, degree;  knots_type= knots_type, δt=δt)
    end

    """
        BSpline(nodes::AbstractArray, weights::AbstractArray; knots_type::Symbol=:uniform δt=0.025)::BSpline

    Constructor for a rational B-Spline curve (N=1).
    """
    function BSpline(nodes::AbstractArray, weights::AbstractArray, degree::Int; knots_type::Symbol=:uniform, δt=0.025)::BSpline
        # get knots
        η = size(nodes, 2) - 1
        knots = eval(:($knots_type($η, $degree)))

        return BSpline(nodes, weights, knots, degree;  δt= δt)
    end

    """
        BSpline(nodes::AbstractArray, weights::AbstractArray, knots::AbstractArray; knots_type::Symbol=:uniform δt=0.025)::BSpline

    Constructor for a rational B Spline given a set of knot
    """
    function BSpline(nodes::AbstractArray, weights::AbstractArray, knots::AbstractArray, degree::Int; δt=0.025)::BSpline
        # get params
        τ, N, η, d = get_curve_parameters(nodes, δt)

        # compute coordinates
        numerator(i) = N_basis(τ; k=knots, i=i, d=degree)' .* nodes[:, i+1] .* weights[i+1]
        denominator(i) =  N_basis(τ; k=knots, i=i, d=degree)' .* weights[i+1]

        coordinates = ∑(numerator, 0:η) ./ ∑(denominator, 0:η)

        return BSpline(
            nodes, coordinates, weights, knots, degree, η, d, N, nothing
        )
    end

    """
        BSpline(t)

    Function to evaluate  a B Spline curve at a parameter `t ∈ [0, 1]`
    """
    function (bspline::AbstractBSpline)(t)
        numerator(i) = N_basis(τ; k=bspline.knots, i=i, d=bspline.degree)' .* bspline.nodes[:, i+1] .* bspline.weights[i+1]
        denominator(i) =  N_basis(τ; k=bspline.knots, i=i, d=bspline.degree)' .* bspline.weights[i+1]

        ∑(numerator, 0:bspline.η) ./ ∑(denominator, 0:bspline.η)

    end

    # ---------------------------------------------------------------------------- #
    #                                    BEZIER                                    #
    # ---------------------------------------------------------------------------- #

    """
        Bezier

    Struct of type Bezier <: AbstractBezier.
    Contains all information about an N-dimensional Bezier surface.

    For `deg=1` we have a Bezier curve and:
        `nodes`: `dxN` array of nodes coordinates
        `η`: dimensionality of the nodes (`N-1`)
        `coordinates`: d x K points with the coordinates of the Bezier curve.

    for `deg>1` the type represents `deg` dimensional Bezier surfaces
    """
    struct Bezier <: AbstractBezier
        nodes::AbstractArray            # coordinates of nodes (d x N (x M ...))
        coordinates::AbstractArray      # coordinates of points along the spline
        weights::AbstractArray          # weights, shape N x M ...
        η                               # shape of the nodes -1  | degree (linera, quadratic, cubic..)
        d::Int                          # dimensionality of the Euclidean space
        N::Int                          # spline dimensionality: 1 - curve, 2 - surface...
        inner_fn                        # a collection of lower dimensionality bezier curves
    end

    # ------------------------------- constructors ------------------------------- #
    """
        Bezier(nodes; δt=0.05)

    Constructor for a degree 1 (`n=1`) `Bezier` curve defined by a set of nodes
    """
    function Bezier(nodes::AbstractArray; δt=0.025)::Bezier
        weights = ones(size(nodes)...)
        return Bezier(nodes, weights; δt=δt)
    end

    """
        Bezier(nodes::AbstractArray, weights::AbstractArray; δt = 0.025)

        Constructor for a degree 1 rational bezier.
    """
    function Bezier(nodes::AbstractArray, weights::AbstractArray; δt = 0.025)::Bezier
        τ, N, η, d = get_curve_parameters(nodes, δt)

        numerator(i) = bernstein(τ; i=i, n=η)' .* nodes[:, i+1] .* weights[i+1]
        denominator(i) = bernstein(τ; i=i, n=η)' .* weights[i+1]
        coordinates = ∑(numerator, 0:η) ./ ∑(denominator, 0:η)

        return Bezier(
            nodes, coordinates, weights, η, d, N, nothing
        )
    end

    # ---------------------------- function evaluation --------------------------- #∂
    """
        Bezier(t)

    Evaluate a `Bezier` curve (N=1) at parameter value `t ∈ [0,1]`
    """
    function (bez::AbstractBezier)(t)
        numerator(i) = bernstein(t; i=i, n=bez.η)' .* bez.nodes[:, i+1] .* bez.weights[i+1]
        denominator(i) = bernstein(t; i=i, n=bez.η)' .* bez.weights[i+1]

        ∑(numerator, 0:bez.η) ./ ∑(denominator, 0:bez.η)
    end

    # ---------------------------------------------------------------------------- #
    #                                Tensor Product                                #
    # ---------------------------------------------------------------------------- #

    """
        ×(x::AbstractArray, y::AbstractArray)

    Cartesian product between two arrays, without the 3 Euclidean coordinates.
    Given a N x M... array `x` and a P x Q array `y` it returns an array of shape
    N x M... x P x Q... with the product of the two arrays.
    """
    function  ×(x::AbstractArray, y::AbstractArray)::AbstractArray
        if ndims(x) == ndims(y) == 1
            return x .* y'
        else
            # initialize an array of the appropriate shape
            N = zeros(size(x)..., size(y)...)

            for Ix in CartesianIndices(x), Iy in CartesianIndices(y)
                N[Ix, Iy] = x[Ix] .+ y[Iy]
            end
        end
        return N
    end

    """
        ⊗(x::AbstractArray, y::AbstractArray)

    Tensor product between two arrays.
    Given x as a d x N x M... and y as a d x P x Q...
    it returns an array of shape d x N x M .. x P x Q x ...
    with the product of the two input arrays.

    This is done by iterating over the dimension d using Einstein summatin.
    Then each N x M... and P x Q... array is repeated/permuted to
    to create two N x M... x P x Q... arrays which are then summed.
    """
    function ⊗(x::AbstractArray, y::AbstractArray)::AbstractArray
        if ndims(x) == ndims(y) == 1
            return x .* y'
        else
            # initialize an array of the appropriate shape
            N = zeros(3, dimensions(x) ..., dimensions(y)...)

            for Ix in CartesianIndices(x[1, ..]), Iy in CartesianIndices(y[1, ..])
                N[:, Ix, Iy] = x[:, Ix] .+ y[:, Iy]
            end
        end
        return N
    end

    """
        ⊗(b1::AbstractBezier, b2::AbstractBezier)::AbstractBezier

    Tensor product of two beziers
    """
    function ⊗(b1::T, b2::T) where T <: AbstractBezier
        # get nodes
        nodes = b1.nodes ⊗ b2.nodes
        n, m = nnodes(nodes)

        # get weights
        weights = b1.weights ⊗ b2.weights

        # get coordinates
        coordinates = b1.coordinates ⊗ b2.coordinates

        # concatenate inner functions
        β1 = b1.inner_fn == nothing ? [b1] : b1.inner_fn
        β2 = b2.inner_fn == nothing ? [b2] : b2.inner_fn

        return T(nodes, coordinates, weights, (n, m), b1.d, b1.N+b2.N, hcat(β1, β2))
    end

    """
        ⊗(b1::AbstractBSpline, b2::AbstractBSpline)::AbstractBSpline

    Tensor product of two beziers
    """
    function ⊗(b1::T, b2::T) where T <: AbstractBSpline
        if b1.degree != b2.degree
            @warn "BSpline ⊗ BSpline, '⊗' is only defined for B Splines with the same degree"
        end

        # get nodes
        nodes = b1.nodes ⊗ b2.nodes
        n, m = nnodes(nodes)

        # get knots
        knots = b1.knots × b2.knots

        # get weights
        weights = b1.weights ⊗ b2.weights

        # get coordinates
        coordinates = b1.coordinates ⊗ b2.coordinates

        # concatenate inner functions
        β1 = b1.inner_fn == nothing ? [b1] : b1.inner_fn
        β2 = b2.inner_fn == nothing ? [b2] : b2.inner_fn

        return T(nodes, coordinates, weights, knots, b1.degree, (n, m), b1.d, b1.N+b2.N, hcat(β1, β2))
    end

end

max_stars_repo_name: FedeClaudi/Splines

max_issues_repo_name: FedeClaudi/Splines

max_forks_repo_name: FedeClaudi/Splines


text:
__precompile__()

module MarkableIntegers

export Marked, Unmarked,
       @mark!, @unmark!,
       ismarked, isunmarked,
       find_marked, find_unmarked, all_marked, all_unmarked,
       MarkInt128, MarkInt64, MarkInt32, MarkInt16, MarkInt8, MarkInt,
       MarkUInt128, MarkUInt64, MarkUInt32, MarkUInt16, MarkUInt8, MarkUInt,
       MarkableSigned, MarkableUnsigned, MarkableInteger

import Base: @pure, promote_type, promote_rule, convert, mark, unmark, ismarked,
    Integer, Signed, Unsigned, signed, unsigned,
    leading_zeros, trailing_zeros, leading_ones, trailing_ones,
    convert, promote_rule, string, show,
    (<=), (<), (==), (!=), (>=), (>), isless, isequal,
    (~), (&), (|), (⊻),
    findall

import Base.Math: zero, one, iszero, isone, isinteger, typemax, typemin,
    isodd, iseven, sign, signbit, abs, copysign, flipsign,
    (+), (-), (*), (/), (\), (^), div, fld, cld, mod, rem,
    sqrt, cbrt

import Base.Checked:  add_with_overflow, sub_with_overflow, mul_with_overflow,
    checked_neg, checked_abs, checked_add, checked_sub, checked_mul,
    checked_div, checked_fld, checked_cld, checked_rem, checked_div

include("type.jl")
include("promote.jl")
include("notation.jl")
include("bits.jl")
include("compare.jl")
include("ops.jl")

end # module MarkableIntegers


text:
"""
    simplify(expr::Expr)

Simplify expression.
"""
function simplify(ex::Expr)
    @assert ex.head == :call
    return simplify(Val{ex.args[1]}, ex)
end

"""
    simplify(f::Union{Number, Symbol})
"""
function simplify(f::Union{Number, Symbol})
    return f
end

"""
    simplify(Val{:*}, ex::Expr)
"""
function simplify(::Type{Val{:*}}, ex::Expr)
    args = simplify.(ex.args[2:end])
    0 in args && return 0
    filter!(k -> !(isa(k, Number) && k == 1), args)
    length(args) == 1 && return first(args)
    return Expr(:call, :*, args...)
end

"""
    simplify(Val{:+}, ex::Expr)
"""
function simplify(::Type{Val{:+}}, ex::Expr)
    args = simplify.(ex.args[2:end])
    filter!(k -> !isa(k, Number) || k != 0, args)
    length(args) == 1 && return first(args)
    return Expr(:call, :+, args...)
end

"""
    simplify(Val{:-}, ex::Expr)
"""
function simplify(::Type{Val{:-}}, ex::Expr)
    args = simplify.(ex.args[2:end])
    filter!(k -> !isa(k, Number) || k != 0, args)
    length(args) == 1 && return first(args)
    return Expr(:call, :-, args...)
end

"""
    simplify(Val{:^}, ex::Expr)
"""
function simplify(::Type{Val{:^}}, ex::Expr)
    args = simplify.(ex.args[2:end])
    return Expr(:call, :^, args...)
end

"""
    simplify(Val{:/}, ex::Expr)
"""
function simplify(::Type{Val{:/}}, ex::Expr)
    args = simplify.(ex.args[2:end])
    return Expr(:call, :/, args...)
end


text:
# This file was generated, do not modify it. # hide
ŷ = predict(tuned_ensemble, rows=test)
rms(ŷ, y[test])

PyPlot.close_figs() # hide


text:
using DSP, Test

# Naive rational resampler
function naivefilt(h::Vector, x::Vector, resamplerate::Union{Integer, Rational}=1)

    upfactor     = numerator(resamplerate)
    downfactor   = denominator(resamplerate)
    xLen         = length(x)
    xZeroStuffed = zeros(eltype(x), length(x) * upfactor)

    for n in 0:length(x)-1
        xZeroStuffed[n*upfactor+1] = x[n+1]
    end

    y = filt(h, one(eltype(x)), xZeroStuffed)
    y = [y[n] for n = 1:downfactor:length(y)]
end

# Naive arbitrary resampler
function naivefilt(h::Vector, x::Vector, resamplerate::AbstractFloat, numfilters::Integer=32)
    xLen          = length(x)
    xInterpolated = naivefilt(h, x, numfilters//1)
    xLen          = length(xInterpolated)
    yLen          = ceil(Int, xLen * resamplerate)
    y             = similar(x, yLen)
    yIdx          = 1
    xIdx          = 1
    α             = 0.0
    (δ, ϕStride)  = modf(numfilters/resamplerate)
    ϕStride       = convert(Int, ϕStride)

    while xIdx < xLen
        yLower  = xInterpolated[xIdx]
        yUpper  = xInterpolated[xIdx+1]
        y[yIdx] = yLower + α*(yUpper - yLower)
        yIdx   += 1
        α      += δ
        xIdx   += floor(Int, α) + ϕStride
        α       = mod(α, 1.0)
    end

    resize!(y, yIdx-1)

    return y
end

# Disable time and printf macros when not running interactivly ( for travis )
if isinteractive()
    macro timeifinteractive(ex)
        :(@time $(esc(ex)))
    end

    macro printfifinteractive(args...)
        :(@printf($(map(esc, args)...)))
    end
else
    macro timeifinteractive(ex)
        esc(ex)
    end

    macro printfifinteractive(args...)
    end
end

#
# Single rate filtering
#

function test_singlerate(h, x)
    xLen       = length(x)
    hLen       = length(h)
    pivotPoint = min(rand(50:150), div(xLen, 4))
    x1         = x[1:pivotPoint]
    x2         = x[pivotPoint+1:end]

    @printfifinteractive( "\n\n" )
    @printfifinteractive( "____ _ _  _ ____ _    ____    ____ ____ ___ ____\n" )
    @printfifinteractive( "[__  | |\\ | | __ |    |___    |__/ |__|  |  |___\n" )
    @printfifinteractive( "___] | | \\| |__] |___ |___    |  \\ |  |  |  |___\n" )
    @printfifinteractive( "\nTesting single-rate fitering, h is %s, x is %s. xLen = %d, hLen = %d", string(eltype(h)), string(eltype(x)), xLen, hLen )

    @printfifinteractive( "\n\tBase.filt\n\t\t")
    @timeifinteractive naiveResult = filt(h, 1.0, x)

    @printfifinteractive( "\n\tDSP.filt( h, x, 1//1 )\n\t\t" )
    @timeifinteractive statelesResult = DSP.filt( h, x )
    @test naiveResult ≈ statelesResult

    @printfifinteractive( "\n\tDSP.filt. length(x1) = %d, length(x2) = %d\n\t\t", length(x1), length(x2) )
    myfilt = DSP.FIRFilter(h, 1//1)
    @timeifinteractive begin
        y1 = DSP.filt(myfilt, x1)
        y2 = DSP.filt(myfilt, x2)
    end
    statefulResult = [y1; y2]
    @test naiveResult ≈ statefulResult

    @printfifinteractive( "\n\tDSP.filt filt. Piecewise for first %d inputs\n\t\t", length(x1) )
    DSP.reset!(myfilt)
    @timeifinteractive begin
        for i in 1:length(x1)
            y1[i] = DSP.filt(myfilt, x1[i:i])[1]
        end
        y2 = DSP.filt(myfilt, x2)
    end
    piecewiseResult = [y1; y2]
    @test naiveResult ≈ piecewiseResult

    DSP.reset!(myfilt)
    @test inputlength(myfilt, length(piecewiseResult)) == xLen
end

#
# Decimation
#

function test_decimation(h, x, decimation)
    xLen       = length(x)
    hLen       = length(h)
    pivotPoint = min(rand(50:150), div(xLen, 4))
    x1         = x[1:pivotPoint]
    x2         = x[pivotPoint+1:end]

    @printfifinteractive( "\n\n" )
    @printfifinteractive( "___  ____ ____ _ _  _ ____ ___ _ ____ _  _ \n" )
    @printfifinteractive( "|  \\ |___ |    | |\\/| |__|  |  | |  | |\\ | \n" )
    @printfifinteractive( "|__/ |___ |___ | |  | |  |  |  | |__| | \\| \n" )
    @printfifinteractive( "\nTesting decimation. h::%s, x::%s. xLen = %d, hLen = %d, decimation = %d", string(typeof(h)), string(typeof(h)), xLen, hLen, decimation )

    @printfifinteractive( "\n\tNaive decimation\n\t\t")
    @timeifinteractive naiveResult = naivefilt(h, x, 1//decimation)

    @printfifinteractive( "\n\tDSP.filt( h, x, 1//%d)\n\t\t", decimation )
    @timeifinteractive statelesResult = DSP.filt(h, x, 1//decimation)
    @test naiveResult ≈ statelesResult

    @printfifinteractive( "\n\tDSP.filt decimation. length(x1) = %d, length(x2) = %d\n\t\t", length(x1), length(x2) )
    myfilt = DSP.FIRFilter(h, 1//decimation)
    @timeifinteractive begin
        y1 = DSP.filt(myfilt, x1)
        y2 = DSP.filt(myfilt, x2)
    end
    statefulResult = [y1; y2]
    @test naiveResult ≈ statefulResult

    @printfifinteractive( "\n\tDSP.filt decimation. Piecewise for first %d inputs.\n\t\t", length(x1) )
    DSP.reset!(myfilt)
    y1 = similar( x, 0 )
    @timeifinteractive begin
        for i in 1:length(x1)
            append!(y1, DSP.filt(myfilt, x1[i:i]))
        end
        y2 = DSP.filt(myfilt, x2)
    end
    piecewiseResult = [y1; y2]
    @test ≈(naiveResult, piecewiseResult, atol=sqrt(eps(real(one(eltype(x))))))

    DSP.reset!(myfilt)
    @test inputlength(myfilt, length(piecewiseResult)) == xLen
end

#
# Interpolation
#

function test_interpolation(h, x, interpolation)
    xLen       = length(x)
    hLen       = length(h)
    pivotPoint = min(rand(50:150), div(xLen, 4))
    x1         = x[1:pivotPoint]
    x2         = x[pivotPoint+1:end]

    @printfifinteractive( "\n\n" )
    @printfifinteractive( "_ _  _ ___ ____ ____ ___  _    ____ ____ ___ _ ____ _  _ \n" )
    @printfifinteractive( "| |\\ |  |  |___ |__/ |__] |    |  | |__|  |  | |  | |\\ | \n" )
    @printfifinteractive( "| | \\|  |  |___ |  \\ |    |___ |__| |  |  |  | |__| | \\| \n" )
    @printfifinteractive( "\nTesting interpolation, h::%s, x::%s. xLen = %d, hLen = %d, interpolation = %d", typeof(h), typeof(x), xLen, hLen, interpolation )

    @printfifinteractive( "\n\tNaive interpolation with Base.filt\n\t\t")
    @timeifinteractive begin
        xZeroStuffed = zeros(eltype(x), xLen * interpolation)
        for n = 0:xLen-1;
            xZeroStuffed[n*interpolation+1] = x[n+1]
        end
        naiveResult = filt(h, one(eltype(h)), xZeroStuffed)
    end

    @printfifinteractive( "\n\tDSP.filt( h, x, %d//1 )\n\t\t", interpolation )
    @timeifinteractive statelesResult = DSP.filt( h, x, interpolation//1 )
    @test naiveResult ≈ statelesResult

    @printfifinteractive( "\n\tDSP.filt interpolation. length(x1) = %d, length(x2) = %d\n\t\t", length(x1), length(x2) )
    myfilt = DSP.FIRFilter( h, interpolation//1 )
    @timeifinteractive begin
        y1 = DSP.filt(myfilt, x1)
        y2 = DSP.filt(myfilt, x2)
    end
    statefulResult = [y1; y2]
    @test naiveResult ≈ statefulResult

    DSP.reset!(myfilt)
    @test inputlength(myfilt, length(statefulResult)) == xLen

    @printfifinteractive( "\n\tDSP.filt interpolation. Piecewise for first %d inputs\n\t\t", length(x1) )
    DSP.reset!(myfilt)
    y1 = similar(x, 0)
    @timeifinteractive begin
        for i in 1:length(x1)
            append!(y1, DSP.filt(myfilt, x1[i:i]))
        end
        y2 = DSP.filt(myfilt, x2)
    end
    piecewiseResult = [y1; y2]
    @test ≈(naiveResult, piecewiseResult, atol=sqrt(eps(real(one(eltype(x))))))
end

#
# Rational resampling
#

function test_rational(h, x, ratio)
    xLen       = length(x)
    hLen       = length(h)
    pivotPoint = min(rand(50:150), div(xLen, 4))
    x1         = x[1:pivotPoint]
    x2         = x[pivotPoint+1:end]
    upfactor   = numerator(ratio)
    downfactor = denominator(ratio)
    resultType = promote_type(eltype(h), eltype(x))

    @printfifinteractive( "\n\n" )
    @printfifinteractive( "      ____ ____ ___ _ ____ _  _ ____ _    \n" )
    @printfifinteractive( "      |__/ |__|  |  | |  | |\\ | |__| |    \n" )
    @printfifinteractive( "      |  \\ |  |  |  | |__| | \\| |  | |___ \n" )
    @printfifinteractive( "                                          \n" )
    @printfifinteractive( "____ ____ ____ ____ _  _ ___  _    _ _  _ ____\n" )
    @printfifinteractive( "|__/ |___ [__  |__| |\\/| |__] |    | |\\ | | __\n" )
    @printfifinteractive( "|  \\ |___ ___] |  | |  | |    |___ | | \\| |__]\n" )
    @printfifinteractive( "\n\nTesting rational resampling, h::%s, x::%s. xLen = %d, hLen = %d, ratio = %d//%d", string(typeof(h)), string(typeof(x)), xLen, hLen, upfactor, downfactor )

    @printfifinteractive( "\n\tNaive rational resampling\n\t\t")
    @timeifinteractive naiveResult = naivefilt(h, x, ratio)

    @printfifinteractive( "\n\tDSP.filt( h, x, %d//%d )\n\t\t", upfactor, downfactor )
    @timeifinteractive statelesResult = DSP.filt(h, x, ratio)
    @test naiveResult ≈ statelesResult

    @printfifinteractive( "\n\tDSP.filt rational resampling. length(x1) = %d, length(x2) = %d\n\t\t", length(x1), length(x2) )
    myfilt = DSP.FIRFilter(h, ratio)
    @timeifinteractive begin
        s1 = DSP.filt(myfilt, x1)
        s2 = DSP.filt(myfilt, x2)
    end
    statefulResult = [s1; s2]
    @test naiveResult ≈ statefulResult

    @printfifinteractive( "\n\tDSP.filt rational. Piecewise for all %d inputs\n\t\t", length( x ) )
    reset!(myfilt)
    y1 = similar(x, 0)
    @timeifinteractive begin
        for i in 1:length(x)
            append!(y1, DSP.filt(myfilt, x[i:i]))
        end
    end
    piecewiseResult = y1
    @test naiveResult ≈ piecewiseResult

    DSP.reset!(myfilt)
    @test inputlength(myfilt, length(piecewiseResult)) == xLen
end

#
# Arbitrary resampling
#

function test_arbitrary(Th, x, resampleRate, numFilters)
    cutoffFreq      = 0.45
    transitionWidth = 0.05
    h               = digitalfilter(Lowpass(cutoffFreq, fs=numFilters), FIRWindow(transitionwidth=transitionWidth/numFilters)) .* numFilters
    h               = convert(Vector{Th}, h)
    myfilt          = DSP.FIRFilter(h, resampleRate, numFilters)
    xLen            = length(x)

    @printfifinteractive( "____ ____ ___      ____ ____ ____ ____ _  _ ___  _    _ _  _ ____\n" )
    @printfifinteractive( "|__| |__/ |__]     |__/ |___ [__  |__| |\\/| |__] |    | |\\ | | __\n" )
    @printfifinteractive( "|  | |  \\ |__] .   |  \\ |___ ___] |  | |  | |    |___ | | \\| |__]\n" )
    @printfifinteractive( "\n\nh::%s, x::%s, rate = %f, Nϕ = %d, xLen = %d, ", string(typeof(h)), string(typeof(x)), resampleRate, numFilters, length(x) )

    @printfifinteractive( "\n\tNaive arbitrary resampling\n\t\t" )
    @timeifinteractive naiveResult = naivefilt(h, x, resampleRate, numFilters)

    @printfifinteractive( "\n\tStateless arbitrary resampling\n\t\t" )
    @timeifinteractive statelessResult = DSP.filt(h, x, resampleRate, numFilters)

    @printfifinteractive( "\n\tStateful arbitrary resampling\n\t\t" )
    @timeifinteractive statefulResult = DSP.filt(myfilt, x)

    # DSP.reset!(myfilt)
    # TODO: figure out why this fails
    # @test inputlength(myfilt, length(statefulResult)) == xLen

    @printfifinteractive( "\n\tPiecewise arbitrary resampling\n\t\t" )
    reset!(myfilt)
    piecwiseResult = eltype(x)[]
    sizehint!(piecwiseResult, ceil(Int, length(x)*resampleRate))
    @timeifinteractive for i in 1:length(x)
        thisY = filt(myfilt, x[i:i])
        append!(piecwiseResult, thisY)
    end

    commonLen = min(length(naiveResult), length(statelessResult), length(statefulResult), length(piecwiseResult))
    resize!(naiveResult, commonLen)
    resize!(statelessResult, commonLen)
    resize!(statefulResult, commonLen)
    resize!(piecwiseResult, commonLen)

    @test naiveResult ≈ statelessResult
    @test naiveResult ≈ statefulResult
    @test naiveResult ≈ piecwiseResult
end

#
# Run the tests
#

@testset "interp.=$interpolation, dec.=$decimation, Th=$Th, Tx=$Tx" for interpolation in [1, 5, 14, 23],
            decimation in [1, 9, 17, 21],
                Th in [Float32, Float64],
                    Tx in [Float32, Float64, ComplexF32, ComplexF64]

    h     = rand(Th, rand(16:128))
    xLen  = rand(200:300)
    xLen  = xLen-mod(xLen, decimation)
    x     = rand(Tx, xLen)
    ratio = interpolation//decimation
    if ratio == 1
        test_singlerate(h, x)
    end
    if decimation != 1
        test_decimation(h, x, decimation)
    end
    if interpolation != 1
        test_interpolation(h, x, interpolation)
    end
    if numerator(ratio) == interpolation && denominator(ratio) == decimation && numerator(ratio) != 1 && denominator(ratio) != 1
        test_rational(h, x, ratio)
        if Tx in [Float32, ComplexF32]
            test_arbitrary(Th, x, convert(Float64, ratio)+rand(), 32)
        end
    end
    if decimation == 1
        test_rational(h, x, interpolation)
    end
end